K M M Yunusova

S L L Bayliss

T Chanelière

V Derkach

J E E Anthony

Alexei D Chepelianskii

L R R Weiss

Spin fine-structure

published or not. The documents may come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

Bi-exciton states are key excited-state species in a range of nano-structured materials from quantumconfined inorganic systems [START_REF] Bryant | Biexciton binding in quantum boxes[END_REF][START_REF] Chen | Biexciton quantum coherence in a single quantum dot[END_REF][START_REF] Hu | Biexcitons in semiconductor quantum dots[END_REF] to synthetic molecular structures [START_REF] Baldo | Transient analysis of organic electrophosphorescence. ii. transient analysis of triplet-triplet annihilation[END_REF][START_REF] Klimov | Biexcitons in π-conjugated oligomers: Intensity-dependent femtosecond transient-absorption study[END_REF][START_REF] Masui | Analysis of exciton annihilation in high-efficiency sky-blue organic lightemitting diodes with thermally activated delayed fluorescence[END_REF]. In organic semiconductors the excitonpair is an intermediate in both the process of singlet fission [START_REF] Smith | Singlet fission[END_REF][START_REF] Singh | Laser generation of excitons and fluorescence in anthracene crystals[END_REF][START_REF] Swenberg | Bimolecular radiationless transitions in crystalline tetracene[END_REF] -the formation of a pair of spin-1 (triplet) excitons from an initial spin-0 (singlet) excitation -and its reverse process, triplet-triplet annihilation [START_REF] Singh-Rachford | Photon upconversion based on sensitized triplet-triplet annihilation[END_REF]. While singlet fission is of particular interest for photovoltaics [START_REF] Rao | Harnessing singlet exciton fission to break the shockley-queisser limit[END_REF][START_REF] Hanna | Solar conversion efficiency of photovoltaic and photoelectrolysis cells with carrier multiplication absorbers[END_REF][START_REF] Tayebjee | Thermodynamic limit of exciton fission solar cell efficiency[END_REF] where it has been shown to increase efficiencies of solar energy harvesting beyond traditional limits [START_REF] Einzinger | Sensitization of silicon by singlet exciton fission in tetracene[END_REF][START_REF] Congreve | External quantum efficiency above 100% in a singlet-exciton-fission-based organic photovoltaic cell[END_REF], triplet-triplet annihilation enables the up-conversion of light via the conversion of two low-energy triplet excitons to one higher energy singlet excitation [START_REF] Singh-Rachford | Photon upconversion based on sensitized triplet-triplet annihilation[END_REF][START_REF] Cheng | On the efficiency limit of triplet-triplet annihilation for photochemical upconversion[END_REF]. This process is of particular interest for catalysis [START_REF] Ravetz | Photoredox catalysis using infrared light via triplet fusion upconversion[END_REF][START_REF] Khnayzer | Upconversion-powered photoelectrochemistry[END_REF], photovoltaics [START_REF] Cheng | On the efficiency limit of triplet-triplet annihilation for photochemical upconversion[END_REF][START_REF] Gray | Triplet-triplet annihilation photonupconversion: towards solar energy applications[END_REF], and bio-imaging [START_REF] Liu | Blue-emissive upconversion nanoparticles for low-power-excited bioimaging in vivo[END_REF][START_REF] Liu | Highly photostable near-ir-excitation upconversion nanocapsules based on triplet-triplet annihilation for in vivo bioimaging application[END_REF]. In both fission and annihilation the bi-exciton state mediates the transition between one excited-state and two dissociated excited-states [START_REF] Burdett | Quantum beats in crystalline tetracene delayed fluorescence due to triplet pair coherences produced by direct singlet fission[END_REF][START_REF] Merrifield | Theory of magnetic field effects on the mutual annihilation of triplet excitons[END_REF][START_REF] Zimmerman | Singlet fission in pentacene through multi-exciton quantum states[END_REF].

The wavefunction of these transient, intermediate pairs remains challenging to probe. Purely optical characterization of bi-excitons can be ambiguous as their optical signatures typically overlap with those of singly-excited states. By contrast, light-induced spin resonance has played a key role in showing, unexpectedly, that in several molecular systems, singlet fission produces a long-lived, bound spin-2 bi-exciton, termed the quintet state [START_REF] Tayebjee | Quintet multiexciton dynamics in singlet fission[END_REF][START_REF] Weiss | Strongly exchange-coupled triplet pairs in an organic semiconductor[END_REF][START_REF] Lubert-Perquel | Multiple quintets via singlet fission in ordered films at room temperature[END_REF]. The observation of the quintet state reflects a dominant exchange interaction between the triplets within a pair. This exchange interaction (parametrized by exchange energy J) separates the pure singlet (S = 0) from the triplet (S = 1), and quintet (S = 2) pairings by J and 3J respectively [START_REF] Bayliss | Site-selective measurement of coupled spin pairs in an organic semiconductor[END_REF]. The spin signatures of this strongly coupled state provide an unambiguous marker of the bi-exciton state [START_REF] Weiss | Strongly exchange-coupled triplet pairs in an organic semiconductor[END_REF][START_REF] Bayliss | Site-selective measurement of coupled spin pairs in an organic semiconductor[END_REF].

Where are such pairs located in a solid-state system with many possible pair configurations? Here the sensitivity to inter-spin distance and orientation of the spin dipole-dipole interaction provides a window into the geometry of the multi-exciton state. Indeed, spin-spin dipolar interactions are used as a ruler for nano-scale distances in widespread applications of spin-labels to biological and synthetic systems [START_REF] Berliner | Spin labeling: theory and applications[END_REF]. Typically non-native spin-probes are required to take advantage of the dipolar interaction between synthetic radical-bearing ligands attached to biologically relevant structures. Conveniently, the dipolar interactions between the two triplet states in the quintet Hamiltonian provides a native probe of its spatial confinement and orientation [START_REF] Weil | Electron paramagnetic resonance: elementary theory and practical applications[END_REF][START_REF] Benk | Theory of two coupled triplet states: Application to bicarbene structures[END_REF].

These dipolar interactions are reflected in the fine structure in the quintet (and triplet) spin manifold: the energetic splitting within the spin-manifold even in the absence of a magnetic field. The zero-field splitting (ZFS) Hamiltonian in its principal axes takes the form

Ĥzfs /h = S T •D•S = D( Ŝ2 z - 1 3 S(S+1))+E( Ŝ2 x -Ŝ2 y ) (1)
where D is the dipolar tensor (D-tensor) and S is the relevant vector of spin operators ( with total spin S = 1, 2 for triplet, quintet states) defined along the principal axes (x, ŷ, ẑ) of the D-tensor. These directions are defined by the symmetry of the spin distribution with the ẑ-axis defined by the axis of maximum absolute magnitude of dipolar coupling. The strength of the dipolar coupling parameters (D, E) reflect the degree of confinement and asymmetry of the spin-system as we describe in detail in the next section.

Experimentally these parameters (D, E and x, y, ẑ) can be extracted from magnetic resonance spectra, which can be detected via microwave absorption in a cavity or via spin-dependent electrical or optical transitions. In this work we use optically detected magnetic resonance (ODMR) as shown in Fig. 1a, which combines three key advantages for probing the geometry of the transient bi-exciton state: (1) sensitivity of photoluminescencedetection to the recombination of intermediate states,

(2) detection sensitivity from optical (versus microwave) read-out, and (3) the sensitivity to dipolar interactions (on the MHz-GHz scale) of magnetic resonance.

II. THEORY

In this section we establish the relationship between the triplet pair geometry and the resulting measurable S = 2 bi-exciton fine structure (D Q , E Q and the principal axes xq , ŷq , ẑq ), which we then go on to determine experimentally. Previous measurements of quintet states in various materials have found that D Q ∼ D T /3 [START_REF] Tayebjee | Quintet multiexciton dynamics in singlet fission[END_REF][START_REF] Weiss | Strongly exchange-coupled triplet pairs in an organic semiconductor[END_REF], yet this relation holds only in the case of co-linear triplet D-tensors and in the absence of dipolar coupling between the excitons, and so acts as an approximate guide to identifying quintet spectra [START_REF] Benk | Theory of two coupled triplet states: Application to bicarbene structures[END_REF].

For triplet excitons in organic semiconductors D T and E T are set by the spatially averaged spin-spin interaction of electron and hole, reflecting triplet states localized on a single molecule in the material used here, TIPS-tetracene [START_REF] Bayliss | Localization length scales of triplet excitons in singlet fission materials[END_REF][START_REF] Bayliss | Geminate and nongeminate recombination of triplet excitons formed by singlet fission[END_REF]. As shown in Fig. 1b, the ẑ-axis of maximal spinspin coupling is out of the plane of the TIPS-tetracene molecule and the x-and ŷ-axes are in the plane of the molecule along the long and short axes.

The S = 2 triplet-pair fine-structure depends on the underlying orientation of two such triplets localized on two molecules (labeled here A and B). We assume that each triplet has the same zero-field parameters (D T , E T ) and differ only in orientation and position. We define the principal axes of the triplet state on molecule A as (x a , ŷa , ẑa ) and on molecule B as (x b , ŷb , ẑb ) with the unit vector between them given by ûab = r ab /| r ab |. The zerofield Hamiltonian of the pair in the uncoupled basis is then given by

Ĥ(1⊗1) zf s /h = i=a,b S T i •D i T •S i -Γ(û ab •S a )(û ab •S b )+JS a •S b (2) where Γ = 3µ0µ 2 B g 2 4π| r ab | 3
gives the strength of the dipolar interaction with µ 0 the magnetic permeability of free-space, µ B the Bohr magneton, g the g-factor, and | r ab | the intertriplet distance. In the limit of strong exchange coupling (J D T ) the Hamiltonian is approximately diagonal in the coupled spin-basis defined by the states of pure total spin [START_REF] Bayliss | Site-selective measurement of coupled spin pairs in an organic semiconductor[END_REF][START_REF] Bayliss | Spin signatures of exchange-coupled triplet pairs formed by singlet fission[END_REF][START_REF] Yago | Magnetic field effects on triplet pair generated by singlet fission in an organic crystal: application of radical pair model to triplet pair[END_REF]. Projecting the above Hamiltonian into the S = 2 subspace and converting to the coupled basis gives the quintet zero-field Hamiltonian as

