Marin Fouchier
email: marin.fouchier@fresnel.fr

Myriam Zerrad

Michel Lequime

Claude Amra

A breakthrough instrument for wide-range wavelength and angle resolved light scattering measurements of high performances optical coatings

Keywords: light scattering, optical densities, metrology, optical thin-films

Due to the diversity and complexity of optical functions they address, optical coatings are widely used in space optics but the level of both their specifications and performances is now so high that the characterization of their optical properties requires the development of extreme metrology. In this context we present in this paper the last state-of-the-art instrument for thin-films characterization. This apparatus is unique by the performances and the levels of detectivity achieved. It operates spectrally and angularly resolved characterization of Angle Resolved Scattering (ARS) over a wide and continuous spectral range from 400 nm to 1.7 μm. The spectral resolution is lower than 1 nm, while the angular resolution is lower than 1° with a precision better than 1%. The lowest levels of scattering measured are close to 10 -8 sr -1 , near to the limit of air particles scattering. The large angular range available offers the opportunity to evaluate scattering functions both in the forward and backward half-spaces. In specific configurations, this same setup records specular properties, achieving again unreached levels of detection and precision over the same spectral range and with the same resolution. All measurements have shown great agreement with theoretical designs and numerical predictions, for both scattering and transmittance records. This instrument is an upgraded version of an apparatus developed in the Scattering group for more than five years formerly operating in the visible wavelength spectrum and now extended to near-infrared with the same level of performances.

INTRODUCTION

Optical thin-film components are widely used for space applications from Earth imaging [1] to green house gas monitoring [START_REF] Pasternak | The microcarb instrument[END_REF] and gravitational waves detection [START_REF] Audley | Laser Interferometer Space Antenna[END_REF]. The levels of both their specifications and performances are so high that they require very strong and accurate metrology. Light scattering phenomena severely affect components performances by reducing or distorting light signal or altering image quality. Accurate and performant metrology tools are needed to qualify these space dedicated components. SALSA (Spectral and Angular Light Scattering Apparatus) is a spectrophotometric instrument, spectrally and angularly resolved, which performs such a characterization. The instrument emerges from a more than five years development process. Previous publications have presented its operation in the visible range (400 -1000 nm) for both scattered light and specular functions estimation [START_REF] Liukaityte | Broadband spectral transmittance measurements of complex thin-film filters with optical densities of up to 12[END_REF]- [START_REF] Fouchier | Wide-range wavelength and angle resolved light scattering measurement setup[END_REF]. This paper presents an upgraded version of the SALSA instrument carrying out measurements in the visible and the near-infrared spectra, from 400 to 1700 nm. Its detection limit is closed to the air particle scattering level ( ) while it measures up to 13 optical densities (ODs) on specular functions.

Many instruments have been proposed for this metrology but most of them use laser lines to reach sufficient power and measure scattering at discrete wavelengths in the UV, visible or infrared [START_REF] Finck | International round-robin experiment for angle-resolved light scattering measurement[END_REF], [START_REF] Hauptvogel | Light scattering characterization of optical components for space applications[END_REF]. The configuration of the setup presented in this paper overcomes this limitation and owns a continuous spectral range.

ICSO 2020 International Conference on Space Optics

Virtual Conference 30 March-2 April 2021

PRESENTATION OF SALSA INSTRUMENT

Design of the instrument

SALSA sits in an ISO6 clean room facility that preserves it from general environmental variations and dusts. Figure 1 shows a functional representation of the SALSA instrument composed of different subassemblies. The emission subassembly is based on a supercontinuum laser source (FIANIUM WL-SC-400-8) which outputs a very wide spectrum light signal from 400 nm to 2.4 μm with a total integrated power of 8 W. A dichroic beam splitter (FIANIUM Splitter-950-01) splits the broadband spectrum into two channels: a visible channel (from 400 to 975 nm) and an infrared channel (from 975 nm to 2.4 μm).

The optical components used on both channels are functionally identical but differ on their spectral operation ranges. The two channels are used independently and asynchronously, if one is activated then the other one stay closed. Two tunable volume hologram filters (Photon Etc. LLTF Contrast VIS HP8 and LLTF Contrast SWIR HP8 respectively) set the working wavelength and reduce the spectral width of the signal between 2 and 3.5 nm. Electro-mechanical shutters (SH1 and SH2) block the light beam and transmit it only during acquisitions. Order sorting filters (OSF1 and OSF2) clean the signal from parasitic harmonic orders induced by the volume hologram filters. Mirrors (M1, M2 and M3) spatially align the optical beams and beam splitters (BS1 and BS2) collect a small portion of the signal, around 10%, as reference signals. Then the main signal goes through a set of four neutral optical density stages (ODs). The combination of these densities results in a total OD between 0 and 10, adjustable by a pitch of 1. The only structural difference between the visible channel and the infrared channel is a removable recombination mirror (RM) which reflects the visible channel towards the optical densities. This mirror is mounted on a translation stage, when performing acquisitions in the visible spectrum the mirror reflects the light beam towards the output of the emission subassembly, while it clears the way for the infrared channel when performing infrared acquisitions. The light beam is finally collected by a parabolic off-axis reflective collimator (ORC).

ICSO 2020 International Conference on Space Optics

Virtual Conference 30 March-2 April 2021

The measurement subassembly is connected to the emission subassembly with a multimode fiber optical link (FOL1, core diameter 50 μm). A reflective collimator (TRC) outputs the optical beam which illuminates the sample under characterization (TFF). The incidence angle is set by a rotation stage that holds the sample. After being transmitted, reflected, or scattered, light is collected by a second reflective collimator (RRC) mounted on a rotation arm centered on the sample. Therefore, by adjusting the angular position of the mechanical arm, it is possible to measure scattered light at any angle around the sample and in the plane of incidence. The measurement subassembly is connected to the detection subassembly with a multimode fiber optical link (FOL2, core diameter 105 μm).

The detection subassembly includes a Czerny-Turner monochromator (ANDOR Shamrock SR-193i-B1-SIL), two gratings are mounted back-to-back on a turret in the spectrograph. The first one is an 1800-lines-per-milimeters holographic grating optimized for a range from 360 to 1100 nm. The second one is a 760-lines-per-milimeters optimized for 850 nm to 2 μm. The spectrograph owns two output ports, on each is mounted a scientific grade camera: a CCD detector on port1 (ANDOR Newton DU920P-BVF) with 255 x 1024 pixels of 26 x 26 μm each and an InGaAs linear detector on port2 (ANDOR DU491A) with 1 x 1024 pixels of 500 x 25 μm each. A sliver-coated flip mirror (FM) leads the light beam to one of the two ports. The CCD camera detects visible light, while InGaAs is sensitive to the infrared spectrum.

Operation of the instrument

The measurement procedure has been already presented in our previous papers [START_REF] Lequime | Breakthrough spectrophotometric instrument for the ultra-fine characterization of the spectral transmittance of thin-film optical filters[END_REF]- [START_REF] Fouchier | Wide-range wavelength and angle resolved light scattering measurement setup[END_REF]. For measuring the amount of scattered light by any sample at a wavelength and an angle , the rotating arm is set at the chosen angle. The value of determines which channel (visible or infrared) to use in the emission and detection subassemblies. The involved LLTF filter, the recombination mirror, the spectrometer grating, and the detection camera are tuned to the chosen value. Initial values are set for the (exposition time) and (overall optical densities) parameters. Then the electromechanical shutter is opened to let the light beam illuminate the sample. The acquisitions on the reference photodiode and the camera are synchronized thanks to this single shutter. The signal on the camera is expected to be around 70% of its saturation level. The effective level is adjusted by tuning the and parameters simultaneously. Several acquisitions might be executed until the signal achieves a convenient level. The CDD camera has numerous lines of pixels and the light beam illuminates only a small area on the detector. The illuminated area is small enough to acquire simultaneously the signal and the noise using different zones of the detector matrix. On the other hand, the InGaAs camera has only one line of pixels which collect the signal. transmittance and the reflectance. In any case a calibration is required. A Spectralon white diffuse reflectance sample defines the standard levels for the scattering measurements, while a baseline (with no sample) acquisition defines the reference level for specular functions measurements.

SAMPLE CHARACTERIZATION AND COMPARISON TO MODELS

Light scattering characterization

A dedicated bandpass filter centered at 1273 nm has been used for the qualification of the instrument, this sample is called F0 in this paper. This filter has been designed and manufactured by CNES and CILAS ArianeGroup and is intended for a space application. It is composed of around 100 alternating layers of SiO2 and Nb2O5. The calibration of the setup requires the measurement of the ARS of a perfect white reflectance Lambertian diffuser (LABSPHERE Spectralon). The ARS analytical formulation for this material is the following:

Where is the albedo of the sample at wavelength , the angular cosine dependance is a specific feature used to check the correct alignment of the setup. The ARS measured on this sample is compared to a cosine shape, a good agreement assesses a right optical alignment of the sample and the rotating arm. The angular resolution is granted by the constitution of the setup and the use of a high-precision rotation stage. The angular response presents a relative decrease of 8 decades at 1° from the specular beam position. Therefore, we can confidently assert that the angular resolution of the setup is below 1°. ARS measurements on the sample F0 show both spectral and angular variations. Recent works carried out in our group have led to a complete and powerful theoretical model implementation which describe such scattering from optical thinfilms [START_REF] Amra | Electromagnetic Optics of Thin-Film Coatings: Light Scattering, Giant Field Enhancement, and Planar Microcavities[END_REF]. The comparison between the model and the measurements presents an excellent and unprecedented agreement. Not to forget, the sample under characterization is very complex which makes it difficult to characterize experimentally and numerically. The model also highlights the spectral and angular peaks observed experimentally. An acquisition without sample depicts the Signature of the instrument, i.e., the lowest detectable values. The level of the Signature shows that the setup should be able to measure down to , closed to the air particles scattering limit. F0 has been characterized at some specific wavelengths, Figure 5 and Figure 6 depict the measurement and model data at 1115 nm and 1270 nm with a fine angular resolution (1 nm). The same comparison has been carried along all the near-infrared accessible spectrum with a lower angular resolution and the agreement with the model remains excellent.

Specular transmittance characterization

The same F0 filter has been subjected to specular transmittance characterization, illustrated on Figure 7. The measurement data shows an excellent agreement with the design data: the high rejection levels are well defined, as well as the narrow spectral peaks and the steep edges. These results validate the asserted that the spectral resolution is lower than 1nm. Once again the Signature of the instrument is very low, around , which is 6 decades below the state-of-the-art [START_REF] Johansen | THIN-FILM INTERFERENCE FILTERS FOR LIDAR[END_REF]. The rejection of the above-mentioned optical filter is considerably larger than the noise level. To validate the detection limits, a 60 nm thick silver monolayer neutral density, called Ag2, is superimposed on F0. This optical density has a transmittance around in the near infrared range. The total transmittance of the two components superposition is measured between 1025 nm and 1275 nm, roughly centered on one of the two main deep rejection bands. The resulting transmittance value is the one of the F0 filter lowered by a factor of . The lowest level achieved in this configuration is , close to the Signature limit.

