Ernesto Garcia

Fernando Pirani

Antonio Laganà

Carles Martí

The role of the longrange tail of the potential in O 2 +N 2 collisional inelastic vibrational energy transfers

In the study of nonreactive energy transfer between the O 2 and N 2 molecules bearing different vibrationally excited states we have faced the problem of selecting a proper formulation of the in teraction. To this end we have compared the values of the related observables computed either on a potential energy surface globally fitted to a very large of ab initio potential energy values [Varga et al., J. Chem. Phys., 2016, 144, 024310] and on a more traditional one formulated as a combina tion of an intra and intermolecular model component of the interaction (both based on combined use of experimental and ab initio information) [Garcia et al., J. Phys. Chem. A, 2016, 120, 5208] in order to enforce an appropriate modelling of the longrange tail of the potential, crucial for the description of the inelastic vibrational energy transfer. The detailed graphical analysis of the po tential plus the quantitative analysis of the computed opacity functions, of the statetostate rate coefficients, of the second virial coefficient and of the integral nonreactive cross section allowed us to conclude that the model formulation of the interaction has to be preferred for nonreactive studies of the O 2 + N 2 energy transfer processes in thermal and subthermal regimes.

Introduction

In the study of the elementary molecular processes of key impor tance for gas phase simulations (like plasma chemistry, aerother modynamics, etc.) [START_REF] Park | Nonequilibrium Hypersonic Aerothermodynamics[END_REF][START_REF] Capitelli | Plasma Kinetics In Atmospheric Gases[END_REF][START_REF] Da Silva | [END_REF][4][5][6][7][8] our interest has recently focused both on the peculiar role played by nonreactive quasiresonant (the vari ation of the vibrational energy of one partner is almost exactly compensated by the opposite sign variation of the vibrational en ergy of the other partner) collisions and on the computational prediction of the experimental nonArrhenius behaviour at sub thermal conditions. [9][10][11][12][13] Due to the strong entrance channel lo calization of the associated energy transfer mechanisms of these processes, a basic ingredient of the related studies is the accurate determination of that part of the Potential Energy Surface (PES) to which dynamical calculations of the detailed statetostate rate coefficients are highly sensitive. Accordingly, before undertaking a systematic investigation of such quantities for the title system we devoted specific attention to the formulation of the related PES and to the calculation of quantities amenable to the integral I T defined as with R being the intermolecular distance and V being the value of the potential energy properly averaged (or integrated) over the angular variables and the internuclear distances. In particu lar, we considered the following quantities: the value of the to tal integral (elastic and inelastic) nonreactive cross section Q T , its glory interference pattern and the temperature dependence of the second virial coefficient B(T ) which are all observables sen sitive to the potential features via the above mentioned integral I T . [START_REF] Hirst | A Computational Approach to Chemistry[END_REF][START_REF] Stogryn | [END_REF] As is well known, in fact, and will be further discussed later, typically the long and intermediate range intermolecular in teraction of diatomic apolar molecules (as is the case of the title system especially at low temperature), is attractive and the area of its negative (lower than the asymptotic value) part contributes significantly to the value of the integral of eq. 1.

I T = Z  0 R 2 exp - V k B T -1 dR, (1
For this reason in the past we have formulated the interac tion of nonreactive diatomdiatom systems using a semiempiri cal functional form, whose parameters associated with the elec tronic polarizability and permanent multipoles are at the same time tested on the measured Q T and B(T ) values and fine tuned by means of ab initio calculations, specifically tailored to repro duce the low temperature dependence of nonreactive processes. In this way we have been able to perform extended studies of the properties of the N 2 + N 2 system (see refs. 10,16-18) and to in vestigate some aspects of the O 2 + N 2 one (see refs. 9,12,13,19) on PESs referring either to diatoms at internuclear (or intramolec ular) distances ( r) close to the equilibrium value r e (MF1) 12 or taking into consideration also their configuration distortions oc curring away from equilibrium (MF2). 13 Recently, thanks to the massive use of fast High Performance Computer (HPC) architectures integrated with High Throughput Computer (HTC) networked (distributed) infrastructures, large quantities of molecular geometry electronic structure values have become more easily obtainable. Accordingly, the alternative ap proach of determining the value of the ab initio potential energy for a large grid of molecular geometries has become popular de spite its much heavier computational cost. Then the large mass of computed ab initio potential energy values are globally fitted using a permutationally invariant polynomial in either bond or der variables [20][21][22] or polynomials in mixed exponentialgaussian variables 6,23 or using a local moving least squares interpolating functions. 24,25 This method makes internally consistent the for mulation of two, three and four body components of the interac tion and has the advantage of providing a valid description of the collision induced fragmentation of the involved molecules. 6,26 In the progress of our studies aimed at characterizing the non reactive diatomdiatom behaviour of O 2 + N 2 collisions, we first started using an older semiempirical PES (GB1) 27 , and we moved later to the more recently formulated semiempirical PESs MF1 and MF2. 12,13 In particular, Quantum Classical (QC) values com puted on MF1 12 agreed with measured vibrationto vibration (VV) and vibrationtotranslation (VT) rate coefficients 28 and with VV rate coefficients for transitions involving the lowest vibra tional states. 29 The recent publication of a new PES (called here MN) 23 based only on ab initio potential energy values allowed a comparative study of the rate coefficients computed on PESs rang ing fom empirical (GB1) to completely ab initio (MN) ones using the well established QuasiClassical Trajectory (QCT) technique. [START_REF] Andrienko | 55th AIAA Aerospace Sciences Meeting[END_REF] The comparison confirmed, however, as already mentioned in ref. 23, that the MN PES is not suited for thermal and subthermal in vestigations due to its inaccurate formulation of the longrange weak interaction.

A more quantitative analysis of the impact of such inaccuracy of the computed value of the O 2 + N 2 rate coefficients for in elastic vibrational energy transfer in conditions in which neither atom exchange nor fragmentation occurs is the central goal of the present paper and is instrumental to ground our ongoing QCT and Quantum Classical studies of thermal and subthermal pro cesses. [START_REF] Kurnosov | [END_REF] Accordingly, the paper is articulated as follows: In sec tion 2 the key differences of the nonreactive outcomes of QCT calculations performed on both the MF2 and the MN O 2 + N 2 PESs are discussed. In section 3 a detailed analysis of the features of the PESs used for the calculations with respect to some exper imental and theoretical information discriminant for the class of processes of our interest is performed. In section 4 a rational ization of the different dynamical behaviours of the system on the considered PESs is given in terms of opacity functions and vibrational statetostate rate coefficients. In section 5 some con clusions are drawn.

A comparison of MN versus MFenergy transfer QCT rate coefficients

As already mentioned, for the systematic investigation of the en ergy transfer properties of diatomdiatom nonreactive collisions we usually formulate the PES by partitioning the overall interac tion V in an intramolecular (intra) and an intermolecular ( inter) component as follows

V = V intra +V inter . ( 2
)
In eq. 2 V intra is usually described as a sum of the separated diatomic potentials while V inter is semiempirically formulated in terms of isotropic and anisotropic components of the interaction between the two separated diatoms whose parameter values are tested against experimental and ab initio data. 12 . For this pur pose each molecule i is either set at its equilibrium configuration (the intramolecular distancer i is frozen at its equilibrium valuer ei in the MF1 PES) or assumed to progressively distort (as is for the MF2 PES) while approaching the close distance configurations. This formulation of V allows a smooth evolution of the PES from a pure longrange asymptotic form into a shortrange structured one if at short intermolecular distance R (ie. when the value of R becomes comparable with that of the involved intramolecular distances r i ) the V inter formulation is modulated to account for the failure of the initial diatomdiatom model. In our approach this flexibility is gained by formulating V inter in terms of the in ternuclear distance dependence of both the molecular bond po larizability (r i ) and the electric quadrupole moment Q(r i ) eval uated using multireference ab initio techniques 9 and then fit ted to a polynomial in r i . V intra is made, instead, of the two di atomic potentials which can be formulated as Morse potentials or more in general as BondOrder (BO) polynomials, V BO (n i ), with n i = exp[- (r i -r ei )] and  being a parameter best fitted to the spectroscopic and scattering properties of the diatom. [START_REF] Garcia | [END_REF] Attempts to unify the formulation of V intra and V inter in terms of bond order variables are also being carried out. [START_REF] Kurnosov | [END_REF] On the contrary, the MN PES is based on the massive calcu lation of a dense grid of 55,000 high level ab initio points and on their fitting using polynomials in mixed exponentialgaussian variables. 23 As already mentioned, while this allows the simul taneous description of all the possible rearrangements and frag mentations of the O 2 + N 2 system, does not provide sufficient information to accurately describe the longrange interaction and enable a proper evaluation of nonreactive rate coefficients at thermal and subthermal conditions.

In order to quantify the difference in the dynamics and kinetics properties of the O 2 + N 2 computed on the MF2 and MN PESs using the QCT popular VENUS code [START_REF] Hase | VENUS96: A General Chemical Dynamics Computer Program[END_REF] we used the same initial conditions adopted in refs. 12 and 13.

Typical outcomes of the QCT calculations performed for col lisions occurring with transfer of VT energy are illustrated in Figure 1. In Figure 1 the value of the normalized to the maxi mum (absolute values of this type of quantities will be discussed later in the paper) difference between the inelastic rate coeffi cients computed on the MF2 and on the MN PES for the process

O 2 (v = 13,19,25) + N 2 (u = 0) → O 2 (v 0 ) + N 2 (u 0 = 0) at 1000
and 5000 K is plotted as a function of E v 0 -v where E v is the en ergy associated with the vibrational state v. As apparent from the figure, the largest differences in the value of the rate coef ficients computed on the two PESs occur for transitions being (v 0v) = ±1 emphasizing the importance of relying on an ac curate formulation of the entrance channel of the PES to prop erly describe the dominant nonreactive processes and under stand the related energy transfer mechanisms. Strictly speaking the comparison between the results obtained on the MF2 and the MN PESs is not on exactly corresponding processes due to the fact that V intra is formulated in a different way (a Morse curve, )) i (r -r e )]) [START_REF] Capitelli | Plasma Kinetics In Atmospheric Gases[END_REF] -D e , for MN, being D e the dissociation energy of the diatom, that rep resents the intramolecular motion in a more hard and/or har monic fashion respect to the canonical Morse) in the two PESs. This makes the vibrational energy E v associated with the same vi brational numbers not identical and this is the reason why in the figure we plot the results as a function of the corresponding vi brational energy jump E v 0 -v rather than the vibrational number difference (v 0 -v). As is appreciable by the comparison of the plots of the upper and the lower panel of Fig. 1, the difference between the MF2 and MN rate coefficients increases and extends to larger E v 0 -v as the temperature increases.