Ĥ(2) zf s /h = S T • D Q • S (3)
where S = ( Ŝx , Ŝy , Ŝz ) are the Pauli spin operators for total spin-2. The quintet zero-field tensor D Q in terms of the underlying triplet fine structure, inter-triplet distance, and dipolar interaction is given by

D Q = D T 6 ( i=a,b ẑi ẑT i - 2 3 Î3 ) (4)
+ E T 6 i=a,b (x i xT i -ŷi ŷT i ) - Γ 3 (û ab ûT ab - 1 3 Î3 )
where the full derivation is given in the Supporting Information (Appendix A). The eigenvalues of

D Q are or- dered |D y | < |D x | < |D z | , which gives D Q = 3 2 D z , E Q = 1 2 (D y -D x )
, and the principal axes xq , ŷq , ẑq (the eigenvectors of D Q ). In the following sections we will apply this framework to link the measured values and orientation of the D-tensors of the triplet and quintet states to their localization in the crystal structure of TIPS-tetracene.

III. FINE STRUCTURE TENSOR PARAMETERS OF THE QUINTET STATE

TIPS-tetracene is a solution-processable singlet fission material of interest for its high singlet fission efficiency [START_REF] Stern | Vibronically coherent ultrafast tripletpair formation and subsequent thermally activated dissociation control efficient endothermic singlet fission[END_REF][START_REF] Stern | Identification of a triplet pair intermediate in singlet exciton fission in solution[END_REF]. Structured with side-chain modification of the canonical fission molecule, tetracene, as shown in Fig. 1b, TIPS-tetracene crystallizes with a distinct molecular packing shown in Fig. 1c. The unit cell contains four orientationally inequivalent molecules yielding six possible nearest-neighbor pair sites within a unit cell in addition to non-crystalline defect sites. Both biexcitons and free triplet excitons have been identified in this material with transient electron spin resonance, but the nature of the local site housing the bound pair remains a key question. The theoretical framework above demonstrates how the parameters of the quintet fine structure (its principal values D Q , E Q and the orientation of the principal axes xq , ŷq , ẑq ) are specific to the geometry of the underlying triplet pair and so each distinct potential pair site can be identified by its unique fine structure parameters in a single crystal. (Note that these parameters are calculated explicitly from the crystal structure for each of the six possible pair states and reported in Section V).

With this motivation we use a macroscopic crystal (∼mm-scale single-crystalline domain) and measure the principal values of the D-tensors of the triplet and quintet states using broadband zero-field ODMR and then apply a magnetic field to find the orientation of their principal axes in the laboratory frame in Section IV. The experimental setup is shown in Fig. 1a and includes 532 nm continuous-wave (CW) light excitation, microwave radiation (B 1 ) with variable frequency delivered by a broadband copper strip-line and a static magnetic field (B 0 ). The ODMR signal is measured by lock-in detection of microwave-induced changes in photoluminescence (see the Methods section for further details).

We first perform standard fixed-frequency (9 GHz) field-swept ODMR to confirm that, in agreement with previous measurements using transient electron spin resonance, we observe two pairs of spin-transitions consistent with the ∆m = ±1 transitions of the S = 1 triplet exciton (T ± ) and the ∆m = ±1 transitions of the S = 2 quintet state (Q ± ), as shown in Fig. 2b [START_REF] Weiss | Strongly exchange-coupled triplet pairs in an organic semiconductor[END_REF][START_REF] Lubert-Perquel | Multiple quintets via singlet fission in ordered films at room temperature[END_REF]. The corresponding energy-level diagrams and transitions fields are shown in Fig. 2a.

These observed high-field transitions are correlated with zero-field transitions (magnetic field strength B 0 = 0, Fig. 2c,d), measured for the first time here. These zero-field spectra allow us to extract the zero-field splitting parameters for the triplet and quintet states directly.

In the absence of a magnetic field, the energy eigenstates of triplets and quintets are determined by their spin-spin dipolar interactions. The energy eigenstates are given by the B 0 = 0 spin-eigenstates. Microwave transitions at zero-field are determined by the energy-splitting between these zero-field states. As shown in the theoretical sketch in Fig. 2c, the triplet energy levels are separated in energy by |D T | and further split by |2E T | between the two upper eigenstates, where D T and E T are the zero-field splitting parameters for the triplet. Zero-field ODMR resonances therefore occur at microwave frequencies ν = D T ± E T .

The quintet eigenstates in zero-field are likewise defined by the eigenstates of the quintet fine-structure tensor. As shown in Fig. 2c, the three lowest quintet levels are split by |D Q | from the ground state to the first two states with a further splitting of |6E Q | between those two upper levels. This leads to ODMR transition frequencies at ν = D Q ± 3E Q . Note that the previously reported D-parameters for TIPS-tetracene are D T ∼ 1.4 GHz and D Q ∼ D T /3 [START_REF] Weiss | Strongly exchange-coupled triplet pairs in an organic semiconductor[END_REF][START_REF] Bayliss | Geminate and nongeminate recombination of triplet excitons formed by singlet fission[END_REF]. However previous experiments did not determine E-parameters and molecular frame orientations which, as we will show, are needed to identify the contribution of triplet dipole-dipole interaction to the quintet fine structure tensor.

The magnetic resonance spectra of triplets and biexcitons can be distinguished in ODMR due to the difference in lifetime of the two species [START_REF] Weiss | Strongly exchange-coupled triplet pairs in an organic semiconductor[END_REF][START_REF] Bayliss | Geminate and nongeminate recombination of triplet excitons formed by singlet fission[END_REF]. The microwave amplitude modulation frequency (137 Hz) is chosen to correspond with the inverse lifetime of the longer-lived species (triplets), which then appear with equal amplitude on the in-phase (X-channel) and out-of-phase (Ychannel) lock-in channels. The signal from shorter lived bi-excitons appears only on the X-channel and can be isolated by subtracting X and Y channels. The zero-field X-and Y-channel ODMR spectra are plotted in Fig. 2d in black (X-channel) and blue (Y-channel). The transitions on the Y-channel are consistent with triplets with |D T | = 1.4 GHz and |E T | = 14 MHz (Fig. 2d). Transitions in the frequency region expected for the quintet only appear on the X-channel and give |D Q | = 477 MHz and |E Q | = 22 (Fig. 2d inset). The sign of D T was determined in Ref. [START_REF] Weiss | Strongly exchange-coupled triplet pairs in an organic semiconductor[END_REF] to be positive, which in turn determines E T to be negative. The resolution of the splitting afforded by this technique allows extraction of E Q and E T , highlighting the sensitivity of the zero-field approach. The measurement of the E-parameters here is made possible by the reduced linewidths observed at zero-field relative to previous measurements under non-zero magnetic field.

IV. EXTRACTION OF BI-EXCITON FINE STRUCTURE ORIENTATION

Having extracted the principal components of the triplet and quintet fine-structure tensors at zero-field, we now map the resonant frequencies as a function of mag- netic field to determine the corresponding orientations of the principal axes. The experimental ODMR maps for quintet and triplet states are shown in Fig. 3c,f (separated by lock-in detection phase as described in the previous section).

The orientation of the quintet fine-structure axes in the lab frame is obtained by fitting these maps with the spin transitions predicted by the fine structure parameters determined at zero-field (see Sec. III) with the orientation as input. The polar angle θ and azimuthal angle φ parametrize the orientation of the principal axes relative to the magnetic field in the laboratory frame as shown in Fig. 3b,e. There are 10 possible transitions between the five quintet spin sub-levels (Fig. 3a), which are overlaid on the quintet ODMR map (Fig. 3c). It should be noted that the visibility of transitions depends on populations and selection rules, and transitions 3,4,9, and 6 are not observed experimentally. The observed resonances cannot be fit by a spin-1 state, which further confirms the assignment (see SI). The quintet state is oriented with fixed θ q = 90 ± 5 • between ẑq and B 0 and φ q = 30 ± 5 • between xq and B 0 .

The evolution of the ∼ 1.4 GHz zero-field transitions with field, shown in Fig. 3f are consistent with θ t, = 0 • and θ t = 90 • with the latter perpendicular component more prominent (calculated transitions for θ t = 90 • are plotted in black) and the θ t = 0 • consistent with the presence of a weak powder background. The dominant θ t = 90 • triplet orientation correlates with the high-field spectrum (Fig. 2b): triplet peaks are separated by ∆B 0 ∼ D/gµ B , which occurs when θ t ∼ 90 • , whereas no peaks are observed for θ t ∼ 0 (separation ∆B 0 ∼ 2D/gµ B ) . [Note that φ t could not be extracted reliably as E T D T , and this parameter was not required for subsequent analysis because the triplet states are nearly axially symmetric (i.e., as E T ≈ 0, E Q and D Q do not depend significantly on the relative orientation of xt and ŷt ).] As the D-tensor values and orientations in the laboratory frame are obtained from a bulk crystalline sample, we can now compare them with the theoretically predicted D-tensors in the TIPS-tetracene crystal structure.

V. PAIR SITES IN THE TIPS-TETRACENE CRYSTAL STRUCTURE

We now compare the experimental observations to the quintet fine-structure theoretically predicted in the TIPS-tetracene crystal structure. As shown in Fig. 1c, in the TIPS-tetracene unit-cell molecules are packed with two approximately parallel dimers each rotated with respect to the other by ∼ 90 • about the ĉ axis. There are six potential nearest-neighbor dimer configurations and for each we can calculate the expected D Q , E Q parameters using Eq. ( 4), the molecular triplet axes shown in Fig. 1b and the intermolecular distances extracted from the crystal structure determined by X-ray diffraction [START_REF] Eaton | Experimental crystal structure determination[END_REF]. These values are summarized in Table I where we have taken nearest-neighbor dimer configurations. (We note that we have taken a point-dipole approximation for each triplet localized at the center of each molecule in the pair which is likely the main source of error in the predicted fine structure parameters.)