Improve the spectral resolution

The measurement campaign for F0 characterization over the whole available spectrum counts 1301 data points, from 400 nm to 1.7 μm by step of 1 nm. The resulting resolution is 1 nm, but it can be improved with an adequate data postprocessing. The combination of a monochromator with an array-based detector shows a great interest in terms of spectral discrimination. The light beam arriving at the monochromator has a small spectral bandwidth (between 2 and 3.5 nm) dispersed spatially by the grating. The signal acquired by the camera is spread over a horizontal line of pixels Considering the dispersion of the grating it is possible to discretize the wavelengths closed to the working one using different pixels of the camera. This configuration cleans the signal from non-desired spectral contributions and increases the spectral resolution. The instrument can consequently resolve narrow spectral peaks (<1 nm). Figure 9 illustrates this spectral dispersion on the raw signals detected by the cameras, the blue points represent baseline acquisitions (no sample), the red points represent acquisitions with sample. The presence of a sample under characterization modifies the signals and provides information about the spectral function locally. The spatial information along the pixel line is transcribed into spectral information. Figure 10 depicts a high-resolution processing, resulting in a transmittance resolution of 0.1 nm obtained from a 1 nm resolution acquisition. The comparison with the high-resolution processing is excellent and suggests that the narrow spectral peaks are even better resolved. The transmittance levels are displayed only between and on vertical axis and from 1075 nm to 1300 nm for ease of reading. processing. The high-resolution processing is coherent with the classical processing and improves the spectral information.

CONCLUSION

We have presented an angularly and spectrally resolved scatterometer that performs optical thin-film components characterization. The setup outperforms other existing instruments and offers very broad and continuous spectral and angular ranges. The reached detection limits are for scattered light and for transmittance measurements in the SWIR. This instrument developed over the last years is now operating in the visible as in the near-infrared and it is part of the DIFFUSIF metrological platform [START_REF]Diffusif[END_REF]. The performances have been qualified with the characterization of a dedicated filter intended for space applications.

  The reference signals collected by the above-mentioned beam splitters (BS1 and BS2) are also sent to reflective collimators (RC R1 and RC R2) connected to photodiodes (VIS_PD and IR_PD, respectively FEMTO DE-OE-200-SI and DE-OE-200-In1) via multimode optical fibers (FOL R1 and FOL R2).

Figure 1

 1 Figure 1 Spectral and Angular Light Scattering Apparatus (SALSA) is composed of three subassemblies: Emission, Measurement and Detection.

Figure 2

 2 Figure 2 SALSA instrument in ISO6 clean room in the Espace Photonique of Institut Fresnel. Measurement subassembly: a reflective collimator (right) outputs the optical beam which illuminates the sample (middle), the scattered light is collected by a second reflective collimator mounted on a rotating arm (left).

  It is not possible to acquire the signal and the noise simultaneously. The noise level is measured immediately after the signal with another acquisition of the same duration but with no illumination of the sample. Depending on the position of the rotating arm, SALSA measures either the scattered light or the specular transmitted or reflected light. The scattering function estimated is the Angle Resolved Scattering (ARS), a spectrally and angularly resolved function which describes the angular distribution of the scattered light around the sample. The specular functions estimated are the

Figure 3

 3 Figure 3 Comparison between the ARS cosine model and the normalized ARS measured on the Spectralon sample at 1490 nm.

Figure 4

 4 Figure 4 Angular response of the instrument to the direction defined by the transmitted specular beam.

Figure 5

 5 Figure 5 ARS measured on F0 at 1115 nm. Comparison with Model and Signature.

Figure 6

 6 Figure 6 ARS measured on F0 at 1270 nm. Comparison with Model and Signature.

Figure 7

 7 Figure 7 Transmittance of the dedicated F0 filter from 400 to 1650 nm. Comparison between the design, the measurement, and the Signature of the instrument.

Figure 8

 8 Figure 8 Transmittance of F0 superposed with a silver monolayer density (Ag2) compared with single F0 and single Ag2 on the [1025-1275 nm] band. The resulting transmittance is close to the Signature level.

Figure 9 Figure 10

 910 Figure 9 Raw signals acquired on the 1024 pixels line. The spatial position corresponds to a spectral information around the working wavelength (left 850 nm, right 1350 nm).

Conference on Space Optics Virtual Conference 30 March-2 April 2021

	ICSO 2020
	International
				Design	Measurement	Signature
		1,0E+00
		1,0E-02
		1,0E-04
	tansmittance	1,0E-08 1,0E-06
		1,0E-10
		1,0E-12
		1,0E-14
			400	600	800	1000	1200	1400	1600
					wavelength (nm)
				F0	Ag2	F0 & Ag2	Signature
			1,0E+00
			1,0E-02
		transmittance	1,0E-10 1,0E-08 1,0E-06 1,0E-04
			1,0E-12
			1,0E-14
			1025	1075	1125	1175		1225	1275
					wavelength (nm)

Proc. of SPIE Vol. 11852 118521C-2 Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 12 Jul 2021 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use

Proc. of SPIE Vol. 11852 118521C-3 Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 12 Jul 2021 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use

Proc. of SPIE Vol. 11852 118521C-4 Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 12 Jul 2021 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use

Proc. of SPIE Vol. 11852 118521C-5 Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 12 Jul 2021 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use

Proc. of SPIE Vol. 11852 118521C-6 Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 12 Jul 2021 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use

Proc. of SPIE Vol. 11852 118521C-7 Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 12 Jul 2021 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use

Proc. of SPIE Vol. 11852 118521C-8 Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 12 Jul 2021 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use

Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 12 Jul 2021 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use

Proc. of SPIE Vol. 11852 118521C-10 Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 12 Jul 2021 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use

ACKNOWLEDGEMENT

The authors would like to thank CNES (Centre National d'Etudes Spatiales), CILAS ArianeGroup, DGA (Direction Générale de l'Armement) and LabTOP (joint laboratory between CILAS ArianeGroup and Institut Fresnel) for their contribution to this work.


B João

Carlos De O Souza

João S Andrade

Jurandir De O Lopes
email: jurandir@ufpi.edu.br

João Carlos

O Souza

An inertial proximal point method for difference of maximal monotone vector fields in Hadamard manifolds

Keywords: Variational inclusion, Proximal point method, Monotone vector fields, DC functions, Hadamard manifolds Mathematics Subject Classification 65K05, 90C26, 90C48

We propose an inertial proximal point method for variational inclusion involving difference of two maximal monotone vector fields in Hadamard manifolds. We prove that if the sequence generated by the method is bounded, then every cluster point is a solution of the non-monotone variational inclusion. Some sufficient conditions for boundedness and full convergence of the sequence are presented. The efficiency of the method is verified by numerical experiments comparing its performance with classical versions of the method for monotone and non-monotone problems.

Introduction

The problem of finding a zero of maximal monotone operators T , i.e., 0 ∈ T (x) [START_REF] Ansari | Existence and boundedness of solutions to inclusion problems for maximal monotone vector fields in Hadamard manifolds[END_REF] includes, as special cases, optimization and min-max problems, complementarity problems, and variational inequalities. It finds many important applications in scientific fields such as image processing, computer vision, machine learning and signal processing. For this reason, in recent years much attention has been given to develop efficient and implementable numerical methods for solving this problem in different contexts; see for instance [START_REF] Ansari | Existence and boundedness of solutions to inclusion problems for maximal monotone vector fields in Hadamard manifolds[END_REF][START_REF] Ansari | Inexact proximal point algorithms for inclusion problems on Hadamard manifolds[END_REF][START_REF] Ansari | Proximal point algorithm for inclusion problems in Hadamard manifolds with applications[END_REF] and references therein. One of the fundamental approaches for solving [START_REF] Ansari | Existence and boundedness of solutions to inclusion problems for maximal monotone vector fields in Hadamard manifolds[END_REF] is the proximal point method which its origin can be traced back to Martinet [4] in the context of convex minimization and Rockafellar [START_REF] Rockafellar | Characterization of the subdifferentials of convex functions[END_REF] in the general setting of maximal monotone operators in Hilbert space. In the Riemannian setting this method was studied by Li et al. [START_REF] Li | Monotone vector fields and the proximal point algorithm on Hadamard manifolds[END_REF] based on first extension of the proximal point method to the Riemannian context proposed by Ferreira and Oliveira [START_REF] Ferreira | Proximal point algorithm on Riemannian manifolds[END_REF] for convex minimization problems.

The proximal point method computes at each iteration the well known resolvent operator introduced by Moreau [START_REF] Moreau | Proximité et dualité dans un espace Hilbertien[END_REF]. Unfortunately, in many interesting cases, the evaluation of the resolvent operator is as difficult as solving the original problem. On the other hand, in many problems, the operator T can be written as the sum of two maximal monotone operators, namely T = A + B, such that the resolvent operator of each component is much easier to compute than the original operator T . The so-called forward-backward methods overcome this drawback combining the resolvents of each component to find a solution of the original problem; for instance the Douglas-Rachford algorithm [START_REF] Douglas | On the numerical solution of heat conduction problems in two and three space variables[END_REF] among others.

It is well known that the sum of two monotone operators is a monotone operator, whereas the difference of two monotone operators is not necessarily a monotone operator. Therefore, the problem of finding a zero of the difference of two monotone operators can be very difficult. It generalizes the problem of finding the critical points of the difference of two convex functions (DC functions) and was not studied extensively yet (even in the linear setting); see for instance [START_REF] Moudafi | On the difference of two maximal monotone operators: regularization and algorithmic approach[END_REF][START_REF] Alimohammady | Notes on the difference of two monotone operators[END_REF][START_REF] Alimohammady | Inertial proximal algorithm for difference of two maximal monotone operators[END_REF][START_REF] Moudafi | On critical points of the difference of two maximal monotone operators[END_REF][START_REF] Noor | On difference of two monotone operators[END_REF]. In this direction, Souza and Oliveira [START_REF] Souza | A proximal point method for DC functions on Hadamard manifolds[END_REF] proposed the first extension of the proximal point method for DC functions in Hadamard manifolds.

