# Inertial Bregman Proximal Gradient under Partial Smoothness

Jean-Jacques Godeme*

Contributing authors: e-mail: jeanjacquesgodeme@gmail.com;

## Abstract

This work considers an Inertial version of Bregman Proximal Gradient algorithm (IBPG) for minimizing the sum of two single-valued functions in finite dimension. We suppose that one of the functions is proper, closed, and convex but non-necessarily smooth whilst the second is a smooth enough function but not necessarily convex. For the latter, we ask the smooth adaptable property (smad) with respect to some kernel or entropy which allows to remove the very popular global Lipschitz continuity requirement on the gradient of the smooth part. We consider the IBPG under the framework of the triangle scaling property (TSP) which is a geometrical property for which one can provably ensure acceleration for a certain subset of kernel/entropy functions in the convex setting. Based on this property, we provide global convergence guarantees when the entropy or kernel is strongly convex under the framework of the Kurdyka-Łojasiewicz (KL) property. Turning to the local convergence properties, we show that when the nonsmooth part is partly smooth relative to a smooth submanifold, IBPG has a finite activity identification property before entering a local linear convergence regime for which we establish a sharp estimate of the convergence rate. We report numerical simulations to illustrate our theoretical results on low-complexity regularized phase retrieval.

**Keywords:** Inertial scheme, Bregman proximal gradient, Partial smoothness, Trap avoidance.

# 1 Introduction

## 1.1 Problem statement

In this work, we study the following, non-necessarily smooth nor convex optimization problem

$$ \inf_{x \in \mathbb{R}^n} \left\{ \Phi(x) \stackrel{\text{def}}{=} F(x) + G(x) \right\}, \quad (\mathcal{P}) $$

under the following assumptions
**Assumption A.**

(A.1) $F: \mathbb{R}^n \to \mathbb{R}$ is $C^1$-smooth,

(A.2) $G: \mathbb{R}^n \to \mathbb{R}$ is a proper, lower semi-continuous and convex,

(A.3) $\Phi$ is bounded from below, i.e. $\inf \Phi > -\infty$.

In many applications like automatic differentiation, inverse problems in particular phase retrieval machine learning for Artificial Intelligence and physics applications, $G$ plays the central role of a penalty or regularization term. It is intended to encode structural properties or prior information about the set of solutions. On the other hand, $F$ is a data fidelity term between the true vector and the estimated vector. Throughout this paper, $F$ can be nonconvex which allows us to cover a various problems of these fields. Before, we describe the method that we propose to solve this problem, we recall the definition of Bregman divergence, an important notion that will be at the heart in both our algorithm and his theoretical analysis.

**Definition 1.1. (Bregman divergence)** The Bregman divergence with respect to a convex function $\phi$ is given by the following expression

$$D_{\phi}^{v}(x, y) \stackrel{\text{def}}{=} \begin{cases} \phi(x) - \phi(y) - \langle v; x - y \rangle, & \text{if } (x, y) \in (\text{dom}\phi \times \text{intdom}\phi), v \in \partial\phi(y), \\ + \infty & \text{otherwise.} \end{cases} \quad (1.1)$$

Here, $\phi$ is the Bregman entropy or Kernel which generates the Bregman divergence $D_{\phi}^{v}$.

To solve the optimization problem (P), we will link $\Phi$ to a Bregman divergence $\psi$ with respect to which, the smooth part $F$ is relatively smooth: a notion that will be properly explained in Section A. In this work, we propose Algorithm 1 which we name Inertial Bregman Proximal Gradient (IBPG).

<caption>**Algorithm 1:** Inertial Bregman Proximal Gradient</caption>

**Parameters:** $\kappa \in ]1, 2]$;
**Initialization:** $z_{-1} = x_{-1}, z_0 = x_0 \in \mathbb{R}^n, a_{-1} = a_0 = 1$, and $0 < \underline{a} < \bar{a} \le 1$;
**for** $k = 0, 1, \dots$ **do**
    $y_k = z_k + a_k(x_k - z_k)$;
    $x_{k+1} = (\nabla\psi + \gamma_k \partial G)^{-1} (\nabla\psi(y_k) - \gamma_k \nabla F(y_k))$, $\gamma_k = \frac{a_k^{\kappa-1}}{L}$;
    $z_{k+1} = x_k + a_k(x_{k+1} - x_k)$;
    Choose $a_{k+1} \in [\underline{a}, \bar{a}]$ s.t. $(a_{k+1}^{1-\kappa} + 1)^{1/\kappa} (1 - a_{k+1}) < a_k^{1/\kappa-1} / (1 - a_k)$.
    (IBPG)
It can be easily shown that $y_k$ can be written as a recursion of $(x_k, x_{k-1})$ with inertial parameters that depend on $(a_k, a_{k-1})$. When $a_k \equiv 1$, we recover the Bregman Proximal Gradient (BPG without inertia) whose global convergence was already established in [20].

In [35], a slightly different algorithm called Accelerated Bregman Proximal Gradient, was considered in the convex case. This scheme was inspired by the Improved Interior Gradient Algorithm for conic optimization as developed in [7, Section 5]. It was shown in [35] that Accelerated Bregman Proximal Gradient, indeed, provides acceleration when the entropy $\psi$ satisfies the triangle scaling property (TSP) (see Definition A.9) with a convergence rate of $O(k^{-\kappa})$, where $\kappa \in ]1, 2]$ is the triangle scaling exponent (TSE). For $\kappa = 2$, one recovers the famous Nesterov-like accelerated rate.

In the nonconvex case, another inertial scheme was analyzed in [56] when the Bregman entropy is also strongly convex and $G$ is weakly convex. Global convergence of the iterates under the KL property was proved there using line search on both the extrapolation (inertial) parameter and the descent step-size. Adding more structure of the choice of the Bregman entropy, (i.e. TSP) will allow to have a sharper analysis and results.

## 1.2 Prior work

### Inertial Bregman proximal gradient algorithms

Inertial methods emerged from the quest of accelerating the convergence of first-order optimization methods such as gradient descent. This starts with the Heavy-ball method with friction [60] for gradient descent ($G \equiv 0$). This approach can be interpreted as a discretization of a nonlinear second-order dynamical system, specifically an oscillator with viscous damping. This idea permeates now all optimization techniques and has been applied for instance to the proximal point method [1, 2] and to the inertial Forward-Backwards type methods [55, 5, 53]. In terms of acceleration for convex programming, the accelerated FISTA method [12, 57] achieves a convergence rate of $O(k^{-2})$ for the sequence of objective functions which has been improved to $o(k^{-2})$ in [4], with convergent iterates in [26].

For Bregman-based methods, first-order methods achieve the convergence $O(k^{-1})$ for the sequence of objective [14, 9, 20, 54, 66] as in the Euclidean case. A natural question is whether the Bregman proximal gradient algorithm can be accelerated in the relative smooth framework. One can notice many attempts to answer this question within [9, 54], and the survey paper [66, Section 6]. Positive answers have already been provided under strict additional regularity assumptions. When the Bregman entropy $\psi$ is a strongly convex Legendre kernel and the smooth part of the objective has a Lipschitz continuous gradient, the Improved Interior Gradient Algorithm in [7] admits an accelerated $O(1/k)$ convergence rate on the objective, by using the same inertial technique as Nesterov-type methods. For a subset of relatively smooth functions, [35] shows that the convergence rate of the objective can be improved from $O(1/k)$ to $O(1/k^\kappa)$ where $\kappa \in [1, 2]$ is determined by some crucial triangle scaling property of the Bregman distance, whose genericity is unclear. The general case was still opened until the work of [29] which showed that the $O(1/k)$ rate is optimal for first-order algorithms over the subset of relatively smooth functions, and this cannot be improved in general. As far as the nonconvex case is concerned, inertial versions of the proximal gradient method were analyzed in [56, 70, 39, 72] when the entropy is strongly convex. All these works use either backtracking or line search on both the extrapolation (inertial) parameter and the descent step-size.
Global convergence of the iterates under KL were proved in [56, 72] while [70] showed linear convergence under certain error bound condition.

## Activity identification

Finite activity identification of underlying manifolds is an important phenomenon occurring when different types of algorithms are used to solve structured nonsmooth minimization problems. Such analysis can be traced back to the work of [24] and [32], where a nondegeneracy condition is used to ensure that the optimal active constraints are identified after finitely many iterations. They showed how polyhedral faces of convex sets could be identified finitely. Later [71] extended these results by introducing the concept of smooth identifiable surface and providing an algorithm that identifies the active constraints in convex problems in a finite number of iterations. This work was then generalized in [44] to the notion of partial smoothness for general non necessarily convex problems. Among algorithms that identify active manifolds of partly smooth functions, one can cite the (sub)gradient projection method, Newton-like methods, and the proximal point algorithm as shown in [38, 37, 36] have shown for the (sub)gradient projection method, Newton-like methods, and the proximal point algorithm. A comprehensive study of finite activity identification as well as sharp local linear convergence has been established by Liang et al. in a series of papers for a variety of operator splitting algorithms: forward-backward-type algorithms including accelerated ones [46, 48, 47], for Douglas-Rashford/ADMM [50] and for the Primal-Dual splitting [49].

## Strict saddle avoidance

A driving theme in nonconvex optimization, supported by empirical evidence, is that simple algorithms often work well in highly nonconvex and even nonsmooth settings by avoiding “bad” critical points. A growing body of literature provides one compelling explanation for this escape property or trap avoidance phenomenon. Namely, typical smooth objective functions provably satisfy the strict saddle property, meaning each critical point is either a local minimizer or has a direction of strictly negative curvature. For such functions, either randomly initialized gradient-type methods [34, 58], or stochastically perturbed gradient methods [59, 23] provably escape all strict saddle points, generically on initialization or on the noise. At the heart of the analysis, in all these work is the use of the Center Stable Manifold Theorem [63, Theorem III.7] which finds its roots in the work of Poincaré. In [34], it was proved that the heavy ball method with friction, applied to a $C^2$-Morse function, provably converges to a local minimizer generically with any initialization. Morse functions are known to be generic in the Baire sense in the space of $C^2$-functions (see [6]). In [43], it has been shown that when the objective function is sufficiently smooth with a Lipschitz continuous gradient, a broader class of first-order methods avoid strict saddle points. In the nonsmooth setting, Euclidean proximal methods as explored in [28], are effective in avoiding a nonsmooth version of strict saddle points: “active” strict saddle points, which are strict saddle points with respect to an underlying activity manifold. In parallel to [28], we would like to mention the recent work of [13] which shows that stochastic subgradient descent also escapes “active” strict saddle points. The last two papers rely on important geometric conditions that turn out to be generic in the space of tame weakly convex functions.
## 1.3 Contributions

In this work, we study the global and local convergence properties of Algorithm 1 for a subclass of Bregman kernels that satisfies the TSP property. The main contributions of this work are the following

### Global convergence of the scheme

Under the KL-property, we show that bounded iterates generated by Algorithm 1 converge to a critical point of the objective function when the inertial parameters satisfy the condition that $a_k \in [\bar{a}, \underline{a}] \subset ]0, 1[$. We also show that starting near an optimal solution, the generated sequence converges to it.

### Finite activity identification

We establish that sequence generated by IBPG enjoy a finite activity identification property. More precisely, we show that when the nonsmooth part $G$ is partly smooth with respect to an underlying manifold $M_{x_*}$ near some critical point $x_*$, the iterates will identify the manifold under an appropriate nondegeneracy condition. This identification property implies the existence of some large integer $K$ such that the iterates generated after this number lie in the manifold $M_{x_*}$. The nondegeneracy condition cannot be relaxed in general as shown in [38, Example 4.1, 4.2, 4.3] respectively for the projected gradient descent, Newton method, and the proximal point algorithm.

### Local linear convergence analysis

After activity identification, we show that locally along the active manifold $M_{x_*}$, we can linearize the iterates generated by IBPG when $F$ is locally $C^2$. In this case, we provide a spectral analysis of the linearized system, and under appropriate restricted injectivity, we exhibit a linear convergence regime for proper choice of the inertial parameter. This choice depends in particular on the TSE parameter $\kappa$, which generalizes the Euclidean case for which this parameter is just 2. Our work hence extends that of [47] to non-euclidean, Bregman-based geometry.

### Escape property in the smooth case

Equipped with the center stable manifold theorem, we study the trap avoidance of the inertial Bregman gradient method where $G = 0$. We show that the scheme generically avoids strict saddle points, which are critical points where the function has at least one direction of negative curvature.

Several numerical results are reported that confirm all our theoretical findings.

## 1.4 Paper organization

The rest of the paper is organized as follows. In section 2, we collect the notations and general definitions needed throughout this work. Section 3 contains the global convergence analysis of the Inertial Bregman Proximal Gradient scheme with the Kurdyka-Łojasiewicz property. Section 4 is a self-contain local analysis of the scheme that includes a finite activity identification of underline manifolds with the notion of partial smoothness, a local linearization of the iterates, and finally a local linear convergence. Section 5 contains a study of trap avoidance
in the smooth case. In Section 6, we highlight our results with a nonlinear inverse problems that arise in many areas of imaging science with applications that include diffraction imaging, astronomical imaging, and microscopy namely phase retrieval with prior information on the signal that we want to recover such as sparsity, low-rank, etc. All the technical proofs, the Bregman and KL toolbox are postponed in the Appendices.

# 2 Preliminaries

## 2.1 Notations

In this work, for any nonempty convex set $\Omega \subseteq \mathbb{R}^n$, the following subsets $\text{int}(\Omega)$, $\text{ri}(\Omega)$, $\text{bd}(\Omega)$, $\text{rbd}(\Omega)$ and $\text{aff}(\Omega)$ denote respectively its interior, relative interior, boundary, relative boundary, and its affine hull. We also denote by $P_\Omega$ the orthogonal projector onto $\Omega$. The subset parallel to the set $\Omega$ passing through the origin 0 is denoted $\text{par}(\Omega) = \text{aff}(\Omega) - \Omega$. We denote by $\langle \cdot, \cdot \rangle$ the usual scalar product and $\|\cdot\|$ his corresponding norm. $B(x, r)$ is the corresponding ball of radius $r$ centered at $x$ and $\mathbb{S}^{n-1}$ is the corresponding unit sphere. $\Gamma_0(\mathbb{R}^n)$ is the usual set of proper, closed, and convex functions. $\mathbb{N}$ is the set of non-negative integers. For $m \in \mathbb{N}^*$, $[m]$ is the discrete set $\{1, \cdots, m\}$. For any matrix $A$, $A^\top$ denote the transpose of $A$ and $A^*$ the adjoint (transpose conjugate) when $A$ is a complex matrix. We used the script $\ll$ for the partial order of Loewner between positive semidefinite matrices. $\text{dom}(f)$ is the domain of the function $f$. $f^*$ denotes the Legendre-Fenchel conjugate of $f$.

## 2.2 General Definitions

Throughout this section, let $g: \mathbb{R}^n \to \bar{\mathbb{R}} = \mathbb{R} \cup \{\infty\}$ be a function such that $g$ is proper and lower-semicontinuous.

**Definition 2.1. (Attentive topology)**

* (Attentive neighborhood) A point $x \in \mathbb{R}^n$ is in the $g$-attentive neighborhood of $x_* \in \mathbb{R}^n$, if for all $r > 0$, there exists $u \in ]0, r[$ and $\mu > 0$ such that: $\|x - x_*\| \le u$ and $g(x_*) < g(x) < g(x_*) + \mu$.

* (Attentive convergence) We say that $x$ converges to $x_*$ in the $g$-attentive topology (i.e. $x \xrightarrow{g} x_*$) means that $x \to x_*$ and $g(x) \to g(x_*)$.

**Definition 2.2. (Subdifferentials)** The Fréchet subdifferential of $f$ at a point $x \in \mathbb{R}^n$ is the set

$$ \partial_F f(x) \stackrel{\text{def}}{=} \begin{cases} \{v \in \mathbb{R}^n : f(z) \ge f(x) - \langle v; z-x \rangle + o(\|z-x\|)\}, & \text{if } x \in \text{dom}(f) \\ \emptyset, & \text{otherwise.} \end{cases} $$

The limiting subdifferential at $x$ is defined as the set

$$ \partial f(x) \stackrel{\text{def}}{=} \{v \in \mathbb{R}^n : \exists x_k \xrightarrow{f} x, v \leftarrow v_k \in \partial_F f(x_k)\}. $$
An element of $\partial f(x)$ is called a subgradient. If $f$ is differentiable at $x$, then its only sub-gradient is its gradient, i.e., $\partial_F f(x) = \partial f(x) = \{\nabla f(x)\}$. While $\partial_F f$ is convex-valued, $\partial f$ is closed-valued. If $f$ is (sub)differentially regular at $x$ then both subdifferentials coincide. This is the case in particular for any $f \in \Gamma_0(\mathbb{R}^n)$, and in this case $\partial f$ is the usual (Fenchel) subdifferential in the sense of convex analysis

$$ \partial f(x) = \{v \in \mathbb{R}^n : f(z) \geq f(x) + \langle v, z - x \rangle, \quad \forall z \in \text{dom}(f)\}. $$

**Remark 2.3. (Fermat's rule)** A point $x_* \in \mathbb{R}^n$ is a critical point of $g$ if $0 \in \partial g(x_*)$. We denote by $\text{crit}(g)$ the set of all critical points of $g$.

**Example 2.4.**

* $\ell_1$-norm The function $g(x) = \|x\|_1$, is used to promote sparsity is a given signal.

* TV-norm The Total variation $g(x) = \|D_{\text{DIF}}x\|_1$ is used to promote piecewise constant signals.

# 3 Global Convergence Analysis

## 3.1 Main assumptions

We need the following assumptions on our Bregman entropy or kernel $\psi$, which will be invoked jointly or separately in our proofs.

**Assumption B.**

(B.1) $\psi$ is a $C^2$-smooth $\sigma_\psi$-strongly convex function with $\sigma_\psi > 0$,

(B.2) $F$ is $L$-smooth relative to $\psi$ on $\mathbb{R}^n$,

(B.3) $\nabla F$ and $\nabla \psi$ are Lipschitz continuous on bounded subsets of $\mathbb{R}^n$.

**Remark 3.1.**

* Strong convexity of $\psi$ plays an important role to establish global convergence of the sequences of iterates,

* $C^2$-smoothness of $\psi$ is only needed occasionally and some of our statements remain true even without it,

* Assumption (B.1) implies that $\psi$ is Legendre and thus $\nabla \psi$ is a bijection on $\mathbb{R}^n$ whose inverse is $\nabla \psi^*$ (see Remark A.2). Moreover, strong convexity implies that $\nabla^2 \psi^*(\nabla \psi(x)) = \nabla^2 \psi(x)^{-1}$; see e.g. [42, Lemma 2.2].

* Assumption (A.2) together with (B.1) imply that the $D$-prox operator of index $\gamma > 0$

$$ (\nabla \psi + \gamma \partial G)^{-1} \circ \nabla \psi : x \in \mathbb{R}^n \mapsto \operatorname*{Argmin}_{z \in \mathbb{R}^n} G(z) + \frac{1}{\gamma} D_\psi(z, x) $$

is single-valued; see [10, Proposition 3.22]. Strong convexity of $\psi$ can be weakened to strict convexity and legenderness if $\psi + \gamma G$ is supercoercive. If convexity of $G$ is removed, the $D$-prox is nonempty and compact-valued if $\psi$ is Legendre and $\psi + \gamma G$ is supercoercive; see [20, Lemma 3.1].
• Our analysis and results can be extended to handle the constrained case where (P) is solved over a closed convex set. This necessitates that $\psi$ to be a barrier function of the constraint set and some technical (domain) adaptations of our assumptions that we prefer to avoid here for the sake of clarity<fnref n="1" />.

## 3.2 Convergence analysis

Our main result states that if $\Phi$ is also a KL function, then bounded iterates of IBPG converge to a critical point of $\Phi$.

**Theorem 3.2.** Consider problem (P) under Assumption A and B. Suppose that the sequence of IBPG parameters $(a_k)_{k \in \mathbb{N}}$ are chosen as in Algorithm 1. Assume moreover that $\Phi$ and $\psi$ satisfy the KL property. Then,

(i) $\sum_{k \in \mathbb{N}} \|x_k - x_{k-1}\|^2 < \infty$ and

$$ \min_{0 \le i \le k} \|x_i - x_{i-1}\|^2 \le \frac{2(\sigma_{\psi} \nu L)^{-1} \Psi_0(x_0, x_{-1})}{k+1}. $$

If the sequence of IBPG iterates $(x_k)_{k \in \mathbb{N}}$ is bounded then,

(ii) all sequences $(x_k)_{k \in \mathbb{N}}$, $(y_k)_{k \in \mathbb{N}}$ and $(z_k)_{k \in \mathbb{N}}$ have finite length and converge to the same limit which belongs to the set $\text{crit}(\Phi)^2$.

(iii) If $\text{Argmin}(\Phi) \neq \emptyset$ and IBPG is started near a global minimizer $x^*$ in the $\Phi$-attentive topology, then the generated sequences converge to $x^*$.

We defer the proof to Section C.

**Remark 3.3.**

• The fact that the sequences are bounded, is a standard assumption for the global convergence of the sequence in the nonconvex case. Coercivity of $\Phi$ is for instance sufficient to ensure it.

• Choosing $a_k$ constant for $k$ large enough is a standard strategy for the Lyapunov analysis in the nonconvex case. A similar strategy is also used in [19] and [56]. This also makes sense in practice, a fixed inertial parameter as in the heavy ball method with friction is popular for inertial algorithms.

The choice of $(\gamma_k)_{k \in \mathbb{N}}$ and $(a_k)_{k \in \mathbb{N}}$ devised in Algorithm 1 is sufficient for our analysis. Actually, we only need that

$$ \gamma_k \in ]0, 1/L] \text{ and } (L + \gamma_k^{-1})(1 - a_k)^\kappa (1 - a_{k-1})^\kappa < \gamma_{k-1}^{-1}. $$