The observed quintet and triplet signatures and extracted angles of θ q ∼ 90 • and φ q ∼ 30 • are consistent with exchange-coupled triplets localized on dimers 1,2 and 3,4 (highlighted in blue in Table I) where molecules in the unit cell are labeled 1 -4 in Fig. 1c and the calculated angles are for B 0 ĉ as this orientation is consistent with known crystal facets, both the quintet and triplet transition maps, and the 9 GHz field-swept spectrum presented in Fig. 1c (see SI). Interestingly, these dimer configurations correspond to the most closely πstacked dimers in the crystal suggesting that this geometry is favorable for binding of the multi-exciton state. The extracted local quintet fine structure is visualized in Fig. 4 where the full quintet and triplet dipolar interactions are shown with respect to the magnetic field in the lab frame and crystallographic axes, summarizing the local structure of quintet and triplet states in TIPStetracene and their relation to intermolecular geometry.

VI. CONCLUSIONS

We have developed the theoretical and experimental approaches to connect inter-molecular structure and biexciton geometry. In the context of organic semiconductors, understanding this connection presents a path toward molecular engineering of bi-excitons. We find that the triplet pair is localized on the closest π-stacked dimers of the crystal structure. As the fine structure is consistent with minimal modification from the ground-state crystal structure, these results suggest that geometric reorganization is negligible in the quintet bi-exciton excitedstate. This description of the geometry of the triplet pair therefore sets the foundation for orientation-selective, alloptical experiments where time-resolved measurements would allow investigation of the transient localization and molecular reorganization of the pair-state.

Further, we have shown that the energy-level structure, determined by the inter and intra-triplet dipolar interactions, is highly sensitive to the underlying molecular structure. This sensitivity of the effective spin Hamiltonian to intermolecular geometry further opens up opportunities to tune spin parameters via crystal structure engineering of spin-based organic devices. Tuning spin energy levels and zero-field interactions is critical in spinbased applications, for example in tuning the amplification frequency in the newly emerging class of organic molecular masers [START_REF] Oxborrow | Roomtemperature solid-state maser[END_REF][START_REF] Bogatko | Molecular design of a roomtemperature maser[END_REF] or in minimizing spin relaxation and dephasing [START_REF] Solomon | Relaxation processes in a system of two spins[END_REF][START_REF] Shiddiq | Enhancing coherence in molecular spin qubits via atomic clock transitions[END_REF]. Dipolar spin interactions thereby provide a route toward engineering the spin levelstructure in organic materials in future spin-based quantum technologies.

Figure 1 .

 1 Figure 1. Broadband ODMR of the triplet-pair state. (a) Schematic of the experimental set-up. Crystalline samples of TIPS-tetracene were optically illuminated (532 nm, continuous-wave (CW) excitation) under amplitudemodulated microwave excitation (B1) using a broadband strip-line (0 -9 GHz) in liquid helium (4 K). Photoluminescence was collected (via 1 mm core optical fiber) to detect the microwave-induced change in photoluminescence (ODMR) as a function of both microwave frequency and static magnetic field (B0) with B0 ⊥ B1. (b) Molecular structure of TIPStetracene and corresponding principal axes of the intra-triplet dipolar interaction (xt, ŷt, ẑt). (c) Solid-state crystal structure of TIPS-tetracene with four rotationally inequivalent molecules per unit cell labelled 1-4 and unit cell axes (â, b, ĉ).

Figure 2 .

 2 Figure 2. Field-swept and zero-field ODMR of the triplet pair state. (a) Energy levels of the triplet and quintet m = 0, ±1 sub-levels as a function of field with transitions at 9 GHz marked with arrows to correspond to experimentally observed transitions in (b). (b) ODMR spectrum at 9 GHz showing inner quintet transitions (Q ± ) and the outer triplet transitions (T ± ). (c) Energy level diagram of triplet and quintet zero-field spin sub-levels. (d) Zero-field ODMR spectrum where the triplet transitions are split by the zero-field ET -parameter and with centre frequency set by the DT -parameter. Inset: Low frequency quintet transitions marked with corresponding zero-field splitting parameters (DQ, EQ) of the quintet pair-state.

Figure 3 .

 3 Figure 3. Fine-structure tensor from from broadband ODMR. (a,d) Energy level diagram for the quintet (a) and triplet (d) states as a function of magnetic field. Arrows indicate potential transitions, corresponding to lines on Fig. 3c and 3f respectively. (b,e) Schematic representation of the orientation of B in the quintet (b) and triplet (e) fine-structure axes. (c) ODMR transitions associated with the quintet state with overlay of simulated transitions. Signal has been isolated by subtracting a scaled out-of-phase (Y)-channel signal from the in-phase signal to remove triplet contributions. Black lines show simulations given for θq = 90 • and φq = 30 • with an uncertainty of ±5 • . (f) Y-channel (out-of-phase) ODMR map of ODMR transitions attributed to the triplet state with overlay of calculated transitions for the triplet state in black with θt = 90 • .

Figure 4 .

 4 Figure 4. Local geometry of the quintet fine-structure tensor in the TIPS-tetracene crystal with an angle of θq = 91.6 • and φq = 30.6 • , where the crystal is oriented with B0 ĉ with two co-planar triads of vectors, (xq, ŷq, B0) and (ẑq, r ab , B0).

Acknowledgments We are thankful for insightful discussions with H. Bouchiat and R.H. Friend and acknowledge support from ANR SPINEX and Labex ANR-10-LABX-0039-PALM. L.R. Weiss acknowledges support from the Clare College Junior Research Fellowship.

Author Contributions The data were measured and analyzed by K.M.Y., L.R.W., and A.D.C. J.E.A. provided materials. All authors contributed to planning the experiment, discussions and the preparation of the manuscript.

Author Information The authors declare no competing financial interests. Correspondence should be addressed to A.D.C. (chepelianskii@lps.u-psud.fr) and L.R.W (lrweiss@uchicago.edu).

Methods Optically detected magnetic resonance is performed 532 nm CW laser illumination via 1 mm-diameter optical fiber (Fig. 1a) in liquid helium at 4 K. An amplitude modulated microwave field (B 1 ) in a frequency range 0 -9 GHz was applied using a broadband copper stripline. The sample holder was fixed inside a superconducting solenoid. Photoluminescence was collected (via the same optical fiber and avalanche photodiode) to detect the microwave-induced change in photoluminescence (ODMR) as a function of both microwave frequency and static magnetic field (B 0 ) with B 0 ⊥ B 1 .


Kimon Antonakopoulos
email: kimon.antonakopoulos@inria.fr

E Veronica Belmega
email: belmega@ensea.fr

Panayotis Mertikopoulos
email: panayotis.mertikopoulos@imag.fr

An Adaptive Mirror-Prox Algorithm for Variational Inequalities with Singular Operators