On the other hand, Polyak [START_REF] Polyak | Some methods of speeding up the convergence of iteration methods[END_REF] introduced the so-called heavy ball method for minimizing a smooth convex function. The difference compared to the gradient method is that each iteration, an extrapolation point (which combines the current and the previous iterates) is used instead of the current iterate. This minor change improves the performance of the gradient method. Alvarez and Attouch [START_REF] Alvarez | An inertial proximal method for maximal monotone operators via discretization of a nonlinear oscillator with damping[END_REF] adapted this idea (called inertial method) to the proximal point method for maximal monotone operators. Inertial methods were considered in the particular case of DC functions in the linear context; see for instance [START_REF] Maingé | Convergence of new inertial proximal methods for DC programming[END_REF][START_REF] Oliveira | Level an inertial algorithm for DC programming[END_REF]. However, to the best of our knowledge, there does not exist any study of inertial proximal method for variational inclusion involving the difference two monotone vector fields in Hadamard manifolds and, in particular, in the linear setting. The present work is a contribution towards this goal.

The aim of this paper is to proposed an inertial version of the proximal point method for the variational inclusion problem 0 ∈ A(x) -B(x), [START_REF] Ansari | Inexact proximal point algorithms for inclusion problems on Hadamard manifolds[END_REF] where A and B are maximal monotone vector fields on finite dimensional Hadamard manifold M. It is equivalent to the problem of finding x ∈ M such that

A(x) ∩ B(x) = ∅. (3)
We illustrate with some preliminaries numerical experiments that the inertial version of the proximal point method has a better performance compared to its classical version.

Although we propose an inertial version of the proximal point method for variational inclusion involving the difference of two maximal monotone vector fields in Hadamard manifolds, some results related to the proposed method are new for DC functions in the linear setting as well as in some particular instances in Hadamard manifolds, namely, if B(x) = 0 in [START_REF] Ansari | Inexact proximal point algorithms for inclusion problems on Hadamard manifolds[END_REF] we have an inertial version of the proximal point method considered by Li et al. [START_REF] Li | Monotone vector fields and the proximal point algorithm on Hadamard manifolds[END_REF] for finding a singularity of a maximal monotone vector field; if A(x) = ∂ g(x) and B(x) = ∂h(x) in [START_REF] Ansari | Proximal point algorithm for inclusion problems in Hadamard manifolds with applications[END_REF], where g, h are convex functions, then our method becomes an inertial version of the proximal point method for DC functions ( f (x) = g(x)h(x)) proposed by Souza and Oliveira [START_REF] Souza | A proximal point method for DC functions on Hadamard manifolds[END_REF]; if A(x) = ∂ f (x) and B(x) = 0 in [START_REF] Ansari | Inexact proximal point algorithms for inclusion problems on Hadamard manifolds[END_REF], where f is a convex function, then we have an inertial version of the proximal point method proposed by Ferreira and Oliveira [START_REF] Ferreira | Proximal point algorithm on Riemannian manifolds[END_REF] for convex minimization problems; and finally, if the parameter which appears in the direction of the inertial method is zero, we have the (non-boosted) proximal point method for difference of maximal monotone operators extending the algorithm proposed by Souza and Oliveira [START_REF] Souza | A proximal point method for DC functions on Hadamard manifolds[END_REF].

In our convergence analysis, we prove that every cluster point of the sequence generated by the inertial method, if any, is a solution of the variational inclusion problem. Furthermore, we present some sufficient conditions for boundedness and full convergence of the method which are new even for DC functions in Hadamard manifolds. To show the efficiency of the method, we provide some numerical experiments solving minimization problems in a genuine (with non-constant curvature different from zero) Hadamard manifold involving convex and non-convex functions as well as variational inclusion involving operators which are not the subdifferential of convex functions in Euclidean space.

The study of algorithms for solving non-monotone variational inclusion is a very difficult problem and interesting by itself. The variational inclusion problem involving the difference of maximal monotone vector fields generalizes the very important problem of minimizing DC functions which has a lot of applications. Another interesting application of difference of monotone operators in linear spaces is given by Attouch and Théra [START_REF] Attouch | A general duality principle for the sum of two operators[END_REF]. They mention that most of the equations arising in physics, economics, among other, can be written as follows

0 ∈ C(x) = A(x) + B(x), ( 4
)
where A, B are possibly multivalued operators. The splitting of the operator C into the sum of two operators A and B has usually a deep physical or economical meaning since A and B may have very distinct properties. Instead of directly study (4), the authors in [START_REF] Attouch | A general duality principle for the sum of two operators[END_REF] proposed an equivalent duality transformation and they consider the following difference of multivalued operators

0 ∈ A -1 (y) -B -1 (-y)
offering a new duality approach to some central questions in the theory of variational inequalities and maximal monotone operators. The remainder of this paper is organized as follows. In Sect. 2, we recall some notations, definitions and preliminary results in Riemannian manifolds, convexity and vector fields in Hadamard manifolds which will be used for further analysis. In Sect. 3, we present our algorithm as well as its well definition. In Sect. 4, we establish convergence analysis of the proposed algorithm under some mild conditions for the non-monotone case. In Sect. 5, we present the convergence analysis of the method for the monotone case. In Sect. 6 some numerical experiments are reported to support the theoretical results obtained and illustrate the feasibility and efficiency of the proposed algorithm comparing its performance with the classical version of the proximal point method in both Euclidean space and Hadamard manifolds. Some concluding remarks are given in the last section.

Notation and basic concepts

The standard notations, results and preliminary concepts of Riemannian geometry used throughout the paper can be found, for instance, in Sakai [START_REF] Sakai | Riemannian geometry. Translations of mathematical monographs[END_REF], Udriste [START_REF] Udriste | Convex functions and optimization algorithms on Riemannian manifolds[END_REF] and do Carmo [START_REF] Do Carmo | Riemannian geometry[END_REF]. We follow the notation, terminology and results of [START_REF] Souza | A proximal point method for DC functions on Hadamard manifolds[END_REF] and [START_REF] Almeida | A modified proximal point method for DC functions on Hadamard manifolds[END_REF].

Throughout this paper, we will assume that M a finite dimensional Hadamard manifold. We denote by T x M the tangent space of M at x. Recall that the parallel transport along the geodesic γ from

P γ,γ (b 2 ),γ (b 1 ) • P γ,γ (b 1 ),γ (a) = P γ,γ (b 2 ),γ (a) and P -1 γ,γ (b),γ (a) = P γ,γ (a),γ (b) .
If γ (a) = p and γ (b) = q, we will write P q, p instead of P γ,γ (b),γ (a) in the case when γ a minimal geodesic joining p to q and no confusion arises. The restriction of a geodesic to a closed bounded interval is called a geodesic segment. Given points x, y ∈ M, we denote the geodesic segment from x to y by [x, y]. We usually do not distinguish between a geodesic and its geodesic segment, as no confusion can arise. We denote by exp x : T x M → M the exponential map. For any x ∈ M we can define the exponential inverse mapping exp -1 x :

M → T x M which is C ∞ . Since d(x, x ) = exp -1 x (x) , then the map ρ x : M → R defined by ρ x (x) = 1 2 d 2 (x, x
) is C ∞ and its gradient at x, denoted by gradρ x (x), is given by gradρ x (x) =exp -1

x (x ); see for instance [START_REF] Ferreira | Proximal point algorithm on Riemannian manifolds[END_REF]Proposition 3.3]. Using the properties of the parallel transport and the exponential map, we obtain the following proposition that will be used in the next sections.

Proposition 1 Let M be a Hadamard manifold. Let x ∈ M and {x k } ⊂ M be such that x k → x. Then the following assertions hold.

For any y ∈ M, we have

exp -1 x k y → exp -1 x y and exp -1

y x k → exp -1 y x. 2. If v k ∈ T x k M and v k → v, then v ∈ T x M. 3. Given u k , v k ∈ T x k M and u, v ∈ T x M, if u k → u and v k → v, then u k , v k → u, v . 4. For any u ∈ T x M, the function F : M → T M defined by F(z) = P z,x u for each z ∈ M is continuous on M. Proof See [5, Lemma 2.4].
Let X (M) denote the set of all multivalued vector fields A : M -→ 2 T M such that A(x) ⊂ T x M, for each x ∈ M, T M = ∪ x∈M T x M and the domain of A is closed and convex, where the domain D(A) of A is defined by

D(A) = {x ∈ M; A(x) = ∅}. Definition 1 Let A ∈ X (M).
Then A is said to be: (i) Monotone if the following condition holds for any x, y ∈ D(A):

u, exp -1 x y ≤ v, -exp -1 y x , ∀u ∈ A(x), v ∈ A(y); (5)
(ii) Strictly monotone if (5) holds with strict inequality for any x, y ∈ D(A) with x = y, that is, u, exp -1 x y < v,exp -1 y x , ∀u ∈ A(x), v ∈ A(y); (iii) Strongly monotone if there exists ρ > 0 such that, for any x, y ∈ D(A), we have u, exp -1

x y -v,exp -1 y x ≤ -ρd 2 (x, y), ∀u ∈ A(x), v ∈ A(y); (iv) Maximal monotone if it is monotone and the following implication holds for any x ∈ D(A) and u ∈ T x M: Definition 3 Given λ > 0 and A ∈ X (M), the resolvent of order λ is the set-valued mapping J A λ : M -→ 2 M defined by [START_REF] Li | Resolvents of set valued monotone vector fields in Hadamard manifolds[END_REF]Remark 3] by the definition of the resolvent of a vector field, the domain of the resolvent J A λ is the range of the vector field defined by x → exp x λA(x). We will denote this range as R(exp (•) λA(•)). Then we have that

u, exp -1 x y ≤ v, -exp -1 y x , ∀y ∈ D(A), v ∈ A(y) ⇒ u ∈ A(x). Definition 2 Let A ∈ X (M)
J A λ (x) = {z ∈ M; x ∈ exp z λA(z)}, ∀x ∈ M. Remark 2 As mentioned in
D(J A λ ) = R(exp (•) λA(•)
). Definition 4 Let K ⊆ M be a non-empty, closed and convex set. Given a mapping T : K ⊆ M → M, we say that T is firmly non-expansive if for any x, y ∈ K , the function

: [0, 1] → R + defined by (t) = d(exp x t exp -1 x T x, exp y t exp -1 y T y), ∀t ∈ [0, 1], is non-increasing.
Next, we state a well known result on firmly non-expansive mappings.

Proposition 3 Let K ⊆ M be a non-empty, closed and convex set and T : K ⊆ M → M. The following assertions are equivalent. 1. T is firmly non-expansive; 2. For any x, y ∈ K and t ∈ [0, 1], one has d(T (x), T (y)) ≤ d(exp x t exp -1

x T x, exp y t exp -1 y T y).

Proof See [START_REF] Li | Resolvents of set valued monotone vector fields in Hadamard manifolds[END_REF]Proposition 5].

The next result establishes the relation between the firm non-expansivity of the resolvent and the monotonicity of the corresponding vector field. Proposition 4 Let A ∈ X (M). The following assertions hold for any λ > 0.

(i) The vector field A is monotone if and only if J A λ is single-valued and firmly non-expansive. (ii) if D(A) = M, the vector field A is maximal monotone if and only if J A