D e (1 -exp[- (r -r e )]) 2 -D e ,
In order to single out the importance of this effect for quasi resonant processes, we report in Table 1 the ratios of the rate coefficients computed with both PESs for VV O

2 (v) + N 2 (u) → O 2 (v 0 = v + 1) + N 2 (u 0 = u -1)
energy transfer transitions at the temperature of T = 750, 1500 and 5000 K, v = 15 and u = 15,25,35,45 (with (15,15|16,14), (15,25|16,24), (15,35|16,34) and (15,45|16,44) being their more compact notation). The ratio for k (15,15|16,14) at T = 750 K is not reported because it is zero on the MF2 PES. It is important to point out here that the effect is usually not smaller than one order of magnitude and increases as T lowers (e.g. for the (15,35|16,34) transition that occurs with an energy gap of a few thousandth of eV is 13 at T = 1500 K and 40 at T = 750 K) while it less predictably depends on u.

(v) + N 2 (u) → O 2 (v 0 = v +1) + N 2 (u 0 = u -1)

A detailed analysis of the intermolecular interaction

In order to single out the specific origin of the above commented differences in the rate coefficients computed on the MF2 and MN PESs, we have carried out a detailed comparison of the related entrance channels.

A first appreciation of the differences between the MN and the MF2 PESs is provided by the already mentioned formulation of V intra for which the MF2 PES adopts a Morse while for MN PES adopts a generalized Morse having a more harmonic nature (sharper repulsive walls) resulting in a larger spacing of the vibra tional levels and the consequent displacement in energy of levels corresponding to the same quantum number in Fig. 1. What marks the real difference for nonreactive scattering is, as already mentioned in the Introduction, the shape of the en trance channel longrange potential and in particular the features related to V inter that in MF2 is accurately modelled via quantum calculations to reproduce B(T ) and Q T experimental information and therefore amenable to the value of the I T integral of eq. 1. In particular, the second virial coefficient available for a wide range of temperatures (in which its value switches from negative to pos itive) and the interference oscillations modulating the velocity de pendence of the total integral nonreactive cross section provide a detailed probe of some critical features of V inter such as the depth of the well, the position of the repulsive wall and the anisotropy of the negative portion of the interaction. On the contrary, as al ready commented before, the MN PES although based on the fit ting of a large number of ab initio values is preferentially referred to the strong interaction configurations (up to a few angstroms) with very limited reference to the longrange ones.

Clear examples of the difference between the MN and the MF2 potentials are illustrated in the plots of Figure 2, where V inter of MF1 is also reported. In the figure a selection of stretched (rhs column), equilibrium (central column) and compressed (lhs col umn) diatomic molecules configurations of the O 2 + N 2 system are shown for the various H, X, T a , T b and I arrangements (see refs. 12 and 13 for the specification of the various configura tions). The plots clearly show that: • the curves of the MF1 PES differ from those of the MF2 one by small percentages both around the minimum and in the longrange attraction (by definition the MF1 and MF2 curves coincide when both molecules are at equilibrium, see central column panels);

• outside the potential energy regions in which the diatoms are in equilibrium (see lhs and rhs column panels), the long range MN potential is either flat or repulsive (i.e. its wells, if any, are negligible) and the repulsive wall is displaced to larger distances (of about 1.0 Å for compressed diatoms and of about 0.5 Å for stretched diatoms);

• in the potential energy regions in which the diatoms are in equilibrium the MF1 and MF2 PESs are identical and the MF and MN PESs are the most similar possible though with the MF wells being definitely deeper than the MN ones;

• for the H and X equilibrium geometries the wells of the MN PESs are displaced deeper inside the shorter range region while for the corresponding I one the displacement is to larger distances.

All in all, therefore, despite a certain similarity of the plots when the diatoms are at equilibrium, the MF and MN PESs differ significantly. In order to put the above made considerations on a more quantitative ground we have computedQ T , its glory pattern and the second virial coefficient B(T ) on the MN PES (details of the procedure and of its application to the MF PES are given in ref. 12) obtaining that:

• the absolute value of Q T computed on MN is at least a factor of 3 lower than the one obtained on the MF2 PES telling us that the longrange attraction of the former is too weak;

• the glory interference pattern (highly sensitive to the depth, position and shape of the longrange potential well) com puted on the MN PES falls completely outside the velocity range probed by the scattering experiment telling us that the area associated with the negative part of V inter is too small;

• the Boyle temperature of B(T ) computed on the MN PES is at least 300 K lower than the experimental value (400 K) telling us that the balance between the attraction and the repulsive part of the MN potential is inappropriate.

From the opacity function to rate coeffi cients

In order to investigate in detail the dynamical behaviour of the system on the MF2 and MN PESs we have carried out extended QCT calculations and analysed the behaviour of several batches of trajectories by plotting them as a function of the (boxed) im pact parameter b (opacity functions). Four families of opacity functions for the (v,u = 0|v 0 = v + ,u 0 = 0) VT transitions are plot ted in Figure 3 definitely larger than that computed for the MF2 PES (in accord with the fact that the repulsive wall of the MN PES protrudes to larger distances and therefore can involve larger values of b) even when the low impact parameter value of the probability of the former is smaller than that of the latter. At the same time, however, when increasing the temperature (that is in going from the lhs to the rhs column) the increase of the probability spreads almost equally well on both the height of the plateau and on its extension on b implying a kind of increased balancing between large b mechanisms and small b ones for which the opacity func tions computed on MF2 and MN tend on the average to converge thanks to a better balancing of the experiencing of V intra and V inter . The balancing between the effects of the intramolecular and inter molecular components of V offers a rationale for the fact that for v = 19 excitation and deexcitation of O 2 of one vibrational quan tum computed on the MN (that is more harmonic, i.e. harder) and on the MF2 (that is less harmonic, i.e. softer) PESs turn out to make the behaviour of the related opacity function quite differ ent (see the central panels of Fig. 3). In fact, for sufficiently de formed molecules, as is the case of O 2 in v = 19, the pronounced shift to larger distances of the MN V inter repulsive wall (see the rhs column of Fig. 2) makes, as discussed already in ref. On the ground of the additional information gained thanks to the more detailed nature of the opacity function, we are now bet ter equipped to discuss the absolute value (rather that the nor malized difference as done in Fig. 1) of the vibrational stateto state rate coefficients. For this reason we plot in Figure 4 the QCT k(v,u|v 0 ,u 0 ) rate coefficients computed on MF2 and MN. As a countercheck, we consider here for comparison also the hybrid PES (MF2MN) obtained by replacing in MF2 the diatomic poten tial of O 2 and N 2 with that of the MN PES (this replacement of the V intra with that of MN having a more harmoniclike nature should tends to overwhem the effect of the MF2 V inter . The values of the rate coefficients of the O 2 (v) + N 2 (u = 0) → O 2 (v 0 ) + N 2 (u 0 = 0) processes are plotted as a function of thev 0 -v jump at T = 1000, 3000 and 5000 K forv = 13, 19 and 25. The fig ure shows, indeed, that the rates computed on the hybrid MF2MN PES are always lower than those computed on the MF2 and MN PESs. The comparison between the MF2MN and the MF2 rates confirms that, as suggested by the analysis of the opacity function of Fig. 3, both the more repulsive longrange tail of V inter and the sharper shape of V intra differentiate the efficiency of the MN and the MF based PESs. However, while such differentiation is mini mal at temperatures higher than 1000 K (due to the negligibility of the negative part of V inter with respect to the high value of the energies involved) at lower temperatures (including T = 1000 K) the differentiation increases and this leads, as temperature de creases, to more than one order of magnitude difference in the value of rate coefficients.

Conclusions

The need for consolidating our trust on the Potential Energy Sur face to be used for massive calculations of thermal and sub thermal properties of the O 2 (v) + N 2 (u) nonreactive collision processes and in particular for the study of quasiresonant colli sion and their nonArrhenius behaviour at low temperature using quantum classical means, has led us to investigate in detail the characteristics of the entrance channel of the available PESs. For this reason we started by comparing VV and VT rate coefficients computed on our more traditional MF2 formulation of the PES based on its separation into a V intra and a V inter component and those computed on the recent MN one formulated as a polynomial of exponentialgaussian variables based on the fitting of a large number of ab initio points. Despite the fact that the MN PES bears the advantage of describing in a homogeneous formalism all the channels involved in the O 2 (v) + N 2 (u) collisions (including dis sociations that cannot be completely addressed by the MF PESs), presently it does not accurately describe thermal and subthermal nonreactive and nondissociative energy transfer processes. Ac cordingly, while the nonreactive rate coefficients computed on MF and MN PESs agree within one order of magnitude for tem peratures higher than 1000 K, when temperature decreases going from T = 1000 K down the disagreement increases to the point of preventing the use of the MN potential for our planned compu tational quantumclassical versus quasiclassical study of the non Arrhenius behaviour of the related rate coefficients at low tem perature.

)a

  Departamento de Quimica Fisica, Universidad del Pais Vasco (UPV/EHU), 01006 Vi toria, Spain. Fax: 34 945 013014; Tel: 34 945 013063; Email: e.garcia@ehu.es b Dipartimento di Chimica, Biologia e Biotecnologie, Università di Perugia, 06123 Pe rugia, Italy

Fig. 1

 1 Fig. 1 Normalized to the maximum value differences of QCT rate coefficients k(v,u|v 0 ,u 0 ) computed on MN and MF2 for the O 2 (v = 13,19,25) + N 2 (u = 0) → O 2 (v 0 ) + N 2 (u 0 = 0) process plotted as a function of E v 0 -v at T = 1000 and 5000 K.