Lipschitz continuity is a central requirement for achieving the optimal O(1/T ) rate of convergence in monotone, deterministic variational inequalities (a setting that includes convex minimization, convex-concave optimization, nonatomic games, and many other problems). However, in many cases of practical interest, the operator defining the variational inequality may exhibit singularities at the boundary of the feasible region, precluding in this way the use of fast gradient methods that attain this optimal rate (such as Nemirovski's mirror-prox algorithm and its variants). To address this issue, we propose a novel regularity condition which we call Bregman continuity, and which relates the variation of the operator to that of a suitably chosen Bregman function. Leveraging this condition, we derive an adaptive mirror-prox algorithm which attains the optimal O(1/T ) rate of convergence in problems with possibly singular operators, without any prior knowledge of the degree of smoothness (the Bregman analogue of the Lipschitz constant). We also show that, under Bregman continuity, the mirror-prox algorithm achieves a O(1/ √ T ) convergence rate in stochastic variational inequalities.

Local Bregman functions and the associated divergence. The notion of a dual local norm presented above will be our principal measure of distance in V * . To proceed, we will also need to adapt the notion of a Bregman (or distance-generating) function on X :

Definition 2. Let • x be a local norm on X . We say that h : V → R is a Bregman function on X if:

1. h is proper, l.s.c., convex, and cl(dom h) = cl(X ).

2. The subdifferential of h admits a continuous selection, i.e., a continuous function ∇h such that ∇h(x) ∈ ∂h(x) for all x ∈ X • ≡ dom ∂h.

Introduction

The seminal introduction of generative adversarial networks (GANs) [START_REF] Goodfellow | Generative adversarial nets[END_REF] has ushered in a new optimization paradigm in deep learning: instead of focusing on the minimization of an empirical loss function, GAN training hinges on a zero-sum game between a generator and a discriminator. In fact, in many cases GAN training goes even beyond the min-max setting, either because there are more than two networks involved, or because the objectives of the generator and the discriminator are not entirely opposed -e.g., as in the widely used ACGAN framework of Odena et al. [START_REF] Odena | Conditional image synthesis with auxiliary classifier GANs[END_REF]. In these cases, the most compact way of representing the problem's training landscape is by means of a variational inequality (VI).

Tracing their origins to the work of Stampacchia [START_REF] Stampacchia | Formes bilineaires coercitives sur les ensembles convexes[END_REF] on the Signorini problem, variational inequalities have since found a broad range of applications in physics, engineering, economics -and, more recently, machine learning. One of the main reasons for their extensive applicability is that they comprise a flexible optimization framework which can simultaneously account for loss function minimization, saddle-point, game-theoretic, and fixed point problems. As a result, there has been considerable interest in the literature to develop optimal algorithms for solving VI problems; for an appetizing introduction, see [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF] and references therein.

One of the most widely studied methods for this purpose is ordinary gradient descent -also known as the forward-backward (FB) algorithm in operator theory [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]. 1 In monotone, deterministic variational inequalities, the convergence of the method is guaranteed under a condition known as cocoercivity. By the Baillon-Haddad theorem, if the operator defining the variational inequality is a gradient field (i.e., in loss minimization problems), this condition is equivalent to Lipschitz smoothness of the associated loss function [START_REF] Baillon | Quelques propriétés des opérateurs angle-bornés et n-cycliquement monotones[END_REF][START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]. However, cocoercivity may fail to hold even in simple, bilinear min-max problems, in which case gradient descent provably fails to converge -see e.g., [START_REF] Gidel | Negative momentum for improved game dynamics[END_REF][START_REF] Mertikopoulos | Cycles in adversarial regularized learning[END_REF][START_REF] Mertikopoulos | Optimistic mirror descent in saddle-point problems: Going the extra (gradient) mile[END_REF] for a precise statement.

The first algorithm achieving convergence in (pseudo-)monotone variational inequalities without cocoercivity is the extra-gradient (EG) algorithm of Korpelevich [START_REF] Korpelevich | The extragradient method for finding saddle points and other problems[END_REF], which only requires Lipschitz continuity of the underlying operator. 2 The asymptotic convergence result of Korpelevich [START_REF] Korpelevich | The extragradient method for finding saddle points and other problems[END_REF] was subsequently extended by Nemirovski [START_REF] Nemirovski | Prox-method with rate of convergence O(1/t) for variational inequalities with Lipschitz continuous monotone operators and smooth convex-concave saddle point problems[END_REF] who introduced the mirror-prox (MP) algorithm, a Bregman variant of the EG algorithm with ergodic averaging. As was shown in [START_REF] Nemirovski | Prox-method with rate of convergence O(1/t) for variational inequalities with Lipschitz continuous monotone operators and smooth convex-concave saddle point problems[END_REF], the mirror-prox algorithm attains a O(1/T ) ergodic convergence rate in monotone variational inequalities with Lipschitz continuous operators, and this rate cannot be improved without further assumptions.

However, in many applications and problems of practical interest, Lipschitz continuity may also fail to hold, either because the loss profile of the problem grows too rapidly (e.g., as in support vector machines or GAN models with Kullback-Leibler losses), or because the problem exhibits singularities near the boundary of the feasible region (e.g., as in resource allocation and inverse problems). In these cases, one would still want to apply a fast method like mirror-prox, but the lack of smoothness means that there are no convergence guarantees -asymptotically, ergodically, or otherwise.

Our contributions. Our starting point is the observation that this failure stems from the fact that Lipschitz continuity of the operator is defined relative to a global norm. Because of this, the standard Lipschitz framework is not well-suited to problems with singularities or rapid growth: a global norm is oblivious to the geometry of the feasible region (and, in particular, its boundary), so it cannot capture the finer features of the problem's loss landscape.

To overcome this limitation, we introduce a novel regularity condition, which we call Bregman continuity, and which is made-to-order for the singularity landscape of the problem at hand. Specifically, instead of defining Lipschitz continuity relative to a global norm, we define it in terms of a family of local norms and a suitably chosen Bregman function. This leads to an intricate interplay between different geometric notions of distance (the Bregman divergence and the local norm), but it also introduces the flexibility required to tackle variational inequalities with singular operators.

Under this assumption, we show that the mirror-prox algorithm attains the optimal O(1/T ) convergence rate in variational inequalities with (possibly) singular operators, provided that the method is run with the same Bregman function that is used to define Bregman continuity. As in the standard Lipschitz framework, the method's convergence requires a step-size of the form γ < 1/β, where β is the Bregman constant of the operator (i.e., the Bregman analogue of the Lipschitz constant). Estimating this constant can be fairly challenging in practice (if not downright impossible), so we also introduce an adaptive mirror-prox (AMP) method which attains the same O(1/T ) rate without requiring any a priori estimation of β -essentially, the Bregman constant is learned at the same time as the problem's landscape. Finally, we provide a variant of the method for stochastic variational inequalities, and we establish a O(1/ √ T ) convergence rate in this setting. To the best of our knowledge, these are the first results of this kind in the literature.

Related work.

Owing to their optimal rate guarantees, the extra-gradient and mirror-prox algorithms have been at the forefront of an extensive literature which is impossible to adequately review here. As a purely indicative list of contributions in the Lipschitz continuous setting (and with no illusion of being comprehensive), we refer the reader to Juditsky et al. [START_REF] Juditsky | Solving variational inequalities with stochastic mirror-prox algorithm[END_REF], Chambolle and Pock [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF], Malitsky [START_REF] Malitsky | Projected reflected gradient methods for monotone variational inequalities[END_REF], Iusem et al. [START_REF] Iusem | Extragradient method with variance reduction for stochastic variational inequalities[END_REF] and Mokhtari et al. [START_REF] Mokhtari | A unified analysis of extra-gradient and optimistic gradient methods for saddle point problems: proximal point approach[END_REF] for some recent developments. Especially in learning theory, there has been a surge of interest motivated by the application of EG/MP methods to GAN training, see e.g., [START_REF] Daskalakis | Training GANs with optimism[END_REF][START_REF] Gidel | A variational inequality perspective on generative adversarial networks[END_REF][START_REF] Mertikopoulos | Optimistic mirror descent in saddle-point problems: Going the extra (gradient) mile[END_REF][START_REF] Yadav | Stabilizing adversarial nets with prediction methods[END_REF] and references therein. 3Going beyond the Lipschitz regime, Bauschke et al. [START_REF] Bauschke | A descent lemma beyond Lipschitz gradient continuity: First-order methods revisited and applications[END_REF] recently introduced a "Lipschitz-like" smoothness condition for convex minimization problems and used it to establish a O(1/T ) value convergence rate for mirror descent methods (as opposed to mirror-prox). Always in the context of loss minimization problems, Bolte et al. [START_REF] Bolte | First order methods beyond convexity and Lipschitz gradient continuity with applications to quadratic inverse problems[END_REF] subsequently extended the results of Bauschke et al. [START_REF] Bauschke | A descent lemma beyond Lipschitz gradient continuity: First-order methods revisited and applications[END_REF] to unconstrained non-convex problems that satisfy the Kurdyka-Łojasiewicz (KL) inequality, while Lu et al. [START_REF] Lu | Relatively-smooth convex optimization by first-order methods and applications[END_REF] considered functions that are also relatively strongly convex and showed that mirror descent achieves a geometric convergence rate in this context. Finally, in a very recent preprint, Hanzely et al. [START_REF] Hanzely | Accelerated Bregman proximal gradient methods for relatively smooth convex optimization[END_REF] examined the rate of convergence of an accelerated variant of mirror descent under the same Lipschitz-like smoothness assumption.

The condition of Bauschke et al. [START_REF] Bauschke | A descent lemma beyond Lipschitz gradient continuity: First-order methods revisited and applications[END_REF] is remarkably simple as it only posits that the problem's loss function f is such that βh -f is convex for some reference Bregman function h and some β > 0. A straightforward extension of this condition to an operator setting would be to require the monotonicity of β∇h -A, where A is the operator defining the variational inequality under study. However, the cornerstone of this "Lipschitz-like" condition is a descent lemma which does not carry over to variational inequalities, so it does not seem possible to extend the analysis of Bauschke et al. [START_REF] Bauschke | A descent lemma beyond Lipschitz gradient continuity: First-order methods revisited and applications[END_REF] to an operator setting. Lu [START_REF] Lu | Relative-continuity" for non-Lipschitz non-smooth convex optimization using stochastic (or deterministic) mirror descent[END_REF] also considered a "relative continuity" condition for loss minimization problems positing that ∇f (x) ≤ M inf x 2D(x , x)/ x -x (where f is the problem's objective and D is the Bregman divergence of h). Written this way, the condition of Lu [START_REF] Lu | Relative-continuity" for non-Lipschitz non-smooth convex optimization using stochastic (or deterministic) mirror descent[END_REF] can also be extended to an operator setting, but this would provide a surrogate for operator boundedness, not Lipschitz continuity (since A = ∇f in minimization problems). Since the optimal O(1/T ) convergence rate of the mirror-prox algorithm is tied to the regularity of A -as opposed to its boundedness -the condition of Lu [START_REF] Lu | Relative-continuity" for non-Lipschitz non-smooth convex optimization using stochastic (or deterministic) mirror descent[END_REF] does not seem applicable to the setting under study. Accordingly, there is no overlap in results or methodology with this particular strand of the literature.

Finally, in a very recent paper, Bach and Levy [START_REF] Bach | A universal algorithm for variational inequalities adaptive to smoothness and noise[END_REF] introduced a universal variant of the mirror-prox algorithm which is model-agnostic and achieves an optimal convergence rate in stochastic and/or smooth settings. Achieving optimal rates in the setting of Bach and Levy [START_REF] Bach | A universal algorithm for variational inequalities adaptive to smoothness and noise[END_REF] relies crucially on the operator being Lipschitz continuous (albeit with a possibly unknown constant) and the feasible region having a finite Bregman diameter. The algorithm we propose in this work is not universal but it is adaptive, and it does not require either Lipschitz continuity or a finite Bregman diameter. In this manner, our work also provides an important first step towards extending the universal analysis of Bach and Levy [START_REF] Bach | A universal algorithm for variational inequalities adaptive to smoothness and noise[END_REF] to VI problems with singularities.

Preliminaries

Let X be a convex -but not necessarily closed or compact -subset of a d-dimensional normed space V, and let V * denote the dual space of V. The variational inequality (VI) problem associated to a continuous operator A : X → V * consists of finding x * ∈ X such that A(x * ), x -x * ≥ 0 for all x ∈ X .

(VI)

Following [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF], we will refer to this problem as VI(X , A) and we will write X * ≡ Sol(X , A) for its set of solutions. Note also that, if X is not closed, A may exhibit a singularity at a residual point x ∈ bd(X ) \ X in the sense that A does not admit a continuous extension to x.

In the literature, this formulation of the problem is often reffered to as a Stampacchia variational inequality (SVI) [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF] or a "strong" variational inequality [START_REF] Juditsky | Solving variational inequalities with stochastic mirror-prox algorithm[END_REF][START_REF] Nesterov | Dual extrapolation and its applications to solving variational inequalities and related problems[END_REF]. For illustration purposes, we present some archetypal examples of such problems below: Example 2.1 (Loss minimization). If A = ∇f for some convex loss function f on X = R d , solutions of (VI) coincide with the global minimizers of f .

Example 2.2 (Min-max optimization). Suppose that