λ is single-valued, firmly non-expansive and the domain D(J A λ ) = M.

Proof See [START_REF] Li | Resolvents of set valued monotone vector fields in Hadamard manifolds[END_REF]Theorem 4].

From Proposition 4 and Remark 2, we have the following result which constitutes a counterpart in the setting of Hadamard manifolds of the well known Minty's theorem.

Corollary 1 Let A ∈ X (M) be monotone such that D(A) = M, and let λ > 0. Then A is maximal monotone if and only if R(ex p (•) λA(•)) = M.

An inertial proximal point method

Throughout this paper, we assume that M is a finite dimensional Hadamard manifold and A, B ∈ X (M) are maximal monotone vector fields with

D(A) = D(B) = M.
We are interested in solving the following problem:

find x * ∈ M such that A(x * ) ∩ B(x * ) = ∅. (6)
A point x ∈ M satisfying ( 6) is said to be the critical point of the difference A -B. The set of critical points of A -B is defined by

S = {x ∈ M; A(x) ∩ B(x) = ∅}.
It is well known that the sum of two monotone vector fields remains monotone. However, the difference of two monotone vector fields is possibly non-monotone even in the Euclidean setting as we can see in the following simple example.

Example 1 Let g, h : R → R be convex functions given by g(x) = 1 4 x 4 and h(x) = 1 2 x 2 . Recall that the subdifferential of a convex function is a monotone operator, i.e., A(x) = ∇g(x) and B(x) = ∇h(x) are monotone operators. However, we can easily see that the function f given by f (x) = g(x)h(x) is not convex and its subdifferential is not monotone. On the other hand, the set of critical points of f satisfies [START_REF] Ferreira | Proximal point algorithm on Riemannian manifolds[END_REF].

In this paper, we suppose that S = ∅ and A or B is ρ-strongly monotone.

Remark 3

Although in our proof we need the strongly monotonicity of only one of the vector fields A and B, it is worth to mention that the strongly monotonicity of both A and B is not a restrictive assumption. In fact, if it does not hold for A and B, then we can obtain another decomposition satisfying this condition as follows

T (x) = A(x) -B(x) = [A(x) -ρ exp -1 x y] -[B(x) -ρ exp -1 x y],
for ρ > 0 arbitrary and y ∈ M fixed. One has that Ã

(x) = A(x) -ρ exp -1 x y and B(x) = B(x) -ρ exp -1
x y are ρ-strongly monotone; see [START_REF] Li | Monotone vector fields and the proximal point algorithm on Hadamard manifolds[END_REF]Remark 4.4].

Inertial proximal point method (IPPM)

Step 1: Given an initial point

x 0 ∈ M, γ k ∈ [0, ρ 2 )
and a bounded sequence of positive numbers {μ k } (to be specified latter). Define x -1 = x 0 .

Step 2: Given

x k ∈ M, define d k = γ k exp -1 x k x k-1 . Find w k ∈ B(x k ) and set y k = exp x k μ k (w k + d k ). ( 7
)
Step 3:

Compute x k+1 ∈ M such that 0 ∈ A(x k+1 ) - 1 μ k exp -1 x k+1 y k . ( 8
) If x k+1 = x k and d k = 0, stop. Otherwise, set k = k + 1 and return to Step 2.

Remark 4

The well definition of the sequences {y k } and {x k } directly follows from the fact the the exponential map is a global diffeomorphism in Hadamard manifolds and the vector field defined by

x → A(x) -1 μ k exp -1
x y k is strongly monotone, and hence, it has a unique singularity, respectively; see [START_REF] Li | Monotone vector fields and the proximal point algorithm on Hadamard manifolds[END_REF]Theorem 4.3 and Remark 4.4]. Moreover, from ( 7) and ( 8), we have

1 μ k exp -1 x k y k -d k ∈ B(x k ) and 1 μ k exp -1 x k+1 y k ∈ A(x k+1 ).
Thus, if x k+1 = x k and d k = 0, we obtain

1 μ k exp -1 x k y k ∈ B(x k ) and 1 μ k exp -1 x k y k ∈ A(x k ). Therefore, 1 μ k exp -1 x k y k ∈ A(x k )∩ B(x k ), i.e., A(x k )∩ B(x k ) = ∅ and hence x k is a solution of (6).
Remark 5 In the linear setting, Alimohammady and Ramazannejad [START_REF] Alimohammady | Inertial proximal algorithm for difference of two maximal monotone operators[END_REF] proposed an inertial proximal point algorithm for the difference of maximal monotone operators based on the method considered by Moudafi [START_REF] Moudafi | On critical points of the difference of two maximal monotone operators[END_REF]. These methods are based on the resolvent operator and the Yosida approximate. It is worth to mention that in the linear setting our method is different from the algorithms proposed in [START_REF] Moudafi | On the difference of two maximal monotone operators: regularization and algorithmic approach[END_REF][START_REF] Alimohammady | Inertial proximal algorithm for difference of two maximal monotone operators[END_REF][START_REF] Moudafi | On critical points of the difference of two maximal monotone operators[END_REF]. Note that if γ k = 0, for all k ∈ N in IPPM, then d k = 0 and (7) becomes

w k ∈ B(x k ) and y k = exp x k μ k w k . Since, x k+1 ∈ M is defined as (8), if B(x) = ∂h(x) and A(x) = ∂ g(x)
, where g, h : M → R are convex functions, then ( 8) is equivalent to

x k+1 = arg min x∈M {g(x) + 1 μ k d 2 (x, y k )}
which is the proximal point method for DC functions proposed by [START_REF] Souza | A proximal point method for DC functions on Hadamard manifolds[END_REF]. Therefore, IPPM with γ k = 0, for all k ∈ N, ca be viewed as an extension to difference of maximal monotone vector fields of the method considered in [START_REF] Souza | A proximal point method for DC functions on Hadamard manifolds[END_REF] for difference of convex functions. In this case, our method is different from the regularized methods for the difference of maximal monotone operators proposed by Moudafi [START_REF] Moudafi | On the difference of two maximal monotone operators: regularization and algorithmic approach[END_REF][START_REF] Moudafi | On critical points of the difference of two maximal monotone operators[END_REF] in the linear setting.

It is worth to mention that IPPM is new, in the Hadamard setting, even for solving monotone vector fields. Therefore, we consider a monotone version of IPPM by doing B(x) = 0. Then, IPPM becomes the following: IPPM: monotone version (mIPPM)

Step 1: Given an initial point

x 0 ∈ M, γ k ∈ [0, ρ
2 ) and a bounded sequence of positive numbers {μ k } (to be specified latter). Define

x -1 = x 0 . Step 2: Given x k ∈ M, define d k = γ k exp -1 x k x k-1 .
Find

y k = exp x k μ k d k . ( 9
)
Step 3:

Compute x k+1 ∈ M such that 0 ∈ A(x k+1 ) - 1 μ k exp -1 x k+1 y k . ( 10
) If x k+1 = x k and d k = 0, stop. Otherwise, set k = k + 1 and return to Step 2.
In the sequel, we present a convergence analysis for both algorithms IPPM for the difference of maximal monotone vector fields (which is a possibly non-monotone vector field) and mIPPM for finding a zero of a monotone vector field. To this end, we will consider their convergence analysis separately.

Convergence analysis: possibly non-monotone case

From now on, we consider {x k } the sequence generated by IPPM and we assume that x k+1 = x k , for all k ∈ N, otherwise the algorithm returns a solution of the problem. Now we shall establish its convergence properties. Proposition 5 Suppose that A or B is ρ-strongly monotone. Then,

ρ + 1 b d(x k+1 , x k ) ≤ P x k ,x k+1 u k+1 -v k + ρ 2 d(x k , x k-1 ), ( 11
)
for any v k ∈ A(x k ) and u k+1 ∈ B(x k+1 ).
Proof We will suppose that A is monotone and B is ρ-strongly monotone. The other case is analogous. By ( 7) e (8), we have

1 μ k exp -1 x k y k -d k ∈ B(x k ) and 1 μ k exp -1 x k+1 y k ∈ A(x k+1 ). Given x ∈ M, since B is ρ-strongly monotone, then 1 μ k exp -1 x k y k -d k , exp -1 x k x -u, -exp -1 x x k ≤ -ρd 2 (x, x k ), ∀ u ∈ B(x). ( 12
)
It follows from monotonicity of A that

1 μ k exp -1 x k+1 y k , exp -1 x k+1 x ≤ v, -exp -1 x x k+1 , ∀ v ∈ A(x). ( 13
)
Taking x = x k+1 in (12), x = x k in (13) and summing up these inequalities, we have

1 μ k [ exp -1 x k y k -d k , exp -1 x k x k+1 + ex p -1 x k+1 y k , exp -1 x k+1 x k ] ≤ u k+1 , -exp -1 x k+1 x k + v k , -exp -1 x k x k+1 -ρd 2 (x k+1 , x k ) + d k , exp -1 x k x k+1 , ( 14
)
for any v k ∈ A(x k ) and u k+1 ∈ B(x k+1 ). Let (y k , x k , x k+1 ) be the geodesic triangle with θ = ∠(exp -1

x k y k , exp -1 x k x k+1
). Then, by the Comparison Theorem for Triangles [20, Proposition 4.5], we have

d 2 (y k , x k ) + d 2 (x k , x k+1 ) -2 exp -1 x k y k , exp -1 x k x k+1 ≤ d 2 (y k , x k+1 ). ( 15
)
Similarly, to the geodesic triangle

(y k , x k+1 , x k ) with θ = ∠(exp -1 x k+1 y k , exp -1 x k+1 x k ), we have d 2 (y k , x k+1 ) + d 2 (x k+1 , x k ) -2 exp -1 x k+1 y k , exp -1 x k+1 x k ≤ d 2 (y k , x k ). ( 16
)
Adding ( 15) and ( 16), we obtain

d 2 (x k+1 , x k ) ≤ exp -1 x k y k , exp -1 x k x k+1 + exp -1 x k+1 y k , exp -1 x k+1 x k . ( 17
)
Since μ k > 0, for all k ∈ N, and using ( 17) in ( 14), we have

1 μ k d 2 (x k+1 , x k ) ≤ u k+1 , -exp -1 x k+1 x k + v k , -exp -1 x k x k+1 -ρd 2 (x k+1 , x k ) + d k , exp -1 x k x k+1 . ( 18
)
Thus, using the fact that 0 < μ k ≤ b, we obtain

ρ + 1 b d 2 (x k+1 , x k ) ≤ u k+1 , -exp -1 x k+1 x k + v k , -exp -1 x k x k+1 + d k , exp -1 x k x k+1 , ( 19
)
for any v k ∈ A(x k ) and u k+1 ∈ B(x k+1 ).

On the other hand, using the parallel transport properties, we have that exp -1