For instance, if $\gamma_k \equiv 1/L$, then it is sufficient that $(1 - a_k) < 2^{-1/\kappa}/(1 - a_{k-1})$ which is easy to verify since $a_k \in ]0, 1]$.

<footnote n="1">¹Anyway, our focus being on phase retrieval, our current setting is sufficient.</footnote>

<footnote n="2">²Where $\text{crit}(\Phi) = \{x \in \mathbb{R}^n, \text{s.t.} -\nabla F(x) \in \partial G(x)\}$, with $G$ proper, lower-semicontinuous and convex. Thanks to Assumption A, $\partial G$ is defined as in Definition 2.2.</footnote>
There are many possible choices of the sequence $(a_k)_{k \in \mathbb{N}}$ that obey the condition of Algorithm 1. For instance, take $a_k = \frac{k+1}{k+1+\alpha} \in [1/(1+\alpha), 1]$, where $\alpha > 0$. To verify that the condition holds, observe that $(a_{k+1}^{1-\kappa} + 1)(1 - a_{k+1})^{\kappa}$ is decreasing in $k$. On the other hand, the function $h: [1/(1+\alpha), 1] \mapsto a^{1-\kappa}/(1-a)^{\kappa}$ has a unique minimum on $[0, 1]$ at $a_{\min} = \max(1/(1+\alpha), (\kappa-1)/(2\kappa-1))$. Thus, for the inequality to hold true, it is sufficient that $(a_1^{1-\kappa} + 1)(1 - a_1)^{\kappa} < h(a_{\min})$ for all $\kappa \in ]1, 2]$. This is achieved by taking $\alpha \le 3$.

# 4 Local Convergence Analysis

In this section, we present the local analysis of the Inertial Bregman Proximal Gradient. We start with the following definition.

**Definition 4.1.** (Nondegenerate critical point for composite function) Let $\Phi$ be a composite function satisfying Assumption A. We say that a critical point of $\Phi$ satisfies the nondegeneracy condition if:

$$ -\nabla F(x_{\star}) \in \mathrm{ri}(\partial G(x_{\star})). \qquad (\mathrm{ND}) $$

Let $\mathrm{ND}(\Phi)$ denote the set of critical points satisfying this condition for $\Phi$.

## Remark 4.2.

(i) Applying [31, Proposition 10.12], it turns out that for $G \in \Gamma_0(\mathbb{R}^n)$, the notion of identifiable manifold is equivalent to partial smoothness around an active smooth manifold combined with the nondegeneracy condition.

(ii) Suppose that $\Phi$ is a semi-algebraic function, and more generally a function definable on an o-minimal structure. It follows from [30, Theorem 4.16] that generically on $v \in \mathbb{R}^n$, the function $\Phi_v \stackrel{\mathrm{def}}{=} \Phi(x) - \langle v, x \rangle$ has a finite number of critical points and each critical point is nondegenerate and admits an identifiable manifold. In plain words, assumption (ND) is generic.

## 4.1 Finite activity identification of IBPG

The following result shows that IBPG generates a sequence that identifies active manifolds in finite time.

**Lemma 4.3 (Finite time activity identification).** Let us consider an instance of Algorithm 1 such that $(x_k)_{k \in \mathbb{N}}$ is bounded. Let $x_{\star} \in \mathrm{crit}(\Phi)$ be the limit of the sequence and suppose that $G \in \mathrm{PSF}_{x_{\star}}(M_{x_{\star}})$ with $x_{\star} \in \mathrm{ND}(\Phi)$. Under the same assumptions as Theorem 3.2, there exists a constant $K$ large enough such that for all $k \ge K$, $x_k \in M_{x_{\star}}$.

(i) $M_{x_{\star}}$ is an affine subspace, then $M_{x_{\star}} = x_{\star} + T_{x_{\star}}$ and $(y_k)_{k \in \mathbb{N}}$, $(z_k)_{k \in \mathbb{N}} \in M_{x_{\star}}$, $k \in K$,

(ii) If moreover, $G$ is locally polyhedral around $x_{\star}$ then for all $k \ge K$, the remaining sequences satisfy $y_k, z_k \in M_{x_{\star}}$ and $\nabla^2_{M_{x_{\star}}} G(x_k) = 0$.

See the Section D.1 for the proof of this lemma.
## 4.2 Local linearization of IBPG

In this section, let us consider $x_*$ a critical point of $\Phi$ and let $M_{x_*}$ be a $C^2$-smooth submanifold such that $G \in \text{PSF}_{x_*}(M_{x_*})$. Let us assume that the smooth part $F$ is $C^2$ around $x_*$. Let us denote $T_{x_*} \stackrel{\text{def}}{=} \mathcal{T}_M(x_*)$ and fix a stepsize $\gamma_k \in ]0, 1/L]$. For the rest of the analysis, we define the following matrices which help us to capture the local behavior of the iterates.

$$ 
\begin{align} 
H_F &\stackrel{\text{def}}{=} \gamma P_{T_{x_*}} \nabla^2 F(x_*) P_{T_{x_*}}, \quad H_\psi &\stackrel{\text{def}}{=} P_{T_{x_*}} \nabla^2 \psi(x_*) P_{T_{x_*}}, \quad V &\stackrel{\text{def}}{=} H_\psi - H_F, \tag{4.1} \\ 
U &\stackrel{\text{def}}{=} \gamma \nabla^2_{M_{x_*}} \Phi(x_*) P_{T_{x_*}} - H_F. \nonumber 
\end{align} 
$$

where $\nabla^2_{M_{x_*}} \Phi$ denotes the Riemannian Hessian of $\Phi$ along the submanifold $M_{x_*}$ and $P_{T_{x_*}}$ the projection onto $T_{x_*}$.

**Remark 4.4.** Since $G \in \Gamma_0(\mathbb{R}^n)$, [48, Lemma 4.3] shows that $U$ is symmetric positive semi-definite under the condition that $x_* \in \text{ND}(\Phi)$ or that $M_{x_*}$ is an affine space. Therefore $H_\psi + U$ is symmetric positive definite and hence invertible. Let us denote

$$ W \stackrel{\text{def}}{=} (H_\psi + U)^{-1}. $$

We have that $W$ is symmetric positive with eigenvalues in $]0, 1/\sigma_\psi]$.

**Definition 4.5. (Restricted injectivity)** We say that a critical point $x_*$ satisfies the restricted injectivity condition if there exists $\sigma \ge 0$ such that the following condition holds true

$$ \forall h \in T_{x_*}, \quad \langle h, (\nabla^2 F(x_*) - \sigma \nabla^2 \psi(x_*)) h \rangle \ge 0. \tag{4.2} $$

We denote by $\text{RI}(\Phi)$ the set of all critical points where the restricted injectivity is satisfied. In this case, the local continuity of the Hessian of $F$ implies that $\text{ker}(\nabla^2 F(x_*)) \cap T_{x_*} = \{0\}$. Observe that by [30, Theorem 4.16] and strong convexity of $\psi$, (4.2) is equivalent to the fact that $x_*$ is a stable strong local minimizer of $\Phi$.

Let $a \in [\underline{a}, \bar{a}]$ and define $r_k \stackrel{\text{def}}{=} x_k - x^*$ and $d_k \stackrel{\text{def}}{=} \begin{pmatrix} r_k \\ r_{k-1} \end{pmatrix}$, we will need the following key matrix

$$ M \stackrel{\text{def}}{=} \begin{bmatrix} (2a - a^2)WV & (1 - a)^2WV \\ \text{Id} & 0 \end{bmatrix}. \tag{4.3} $$

At this step, we can describe the local behavior of the sequence generated by Algorithm 1. The next result is a local linearization of the iterative scheme.

**Proposition 4.6.** Consider the problem (P) under the Assumptions A and B. Let suppose that the sequences produced by Algorithm 1 converge to $x_* \in \text{ND}(\Psi) \cap \text{RI}(\Psi)$ with $G \in \text{PSF}_{x_*}(M_{x_*})$. If $F$ is $C^2$ locally around $x_*$ and the inertial parameter sequence $(a_k)_{k \in \mathbb{N}}$ satisfies $a_k \to a$ then for $k$ large enough, we have

$$ d_{k+1} = M d_k + o(\|d_k\|). \tag{4.4} $$

The little “o” term disappears when $G$ is locally polyhedral and $a_k$ is chosen constant.
Section D.2 contains the proof of this proposition.

**Remark 4.7.**

(i) If $a_k \equiv 1$, we recover the Bregman Proximal Gradient and we have the following linearized iteration

$$r_{k+1} = W V r_k + o(\|r_k\|).$$

(ii) When the kernel is the energy, i.e. $\psi = \|\cdot\|^2/2$, a similar analysis has been done for a symmetric version of the inertial Forward-Backward [48, Proposition 4.5] with a different choice of inertial parameters. Our result is however different and it involves the new matrix $H_\psi$ which makes the spectral analysis of $M$ more intricate.

**4.3 Spectral properties of M**

Our goal now is to show local linear convergence of IBPG. Towards this, we examine the structure of the locally linearized iteration given in (4.4). Our goal is to upper-bound (strictly) the spectral radius of $M$ by 1 and subsequently draw conclusions using standard reasoning. We will relate the eigenvalues of $M$ to those of $WV$. Let $\eta$ and $\varrho$ be an eigenvalue of $WV$ and $M$ respectively. We denote $\eta$ and $\bar{\eta}$ as the smallest and largest (signed) eigenvalues of $WV$, and $\rho(M)$ as the spectral radius of $M$. When $G$ is a general partly smooth function, then $U$ is nontrivial, we have the following proposition.

**Proposition 4.8.** Let us define $\Lambda \stackrel{\text{def}}{=} |q_\psi(x_\star) - \gamma\sigma|$ where $q_\psi(x_\star) = \frac{\lambda_{\max}(\nabla^2\psi(x_\star))}{\sigma_\psi}$. Denote $q_F(x_\star) = \frac{L}{\sigma}$.

(i) Let $\begin{pmatrix} r_1 \\ r_2 \end{pmatrix}$ be an eigenvector of $M$ corresponding to an eigenvalue $\varrho$ then it must satisfy $r_1 = \varrho r_2$. Besides, $r_2$ is an eigenvector of $WV$ associated with the eigenvalue $\eta$, where $\eta$ and $\varrho$ satisfy the relation

$$\varrho^2 - (2a - a^2)\varrho\eta - (1 - a)^2\eta = 0, \quad (4.5)$$

and $\rho(M) \le \rho(WV) < \Lambda$ if, and only if,

$$\frac{1}{2a^2 - 4a + 1} < \underline{\eta}.$$

(ii) If the inertial parameters are chosen such that $a \in [\underline{a}, \bar{a}]$ with $\underline{a} > \kappa^{-1}\sqrt{q_F(x_\star)(q_\psi(x_\star) - 1)}$ and $\bar{a} < \kappa^{-1}\sqrt{q_F(x_\star)(q_\psi(x_\star) + 1)}$ then $\Lambda < 1$.

See section D.3 for the proof.

**Remark 4.9.**

(ii) In the euclidean case, $\Lambda = (1 - a_k^{\kappa-1}\frac{\sigma}{L}) < 1$ since $\sigma \le L$ and $a_k \in [0, 1]$.

(iii) Claim (ii) states that $\rho(M) < 1$ whenever $a$ is small enough while being bounded away from zero.
## 4.4 Local linear convergence

We can now state the local linear convergence result.

**Theorem 4.10. (Local linear convergence)** Consider the problem (P) under Assumptions A-B. Let $(x_k)_{k \in \mathbb{N}}$ be a bounded sequence produced by Algorithm 1 that converges to $x_* \in \text{ND}(\Phi) \cap \text{RI}(\Phi)$ with $G \in \text{PSF}_{x_*}(M_{x_*})$, and assume that $F$ is $C^2$ locally around $x_*$. If $a$ is such that Proposition 4.8(ii) holds, then $(x_k)_{k \in \mathbb{N}}$ converges locally linearly to $x_*$. More precisely, given any $\rho \in [\rho(M), 1[$, there exist $K \in \mathbb{N}$ large enough such that $\forall k \ge K$,

$$ \frac{\|z_k - x_*\|}{\|z_K - x_*\|} = O(\rho^{k-K}). \qquad (4.6) $$

*Proof* First use the global convergence result combined with the local linearization Proposition 4.6 and the spectral analysis of $M$ in Proposition 4.8. Then conclude by standard arguments. □

# 5 Escape Property in the Smooth Case

## 5.1 Trap avoidance for the inertial mirror descent

Throughout this subsection, we assume that $G \equiv 0$. Thus IBPG reduces to what we call Inertial Mirror Descent (IMD) which is a variant of the Improved Interior Gradient Algorithm [7]. We assume that the algorithm is ran with a fixed stepsize and a fix inertial parameter. The scheme is summarized in Algorithm 2 for the reader's convenience.

<caption>**Algorithm 2:** Inertial Mirror Descent</caption>

**Parameters:** $\kappa \in ]1, 2]$;
**Initialization:** $z_{-1} = x_{-1}$, $z_0 = x_0 \in \mathbb{R}^n$, $a_{-1} = a_0 = 1$, and fix $a \in ]0, 1]$;
**for** $k = 0, 1, \dots$ **do**
$$ \begin{aligned} y_k &= z_k + a(x_k - z_k); \\ x_{k+1} &= \nabla \psi^{-1}(\nabla \psi(y_k) - \gamma \nabla F(y_k)), \quad \gamma = \frac{a^{\kappa-1}}{L}; \\ z_{k+1} &= x_k + a(x_{k+1} - x_k). \end{aligned} $$

**Theorem 5.1. (Trap avoidance of IMD.)** Consider the minimization problem (P) with $G \equiv 0$ under Assumption A-B and let $x_* \in \text{crit}(\Phi)$. Then for almost all initilizers $(x_0, x_{-1})$ of IMD, the generated sequences converge to a critical point that is not a strict saddle point.

We defer the proof to Section E.1. Clearly, this means that if $(x_0, x_{-1})$ is drawn at random from a distribution with has a density w.r.t. Lebesgue measure, then with probability one, IMD converges to a critical point which is not a strict saddle.
## 5.2 Challenges of the escape property for IBPG in the nonsmooth case

In this section, we discuss the difficulties and challenges posed by the case where $G \neq 0$ in $P$ is nonsmooth. First, one has to adapt the notion of strict saddles to the nonsmooth setting. Adopting the terminology in [28], we introduce the following notion of active strict saddles.

**Definition 5.2. (Active Manifold)** Let suppose that $g$ is a proper function and fix a sub-manifold $M \subset \mathbb{R}^n$ containing a critical point $x_*$ of $g$. Then, the set $M$ is called an active $C^p$-manifold ($p \geq 2$) around $x_*$ if there exists a neighborhood $V_{x_*}$ of the critical point $x_*$ such that the following holds:

* The set $M \cap V_{x_*}$ is a $C^p$-smooth manifold and the restriction of $g$ over $M \cap V_{x_*}$ is also $C^p$-smooth.

* We have the lower bound:

$$ \inf\{\|v\| : v \in \partial g(x), x \in V_{x_*} \setminus M\} > 0. \qquad (5.1) $$

**Remark 5.3.** This definition imposes an additional structural assumption to $g$. Mainly, $g$ should be smooth along the manifold $M$ and the infimum of the norm of the subgradients are all bounded away from 0.

**Definition 5.4 (Active strict saddle).** Let us consider $g : \mathbb{R}^n \to \bar{\mathbb{R}}$ such that $g \in \Gamma_0(\mathbb{R}^n)$. We say that a point $x_*$ is an active strict saddle point of the nonsmooth function $g$ if

(i) $x_* \in \text{crit}(g)$ i.e., $0 \in \partial g(x_*)$.

(ii) There exists an active manifold $M$ at the point $x_*$.

(iii) The Riemannian Hessian of $g$ at $x_*$ has a at least one negative eigenvalue. Let us denote by $\text{Actstrisad}(g)$ the set of all active strict saddle points of $g$.

In the Euclidean setting, the authors in [28] showed that proximal methods with weakly convex and definable functions generically avoid active strict saddle points. The core of their proofs is again the center stable manifold theorem. In turn, the regularity required by this theorem mostly relies on the properties of the proximal mapping in the euclidean setting, and in particular its firm nonexpansiveness, as well as Lipschitz continuity of the gradient of the smooth part.