A = (∇ x1 f, -∇ x2 f ) for some real-valued function f (x 1 , x 2 ) with x 1 ∈ X 1 , x 2 ∈ X 2 . If f is convex-concave (i.e., convex in x 1 and concave in x 2 ), any solution x * = (x * 1 , x * 2 ) of (VI) is a global saddle-point of f , i.e., f (x * 1 , x * 2 ) ≤ f (x 1 , x * 2 ) and f (x * 1 , x * 2 ) ≥ f (x * 1 , x 2 ) (2.1) for all x 1 ∈ X 1 , x 2 ∈ X 2 .
Problems of this type have attracted considerable interest in the fields of machine learning and artificial intelligence because they constitute the basic optimization framework for GANs [START_REF] Goodfellow | Generative adversarial nets[END_REF]. For a series of recent papers focusing on the interplay between GAN and saddle-point problems / variational inequalities, see [START_REF] Daskalakis | Training GANs with optimism[END_REF][START_REF] Gidel | A variational inequality perspective on generative adversarial networks[END_REF][START_REF] Liang | Interaction matters: A note on non-asymptotic local convergence of generative adversarial networks[END_REF][START_REF] Mertikopoulos | Optimistic mirror descent in saddle-point problems: Going the extra (gradient) mile[END_REF][START_REF] Yadav | Stabilizing adversarial nets with prediction methods[END_REF] and references therein.

Example 2.3 (Resource sharing problems). Consider a set of resources r ∈ R = {1, . . . , R} serving a stream of demands that arrive at a rate of ρ per unit of time (for instance, a GPU cluster or a computing grid processing a stream of jobs). If the load on the r-th resource is x r , the expected service time in the standard Kleinrock model [START_REF] Kleinrock | Queueing Systems[END_REF] is given by the M/M/1 loss function

r (x r ) = 1 c r -x r , (2.2)
where c r denotes the capacity of the resource. In this setting, the set of feasible resource allocations is 4 and we say that a resource allocation profile x * ∈ X * is at Nash/Wardrop equilibrium [START_REF] Nisan | Algorithmic Game Theory[END_REF][START_REF] Wardrop | Some theoretical aspects of road traffic research[END_REF] if

X ≡ {(x 1 , . . . , x R ) : 0 ≤ x r < c r , x 1 + • • • + x R = ρ},
r (x * r ) ≤ r (x r ) for all x ∈ X and all r ∈ R such that x * r > 0 (2.
3) i.e., when no job would be better served by transferring it to a different priority queue. In this case, if we let A(x) = ( 1 (x 1 ), . . . , R (x R )), a standard calculation shows that x * is an equilibrium allocation if and only if it solves the associated variational inequality problem for A.

The most widely used assumption in the literature for solving VI problems is monotonicity, i.e., A(x [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]; finally, by direct calculation, it is straightforward to see that the operator defined in Example 2.3 is monotone. For an introduction to the theory of monotone operators, we refer the reader to Facchinei and Pang [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF] and Bauschke and Combettes [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]. Now, drawing on Nesterov [START_REF] Nesterov | Dual extrapolation and its applications to solving variational inequalities and related problems[END_REF][START_REF] Nesterov | Primal-dual subgradient methods for convex problems[END_REF] and Juditsky et al. [START_REF] Juditsky | Solving variational inequalities with stochastic mirror-prox algorithm[END_REF], if A is monotone, the quality of a candidate solution x ∈ X can be assessed via the restricted gap (or merit) function

) -A(x), x -x ≥ 0 for all x, x ∈ X . (2.4) When A = ∇f , this condition is equivalent to f being convex; likewise, when A = (∇ x1 f, -∇ x2 f ) as in Example 2.2, monotonicity is equivalent to f being convex-concave
Gap C (x) = sup x∈C A(x), x -x , (2.5a)
where C is a nonempty convex subset of X . The rationale behind this definition is that, if x * solves (VI), monotonicity gives A(x), x * -x ≤ A(x * ), x * -x ≤ 0, so the quantity being maximized in (2.5) is small if x is an approximate solution of (VI). Formally, we have:

Lemma 1. Suppose that A is monotone. If x * solves (VI), we have Gap C (x * ) = 0 whenever x * ∈ C. Conversely, if Gap C (x) = 0 and C contains a neighborhood of x in X , x is a solution of (VI).
This lemma extends a similar result by Nesterov [START_REF] Nesterov | Dual extrapolation and its applications to solving variational inequalities and related problems[END_REF], so we defer its proof to the paper's supplement.

In view of all this, we will employ the gap function Gap C (x) as our main figure of merit and we will use it to state our convergence rate guarantees in the sequel.

Bregman continuity

In addition to monotonicity, a standard assumption for solving variational inequalities is Lipschitz continuity, i.e., A(x ) -A(x) * ≤ β x -x (Lip) for some β > 0 and for all x, x ∈ X . This definition involves two distinct (but related) measures of distance: (i) the primal norm on V which measures distances between the primal points x, x ∈ X ; and (ii) the dual norm on V * which measures the distance between the dual vectors A(x), A(x ) ∈ V * . 5Importantly, both of these notions are global, i.e., they do not depend on the point in space at which they are calculated; as such, Lipschitz continuity is oblivious to the geometry of X (and, in particular, its boundary). In the sequel, we describe a way to overcome this limitation by introducing two distinct notions of distance that are tailored to the geometry of X and the singularity landscape of A.

Local norms. The first measure of distance that we define is that of local norm on X : Definition 1. Let Z = span(X -X ) denote the tangent hull of X , i.e., the subspace of V spanned by all possible displacement vectors of the form z = x -x, x, x ∈ X . A local norm on X is a continuous assignment of a norm • x on Z at each x ∈ X . 6 The induced dual local norm is then defined as

v x, * = max z∈Z {| v, z | : z x ≤ 1} for all v ∈ V * . (3.1)
For ease of presentation, we tacitly assume in what follows that z x ≥ µ z for some µ > 0 and all x ∈ X , z ∈ Z. This can always be achieved by taking • x ← • x + µ • so there is no loss of generality. Note in particular that this implies that v x, * ≤ (1/µ) v for all x ∈ X and all v ∈ Z * .

For intuition, we present some key examples below:

Example 3.1 (Euclidean geometry). Let X = R d so Z = R d .
The Euclidean norm on X is given by the standard expression z 2 2 = d j=1 z 2 j , and the associated dual norm is the same.

Example 3.2 (Shahshahani p-norm). Let X = R d ++ so, again, Z = R d .
The Shahshahani p-norm on X is defined for all p > 1 as

z x = |z 1 | p /x 1 + • • • + |z d | p /x d 1/p for all x ∈ X , z ∈ Z. (3.2)
By a straightforward application of Hölder's inequality, the corresponding dual norm is given by

v x, * = x q-1 1 |v 1 | q + • • • + x q-1 d |v d | q 1/q (3.3)
with the usual convention p -1 + q -1 = 1. In particular, for p → 1 + , we get the limiting expression

v x, * = max{x 1 |v 1 |, . . . , x d |v d |}. (3.4)
This metric plays a major role in, among others, game theory, optimal transport, machine learning, information theory, and many other fields -see e.g., [START_REF] Absil | Optimization Algorithms on Matrix Manifolds[END_REF][START_REF] Alvarez | Hessian Riemannian gradient flows in convex programming[END_REF][START_REF] Kakade | Regularization techniques for learning with matrices[END_REF][START_REF] Mertikopoulos | Riemannian game dynamics[END_REF][START_REF] Mertikopoulos | Distributed stochastic optimization via matrix exponential learning[END_REF][START_REF] Shahshahani | A New Mathematical Framework for the Study of Linkage and Selection[END_REF][START_REF] Sra | Optimization for Machine Learning[END_REF]] and references therein.

3. h is strongly convex relative to the underlying local norm, i.e.,

h(p) ≥ h(x) + ∇h(x), p -x + 1 2 K p -x 2 x (3.5)
for some K > 0 and all p ∈ X , x ∈ X • .

The Bregman divergence induced by h is then defined for all p ∈ X , x ∈ X • as

D(p, x) = h(p) -h(x) -∇h(x), p -x . (3.6)
As an immediate consequence of the above, we have:

Lemma 2. A Bregman function h is K-strongly convex relative to • x if and only if D(p, x) ≥ 1 2 K p -x 2
x for all p ∈ X and all x ∈ X • .

(3.7)

The main difference between Definition 2 and the standard assumptions in the literature [START_REF] Bauschke | A descent lemma beyond Lipschitz gradient continuity: First-order methods revisited and applications[END_REF][START_REF] Juditsky | Solving variational inequalities with stochastic mirror-prox algorithm[END_REF][START_REF] Mertikopoulos | Learning in games via reinforcement and regularization[END_REF][START_REF] Mertikopoulos | On the convergence of gradient-like flows with noisy gradient input[END_REF][START_REF] Mertikopoulos | Learning in games with continuous action sets and unknown payoff functions[END_REF][START_REF] Nemirovski | Robust stochastic approximation approach to stochastic programming[END_REF][START_REF] Nesterov | Dual extrapolation and its applications to solving variational inequalities and related problems[END_REF][START_REF] Nesterov | Primal-dual subgradient methods for convex problems[END_REF] is the strong convexity requirement relative to the local norm • x (whose choice, in turn, is aimed to capture the singularity landscape of the operator). We illustrate this with two examples below:

Example 3.3. Suppose that X = R d is
2 x = d i=1 |z| 2 i /(1 -x i ) 2 for x ∈ X , z ∈ R d (cf. Example 3.2 above). If we set h(x) = d i=1 1/(1 -x i ) (3.8)
a straightforward calculation gives

D(p, x) = d i=1 (p i -x i ) 2 (1 -p i )(1 -x i ) 2 ≥ d i=1 (p i -x i ) 2 (1 -x i ) 2 = p -x 2 x , (3.9)
i.e., h is strongly convex relative to • x . Importantly, since • x ≥ • 2 , this Bregman function is also strongly convex relative to the standard Euclidean norm. However, even though the Euclidean regularizer of Example 3.3 is strongly convex relative to any global norm on X , it cannot be strongly convex relative to the local norm • x because of the singularity of the latter when x i → 1 -.