x k x k+1 = -P x k ,x k+1 exp -1
x k+1 x k , and hence

u k+1 , -exp -1 x k+1 x k + v k , -exp -1 x k x k+1 = P x k ,x k+1 u k+1 -v k , exp -1 x k x k+1
and hence, using this fact in [START_REF] Attouch | A general duality principle for the sum of two operators[END_REF], one has

(ρ + 1 b )d 2 (x k+1 , x k ) ≤ P x k ,x k+1 u k+1 -v k + d k , exp -1 x k x k+1 , ( 20
)
for any v k ∈ A(x k ) and u k+1 ∈ B(x k+1 ). Therefore, the desired result follows applying the Cauchy-Schwarz inequality and triangular inequality in [START_REF] Sakai | Riemannian geometry. Translations of mathematical monographs[END_REF] and using the fact that, in Hadamard manifolds, exp -1

x k x k+1 = d(x k+1 , x k ) and d k = γ k d(x k , x k-1 ). Proposition 6 Suppose that {x k } is bounded, then there exist constants L, M ≥ 0 such that ||P x k ,x k+1 u k+1 -v k || ≤ L and lim sup k→∞ d(x k+1 , x k ) = M, for any v k ∈ A(x k ), u k+1 ∈ B(x k+1
) and k large enough. Moreover, if A or B is ρ-strongly monotone and

ρ 2 + 1 b M > L, ( 21
)
then lim k→∞ d(x k+1 , x k ) = 0. Consequently lim k→∞ d k = 0.
Proof Since {x k } is bounded it follows from the maximal monotonicity of A and B (see Proposition 2) that there exist constants

K 1 , K 2 > 0 such that ||v k || < K 1 and ||u k || < K 2
for k large enough. Thus, from triangular inequality, we have

||P x k ,x k+1 u k+1 -v k || ≤ ||P x k ,x k+1 u k+1 + v k ||.
Using the fact that the parallel transport mapping is an isometry, we have

||P x k ,x k+1 u k+1 -v k || < L = K 1 + K 2 , ( 22
)
for k large enough. On the other hand, since {x k } is bounded, we have that {d(x k+1 , x k )} is also bounded, and hence, there exists M = lim sup k→∞ d(x k+1 , x k ) and the first inequality is proved. Now, if M = 0, then the second assertion directly follows. Otherwise, suppose that M > 0 and, by assumption, we have that A or B is ρ-strongly monotone with ρ > 0 satisfying [START_REF] Udriste | Convex functions and optimization algorithms on Riemannian manifolds[END_REF]. Thus, combining [START_REF] Alimohammady | Inertial proximal algorithm for difference of two maximal monotone operators[END_REF] with ( 22), we have

ρ + 1 b d(x k+1 , x k ) ≤ P x k ,x k+1 u k+1 -v k + ρ 2 d(x k , x k-1 ) ≤ L + ρ 2 d(x k , x k-1 ).
Taking the lim sup k→∞ in the above inequality, this implies that ρ 2 + 1 b M ≤ L which is a contradiction. Therefore, the assertion is proved.

Remark 6 If A(•) = ∂ g(•)
and B(•) = ∂h(•) with g and h convex functions (and hence, maximal monotone vector fields) and f (x) = g(x)h(x), then (21) can be replaced by the lower boundedness of f . However, the key of the proof still follows from an inequality like [START_REF] Alimohammady | Inertial proximal algorithm for difference of two maximal monotone operators[END_REF], more precisely, Proposition 5]. Dealing with general maximal monotone vector fields, it is worth to mention that condition [START_REF] Udriste | Convex functions and optimization algorithms on Riemannian manifolds[END_REF] is not restrictive due to the fact that the parameter ρ can be taken large enough so that (21) holds; see Remark 3.