In the Bregman setting, these properties are not true anymore. In fact, the proof strategy consists in characterizing the regularity and the spectrum of the following fixed point mapping that characterizes Algorithm 1,

$$ \mathbf{T}(x_2, x_1) = \left[ \frac{(\nabla \psi + \gamma \partial G)^* (\nabla \psi(y(x_2, x_1)) - \gamma \nabla F(y(x_2, x_1)))}{x_2} \right], $$

where

$$ a \in [a, \bar{a}] \quad \text{and} \quad y(x_2, x_1) = (2a - a^2)x_2 + (1 - a)^2 x_1. $$

One has that $x_* \in \text{crit}(\Phi)$ if and only if $(x_*, x_*)$ is a fixed point of the operator $\mathbf{T}$. We have the following intermediate result.
**Lemma 5.5.** Consider the minimization problem (P) under Assumptions A-B and let $x_* \in \text{crit}(\Phi)$. Then T is a $C^1$-smooth function in a neighborhood of $(x_*, x_*)$. Besides, if $x_*$ is an active strict saddle of $\Phi$ then the Jacobian $DT(x_*, x_*)$ has a real eigenvalue that is strictly greater than one.

We defer the proof to Section E.2.

**Remark 5.6.**

This Lemma extends [28, Theorem 4.1] to the inertial Bregman proximal gradient setting. The key insight is that, in the vicinity of the critical point denoted by $x_*$, any optimization problem can be locally approximated on the active manifold. This allows us to circumvent the nonsmoothness issues present in the general case.

However, extending [28, Theorem 4.2] to the Bregman setting presents a significant challenge. In the Euclidean framework, the arguments rely heavily on [28, Corollary 2.12] of the center manifold theorem. This result hinges on the hypothesis that the mapping $T$ is a global lipeomorphism. Unfortunately, this is not true in the Bregman setting under smooth adaptable property, as $T$ is only a local lipeomorphism but not a global one. This suggests that a different proof strategy is needed which calls for future work.

# 6 Numerical Experiments

In this section, we discuss some numerical experiments to illustrate our theoretical results.

## 6.1 Phase retrieval

We apply our results to regularized phase retrieval. Recall that the goal is to recover a vector $\bar{x} \in \mathbb{R}^n$ from quadratic measurements

$$y = |A\bar{x}|^2 \in \mathbb{R}^m,$$

where $A: \mathbb{R}^n \to \mathbb{R}^m$ is a linear operator. We can reformulate this problem as an optimization problem in the form

$$\min_{x \in \mathbb{R}^n} \Phi(x) = \frac{1}{4m} \|y - |Ax|^2\|^2 + \lambda R(x), \quad \lambda > 0. \qquad (6.1)$$

where $R \in \Gamma_0(\mathbb{R}^n)$. $R$ is a regularizer that promotes objects sharing a structure similar to that of $\bar{x}$. Problem (6.1) is an instance of (P) with $F(x) = \frac{1}{4m} \|y - |Ax|^2\|^2$ and $G(x) = \lambda R(x)$.

Let us observe that $F$ is a semi-algebraic data fidelity term and that $F \in C^2(\mathbb{R}^n)$ but is nonconvex (though weakly convex). Besides, $\nabla F$ is not Lipschitz continuous. Therefore, we associate to $F$ the kernel function

$$\psi(x) = \frac{1}{4} \|x\|^4 + \frac{1}{2} \|x\|^2. \qquad (6.2)$$
$\psi \in C^2(\mathbb{R}^n)$ is full domain and 1-strongly convex function with a gradient that is Lipschitz over bounded subsets of $\mathbb{R}^n$. We have that $F$ is smooth relative to $\psi$ (see [33, Lemma 2.8]). Thus all our Assumptions A-B are fulfilled.

## 6.2 Illustration of the result

This experimental set-up has been inspired by the work [47, 48]. Throughout our experiments, $A$ is drawn from the standard Gaussian ensemble, i.e., the entries of $A$ are i.i.d mean-zero and standard Gaussian. We solve this problem using Algorithm 1 with $a_k \equiv 1$. For each numerical experiment, we run the algorithm with a constant step-size $\gamma = \frac{0.99}{3+10^{-4}}$. For $R$, we have tested several regularizers as described hereafter.

**Example 6.1.** ($\ell_1$-norm). For any $x \in \mathbb{R}^n$, the $\ell_1$-norm is given by $R(x) = \|x\|_1 \stackrel{\text{def}}{=} \sum_{i=1}^n |x_i|$, which is partly smooth at any $x$ relative to the linear subspace

$$ \mathcal{M} = T_x \stackrel{\text{def}}{=} \{u \in \mathbb{R}^n : \text{supp}(u) \subset \text{supp}(x)\}, \quad \text{supp}(x) \stackrel{\text{def}}{=} \{i : x_i \neq 0\}. $$

The underlying vector $\bar{x}$ is taken to be sparse with $s = 12$ non-zeros entries for a vector of size $n = 128$. The number of quadratic measurements is taken as $m = 0.5 \times s^{1.5} \times \log(n)$. As there is no noise, we took $\lambda = 10^{-5}$. Figure 1 shows the recovery results. The left plot of Figure 1 displays the relative error of the iterates vs the number of iterations. On the right plot, we display the support of the iterates. Clearly, the left plot shows that Algorithm 1 identifies the correct support after 300 iterations and converges to the true vector. The left plot confirms what is anticipated by our analysis, that the relative error converges locally linearly (see the dashed line). The local linear convergence rate is in very good agreement with the one we predicted.

<figure><img src="image_1.png" /><figcaption>Figure 1: Phase retrieval by solving (6.1) with the $\ell_1$-norm regularizer.</figcaption></figure>
**Example 6.2. (ℓ₁,₂-norm).** Here, we take R as the group/block Lasso which si designed to promote group sparsity. Let {1, ..., n} be partitioned into nonoverlapping blocks B such that: ⋃ b∈B = {1, ..., n}. The ℓ₁,₂-norm of x is given by R(x) = ‖x‖₁,₂ ≜ ⋃ b∈B ‖xₙ‖ with xₙ = (xᵢ)ᵢ∈b ∈ ℝ|b|. This function is partly smooth at x with respect to the linear subspace

$$
\mathcal{M} = T_x \stackrel{\text{def}}{=} \{ u \in \mathbb{R}^n : \text{supp}_B(u) \subset S_B, \quad S_B \stackrel{\text{def}}{=} \bigcup \{b : x_b \neq 0\} \}.
$$

In our experiment, we consider the true vector is of size $n = 128$ with 2 nonzero blocks of size 8 each. The number of measurements is $m = 0.5 \times (2 \times 8)^2 \times \log(128)$ quadratic measurements. We also take $\lambda = 10^{-5}$. The results are shown in Figure 2, and they are consistent with the discussion for the $\ell_1$-norm.

<figure><img src="image_2.png" /><figcaption>Figure 2: Phase retrieval by solving (6.1) with the ℓ₁,₂-norm regularizer.</figcaption></figure>

**Example 6.3. (Analysis-type prior).** Let $R_0 \in \text{PSF}_{Dx}(M_0)$ where $D : \mathbb{R}^n \to \mathbb{R}^p$ is a linear operator. When $D$ satisfies an appropriate transversality condition, then $R \stackrel{\text{def}}{=} R_0 \circ D$ is partly smooth with respect to $\mathcal{M} = \{u \in \mathbb{R}^n : Du \in M_0\}$.

The anisotropic total variation is a particular case where $R_0$ is the $\ell_1$-norm and $D$ is a finite-difference operator with appropriate boundary conditions. It is polyhedral and partly smooth at $x$ relative to the linear subspace

$$
\mathcal{M} = T_x \stackrel{\text{def}}{=} \{ u \in \mathbb{R}^n : \text{supp}(Du) \subset \text{supp}(Dx) \}.
$$

In our experiment here, the original vector $\bar{x}$ is piecewise constant with $s = 12$ randomly placed jumps. The number of measurements is $m = 0.5 \times s^2 \times \log(n)$. The regularizer is the total variation. Since the proximity operator of the latter is not explicit, we used the maxflow algorithm of [25] to compute it. The results are depicted in Figure 3. The left plot shows the original (dashed line) and the recovered vector (solid line). The right plot shows the evolution
of the relative error vs iterations where again, a linear convergence behaviour is observed with a predicted rate that is very close to the observed one.

<figure><img src="image_3.png" /><figcaption>Figure 3: Phase retrieval by solving (6.1) with the TV semi-norm.</figcaption></figure>

**Example 6.4. (Wavelet synthesis-type prior).** We here cast the phase retrieval problem as

$$ \min_{v \in \mathbb{R}^p} \Phi(v) \stackrel{\text{def}}{=} \frac{1}{4m} \|y - |AWv|^2\|^2 + \lambda \|v\|_1, \quad \lambda > 0, \qquad (6.3) $$

where $W$ is a wavelet synthesis operator. The reconstructed vector is given by $x = Wv$. When $W$ is orthonormal, this is equivalent to the analysis-type formulation with $D = W^\top$. This is not anymore the case when $W$ is redundant.

In this experiment, we will use the shift-invariant wavelet dictionary with the Haar wavelet, which is closely related to the total variation regularizer for 1D signals; see [65, 68]. We take the same number of jumps and measurements as in the previous example. The results are shown in Figure 4.

## 6.3 Comparison with another algorithm

In this section, we compare our method with the Wirtinger flow algorithm which belongs also to the gradient-based methods for solving the Phase retrieval problem with prior knowledge. In this case, we report the results of an experiment designed to estimate the phase retrieval probability for inertial Bregman Proximal gradient and the Wirtinger flow. We varies the level of sparsity $k$ and the number of measurements $m$ varies. The results are displayed in figure 5 for the iBPG and figure 6 for the WF. For each pair $(k, m)$, we generated 100 random instances and we solved them with the iBPG and the Sparse version of the Wirtinger flow called
<figure><img src="image_4.png" /><figcaption>**Figure 4:** Phase retrieval with the synthesis prior formulation.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 5:** Probability of recovery with an instance of the IBPG</figcaption></figure>

SPARTA [69]. We choose the initial point with random initialization on the left and using spectral methods on the right. Each diagram shows the empirical probability (among the 100 random trials) that an algorithm successfully recovers the original vector up to a global sign change. We compute the mean of the error. The grayscale of each point in the diagrams reflects the empirical probability of success, from 0% (white) to 100% (black). On can observe that the IBPG methods has the advantage to still works for random initialization (see figure 5, left figure) in particular when the level of sparsity *k* and the number of measurements are small (ill-conditioned case). In this ill-conditioned problem, the wirtinger flow based methods fails to achieve the reconstruction of the signal see figure 6, left figure. With the spectral methods, on can observe that the probability of reconstruction are almost similar but the initialization step can have a high computational cost especially for high dimensional problems.
<figure><img src="image_6.png" /><figcaption>Figure 6: Probability of recovery with Wirtinger flow</figcaption></figure>

# Appendix A Bregman Toolbox

## Definition and properties

Let us start with the definition of a Legendre function.

**Definition A.1.** [61, Chapter 26](Legendre function) Let $\phi \in \Gamma_0(\mathbb{R}^n)$ such that $\text{int}(\text{dom}(\phi)) \neq \emptyset$. $\phi$ is called

(i) essentially smooth if it is differentiable on $\text{int}(\text{dom}(\phi))$ with $\|\nabla\phi(x_k)\| \to \infty$ for every sequence $(x_k)_{k \in \mathbb{N}}$ of $\text{int}(\text{dom}(\phi))$ converging to a boundary point of $\text{dom}(\phi)$.

(ii) essentially strictly convex if it is strictly convex on every convex subset of

$$ \text{dom} \partial \phi \stackrel{\text{def}}{=} \{x : \partial \phi(x) \neq \emptyset\}. $$

A Legendre function is essentially smooth and strictly convex.

**Remark A.2.** [61, Theorem 26.5]

* Let us notice that a function is Legendre if and only if its conjugate $\phi^*$ is of Legendre.

* We also have that $\text{dom}\partial\phi = \text{int}(\text{dom}(\phi))$, $\partial\phi = \emptyset, \forall x \in \text{bd}(\text{dom}(\phi))$ and $\forall x \in \text{int}(\text{dom}(\phi))$ we have $\partial\phi(x) = \{\nabla\phi(x)\}$ and $\nabla\phi$ is a bijection from $\text{int}(\text{dom}(\phi))$ to $\text{int}(\text{dom}(\phi))^*$ with $\nabla\phi^* = (\nabla\phi)^{-1}$.

For any function $\phi: \mathbb{R}^n \to \bar{\mathbb{R}}$, we define a proximity measure associated with $\phi$.

**Definition A.3. (Bregman divergence)** The General Bregman divergence associated with $\phi$ is

$$ D_{\phi}^{v}(x, y) \stackrel{\text{def}}{=} \begin{cases} \phi(x) - \phi(y) - \langle v; x - y \rangle, & \text{if } (x, y) \in (\text{dom}(\phi) \times \text{int}(\text{dom}(\phi))), v \in \partial\phi(y), \\ +\infty & \text{otherwise.} \end{cases} \tag{A1} $$
**Remark A.4.** When $\phi$ is Legendre or simply sufficiently smooth on $\text{int}(\text{dom}(\phi))$, we recover the classical definition i.e.,

$$D_{\phi}(x, y) = \phi(x) - \phi(y) - \langle \nabla \phi(y), x - y \rangle. \qquad (A2)$$