Bregman continuity. We are now in a position to introduce the notion of Bregman continuity: Definition 3. Let h be a local Bregman function relative to some local norm • x on X . We say that the operator A :

X → V * is β-Bregman continuous if A(x ) -A(x) x, * ≤ β 2D(x, x ) for all x, x ∈ X . (
A(x ) -A(x) 2 x, * = d i=1 (x i -x i ) 2 (1 -x i ) 2 ≤ d i=1 (x i -x i ) 2 (1 -x i )(1 -x i ) 2 = D(x, x ) (3.10)
i.e., A is (1/ √ 2)-Bregman continuous relative to h. However, given the singularity of A(x) as x i → 1 -, we see that A cannot be Lipschitz continuous relative to any global norm on X .

Importantly, this example suggests the following rule of thumb: if the Jacobian of A exhibits a singularity of the form O(φ(x)) near the residual set cl(X ) \ X of X , taking • x = Θ(φ(x)) and h(x) = Θ(φ(x)) allows A to be Bregman continuous, despite this singularity. This heuristic provides a principled choice of Bregman data under which A satisfies (BC).

The mirror-prox algorithm

In this section, we present the main algorithmic method that we will use to solve (VI) under Bregman continuity. Our core assumptions in that regard will be: Assumption 1. The solution set X * ≡ Sol(X , A) of (VI) is nonempty. Assumption 2. A is monotone and β-Bregman continuous.

In addition to the above, we assume that the optimizer gains access to A via an oracle which, when called at the t-th stage of a sequence X t ∈ X , returns (possibly imperfect) feedback of the form

V t = A(X t ) + U t , (4.1)
where U t ∈ V * is an additive noise variable. The two cases of interest that we consider here are (i) when U t = 0 for all t; and (ii) when U t satisfies the statistical hypotheses:

a) Zero-mean: E[U t | F t ] = 0. (4.2a) b) Finite variance: E[ U t 2 * | F t ] ≤ σ 2 .
(4.2b) with F t denoting the history (natural filtration) of X t . For obvious reasons, we will refer to the first case (U t = 0) as a perfect oracle, and to the second one as a stochastic oracle.

Following Nemirovski [START_REF] Nemirovski | Prox-method with rate of convergence O(1/t) for variational inequalities with Lipschitz continuous monotone operators and smooth convex-concave saddle point problems[END_REF] and Juditsky et al. [START_REF] Juditsky | Solving variational inequalities with stochastic mirror-prox algorithm[END_REF], the mirror-prox (MP) algorithm can be stated in recursive form as follows:

X t+1/2 = P Xt (-γ t V t ) X t+1 = P Xt (-γ t V t+1/2 ) (MP)
where γ t > 0 is a variable step-size sequence (discussed in detail below), and the so-called "proxmapping" P : X • × V * → X is defined as

P x (y) = arg min x ∈X { y, x -x + D(x , x)} (4.3)
with D(•, •) denoting the divergence of an underlying Bregman function h : X → R. For concreteness, we also assume in what follows that (MP) is initialized at the so-called "prox-center" of X , i.e.,

X 1 = x c ≡ arg min x∈X h(X ). (4.4)
Remark 1. In general, calculating mirror steps can be computationally expensive -just like Euclidean projections in several cases. In what follows, we tacitly assume that our setting is "prox-friendly" [START_REF] Juditsky | Solving variational inequalities with stochastic mirror-prox algorithm[END_REF][START_REF] Nemirovski | Prox-method with rate of convergence O(1/t) for variational inequalities with Lipschitz continuous monotone operators and smooth convex-concave saddle point problems[END_REF][START_REF] Nesterov | Dual extrapolation and its applications to solving variational inequalities and related problems[END_REF] in the sense that the update (4.3) can be computed efficiently. For instance, it is straightforward to check that Example 3.4 is prox-friendly.

Heuristically, the main idea behind (MP) is that, at each t = 1, 2, . . . , the oracle is called at the algorithm's base state X t to generate an intermediate, leading state X t+1/2 ; subsequently, the base state is updated with oracle information from the leading state X t+1/2 and the process repeats. In this way, (MP) essentially tries to "anticipate" the change of A along a prox-step, and to exploit this "forward" information in order to achieve a faster convergence rate than ordinary forwardbackward/gradient descent schemes. For this anticipatory scheme to work, the variation of the operator A must be sufficiently gradual, hence the need for Lipschitz continuity in the classical analysis of the algorithm [START_REF] Juditsky | Solving variational inequalities with stochastic mirror-prox algorithm[END_REF][START_REF] Nemirovski | Prox-method with rate of convergence O(1/t) for variational inequalities with Lipschitz continuous monotone operators and smooth convex-concave saddle point problems[END_REF][START_REF] Nesterov | Dual extrapolation and its applications to solving variational inequalities and related problems[END_REF]. If this variation is unbounded (e.g., if A exhibits singularities), this look-ahead mechanism could break down completely and the algorithm might fail to converge altogether. Our first result below is that, despite such singularities, Bregman continuity allows us to recover the optimal convergence rate of (MP):

Theorem 1. Assume that A satisfies Assumptions 1 and 2, and let Gap H denote the restricted gap function for the Bregman zone C H = {x ∈ X : D(x, x c ) ≤ H}. Suppose further that (MP) is run with a K-strongly convex Bregman function and oracle feedback of the form (4.1). Then, for all H > 0, the averaged sequence XT = T t=1 γ t X t+1/2 T t=1 γ t enjoys the following gap bounds: a) If σ 2 = 0 and the algorithm's step-size satisfies

0 < γ min ≡ inf t γ t ≤ sup t γ t ≡ γ max ≤ √ K/β, (4.5)
we have

Gap H ( XT ) ≤ H γ min 1 T (4.6) b) Otherwise, if σ 2 > 0 and γ t ≤ K/2/β, we have E[Gap H ( XT )] = O H+σ 2 T t=1 γ 2 t T t=1 γt (4.7) In particular, if γ t ∝ 1/ √ T , we get E[Gap H ( Xt )] = O(1/ √ T ).
Remark. We should stress here that the above rates are attained even if X is unbounded and/or the "Bregman content"

H X ≡ sup x D(x, x c ) = sup h -inf h of X is infinite.
As we show in the supplement, the key step in the proof of the deterministic part of Theorem 1 is the following energy inequality for an arbitrary target point p ∈ C H :