1 b d 2 (x k+1 , x k ) ≤ f (x 0 ) -f * , where f * = inf x∈M f (x); see [14,
Next, we present partial and full convergence results for IPPM for the more general case of difference of monotone vector fields in Hadamard manifolds. Furthermore, we introduce an assumption on the proximal parameter {μ k } and a Lipschitz continuity of B in order to obtain full convergence of the sequence to a critical point of A -B.

Partial convergence analysis

In this subsection, we suppose that the assumptions in Proposition 6 hold, i.e., we assume that {x k } is bounded, A or B is ρ-strongly monotone and (21) holds. Furthermore, we consider a, b > 0 and {μ k } such that

(23) a ≤ μ k ≤ b. Theorem 1 Every cluster point of {x k } is a critical point of A -B.
Proof Let μ, x and y be cluster points of {μ k }, {x k } and {y k }, respectively. Without loss of generality we can take subsequences {μ k j }, {x k j } and {y k j } converging respectively to μ, x and y (we can extract another subsequence if necessary). From Proposition 6, we have that x k j +1 → x. It follows from ( 7) that

1 μ k j exp -1 x k j y k j -d k j ∈ B(x k j ).
By the monotonicity of B, we have

1 μ k j exp -1 x k j y k j -d k j , exp -1 x k j z ≤ u, -exp -1 z x k j , ∀u ∈ B(z), z ∈ M.
Letting j → +∞ in last inequality and using Proposition 1, we obtain

1 μ exp -1 x y, exp -1 x z ≤ u, -exp -1 z x , ∀u ∈ B(z), z ∈ M,
and hence, by the maximality of B, we have

1 μ exp -1 x y ∈ B(x). ( 24
)
Now, from (8), we have

1 μ k j exp -1 x k j +1 y k j ∈ A(x k j +1 ).
Similarly, by the maximal monotonicity of A, we have

1 μ k j exp -1 x k j +1 y k j , exp -1 x k j +1 z ≤ u, -exp -1 z x k j +1 , ∀u ∈ A(z), ∀z ∈ M,
and letting j → +∞, we show that

1 μ exp -1 x y ∈ A(x). ( 25
)
Therefore, from ( 24) and ( 25), we have that 1 μ exp -1 x y ∈ A(x) ∩ B(x) which means that (6) holds, i.e., x is critical point of A -B.

Sufficient conditions for boundedness

Now, for the sake of completeness, we state the following technical tool which will be used in the sequel.

Lemma 1 Let {α k }, {β k } and { k } three sequences of non-negative numbers satisfying α k+1 ≤ (1 + β k )α k + k . If ∞ k=0 β k < +∞ and ∞ k=0 k < +∞, then {α k } is convergent.
Proof See Polyak [START_REF] Polyak | Introduction to optimization[END_REF].

Recently, some works have proved the convergence of the whole sequence for several methods applied to DC functions supposing the Kurdyka-Łojasiewicz property of the objective function; see [START_REF] Aragon Artacho | The boosted difference of convex functions algorithm for non-smooth functions[END_REF][START_REF] Neto | A generalized proximal linearized algorithm for DC functions with application to the optimal size of the firm problem[END_REF][START_REF] Ferreira | Boosted scaled subgradient method for DC programming[END_REF][START_REF] Thi | Convergence analysis of difference-of-convex algorithm with sub-analytic data[END_REF]. However, even in this case, it is supposed that the sequence generated by the method is bounded. The aim of this subsection is to introduce sufficient conditions in order to guarantee the boundedness of the sequence generated by our method. It is worth to mention that this result is new even for DC problems in the Euclidean setting.

In the remain of this section, we consider Algorithm IPPM with the following condition instead of (23):

0 < μ k ≤ b, k ∈ N. (26)
Furthermore, we consider the following assumptions in the vector field B, the proximal parameters μ k and to sequence of non-negative numbers {γ k d(x k , x k-1 )}. (A1) Given x, y ∈ M, there exists a constant L > 0 such that u -P x,y v ≤ Ld(x, y), for any u ∈ B(x) and v ∈ B(y).

(A2) The sequence of positive numbers {μ k } satisfies

∞ k=0 μ k < ∞. (A3) The sequence {γ k d(x k , x k-1 )} is bounded.
Remark 7 Note that (A1) is Lipschitz type assumption for vector fields. It is worth to mention that assumption (A1) is a natural extension of the Lipschitz continuity of the gradient function from DC function to difference of maximal monotone vector fields. In the literature of algorithm for DC functions some papers have considered the non-smooth case where f (x) = g(x)h(x) with g, h convex function, g possibly non-smooth and h a function C 1,1 , i.e., h is differentiable and its gradient is Lipschitz continuous; see for instance [START_REF] Aragon Artacho | The boosted difference of convex functions algorithm for non-smooth functions[END_REF][START_REF] Neto | A generalized proximal linearized algorithm for DC functions with application to the optimal size of the firm problem[END_REF][START_REF] Ferreira | Boosted scaled subgradient method for DC programming[END_REF][START_REF] Thi | Convergence analysis of difference-of-convex algorithm with sub-analytic data[END_REF]. In the context of vector fields, condition (A1) implies that the vector field B is single-valued. On the other hand, (A2) is a classical assumption on proximal point methods; see for instance [START_REF] Ferreira | Proximal point algorithm on Riemannian manifolds[END_REF]. Dealing with DC functions it has been natural to suppose boundedness assumption of the sequence for different methods; see for instance [START_REF] Souza | A proximal point method for DC functions on Hadamard manifolds[END_REF][START_REF] Maingé | Convergence of new inertial proximal methods for DC programming[END_REF][START_REF] Oliveira | Level an inertial algorithm for DC programming[END_REF][START_REF] Almeida | A modified proximal point method for DC functions on Hadamard manifolds[END_REF][START_REF] Aragon Artacho | The boosted difference of convex functions algorithm for non-smooth functions[END_REF][START_REF] Neto | A generalized proximal linearized algorithm for DC functions with application to the optimal size of the firm problem[END_REF][START_REF] Ferreira | Boosted scaled subgradient method for DC programming[END_REF][START_REF] Thi | Convergence analysis of difference-of-convex algorithm with sub-analytic data[END_REF][START_REF] Bento | Proximal point method for a special class of non-convex functions on Hadamard manifolds[END_REF]. Some works have replaced this assumption by the concept of Kurdyka-Łojasiewicz inequality; see for instance [START_REF] Aragon Artacho | The boosted difference of convex functions algorithm for non-smooth functions[END_REF][START_REF] Neto | A generalized proximal linearized algorithm for DC functions with application to the optimal size of the firm problem[END_REF][START_REF] Ferreira | Boosted scaled subgradient method for DC programming[END_REF][START_REF] Thi | Convergence analysis of difference-of-convex algorithm with sub-analytic data[END_REF]. Here, we consider the boundedness assumption (A3) in order to obtain that {x k } is bounded. It is quite natural because γ k ∈ [0, ρ 2 ), for all k ∈ N, and the proximal point method and its variants usually have the property that {d(x k , x k-1 )} converges to zero.

Theorem 2 Suppose that assumptions (A1), (A2) and (A3) hold. Then, {d(x k , x * )} is convergent for every x * ∈ S. In particular, {x k } is bounded.

Proof By (8), we have that exp -1 x k+1 y k ∈ μ k A(x k+1 ) and hence

y k ∈ exp x k+1 μ k A(x k+1 ).
This means that x k+1 = J A μ k (y k ) and from (7), we have

x k+1 = J A μ k (exp x k μ k (w k + d k )),
where w k = B(x k ). Let x * ∈ S be fixed and take w * ∈ A(x * ) ∩ B(x * ) such that

x * = J A μ k (exp x * μ k w * ).
Since A is maximal monotone from Proposition 4 its resolvent is firmly nonexpansive. Thus, from Proposition 3, we have

d(x k+1 , x * ) = d(J A μ k (exp x k μ k (w k + d k )), J A μ k (exp x * μ k w * )) ≤ d(exp x k μ k (w k + d k ), exp x * μ k w * ) ≤ d(exp x k μ k (w k + d k ), x * ) + d(x * , exp x * μ k w * ) ≤ d(x k , x * ) + d(exp x k μ k (w k + d k ), x k ) + μ k w * = d(x k , x * ) + μ k w k + d k + μ k w * ≤ d(x k , x * ) + μ k ( w k -P x k ,x * w * + P x k ,x * w * + d k ) + μ k w * ≤ d(x k , x * ) + μ k (Ld(x k , x * ) + w * + d k ) + μ k w * = (1 + μ k L)d(x k , x * ) + μ k γ k d(x k , x k-1 ) + 2μ k w * ,
where the triangular inequality was successively applied and in last inequality we used (A1).

Therefore, applying Lemma 1 with

α k = d(x k , x * ), β k = Lμ k and k = μ k γ k d(x k , x k-1
) + 2μ k w * taking into account that assumptions (A2) and (A3) hold, then we conclude that {d(x k , x * )} is convergent and hence {x k } is bounded.

In Proposition 6, we proved that lim k→∞ d(x k+1 , x k ) = 0 under the assumption [START_REF] Udriste | Convex functions and optimization algorithms on Riemannian manifolds[END_REF]. In the absence of this assumption, we prove in the sequel that this assertion still holds under the conditions (A1), (A2) and (A3) for the choice of μ k as in [START_REF] Aragon Artacho | The boosted difference of convex functions algorithm for non-smooth functions[END_REF].

Corollary 2 If assumptions (A1), (A2) and (A3) hold, then

lim k→+∞ d(x k+1 , x k ) = 0.
Proof From (18), we have

1 μ k d 2 (x k+1 , x k ) ≤ u k+1 , -exp -1 x k+1 x k + v k , -exp -1 x k x k+1 + d k , exp -1 x k x k+1 , u k+1 = B(x k+1 ) and v k ∈ A(x k ).
Using the parallel transport, we have that exp -1

x k x k+1 = -P x k ,x k+1 exp -1 x k+1 x k and then, 1 μ k d 2 (x k+1 , x k ) ≤ P x k ,x k+1 u k+1 -v k + d k exp -1 x k x k+1 .
This implies that

0 ≤ d(x k+1 , x k ) ≤ μ k ( u k+1 + v k + γ k d(x k , x k-1 )) ( 27
)
taking into account that exp -1 x k x k+1 = d(x k+1 , x k ) > 0. Since, from Theorem 2, {x k } is bounded, we have from Proposition 2 and (A3) that the sequences {u k }, {v k } and {γ k d(x k , x k-1 )} are bounded. Therefore, letting k → +∞ in [START_REF] Neto | A generalized proximal linearized algorithm for DC functions with application to the optimal size of the firm problem[END_REF] taking into account that from (A2) we have that lim k→+∞ μ k = 0, we obtain that lim k→+∞ d(x k+1 , x k ) = 0 and the assertion is proved.

Sufficient condition for full convergence

In this section, we provide a sufficient condition to obtain the convergence of the whole sequence generated by our method. Under the assumption (A1), (A2) and (A3) we proved that lim k→+∞ d(x k+1 , x k ) = 0 which is a classical behaviour of the proximal point method. It follows from (A2) that lim k→+∞ μ k = 0. Since this auxiliary sequence of parameters is freely chosen satisfying some suitable conditions we will take it in such a way that it does not go to zero faster than d(x k+1 , x k ). This leads us to the following condition:

(A3 * ) Suppose that lim k→∞ d(x k+1 , x k ) μ k = 0.
Remark 8 Note that under assumption (A2), condition (A3 * ) implies (A3). Indeed, (A3 * ) implies that lim k→∞ d(x k+1 , x k ) = 0 due to the fact that lim k→∞ μ k = 0 from (A2).

Since 0 ≤ γ k < ρ 2 , we have that lim k→∞ γ k d(x k , x k-1 ) = 0 and hence, it is bounded. Therefore, replacing assumption (A3) by (A3 * ) the results of Sect. 4.2 remain true and we additionally obtain the convergence of the method as stated in the next result. Recall that we are considering algorithm IPPM with [START_REF] Aragon Artacho | The boosted difference of convex functions algorithm for non-smooth functions[END_REF]. Proof In view of Theorem 2 it is enough to show that every cluster point of {x k } belongs to S. Indeed, if x * is a cluster point of {x k } which belongs to S, then there exists a subsequence {x k j } such that {d(x k j , x * )} converges to zero. Since from Theorem 2, we have that {d(x k , x * )} is convergent, thus it converges to zero, and hence, {x k } converges to x * ∈ S. Now, we shall prove that an arbitrary cluster point x of {x k } belongs to S. Let {x k j } be a subsequence of {x k } such that x k j → x. From Corollary 2, we have that lim j→∞ x k j +1 = x. Moreover, combining [START_REF] Moreau | Proximité et dualité dans un espace Hilbertien[END_REF] with Theorem 2 and Proposition 2, we have that {w k } is bounded. Therefore, without loss of generality, we may assume that there exists a subsequence {w k j } converging to w. Since w k j = B(x k j ) from the monotonicity of B, we have

w k j , exp -1 x k j z ≤ u, -exp -1 z x k j , ∀u = B(z), z ∈ M.
Letting j → ∞, we obtain that w = B(x). On the other hand, from (8), we have

1 μ k j exp -1 x k j +1 y k j ∈ A(x k j +1 ),
where y k j = exp x k j μ k j w k j . We claim that lim j→∞

1 μ k j exp -1
x k j +1 y k j = w, and thus, w ∈ A(x). Indeed,

1 μ k j exp -1 x k j +1 y k j -P x k j +1 ,x w ≤ 1 μ k j exp -1 x k j +1 y k j -P x k j +1 ,x k j (w k j + d k j ) + P x k j +1 ,x k j (w k j + d k j ) -P x k j +1 ,x w . ( 28
)
Let us prove that the right-hand side of the above inequality vanishes as j → ∞. First, note that

1 μ k j exp -1 x k j +1 y k j -P x k j +1 ,x k j (w k j + d k j ) 2 = 1 μ 2 k j exp -1 x k j +1 y k j 2 + 1 μ 2 k j exp -1 x k j y k j 2 -2 1 μ k j exp -1 x k j +1 y k j , P x k j +1 ,x k j (w k j + d k j ) = 1 μ 2 k j d 2 (x k j +1 , y k j ) + d 2 (y k j , x k j ) -2 1 μ k j exp -1 x k j +1 y k j , P x k j +1 ,x k j (w k j + d k j ) = 1 μ 2 k j d 2 (x k j +1 , y k j ) + d 2 (y k j , x k j ) - 2 μ k j P y k j ,x k j +1 (exp -1 x k j +1 y k j ), P y k j ,x k j +1 (P x k j +1 ,x k j (w k j + d k j )) = 1 μ 2 k j d 2 (x k j +1 , y k j ) + d 2 (y k j , x k j ) - 2 μ k j -exp -1 y k j x k j +1 , P y k j ,x k j ( 1 μ k j exp -1 x k j y k j ) = 1 μ 2 k j d 2 (x k j +1 , y k j ) + d 2 (y k j , x k j ) -2 1 μ k j -exp -1 y k j x k j +1 , - 1 μ k j exp -1 y k j x k j = 1 μ 2 k j d 2 (x k j +1 , y k j ) + d 2 (y k j , x k j ) -2 1 μ 2 k j exp -1 y k j x k j +1 , exp -1 y k j x k j . ( 29
)
From the geodetic triangle (x k j +1 , y k j , x k j ), we have

d 2 (x k j +1 , y k j ) + d 2 (y k j , x k j ) -2 exp -1 y k j x k j +1 , exp -1 y k j x k j ≤ d 2 (x k j +1 , x k j ).
From the geodetic triangle (x k j +1 , y k j , x k j ), we have

d 2 (x k j +1 , y k j ) + d 2 (y k j , x k j ) -2 exp -1 y k j x k j +1 , exp -1 y k j x k j ≤ d 2 (x k j +1 , x k j ).
Using last inequality in ( 29), we have

1 μ k j exp -1 x k j +1 y k j -P x k j +1 ,x k j (w k j + d k j ) ≤ 1 μ k j d(x k j +1 , x k j ) which from (A4) implies lim j→∞ 1 μ k j exp -1 x k j +1 y k j -P x k j +1 ,x k j (w k j + d k j ) = 0. ( 30
)
Now, by (A1) and triangular inequality, we have that

P x k j +1 ,x k j (w k j + d k j )P x k j +1 ,x w = P x k j +1 ,x k j (w k j + d k j ) -P x k j +1 ,x k j (P x k j ,x w) = w k j + d k j -P x k j ,x w ≤ Ld(x k j , x) + γ k j d(x k j , x k j -1 ).
From Corollary 2 and x k j → x, letting j → ∞ in the above equality we obtain

lim j→∞ P x k j +1 ,x k j (w k j + d k j ) -P x k j +1 ,x w = 0. ( 31
)
Therefore, using ( 30) and ( 31) in ( 28), we have

lim j→∞ 1 μ k j exp -1 x k j +1 y k j = w.
Using now the monotonicity of A, we have

1 μ k j exp -1 x k j +1 y k j , exp -1 x k j +1 z ≤ u, -exp -1 z x k j +1 , ∀u ∈ A(z), z ∈ M.
Letting j → ∞ and using the maximality of A, we have that w ∈ A(x) and thus w ∈ A(x) ∩ B(x) implying that x ∈ S. This completes the proof.

Convergence analysis: monotone case

As we mention before, IPPM is new even for finding a zero of a maximal monotone vector field. In this section, we consider the variational problem (2) with B(x) = 0 and the algorithm mIPPM. Here, we will not need to suppose assumption (A1).

In order to prove the convergence of mIPPM let us consider A -1 (0) ⊂ M the set of singularities of A and assume that A -1 (0) = ∅, i.e., A -1 (0) = {x ∈ M; 0 ∈ A(x)}.

In this section, we will consider two possible choices for {μ k } (26) and ( 23) as in the non-monotone case, i.e., 0 < μ k ≤ b and 0 < a ≤ μ k ≤ b. In the first case, we assume that assumptions (A2) and (A3 * ) hold. In the second one, we replace these assumptions by

∞ k=0 γ k d(x k , x k-1 ) < ∞. ( 32
)
Remark 9 In the linear setting, condition ( 32) is considered for instance by Alimohammady and Ramazannejad [START_REF] Alimohammady | Inertial proximal algorithm for difference of two maximal monotone operators[END_REF]. In Alvarez and Attouch [START_REF] Alvarez | An inertial proximal method for maximal monotone operators via discretization of a nonlinear oscillator with damping[END_REF], a similar condition is proposed in Hilbert spaces. More precisely, it is assumed that ∞ k=0 γ k x kx k-1 2 < ∞; see [START_REF] Alvarez | An inertial proximal method for maximal monotone operators via discretization of a nonlinear oscillator with damping[END_REF]Theorem 2.1]. Furthermore, such a condition holds for some special cases which can be verified a priori; see [START_REF] Alvarez | An inertial proximal method for maximal monotone operators via discretization of a nonlinear oscillator with damping[END_REF]Proposition 2.1].