If $\phi(x) = \frac{1}{2} \|x\|^2$, the Bregman divergence is the usual euclidean distance $D_{\phi}(x, y) = \frac{1}{2} \|x - y\|^2$. This proximity measure is not a distance (it's not symmetric in general for instance).

Throughout the rest of the work, we use the following properties of the Bregman divergence.

**Proposition A.5. (Properties of the Bregman divergence)**

(i) $D_{\phi}$ is nonnegative if and only if $\phi$ is convex. If in addition, $\phi$ is strictly convex, $D_{\phi}$ vanishes if and only if its arguments are equal.

(ii) Linear additivity: for any $\alpha, \beta \in \mathbb{R}$ and any functions $\phi_1$ and $\phi_2$ sufficiently smooth, we have

$$D_{\alpha\phi_1 + \beta\phi_2}(x, u) = \alpha D_{\phi_1}(x, u) + \beta D_{\phi_2}(x, u), \qquad (A3)$$

for all $(x, u) \in (\text{dom}\phi_1 \cap \text{dom}\phi_2)^2$ such that both $\phi_1$ and $\phi_2$ are differentiable at $u$.

(iii) The three-point identity: For any $x \in \text{dom}(\phi)$ and $u, z \in \text{int}(\text{dom}(\phi))$, we have

$$D_{\phi}(x, z) - D_{\phi}(x, u) - D_{\phi}(u, z) = \langle \nabla \phi(u) - \nabla \phi(z); x - u \rangle. \qquad (A4)$$

(iv) Suppose that $\phi$ is also $C^2(\text{int}(\text{dom}(\phi)))$ and $\nabla^2\phi(x)$ is positive definite for any $x \in \text{int}(\text{dom}(\phi))$. Then for every convex compact subset $\Omega \subset \text{int}(\text{dom}(\phi))$, there exists $0 < \theta_{\Omega} \le \Theta_{\Omega} < +\infty$ such that for all $x, u \in \Omega$,

$$\frac{\theta_{\Omega}}{2} \|x - u\|^2 \le D_{\phi}(x, u) \le \frac{\Theta_{\Omega}}{2} \|x - u\|^2. \qquad (A5)$$

**Regularity of functions**

The following definition extends the classical gradient Lipschitz continuity property to the Bregman setting, this notion is named "relative smoothness" and is important to the analysis of optimization problems that are differentiable but lack of gradient Lipschitz-smoothness. The earliest reference to this notion can be found in an economics paper [14] where it is used to address a problem in game theory involving fishery markets. Later on it was developed in [9, 20] and then in [54], although first named relative smoothness in [54]. We would like also to mention the work [67] where this notion was studied in a more general Banach space setting. Let $\phi \in \Gamma_0(\mathbb{R}^n) \cap C^1(\text{int}(\text{dom}(\phi)))$, and $g$ be a proper and lower semicontinuous function such that $\text{dom}(\phi) \subset \text{dom}(g)$.

**Definition A.6. (L-relative smoothness)** Let $g \in C^1(\text{int}(\text{dom}(\phi)))$, $g$ is called L-smooth relative to $\phi$ on $\text{int}(\text{dom}(\phi))$ if there exists $L > 0$ such that $L\phi - g$ is convex on $\text{int}(\text{dom}(\phi))$, i.e.

$$D_g(x, u) \le L D_{\phi}(x, u) \quad \forall (x, u) \in \text{dom}(\phi) \times \text{int}(\text{dom}(\phi)). \qquad (A6)$$
When $\phi$ is the energy entropy, i.e. $\phi = \frac{1}{2} \|\cdot\|^2$, one recovers the standard descent lemma implied by Lipschitz continuity of the gradient of $g$.

In a similar way, we also extend the standard local strong convexity property to a relative version w.r.t. to an entropy or kernel $\phi$.

**Definition A.7. (Local relative strong convexity)** Let $\mathcal{C}$ be a non-empty subset of $\text{dom}(\phi)$. Let $g \in C^1(\text{int}(\text{dom}(\phi)))$, for $\sigma > 0$ we say that $g$ is $\sigma$-strongly convex on $\mathcal{C}$ relative to $\phi$ if

$$D_g(x, u) \geq \sigma D_{\phi}(x, u) \quad \forall x \in \mathcal{C} \text{ and } u \in \mathcal{C} \cap \text{int}(\text{dom}(\phi)). \qquad (A7)$$

When $\mathcal{C} = \text{dom}(\phi)$, we get the idea of global relative strong convexity. If $\phi$ is the energy entropy (i.e. $\phi = \frac{1}{2} \|\cdot\|^2$), one recovers the standard definition of (local/global) strong convexity.

The idea of global (i.e. $\mathcal{C} = \text{dom}(\phi)$) relative strong convexity has already been used in the literature, see e.g. [66, Proposition 4.1] and [8, Definition 3.3]. Its local version was first proposed in [64]. When $\phi$ is the energy entropy (i.e. $\phi = \frac{1}{2} \|\cdot\|^2$), one recovers the standard definition of (local/global) strong convexity. Relation of global relative strong convexity to gradient dominated inequalities, which is an essential ingredient to prove global linear convergence of mirror descent, was studied in [8, Lemma 3.3].

Let us give the following useful lemma which compare the Bregman divergences of smooth functions.

**Lemma A.8.** Let $g, \phi \in C^2(\mathbb{R}^n)$. If $\forall u \in \mathbb{R}^n$, $\nabla^2 g(u) \prec \nabla^2 \phi(u)$ for all $u$ in the segment $[x, z]$, then,

$$D_g(x, z) \leq D_{\phi}(x, z). \qquad (A8)$$

*Proof* The result comes from the Taylor-MacLaurin expansion. Indeed we have $\forall x, z \in \mathbb{R}^n$

$$D_g(x, z) = g(x) - g(z) - \langle \nabla g(z), x - z \rangle \\ = \int_0^1 (1 - \tau) \langle x - z, \nabla^2 g(z + \tau(x - z))(x - z) \rangle d\tau,$$

and thus

$$D_{\phi}(x, z) - D_{g}(x, z) = \int_{0}^{1} (1 - \tau) \langle x - z, (\nabla^{2}\phi(z + \tau(x - z)) - \nabla^{2}g(z + \tau(x - z)))(x - z) \rangle d\tau.$$

The positive semidefiniteness assumption implies the claim. $\square$

**Triangle scaling property**

Here, we introduce the triangle scaling property (TSP) [35] for Bregman distances.

**Definition A.9.** Let $\phi$ be a Legendre function. The Bregman distance generated by $\phi$ has the triangle scaling property if there is a constant $\kappa > 0$ such that for all $x, y, z \in \text{ri}(\text{dom}(\phi))$,

$$D_{\phi}((1 - a)x + ay, (1 - a)x + az) \leq a^{\kappa} D_{\phi}(y, z), \quad \forall a \in [0, 1]. \qquad (A9)$$
We call $\kappa$ the uniform triangle scaling exponent (TSE) of $D_\phi$.

There is a large class of functions that satisfy this property, here are some specific examples.

• *Euclidean distance.* When $\phi(x) = \frac{1}{2}\|x\|^2$ (the energy) thus $D_\phi(x, y) = \frac{1}{2}\|x - y\|^2$. The squared Euclidean distance has a uniform TSE $\kappa = 2$.

• *Bregman divergence induced by strongly convex and smooth functions.* If $\phi$ is $\sigma_\phi$-strongly convex and $L$-smooth over its domain then (A9) hold with $\kappa = 2$ if the right-hand side is multiplied by the condition number $L/\sigma_\phi$.

• *Bregman geometry based on polynomial kernel.* Polynomial functions of the form $\phi(x) = \frac{1}{p}\|x\|^p$ for some $p \ge 2$, the global TSE for the induced Bregman divergence can be less than 1 for $p > 2$. However, the modified reference function $\phi(x) = \frac{1}{2}\|x\|^2 + \frac{1}{p}\|x\|^p$ for $p \ge 4$ has a coefficient $\kappa > 1$, or $\kappa = 2$ with an additional factor on the right-hand side of (A9), over a bounded domain. It turns that this choice of $\phi$ is precisely the one that we will make for phase retrieval.

# Appendix B KL Functions

This section encompasses all the essential components required for the axiomatization of convergence for KL functions and therefore can be skipped by an experienced reader. We start by defining the non-smooth KL property which is an additional assumption on the class of functions that we consider. This property gives a hint about the geometric bearing of the function near the point where it is satisfied.

**Definition B.1. (Non-smooth KL property)** A proper and lower semicontinuous function $g: \mathbb{R}^n \to \bar{\mathbb{R}}$ has the KL property at a point $x_* \in \text{dom}(g)$ if there exists a neighborhood $U_{x_*}$, $\eta > 0$ and a concave real-valued function $\varphi \in C^1([0, \eta])$, with $\varphi(0) = 0$ and $\varphi' > 0$, such that

$$ \varphi'(g(x) - g(x_*)) \text{dist}(0, \partial g(x)) \ge 1, \quad \forall x \in U_{x_*} \cap \{x \in \mathbb{R}^n : g(x_*) < g(x) < g(x_*) + \eta\}. $$

If $g$ has the KL property at each point of $\text{dom}(g)$, $g$ is called a KL function.

$\varphi$ is known as the desingularizing function. The KL property is also closely related to error bounds and the broader notion of “(sub)metric regularity”. We refer the reader to [40] and [17] for a detailed study of these notions. In general, it is not obvious to check whether a given function is KL or not. Actually, this is a very deep question that has been studied at the interface of analysis and algebraic geometry. For smooth functions, it has been shown that semi-algebraic and sub-analytic are Łojaciewicz in the seminal works of [51, 52, 41]. This has been extended to the nonsmooth case and then widely studied in the scope of optimization in [16, 15, 17]. For instance, functions definable on $0$-minimal structures are KL. This covers most functions studied in practice, and for instance those in this manuscript.

The following uniformization of the KL property will be very useful; see [18, Lemma 6].

**Lemma B.2 (Uniformized KL property).** Let $\Omega$ be a compact set, and let $g: \mathbb{R}^n \to \bar{\mathbb{R}}$ be a proper and closed function. Assume that $g$ is constant on $\Omega$ and satisfies the KL property at
each point of $\Omega$. Then, there exist $\varepsilon > 0$, $\eta > 0$, if there exists $\eta > 0$ and a concave real-valued function $\varphi \in C^1([0, \eta])$ with $\varphi(0) = 0$, $\varphi' > 0$ such that for all $x_* \in \Omega$, one has

$$ \varphi'(g(x) - g(x_*))\mathrm{dist}(0, \partial g(x)) \ge 1, $$

and all $x \in \mathbb{R}^n$ such that $\mathrm{dist}(x, \Omega) < \varepsilon$ and $g(x_*) < g(x) < g(x_*) + \eta$.

For the convergence of our inertial Bregman proximal gradient (IBPG) algorithm, we will use a general convergence mechanism as first axiomatized in [3] for descent algorithms and generalized in [18] on the so-called PALM algorithm, so that it can be used and applied to any given algorithm such as ours (see also, e.g. [20, Appendix 6] for a self-contained presentation). The main goal is to prove that the whole sequence $(x_k)_{k \in \mathbb{N}}$ generated by IBPG, converges to a critical point. For that purpose, considering a Lyapunov function $\Psi$ associated to IBPG, it has to satisfy the following three key conditions.

**Definition B.3 (Descent-like method).** A sequence $(x_k)_{k \in \mathbb{N}}$ is called descent-like for the function $\Psi$ if the following conditions hold:

(C.1) *Sufficient decrease condition.* There exists a positive scalar $\rho_1$ such that

$$ \rho_1 \|x_k - x_{k-1}\|^2 \le \Psi(x_k, x_{k-1}) - \Psi(x_{k+1}, x_k), \quad \forall k \in \mathbb{N}. $$

(C.2) *Relative error condition.* There exists $K \in \mathbb{N}$ and $\rho_2 > 0$ such that $\forall k \ge K$, there exists $v_{k+1} \in \partial\Psi(x_{k+1}, x_k)$ such that

$$ \|v_{k+1}\| \le \rho_2 (\|x_{k+1} - x_k\| + \|x_k - x_{k-1}\|). $$

(C.3) *Continuity condition.* Let $x^*$ be a limit point of a subsequence $(x_k)_{k \in \mathcal{K} \subset \mathbb{N}}$ then we have that $\limsup_{k \in \mathcal{K} \subset \mathbb{N}} \Psi(x_k, x_{k-1}) \le \Psi(x^*) \stackrel{\mathrm{def}}{=} \Psi(x^*, x^*)$.

Condition (C.1) is intended to model a descent property of the Lyapunov function, and hence a dissipation of the energy $\Psi$. (C.2)<fnref n="3" /> originates from the well-known fact that most algorithms in optimization generate sequences via exact or inexact minimization of subproblems and condition (C.2) reflects relative inexact optimality conditions for such minimization subproblems. Condition (C.3) is a weak requirement which, in particular, holds when $\Psi$ is continuous. However, the latter is not mandatory in general as the nature of the algorithm (IBPG here) will force the sequences to comply with (C.3) under a simple lower semicontinuity assumption.

Equipped with Definition B.3, and when $\Psi$ satisfies the KL property, the following global convergence result holds true.

**Theorem B.4 (Global convergence).** Let $(x_k)_{k \in \mathbb{N}}$ be a bounded sequence generated by a descent-like method for $\Psi$. If $\Psi$ satisfies the KL property, then the sequence $(x_k)_{k \in \mathbb{N}}$ has finite length, i.e., $\sum_{k \in \mathbb{N}} \|x_{k+1} - x_k\| < +\infty$ and it converges to $x_* \in \mathrm{crit}(\Psi)$.

<footnote n="3">³The original version of this condition in [3] involves only the first term in the bound. The reasoning however remains the same with this version of the inequality; see e.g. [21, 47, 56].</footnote>
## Appendix C Proof of Global Convergence

**Lemma C.1.** Let $(x_k)_{k \in \mathbb{N}}$, $(y_k)_{k \in \mathbb{N}}$ and $(z_k)_{k \in \mathbb{N}}$ be generated by IBPG (Algorithm 1). Then the following holds true

(i) If $(x_k)_{k \in \mathbb{N}}$ is bounded, then the other sequences are also bounded.

(ii) If $(x_k)_{k \in \mathbb{N}}$ has a limit, then the remaining sequences also converge to the same limit.

*Proof* (i) $z_k$ being a convex combination of $x_k$ and $x_{k-1}$, the conclusion is immediate. The same holds also for $y_k$.

(ii) Suppose that $x_k \to x^*$. We have by definition

$$ \|z_k - x^*\| \le \|x_k - x^*\| + \|x_k - x_{k-1}\|. $$

Passing to the limit we get that $z_k \to x^*$. Moreover,

$$ \|y_k - x^*\| \le (1 - a_k) \|z_k - x^*\| + a_k \|x_k - x^*\| \le \|z_k - x^*\| + \|x_k - x^*\|. $$

and thus $y_k \to x^*$.

□

Our global convergence analysis will be based on a Lyapunov analysis with the energy function $\Psi_k$ on $\mathbb{R}^{3n}$ defined as

$$ \Psi_k(x_k, x_{k-1}) = \Phi(x_k) - \inf \Phi + a_{k-1}^{1-\kappa} LD_{\psi}(x_k, x_{k-1}), $$

where $\kappa \in [1, 2]$ is the TSE parameter of $\psi$. The subscript $k$ underscores the fact that $\Psi_k$ depends on $a_{k-1}$. Observe that $\Psi_k$ is non-negative. The first part of $\Psi_k$ corresponds to the potential energy of IBPG seen as a dissipative (discrete) dynamical system. The second term, which captures how the iterates remain close to each other, can be interpreted as a discrete Bregman version of the kinetic energy (involving the discrete velocity) of the system. The following lemma shows that $\Psi_k$ is indeed a Lyapunov function for IBPG.

**Lemma C.2.** Under Assumptions A and (B.2), there exists $\nu \in [0, 1]$ such that the sequences generated by Algorithm 1 satisfy $\forall k \ge 0$

$$ \Psi_{k+1}(x_{k+1}, x_k) \le \Psi_k(x_k, x_{k-1}) - \nu LD_{\psi}(x_k, x_{k-1}) - (\bar{a}^{1-\kappa} - 1) LD_{\psi}(x_{k+1}, y_k). \quad (C1) $$

*Proof* By L-smoothness of $F$ relative to $\psi$, we have

$$ \Phi(x_{k+1}) \le F(y_k) + \langle \nabla F(y_k), x_{k+1} - y_k \rangle + G(x_{k+1}) + LD_{\psi}(x_{k+1}, y_k). \quad (C2) $$

Observe that

$$ x_{k+1} \in \underset{x \in \mathbb{R}^n}{\operatorname{Argmin}} F(y_k) + \langle \nabla F(y_k), x - y_k \rangle + G(x) + \frac{1}{\gamma_k} D_{\psi}(x, y_k). \quad (C3) $$

Then convexity of $G$ gives that $\forall x \in \mathbb{R}^n$ (see e.g. [27, Lemma 3.2])

$$ F(y_k) + \langle \nabla F(y_k), x_{k+1} - y_k \rangle + G(x_{k+1}) + \frac{1}{\gamma_k} D_{\psi}(x_{k+1}, y_k) \le F(y_k) + \langle \nabla F(y_k), x - y_k \rangle + G(x) \\ + \frac{1}{\gamma_k} D_{\psi}(x, y_k) - \frac{1}{\gamma_k} D_{\psi}(x_{k+1}, x). \quad (C4) $$
Inserting this into (C2), we get

$$
\begin{align*}
\Phi(x_{k+1}) &\le F(y_k) + \langle \nabla F(y_k), x - y_k \rangle + G(x) + \frac{1}{\gamma_k} D_\psi(x, y_k) - \frac{1}{\gamma_k} D_\psi(x_{k+1}, x) - \left(\frac{1}{\gamma_k} - L\right) D_\psi(x_{k+1}, y_k) \\
&= \Phi(x) - D_F(x, y_k) + \frac{1}{\gamma_k} D_\psi(x, y_k) - \frac{1}{\gamma_k} D_\psi(x_{k+1}, x) - \left(\frac{1}{\gamma_k} - L\right) D_\psi(x_{k+1}, y_k) \\
&\le \Phi(x) + \left(\frac{1}{\gamma_k} + L\right) D_\psi(x, y_k) - \frac{1}{\gamma_k} D_\psi(x_{k+1}, x) - \left(\frac{1}{\gamma_k} - L\right) D_\psi(x_{k+1}, y_k),
\end{align*}
$$

where we used again L-smoothness of F relative to ψ. Applying this inequality at x = x_k, we obtain

$$
\Phi(x_{k+1}) + \frac{1}{\gamma_k} D_\psi(x_{k+1}, x_k) \le \Phi(x_k) + \left(\frac{1}{\gamma_k} + L\right) D_\psi(x_k, y_k) - \left(\frac{1}{\gamma_k} - L\right) D_\psi(x_{k+1}, y_k).
$$

We now use the TSP property twice to get

$$
\begin{align*}
D_{\psi}(x_k, y_k) &= D_{\psi}((1-a_k)x_k + a_k x_k, (1-a_k)z_k + a_k x_k) \le (1-a_k)^{\kappa} D_{\psi}(x_k, z_k) \quad \text{and} \\
D_{\psi}(x_k, z_k) &= D_{\psi}((1-a_{k-1})x_k + a_{k-1} x_k, (1-a_{k-1})x_{k-1} + a_{k-1} x_k) \le (1-a_{k-1})^{\kappa} D_{\psi}(x_k, x_{k-1}).
\end{align*}
$$

Combining the above inequalities, we arrive at

$$
\begin{align*}
\Phi(x_{k+1}) + \frac{1}{\gamma_k} D_\psi(x_{k+1}, x_k) \le \Phi(x_k) &+ \left(\frac{1}{\gamma_k} + L\right) (1-a_k)^\kappa (1-a_{k-1})^\kappa D_\psi(x_k, x_{k-1}) \\
& - \left(\frac{1}{\gamma_k} - L\right) D_\psi(x_{k+1}, y_k).
\end{align*}
$$

Thus, in view of the choice of the parameters, there exists $\nu[0, 1]$ such that

$$
\begin{align*}
\Phi(x_{k+1}) + a_k^{1-\kappa} LD_\psi(x_{k+1}, x_k) &\le \Phi(x_k) + a_{k-1}^{1-\kappa} LD_\psi(x_k, x_{k-1}) - \nu a_{k-1}^{1-\kappa} LD_\psi(x_k, x_{k-1}) \\
&\quad - (a_k^{1-\kappa} - 1) LD_\psi(x_{k+1}, y_k) \\
&\le \Phi(x_k) + a_{k-1}^{1-\kappa} LD_\psi(x_k, x_{k-1}) - \nu LD_\psi(x_k, x_{k-1}) \\
&\quad - (\bar{a}^{1-\kappa} - 1) LD_\psi(x_{k+1}, y_k).
\end{align*}
$$

Subtracting inf Φ on both sides, we get the claimed inequality.

<figure><img src="image_7.png" /></figure>

Capitalizing on the above descent property of the Lyapunov function Ψₖ, we get some
preliminary convergence properties<fnref n="4" />.

**Proposition C.3.** Under Assumptions A and (B.1)-(B.2), the sequences generated by Algorithm 1 (IBPG) are such that:

(i) The sequence $\{\Psi_k(x_k, x_{k-1})\}_{k \in \mathbb{N}}$ is decreasing and thus $\lim_{k \to +\infty} \Psi_k(x_k, x_{k-1}) \ge 0$ exists.

(ii) $\sum_{k \in \mathbb{N}} \|x_k - x_{k-1}\|^2 < \infty$ and thus $\lim_{k \to \infty} \|x_k - x_{k-1}\| = 0$.

(iii) We have the rate

$$
\min_{0 \le i \le k} \|x_i - x_{i-1}\|^2 \le \frac{2(\sigma_{\psi} \nu L)^{-1} \Psi_0(x_0, x_{-1})}{k+1}.
$$

<footnote n="4">⁴Only convexity of $\psi$ is needed for claim (i) to hold.</footnote>
*Proof* (i) The first claim comes from Lemma C.2 and non-negativity of $D_\psi$ thanks to convexity of $\psi$. The existence of the limit then follows as $\Psi_k$ is non-negative and decreasing (recall that $\Phi$ is bounded from below, $a_k \ge a > 0$ and $\kappa \in ]0, 1]$).

(ii) We use $\sigma_\psi$-strong convexity of $\psi$ and then sum (C1) dropping the last non-negative term to get

$$ \frac{\sigma_\psi \nu L}{2} \sum_{i=0}^{k} \|x_i - x_{i-1}\|^2 \le \nu L \sum_{i=0}^{k} D_\psi(x_i, x_{i-1}) \\ \le \Psi_0(x_0, x_{-1}) - \Psi_k(x_k, x_{k-1}) \le \Psi_0(x_0, x_{-1}). \quad (C5) $$

Passing to the limit as $k \to +\infty$ we get the summability claim.

(iii) We have

$$ (k+1) \min_{0 \le i \le k} \|x_i - x_{i-1}\|^2 \le \sum_{i=0}^{k} \|x_i - x_{i-1}\|^2. $$

Combining this with (C5), we conclude.

□

To prove global convergence, it is sufficient to show that IBPG is a descent-like method according to Definition B.3 and then to invoke Theorem B.4. For this we need to construct an appropriate function $\Psi$ that verifies the conditions of Definition B.3. The sequence of functions $\Psi_k$ could do the job if $a_k$ is taken fixed, say equal to $a \in [a, \bar{a}]$, for all $k \ge K$, where $K$ is arbitrarily large (see the discussion in Remark 3.3). We therefore consider the energy function

$$ \Psi(x_k, x_{k-1}) = \begin{cases} \Psi_k(x_k, x_{k-1}) & \text{if } k < K, \\ \Phi(x_k) - \inf \Phi + a^{1-\kappa} L D_\psi(x_k, x_{k-1}) & \text{otherwise.} \end{cases} $$

Observe first that $\Psi$ is KL since both $\Phi$ and $\psi$ are. The following proposition shows that the sequence generated by IBPG is a descent-like sequence for the new Lyapunov function $\Psi$.

**Proposition C.4.** Assume that Assumptions A and Assumptions B hold. Let $(x_k)_{k \in \mathbb{N}}$ be a bounded sequence generated by Algorithm 1. Then $(x_k)_{k \in \mathbb{N}}$ is a gradient-like descent sequence. Moreover the set of cluster points of $(x_k)_{k \in \mathbb{N}}$ is a nonempty compact set of $\text{crit}(\Phi)$.

*Proof*

* Sufficient decrease condition. From Lemma C.2 $\sigma_\psi$-strong convexity of $\psi$, we have $\forall k \in \mathbb{N}$