Fig. 2

 2 Fig. 2 Plot of the longrange interaction at different values of the O 2 and N 2 internuclear distances (r O 2 = 1.007 Å, r N 2 = 0.898 Å, left hand side panels; r O 2 = r e,O 2 = 1.207 Å, r N 2 = r e,N 2 = 1.098 Å, central panels; r O 2 = 1.407 Å, r N 2 = 1.298 Å, right hand side panels) as a function of the intermolecular distance R for the MF1 (red line), MF2 (black line) and MN (blue line) PESs and the various H, X, T a , T b and I configurations.

  by assigning to  the values 2, 1, 1, 2 (in going from the top to the bottom panel) at the temperatures of 3000 K (left hand side panels) and 5000 K (right hand side panels) for v = 13, 19 and 25. The figure shows that in general the opac ity function computed on the MN PES dies at values of b (b max )

Fig. 3

 3 Fig. 3 Opacity function for the O 2 (v = 13,19,25) + N 2 (u = 0) → O 2 (v 0 = v ±1, ± 2) + N 2 (u 0 = 0) process plotted as a function of the impact parameter b at T = 3000 and 5000 K.

Fig. 4

 4 Fig. 4 Vibrational statetostate rate coefficients for the O 2 (v = 13,19,25) + N 2 (u = 0) → O 2 (v 0 ) + N 2 (u 0 = 0) process computed at T = 1000, 3000 and 5000 K on the MN, MF2 and MF2MN PESs plotted as a function of the vibrational jump v 0 -v.

Table 1

 1 Ratio of the MN over the MF2 VV thermal rate coefficients for the O 2

  10, the maximum value of b for which the opacity function is not zero, about 50% larger than that on MF2. The same does not occur when O 2 is either at equilibrium or scarcely deformed (MN and MF2 V inter of the central column of Fig.2for O 2 in v = 13 are quite similar) or largely deformed because highly excited (O 2 in v = 25) in which case in which case the effect of a V inter variation is of opposite sign and comparable with that of V intra .

	Probability	0.10 0.15 0.20			v=13 MF2 v=13 MN v=19 MF2 v=19 MN v=25 MF2 v=25 MN T=3000 K (v,0|v-2,0)				T=5000 K v=25 MN v=25 MF2 v=19 MN v=19 MF2 v=13 MN v=13 MF2 (v,0|v-2,0)	0.10 0.15 0.20 Probability
		0.05										0.05
	Probability	0.00 0.10 0.15 0.20			v=13 MF2 v=13 MN v=19 MF2 v=19 MN v=25 MF2 v=25 MN T=3000 K (v,0|v-1,0)				T=5000 K v=25 MN v=25 MF2 v=19 MN v=19 MF2 v=13 MN v=13 MF2 (v,0|v-1,0)	0.00 0.10 0.15 0.20 Probability
		0.05										0.05
	Probability	0.00 0.10 0.15 0.20			v=13 MF2 v=13 MN v=19 MF2 v=19 MN v=25 MF2 v=25 MN T=3000 K (v,0|v+1,0)				T=5000 K v=25 MN v=25 MF2 v=19 MN v=19 MF2 v=13 MN v=13 MF2 (v,0|v+1,0)	0.00 0.10 0.15 0.20 Probability
		0.05										0.05
	Probability	0.00 0.10 0.15 0.20			v=13 MF2 v=13 MN v=19 MF2 v=19 MN v=25 MF2 v=25 MN T=3000 K (v,0|v+2,0)				T=5000 K v=25 MN v=25 MF2 v=19 MN v=19 MF2 v=13 MN v=13 MF2 (v,0|v+2,0)	0.00 0.10 0.15 0.20 Probability
		0.05										0.05
		0 0.00	1	2 Impact parameter /Å 3	4	5	0	1	2 Impact parameter /Å 3	4	5	0.00

Journal Name, [year], [vol.], 1-7 | 1

Journal Name, [year], [vol.], 1-7 | 3

| 1-7 Journal Name, [year], [vol.],

Journal Name, [year], [vol.], 1-7 | 5

| 1-7 Journal Name, [year], [vol.],

Journal Name, [year], [vol.], 1-7 | 7

Acknowledgements

We thank the authors of ref. 23 for making available their po tential energy surface through the POTLIB repository. 34 Thanks are also due to the European Grid Infrastructure (EGI) for sup porting the COMPCHEM Virtual Organization and the CMMST Virtual Research Community, and the OU Supercomputing Center for Education & Research (OSCER) at the University of Oklahoma (OU) for providing computing resources and services. E.G. ac knowledges financial support from the MINECO/FEDER of Spain under grant CTQ201565033P. F.P. acknowledges the Diparti mento di Chimica, Biologia e Biotecnologie dell'Università di Pe rugia and the "Fondazione Cassa Risparmio Perugia" (Codice Pro getto: 2015.0331.021 Ricerca Scientifica e Tecnologica). C.M. thanks The European Joint Doctorate on TCCM project ITNEJD 642294Theoretical Chemistry and Computational Modeling.


Heinz H Bauschke
email: heinz.bauschke@ubc.ca.

Patrick L Combettes

Dominikus Noll
email: noll@mip.ups-tlse.fr.

Joint minimization with alternating Bregman proximity operators

Keywords: 2000 Mathematics Subject Classification: Primary 90C25; Secondary 49M05, 49M45, 65K05, 90C30 Alternating minimization, alternating projections, Bregman distance, Bregman projection, left proximity operator, Moreau envelope, proximal point algorithm, prox operator, right proximity operator

A systematic study of the proximity properties of Bregman distances is carried out. This investigation leads to the introduction of a new type of proximity operator which complements the usual Bregman proximity operator. We establish key properties of these operators and utilize them to devise a new alternating procedure for solving a broad class of joint minimization problems. We provide a comprehensive convergence analysis of this algorithm. Our framework is shown to capture and extend various optimization methods.

Introduction 1.Standing assumptions

Throughout, X = R J is the standard Euclidean space with inner product •, • and induced norm

• , and

(1) f : X → ]-∞, +∞] is convex and differentiable on U = int dom f = ∅.

Recall that (see [START_REF] Censor | Parallel Optimization: Theory, Algorithms, and Applications[END_REF])

(2) D f : X × X → [0, +∞] : (x, y) → f (x) -f (y) -f (y), x -y , if y ∈ U ; +∞, otherwise
is the Bregman distance associated with f , also denoted by D for brevity. Let Γ 0 (X) be the set of all proper lower semicontinuous convex functions from X to ]-∞, +∞]. In addition, f satisfies the following standard properties:

A1 f ∈ Γ 0 (X) is a convex function of Legendre type, i.e., f is essentially smooth and essentially strictly convex in the sense of [START_REF] Rockafellar | Convex Analysis[END_REF]Section 26];

A2 f exists and is continuous on U ;

A3 D is jointly convex, i.e., convex on X × X;

A4 (∀x ∈ U ) D(x,
•) is strictly convex on U ;

A5 (∀x ∈ U ) D(x, •) is coercive, i.e., the lower level set {y ∈ X : D(x, y) ≤ η} is bounded, for every η ∈ R.

These assumptions allow us to encompass several important scenarios, see Example 2.5. Finally, ϕ and ψ are two functions such that

(3)

     ϕ ∈ Γ 0 (X), (∀y ∈ U ) ϕ(•) + D(•, y) is coercive, dom ϕ ∩ U = ∅, and
     ψ ∈ Γ 0 (X), (∀x ∈ U ) ψ(•) + D(x, •) is coercive, dom ψ ∩ U = ∅.

Problem statement

Bregman distances were introduced in [START_REF] Bregman | The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming[END_REF] as an extension to the usual discrepancy measure (x, y) → x -y 2 and have since found numerous applications in optimization, convex feasibility, convex inequalities, variational inequalities, monotone inclusions, equilibrium problems; see [START_REF] Bauschke | Bregman monotone optimization algorithms[END_REF][START_REF] Butnariu | Totally Convex Functions for Fixed Points Computation and Infinite Dimensional Optimization[END_REF][START_REF] Censor | Parallel Optimization: Theory, Algorithms, and Applications[END_REF] and the references therein. The problem under consideration in the present paper is the joint minimization problem (4) minimize Λ : (x, y) → ϕ(x) + ψ(y) + D(x, y) over U × U.

The optimal value of (4) and its set of solutions will be denoted by [START_REF] Bauschke | Legendre functions and the method of random Bregman projections[END_REF] p = inf Λ(U × U ) and S = (x, y) ∈ U × U : Λ(x, y) = p , respectively.

The objective function Λ in (4) consists of a separable term (x, y) → ϕ(x) + ψ(y) and of a coupling term D. This structure arises explicitly or implicitly in a variety of problems, for instance in the areas of image processing [START_REF] Aujol | Image decomposition into a bounded variation component and an oscillating component[END_REF][START_REF] Vese | Modeling textures with total variation minimization and oscillating patterns in image processing[END_REF], signal recovery [START_REF] Combettes | Inconsistent signal feasibility problems: Least-squares solutions in a product space[END_REF], statistics [START_REF] Byrne | Proximity function minimization using multiple Bregman projections, with applications to split feasibility and Kullback-Leibler distance minimization[END_REF][START_REF] Csiszár | Information geometry and alternating minimization procedures[END_REF][START_REF] Iusem | A short convergence proof of the EM algorithm for a specific Poisson model[END_REF], mechanics [START_REF] Mercier | Inéquations Variationnelles de la Mécanique[END_REF], and wavelet synthesis [START_REF] Pesquet | Wavelet synthesis by alternating projections[END_REF]. Further applications will be described in Section 5.

Let ∆ = {(x, x) : x ∈ X}. Then it follows from Lemma 2.4(i) and A1 that [START_REF] Bauschke | Bregman monotone optimization algorithms[END_REF] (∀(x, y) ∈ U × U ) D(x, y) = 0 ⇔ (x, y) ∈ ∆.

Therefore, Problem (4) can be viewed as a relaxation of [START_REF] Bauschke | Iterating Bregman retractions[END_REF] minimize (x, y) → ϕ(x) + ψ(y) + ι ∆ (x, y) over U × U, which, in turn, is equivalent to the standard problem [START_REF] Bauschke | The asymptotic behavior of the composition of two resolvents[END_REF] minimize ϕ + ψ over U.