Next results ensure that the sequence {x k } generated by mIPPM converges to a singularity of A. We will consider both possible choices of μ k as mentioned before.

Theorem 4 The sequence {d(x k , x * )} is convergent for every x * ∈ A -1 (0) and, in particular, {x k } is bounded.

Proof By [START_REF] Alimohammady | Notes on the difference of two monotone operators[END_REF], we have that exp -1 x k+1 y k ∈ μ k A(x k+1 ) and hence

y k ∈ exp x k+1 μ k A(x k+1 ).
This means that x k+1 = J A μ k (y k ) and from (9), we have

x k+1 = J A μ k (exp x k μ k (d k )).
Let x * ∈ A -1 (0) be fixed, then 0 ∈ A(x * ) such that

x * = J A μ k (exp x * μ k 0) = J A μ k (x * ).
Since A is maximal monotone from Proposition 4 its resolvent is firmly nonexpansive. Thus, from Proposition 3, we have

d(x k+1 , x * ) = d(J A μ k (exp x k μ k d k ), J A μ k (x * )) ≤ d(exp x k μ k d k , x * ) ≤ d(x k , x * ) + d(exp x k μ k d k , x k ) = d(x k , x * ) + μ k d k = d(x k , x * ) + μ k γ k d(x k , x k-1 ),
where the triangular inequality was applied. Therefore, the desired result follows from Lemma 1 with

α k = d(x k , x * ), β k = 0 and k = μ k γ k d(x k , x k-1 ) taking into account that ∞ k=0 k < ∞ for both possible choices of μ k . Indeed, if 0 < μ k ≤ b, then assump- tions (A2) and (A3 * ) hold, and hence, {γ k d(x k , x k-1 )} is bounded and ∞ k=0 μ k < ∞ which implies that ∞ k=0 k < ∞. If 0 < a ≤ μ k ≤ b, then {μ k } is bounded and from (32) we obtain that ∞ k=0 k < ∞.
This completes the proof.

As mentioned in Remark 8, under assumption (A2) we have that condition (A3 * ) implies lim k→∞ d(x k+1 , x k ) = 0. This is the case if {μ k } is given by ( 26), namely, 0 <μ k ≤ b. On the other hand, if {μ k } is givenby [START_REF] Almeida | A modified proximal point method for DC functions on Hadamard manifolds[END_REF], i.e., 0 < a ≤ μ k ≤ b, we prove next that this condition still holds under assumption [START_REF] Nesterov | On the Riemannian geometry defined by self-concordant barriers and interiorpoint methods[END_REF].

Corollary 3 Let {x k } be a sequence generated by mIPPM with {μ k } given by [START_REF] Almeida | A modified proximal point method for DC functions on Hadamard manifolds[END_REF]. Then, lim k→+∞ d(x k+1 , x k ) = 0.

Proof From (10), we have that exp -1 x k+1 y k ∈ μ k A(x k+1 ). Take x * ∈ A -1 (0). By monotonicity of A, we have

1 μ k exp -1 x k+1 y k , exp -1 x k+1 x * ≤ 0, -exp -1 x * x k+1 = 0. ( 33
)
From the geodetic triangle (y k , x k+1 , x * ) with θ = ∠(exp -1 x k+1 y k , exp -1 x k+1 x * ), we obtain

d 2 (y k , x k+1 ) + d 2 (x k+1 , x * ) -2 exp -1 x k+1 y k , exp -1 x k+1 x * ≤ d 2 (y k , x * ). ( 34
)
Combining ( 33) with (34), we have

d 2 (y k , x k+1 ) + d 2 (x k+1 , x * ) ≤ d 2 (y k , x * )
which applying the triangular inequality leads to

d 2 (y k , x k+1 ) ≤ d 2 (y k , x k ) + 2d(y k , x k )d(x k , x * ) + [d 2 (x k , x * ) -d 2 (x k+1 , x * )]. ( 35
)
Since

y k = exp x k μ k d k , one has d(y k , x k ) = μ k γ k d(x k , x k-1
). Thus, we obtain that [START_REF] Nesterov | On the Riemannian geometry defined by self-concordant barriers and interiorpoint methods[END_REF]. Therefore, letting k → ∞ in [START_REF] Bhatia | Positive definite matrices[END_REF] and using the fact that {d(x k , x * )} is convergent, we have that lim k→∞ d(y k , x k+1 ) = 0. Thus, from triangular inequality, we have

lim k→∞ d(y k , x k ) = 0 taking into account that {μ k } is bounded and lim k→∞ γ k d(x k , x k-1 ) = 0 from
0 ≤ lim k→∞ d(x k+1 , x k ) ≤ lim k→∞ [d(x k+1 , y k ) + d(y k , x k )] = 0,
and hence, lim

k→∞ d(x k+1 , x k ) = 0.
Theorem 5 The sequence {x k } converges to a singularity of A.

Proof In view of Theorem 4 it is enough to show that every cluster point of {x k } belongs to S. Indeed, if x * is a cluster point of {x k } which belongs to A -1 (0), then there exists a subsequence {x k j } such that {d(x k j , x * )} converges to zero. Since from Theorem 4, we have that {d(x k , x * )} is convergent, thus it converges to zero, and hence, {x k } converges to x * ∈ A -1 (0). Now, we shall prove that an arbitrary cluster point x of {x k } belongs to A -1 (0). Let {x k j } be a subsequence of {x k } such that x k j → x. From assumption (A3 * ) and Corollary 3, for both possible choices of μ k , we have that lim k→∞ d(x k+1 , x k ) = 0 and hence lim j→∞ x k j +1 = x. From the definition of mIPPM, we have that

d k j = γ k j exp -1 x k j x k j -1 and y k j = exp x k j μ k j d k j . Thus, 1 μ k j exp -1 x k j +1 y k j = 1 μ k j d(x k j +1 , y k j ) ≤ 1 μ k j d(x k j +1 , x k j ) + 1 μ k j d(x k j , y k j ) ≤ 1 μ k j d(x k j +1 , x k j ) + ||d k j ||,
for both choices of μ k . The right-hand side of the above inequality goes to zero as j → ∞ from (A3 * ), for 0 < μ k ≤ b, and from 1

μ k j d(x k j +1 , x k j ) ≤ 1 a d(x k j +1 , x k j ), for 0 < a ≤ μ k ≤ b,
taking into account that lim k→∞ d(x k+1 , x k ) = 0 in both cases. Now using the monotonicity of A for 1

μ k j exp -1 x k j +1 y k j ∈ A(x k j +1 ), we have 1 μ k j exp -1 x k j +1 y k j , exp -1 x k j +1 z ≤ u, -exp -1 z x k j +1 , ∀u ∈ A(z), z ∈ M.
Letting j → +∞ in last inequality and using the fact that 1

μ k j exp -1
x k j +1 y k j → 0 ∈ T x M together with Proposition 1 and the maximality of A, we obtain that 0 ∈ A(x) and the proof is completed.

Numerical experiments

In this section, we solve unconstrained minimization DC problems in a genuine Hadamard manifold considered by Almeida et al. [START_REF] Almeida | A modified proximal point method for DC functions on Hadamard manifolds[END_REF]. We run IPPM and mIPPM comparing its performance with proximal point methods [START_REF] Ferreira | Proximal point algorithm on Riemannian manifolds[END_REF] (PPM in the convex case) and [START_REF] Souza | A proximal point method for DC functions on Hadamard manifolds[END_REF] (DCPPM in the non-convex DC case). We also present some experiments for monotone and non-monotone operators including an operator which its components are not the subdifferential (or gradient) of a convex function and an example of a point-to-set operator in the Euclidean setting.

Let P n be the set of the symmetric matrices, P n + be the cone of the symmetric positive semi-definite matrices and P n ++ be the cone of the symmetric positive definite matrices both n × n. Let M = (P n ++ , •, • ) be the Riemannian manifold endowed with the Riemannian metric given by

U , V = tr(V X -1 U X -1 ), X ∈ M, U , V ∈ T X M, ( 36
)
where tr(X ) denotes the trace of X ∈ P n and T X M ≈ P n , with the corresponding norm denoted by . ; see Rothaus [START_REF] Rothaus | Domains of positivity[END_REF]. In this case, for any X , Y ∈ M the unique geodesic joining those two points is given by:

γ (t) = X 1/2 X -1/2 Y X -1/2 t X 1/2 , γ (0) = X 1/2 ln(X -1/2 Y X -1/2 )X 1/2 , t ∈ [0, 1],
see, for instance, Nesterov and Todd [START_REF] Nesterov | On the Riemannian geometry defined by self-concordant barriers and interiorpoint methods[END_REF] and Bento et al. [START_REF] Bento | Proximal point method for a special class of non-convex functions on Hadamard manifolds[END_REF]. More precisely, M is a Hadamard manifold; see, for instance, [33, Theorem 1.2., page 325]. One can compute the curvature of M and verify that it has non-constant curvature; see Lenglet et al. [34,page 428]. We consider the exponential map and its inverse as follows

exp X V = X 1/2 e X -1/2 VX -1/2 X 1/2 , exp X -1 V = X 1/2 ln(X -1/2 VX -1/2 )X 1/2 ;
see Bhatia [START_REF] Bhatia | Positive definite matrices[END_REF]Chapter 6].