$$ \Psi(x_{k+1}, x_k) \le \Psi(x_k, x_{k-1}) - \nu L D_\psi(x_k, x_{k-1}) \\ \le \Psi(x_k, x_{k-1}) - \frac{\sigma_\psi \nu L}{2} \|x_k - x_{k-1}\|^2 $$

which shows (C.1) in Definition B.3.

* Relative error condition. We have to show (C.2) in Definition B.3.

$\psi$ is $C^2$ hence $D_\psi(\cdot, \cdot)$ is $C^1$ jointly in its arguments. The sum rule of the limiting subdifferential applies in this case and tells us that for $k \ge K$, we have

$$ \partial \Psi(x_{k+1}, x_k) = \left( \begin{array}{l} \nabla F(x_{k+1}) + \partial G(x_{k+1}) + a^{1-\kappa} L (\nabla \psi(x_{k+1}) - \nabla \psi(x_k)) \\ -a^{1-\kappa} L \nabla^2 \psi(x_k) (x_{k+1} - x_k) \end{array} \right). \quad (C6) $$

From the update equation of $x_{k+1}$ by IBPG, we have

$$ \nabla \psi(y_k) - \gamma \nabla F(y_k) - \nabla \psi(x_{k+1}) \in \gamma \partial G(x_{k+1}), \quad (C7) $$
where $\gamma = a^{\kappa-1}/L$. Set $v_{k+1} \stackrel{\text{def}}{=} (v_{k+1}^1, v_{k+1}^2)$ where

$$v_{k+1}^1 = (\nabla F(x_{k+1}) - \nabla F(y_k)) + a^{1-\kappa} L (\nabla \psi(y_k) - \nabla \psi(x_k))$$

$$\text{and } v_{k+1}^2 = -a^{1-\kappa} L \nabla^2 \psi(x_k)(x_{k+1} - x_k).$$

In view of (C6) and (C7), we have

$$v_{k+1} \in \partial\Psi(x_{k+1}, x_k).$$

We shall now bound $\|v_{k+1}\|^2 = \|v_{k+1}^1\|^2 + \|v_{k+1}^2\|^2$. Since $(x_k)_{k \in \mathbb{N}}$ is bounded and $\nabla \psi$ is Lipschitz continuous on any bounded subset by (B.3), there exists $L_\psi$ such that

$$\|v_{k+1}^2\| \le a^{1-\kappa} LL_\psi \|x_{k+1} - x_k\| \le \underline{a}^{1-\kappa} LL_\psi \|x_{k+1} - x_k\|.$$

Moreover, $(x_k)_{k \in \mathbb{N}}$ is also bounded by Lemma C.1. Assumption (B.3) then entails that there exist $L_F > 0$ such that

$$\|v_{k+1}^1\| \le L_F \|x_{k+1} - x_k\| + L_F \|x_k - y_k\| + \underline{a}^{1-\kappa} LL_\psi \|x_k - y_k\|.$$

where we used that $a \in [0, 1]$ and $\kappa \in [1, 2]$. Now, by definition of the iterates

$$x_k - y_k = x_k - z_k - a(x_k - z_k) = (1-a)(x_k - x_{k-1} - a(x_k - x_{k-1})) = (1-a)^2(x_k - x_{k-1}).$$

Therefore

$$\begin{aligned} \|v_{k+1}^1\| &\le L_F \|x_{k+1} - x_k\| + (L_F + \underline{a}^{1-\kappa} LL_\psi)(1-a)^2 \|x_k - x_{k-1}\| \\ &\le L_F \|x_{k+1} - x_k\| + (L_F + \underline{a}^{1-\kappa} LL_\psi)(1-\underline{a})^2 \|x_k - x_{k-1}\|. \end{aligned}$$

Taking $\rho_2 = L_F + (L_F + \underline{a}^{1-\kappa} LL_\psi)(1-\underline{a})^2 + \underline{a}^{1-\kappa} LL_\psi$, we get the claim.

* Continuity condition. Since $(x_k)_{k \in \mathbb{N}}$ is bounded, its set of cluster points is a nonempty set. It is also compact set as the intersection of compact sets. Let us consider a subsequence $(x_{k_j})_{j \in \mathbb{N}}$ that converges to some limit $x^*$. From Proposition C.3(i) and Lemma C.1, we have that $(x_{k_j+1})_{k \in \mathbb{N}}$ and $(y_{k_j})_{k \in \mathbb{N}}$ converge to the same limit. Now arguing as in (C3)-(C4), we get

$$\begin{aligned} G(x_{k_j+1}) &\le G(x^*) + \langle \nabla F(y_{k_j}), x^* - x_{k_j+1} \rangle - \frac{1}{\gamma} D_\psi(x_{k_j+1}, y_{k_j}) + \frac{1}{\gamma} D_\psi(x^*, y_{k_j}) - \frac{1}{\gamma} D_\psi(x_{k_j+1}, x^*) \\ &\le G(x^*) + \langle \nabla F(y_{k_j}), x^* - x_{k_j+1} \rangle + \frac{1}{\gamma} D_\psi(x^*, y_{k_j}). \end{aligned}$$

Passing to the limit as $j \to +\infty$ and using continuity of $D_\psi$ we get that

$$\limsup_{j \to +\infty} G(x_{k_j+1}) \le G(x^*).$$

Combining this with continuity of $F$ proves (C.3) in Definition B.3.

Let $x_{k_j} \to x^*$. We argued above that $y_{k_j} \to x^*$ and $x_{k_j+1} \to x^*$. We then have from (C7), continuity of $\nabla \psi$ and $\nabla F$, and sequential closedness of $\nabla G$ that

$$-\nabla F(x^*) \in \partial G(x^*),$$

i.e. $x^* \in \text{crit}(\Phi)$.

□
Proof of Theorem 3.2

Proof (i) This claim comes from Proposition C.3.

(ii) This is a consequence of Proposition C.4 and Theorem B.4 after observing from (C6) that $\mathrm{crit}(\Psi) = \{(x_*, x_*): x_* \in \mathrm{crit}(\Phi)\}$.

(iii) We now turn to proving convergence to a global minimizer. We introduce the following extended variable, $\tilde{x}_k = (x_k, x_{k-1}) \in \mathbb{R}^n \times \mathbb{R}^n$ such that $\Psi(x_k, x_{k-1}) = \Psi(\tilde{x}_k)$ for all $k \in \mathbb{N}$. Let us choose radius $r > \rho_2 > 0$ such that $\eta < \rho_1(r - \rho_2)^2$. Let us suppose that the initial point $x_0$ is chosen such that the following conditions hold,

$$ \Phi(x^*) = \Psi(\tilde{x}^*) \le \Psi(\tilde{x}_0) < \Psi(\tilde{x}^*) + \eta = \Phi(x^*) + \eta \quad (C8) $$

$$ \|x_0 - x^*\| + 2\sqrt{\frac{\Psi(\tilde{x}_0) - \Psi(\tilde{x}^*)}{\rho_1}} + \frac{\rho_2}{\rho_1}\varphi(\Psi(\tilde{x}_0) - \Psi(\tilde{x}^*)) < \rho. \quad (C9) $$

The condition (C.1) combined with (C8) imply that for any $k \in \mathbb{N}$, $\Psi(\tilde{x}^*) \le \Psi(\tilde{x}_{k+1}) \le \Psi(\tilde{x}_0) < \Psi(\tilde{x}^*) + \eta$, and moreover

$$ \|x_{k+1} - x_k\| \le \sqrt{\frac{\Psi(\tilde{x}_{k+1}) - \Psi(\tilde{x}_{k+2})}{\rho_1}} \le \sqrt{\frac{\Psi(\tilde{x}_{k+1}) - \Psi(\tilde{x}^*)}{\rho_1}}. \quad (C10) $$

We deduce that if for any $k \in \mathbb{N}$, $x_k \in B(x^*, \rho)$ then $x_{k+1} \in B(x^*, r)$. Indeed, by the triangle inequality

$$ \|x_{k+1} - x^*\| \le \|x_k - x^*\| + \sqrt{\frac{\Psi(\tilde{x}_{k+1}) - \Psi(\tilde{x}^*)}{\rho_1}} = \rho + (r - \rho) = r. \quad (C11) $$

It remains to show that $\forall k \in \mathbb{N}$, $x_k \in B(x^*, \rho)$. We argue by induction. The triangle inequality gives

$$ \|x_1 - x^*\| \le \|x_0 - x^*\| + \sqrt{\frac{\Psi(\tilde{x}_1) - \Psi(\tilde{x}^*)}{\rho_1}} \le \|x_0 - x^*\| + \sqrt{\frac{\Psi(\tilde{x}_0) - \Psi(\tilde{x}^*)}{\rho_1}} < \rho, $$

which means that $x_1 \in B(x^*, \rho)$. We also have

$$ \|x_{k+1} - x^*\| \le \|z_0 - x^*\| + 2 \|x_1 - x_0\| + \sum_{j=1}^{k} \|x_k - x_{k-1}\|, $$

Standard arguments with the KL inequality show that

$$ \sum_{i=l+1}^{k} \|x_k - x_{k-1}\| \le \frac{\rho_2}{\rho_1} \varphi(\Psi(x_{l+1}, x_l) - \Psi(\tilde{x}^*)). \quad (C12) $$

Applying this bound (C12) with $l = 0$ and combining with (C9) yields

$$ \|x_{k+1} - x^*\| \le \|x_0 - x^*\| + 2\sqrt{\frac{\Psi(\tilde{x}_0) - \Psi(\tilde{x}^*)}{\rho_1}} + \frac{\rho_2}{\rho_1}\varphi(\Psi(\tilde{x}_0) - \Psi(\tilde{x}^*)) < \rho, $$

which implies that $x_{k+1} \in B(x^*, \rho)$.

If we start close enough to $x^*$ so that (C8)-(C9) holds the sequence $(x_k)_{k \in \mathbb{N}}$ will remain in the neighborhood $B(x^*, \rho)$ and converges to a critical point, say $x^*$. Moreover $\Psi(x_k) \to \Phi(x^*) \ge \Phi(x^*)$. Let us assume that $\Phi(x^*) > \Phi(x^*)$. Since $\Phi$ has the KL property at $x^*$ and thus

$$ \varphi'(\Phi(x^*) - \Phi(x^*)) \mathrm{dist}(0, \partial\Phi(x^*)) \ge 1. $$

This is a contradiction since $\varphi'(s) > 0$ for $s \in ]0, \eta[$ and $\mathrm{dist}(0, \partial\Phi(x^*)) = 0$ since $x^*$ is a critical point. We conclude that $x^*$ is indeed a global minimizer. $\square$
# Appendix D Proofs of Local Convergence

At this juncture, we pause to present the following Lemma, which forms the essence of our framework. It is an extension of [11, Proposition 4.4] to the Bregman setting when the Bregman entropy is strongly convex.

**Lemma D.1. (A Firmly nonexpansive map)** Let $G \in \Gamma_0(\mathbb{R}^n)$ and $\psi$ be a strongly convex function with full domain then the following map $J \stackrel{\text{def}}{=} (\nabla\psi + \partial G)^{-1}$ and $\text{Id} - J$ are firmly non-expansive.

**Remark D.2.** Let us notice that if $\psi$ is not strongly convex but just convex, then from (D1) one has that $J$ is only monotone.

*Proof* Let $x, y \in \mathbb{R}^n$ such that $p = J(x)$ and $q = J(y)$. Since $G \in \Gamma_0(\mathbb{R}^n)$ we have the following statements

$$D_G^{-\nabla\psi(p)+x}(q,p) \geq 0 \text{ and } D_G^{-\nabla\psi(q)+y}(p,q) \geq 0.$$

If we sum up both inequalities, we find that

$$D_G^{-\nabla\psi(p)+x}(q,p) + D_G^{-\nabla\psi(q)+y}(p,q) \geq 0 \iff \langle Jx - Jy, x - y \rangle \geq D_\psi(q,p) + D_\psi(p,q). \quad (\text{D1})$$

Therefore we get the desired result through strong convexity of $\psi$ and [11, Proposition 4.4]. □

## D.1 Proof of Lemma 4.3

*Proof* We have $0 \in \partial\Phi(x_*)$, since $G \in \text{PSF}_{x_*}(M_{x_*})$ and $F \in C^1$ thanks to [44, Corollay 4.7] (smooth perturbation of partly smooth functions), we have that $\Phi \in \text{PSF}_{x_*}(M_{x_*})$. From the global convergence Theorem 3.2, we have $x_k \xrightarrow{\Phi} x_*$. Let us consider the iteration of Algorithm 1 and define

$$v_{k+1} = \frac{1}{\gamma_k} (-\nabla\psi(x_{k+1}) + \nabla\psi(y_k)) + \nabla F(x_{k+1}) - \nabla F(y_k),$$

we have that $\forall k \in \mathbb{N}, v_k \in \partial\Phi(x_k)$. Besides,

$$\|v_{k+1}\| \leq \frac{1}{\gamma_k} \|\nabla\psi(x_{k+1}) - \nabla\psi(y_k)\| + \|\nabla F(x_{k+1}) - \nabla F(y_k)\|$$

Since $(x_k)_{k \in \mathbb{N}}$ is assumed to be bounded and so thus $(y_k)_{k \in \mathbb{N}}$ by Lemma C.1, we deduce from Assumption (B.3) that there exists a positive scalar $M_1, M_2 > 0$ such that

$$\begin{align*}
\|v_{k+1}\| &\leq \left(\frac{M_1}{\gamma_k} + M_2\right) \|x_{k+1} - y_k\|, \\
&= \left(\frac{M_1}{\gamma_k} + M_2\right) \|x_{k+1} - (1 - a_k)(1 - a_{k-1})x_{k-1} - (1 - a_k)a_{k-1}x_k - a_kx_k\|, \\
&\leq \left(\frac{M_1}{\gamma_k} + M_2\right) (\|x_{k+1} - x_k\| + (1 - a_k)(1 - a_{k-1}) \|x_k - x_{k-1}\|), \\
&\leq \left(\frac{M_1}{\gamma_k} + M_2\right) \|x_{k+1} - x_k\| + \left(\frac{M_1}{\gamma_k} + M_2\right) (1 - a_k)(1 - a_{k-1}) \|x_k - x_{k-1}\|,
\end{align*}$$

Therefore, we obtained that for $k \to \infty$ we get that $a_k \to a \in [a, \bar{a}]$ implying that $(\gamma_k)_{k \in \mathbb{N}}$ converges. The latter combined with the fact that $x_k \to x_*$, we deduce that $\|v_k\| \to 0$ as $k \to \infty$. We conclude that $x_k$ identifies $M_{x_*}$ thanks to [31, Proposition 10.12].
(i) If the active manifold $\mathcal{M}_{x^*}$ is an affine subspace, then $\mathcal{M}_{x^*} = x^* + T_{x^*}$ due to the normal sharpness property, and the claim follows immediately.

(ii) When $G$ is locally polyhedral around $x^*$, $\mathcal{M}_{x^*}$ is an affine subspace, and the identification of $(y_k)_{k \in \mathbb{N}}$ and $(z_k)_{k \in \mathbb{N}}$ follows from (i). For the rest, it is sufficient to observe that, by polyhedrality, for any $x$ in $\mathcal{M}_{x^*}$ near $x^*$, $\partial G(x) = \partial G(x^*)$, combining Fact F.4 and Fact F.5, we arrive at the second conclusion.
□

The next Lemma gives the spectral properties of the matrices defined (4.1).

**Lemma D.3.** Under the Assumption A-B, let $x_* \in \text{RI}(\Phi) \cap \text{ND}(\Phi)$ such that $F$ is locally $C^2$ around $x_*$. Then for any stepsize $\gamma \in ]0, 1/L]$, we have

(i) $H_F$ is symmetric positive definite with eigenvalues in $]\gamma\sigma\sigma_\psi, L\lambda_{\max}(\nabla^2\psi(x_*))\gamma]$.

(ii) $V$ has eigenvalues in

$$[\lambda_{\max}(\nabla^2\psi(x_*))(1 - \gamma L), \lambda_{\max}(\nabla^2\psi(x_*)) - \gamma\sigma_\psi\sigma, ]$$

hence $H_{\psi}^{-1}V$ has eigenvalues in

$$[1 - \gamma L \frac{\lambda_{\max}(\nabla^2\psi(x_*))}{\sigma_\psi}, 1 - \gamma\sigma \frac{\sigma_\psi}{\lambda_{\max}(\nabla^2\psi(x_*))}, ]$$

(iii) $WH_\psi$ has eigenvalues in $]0, q_\psi(x_*)]$.

(iv) If either $x_* \in \text{ND}(\Psi)$ or $\mathcal{M}_{x_*}$ is affine then $WV$ has eigenvalues in $]-\Lambda, \Lambda]$ where we recall that $\Lambda = |q_\psi(x_*) - \gamma\sigma|$ with $q_\psi(x_*) = \frac{\lambda_{\max}(\nabla^2\psi(x_*))}{\sigma_\psi}$.

### Proof

(i) Combining the fact that $F$ is locally $C^2$ and $x_* \in \text{RI}(\Psi)$ we get that $\exists\sigma > 0, \forall h \in T_{x_*}$,

$$\langle h; \nabla^2 F(x_*) h \rangle \geq \langle h; \sigma \nabla^2 \psi(x_*) h \rangle \geq \sigma \sigma_\psi \|h\|^2$$

where the last part comes from the strong convexity of $\psi$. This implies that $\lambda_{\min}(H_F) \geq \gamma\sigma\sigma_\psi$. Since $F$ is $L$-smooth relative to $\psi$, $\forall h \in T_{x_*}$,

$$\langle h; \nabla^2 F(x_*) h \rangle \leq L \langle h; \nabla^2 \psi(x_*) h \rangle \leq L \lambda_{\max}(\nabla^2 \psi(x_*)) \|h\|^2.$$

(ii) We have that $V = H_\psi - H_F$ and that $H_\psi$ has eigenvalues in $[\sigma_\psi, \lambda_{\max}(\nabla^2\psi(x_*))]$. We combine this with the previous claim on $H_F$ and the eigenvalues of the difference of two positive definite matrices.

For the second claim, let us observe that $H_{\psi}^{-1}V$ is similar to the following matrix $\text{Id} - H_{\psi}^{-1/2}H_F H_{\psi}^{-1/2}$. Besides, we have that $H_{\psi}^{-1/2}H_F H_{\psi}^{-1/2}$ has eigenvalues in $[\gamma\sigma \frac{\sigma_\psi}{\lambda_{\max}(\nabla^2\psi(x_*))}, \gamma L \frac{\lambda_{\max}(\nabla^2\psi(x_*))}{\sigma_\psi}]$ thus, the difference of two positive definite matrices.

(iii) From Remark 4.4, we have that $W$ has eigenvalues in $]0, 1/\sigma_\psi]$ which implies that $WH_\psi$ has eigenvalues in $]0, \frac{\lambda_{\max}(\nabla^2\psi(x_*))}{\sigma_\psi}]$.
(iv) Let us observe that $WV = W^{1/2}(W^{1/2}VW^{1/2})W^{-1/2}$, thus $WV$ is similar to $W^{1/2}VW^{1/2}$. We have