For the sake of illustration, let us consider the case when f = 1 2 • 2 , so that U = X and D : (x, y) → 1 2 x -y 2 . If ϕ and ψ are the indicator functions of two nonempty closed convex sets A and B, respectively, then [START_REF] Bauschke | The asymptotic behavior of the composition of two resolvents[END_REF] corresponds to the convex feasibility problem of finding a point in A ∩ B. When no such point exists, a sensible alternative is to look for a pair (x, y) ∈ A × B such that x -y = inf A -B . This formulation, which corresponds to (4), was proposed in [START_REF] Cheney | Proximity maps for convex sets[END_REF] and has found many applications in engineering [START_REF] Combettes | Inconsistent signal feasibility problems: Least-squares solutions in a product space[END_REF][START_REF] Mercier | Inéquations Variationnelles de la Mécanique[END_REF][START_REF] Pesquet | Wavelet synthesis by alternating projections[END_REF]. The algorithm devised in [START_REF] Cheney | Proximity maps for convex sets[END_REF] to solve this joint best approximation problem is the alternating projections method [START_REF] Bauschke | The method of forward projections[END_REF] fix x 0 ∈ X and set (∀n ∈ N) y n = P B (x n ) and x n+1 = P A (y n ).

More generally, let prox θ : x → argmin y θ(y)+ 1 2 x-y 2 be the proximity operator [START_REF] Moreau | Propriétés des applications prox[END_REF][START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF] associated with a function θ ∈ Γ 0 (X). In [START_REF] Acker | Convergence d'un schéma de minimisation alternée[END_REF], [START_REF] Bauschke | The method of forward projections[END_REF] was extended to the algorithm [START_REF] Bertsekas | Parallel and Distributed Computation: Numerical Methods[END_REF] fix x 0 ∈ X and set (∀n ∈ N)

y n = prox ψ (x n ) and x n+1 = prox ϕ (y n ) in order to solve (11) minimize (x, y) → ϕ(x) + ψ(y) + 1 2 x -y 2 over X × X.
The purpose of this paper is to introduce and analyze a proximal-like method to solve (4) under the assumptions stated above. The lack of symmetry of D prompts us to consider two single-valued operators defined on U , namely [START_REF] Butnariu | Iterative averaging of entropic projections for solving stochastic convex feasibility problems[END_REF] ← --

prox ϕ : y → argmin x∈U ϕ(x) + D(x, y) and --→ prox ψ : x → argmin y∈U ψ(y) + D(x, y).
The operators ← -prox ϕ and --→ prox ψ will be called the left and the right proximity operator, respectively. While left proximity operators have already been used in the literature (see [START_REF] Bauschke | Bregman monotone optimization algorithms[END_REF] and the references therein), the notion of a right proximity operator at this level of generality appears to be new. We note that [27, p. 26f] observes (but does not exploit) a superficial similarity between the iterative step of a multiplicative algorithm and the application of the right proximity operator --→ prox ψ in the Kullback-Leibler divergence setting (see Example 2.5(ii)), where ψ is assumed to be the sum of a continuous convex function and the indicator function of the nonnegative orthant in X.

In this paper, we shall provide a detailed analysis of these operators and establish key properties. With these tools in place, we shall be in a position to tackle (4) by alternating minimizations of Λ. We thus obtain the following algorithm [START_REF] Butnariu | Iterative methods of solving stochastic convex feasibility problems and applications[END_REF] fix x 0 ∈ U and set (∀n ∈ N)

y n = --→ prox ψ (x n ) and x n+1 = ← -- prox ϕ (y n ).
It is important to realize that it is quite nontrivial to see that this iteration is even well defined. The difficulty lies in guaranteeing that every iterate formally defined in [START_REF] Butnariu | Iterative methods of solving stochastic convex feasibility problems and applications[END_REF] lies again in U , so that the iterative update can be carried out. The crucial details of our analysis rely on various results on the interplay between the Bregman distance and the assumptions A1-A5 imposed on f . Armed with those results, we shall analyze the asymptotic behavior of this algorithm and, in particular, we shall establish convergence to a solution of (4). In the special case when ψ = 0, we recover variants and particular versions of the classical Bregman proximal method proposed in [START_REF] Censor | Proximal minimization algorithm with D-functions[END_REF] (see also [START_REF] Butnariu | Totally Convex Functions for Fixed Points Computation and Infinite Dimensional Optimization[END_REF] and [START_REF] Censor | Parallel Optimization: Theory, Algorithms, and Applications[END_REF]). Moreover, if we let ϕ = 0, we obtain a completely new proximal point method. We shall also extend and recover special cases of various known parallel decomposition algorithms, including least-squares techniques for inconsistent feasibility problems with finitely many sets. Let us also note that if ϕ is an indicator function and ψ is the sum of an indicator function and a differentiable convex function, then problem (4) reduces to a setting discussed in [START_REF] Eggermont | On EM-like algorithms for minimum distance estimation[END_REF]Remark 2.18]. However, the proofs in that manuscript are somewhat sketchy as several details are omitted. For instance, [START_REF] Eggermont | On EM-like algorithms for minimum distance estimation[END_REF] does not explain why the iteration ( 13) is well defined. Algorithm (13) may also be interpreted as a cyclic descent or nonlinear Gauss-Seidel method. However, the typical general convergence results for the latter methods (see, e.g., [10, Proposition 3.3.9]) fail to cover our main result (Theorem 4.4).

The paper is organized as follows. In Section 2, we collect the technical results required by our analysis. Left and right Bregman proximity operators are introduced and studied in Section 3. The asymptotic properties of Algorithm [START_REF] Butnariu | Iterative methods of solving stochastic convex feasibility problems and applications[END_REF] are investigated in Section 4. Finally, various applications and connections with previous works are described in Section 5.

Notation and conventions. Given a function g, denote its subdifferential map (resp. gradient map, conjugate function, and domain) by ∂g (resp. g or ∇g, g * and dom g). If C is a set, then we write ι C (resp. int C and cl C) for its indicator function (resp. interior and closure). We write N = {0, 1, 2, . . .} for the nonnegative integers. When dealing with the Boltzmann-Shannon entropy, it will be convenient to define 0 • ln(0) = 0, and to allow expressions such as x ≤ y, x • y, and x/y, which are understood coordinate-wise, for two vectors x and y in R J .

Auxiliary results

To make the paper self contained and to improve the presentation of the proofs of the main results in the later sections, we collect in this section several technical results. Lemma 2.1 Let g : X → ]-∞, +∞] be a proper convex function with V = int dom g. Proof. The inequality inf g(C) ≥ inf g(C) is clear. Since C ∩ dom g = ∅, the convexity of dom g and [START_REF] Rockafellar | Convex Analysis[END_REF]Corollary 6.3.2] imply that there exists c ∈ C ∩ ri dom g. Now fix x ∈ C ∩ dom g and note that, by [START_REF] Rockafellar | Convex Analysis[END_REF]Theorem 6.1], [START_REF] Butnariu | Totally Convex Functions for Fixed Points Computation and Infinite Dimensional Optimization[END_REF] ]x, c] ⊂ C ∩ ri dom g.

(i) If V = ∅ and g is differentiable on V , then g is continuous on V .
Next, we define, for every α ∈ ]0, 1], (i) Energy:

x α = (1 -α)x + αc ∈ C ∩ ri dom g. It follows from the segment continuity property [40, Theorem 7.5] that g(x) = lim α→0 + g(x α ). Thus, g(x) ≥ inf g(C). We conclude that inf g(C) ≥ inf g(C). Lemma 2.4 (i) (∀x ∈ X)(∀y ∈ U ) D(x, y) = 0 ⇔ x = y. (ii) (∀y ∈ U ) D(•, y) is coercive. (iii) If x ∈ U and (y n ) n∈N is a sequence in U such that y n → y ∈ bdry U , then D(x, y n ) → +∞.
If f : x → 1 2 x 2 , then U = X and D(x, y) = 1 2 x -y 2 .
(ii) Boltzmann-Shannon entropy:

If f : x → J j=1 ξ j ln(ξ j ) -ξ j , then U = {x ∈ X :
x > 0} and one obtains the Kullback-Leibler divergence

D(x, y) = J j=1 ξ j ln(ξ j /η j ) -ξ j + η j , if x ≥ 0 and y > 0; +∞, otherwise. (iii) Fermi-Dirac entropy: If f : x → J j=1 ξ j ln(ξ j ) + (1 -ξ j ) ln(1 -ξ j ), then U = {x ∈ X : 0 < x < 1} and D(x, y) = J j=1 ξ j ln(ξ j /η j ) + (1 -ξ j ) ln (1 -ξ j )/(1 -η j ) , if 0 ≤ x ≤ 1 and 0 < y < 1; +∞, otherwise.
Lemma 2.6 Suppose that x ∈ X and {u, v} ⊂ U . Then:

(15) D(x, v) = D(x, u) + D(u, v) + f (v) -f (u), u -x .
Moreover, D is continuous on U × U and D(u, •) ∈ Γ 0 (X).

Proof. The proof of the identity ( 15) is clear from [START_REF] Aujol | Image decomposition into a bounded variation component and an oscillating component[END_REF]. The continuity of D on U × U follows from Lemma 2.1(i) The function D(u, •) is convex by A3, and proper since u ∈ U . To verify lower semicontinuity of D(u, •), it suffices -in view of (2) -to take a sequence (y n ) n∈N in U that converges to y ∈ cl(U ) and to show that

D(u, y) ≤ lim D(u, y n ). If y ∈ U , then D(u, y n ) → D(u, y) by continuity of D on U × U . If y ∈ bdry(U ), then Lemma 2.4(iii) implies that D(u, y n ) → +∞ = D(u, y).
Identity ( 15) is also known as the "three points identity", see [START_REF] Chen | Convergence analysis of a proximal-like minimization algorithm using Bregman functions[END_REF]Lemma 3.1]. The next result follows by expanding (2) and some calculus.

Lemma 2.7 Take z ∈ U and h ∈ X. Then:

(16) lim t→0 + D(z, z + th) t = 0 = lim t→0 + D(z + th, z) t .
Because of A2 and A3, the function D f conforms to (1) and therefore its Bregman distance D D f , which will play a central role in our analysis, is well-defined.

Lemma 2.8 [9, Lemma 2.9] Take {x, y, u, v} ⊂ U . Then:

(17) D D f (x, y), (u, v) = D f (x, y) + D f (x, u) -D f (x, v) + f (v)(u -v), y -v . Moreover, D D f is continuous on U 4 .
Note that D f itself does not satisfy the counterparts of properties A1-A5; for instance, strict convexity fails as we shall see shortly. However, the expression for D D f becomes simpler when we deal with the energy or the Boltzmann-Shannon entropy (which are defined in Example 2.5):

Example 2.9 [9, Example 2.12] Take {x, y, u, v} ⊂ U . Then:

(i) If f is the energy, then D D f (x, y), (u, v) = D f x, y + (u -v) . (ii) If f is the Boltzmann-Shannon entropy, then D D f (x, y), (u, v) = D f x, yu/v ,
where the product and the quotient is taken coordinate-wise.

We do not know whether a similar simplification can be obtained for the Fermi-Dirac entropy.

Lemma 2.10 Let θ : X → ]-∞, +∞] be convex and x ∈ X be such that dom θ ∩ U = ∅ and

θ(•) + D(x, •) is coercive. Suppose (y n ) n∈N is a sequence in U such that θ(y n ) + D(x, y n ) n∈N is bounded.
Then (y n ) n∈N is bounded and all its cluster points belong to U .

Proof. The coercivity assumption implies the boundedness of (y n ) n∈N . Now let y be a cluster point of (y n ) n∈N , say y kn → y. We argue by contradiction and assume that y ∈ bdry U . By Lemma 2.1(ii), the function θ has an affine minorant, say a. On the other hand, Lemma 2.4(iii) implies that D(x, y kn ) → +∞. Hence -∞ ← θ(y kn ) ≥ a(y kn ) → a(y) > -∞, which is contradictory.

(∀x ∈ C)(∀n ∈ N) D(x, y n+1 ) ≤ D(x, y n ). (18)
Then (y n ) n∈N converges to a point in C if and only if all cluster points of (y n ) n∈N lie in C.

Lemma 2.12 Take θ ∈ Γ 0 (X) such that dom θ ∩ U = ∅. Consider the following properties:

(a) dom θ ∩ U is bounded. (b) inf θ(U ) > -∞.
(c) f is supercoercive, i.e., lim

x →+∞ f (x)/ x = +∞. (d) (∀x ∈ U ) D(x, •) is supercoercive.
Then:

(i) If any of the conditions (a), (b), or (c) holds, then

(19) (∀y ∈ U ) θ(•) + D(•, y) is coercive or, equivalently, (20)
(∀x ∈ U ) θ(•) + D(x, •) is coercive. Proof. By Lemma 2.1(ii), there exists (z * , α) ∈ X × R such that θ(•) ≥ z * , • + α. (a) ⇒ (b): By Cauchy-Schwarz, we have inf θ(U ) = inf θ(dom θ ∩ U ) ≥ -z * • sup dom θ ∩ U + α > -∞ as dom θ ∩ U is bounded. (b) ⇒ (19) (21)
: Suppose to the contrary that there exists a sequence Lemma 2.13 Let g : X → ]-∞, +∞] be proper, coercive, and convex. Then inf g(X) > -∞.

(x n ) n∈N in dom f such that x n → +∞ and (θ(x n ) + D(x n , y)) n∈N is bounded. Now observe that µ = inf θ(dom f ) = inf θ(U ) > -∞
∈ U ) θ + f -f (y), • is coercive ⇔ (∀y ∈ U ) 0 ∈ dom(θ + f -f (y), • ) * ⇔ (∀y ∈ U ) f (y) ∈ int dom (f + θ) * ⇔ (20
θ(x) + D(x, y) ≥ α -f (y) + f (y), y + x f (x) x -z * -f (y) → +∞.
Proof. Set µ = inf g(X) and take a sequence (x n ) n∈N in X such that g(x n ) → µ. Since g is coercive, (x n ) n∈N is bounded and therefore it has a cluster point, say x kn → x. By Lemma 2.1(ii), there exists an affine minorant of g, say a. Then µ ← g(x kn ) ≥ a(x kn ) → a(x) > -∞.

3 Bregman envelopes and proximity operators

Definition 3.1 Take θ : X → ]-∞, +∞]. The left Bregman envelope of θ is (24) ←-env θ : X → [-∞, +∞] : y → inf x∈X θ(x) + D(x, y),
and the right Bregman envelope of θ is

(25) -→ env θ : X → [-∞, +∞] : x → inf y∈X θ(y) + D(x, y).
Let us provide two illustrations of these definitions.

Example 3.2 Suppose f = 1 2 • 2 and take θ : X → ]-∞, +∞]. Then D : (x, y) → 1 2 x -y 2 and ←-env θ = -→ env θ = θ ( 1 2 • 2 ) is the Moreau envelope of θ [41, Section 1.G]. Example 3.3 Let C be a subset of X. The left Bregman distance to C is defined by (26) ← - D C = ←-env ι C : y → inf x∈C D(x, y),
and the right Bregman distance to C is defined by

(27) -→ D C = -→ env ι C : x → inf y∈C D(x, y).
The following propositions collect some basic properties of Bregman envelopes.

Proposition 3.4 Let θ : X → ]-∞, +∞] be such that dom θ ∩ U = ∅. Then:

(i) dom ←env θ = U and (∀y ∈ U ) ←env θ (y) ≤ θ(y).

(ii) dom -→ env θ = dom f and (∀x ∈ U ) -→ env θ (x) ≤ θ(x).

(iii) Suppose that θ is convex. Then ←env θ is convex and continuous on U . If, in addition, (19) holds, i.e., (∀y ∈ U ) θ(•) + D(•, y) is coercive, then ←env θ is proper.

(iv) Suppose that θ is convex. We now provide additional conditions guaranteeing that the infima in Definition 3.1 are uniquely attained in U .

) n∈N in U such that z n → z. Then ← - D C (z n ) ≡ 0, but ← - D C (z) = +∞. Hence ← - D C
is not lower semicontinuous at z. Therefore, left Bregman envelopes need not belong to Γ 0 (X). Left (i.e., classical Bregman) proximity operators have already been used in several works, see e.g., [START_REF] Bauschke | Bregman monotone optimization algorithms[END_REF], [START_REF] Censor | Proximal minimization algorithm with D-functions[END_REF], [19, Chapter 3], [START_REF] Chen | Convergence analysis of a proximal-like minimization algorithm using Bregman functions[END_REF], [START_REF] Eckstein | Nonlinear proximal point algorithms using Bregman functions, with applications to convex programming[END_REF], and [START_REF] Kiwiel | Proximal minimization methods with generalized Bregman functions[END_REF]. On the other hand, the notion of a right proximity operator appears to be new. We note that [27, p. 26f] discusses the formal resemblance between a multiplicative algorithm and the right proximity operator in the Kullback-Leibler divergence setting when θ is the sum of a continuous convex function and the indicator function of the nonnegative orthant in X. Further, we point out below (see Example 3.9) that in the special case of indicator functions, the right proximity operator was previously considered in [START_REF] Bauschke | The method of forward projections[END_REF].

Example 3.8 Suppose that f = 1 2 • 2 and take θ ∈ Γ 0 (X). Since f is supercoercive, it follows from Lemma 2.12(i) that ( 19) is satisfied, and we obtain Moreau's proximity operator [START_REF] Moreau | Propriétés des applications prox[END_REF][START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF][START_REF] Rockafellar | Variational Analysis[END_REF]:

← -- prox θ = --→ prox θ = (Id +∂θ) -1 .
Example 3.9 Let C ⊂ X be a closed convex set such that C ∩ U = ∅. Since ι C is bounded below, Lemma 2.12 guarantees that ( 19) and ( 21) hold; furthermore, ← -prox ι C = ← -P C is the (left, i.e.,) classical Bregman projector onto C [START_REF] Bauschke | Legendre functions and the method of random Bregman projections[END_REF][START_REF] Bregman | The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming[END_REF][START_REF] Byrne | Proximity function minimization using multiple Bregman projections, with applications to split feasibility and Kullback-Leibler distance minimization[END_REF][START_REF] Censor | Block-iterative algorithms with underrelaxed Bregman projections[END_REF] and --→ prox ι C = -→ P C is the right Bregman projector onto C [START_REF] Bauschke | Iterating Bregman retractions[END_REF][START_REF] Bauschke | The method of forward projections[END_REF]. Note that in the last two references, the left and right Bregman projector are called backward and forward Bregman projector. However, because of possible ambiguity in the context of splitting methods, the notions of left and right Bregman projector are preferable.

From now on, we also utilize the notation ∇g to describe the derivative g of a given function g. This increases readability when g is a more complicated expression. The following properties will be needed later. Proposition 3.10 Let θ ∈ Γ 0 (X) be such that dom θ ∩ U = ∅.

(i) Suppose that (19) holds. Then for every (x, y) ∈ U 2 , the following conditions are equivalent:

(a) x = ← -- prox θ (y); (b) 0 ∈ ∂θ(x) + f (x) -f (y); (c) (∀z ∈ X) f (y) -f (x), z -x + θ(x) ≤ θ(z).
Moreover,

(30) ← -- prox θ = (f + ∂θ) -1 • f is continuous on U .
(ii) Suppose that (21) holds. Then for every (x, y) ∈ U 2 , the following conditions are equivalent:

(a) y = --→ prox θ (x); (b) 0 ∈ ∂θ(y) + f (y)(y -x); (c) (∀z ∈ X) f (y)(x -y), z -y + θ(y) ≤ θ(z).
Moreover, --→ prox θ is continuous on U .

Proof. (i): We verify only continuity as the equivalence of (a)-(c), as well as the identity [START_REF] Kiwiel | Proximal minimization methods with generalized Bregman functions[END_REF]

∇(f + θ) * • f = ∂(f + θ) -1 • f = (f + ∂θ) -1 • f = ← -- prox θ is continuous on U . (ii):
The equivalence of (a)-(c) is clear from [START_REF] Iusem | A short convergence proof of the EM algorithm for a specific Poisson model[END_REF] and convex calculus. To establish the continuity of --→ prox θ on U , pick a sequence (x n ) n∈N in U converging to x ∈ U and set y n = --→ prox θ (x n ), for all n ∈ N. Take q ∈ dom θ ∩ U . Then, using Lemma 2.6, Lemma 2.8, and item (ii)(c), we obtain