The numerical experiments are coded in MATLAB R2020b on a machine with a Intel(R) Core(TM) i7, 2.3 GHz CPU and 8 GB memory. We consider two kind of problems. In Problems 1 and 2, we consider minimization problems in the Riemannian manifolds M described above with n × n matrices for n = 5, n = 25, n = 50 and n = 100. We perform the methods performed using "MANOPT" which is a toolbox for optimization on manifolds; see Boumal et al. [START_REF] Boumal | Manopt, a Matlab toolbox for optimization on manifolds[END_REF]. We generate random matrices using "manifold.rand" to start the method. The subproblems are solved using the routine "steepestdescent" with the inner loop stopping if the norm of the gradient drops below = 10 -5 as well as the stopping criterion used to the outer loop was d(x k+1 , x k )< . In Problem 3 and 4, we consider examples in R 2 with operators which are not the subdifferential of a convex function or multi-valued operators. In this context, we use the same stopping rule where the Riemannian distance is the Euclidean norm.

For each problem below, we perform the algorithms 100 times using different random starting points. In all problems, we compare the performance of the methods using random initial starting points but the same at each running for all the methods. We compare the performance of the methods for different choices of the parameters μ k and γ k .

Problem 1

In this problem, we consider the convex minimization problem min

X ∈M f (X ),
where f (X ) = 1 2 ln(det(X )) 2 . Clearly, this problem can be seen as 0 ∈ A(X ) -B(X ) for A(X ) = grad f (X ) and B(X ) = 0, where A is a maximal monotone vector field. The solution of the equation 0 = grad f (X * ) is a matrix such that det(X * ) = 1, and hence, f (X * ) = 0. Note that the Euclidean gradient and hessian of f , denoted by f (X ) and f (X ), are given by f (X ) = ln(det(X ))X -1 and f (X )V = tr (X -1 V )X -1ln(det(X ))X -1 V X -1 . Thus, one can verify that f is not convex in the Euclidean sense while it is in the Riemannian setting using the metric given in [START_REF] Boumal | Manopt, a Matlab toolbox for optimization on manifolds[END_REF] and having in mind that the Riemannian gradient and hessian of f , denoted by grad f (X ) and hess f (X ), are given by grad f (X ) = ln(det(X ))X and hess f (X )V = tr (X -1 V )X .

Problem 2

In this problem, we consider the non-convex DC problem min

X ∈M f (X ) = g(X ) -h(X ),
where f (X ) = ln(det(X )) 4 -2 ln(det(X )) 2 + 1 with g(X ) = ln(det(X )) 4 + 1 and h(X ) = 2 ln(det(X )) 2 . Clearly, this problem can be seen as 0 ∈ A(X ) -B(X ) for A(X ) = grad g(X ) and B(X ) = grad h(X ), where A, B are maximal monotone vector fields. The solution of the equation 0 = grad f (X * ) is a matrix such that det(X * ) = 1, det(X * ) = e or det(X * ) = e -1 with min x∈M f (X ) = 0. It is worth to mention that g and h are not convex functions in the Euclidean sense and become convex in the Riemannian sense using the metric [START_REF] Boumal | Manopt, a Matlab toolbox for optimization on manifolds[END_REF]. Despite the convexity of g and h in M, we have that f is a non-convex function in both Euclidean and Riemannian setting.

Next, we consider some examples where the vector fields are defined in the Euclidean space. In Problem 3, we have an operator which its components are not the subdifferential (or gradient) of a convex function. In Problem 4, we have an example of a point-to-set operator.

Problem 3 Let A, B : R 2 → R 2 given by A(x) = x 1 2 -x 2 , x 1 + x 2 2 and B(x) = (x 2 , -x 1 ) with x = (x 1 , x 2 ).
One has that A is strongly maximal monotone and B is maximal monotone but neither of them is the subdifferential of a convex function because they are not cyclically monotone operators; see [START_REF] Rockafellar | Characterization of the subdifferentials of convex functions[END_REF]Theorem 3]. Let us consider the problem of finding

0 ∈ C(x) = A(x) -B(x)
which is equivalent to show that A(x) ∩ B(x) is non-empty, for some x ∈ R 2 . One can easily show that taking x * = (0, 0), we have that 0 ∈ A(x * ) ∩ B(x * ). Thus, x * = (0, 0) is the solution of the above (non-monotone) variational inclusion.

Problem 4

In this problem, we consider the non-monotone variational problem

0 ∈ A(x) -B(x), where A(x) = ∂ g(x) and B(x) = ∂h(x) with x = (x 1 , x 2 ) and g(x) = ||x|| 2 + ||x||, h(x) = 0.5||x|| 2 + max{-x 1 , 0} + max{-x 2 , 0}.
Clearly, A, B are maximal monotone (point-to-set) operators and this problem can be seen as the non-differentiable and non-convex DC problem min

x∈R 2 f (x) = g(x) -h(x) = 0.5||x|| 2 + ||x|| -max{-x 1 , 0} -max{-x 2 , 0}.
The results are presented in Tables 1,2 4, this last column is considered as max ||C(x k )|| which means the maximal value for ||C(x k )|| because in this example the methods converge to different limit points depending of the initial point, so in this case does not make sense to consider the median.

We also present the behavior of the methods in terms of In Tables 1, the results show that for Problem 1, the method mIPMM outperforms the nonboosted version of the method (the case where γ k = 0) in both number of iterates and CPU time for different choices of μ k and γ k in all the dimensions considered. In this example, our method with γ k = 0 coincides with the classical proximal point method (PPM) considered in [START_REF] Ferreira | Proximal point algorithm on Riemannian manifolds[END_REF].

||d k ||, d(x k+1 ,x k ) μ k and | f (X k )| (or ||C(x k )||). It is shown in Figures 1, 2
In Tables 2, the results show that for Problem 2, the method IPMM outperforms the nonboosted version of the method (the case where γ k = 0) in CPU time for different choices of μ k and γ k in all the dimensions considered but underperform it in number of iterates. In this example, our method with γ k = 0 coincides with the classical proximal point method for DC functions (PPMDC) considered in [START_REF] Souza | A proximal point method for DC functions on Hadamard manifolds[END_REF].

In Tables 3, the results show that for Problem 3, the method IPMM outperforms the nonboosted version of the method (the case where γ k = 0) in both number of iterates and CPU time for different choices of μ k and γ k in all the dimensions considered.

In Tables 4, the results show that for Problem 4, the method IPMM outperforms the nonboosted version of the method (the case where γ k = 0) in both number of iterates and CPU time for different choices of μ k and γ k in all the dimensions considered. In this example, our method with γ k = 0 coincides with the classical proximal point method for DC functions (PPMDC) considered in [START_REF] Souza | A proximal point method for DC functions on Hadamard manifolds[END_REF].

Conclusions

We have considered an inertial version of the proximal point method for finding a zero of (non-monotone) difference of two maximal monotone vector fields in Hadamard manifolds. We have proved that every cluster point of the sequence, if any, is a solution of the problem. Furthermore, we have presented some sufficient conditions for boundedness and full convergence of the proximal point method for difference of monotone vector fields which are new Fig. 6 Running PPMDC [START_REF] Souza | A proximal point method for DC functions on Hadamard manifolds[END_REF] and IPPM for Problem 4

  γ (a) to γ (b) is denoted by P γ,γ (b),γ (a) : T γ (a) M → T γ (b) M. For any a, b, the parallel transport P γ,γ (b),γ (a) is an isometry from T γ (a) M to T γ (b) M. Note that, for any a, b, a 1 , b 1 , we have

Theorem 3

 3 If assumptions (A1), (A2) and (A3 * ) hold, then {x k } converges to a critical point of A -B.

  , 3 and 4 where the column n denotes the size of the n × n matrices, the columns μ k and γ k present the values of these parameters used in the algorithms, the columns min. iter. (k) (resp. min. time), max. iter. (k) (resp. max. time) and med. iter. (k) (resp. med. time) stand to the minimal, maximal and median of iterates (resp. CPU time in seconds) of the methods in 100 runs. The columns ||grad f (X k )|| and | f (X k )| (or ||C(x k )||) show the median of these values in 100 runs. In Table

  , 3, 4, 5 and 6.

Fig. 1 Fig. 2 Fig. 3 Fig. 4

 1234 Fig. 1 Computing ||d k || and d(X k+1 , X k )/μ k with PPM [6] and mIPPM for Problem 1

Fig. 5

 5 Fig. 5 Running IPPM for Problem 3

  and p ∈ D(A). Then A is said to be (i) Upper semicontinuous at p if, for any open set V satisfying A( p) ⊂ V ⊂ T p M, there exists an open neighbourhood U ( p) of p such that P p,x A(x) ⊂ V for any x ∈ U ( p); (ii) Upper Kuratowski semicontinuous at p if, for any sequences {x k } ⊆ D(A) and {u k } ⊂ T M with each u k ∈ A(x k ), the relations lim k→∞ x k = p and lim k→∞ u k = u imply that u ∈ A( p).

	(iii) Locally bounded at p if there exists an open neighbourhood U ( p) of p such that the set
	∪ x∈U ( p) A(x) is bounded;
	(iv) Upper semicontinuous (resp. upper Kuratowski semicontinuous, locally bounded) on M
	if it is upper semicontinuous (resp. upper Kuratowski semicontinuous, locally bounded)
	at each point p ∈ D(A).
	Proposition 2 Suppose that A ∈ X (M) is maximal monotone and D(A) = M. Then A is
	locally bounded on M.
	Proof See [5, Lemma 3.6].
	Remark 1 Suppose that A ∈ X (M) is maximal monotone and D(A) = M. If {x k } is bounded
	and v k ∈ A(x k ) for all k ∈ N, then {v k } is bounded. Indeed by boundedness of {x k } there
	exists C > 0 such that d(x 0 , x k ) ≤ C, for all k ∈ N. Define the open neighbourhood of x 0
	given by V = {x ∈ M; d(x 0 , x) < C + 1}. It follows from Proposition 2 and Definition 2
	(iii) that ∪ x∈V A(x) is bounded. Since x k ∈ V and v k ∈ A(x k ), we conclude that {v k } is
	bounded and the claim is proved.

Table 1

 1 Running 100 times PPM[START_REF] Ferreira | Proximal point algorithm on Riemannian manifolds[END_REF] and m I P P M

		(X k
		||grad f
		Med. time
		Max. time
		Min. time
		Med. iter.(k)
	for Problem 1	Max. iter.(k)
		Min. iter.(k)
		k
		γ
		μ k
		n

Table 4

 4 Running 100 times PPMDC functions[START_REF] Souza | A proximal point method for DC functions on Hadamard manifolds[END_REF] and I P P M

	-09
	1.1912e
	-05
	9.7601e

Acknowledgements

We would like to thank the referees for their constructive remarks which allow us to improve our work. J.C.O. Souza was supported in part by CNPq Grants 424169/2018-5 and 313901/2020-1. The project leading to this publication has received funding from the French government under the "France 2030" investment plan managed by the French National Research Agency (reference: ANR-17-EURE-0020) and from Excellence Initiative of Aix-Marseille University -A*MIDEX.