D(p, X t+1 ) ≤ D(p, X t ) -γ t A(X t+1/2 ), X t+1/2 -p -1 - β 2 γ 2 t K D(X t+1/2 , X t ) (4.8)
There are two points where the Bregman structure of the algorithm can be seen in ( 4 D(X t+1/2 , X t ). In the "vanilla" setting, Lipschitz continuity is used to obtain a comparison of these successive states in terms of a global norm difference of the form X t+1/2 -X t 2 . However, this step also requires A to vary gradually relative to • , which is of course impossible if A exhibits singularities. The key novelty in our setting is the use of the Bregman divergence as a comparator for the algorithm's successive states: it is at this point that the triple interplay between the operator, the local norm and the chosen Bregman function is made manifest, and it is what makes Bregman continuity particularly well-suited for tackling singular problems of this kind. This requires a careful treatment of the various Bregman differences involved, so we defer the details to the supplement.

The adaptive mirror-prox algorithm

A crucial assumption underlying the analysis of the previous section is that the optimizer must know in advance -or be otherwise able to estimate -the Bregman constant β. In practice, this can be difficult to achieve, so it is important to be able to run (MP) with an adaptive step-size policy. Our starting point is the observation that, with accurate feedback, one can estimate β by setting

β t = A(X t+1/2 ) -A(X t ) X t+1/2 , * 2D(X t+1/2 , X t ) (5.1)
whenever X t+1/2 = X t ; obviously, if A is β-Bregman continuous, we have β t ≤ β. However, the fact that the Bregman constant is being under-estimated means that a step-size policy of the form γ t ∝ √ K/β t would over-estimate the inverse Bregman constant 1/β, so the resulting step-size policy would have no reason to satisfy (4.5). To overcome this obstacle, we introduce the following comparison mechanism: first, at each t = 1, 2, . . . , we use the estimation (5.1) to test the step-size γt = √ K/β t . Then, to avoid the growth phenomenon outlined above, we shrink γt by a constant factor of θ and, to avoid running into vanishing step-size issues, we take the previous step-size employed if the shrunk one would be smaller. Formally, we consider the adaptive step-size policy:

γ t+1 = min{γ t , θ √ K/β t } if X t = X t+1/2 , γ t otherwise, (5.2)
with β t defined as in (5.1) and θ ∈ (0, 1) chosen arbitrarily.

For concreteness, we call the resulting algorithm adaptive mirror-prox (AMP) and we provide a pseudocode implementation in Algorithm 1 above. In terms of performance, we have: Theorem 2. Assume that A satisfies Assumptions 1 and 2, and (MP) is run with perfect oracle feedback and the adaptive step-size policy (5.2). Then, with notation as in Theorem 1, the algorithm's ergodic average XT =

T t=1 γ t X t+1/2 T t=1 γ t enjoys the gap bound Gap H ( XT ) = O(1/T ).
We find this result particularly appealing because it yields the optimal O(1/T ) convergence rate of the mirror-prox algorithm, even for possibly singular operators, and even if the operator's Bregman constant is unknown. Its proof relies on using the specific form of the step-size policy (5.2) to control the second term in the energy inequality (4.8); we provide the detailed arguments in the supplement.

○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ ○ �����-�������� (

Numerical experiments

We performed a series of numerical experiments on the resource sharing problem described in Example 2.3 with a set of R = 1000 servers being shared by N = 100 commodities, each with a demand drawn uniformly at random from [0, 1]; the capacity c r of each server r = 1, . . . , R was also drawn randomly from [0, 100]. Subsequently, we ran two variants of the mirror-prox method: (MP) with Euclidean regularization, and (MP) with the Bregman function defined in Example 3.4. For all methods, we ran a range of different constant step-sizes (we present the most representative values, namely γ = 0.001, γ = 0.005, and γ = 0.010). Subsequently, we also ran Algorithm 1 and we plotted the distance from the solution to the induced variational inequality problem as a function of the number of iterations. The main conclusions that can be drawn are as follows:

1. The Euclidean version of the mirror-prox algorithm (i.e., the extra-gradient algorithm) is unstable and does not converge; this is due to the fact that the gradients received are very large (recall that the problem is not Lipschitz continuous), so the algorithm does not exhibit descent or convergence. 2. The MP variant with the non-Euclidean regularizer of Example 3.4 behaves much better and converges (since the VI problem under study is Bregman continuous relative to this Bregman function). However, depending on the method's step-size, the convergence is relatively slow, and there is no easy way to estimate the problem's Bregman constant in order to choose a "good" step-size. 3. By contrast, the adaptive version of the method (Algorithm 1) converges significantly faster than variants with a constant step-size. Importantly, the method's step-size changes significantly from one iteration to the next, at least for the first few thousand iterations, eventually converging to a value close to γ inf ≈ 0.0057. However, running the method with this constant step-size is significantly slower in the warm-up period of the method, and leads to worse results overall. This is due to the fact that, initially, a greedier step-size is able to take larger steps towards the problem's solution.

Concluding remarks

In this work, we introduced a novel regularity condition to account for variational inequalities (both deterministic and stochastic) with possible singularities. This condition, which we call Bregman continuity, is tailored to the operator's singularity landscape and, as such, provides the necessary bedrock to achieve optimal convergence rates via a properly chosen version of the mirror-prox algorithm (with or without knowledge of the problem's Bregman constant). This opens up several interesting research directions: First, an appealing extension would be to develop a "model-agnostic" version of the method (which would concurrently provide optimal rates in stochastic and deterministic settings) or to combine it with backtracking / linesearch to accelerate convergence. Finally, it would also be interesting to examine the method's local convergence properties in non-monotone problems (deterministic or stochastic). We relegate these questions to future work.

B Bregman functions and Bregman continuity B.1 Properties of Bregman functions

In this appendix, we present some basic facts about Bregman functions and proximal mappings. Similar results exist in the literature in different contexts (see e.g., [START_REF] Juditsky | Solving variational inequalities with stochastic mirror-prox algorithm[END_REF][START_REF] Nesterov | Dual extrapolation and its applications to solving variational inequalities and related problems[END_REF][START_REF] Nesterov | Primal-dual subgradient methods for convex problems[END_REF] and references therein), but given that many of our results rely on the use of local -as opposed to global -norms, we provide here complete statements and proofs.

To begin, we introduce two notions that will be particularly useful in the sequel. The first is the convex conjugate of a Bregman function h, i.e., h * (y) = max

x∈V { y, x -h(x)} (B.1)
and the associated primal-dual "mirror map" Q : V * → X :

Q(y) = arg max x∈V { y, x -h(x)} (B.2)
We then have the following basic lemma connecting the above notions: Lemma B.1. Let h be a K-strongly convex Bregman function as above. Then, for all x ∈ X • ≡ dom ∂h, y ∈ V * we have: (B.4) Since, h is strongly convex and y ∈ ∂h(x) due to the first equivalence, it follows that φ(t) ≥ 0 with equality if and only if t = 0. Since, ψ(t) = ∇h(x + t(p -x)) -y, p -x is a continuous selection of subgradients of φ and both φ and ψ are continuous over [0, 1], it follows that φ is continuously differentiable with φ = ψ on [0, 1]. Hence, with φ convex and φ(t) ≥ 0 = φ(0) for all t ∈ [0, 1], we conclude that φ (0) = ∇h(x) -y, p -x ≥ 0 and thus we obtain the result.

1. x = Q(y) ⇐⇒ y ∈ ∂h(x) 2. x + = P x (y) ⇐⇒ ∇h(x) + y ∈ ∂h(x + ) ⇐⇒ x + = Q(∇(x) + y) 3. Finally, if x = Q(y)
The basic ingredient for establishing connections in the Bregman framework is a generalization of the rule of cosines which is known in the literature as the "three-point identity" [START_REF] Chen | Convergence analysis of a proximal-like minimization algorithm using Bregman functions[END_REF] and will be the main tool for deriving the main estimations for our analysis. Being more precise, we have the following lemma: Lemma B.2. Let h be a Bregman function on X . Then, for all p ∈ X and all x, x ∈ X • , we have: D(p, x ) = D(p, x) + D(x, x ) + ∇h(x ) -∇h(x), x -p (B.5)

The proof of this lemma follows as in the classic Bregman case [START_REF] Chen | Convergence analysis of a proximal-like minimization algorithm using Bregman functions[END_REF] so we omit it and proceed to derive some key bounds for the Bregman divergence before and after a mirror step: Proposition B.1. Let h be a local Bregman function with strong convexity modulus K > 0. Fix some p ∈ X and let x + = P x (y) for some x ∈ X • and y ∈ V * . We then have: with Z 1 = x c . We then have -γ t ξ t+1 = γ t U t+ 1 2 , p -X t+

C.3 Adaptive analysis

Our aim here is to prove the O(1/T ) convergence rate of the adaptive mirror-prox algorithm:

Proof of Theorem 2. For simplicity, we will assume in what follows that X t+ 1 2 = X t for all t ≥ 1. Otherwise, if X t+ 1 2 = X t for some t, we would also have γ t+1 = γ t for said value of t by convention; this would not change our arguments below, but it would make them much more cumbersome to write down.

With this caveat in mind, we begin by showing that the adaptive step-size policy γ t+1 = min{γ t , θ √ K/β t } is bounded from below as

γ t ≥ min{γ 1 , θ √ K/β}. (C.23)

  BC)Of course, in the case of a global norm with Bregman function h(x) = (1/2) x 2 (cf. Example 3.3), we recover the standard Lipschitz continuity condition:A(x ) -A(x) * ≤ β 2D(x , x) = β x -x .On the other hand, the example below shows that an operator can be Bregman continuous without being Lipschitz continuous relative to any global norm: Example 3.5. Consider the operator A(x) = (c r • (1 -x r /c r ) -1 ) r∈R defined in Example 2.3. Renormalizing c r to 1 for clarity and using the Bregman data of Examples 3.2 and 3.4, we get:

Figure 1 :

 1 Figure 1: Different variants of the mirror-prox algorithm in the resource sharing problem of Example 2.3. The algorithm labeled "extra-gradient" refers to Euclidean regularization and a constant step size as indicated in the legend; "mirror-prox" was run with the Bregman function of Example 3.4 and step-sizes as in the legend; finally, "adaptive mirror-prox" corresponds to Algorithm 1, i.e., mirror-prox with the adaptive step-size (5.2).

D

  (p, x + ) ≤ D(p, x) -D(x + , x) + y, x + -p (B.6) Proof. By the three-point identity established in Lemma B.2, we get: D(p, x) = D(p, x + ) + D(x + , x) + ∇h(x) -∇h(x + ), x + -p (B.7) By rearranging the terms we get: D(p, x + ) = D(p, x) -D(x + , x) + ∇h(x + ) -∇h(x), x + -p (B.8) Due to (B.3) and the fact that x + = P x (y) so ∇h(x) + y ∈ ∂h(x + ), we get the result.

  endowed with the Euclidean norm as in Example 3.1. Then, setting h(x) = (1/2) x 2 2 , we get the standard expression D(p, x) = (1/2) p -x 2 2 for the associated Bregman divergence. Obviously, h is 1-strongly convex relative to • 2 . Example 3.4. Let X = [0, 1) d (so X is neither open nor closed), and consider the local norm z