$$ \begin{aligned} \|W^{1/2}VW^{1/2}\| & \le \|W^{1/2}\| \|V\| \|W^{1/2}\| \le |\lambda_{\max}(\nabla^2\psi(x_\star)) - \gamma\sigma_\psi\sigma| \|W^{1/2}\|^2. \\ & = \left| \frac{\lambda_{\max}(\nabla^2\psi(x_\star))}{\sigma_\psi} - \gamma\sigma \right|, \end{aligned} $$

where we used the fact that either $x_\star \in \text{ND}(\Psi)$ or $M_{x_\star}$ is affine holds then from Remark 4.4 we have $\|W\| \le \frac{1}{\sigma_\psi}$.

□

Let us define the following matrices i.e.,

$$ H_F^k \stackrel{\text{def}}{=} \gamma_k P_{T_{x_*}} \nabla^2 F(x_*) P_{T_{x_*}}, \quad V^k \stackrel{\text{def}}{=} H_\psi - H_F^k, \quad U^k \stackrel{\text{def}}{=} \gamma_k \nabla^2_{M_{x_*}} \Phi(x_*) P_{T_{x_*}} - H_F. \quad (D2) $$

To enhance readability, we introduce simplified notation for any $k \in \mathbb{N}$,

$$ b_k = (1 - a_k)a_{k-1} + a_k \quad \text{and} \quad c_k = (1 - a_k)(1 - a_{k-1}), \quad (D3) $$

$$ b = 2a - a^2 \quad \text{and} \quad c = (1 - a)^2. \quad (D4) $$

From this notation, we have this obvious Lemma.

**Lemma D.4.** Let us consider the sequences define in (D3), If $a_k \to a$ as $k \to \infty$ then we have $b_k \to b$ and $c_k \to c$.

Following the work of [48, Section B], we also define the matrices

$$ M_1^k = [bW(V^k - V), cW(V^k - V)], \quad M_2^k = [(b_k - b)WV^k, (c_k - c)WV^k]. $$

Therefore, we have the following proposition.

**Proposition D.5.** Under the same assumptions as Proposition 4.6, for $k$ large enough we have

$$ \|y_k - x_*\| = O(\|d_k\|), \|r_{k+1}\| = O(\|d_k\|), \|x_{k+1} - y_k\| = O(\|d_k\|), \quad (D5) $$

$$ \|\nabla F(y_k) - F(x_{k+1})\| = O(\|d_k\|), \quad \|\nabla\psi(y_k) - \nabla\psi(x_{k+1})\| = O(\|d_k\|), \quad (D6) $$

and

$$ \|W(U^k - U)(x_{k+1} - x_*)\| = o(\|d_k\|), \quad \|M_1^k d_k\| = o(\|d_k\|), \quad \|M_2^k d_k\| = o(\|d_k\|). \quad (D7) $$

*Proof*

* From the definition of the sequences we have,

$$ \begin{aligned} \|y_k - x_*\| &= \|(1 - a_k)z_k + a_k x_k - x_*\|, \\ &= \|(1 - a_k)(zk - x_*) + a_k(x_k - x_*)\|. \end{aligned} $$
$$
\begin{align*}
&= \|(1 - a_k)(1 - a_{k-1})r_{k-1} + ((1 - a_k)a_{k-1} + a_k)r_k\|, \\
&\le (1 - a_k)(1 - a_{k-1}) (\|r_{k-1}\| + \|r_k\|), \\
&\le \sqrt{2} \|d_k\|.
\end{align*}
$$

• We recall that $r_{k+1} = x_{k+1} - x_{\star}$ thus,

$$
\begin{align*}
\|r_{k+1}\| &= \left\| (\nabla\psi + \gamma_k \partial G)^{-1} (\nabla\psi(y_k) - \gamma_k \nabla F(y_k)) - (\nabla\psi + \gamma_k \partial G)^{-1} (\nabla\psi(x_\star) - \gamma_k \nabla F(x_\star)) \right\|, \\
&\le \|\nabla\psi(y_k) - \nabla\psi(x_\star)\| + \gamma_k \|\nabla F(y_k) - \nabla F(x_\star)\|, \\
&\le (M_1 + M_2 \gamma_k) \|y_k - x_\star\|, \\
&\le (M_1 + M_2 \gamma_k) \sqrt{2} \|d_k\|,
\end{align*}
$$

where we used the non-expansiveness of the mapping $(\nabla\psi + \gamma_k \partial G)^{-1}$ (see Lemma D.1), the boundedness of the sequence and again Assumption (A.3).

• From Assumption (A.3), there exists $M_1 > 0$ large enough such that

$$
\|\nabla\psi(y_k) - \nabla\psi(x_{k+1})\| \le M_1 \|x_{k+1} - y_k\| = O(\|d_k\|).
$$

Similarly, there exists $M_2$ large enough such that

$$
\|\nabla F(y_k) - \nabla F(x_{k+1})\| \le M_2 \|y_k - x_{k+1}\| \le M_2 \sqrt{2} \|d_k\| = O(\|d_k\|).
$$

• Let us now turn to the proof of (D7), from the definition of $\Phi$, we have that

$$
\lim_{k \to \infty} \frac{\|W(U^k - U)r_{k+1}\|}{\|r_{k+1}\|} \le \lim_{k \to \infty} |\gamma_k - \gamma| \|W\| \|\nabla_{M_{x_*}} \Phi(x_*) P_{T_{x_*}}\| = 0, \quad (D8)
$$

since $\gamma_k \to \gamma$ which means that $\|W(U^k - U)r_{k+1}\| = o(\|r_{k+1}\|) = o(\|d_k\|)$, where we have used (D5).

• We have that,

$$
\begin{align*}
\lim_{k \to \infty} \frac{\|M_1^k d_k\|}{\|d_k\|} &= \lim_{k \to \infty} \frac{\|bW(V^k - V)r_k + cW(V^k - V)r_{k-1}\|}{\|d_k\|} \\
&\le \lim_{k \to \infty} \frac{\max(|b|, |c|) \|W\| \|V^k - V\| (\|r_k\| + \|r_{k+1}\|)}{\|d_k\|} \\
&\le \lim_{k \to \infty} \frac{\max(|b|, |c|) \|W\| |\gamma_k - \gamma| \|P_{T_{x_*}} \nabla^2 F(x_*) P_{T_{x_*}}\| \sqrt{2} \|d_k\|}{\|d_k\|} \\
&= \lim_{k \to \infty} \sqrt{2} |\gamma_k - \gamma| \max(|b|, |c|) \|W\| \|P_{T_{x_*}} \nabla^2 F(x_*) P_{T_{x_*}}\| = 0,
\end{align*}
$$

as the term $\max(|b|, |c|) \|W\| \|P_{T_{x_*}} \nabla^2 F(x_*) P_{T_{x_*}}\|$ is bounded since
$\|W\| \|P_{T_{x_*}} \nabla^2 F(x_*) P_{T_{x_*}}\| \le L \lambda_{\max}(\nabla^2 \psi(x_*))^{-1} \|H_\psi\| \le 1$ and $\max(|b|, |c|) \le 1$.

• To finish the proof, we have

$$
\begin{align*}
\lim_{k \to \infty} \frac{\|M_2^k d_k\|}{\|d_k\|} &= \lim_{k \to \infty} \frac{\|(b_k - b)WV^k r_k + (c_k - c)WV^k r_{k-1}\|}{\|d_k\|}, \\
&\le \lim_{k \to \infty} \frac{\max(|b_k - b|, |c_k - c|) \|W\| \|V^k\| \sqrt{2} \|d_k\|}{\|d_k\|}, \\
&\le \lim_{k \to \infty} \max(|b_k - b|, |c_k - c|) \|W\| \|V^k\| \sqrt{2} = 0,
\end{align*}
$$

where we use (D5), the fact that $\|V^k\| \to \|V\| < \infty$, the boundedness of $\|W\|$ combined with Lemma D.4 to get the result.
□
## D.2 Proof of Proposition 4.6

*Proof*

Since $(z_k)_{k \in \mathbb{N}}$ is a sequence generated by Algorithm 1 converging to $x_* \in \mathrm{crit}(\Phi)$. From the finite identification Lemma 4.3, there exists $K \in \mathbb{N}$ such that $x_k$ is close enough to $x_*$ for $k \ge K$. Let $T_{x_{k+1}}, T_{x_*}$ be their corresponding tangent spaces, and define $\tau_{k+1}: T_{x_*} \to T_{x_{k+1}}$ the parallel translation along the unique geodesic joining $x_{k+1}$ to $x_*$. From the definition of the point $x_{k+1}$ and the fact that $x_* \in \mathrm{crit}(\Phi)$ we have

$$ 
\begin{aligned}
\nabla\psi(y_k) - \nabla\psi(x_{k+1}) - \gamma_k (\nabla F(y_k) - \nabla F(x_{k+1})) &\in \gamma_k \partial\Phi(x_{k+1}) \\
0 &\in \gamma_k \partial\Phi(x_*).
\end{aligned}
$$

We project now this inclusions over the tangents spaces $T_{x_{k+1}}$ and $T_{x_*}$ respectively to get that

$$ 
\begin{aligned}
\tau_{k+1}^{-1} P_{T_{x_{k+1}}} (\nabla\psi(y_k) - \nabla\psi(x_{k+1}) - \gamma_k (\nabla F(y_k) - \nabla F(x_{k+1}))) &= \gamma_k \tau_{k+1}^{-1} \nabla_{M_{x_*}} \Phi(x_{k+1}) \\
0 &= \gamma_k \nabla_{M_{x_*}} \Phi(x_*)
\end{aligned}
$$

where we have used the Fact F.5. After summing both lines and subtracting the value $\tau_{k+1}^{-1} P_{T_{x_{k+1}}} \nabla\psi(x_*)$ we have

$$ 
\begin{aligned}
& \gamma_k (\tau_{k+1}^{-1} \nabla_{M_{x_*}} \Phi(x_{k+1}) - \nabla_{M_{x_*}} \Phi(x_*)) \\
&= \tau_{k+1}^{-1} P_{T_{x_{k+1}}} (\nabla\psi(y_k) - \nabla\psi(x_{k+1})) - \gamma_k \tau_{k+1}^{-1} P_{T_{x_{k+1}}} (\nabla F(y_k) - \nabla F(x_{k+1})).
\end{aligned}
\quad (D9) $$

We combined Lemma F.1, (D6) and Lemma F.2 due to the local $C^2$-smoothness of $\psi$ to obtain that

$$ 
\begin{aligned}
\tau_{k+1}^{-1} P_{T_{x_{k+1}}} (\nabla\psi(y_k) - \nabla\psi(x_{k+1}))
&= P_{T_{x_*}} (\nabla\psi(y_k) - \nabla\psi(x_{k+1})) + o(\|x_{k+1} - y_k\|), \\
&= P_{T_{x_*}} (\nabla\psi(y_k) - \nabla\psi(x_*)) - P_{T_{x_*}} (\nabla\psi(x_{k+1}) - \nabla\psi(x_*)) + o(\|d_k\|), \\
&= P_{T_{x_*}} \nabla^2\psi(x_*) (y_k - x_*) + o(\|y_k - x_*\|) - P_{T_{x_*}} \nabla^2\psi(x_*) (x_{k+1} - x_*) + o(\|r_{k+1}\|) + o(\|d_k\|), \\
&= P_{T_{x_*}} \nabla^2\psi(x_*) (y_k - x_*) - P_{T_{x_*}} \nabla^2\psi(x_*) (x_{k+1} - x_*) + o(\|d_k\|), \\
&= P_{T_{x_*}} \nabla^2\psi(x_*) P_{T_{x_*}} (y_k - x_*) - P_{T_{x_*}} \nabla^2\psi(x_*) P_{T_{x_*}} (x_{k+1} - x_*) + o(\|d_k\|),
\end{aligned}
\quad (D10) $$

where we have used (D6), [50, Lemma 5.1] and the fact that $r_{k+1} = O(\|d_k\|)$. Using similar arguments, we have

$$ 
\begin{aligned}
\tau_{k+1}^{-1} P_{T_{x_{k+1}}} (\nabla F(y_k) - \nabla F(x_{k+1}))
&= P_{T_{x_*}} (\nabla F(y_k) - \nabla F(x_{k+1})) + o(\|d_k\|), \\
&= P_{T_{x_*}} (\nabla F(y_k) - \nabla F(x_*)) - P_{T_{x_*}} (\nabla F(x_{k+1}) - \nabla\psi(x_*)) + o(\|d_k\|), \\
&= P_{T_{x_*}} \nabla^2 F(x_*) (y_k - x_*) + o(\|y_k - x_*\|) - P_{T_{x_*}} \nabla^2 F(x_*) (x_{k+1} - x_*) + o(\|x_{k+1} - x_*\|) + o(\|d_k\|), \\
&= P_{T_{x_*}} \nabla^2 F(x_*) (y_k - x_*) - P_{T_{x_*}} \nabla^2 F(x_*) (x_{k+1} - x_*) + o(\|d_k\|), \\
&= P_{T_{x_*}} \nabla^2 F(x_*) P_{T_{x_*}} (y_k - x_*) - P_{T_{x_*}} \nabla^2 F(x_*) P_{T_{x_*}} (x_{k+1} - x_*) + o(\|d_k\|).
\end{aligned}
\quad (D11) $$

Moreover, we have

$$ \tau_{k+1}^{-1} \nabla_{M_{x_*}} \Phi(x_{k+1}) - \nabla_{M_{x_*}} \Phi(x_*) = \nabla_{M_{x_*}}^2 \Phi(x_*) P_{T_{x_*}} (x_{k+1} - x_*) + o(\|d_k\|). \quad (D12) $$

We replace now the expressions (D10), (D11), (D12), in (D9), we obtain

$$ 
\begin{aligned}
& (P_{T_{x_*}} \nabla^2 \psi(x_*) P_{T_{x_*}} + \gamma_k \nabla_{M_{x_*}}^2 \Phi(x_*) P_{T_{x_*}} - \gamma_k P_{T_{x_*}} \nabla^2 F(x_*) P_{T_{x_*}}) (x_{k+1} - x_*) = (H_\psi + U^k)(x_{k+1} - x_*) \\
& = P_{T_{x_*}} \nabla^2 \psi(x_*) P_{T_{x_*}} (y_k - x_*) - \gamma_k P_{T_{x_*}} \nabla^2 F(x_*) P_{T_{x_*}} (y_k - x_*) + o(\|d_k\|).
\end{aligned}
\quad (D13) $$
We get by factorizing and replacing (D2) that

$$ (H_{\psi} + U^k)(x_{k+1} - x_{\star}) = V^k(y_k - x_{\star}) + o(\|d_k\|). \qquad (D14) $$

Now we replace the expressions of the inertial term in term of $x_k, x_{k-1}$ to obtain that

$$ (H_{\psi} + U^k)r_{k+1} = ((1 - a_k)a_{k-1} + a_k) V^k r_k + (1 - a_k)(1 - a_{k-1})V^k r_{k-1} + o(\|d_k\|), \quad (D15) $$

We can write further

$$ (H_{\psi} + U) r_{k+1} = (U - U^k) r_{k+1} + ((1 - a_k)a_{k-1} + a_k) V^k r_k + (1 - a_k)(1 - a_{k-1})V^k r_{k-1} + o(\|d_k\|), $$

Thanks to Remark 4.4, it is possible to invert the matrix $U + H_{\psi}$ we have,

$$ r_{k+1} = W (U - U^k) r_{k+1} + ((1 - a_k)a_{k-1} + a_k) W V^k r_k + (1 - a_k)(1 - a_{k-1}) W V^k r_{k-1} + o(\|Wd_k\|), $$

Let us use the notation (D3) with the estimates (D7)

$$ d_{k+1} = \left( M + \begin{bmatrix} M_1^k \\ 0 \end{bmatrix} + \begin{bmatrix} M_2^k \\ 0 \end{bmatrix} \right) d_k + o(\|d_k\|) = M d_k + o(\|d_k\|), $$

which concludes the proof. □

## D.3 Proof of Proposition 4.8

Proof We have that

$$ M \begin{pmatrix} r_1 \\ r_2 \end{pmatrix} = \begin{bmatrix} (2a - a^2)H_{\psi}^{-1}V & (a-1)^2H_{\psi}^{-1}V \\ \text{Id} & 0 \end{bmatrix} \begin{pmatrix} r_1 \\ r_2 \end{pmatrix} = \begin{pmatrix} (2a - a^2)H_{\psi}^{-1}Vr_1 + (a-1)^2H_{\psi}^{-1}Vr_2 \\ r_1 \end{pmatrix} = \varrho \begin{pmatrix} r_1 \\ r_2 \end{pmatrix}, $$

therefore $r_1 = \varrho r_2$, and insert it in the first identity to get that

$$ \varrho^2 r_2 = (2a - a^2)H_{\psi}^{-1}V\varrho r_2 + (1 - a)^2 H_{\psi}^{-1}Vr_2 $$

which means that

$$ ((2a - a^2)\varrho + (1 - a)^2) H_{\psi} V r_2 = \varrho^2 r_2 $$

thus there exists $\eta$ such that $H_{\psi}^{-1}Vr_2 = \eta r_2$ moreover $\eta$ satisfies the following equation

$$ \varrho^2 - (2a - a^2)\varrho\eta - (1 - a)^2\eta = 0. \qquad (D16) $$

The rest of the proof follows exactly the same step as the proof of [48, Proposition 4.7, Corollary 4.9] and we conclude with Lemma D.3-(iv). □

# Appendix E Proof of the Escape Property

## E.1 Proof of Theorem 5.1

Let us define the following mapping for $a \in [\underline{a}, \bar{a}]$.

$$ \mathbf{T}(x_2, x_1) \stackrel{\text{def}}{=} \left[ \nabla \psi^{-1} \left( \nabla \psi \left( y(x_2, x_1) \right) - \gamma \nabla F(y(x_2, x_1)) \right) \right], \qquad (E1) $$
where

$$y(x_2, x_1) = (2a - a^2)x_2 + (1 - a)^2 x_1.$$

It is simple to see that $x_* \in \mathrm{crit}(\Phi)$ if and only if $(x_*, x_*)$ is a fixed point of the operator $\mathbf{T}$. We have the following lemma which is an extension of the [33, Lemma A.2] to the inertial case.

**Lemma E.1.** Let $\mathbf{T}$ be defined as in (E1) then,

(a) For all $(x_2, x_1) \in \mathbb{R}^{2n}$, $\det \mathbf{DT}(x_2, x_1) \neq 0$,

(b) The set of strict saddle points is contained in the following set

$$U_{\mathbf{T}} \stackrel{\mathrm{def}}{=} \left\{ (x_2, x_1) \in \mathbb{R}^{2n} : \mathbf{T}(x_2, x_1) = \begin{bmatrix} x_2 \\ x_1 \end{bmatrix}, \max_i |\lambda_i(\mathbf{DT}(x_1, x_2))| > 1 \right\}, \quad (\mathrm{E2})$$

$$= \left\{ (x, x) \in \mathbb{R}^{2n} : \mathbf{T}(x, x) = \begin{bmatrix} x \\ x \end{bmatrix}, \max_i |\lambda_i(\mathbf{DT}(x, x))| > 1 \right\}. \quad (\mathrm{E3})$$