D(x, q) ← D(x, q) + D(x, x n ) = D D f (x, q), (x n , y n ) + D(x, y n ) -f (y n )(x n -y n ), q -y n ≥ D(x, y n ) -f (y n )(x n -y n ), q -y n ≥ D(x, y n ) + θ(y n ) -θ(q). (31)
It follows that θ(y n ) + D(x, y n ) n∈N is bounded. By ( 21) and Lemma 2.10, the sequence (y n ) n∈N is bounded and its cluster points belong to U . Let us extract a converging subsequence, say y kn → y ∈ U . In view of item (ii)(c), we have (32) (∀z ∈ X)(∀n ∈ N) f (y kn )(x kn -y kn ), z -y + θ(y kn ) ≤ θ(z).

We let n tend to +∞ in (32), use continuity of f (see A2) and lower semicontinuity of θ to obtain

(33) (∀z ∈ X) f (y)(x -y), z -y + θ(y) ≤ θ(z).
The equivalence between items (ii)(a) and (ii)(c) now results in y = --→ prox θ (x).

Remark 3.11

The proof of continuity of ← -prox θ presented above extends Lewis' unpublished proof [32] of continuity of ← -P C , where C ⊂ X is a closed convex set such that C ∩ U = ∅. Furthermore, the continuity of ← -P C when f is the Boltzmann-Shannon entropy was first established in [START_REF] Butnariu | Iterative averaging of entropic projections for solving stochastic convex feasibility problems[END_REF].

Proposition 3.12 Let θ ∈ Γ 0 (X) be such that dom θ ∩ U = ∅. (i) If (19) holds, then ←-env θ is differentiable on U and (34) (∀y ∈ U ) ∇ ←-env θ (y) = f (y)(y -← -- prox θ (y)).
(ii) If (21) holds, then -→ env θ is differentiable on U and After dividing this chain of inequalities by t, we take the limit as t → 0 + . Lemma 2.7, A2, and the continuity of P (see Proposition 3.10(i)) imply that the leftmost limit is the same as the rightmost limit, namely f (y)(h), y -P (y) . It follows that

(35) (∀x ∈ U ) ∇ -→ env θ (x) = f (x) -f ( --→ prox θ (x)). Proof. (i): Fix y ∈ U , h ∈ X,
+ ←-env θ (y + th) -←-env θ (y) t = f (y)(h), y -P (y) . (36) lim t→0
(ii): Fix x ∈ U , h ∈ X, and t ∈ ]0, +∞[ such that x + th ∈ U . We set P = --→ prox θ and obtain, using Lemma 2.6 twice, D(x + th, x) + f (x) -f (P (x + th)), th = D(x + th, P (x + th)) -D(x, P (x + th))

= θ(P (x + th)) + D(x + th, P (x + th)) -θ(P (x + th)) -D(x, P (x + th))

≤ θ(P (x + th)) + D(x + th, P (x + th)) -θ(P (x)) -D(x, P (x)) = -→ env θ (x + th) --→ env θ (x) ≤ θ(P (x)) + D(x + th, P (x)) -θ(P (x)) -D(x, P (x)) = D(x + th, P (x)) -D(x, P (x)) = D(x + th, x) + f (x) -f (P (x)), th
Let us divide this chain of inequalities by t, and then take the limit as t → 0 + . Lemma 2.7 and the continuity of P (see Proposition 3.10(ii)) imply that the leftmost limit is the same as the rightmost limit, namely f (x) -f (P (x)), h . Thus [START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF] lim

t→0 + -→ env θ (x + th) --→ env θ (x) t = f (x) -f (P (x)), h .
Remark 3.13 Special cases of Proposition 3.12(i) have been observed previously in the literature; see, for instance, [42, Theorem 4.1(b)] when inf θ(U ) > -∞ (so that [START_REF] Censor | Parallel Optimization: Theory, Algorithms, and Applications[END_REF] holds by Lemma 2.12(i)) and [19, Proposition 3.2.3] when dom θ = X. The differentiability of -→ D C was established in [START_REF] Butnariu | Iterative averaging of entropic projections for solving stochastic convex feasibility problems[END_REF] for the case where f is the Boltzmann-Shannon entropy.

Let us now provide two illustrations of Proposition 3.12.

Example 3.14 If θ ∈ Γ 0 (X) and f = 1 2 • 2 , then Proposition 3.12 reduces to Moreau's gradient formula [START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF]Proposition 7.d

], namely ∇ θ ( 1 2 • 2 ) = Id -(Id +∂θ) -1 .
Example 3.15 Let C ⊂ X be a closed convex set such that C ∩ U = ∅ and take {x, y} ⊂ U . In view of Examples 3.3 and 3.9, setting θ = ι C in Proposition 3.12 yields ∇ ← -

D C (y) = f (y)(y- ← - P C (y)) and ∇ -→ D C (x) = f (x) -f ( -→ P C (x)).
As the following proposition shows, left and right envelopes and prox operators arise naturally in connection with our basic problem (4). Let us introduce the two auxiliary relaxed problems [START_REF] Pesquet | Wavelet synthesis by alternating projections[END_REF] minimize ←env ϕ + ψ over U and (39) minimize ϕ + -→ env ψ over U.

Their solution sets will be denoted by ( 40)

F = y ∈ U | ←-env ϕ (y) + ψ(y) = inf( ←-env ϕ + ψ)(U ) and (41)
E = x ∈ U | ϕ(x) + -→ env ψ (x) = inf(ϕ + -→ env ψ )(U ) ,
respectively. In the standard metric setting, i.e., f = 1 2 • 2 , (38) and ( 39) are the two classical partial Moreau regularizations of ( 8), e.g., [START_REF] Mahey | Partial regularization of the sum of two maximal monotone operators[END_REF]. We now relate the sets E and F to the set S, defined in [START_REF] Bauschke | Legendre functions and the method of random Bregman projections[END_REF], as well as to the operators ← -prox ϕ and --→ prox ψ .

Proposition 3. [START_REF] Byrne | Proximity function minimization using multiple Bregman projections, with applications to split feasibility and Kullback-Leibler distance minimization[END_REF] The following properties hold:

(i) E and F are convex.

(ii) (∀(x, y) ∈ U × U ) (x, y) ∈ S ⇔ x = ← -- prox ϕ (y) and y = --→ prox ψ (x) . (iii) E = Fix ← -- prox ϕ • --→ prox ψ and F = Fix --→ prox ψ • ← -- prox ϕ . (iv) (∀(x, y) ∈ E × F ) x, --→ prox ψ (x) ∈ S and ← -- prox ϕ (y), y ∈ S.
Example 4.12 Take (x, y), ( x, y) ⊂ S. Then:

(i) If f is the energy, then y -x = y -x.

(ii) If f is the Boltzmann-Shannon entropy, then y/x = y/ x.

Proof. Combine Corollary 4.11, Example 2.9, and [START_REF] Bauschke | Bregman monotone optimization algorithms[END_REF].

We now turn to an optimization problem dual to (4), namely, to determine

* ) = f (y) -f (x), f (y)(x -y) . Consequently, if (x n , y n ) n∈N is generated by (46), then f (y n ) -f (x n ), f (y n )(x n -y n ) → (x * , y * ).
Proof. Let (x * , y * ) be a solution of (63) and take (x, y) ∈ S, i.e., (x, y) solves (4). Then, using the Fenchel Duality Theorem The "Consequently" part follows from Theorem 4.4 and the continuity of f and f (see A2).

Remark 4.14 Take (x, y) ∈ S. Proposition 4.13 asserts that (x * , y * ) = f (y)-f (x), f (y)(x-y) is the unique solution of (63).

(i) If f is the energy, then (x * , y * ) = (y -x, x -y) = (x * , -x * ).

(ii) If f is the Boltzmann-Shannon entropy, then (x * , y * ) = ln(y/x), x/y -1 = x * , exp(-x * ) -1 .

(iii) If f is the Fermi-Dirac entropy, then (x * , y * ) = ln y/x (1 -y)/(1 -x) , x -y y(1 -y) .

Note that (i) and (ii) combined with the uniqueness of (x * , y * ) lead to an alternative proof of the identities in Example 4.12.

5 Applications and connections with previous works

Preliminaries

In this section, we discuss various applications of Theorem 4.4 revolving around the basic constrained optimization problem (66) minimize θ over C ∩ U, where θ ∈ Γ 0 (X), dom θ ∩ U = ∅, and C is a closed convex subset of X such that C ∩ U = ∅. We are going to consider increasingly specialized realizations of (66). First, suppose that C = X, that I is an ordered finite index set, and that θ can be decomposed as θ = i∈I ω i θ i , where

(∀i ∈ I) θ i ∈ Γ 0 (X) and dom θ i ∩ U = ∅, (67)
and the weights {ω i } i∈I ⊂ ]0, 1] satisfy i∈I ω i = 1. Then (66) becomes (68) minimize i∈I ω i θ i over U.

In particular, if we set

(69) (∀i ∈ I) θ i = 1 2 max{0, g i } 2
, where g i ∈ Γ 0 (X) and dom g i ∩ U = ∅, then (68) reduces to solving a system of convex inequalities, namely, (70) find x ∈ U such that max

i∈I g i (x) ≤ 0.
Furthermore, if we set (g i ) i∈I = (ι S i ) i∈I , where (S i ) i∈I is a family of closed convex sets such that, for every i ∈ I, S i ∩ U = ∅, then (70) reduces to the basic convex feasibility problem

(71) find x ∈ U ∩ i∈I S i .
We shall employ a product space setup initially introduced in [START_REF] Pierra | Decomposition through formalization in a product space[END_REF] for metric projection methods and revisited in [19, Section 5.9] in the context of feasibility problems with Bregman distances.

Denote the standard Euclidean product space X I by X and write x = (x i ) i∈I , for x ∈ X (hence,

x 2 = i∈I x i 2 ). Now define (72)            ∆ = (x, . . . , x) ∈ X : x ∈ X , f : X → ]-∞, +∞] : x → i∈I ω i f (x i ), U = U I = int dom f , θ : X → ]-∞, +∞] : x → i∈I ω i θ i (x i ).
Then f induces a Bregman distance D on X defined by