1 2

 1 = γ t U t+ 1 2 , Z t -X t+ 1 2 + γ t U t+ 1 2 , p -Z t (C.18)In order to bound the term which depends on p, we have the following:γ t U t+ 1 2 , p -Z t = γ t U t+ 1 2 , p -Z t+1 + γ t U t+ 1 2 , Z t+1 -Z t ≤ µ ∇h(Z t+1 ) -∇h(Z t ), p -Z t+1 + γ 2Hence, by the three-point identity, we obtain:γ t U t+ 1 2 , p -Z t ≤ µ[D(p, Z t ) -D(p, Z t+1 )] -µD(Z t+1 , Z t )where the last inequality is a consequence of the strong convexity of h. Thus, combining all these with the fact that A is monotone, we can telescope and obtainU t+ 1 2 , Z t -X t+ 1 2 ] = E[E[ U t+ 1 2 , Z t -X t+ 1 2 ] | F t+12 ] = 0, taking expectations yields Finally, the RHS of this last inequality is Õ(1/T 1/2 ) if γ t ∝ 1/

								t 2K	U t+ 1 2	2 * ,Zt +	K 2	Z t+1 -Z t	2 Zt
	≤ µ ∇h(Z t+1 ) -∇h(Z t ), p -Z t+1 +	γ 2 t 2µK	U t+ 1 2	2 * +	Kµ 2	Z t+1 -Z t	2 .
									(C.19)
	+	γ 2 t 2µK	U t+ 1 2	2 * +	Kµ 2	Z t+1 -Z t	2
	≤ µ[D(p, Z t ) -D(p, Z t+1 )] +	γ 2 t 2µK	U t+ 1 2	2 *	(C.20)
				T
				γ t U t+ 1 2 , Z t -X t+ 1 2
				t=1
		+	1 µK	T t=1	γ 2 t 2 U t	2 * +	5 2	U t+ 1 2	2 * .
	Hence, after dividing by
	T t=1 γ t U t+ 1 2 , Z t -X t+ 1 2 T t=1 γ t + 1 µK	T t=1 γ 2 t 2 U t	2 * + 5 2 U t+ 1 2	2 *	.
									(C.21)
	Since E[ E[Gap H ( Xt )] ≤ which proves our claim. √ (1 + µ)D + 9σ 2 2µK T t=1 γ 2 t T , t=1 γ t Õ(1/ √ T ) result follows.	(C.22) t, so the

T t=1 γ t A(p), X t+ 1 2 -p ≤ (1 + µ)D(p, x c ) + T t=1

γ t and taking the supremum over p ∈ C H , we get:

Gap H ( XT ) ≤ (1 + µ)H +

33rd Conference on Neural Information Processing Systems (NeurIPS 2019), Vancouver, Canada.

When used to find a zero of a composite operator, the FB algorithm is known as a "splitting" method; see e.g., Bruck Jr.[START_REF] Bruck | On the weak convergence of an ergodic iteration for the solution of variational inequalities for monotone operators in Hilbert space[END_REF], Passty[START_REF] Passty | Ergodic convergence to a zero of the sum of monotone operators in Hilbert space[END_REF],[START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF], and references therein.

An operator A(x) is cocoercive if A(x ) -A(x), x -x ≥ (1/β) A(x ) -A(x)[START_REF] Alvarez | Hessian Riemannian gradient flows in convex programming[END_REF] for some β > 0 and all x, x . Note that Lipschitz continuity is strictly weaker than cocoercivity: the operator A(x1, x2) = (-x2, x1) is Lipschitz continuous over R 2 , but it is not cocoercive; see Section 2 for a detailed discussion.

We note here that the method is sometimes referred to as "optimistic mirror descent". This terminology is due to Rakhlin and Sridharan[43,

[START_REF] Rakhlin | Optimization, learning, and games with predictable sequences[END_REF] and may refer either to the mirror-prox method itself, or to a variant with "gradient extrapolation from the past", as in[START_REF] Gidel | A variational inequality perspective on generative adversarial networks[END_REF].

For posterity, note here that X is convex but it is not necessarily closed.

Recall here that the dual norm of v ∈ V * is defined as v * = maxz∈V {| v, z | : z ≤ 1}.

By that, we have in mind the definition of an absolutely homogeneous Finsler metric[START_REF] Bao | An Introduction to Riemann-Finsler Geometry[END_REF]. Specifically, a local norm is viewed here as continuous nonnegative function F : X × V → R+ with the following propoerties: for all x ∈ X and all z1, z2 ∈ V, we have (i) F (x, z1 + z2) ≤ F (x, z1) + F (x, z2); (ii) F (x, λz) = |λ|z; and (iii) F (x, z) > 0 for all z ∈ V \ {0}. The local norm of z at x is then defined as z x = F (x, z).

Acknowledgments

The authors gratefully acknowledge financial support from the French National Research Agency (ANR) under grant ORACLESS (ANR-16-CE33-0004-01) and the COST Action CA16228 "European Network for Game Theory" (GAMENET).

A Properties of the restricted gap function

In this appendix, we discuss the basic properites of the restricted merit function Gap C introduced in the main text. For completeness, we begin with the proof of Lemma 1, itself an extension of a similar result by [START_REF] Nesterov | Dual extrapolation and its applications to solving variational inequalities and related problems[END_REF]:

Proof of Lemma 1. Let x * ∈ X be a solution of (VI) so A(x * ), x -x * ≥ 0 for all x ∈ X . Then, by monotonicity, we get: A(x), x * -x ≤ A(x) -A(x * ), x * -x + A(x * ), x * -x = -A(x * ) -A(x), x * -x -A(x * ), x -x * ≤ 0, (A.1) so Gap C (x * ) ≤ 0. On the other hand, if x * ∈ C, we also get Gap(x * ) ≥ A(x * ), x * -x * = 0, so we conclude that Gap C (x * ) = 0.

For the converse statement, assume that Gap C (x) = 0 for some x ∈ C and suppose that C contains a neighborhood of x in X . First, we claim that the following inequality holds:

Indeed, assume to the contrary that there exists some x 1 ∈ C such that

This would then give

which is a contradiction. Now, we further claim that x is a solution of (VI),i.e.,:

Indeed, if we suppose that there exists some z 1 ∈ X such that A(x), z 1 -x < 0, then, due to the continuity of A there exists a neighborhood U of x in X such that

Hence, assuming without loss of generality that U ⊂ U ⊂ C (the latter assumption due to the assumption that C contains a neighborhood of x), and taking λ > 0 sufficiently small so that x = x + λ(z 1 -x) ∈ U , we get that A(x), x -x = λ A(x), z 1 -x < 0, in contradiction to (A.2). We conclude that x is a solution of (VI), as claimed.

For intuition, we discuss below the relation of the error function Gap ≡ Gap X to other performance measures that arise in practice (similar considerations also apply to Gap C as well):

(1

(2) In saddle-point problems, the quality of a candidate solution x = (x 1 , x2 ) is often assessed via the Nikaido-Isoda (NI) function

provided of course that the right-hand side is well-posed. Similarly to the minimization framework above, if f is convex-concave, we have Gap(x) ≤ NI(x).

In view of all this, a bound on Gap(x) does not immediately translate to a bound on the value of the loss or the Nikaido-Isoda function (for minimization or min-max problems respectively). However, the same arguments used to obtain the convergence rate of an algorithm relative to Gap can be usually adapted to these measures with minimal effort. Now, using the above estimations, we show below the main inequalities which relate the Bregman divergence between two prox-steps: Proposition B.2. Let h be a Bregman function on X and fix some p ∈ X , x ∈ X • . Letting

x + 1 = P x (y 1 ) and x + 2 = P x (y 2 ), we have:

In particular, we have:

Proof. For the first inequality, by applying (B.1), (B.6) for x + 2 = P x (y 2 ), we get:

For the second inequality, we need to bound

In particular, applying again (B.1), (B.6) for p = x + 2 , we get:

and hence:

(B.13) So, combining the above inequalities we get:

and thus we get the second inequality as well.

B.2 Bregman cocoercivity

For comparison to the operator theory literature, we end this section by introducing a notion close to Bregman continuity, that of Bregman cocoercivity. In particular, we have the following definition: Definition 4. Let A : X → V * be an operator, let • x be a local norm on X and let h be a local Bregman function on X with strong convexity modulus K. We say that A is δ-Bregman cocoercive if

This notion is a straightforward extension of ordinary operator cocoercivity, e.g., as defined in [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF].

As in the base case, we have: Lemma B.3. Let h be a Bregman function with local strong convexity modulus K > 0 and let A : X → V * be a δ-Bregman cocoercive operator. Then, A is 1/(δ √ K)-Bregman continuous.

Proof. By applying the Cauchy-Shwartz inequality to the definition (B.15) of Bregman cocoercivity, we get:

Thus, by simplifying and recalling Lemma 2 in the main paper, we get:

i.e., A is Bregman continuous with modulus of continuity 1/(δ √ K), as claimed.

C Convergence analysis of the mirror-prox algorithm

In this section, we provide the proofs of our main convergence results for the mirror-prox algorithm as presented in Sections 4 and 5 of the main part of the paper.

C.1 Deterministic analysis

The main ingredient of our proof for the deterministic case is the following energy inequality: Proposition C.1. Assume that A satisfies Assumption 2 and (MP) is run with perfect oracle feedback.

Then, for all p ∈ X , we have:

2 ) and x + 2 = X t+1 in Proposition B.2, we readily obtain:

Proceeding line-by-line, the Fenchel-Young inequality applied to the function φ

Thus, by substituting in (C.1), we get

and hence, by Lemma 2, we obtain:

However, the Bregman continuity of A also yields

so our claim follows by combining Eqs. (C.4) and (C.5).

We are now in a position to establish the O(1/T ) convergence rate of (MP) for deterministic problems:

Proof of Theorem 1 -deterministic case. Fix some p ∈ C H . Since γ t ≤ 1/β by assumption, a slight rearrangement of Proposition C.1 readily yields:

Moreover, by the monotonicity of A, we also have:

Thus, combining the two inequalities above, we get

and, proceeding to telescope from t = 1 to T , we obtain:

Then, dividing by T t=1 γ t finally yields

so our result follows by taking the supremum over all p ∈ X such that D(p, x c ) ≤ H (i.e., over all p ∈ C H ).

C.2 Stochastic analysis

We now turn to our results for the convergence of the mirror-prox algorithm in the stochastic case:

Proof of Theorem 1 -stochastic case. Working in the same spirit as for the deterministic case, let x = X t ,

in the first part of Proposition B.2. We then get:

where we used the feedback decomposition

By the second part of Proposition B.2, we also have

Now, by applying the Fenchel-Young inequality to the duality pairing in the above inequality, we get

On the other hand, by the stochastic oracle assumption (4.1), we have:

where the last line follows from the Bregman continuity of A (Assumption 2) and the fact that • x ≥ µ • for some µ > 0 and all x ∈ X (implying in turn that • x, * ≤ µ -1 • * for all x ∈ X ). We thus get

Since γ 2 t ≤ K/(2β 2 ) by assumption, substituting (C.15) in (C.11) and rearranging yields

In order to bound ξ t+ 1 2 , we will need to introduce the auxilliary process

We consider two cases below: First, if γ 1 ≤ θ √ K/β, we will also have

for all t ≥ 1. We then claim that γ t = γ 1 for all t ≥ 1: indeed, assuming inductively that this is the case for some ≥≥ 1 (the claim is trivially true for t = 1), we readily get

which proves our claim in this case.

Assume now that γ 1 > θ √ K/β, in which case we are left to show that γ t ≥ θ √ K/β for all t ≥ 1. Again, the base case t = 1 being true trivially, assume inductively that this holds for some t ≥ 1. Then, one of the following will hold:

In both cases we obtain γ t+1 ≥ θ √ K/β, so the induction is complete.

To proceed, going back to the proof of the basic energy inequality (C.1), the intermediate step (C.4) can be rewritten as

where we used the fact that γ t+1 ≤ θ √ K/β t by construction. However, since γ t+1 is weakly decreasing and bounded from below, it follows that γ t converges to some limit value γ ∞ as t → ∞. In turn, this implies that for all t greater than some (finite) t 0 . Accordingly, rearranging (C.26) and subsequently telescoping as in (C.9), we finally obtain

whenever T > t 0 . Our result then follows by dividing both sides of this last inequality by T t=1 γ t and recalling the fact that γ t ≥ min{γ 1 , θ √ K/β} > 0 for all t.