*Proof*

(a) Since $\psi$ is a $C^2$ function, and thus $\nabla\psi$ is $C^1$, and as $\psi$ is strongly convex, the inverse function theorem ensures that $(\nabla\psi)^{-1}$ is a local diffeomorphism. Moreover, $F$ is $C^2$, therefore we define the following for simplicity

$$\mathbf{A}(x_2, x_1) \stackrel{\mathrm{def}}{=} (2a - a^2)\nabla^2\psi^{-1}(\nabla\psi(y(x_2, x_1)) - \gamma\nabla F(y(x_1, x_2))) \\ \qquad (\nabla^2\psi(y(x_2, x_1)) - \gamma\nabla^2 F(y(x_2, y_1))),$$

and

$$\mathbf{B}(x_2, x_1) \stackrel{\mathrm{def}}{=} (1 - a)^2 \nabla^2 \psi^{-1} (\nabla \psi(y(x_2, x_1)) - \gamma \nabla F(y(x_1, x_2))) \\ \qquad (\nabla^2 \psi(y(x_2, x_1)) - \gamma \nabla^2 F(y(x_2, y_1))),$$

then one can write that

$$\mathbf{DT}(x_2, x_1) = \begin{bmatrix} \mathbf{A}(x_2, x_1) & \mathbf{B}(x_2, x_1) \\ \mathrm{Id} & 0 \end{bmatrix}. \quad (\mathrm{E4})$$

We deduce that for any $(x_2, x_1)$ we have that $\det \mathbf{DT}(x_2, x_1) = \det(-\mathbf{B}(x_2, x_1))$. Therefore, it suffices to show that $\det(-\mathbf{B}(x_2, x_1)) \neq 0$ which hold since $\forall (x_2, x_1)$, $\nabla^2\psi(y(x_2, x_1)) - \gamma\nabla^2 F(y(x_2, x_1))$ is invertible. Indeed we have,

$$\nabla^2\psi(y(x_2, x_1)) - \gamma\nabla^2 F(y(x_2, x_1)) \succ (1 - \gamma L)\nabla^2\psi(y(x_2, x_1)) \succ 0,$$

where we have used the $L$-smooth adaptable property and the strong convexity of $\psi$.

(b) Let $x_*$ be a strict saddle point therefore $(x_*, x_*)$ is a fixed point of $\mathbf{T}$. To have that $(x_*, x_*) \in U_{\mathbf{T}}$ it remains to show that $\mathbf{DT}$ has an eigenvalue of magnitude greater than 1. We have

$$\mathbf{DT}(x_*, x_*) = \begin{bmatrix} (2a - a^2) (\mathrm{Id} - \gamma \nabla^2 \psi(x_*)^{-1} \nabla^2 F(x_*)), & (1 - a)^2 (\mathrm{Id} - \gamma \nabla^2 \psi(x_*)^{-1} \nabla^2 F(x_*)) \\ \mathrm{Id} & 0 \end{bmatrix}.$$

Let us remark that $\mathbf{DT}(x_*, x_*)$ is similar to the following matrix

$$\tilde{\mathbf{T}} = \begin{bmatrix} (2a - a^2) (\mathrm{Id} - \gamma H_{\psi}^{-1/2} \nabla^2 F(x_*) H_{\psi}^{1/2}), & (1 - a)^2 (\mathrm{Id} - \gamma H_{\psi}^{-1/2} \nabla^2 F(x_*) H_{\psi}^{1/2}) \\ \mathrm{Id} & 0 \end{bmatrix},$$
where we have applied the following transformation
$$
\begin{bmatrix}
H_{\psi}^{-1/2} & 0 \\
0 & H_{\psi}^{-1/2}
\end{bmatrix}
D\mathbf{T}(x_{\star}, x_{\star})
\begin{bmatrix}
H_{\psi}^{1/2} & 0 \\
0 & H_{\psi}^{1/2}
\end{bmatrix}.
$$

For $\gamma < 1/L$, the symmetric matrix Id $-\gamma H_{\psi}^{-1/2} \nabla^2 F(x_{\star}) H_{\psi}^{1/2}$ has an eigenvalue of magnitude greater than one, (see Lemma D.3). Let us denote by $v$ the eigenvector associated with this eigenvalue that we denote $\eta > 1$. We claim that the vector $[v_0]$ is an eigenvector associated with the eigenvalue $\eta$ of the matrix $\tilde{\mathbf{T}}$. Indeed, we have

$$
\tilde{\mathbf{T}} \begin{bmatrix} v \\ 0 \end{bmatrix} = \begin{bmatrix} (2a - a^2)v\eta + (1 - a)^2v\eta \\ 0 \end{bmatrix} = \eta \begin{bmatrix} v \\ 0 \end{bmatrix},
$$

which conclude the proof.

<figure><img src="image_8.png" /></figure>

To show our Claim, we combine the global convergence result of Theorem 3.2-(ii) the previous Lemma and the center stable manifold theorem see [43, Corollary 1]. This allows us to say that the set

$$
\left\{ (z_0, z_{-1}) \in \mathbb{R}^{2n} : \lim_{k \to \infty} \mathbf{T}^k((z_0, z_{-1})) \in \text{strisad}(\Phi) \right\}
$$

has measure zero.

E.2 Proof of Lemma 5.5

Let us first observe that since $x_* \in \text{ND}(\Phi)$, using the same arguments as for the finite time identification Lemma 4.3 for any $x = (x_2, x_1) \in \mathbb{R}^{2n}$ near $(x_*, x_*)$, $\mathbf{T}(x_2, x_1) \in \tilde{\mathcal{W}}_{x_*}$. The subsequent portion of the proof faithfully adheres to the perturbation analysis used in [28, Theorem 4.1] with a minor adjustment to the inertial Bregman case.

Let $\tilde{G}: \mathbb{R}^n \rightarrow \mathbb{R}$ be any $C^2$-smooth extension (representative) of $G$ on the neighborhood of $x_*$ in $W_{x_*}$, consider the following problem defining $P_1 \mathbf{T}$ near $(x_*, x_*)$,

$$
\min_{u \in W_{x_*}} \Upsilon(x, u) \stackrel{\text{def}}{=} \left\{ F(y(x)) + \tilde{G}(u) + \langle \nabla F(y(x)), u - y(x) \rangle + \frac{1}{\gamma} D_{\psi}(u, y(x)) \right\}. (\tilde{P}_x)
$$

We can write the map $\mathbf{T}$ as

$$
\mathbf{T}(x) = \left[ \begin{array}{c} \min_{u \in W_{x_*}} \Upsilon(x, u) \\ x \end{array} \right],
$$

To apply the perturbation result [62, Theorem 3.1] to $\Upsilon$. We need a quadratic growth condition which we got using the fact that $\psi$ is strongly convex and thus $u \mapsto \Upsilon(x, u)$ is also $\frac{\sigma_\psi}{\gamma}$- strongly convex. We also need to check a level-boundedness condition. This comes from a sufficient condition see [28, Lemma 2.4]. $u(x)$, the minimizer of $\Upsilon(z, \cdot)$ is a continuous map in near $(x_*, x_*)$ then we apply the perturbation analysis [62, Theorem 3.1] with the following Lagrangian function

$$
\mathcal{L}(x, u, \lambda) \stackrel{\text{def}}{=} \Upsilon(x, u) + \langle W(u), \lambda \rangle,
$$

where $\lambda$ is the vector of Lagrange multipliers.
Since $u \mapsto \Upsilon(x, u)$ is a strongly convex function and $u(x_*, x_*) = x_*$ is the unique minimizer on $W_{x_*}$. Thus there exists optimal Lagrange multipliers $\bar{\lambda}$ such that

$$ \nabla_u \mathcal{L}((x_*, x_*), x_*, \bar{\lambda}) = \nabla \tilde{G}(x_*) + \nabla F(x_*) + \sum_{i=1}^{n} \bar{\lambda}_i \nabla W_i(x_*) = 0. $$

From [62, Theorem 3.1], we get that locally for any $x$ near $(x_*, x_*)$, $u(x) = P_1 T(x)$ is a $C^1$-smooth map and we deduce that for any $h = (h_1, h_2) \in \mathbb{R}^{2n}$ with $\|h\| \to 0$,

$$ \langle \nabla P_1 T(x_*, x_*), h \rangle = \underset{v \in T_{W_{x_*}}}{\operatorname{argmin}} 2 \langle \nabla_{zu}^2 \mathcal{L}((x_*, x_*), x_*, \bar{\lambda}) v, h \rangle + \langle \nabla_{uu}^2 \mathcal{L}((x_*, x_*), x_*, \bar{\lambda}) v, v \rangle. \quad (E5) $$

The expression of the Hessians of the Lagrangian is of the form $2n \times n$

$$ \begin{aligned} \mathcal{H}_{zu} &\stackrel{\text{def}}{=} \nabla_{zu}^2 \mathcal{L}((x_*, x_*), x_*, \bar{\lambda}) \\ &= \left[ (2-a) \left( \nabla^2 F(x_*) - \frac{1}{\gamma} \nabla^2 \psi(x_*) \right), (1-a)^2 \left( \nabla^2 F(x_*) - \frac{1}{\gamma} \nabla^2 \psi(x_*) \right) \right], \end{aligned} $$

and

$$ \mathcal{H}_{uu} \stackrel{\text{def}}{=} \nabla_{uu}^2 \mathcal{L}((x_*, x_*), x_*, \bar{\lambda}) = \nabla^2 \tilde{G}(x_*) + \frac{1}{\gamma} \nabla^2 \psi(x_*) + \sum_{i=1}^{n} \bar{\lambda}_i \nabla_{uu}^2 W_i(x_*). $$

The minimization problem (E5) is a quadratic form over the tangent space $T_{W_{x_*}}$ thus the minimizer $\bar{v}$ must satisfy the following condition

$$ P_{T_{W_{x_*}}} \mathcal{H}_{uu} P_{T_{W_{x_*}}} \bar{v} + P_{T_{W_{x_*}}} P_1 \mathcal{H}_{zu}^{\top} P_{T_{W_{x_*}}} h_1 + P_{T_{W_{x_*}}} P_2 \mathcal{H}_{zu}^{\top} P_{T_{W_{x_*}}} h_2 = 0. \quad (E6) $$

Let us denote $\tilde{\mathcal{H}}_{uu} = P_{T_{W_{x_*}}} \mathcal{H}_{uu} P_{T_{W_{x_*}}}$, $\tilde{\mathcal{H}}_{zu}^2 = P_{T_{W_{x_*}}} P_1 \mathcal{H}_{zu}^{\top} P_{T_{W_{x_*}}}$ and $\tilde{\mathcal{H}}_{zu}^1 = P_{T_{W_{x_*}}} P_2 \mathcal{H}_{zu}^{\top} P_{T_{W_{x_*}}}$. (E6) becomes

$$ \tilde{\mathcal{H}}_{uu} \bar{v} + \tilde{\mathcal{H}}_{zu}^2 h_1 + \tilde{\mathcal{H}}_{zu}^1 h_2 = 0. $$

Since $x_*$ is the unique minimizer of $\Upsilon((x_*, x_*), \cdot)$ and we also observe that $\Upsilon((x_*, x_*), \cdot)$ has the same active manifold $W_{x_*}$ as a sum of $G$ and a smooth function this implies that $\Upsilon((x_*, x_*), \cdot)$ has at least a quadratic growth near $x_*$. Therefore $\tilde{\mathcal{H}}_{uu}$ is symmetric positive definite and invertible. Then we solve (E5) to get that

$$ \langle \nabla P_1 T((x_*, x_*)), h \rangle = -\tilde{\mathcal{H}}_{uu}^{-1} \tilde{\mathcal{H}}_{zu}^2 h_1 - \tilde{\mathcal{H}}_{uu}^{-1} \tilde{\mathcal{H}}_{zu}^1 h_2, \quad \forall h_1, h_2 \in T_{W_{x_*}}. $$

At this point, we get that $T$ is a $C^1$-smooth map. We immediately deduce that

$$ \langle \nabla T((x_*, x_*)), h \rangle = \left[ -\tilde{\mathcal{H}}_{uu}^{-1} \tilde{\mathcal{H}}_{zu}^2 h_2 - \tilde{\mathcal{H}}_{uu}^{-1} \tilde{\mathcal{H}}_{zu}^1 h_1 \right], \quad \forall h_1, h_2 \in T_{W_{x_*}}. $$
It remains now to show that $\nabla T(x_*, x_*)$ has at least one eigenvalue greater than one. Let us first show that when $x_* \in \text{Actstrisad}(\Phi)$, $-\tilde{\mathcal{H}}_{uu}^{-1}\tilde{\mathcal{H}}_{zu}^2$ has an eigenvalue greater than one. Let $\eta$ be an real eigenvalue associated to an eigenvector $v \in \mathcal{T}_{W_{x_*}}$, i.e.

$$-\tilde{\mathcal{H}}_{uu}^{-1}\tilde{\mathcal{H}}_{zu}^2 v = \eta v \iff (\eta\tilde{\mathcal{H}}_{uu} + \tilde{\mathcal{H}}_{zu}^2)v = 0.$$

Let us observe that since $x_*$ is an active strict saddle point for the problem, then $\tilde{\mathcal{H}}_{uu} + \tilde{\mathcal{H}}_{zu}^2$ has a strict negative eigenvalue. Indeed,

$$\tilde{\mathcal{H}}_{uu} + \tilde{\mathcal{H}}_{zu}^2 = P_{\mathcal{T}_{W_{x_*}}} \left( \nabla^2 \tilde{G}(x_*) + \sum_{i=1}^n \bar{\lambda}_i \nabla_{uu}^2 W_i(x_*) + (1-a)^2 \left( \nabla^2 F(x_*) - \frac{1}{\gamma} \nabla^2 \psi(x_*) \right) \right) P_{\mathcal{T}_{W_{x_*}}}.$$

Combining with the fact that $\tilde{\mathcal{H}}_{uu}$ is positive definite means that there exists $\eta > 1$ such that the matrix $\eta\tilde{\mathcal{H}}_{uu} + \tilde{\mathcal{H}}_{zu}$ is singular. By construction, we get that $\eta > 1$ is an eigenvalue of $\nabla T$ associated to the eigenvector $[v \ 0]$ of $T$.

# Appendix F Riemannian Geometry and Partial Smoothness

## F.1 Riemannian geometry

Let $\mathcal{M}$ be a $C^2$-smooth embedded submanifold of $\mathbb{R}^n$ around a point $x \in \mathbb{R}^n$. With some abuse of terminology, we will state $C^2$-manifold instead of $C^2$-smooth embedded submanifold of $\mathbb{R}^n$. The natural embedding of a submanifold $\mathcal{M}$ into $\mathbb{R}^n$ permits to define a Riemannian structure and to introduce geodesics on $\mathcal{M}$, and we simply say $\mathcal{M}$ is a Riemannian manifold. We denote respectively $\mathcal{T}_{\mathcal{M}}(x)$ and $\mathcal{N}_{\mathcal{M}}(x)$ the tangent and normal space of $\mathcal{M}$ at $x \in \mathcal{M}$.

### Exponential map

Geodesics generalize the concept of straight lines from linear spaces to manifolds. It is a smooth curve from an interval of $\mathbb{R}$ to $\mathcal{M}$, with intrinsic acceleration normal everywhere to $\mathcal{M}$. Roughly speaking, it is locally the shortest path between two points on $\mathcal{M}$. Let denote by $\mathfrak{g}(t; x, h)$ the value at $t \in \mathbb{R}$ of the geodesic starting $\mathfrak{g}(0; x, h) = x \in \mathcal{M}$ with velocity $\dot{\mathfrak{g}}(t; x, h) = \frac{d\mathfrak{g}(t; x, h)}{dt} = h \in \mathcal{T}_{\mathcal{M}}(x)$ (uniquely defined). It is important to realize that for every $h \in \mathcal{T}_{\mathcal{M}}(x)$ there exists an interval $I$ around 0 and a unique geodesic $\mathfrak{g}(t; x, h) : I \to \mathcal{M}$ such that $\mathfrak{g}(0; x, h) = x$ and $\dot{\mathfrak{g}}(0; x, h) = t$. The mapping $\text{Exp}_x : \mathcal{T}_{\mathcal{M}}(x) \to \mathcal{M}$, $h \mapsto \text{Exp}_x(h) = \mathfrak{g}(1; x, h)$, is called the Exponential map. Given $x, x' \in \mathcal{M}$, and a direction $h \in \mathcal{T}_{\mathcal{M}}(x)$ we are want such map to fulfill $\text{Exp}_x(h) = x' = \mathfrak{g}(1; x, h)$.

### Parallel translation

Given two points $x, x' \in \mathcal{M}$ let $\mathcal{T}_{\mathcal{M}}(x), \mathcal{T}_{\mathcal{M}}(x')$ be the corresponding tangent spaces. Define $\tau : \mathcal{T}_{\mathcal{M}}(x) \to \mathcal{T}_{\mathcal{M}}(x')$, the parallel translation along the unique geodesic joining $x$ to $x'$, which is an isomorphism and isometry with respect to the Riemannian metric.
### Riemannian gradient and Hessian

For a vector $v \in \mathcal{N}_M(x)$, the Weingarten map of $M$ at $x$ is the operator $\mathfrak{W}_x(\cdot, v) : \mathcal{T}_M(x) \to \mathcal{T}_M(x)$ defined by:

$$ \mathfrak{W}_x(h, v) = -P_{\mathcal{T}_M(x)} dV[h], $$

where $V$ is any local extension of $v$ to a normal vector field on $M$. The definition does not depend of the choice of the extension $V$. The Weingarten map as defined above is a symmetric linear operator which is closely tied to the second fundamental form of $M$ ([22, Definition 5.48]). Let $g$ be a real-valued function which is $C^2$ along $M$ around $x$. The covariant gradient of $g$ at $x' \in M$ is the vector denoted $\nabla_M g(x') \in \mathcal{T}_M(x')$ defined by:

$$ \langle \nabla_M g(x'), h \rangle = \frac{d}{dt} g(P_M(x' + th))|_{t=0}, \quad \forall h \in \mathcal{T}_M(x'), $$

where $P_M$ is the projection onto $M$. The covariant Hessian of $g$ at $x'$ is the symmetric linear mapping $\nabla_M^2 g(x')$ from $\mathcal{T}_M(x')$ to itself which is defined as

$$ \langle \nabla_M^2 g(x') h; h \rangle = \frac{d^2}{dt^2} g(P_M(x' + th))|_{t=0}, \quad \forall h \in \mathcal{T}_M(x'). $$

This definition agrees with the definition using geodesics or connections. Now, assume that $M$ is a Riemannian embedded submanifold of $\mathbb{R}^n$ and $g$ has a $C^2$-smooth restriction on $M$. This can be characterization by the existence of a $C^2$-smooth extension (representative) of $g$ i.e. a $C^2$-smooth function $\tilde{g}$ on $\mathbb{R}^n$ such that $\tilde{g}$ agrees with $g$ on $M$. Thus, the Riemannian gradient $\nabla_M g(x')$ is also given by

$$ \nabla_M g(x') = P_{\mathcal{T}_M(x')} \nabla \tilde{g}(x'), $$

and $\forall h \in \mathcal{T}_M(x')$, the Riemannian Hessian reads

$$ \begin{aligned} \nabla_{\mathcal{T}_M}^2 g(x') h &= P_{\mathcal{T}_M(x')} d(\nabla_{\mathcal{T}_M} g(x')) [h] = P_{\mathcal{T}_M(x')} d(\nabla_{\mathcal{T}_M} \tilde{g}(x')) [h] \\ &= P_{\mathcal{T}_M(x')} \nabla^2 \tilde{g}(x') h + \mathfrak{W}_{x'}(h, P_{\mathcal{N}_{\mathcal{T}_M}(x')} \nabla \tilde{g}(x')). \end{aligned} $$

If $M$ is affine or linear subspace of $\mathbb{R}^n$, then $M = x + \mathcal{T}_M(x)$ and $\mathfrak{W}_{x'}(h, P_{\mathcal{N}_{\mathcal{T}_M}(x')} \nabla \tilde{g}(x')) = 0$ and finally

$$ \nabla_{\mathcal{T}_M}^2 g(x') = P_{\mathcal{T}_M(x')} \nabla^2 \tilde{g}(x') P_{\mathcal{T}_M(x')}. $$

The next two lemmas will be instrumental when analyzing the local convergence behaviour of our inertial algorithm. We refer to [45, Section 2.6] for their proofs.

**Lemma F.1.** Let $x \in M$ and $x_k$ a sequence converging to $x$ in $M$. Denote $\tau_k : \mathcal{T}_M(x_k) \to \mathcal{T}_M(x_k)$ be the parallel translation along the unique geodesic joining $x$ to $x_k$. Then, for any bounded vector $u \in \mathbb{R}^n$, we have:

$$ \left( \frac{1}{\tau_k} P_{\tau_M(x_k)} - P_{\tau_M(x)} \right) u = o (\|u\|). \qquad (F1) $$
**Lemma F.2.** Let $x, x'$ be two close points in $\mathcal{M}$, denote $\tau: \mathcal{T}_{\mathcal{M}}(x) \to \mathcal{T}_{\mathcal{M}}(x')$ the parallel translation along the unique geodesic joining $x$ too $x'$. The Riemannian Taylor expansion of $g \in C^2(\mathcal{M})$ around $x$ reads,