(73) (∀(x, y) ∈ X × X) D(x, y) = i∈I ω i D(x i , y i ).
It is straightforward to verify that f satisfies A1-A5, that dom θ ∩ U = ∅, and that ∆ ∩ U = ∅.

When (66) is not guaranteed to have solutions, one can turn to the Bregman relaxations ( 38) and [START_REF] Pierra | Decomposition through formalization in a product space[END_REF]. Let us now explore these relaxed formulations and derive algorithms to solve them.

Left Bregman relaxation

Setting ϕ = θ and ψ = ι C in ( 8) yields (66). Accordingly, the left relaxation of (66) with respect to θ is derived from [START_REF] Pesquet | Wavelet synthesis by alternating projections[END_REF]

ω i ←-env θ i over U
as a left relaxation of (68) under the assumption that the functions (θ i ) i∈I satisfy (19) (hence so does θ). The convenience of the product space setup of (72) becomes apparent in the following result.

Proposition 5.2 Let (x, y) ∈ U 2 . Then:

(i) ←env θ (y) = i∈I ω i ←env θ i (y i ).

(ii)

-→ P ∆ (x) = (z, . . . , z), where z = i∈I ω i x i .

(iii) ← -prox θ (y) = ← -prox θ i (y i ) i∈I .

(iv) Fix

z ∈ Fix -→ P ∆ • ← -- prox θ ⇔ z = i∈I ω i ← -- prox θ i (z) ⇔ 0 = f (z) i∈I ω i z -← -- prox θ i (z) = i∈I ω i f (z) z -← -- prox θ i (z) ⇔ 0 = i∈I ω i ∇ ←-env θ i (z) = ∇ i∈I ω i ←-env θ i (z)
⇔ z solves (76).

(79)

Important conclusions can be drawn from the above proposition. First, item (i) asserts that Problem (76) in X is equivalent to (80) minimize ←env θ over ∆ ∩ U in X. This is a special case of (74) for which Algorithm (75) becomes (81) y 0 ∈ U and (∀n ∈ N) y n+1 = -→ P ∆ • ← -prox θ (y n ).

A direct application of Propositions 5.1 and 3.16 shows that (y n ) n∈N converges to a fixed point of -→ P ∆ • ← -prox θ , provided that such a point exists. In view of Proposition 5.2(ii)-(iv), we therefore obtain the following proposition.

We conclude with an application of this proposition to the right Bregman relaxation of ( 71 Then (x n ) n∈N converges to a point in G.

Connections with previous works

We conclude by providing links between the results of Sections 5.2 and 5.3 and previous works.

Remark 5.9

(i) When f = 1 2 • 2 , Algorithms (82) and (90) coincide with [39, Algorithm 3.1] (see also [START_REF] Mahey | Partial regularization of the sum of two maximal monotone operators[END_REF]).

(ii) When f = 1 2 • 2 , (83) and (91) reduce to the problem of minimizing a weighted sum of the squares of the distances to the sets whereas (84) and (92) reduce to the method of barycentric metric projections. This framework has been explored from different viewpoints in [START_REF] Bauschke | On the convergence of von Neumann's alternating projection algorithm for two sets[END_REF][START_REF] Combettes | Inconsistent signal feasibility problems: Least-squares solutions in a product space[END_REF][START_REF] De Pierro | A parallel projection method for finding a common point of a family of convex sets[END_REF].

(iii) Algorithm (84) has been studied at various levels of generality in [START_REF] Butnariu | Iterative averaging of entropic projections for solving stochastic convex feasibility problems[END_REF][START_REF] Butnariu | Iterative methods of solving stochastic convex feasibility problems and applications[END_REF][START_REF] Byrne | Proximity function minimization using multiple Bregman projections, with applications to split feasibility and Kullback-Leibler distance minimization[END_REF]. Proposition 5.4 is a particular case of the comparable results in [START_REF] Butnariu | Iterative methods of solving stochastic convex feasibility problems and applications[END_REF][START_REF] Butnariu | Totally Convex Functions for Fixed Points Computation and Infinite Dimensional Optimization[END_REF].

(iv) Algorithm (92) is discussed in [START_REF] Bauschke | Iterating Bregman retractions[END_REF][START_REF] Byrne | Proximity function minimization using multiple Bregman projections, with applications to entropy optimization and Kullback-Leibler distance minimization[END_REF].

(

  ii) The function g admits an affine minorant. Proof. (i): [40, Theorem 25.5]. (ii): [40, Corollary 12.1.2]. Lemma 2.2 [40, Corollary 14.2.2] Let g ∈ Γ 0 (X). Then g is coercive if and only if 0 ∈ int dom g * . Lemma 2.3 Let C be an open convex subset of X and let g ∈ Γ 0 (X) be such that C ∩ dom g = ∅. Then inf g(C) = inf g(C).

Proof. (i): [ 5 ,Example 2 . 5 [ 9 ,

 5259 Theorem 3.7.(iv)]. (ii): [5, Theorem 3.7.(iii)]. (iii): [5, Theorem 3.8.(i)]. Example 2.16] Assumptions A1-A5 hold in the following cases, where x = (ξ j ) 1≤j≤J and y = (η j ) 1≤j≤J are two generic points in R J .

Lemma 2 . 11 [ 9 ,

 2119 Lemma 2.20] or[START_REF] Bauschke | Bregman monotone optimization algorithms[END_REF] Section 4.1] Suppose that ∅ = C ⊂ U and (y n ) n∈N is a sequence in U which is Bregman monotone with respect to C, i.e.,

  ran f ⊂ int dom (f + θ) * . (ii) If any of the conditions (a), (b), or (d) holds, then

  (d) ⇒ (21): Letting y → +∞, we obtain (23) θ(y) + D(x, y) ≥ α + y D(x, y) y -z * → +∞.

Proposition 3 . 5 Remark 3 . 6

 3536 Let θ ∈ Γ 0 (X) such that dom θ ∩ U = ∅. (i) Suppose that y ∈ U and that θ(•) + D(•, y) is coercive. Then there exists a unique point z ∈ U such that ←env θ (y) = θ(z) + D(z, y). (ii) Suppose that x ∈ U and that θ(•) + D(x, •) is coercive. Then there exists a unique point z ∈ U such that -→ env θ (x) = θ(z) + D(x, z). Proof. (i): Apply [6, Proposition 3.21.(ii)], [6, Proposition 3.23.(v)(b)], and [6, Proposition 3.22.(ii)(d)]. (ii): Set µ = -→ env θ (x). By Lemma 2.13, µ ∈ R. Take (z n ) n∈N in U such that θ(z n ) + D(x, z n ) → µ. Then by Lemma 2.10, (z n ) n∈N has a cluster point in U , say z kn → z ∈ U . However, by Lemma 2.6, θ(•) + D(x, •) is lower semicontinuous at z and therefore µ ≤ θ(z) + D(x, z) ≤ lim(θ(z kn ) + D(x, z kn )) = µ. Furthermore, A4 implies that θ(•) + D(x, •) is strictly convex, which secures the uniqueness of z. Suppose that U = X (as happens for the two entropies in Example 2.5) and set C = cl(U ) in Example 3.3. Now pick z ∈ bdry U and (z n

Proposition 3 .Definition 3 . 7

 337 5 allows us to define the following operators on U . Let θ ∈ Γ 0 (X) be such that dom θ ∩ U = ∅. If (19) holds, then the left proximity operator associated with θ is (28) ← -prox θ : U → U : y → argmin x∈X θ(x) + D(x, y). If (21) holds, then the right proximity operator associated with θ is (29) --→ prox θ : U → U : x → argmin y∈X θ(y) + D(x, y).

  and t ∈ ]0, +∞[ such that y + th ∈ U . For the sake of brevity, set P = ← -prox θ . Using Lemma 2.6 twice, we estimate D(y, y + th) + f (y + th) -f (y), y -P (y + th) = D(P (y + th), y + th) -D(P (y + th), y) = θ(P (y + th)) + D(P (y + th), y + th) -θ(P (y + th)) -D(P (y + th), y) ≤ θ(P (y + th)) + D(P (y + th), y + th) -θ(P (y)) -D(P (y), y) = ←env θ (y + th) -←env θ (y) ≤ θ(P (y)) + D(P (y), y + th) -θ(P (y)) -D(P (y), y) = D(P (y), y + th) -D(P (y), y) = D(y, y + th) + f (y + th) -f (y), y -P (y) .

  ,y * )∈X×X ϕ * (x * ) + ψ * (y * ) + D * (-x * , -y * ). This is precisely the standard Fenchel dual of the optimization problem (4) -the minus sign is simply added to ensure that the optimal values of the two problems coincide. Now (3) implies the standard constraint qualification which in the present setting states that int dom D f = U × U and dom (x, y) → ϕ(x) + ψ(y) = dom ϕ × dom ψ possess common points. Consequently, by [40, Corollary 6.3.2 and Theorem 31.1], the minimum in (63) is always attained. More information can be obtained when (4) has solutions: Proposition 4.13 Suppose that S = ∅. Then the minimum in (63) is attained at a unique point (x * , y * ) and, for every (x, y) ∈ S, (x * , y

[ 40 ,

 40 Theorem 31.1] and the Fenchel-Young inequality, we obtain 0 = ϕ(x) + ψ(y) + D(x, y) + ϕ * (x * ) + ψ * (y * ) + D * (-x * , -y * ) ≥ x * , x + y * , y + D(x, y) + D * (-x * , -y * ) x, y) + D * (-x * , -y * ) = (-x * , -y * ), (x, y) and thus, with the help of (42), we obtain (65) (-x * , -y * ) = ∇D(x, y) = f (x) -f (y), f (y)(y -x) .

  -→ P ∆ • ← -prox θ = (z, • • • , z) : z solves (76) . Proof. (i): By definition, ←env θ (y) = inf x∈X θ(x) + D(x, y) = inf x∈X i∈Iω i θ i (y i ) + D(x i , y i ) = i∈I ω i inf x i ∈X θ i (y i ) + D(x i , y i ) = i∈I ω i ←env θ i (y i ). (77) (ii): See [7, Example 3.16(ii)]. (iii): It follows from Proposition 3.10(i) thatx = ← -prox θ (y) ⇔ 0 ∈ ∂θ(x) + f (x) -f (y) = ω i ∂θ i (x i ) i∈I + ω i f (x i ) -ω i f (y i ) i∈I ⇔ (∀i ∈ I) 0 ∈ ∂θ i (x i ) + f (x i ) -f (y i ) ⇔ (∀i ∈ I) x i = ← -prox θ i (y i ). (78) (iv): Since Fix -→ P ∆ • ← -prox θ ⊂ ∆ ∩ U, let us fix z ∈ ∆ ∩ U, say z = (z, • • • , z), where z ∈ U . Then it follows from (ii), (iii), A1, A2, and (34) that

Proposition 5 . 8

 58 Suppose that the solution set G of (91) is nonempty and let (x n ) n∈N be a sequence generated by (92)x 0 ∈ U and(∀n ∈ N) x n+1 = ∇f * i∈I ω i ∇f -→ P S i (x n ) .

  by (b) and Lemma 2.3. Then we arrive at the contradiction +∞ > sup n∈N θ(x

n ) + D(x n , y) ≥ µ + sup n∈N D(x n , y) = +∞

since, by Lemma 2.4(ii), D(•, y) is coercive. (19) ⇔ (20): Using Lemma 2.2, we deduce that (19) ⇔ (∀y

  Then -→ env θ is convex and continuous on U . If, in addition,[START_REF] Cheney | Proximity maps for convex sets[END_REF] holds, i.e., (∀x ∈ U ) θ(•) + D(x, •) is coercive, then -→ env θ is proper.Proof. (i) and (ii) follow at once from Definition 3.1 and (2). (iii): A3 asserts that the function (y, z) → θ(z) + D(z, y) is convex. Hence, it follows from [41, Proposition 2.22.(a)] that the marginal function ←env θ is also convex. The continuity of ←env θ on U then follows from (i) and the fact that every convex function on X is continuous on the interior of its domain[START_REF] Rockafellar | Convex Analysis[END_REF] Theorem 10.1]. It is clear from (i) that ←env θ ≡ +∞. On the other hand it follows from[START_REF] Censor | Parallel Optimization: Theory, Algorithms, and Applications[END_REF] and Lemma 2.13 that

-∞ / ∈ ←env θ (X). (iv): Similar to (iii).

  are known, see e.g.[START_REF] Bauschke | Bregman monotone optimization algorithms[END_REF] Section 3.4]. Since f is Legendre by A1, it is essentially strictly convex and so is f + θ. By[START_REF] Rockafellar | Convex Analysis[END_REF] Theorem 26.3], (f + θ) * is essentially smooth. Now Lemma 2.1(i) implies that ∇(f + θ) * is continuous on int dom (f + θ) * and that f is continuous on U . Therefore,

Acknowledgments H. H. Bauschke's work was partially supported by the Natural Sciences and Engineering Research Council of Canada and by a NATO grant. The work of P. L. Combettes was partially supported by a NATO grant.

Proof. (i): In view of (3) and Proposition 3.4(iii)&(iv), E and F are convex, as sets of minimizers of convex functions. For the remainder of the proof, we fix (x, y) ∈ U × U . (ii): Since [START_REF] Teboulle | Entropic proximal mappings with applications to nonlinear programming[END_REF] ∇D(x, y) = f (x) -f (y), f (y)(y -x) ,

we obtain by standard convex calculus and by invoking items (i)(b) and (ii)(b) of Proposition 3.10 the chain of equivalences (x, y) ∈ S ⇔ (0, 0) ∈ ∂Λ(x, y) = ∂ϕ(x) + f (x) -f (y), ∂ψ(y) + f (y)(y -x)

⇔ 0 ∈ ∂ϕ(x) + f (x) -f (y) and 0 ∈ ∂ψ(y) + f (y)(y -x) [START_REF] Vese | Modeling textures with total variation minimization and oscillating patterns in image processing[END_REF] ⇔ x = ← -prox ϕ (y) and y = --→ prox ψ (x).

(iii): It follows from Proposition 3.12(ii) and Proposition 3.10(i) that

Likewise, it follows from Proposition 3.12(i) and Proposition 3.10(ii) that

(iv) follows at once from (ii) and (iii).

Alternating left and right proximity operators

Recall that the standing assumptions on our problem (4) are described by [START_REF] Bauschke | On the convergence of von Neumann's alternating projection algorithm for two sets[END_REF], that its solution set S and its optimal value p are defined in [START_REF] Bauschke | Legendre functions and the method of random Bregman projections[END_REF].

In [START_REF] Butnariu | Iterative methods of solving stochastic convex feasibility problems and applications[END_REF], we proposed the following algorithm to solve (4).

(46) fix x 0 ∈ U and set (∀n ∈ N)

In view of (3) and Definition 3.7, the sequences (x n ) n∈N and (y n ) n∈N are well-defined. We now study the asymptotic behavior of this algorithm, starting with two key monotonicity properties.

Proposition 4.1 Let (x n , y n ) n∈N be generated by (46). Then:

Proof. This is a direct consequence of Definition 3.7 and (46). Proposition 4.2 Let (x n , y n ) n∈N be generated by (46) and take {x, y} ⊂ U . Then:

Proof. Fix n ∈ N. If x ∈ dom ϕ or y ∈ dom ψ, then (48) is clear. Otherwise, it follows from Lemma 2.8, Lemma 2.6, Proposition 3.10 that

Hence

Corollary 4.3 Let (x n , y n ) n∈N be generated by (46) and suppose that p in (5) is finite. Then

Proof. In view of Proposition 4.1, let λ = lim Λ(x n , y n ) = lim Λ(x n+1 , y n ). Clearly, λ ∈ [p, +∞[. Let us assume that λ > p and we shall obtain a contradiction. Take {x, y} ⊂ U and ε ∈ ]0, +∞[ such that λ = Λ(x, y) + ε. Then Proposition 4.2 yields

Our main convergence result can now be stated and proved.

Theorem 4.4 Let (x n , y n ) n∈N be a sequence generated by algorithm (46) and suppose that S is nonempty (and hence p in (5) is finite). Then

Moreover, (x n , y n ) n∈N converges to a point in S.

Proof. Take (x, y) ∈ S. It follows from (48) that

Therefore, (53) holds. Moreover, (54) and Proposition 3.16 imply that (55) (x n ) n∈N is Bregman monotone with respect to E ⊂ U.

In view of Lemma 2.10 (with θ = 0) and A5, the sequence (x n ) n∈N is bounded and all its cluster points lie in U . Let us consider a convergent subsequence, say (56) Hence ( x, y) ∈ S and thus x ∈ E by Proposition 3.16(ii)&(iii). Therefore, every cluster point of (x n ) n∈N belongs to E. Consequently, utilizing (55) and Lemma 2.11, we conclude that (x n ) n∈N converges a point in E, say x. Set ȳ = --→ prox ψ (x). Proposition 3.16(iv) shows that (x, ȳ) ∈ S. On the other hand, Proposition 3.10 implies that

Let us illustrate Theorem 4.4 by presenting some immediate applications; further examples will be provided in Section 5. converges to a point in S.

Proof. This is a consequence of Example 3.9 and Theorem 4.4.

Remark 4.8 Corollary 4.7 corresponds to Csiszár and Tusnády's classical alternating minimization procedure (see their seminal work [START_REF] Csiszár | Information geometry and alternating minimization procedures[END_REF]) which, in turn, covers the expectation-maximization method for a specific Poisson model [START_REF] Iusem | A short convergence proof of the EM algorithm for a specific Poisson model[END_REF]. For an alternative proof of Corollary 4.7 when A ∩ B = ∅, see [7, Application 5.5].

Corollary 4.9 Take θ ∈ Γ 0 (X) and suppose that its set M of minimizers over U is nonempty. Remark 4.10 Item (i) in Corollary 4.9 goes back to [START_REF] Censor | Proximal minimization algorithm with D-functions[END_REF]. A special case of item (ii) in the context of the Kullback-Leibler divergence appears in [START_REF] Eggermont | Multiplicative iterative algorithms for convex programming[END_REF], see also [START_REF] Eggermont | On EM-like algorithms for minimum distance estimation[END_REF]Remark 2.18]. If f = 1 2 • 2 , then items (i) and (ii) reduce to a classical result of Martinet [START_REF] Martinet | Détermination approchée d'un point fixe d'une application Cas de l'application prox[END_REF].

The next two statements concern the invariance of the solution set S. Proof. Consider the iteration (46) with starting point x 0 = x. Using Proposition 3.16, we see that x n ≡ x and y n ≡ y. Hence (53) yields D D f (x, y), ( x, y) = 0. Proposition 5.3 Suppose that the solution set G of (76) is nonempty and let (y n ) n∈N be a sequence generated by

Then (y n ) n∈N converges to a point in G.

The above result can be applied to the problem of finding relaxed solutions to the inequality problem (70) by choosing (θ i ) i∈I as in (69). This approach is of special interest when (70) has no solution since the standard subgradient projection techniques that are available to solve this problem [START_REF] Bauschke | On projection algorithms for solving convex feasibility problems[END_REF][START_REF] Combettes | Convex set theoretic image recovery by extrapolated iterations of parallel subgradient projections[END_REF][START_REF] Kiwiel | Surrogate projection methods for finding fixed points of firmly nonexpansive mappings[END_REF] all fail in this situation. In the particular case of the convex feasibility problem (71), the relaxed problem ( 76 Then (y n ) n∈N converges to a point in G.

Right Bregman relaxation

The left relaxation techniques developed in Section 5.2 have natural right counterparts. Since the resulting statements have largely similar proofs, we shall only highlight the main aspects of this approach.

The right relaxation of (66) with respect θ is obtained by setting ϕ = ι C and ψ = θ in (39), which yields (85) minimize -→ env θ over C ∩ U.

We derive at once from Theorem 4.4 and Proposition 3.16 the following result.

Proposition 5.5 Suppose that (21) holds and that the solution set E of (85) is nonempty. Then the sequence (x n ) n∈N generated by

converges to a point in E.

Now assume that the functions (θ i ) i∈I in (67) satisfy ( 21) (hence so does θ). Then a right relaxation of ( 68) is (87) minimize

The next two results are the right counterparts of Propositions 5.2 and 5.3.

Proposition 5.6 Let (x, y) ∈ U 2 . Then: Then (x n ) n∈N converges to a point in G.