$$ \frac{1}{\tau} \nabla_{\mathcal{M}} g(x') = \nabla_{\mathcal{M}} g(x) + \nabla_{\mathcal{M}}^2 g(x) P_{\mathcal{T}_{\mathcal{M}}(x)}(x' - x) + o(\|x' - x\|). \quad (F2) $$

## F.2 Partial Smoothness

Introduced by Lewis in [44], “Partial smoothness” captures the characteristics of the geometry of nonsmooth functions. It axiomatizes the notion of active/identifiable submanifold or identifiable surfaces in [71]. A partly smooth function is smooth along the identifiable submanifold and sharp transversally to the manifold. Therefore, the behaviour of the function of its minimizers depends essentially on its restriction to this manifold, hence offering a powerful framework for algorithmic and sensitivity analysis theory.

**Definition F.3 (Partly smooth function).** A function $g \in \Gamma_0(\mathbb{R}^n)$ is $C^2$-partly smooth at a point $x$ relative to the set $\mathcal{M}$ containing $x$, if $\partial g(x) \neq \emptyset$ and $\mathcal{M}$ is an embedded $C^2$-smooth submanifold and there exists a neighborhood $\mathcal{V}_x$ of $x$ such that the following properties hold

(i) (Smoothness) the restriction $g|_{\mathcal{M}}$ is a $C^2$ function in the neighborhood $\mathcal{V}_x$;

(ii) (Sharpness) The affine hull of $\partial g(x)$ is a translation of the space $\mathcal{N}_{\mathcal{M}}(x)$, i.e.

$$ S_x \stackrel{\text{def}}{=} \text{par}(\partial g(x)) = \mathcal{N}_{\mathcal{M}}(x) \Leftrightarrow T_x \stackrel{\text{def}}{=} \text{par}(\partial g(x))^{\perp} = \mathcal{T}_{\mathcal{M}}(x). $$

(iii) (Continuity) The set-valued mapping $\partial g$ is continuous at $x$ relative to $\mathcal{M}$.

Observe that $S_x = T_x^{\perp}$ by definition. Throughout the rest of the work, we denote the class of $C^2$-partly smooth function at $x$ relative to $\mathcal{M}$ by $\text{PSF}_x(\mathcal{M})$.

Owing to the definition of partial smoothness, we have the following facts.

**Fact F.4. (Local normal sharpness)** If $g \in \text{PSF}_x(\mathcal{M})$, then for all point $x' \in \mathcal{M}$ near $x$ we have $\mathcal{T}_{\mathcal{M}}(x') = T_{x'}$. In particular when $\mathcal{M}$ is affine or linear, then $T_{x'} = T_x$.

**Fact F.5.** If $g \in \text{PSF}_x(\mathcal{M})$, then for all $x' \in \mathcal{M}$ near $x$ we have

$$ \nabla_{\mathcal{M}} g(x') = P_{T_{x'}} (\partial g(x')), $$

and this does not depend on the smooth representation of $g$ on $\mathcal{M}$. In turn, for all $h \in T_{x'}$,

$$ \nabla_{\mathcal{M}}^2 g(x')h = P_{T_{x'}} \nabla^2 \tilde{g}(x')h + \mathfrak{W}_{x'}(h, P_{T_{x'}^{\perp}} \nabla \tilde{g}(x')), $$

where $\tilde{g}$ is a smooth representative of $g$ on $\mathcal{M}$ and $\mathfrak{W}_{x'}(\cdot, \cdot)$ is the Weingarten map of $\mathcal{M}$ at $x$.

**Acknowledgments A.** The author would like to thank Prof. Jalal Fadili for the discussions on the topic and the guidance during the author's Ph.D. thesis. Besides, he want to thank the French National Research Agency (ANR) for funding the project FIRST (ANR-19-CE42-0009).
References

[1] F. Alvarez. On the minimizing property of a second order dissipative system in Hilbert spaces. *SIAM J. Control Optim.*, 38(4):1102–1119, 2000.

[2] F. Alvarez and H. Attouch. An inertial proximal method for maximal monotone operators via discretization of a nonlinear oscillator with damping. *Set-Valued Analysis*, 9(1):3–11, 2001.

[3] H. Attouch, J. Bolte, and B. Svaiter. Convergence of descent methods for semi-algebraic and tame problems: Proximal algorithms, forward-backward splitting, and regularized Gauss-Seidel methods. *Math. Program.*, 137(1-2):91–129, 2013.

[4] H. Attouch and J. Peypouquet. The rate of convergence of Nesterov's accelerated forward-backward method is actually faster than $o(k^{-2})$. *SIAM J. Optim.*, 26(3):1824–1834, 2016.

[5] H. Attouch, J. Peypouquet, and P. Redont. A dynamical approach to an inertial forward-backward algorithm for convex minimization. *SIAM J. Optim.*, 24(1):232–256, 2014.

[6] J.-P. Aubin and I. Ekeland. *Applied Nonlinear Analysis*. Elsevier, 1984.

[7] A. Auslender and M. Teboulle. Interior gradient and proximal methods for convex and conic optimization. *SIAM J. Optim.*, 16(3):697–725, 2006.

[8] H. Bauschke, J. Bolte, J. Chen, M. Teboulle, and X. Wang. On linear convergence of non-euclidean gradient methods without strong convexity and Lipschitz gradient continuity. *Journal of Optimization Theory and Applications*, 182(3):1068–1087, 2019.

[9] H. Bauschke, J. Bolte, and M. Teboulle. A descent lemma beyond Lipschitz gradient continuity: First-order methods revisited and applications. *Mathematics of Operations Research*, page 20, 2016.

[10] H. Bauschke, J. Borwein, and P. Combettes. Bregman monotone optimization algorithms. *SIAM Journal on Control and Optimization*, 42(2):596–636, 2003.

[11] H. Bauschke and P. Combettes. *Convex analysis and monotone operator theory in Hilbert spaces*. Springer, 2011.

[12] A. Beck and M. Teboulle. A fast iterative shrinkage-thresholding algorithm for linear inverse problems. *SIAM J. Imaging Sci.*, 2(1):183–202, 2009.

[13] P. Bianchi, W. Hachem, and S. Schechtman. Stochastic subgradient descent escapes active strict saddles on weakly convex Functions. *Mathematics of Operations Research*, 2023.

[14] B. Birnbaum, N. Devanur, and L. Xiao. Distributed algorithms via gradient descent for fisher markets. In *Proceedings of the 12th ACM conference on Electronic commerce*, pages 127–136, 2011.

[15] J. Bolte, A. Daniilidis, and A. Lewis. The Łojasiewicz inequality for nonsmooth subanalytic functions with applications to subgradient dynamical systems. *SIAM J. Optim.*, 17(4):1205–1223, 2007.

[16] J. Bolte, A. Daniilidis, A. Lewis, and M. Shiota. Clarke subgradients of stratifiable functions. *SIAM Journal on Optimization*, 18(2):556–572, 2007.

[17] J. Bolte, A. Daniilidis, O. Ley, and L. Mazet. Characterizations of Łojasiewicz inequalities: Subgradient flows, talweg, convexity. *Trans. Amer. Math. Soc.*, 362(06):3319–3363, December 2009.

[18] J. Bolte, S. Sabach, and M. Teboulle. Proximal alternating linearized minimization for nonconvex and nonsmooth problems. *Math. Program.*, 146(1-2):459–494, August 2014.
[19] J. Bolte, S. Sabach, and M. Teboulle. Nonconvex lagrangian-based optimization: Monitoring schemes and global convergence. Mathematics of Operations Research, 43(4):1210–1232, 2018.

[20] J. Bolte, S. Sabach, M. Teboulle, and Y. Vaisbourd. First order methods beyond convexity and Lipschitz gradient continuity with applications to quadratic inverse problems. SIAM J. Optim., 28(3):2131–2151, 2018.

[21] R. Bøt, R. Çsetnek, and S. László. An inertial forward-backward algorithm for the minimization of the sum of two nonconvex functions. EURO Journal on Computational Optimization, 4(1):3–25, 2016.

[22] N. Boumal. An introduction to optimization on smooth manifolds. Cambridge University Press, 2023.

[23] O. Brandière and M. Duflo. Les algorithmes stochastiques contournent-ils les pièges ? Annales de l'I.H.P. Probabilités et statistiques, 32:395–427, 1996.

[24] J. Burke and J.-J. More. On the identification of active constraints. SIAM Journal on Numerical Analysis, 25(5):1197–1211, 1988.

[25] A. Chambolle and J. Darbon. A parametric maximum flow approach for discrete total variation regularization. In Image Processing and Analysis with Graphs. CRC Press, 2012.

[26] A. Chambolle and C. Dossal. On the convergence of the iterates of the “fast iterative shrinkage/thresholding algorithm”. J Optim Theory Appl, 166(3):968–982, 2015.

[27] G. Chen and M. Teboulle. Convergence analysis of a proximal-like minimization algorithm using Bregman functions. SIAM J. Optim., 3(3):538–543, 1993.

[28] D. Davis and D. Drusvyatskiy. Proximal methods avoid active strict saddles of weakly convex functions. Found Comput Math, 22(2):561–606, 2022.

[29] R. Dragomir, A. Taylor, A. d'Aspremont, and J. Bolte. Optimal complexity and certification of Bregman first-order methods. Math. Program., 194(1):41–83, 2022.

[30] D. Drusvyatskiy, A. Ioffe, and A. Lewis. Generic minimizing behavior in semi-algebraic optimization. SIAM J. Optim., 26(1):513–534, 2016.

[31] D. Drusvyatskiy and A. Lewis. Optimality, identifiability, and sensitivity. Math. Program., 147(1):467–498, 2014.

[32] J. Dunn. On the convergence of projected gradient processes to singular critical points. J Optim Theory Appl, 55(2):203–216, November 1987.

[33] J.-J. Godeme, J. Fadili, X. Buet, M. Zerrad, M. Lequime, and C. Amra. Provable phase retrieval with mirror descent. SIAM J. Imaging Sci., 16(3):1106–1141, 2023.

[34] X. Goudou and J. Munier. The gradient and heavy-ball with friction dynamical systems: the quasiconvex case. Math. Program., 116(1):173–191, 2009.

[35] F. Hanzely, P. Richtárik, and L. Xiao. Accelerated Bregman proximal gradient methods for relatively smooth convex optimization. Comput Optim Appl, 79(2):405–440, 2021.

[36] W. Hare. Identifying active manifolds in regularization problems. In Fixed-Point Algorithms for Inverse Problems in Science and Engineering, Springer Optimization and Its Applications, pages 261–271. Springer, New York, NY, 2011.

[37] W. Hare and A. Lewis. Identifying active constraints via partial smoothness and prox regularity. Journal of Convex Analysis, 2004.

[38] W. Hare and A. Lewis. Identifying active manifolds. Algorithms Operations Research, pages 75–82, 2007.
[39] L. Hien, N. Gillis, and P. Patrinos. Inertial block proximal methods for non-convex non-smooth optimization, 2020.

[40] A. Ioffe. *Variational Analysis of Regular Mappings*. Springer Monographs in Mathematics. Springer International Publishing, 2017.

[41] K. Kurdyka. On gradients of functions definable in o-minimal structures. *Annales de l'institut Fourier*, 48(3):769-783, 1998.

[42] E. Laude, P. Ochs, and D. Cremers. Bregman proximal mappings and Bregman-Moreau envelopes under relative prox-regularity. *J. Optim Theory Appl*, 184(3):724-761, 2020.

[43] J. Lee, I. Panageas, G. Piliouras, M. Simchowitz, M. Jordan, and B. Recht. First-order methods almost always avoid strict saddle points. *Math. Program.*, 176(1-2):311-337, 2019.

[44] A. Lewis. Active sets, nonsmoothness, and sensitivity. *SIAM J. Optim.*, 13(3):702-725, 2002.

[45] J. Liang. *Convergence rates of first-order operator splitting methods*. PhD thesis, Université de Caen, 2016.

[46] J. Liang, J. Fadili, and G. Peyré. Local linear convergence of forward-backward under partial smoothness. *In: NeuIPS 2014*, pages 1970–1978, 2014.

[47] J. Liang, J. Fadili, and G. Peyré. A Multi-step Inertial Forward-Backward Splitting Method for Non-convex Optimization. In *Advances in Neural Information Processing Systems*, volume 29. Curran Associates, Inc., 2016.

[48] J. Liang, J. Fadili, and G. Peyré. Activity identification and local linear convergence of forward-backward-type methods. *SIAM J. Optim.*, 27(1):408–437, 2017.

[49] J. Liang, J. Fadili, and G. Peyré. Local linear convergence analysis of Primal-Dual splitting methods. *Optimization*, 67(6):821–853, June 2018.

[50] J. Liang, J. Fadili, G. Peyré, and R. Luke. Activity identification and local linear convergence of Douglas-Rachford/ADMM under partial smoothness. In *Scale Space and Variational Methods in Computer Vision*, Lecture Notes in Computer Science, pages 642-653. Springer International Publishing, 2015.

[51] S. Łojasiewicz. Une propriété topologique des sous-ensembles analytiques réels. In *Les Équations aux Dérivées Partielles*, pages 87-89. Editions du Centre National de la Recherche Scientifique, 1963.

[52] S. Łojasiewicz. Ensembles semi-analytiques. *Lectures Notes IHES (Bures-sur-Yvette)*, 1965.

[53] D. Lorenz and T. Pock. An inertial forward-backward algorithm for monotone inclusions. *J. Math Imaging Vis*, 51(2):311-325, 2015.

[54] H. Lu, R. Freund, and Y. Nesterov. Relatively smooth convex optimization by first-order methods, and applications. *SIAM Journal on Optimization*, 28(1):333-354, 2018.

[55] A. Moudafi and M. Oliny. Convergence of a splitting inertial proximal method for monotone operators. *Journal of Computational and Applied Mathematics*, 155(2), 2003.

[56] M. Mukkamala, P. Ochs, T. Pock, and S. Sabach. Convex-concave backtracking for inertial Bregman proximal gradient algorithms in nonconvex optimization. *SIAM Journal on Mathematics of Data Science*, 2(3):658-682, 2020.

[57] Y. Nesterov. A method for solving the convex programming problem with convergence rate O(1/k^2). *Proceedings of the USSR Academy of Sciences*, 1983.
[58] I. Panageas and G. Piliouras. Gradient descent only converges to minimizers: non-isolated critical points and invariant regions, 2017.

[59] R. Pemantle. Nonconvergence to unstable points in Urn models and stochastic Approximations. *The annals of probability*, 18(2):698–712, 1990.

[60] B. Polyak. Some methods of speeding up the convergence of iteration methods. *USSR Computational Mathematics and Mathematical Physics*, 4(5):1–17, 1964.

[61] R. Rockafellar. *Convex Analysis*. Princeton University Press, 1970.

[62] A. Shapiro. Second order sensitivity analysis and asymptotic theory of parametrized nonlinear programs. *Mathematical Programming*, 33(3):280–299, 1985.

[63] M. Shub. *Global Stability of Dynamical Systems*. Springer, 1987.

[64] A. Silveti-Falls, C. Molinari, and J. Fadili. A stochastic Bregman primal-dual splitting algorithm for composite optimization. *Pure and Applied Functional Analysis (special issue in honor of L. Bregman)*, 2022.

[65] G. Steidl, J. Weickert, T. Brox, P. Mrázek, and M. Welk. On the equivalence of soft wavelet shrinkage, total variation diffusion, total variation regularization, and sides. *SIAM J. Numer. Anal.*, 42(2):686–713, 2004.

[66] M. Teboulle. A simplified view of first order methods for optimization. *Mathematical Programming*, 170(1):67–96, 2018.

[67] Q. Van Nguyen. Forward-backward splitting with bregman distances. *Vietnam Journal of Mathematics*, 45(3):519–539, 2017.

[68] I. Waldspurger, A. d'Aspremont, and S. Mallat. Phase recovery, MaxCut and complex semidefinite programming. *Math. Program.*, 149(1):47–81, 2015.

[69] G. Wang, Georgios B. Giannakis, J. Chen, and M. Akçakaya. SPARTA: Sparse phase retrieval via truncated amplitude flow. In *2017 IEEE International conference on acoustics, speech and signal processing (ICASSP)*, pages 3974–3978, 2017.

[70] B. Wen, X. Chen, and T. Pong. Linear convergence of proximal gradient algorithm with extrapolation for a class of nonconvex nonsmooth minimization problems. *SIAM Journal on Optimization*, 27(1):124–145, 2017.

[71] S. Wright. Identifiable Surfaces in Constrained Optimization. *SIAM J. Control Optim.*, 31(4):1063–1079, 1993.

[72] Z. Wu, C. Li, M. Li, and A. Lim. Inertial proximal gradient methods with bregman regularization for a class of nonconvex optimization problems. *Journal of Global Optimization*, 79(3):617–644, 2021.