Yoel Drori

George Lan

Yurii Nesterov

Edouard Pauwels

Cristobal Guzman

Grégoire Merci À Mathieu

Thomas E Hadrien

Rémi

Adrien Yann

Yana, Loucas, Raphaël, Oumayma Vivien Francis

Ale, Céline, Hans Justin Bruno

Hélène Pierre-Louis Wilson

En apprentissage statistique et traitement du signal, de nombreuses tâches se formulent sous la forme d'un problème d'optimisation de grande taille. Dans ce contexte, les méthodes du premier ordre, qui utilisent uniquement l'information apportée par le gradient de la fonction objectif, sont privilégiées en raison de leur faible coup par itération et de leur simplicité. Nous étudions dans cette thèse les méthodes du premier ordre à distance de Bregman, qui constituent une généralisation de la célèbre méthode de descente de gradient. Cette généralisation consiste à remplacer la distance euclidienne par une distance plus générale, dite de Bregman, générée par une fonction convexe noyau suffisamment simple. La fonction noyau est choisie de manière à être adaptée à la géométrie de la fonction objectif au travers de la condition de régularité relative, introduite en 2017 par Bauschke, Bolte et Teboulle. Depuis son apparition, cette condition a fait naître de nouvelles perspectives en optimisation du premier ordre.

Tout d'abord, nous appliquons les méthodes de Bregman aux problèmes d'optimisation sur des espaces de matrices de rang faible. En exploitant la structure matricielle et en utilisant la propriété de régularité relative, nous proposons des noyaux de Bregman qui permettent d'améliorer la performance numérique par rapport aux méthodes euclidiennes.

Ensuite, nous nous penchons sur la complexité théorique de ces algorithmes. Un des problèmes les plus importants est de déterminer s'il existe une version accélérée de l'algorithme de gradient de Bregman qui possède un meilleur taux de convergence. Dans le cas général, nous démontrons que la réponse est négative : la complexité de la descente de gradient de Bregman standard ne peux pas être améliorée pour des noyaux génériques. La preuve repose sur un contre-exemple pathologique qui a été découvert au travers de méthodes d'analyses de pire cas par ordinateur. Nous évoquons aussi une tentative pour obtenir des résultats positifs d'accélération en spécialisant cette analyse dans le contexte plus restreint de la géométrie entropique.

Enfin, nous étudions la version stochastique de l'algorithme de Bregman pour minimiser des fonctions sous la forme d'espérance, ainsi que des méthodes de réduction de variance lorsque la fonction objectif est une somme finie. i Je tiens en premier lieu à remercier chaleureusement mes deux directeurs de thèse pour leur encadrement bienveillant et amical. J'ai été marqué, Jérôme, par ton dévouement, ta sincérité, ta philosophie, ta passion pour les belles mathématiques et ton exigence d'une recherche de qualité. Alexandre, ta culture scientifique, ton sens de l'initiative et du relationnel, ton énergie et ton enthousiasme à toute épreuve ont été un moteur essentiel de cette thèse.

Je tiens à exprimer mes sincères remerciements aux membres du jury :

Abstract

We study large-scale optimization problems with applications to signal processing and machine learning. Such problems are typically solved with first-order methods that perform iterative updates using the gradient of the objective function. We focus on the class of Bregman firstorder methods, for which the direction of the gradient step is determined by the Bregman divergence induced by a convex kernel function. The choice of the kernel is guided by the relative smoothness condition, which requires the kernel to be compatible with the objective through a descent inequality. This condition was introduced recently by Bauschke, Bolte and Teboulle in 2017 and has opened new perspectives in first-order optimization.

In the first part, we apply Bregman methods to minimization problems on the space of low-rank semidefinite matrices. By leveraging the matrix structure and using the relative smoothness property, we show that well-chosen Bregman kernels allow to improve performance over standard Euclidean methods.

Then, we study the theoretical complexity of these algorithms. An important question is to determine whether there exists an accelerated version of Bregman gradient descent which achieves a better convergence rate in the same setting. In the general case, we show that the answer is negative as the complexity of the standard Bregman gradient method cannot be improved for generic kernels. The proof relies on a pathological example which was discovered by analyzing the worst-case behavior of Bregman methods with a computer-aided technique called performance estimation. We also detail an attempt towards improving the convergence speed in a more restricted setting, by specializing the performance estimation framework to the entropic geometry.

Finally, we study a stochastic variant of Bregman gradient descent for expectation minimization problems, which are pervasive in machine learning, along with variance reduction methods for finite-sum objectives.

Conclusion and Perspectives

Bibliography ix

Contributions and Thesis Outline

Throughout this thesis, we consider the minimization problem

min x2C f (x)
where C is a convex set and f a differentiable function. We focus on the setting where f satisfies a relative smoothness condition with respect to some kernel function h: f (x)  f (y) + hrf (y), x yi + LD h (x, y), 8x, y 2 int C, where D h is the Bregman divergence induced by h. The basic algorithm for solving this problem is Bregman gradient descent (also known as NoLips or mirror descent).

Chapter 1: we introduce the framework of Bregman gradient methods for relatively-smooth problems, and review related work on theoretical aspects, algorithmic extensions and applications. We also provide a more general historical overview on the use of Bregman divergences in optimization.

Chapter 2: in this chapter, we apply the Bregman gradient method to nonconvex low-rank problems. We wish to minimize a L-smooth function F on the space of semidefinite matrices of size n ⇥ n and rank at most r with the nonconvex Burer-Monteiro formulation, which writes

f (X) = F (XX T ) + g(X)
where X 2 R n⇥r , and g is a (possibly nonsmooth) regularization or constraint function. We show that the factorized function f is not globally smooth with respect to the squared Euclidean norm, but it can be proven to be relatively-smooth with respect to well-chosen quartic kernels, allowing to apply Bregman proximal gradient methods. We study two types of quartic kernels. The norm kernel is defined as

h n (X) = ↵ 4 kXk 4 + 2 kXk 2
and can be used for various functions g, such as the sparsity-inducing `1 norm or the indicator set of the nonnegative orthant. We also introduce the Gram kernel,

h G (X) = ↵ 4 kXk 4 + 4 kX T Xk 2 + 2 kXk 2 ,
which is more complex to use as the Bregman gradient iterates can only be computed easily for problems without penalization, i.e., with g ⌘ 0. We show that, by leveraging the lowrank matrix structure, the Gram kernel allows to improve numerical performance for wellconditionned functions F . We provide numerical experiments on two applications: symmetric nonnegative matrix factorization (SymNMF) and Euclidean distance matrix completion.

Chapter 3: we study the theoretical complexity of relatively-smooth convex minimization. We prove a lower bound stating that Bregman gradient descent is an optimal method in this setting, among the general class of methods that use Bregman divergences and basic linear operations. This result demonstrates that the O(1/k) convergence rate of Bregman gradient descent cannot be improved for generic kernel h, and thus that additional regularity assumptions are needed in order to devise a faster algorithm. The worst-case instances involved in the lower bound were inferred from the numerical solution to a performance estimation problem, which we describe in the next chapter.

Chapter 4: in this chapter, we use computer-aided tools for analyzing the worst-case complexity of Bregman methods. This technique, called performance estimation, was initially designed for the study of Euclidean gradient methods. We extend it to the setting of relativelysmooth problems, for computing theoretical guarantees of Bregman methods with generic kernel h. That is, we search for the worst-case behavior among the class of relatively-smooth problem instances of the form B L (C) = {(f, h) : h is a kernel function on C and f is convex and L-smooth relative to h} .

By establishing interpolation conditions, we reduce this task to a finite-dimensional semidefinite program defined on an open set. Then, we use topological arguments to show that this problem is equivalent to a simpler problem on a larger class B L (C), which can be seen as the closure of B L (C) and involves possibly nonsmooth convex functions.

We showcase this approach on several examples. First, solving the dual performance estimation problem allows to infer analytical proofs for two different convergence results of BGD. Then, we are also able to discover the corresponding worst-case functions f and h. Because of the topological fact mentioned above, these worst-case elements are actually nonsmooth functions from B L (C), and can be approached by a sequence of elements of B L (C). This limiting nonsmooth behavior is at the core of the proof of the lower bound in Chapter 3, which was inferred from these numerically generated worst-case examples.

Chapter 5: in this chapter, we describe an attempt to perform computer-aided analyses of Bregman methods in the more restricted setting of a fixed kernel h. Indeed, results from previous chapters suggest that the class B L (C) of general relatively-smooth problems is too large, as functions can approach pathological nonsmooth instances from B L (C). To consider a more realistic setting, we focus on the particular case of the entropic kernel

h e (x) = d X i=1
x (i) log x (i) x (i) , which is one of the most common kernels used in Bregman methods and presents favorable regularity properties. Therefore, we now search for the worst-case of Bregman methods on the class F he L = {f : f is convex and L-smooth relative to h e } . We show that the corresponding performance estimation problem (PEP) can be reduced to a finite-dimensional problem through interpolation conditions, which we establish using the smoothing properties of Bregman-Moreau envelopes. We then formulate this problem as a Chapter 1

Relative Smoothness and Bregman Optimization Methods

Many tasks in computational science, engineering and operations research can be formulated as constrained optimization problems of the form

min x2C f (x) (P)
where C is a convex subset of R d and f is a differentiable function. We are particularly interested in large-scale applications from signal processing, data analysis and machine learning, where the dimension d typically ranges between 10 3 and 10 9 . Problem (P) is generally solved with iterative algorithms, which successively minimize local approximations of the function f based on its derivatives. As higher-order derivatives are too costly to compute for high dimension d, we consider methods that only rely on the information provided by the gradient rf : they constitute the class of first-order methods.

The simplest and most known first-order method is projected gradient descent (GD), which amounts to iterate

x k+1 = ⇧ C (x k k rf (x k )) , (1.1)
where k 0 is the iteration counter, ⇧ C is the Euclidean projection on C and k > 0 is the step size. Equivalently, this update can be written as

x k+1 = argmin u2C f (x k ) + hrf (x k ), u x k i + 1 2 k ku x k k 2 , (GD)
where h•, •i and k•k denote the Euclidean inner product and norm on R d . Thus, gradient descent successively minimizes a surrogate of f consisting of the first-order approximation penalized with a quadratic term which quantifies the inexactness of the linear model far from x k . The efficiency of the method is connected to how well this surrogate approximates f . The classical condition for analyzing this is to assume that the gradient rf is Lipschitz continuous. This ensures that, for k small enough, the inner objective in (GD) is an upper bound on f , and therefore that every iteration decreases the value of the objective. However, in some problems, the quadratic model is a poor approximation of f and a different geometry might be more adapted. This is the purpose of the Bregman gradient descent method, which performs the update

x k+1 = argmin u2C f (x k ) + hrf (x k ), u x k i + 1 k D h (u, x k ), (BGD)
where the quadratic term has been replaced by the more general Bregman divergence D h (x, y) = h(x) h(y) hrh(y), x yi induced by some convex kernel function h. This divergence is nonnegative by convexity of h. Note that taking h(x) = 1 2 kxk 2 yields the Euclidean distance D h (x, y) = 1 2 kx yk 2 . However, choosing a different kernel fuction h can result in a more accurate local model than the standard Euclidean one. This choice is guided by the so-called relative smoothness condition, introduced in Bauschke et al. [2017], which amounts to assume that the inner objective in (BGD) is an upper bound on f . Naturally, the kernel h also needs to be simple enough so that the subproblem defining x k+1 can be solved easily, that is, in closed form or with efficient numerical schemes.

The Bregman gradient method is also known as Mirror Descent, and has been originally proposed by [START_REF] Nemirovski | Problem Complexity and Method Efficiency in Optimization[END_REF] for minimizing convex, possibly nonsmooth, functions. In their context, the kernel h is chosen according to the set C, in order to minimize its diameter with respect to the corresponding Bregman divergence (see the discussion in Section 1.5). In this thesis, we study the more recent relatively-smooth setting, where h is chosen to fit the curvature of the differentiable function f , and use the Bregman gradient descent terminology, although there is no consensus in the litterature.

Relative smoothness and the Bregman gradient method

In this section, we review the general setting for relatively-smooth optimization as well as the main theoretical results. We use standard notation and definitions from convex analysis; see e.g., [START_REF] Rockafellar | Convex Analysis[END_REF], Bauschke and [START_REF] Heinz | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF].

The first step is the choice of a kernel function h adapted to the optimization set C. While the blanket assumptions on h vary between papers, there is usually a set of minimal requirements. The set C is sometimes called the zone of h. The kernel h induces the Bregman divergence D h (x, y) = h(x) h(y) hrh(y), x yi defined for x 2 dom h, y 2 dom rh. Because of strict convexity of h, the Bregman divergence D h (x, y) is nonnegative and equal to zero if and only if x = y. However, it is not a proper distance, as it is not symmetric in general. If h is sufficiently regular, then D h is locally approximated by a quadratic function when x is close to y: D h (x, y) ⇡ 1 2 hr 2 h(y)(x y), x yi when x ! y, where r 2 h denotes the Hessian matrix of h. Equiped with the geometry induced by h, we can now define for x 2 int C the Bregman gradient map with step size > 0:

T (x) = argmin u2C hrf (x), u xi + 1 D h (u, x),
which we assume to be easily computable. Due to Definition 1.1(iii), T is well defined on int C. We now list the most standard examples of kernels.

• The Euclidean kernel h(x) = 1 2 kxk 2 induces the Euclidean distance D h (x, y) = 1 2 kx yk 2 , and T is the projected gradient step (1.1).

• The quadratic kernel is defined as h(x) = 1 2 hAx, xi, where A 2 R d⇥d is a positive definite matrix. Using this kernel corresponds to performing linear preconditioning.

• The entropy kernel h(x) = P d i=1 x (i) log x (i) defined for x 2 R d + with the convention 0 log 0 = 0, induces the Kullback-Leibler divergence D h = D KL with

D KL (x, y) = d X i=1
x (i) log

x (i) y (i)

!

x (i) + y (i) .

In the case when C is the unit simplex {x 2 R d + :

P d i=1 x (i) = 1}
, the gradient map is j) exp( rf (x) (j) )

T (x) (i) = x (i) exp( rf (x) (i) ) P d j=1 x (
for i = 1 . . . d, see, e.g., [START_REF] Beck | Mirror Descent and Nonlinear Projected Subgradient Methods for Convex Optimization[END_REF]. This formula is also known as exponential weight update, or exponentiated gradient [START_REF] Kivinen | Exponentiated Gradient versus Gradient Descent for Linear Predictors[END_REF]].

• The log kernel h(x) = P d i=1 log x (i) for x 2 R d ++ induces the Itakura-Saito divergence D h = D IS with

D IS (x, y) = d X i=1 log x (i) y (i) ! + x (i) y (i) 1 ! . defined for x, y 2 R d ++ . The corresponding Bregman gradient map for C = R d ++ is T (x) (i) = x (i) 1 + x (i) rf (x (i) )
, i = 1 . . . d.

See Section 1.2 for more examples. It should be emphasized that, while a kernel function is differentiable on the interior of its domain, it is not required to be differentiable on the boundary. For instance, the entropy is continuous but not differentiable on the boundary of R d + . Moreover, the domain of h can be closed, such as for the entropy which is finite on R d + , or open, as for the log kernel.

Convex conjugate and mirror formulation. If h is a kernel function, we define its convex conjugate h ⇤ as h ⇤ (y) = sup u2R d hu, yi h(u).

By Assumption (iii), h ⇤ is finite-valued on R d . Moreover, its gradient satisfies for every u 2 R d [Rockafellar, 1970, Sect. 26] rh ⇤ (y) = argmax u2R d hu, yi h(u).

Hence the Bregman gradient map T can be alternatively written as

T (x) = rh ⇤ [rh(x) rf (x)] ,
which is sometimes called the mirror descent formulation (see Section 1.5 for historical comments on this expression).

We now introduce the fundamental notion of relative smoothness, also known as smooth adaptable property, or Lipschitz-like condition [START_REF] Heinz | A Descent Lemma Beyond Lipschitz Gradient Continuity: First-Order Methods Revisited and Applications[END_REF].

Definition 1.2 (Relative smoothness). We say that the function f is smooth relative to h if it is differentiable on int C and if there exists a constant L > 0 such that f (x)  f (y) + hrf (y), x yi + LD h (x, y) (RelSmooth)

for every x, y 2 int C.

An equivalent characterization of relative smoothness is to assume that Lh f is convex on int C, or, if f and h are twice differentiable, to the condition r 2 f (x) Lr 2 h(x) 8x 2 int C, where denotes the semidefinite order. If h is the Euclidean kernel, then (RelSmooth) corresponds to the standard notion of smoothness in optimization, implied by Lipschitz continuity of rf (which we will refer to as L-smoothness to avoid ambiguity).

Similarly, [START_REF] Lu | Relatively-Smooth Convex Optimization by First-Order Methods, and Applications[END_REF] also defined the notion of relative strong convexity.

Definition 1.3 (Relative strong convexity). We say that the function f is strongly convex relative to h if it is differentiable on int C and if there exists a constant µ > 0 such that f (y) + hrf (y), x yi + µD h (x, y)  f (x)

for every x, y 2 int C.

In the same, way, relative strong convexity is equivalent to the function f µh being convex on int C, or to the second-order characterization µr 2 h(x) r 2 f (x) 8x 2 int C.

The ratio  = L/µ is called the relative condition number : a value close to 1 indicates that f is well approximated by the right-hand side of (RelSmooth). Now, note that this right-hand side provides a simple global majorant of f . A natural approach for solving (P) is then to successively minimize this majorant, and this leads to the Bregman Gradient method (Algorithm 1), also known as NoLips [START_REF] Heinz | A Descent Lemma Beyond Lipschitz Gradient Continuity: First-Order Methods Revisited and Applications[END_REF].

Algorithm 1 Bregman Gradient Descent (BGD) / NoLips

Input: Initial point x 0 2 int C, step size . for k = 0,1,. . . do

x k+1 = argmin u2C f (x k ) + hrf (x k ), u x k i + 1 D h (u, x k )
end for

We now state the theoretical convergence rate of Bregman gradient: as in the Euclidean case, the convergence is sublinear for convex functions and linear for relatively-strongly convex functions.

Proposition 1.4 (BGD convergence rate for convex functions). Assume that f is convex on C and L-smooth relative to a kernel function h. Let {x k } k 0 be the sequence of iterates produced by Algorithm 1 with step size 2 (0, 1/L]. Then, for every u 2 C, we have

f (x k ) f (u)  LD h (u, x 0 ) k , (1.2)
and if f is additionally µ-strongly convex relative to h,

f (x k ) f (u)  L ⇣ 1 µ L ⌘ k D h (u, x 0 ). (1.3)
See [Bauschke et al., 2017, Thm. 1] for the proof of the sublinear rate and [Lu et al., 2018, Thm. 3.1] for the linear convergence result.

Remark 1. If the relative smoothness constant L is unknown or too conservative in theory, one can also use Armijo-like line search techniques to determine the step size dynamically (see Chapter 2).

Remark 2. Let x ⇤ 2 argmin C f . In order to take u = x ⇤ in Equations (1.2) and (1.3) and obtain a bound on the suboptimality gap f (x k ) f (x ⇤ ), we need x ⇤ to belong to the domain of h. In most cases, this condition is trivially satisfied. However, it can fail if x ⇤ lies on the boundary of C and dom h is open, such as for the log kernel.

Benefits of Bregman gradient descent with relative regularity. The convergence rate of the Bregman Gradient method generalizes that of Euclidean gradient descent on L-smooth and strongly convex functions [Nesterov, 2003, Chap. 2]. Using a different geometry than the Euclidean one can benefit in two types of situations:

1. Problems with unbounded curvature: in some problems, the Hessian r 2 f (x) grows unbounded as x reaches some boundary points of C; this is the case for instance with functions involving the Kullback-Leibler divergence (Figure 1.1). Therefore, f is not Lsmooth in the standard sense, but relative smoothness can be established with a well chosen kernel that takes into account this singularity. This allows to apply the Bregman gradient method with a fixed step size and ensure global convergence.

2. Bregman preconditioning: in other situations, even if f is globally L-smooth, the performance of gradient methods can be improved by choosing a kernel that provides a better approximation of the objective, i.e., such that the gap between r 2 h and r 2 f is as small as possible. There is however a tradeoff, as a tighter approximation usually implies that the subproblem defining the Bregman iterate is harder to solve. This is analogous to the effect of preconditioning for solving linear systems.

Examples of relatively-smooth problems

We now provide a list of examples of relatively-smooth problems with various applications.

L-smooth objective and strongly convex kernel

Prior to the introduction of relative smoothness, most analyses of Bregman gradient descent for differentiable objectives made the assumption that the objective f is L-smooth and that the kernel h is -strongly convex with respect to some general norm k • k E [START_REF] Auslender | Interior Gradient and Proximal Methods for Convex and Conic Optimization[END_REF]Teboulle, 2006, Walid et al., 2015]. This is a particular case of a relatively-smooth problem, as we have for x, y 2 int C,

f (x)  f (y) + hrf (y), x yi + L 2 kx yk 2 E  f (y) + hrf (y), x yi + L D h (x, y) (1.4)

Kernels on nonnegative orthant

This category corresponds to kernels defined on the nonnegative orthant C = R d + . These kernels are mainly designed to tackle problems that have unbounded curvature near the boundary of C.

Log kernel: inverse problems with Poisson noise. Consider the problem

min x2R d + f (x) = D KL (b, Ax) (1.5)
where

D KL (x, y) = P m i=1 x (i) log ⇣ x (i) y (i) ⌘ x (i) + y (i) denotes the Kullback-Leibler divergence, A 2 R m⇥d + is a measurement matrix and b 2 R m
+ is a vector of observations. Problem (1.5) models the task of recovering an unknown signal corrupted by Poisson noise, where A encodes the measurement process. It is a fundamental problem in signal processing, with a large range of applications in astronomy, microscopy [START_REF] Bertero | Image Deblurring with Poisson Data: From Cells to Galaxies[END_REF] and medical imaging [START_REF] Ollinger | Maximum-likelihood reconstruction of transmission images in emission computed tomography via the EM algorithm[END_REF][START_REF] Ben-Tal | The Ordered Subsets Mirror Descent Optimization Method with Applications to Tomography[END_REF][START_REF] Avinash | Principles of computerized tomographic imaging[END_REF]. This is a typical situation of an objective function with unbounded curvature, as r 2 f (x) grows to infinity when A j x ! 0 for some j 2 {0 . . . m} (see Figure 1.1). However, when choosing the log-kernel

h(x) = d X i=1 log x (i) ,
we can show that f is smooth relative to h with constant L = kbk 1 [START_REF] Heinz | A Descent Lemma Beyond Lipschitz Gradient Continuity: First-Order Methods Revisited and Applications[END_REF] and thus guarantee that BGD converges with a constant step size  1/L.

Entropy kernel: nonnegative entropy regression. Consider now the problem

min x2R d + f (x) = D KL (Ax, b) where A 2 R m⇥d + , b 2 R m
+ , which differs from (1.5) as D KL is not symmetric. Then, f is smooth relative to the entropy kernel [START_REF] Heinz | A Descent Lemma Beyond Lipschitz Gradient Continuity: First-Order Methods Revisited and Applications[END_REF]. Such a function appears in entropy-regularized optimal transport, see e.g., [START_REF] Chizat | Scaling algorithms for unbalanced optimal transport problems[END_REF], [START_REF] Mishchenko | Sinkhorn Algorithm as a Special Case of Stochastic Mirror Descent[END_REF].

h(x) = d X i=1 x (i) log x (i) with constant L = max j=1...d P m i=1 A ij [

Polynomial kernels

In this situation, the objective function is usually a polynomial of degree higher than 2, whose Hessian then grows unbounded as kxk ! +1. Choosing an appropriate polynomial kernel allows to prove relative smoothness globally, and possibly to improve conditioning. Quartic polynomials. A large class of polynomial optimization problems is composed of quadratic inverse problems, that is, with objectives of the form

min x2R d f (x) = m X j=1 (hx, A j xi b j ) 2
where A j 2 R d⇥d , b j 2 R , j = 1 . . . m. Quadratic inverse problems have a broad range of applications, such as phase retrieval [START_REF] Candès | Phase retrieval via wirtinger flow: Theory and algorithms[END_REF], matrix factorization and unsupervised learning (see Chapter 2 and references therein). As shown by [START_REF] Bolte | First Order Methods Beyond Convexity and Lipschitz Gradient Continuity with Applications to Quadratic Inverse Problems[END_REF], such functions are smooth relative to the quartic kernel

h(x) = 1 4 kxk 4 + 1 2 kxk 2 .
In Chapter 2, we specialize the analysis to low-rank matrix problems and leverage this structure to provide tighter kernels. Another important application is given in [START_REF] Nesterov | Superfast second-order methods for Unconstrained Convex Optimization[END_REF][START_REF] Nesterov | Implementable tensor methods in unconstrained convex optimization[END_REF], [START_REF] Grapiglia | On inexact solution of auxiliary problems in tensor methods for convex optimization[END_REF], where the authors propose to use a refined quartic kernel for solving the subproblem that appears in regularized third-order tensor methods.

Higher-order polynomials. Kernels of higher degree have also been considered for deep linear neural networks [Mukkamala et al., 2019] and structure learning on causal models [START_REF] Romain | A Bregman Method for Structure Learning on Sparse Directed Acyclic Graphs[END_REF].

Other applications

Statistical preconditioning for distributed optimization. Consider the finite-sum minimization problem

min x2R d f (x) = 1 n n X i=1 f i (x)
where f 1 , . . . f n are differentiable functions with similar structure; e.g., the loss function on different parts of the dataset in the setting of empirical risk minimization. Then, statistical preconditioning consists in choosing the kernel h as

h(x) = f i 0 (x) + 2 kxk 2
for some i 0 2 {1 . . . n}, where > 0 is a given constant.Typically, f i 0 is the loss function on a part of the dataset of size n prec . This idea, proposed originally by [START_REF] Shamir | Communication-efficient distributed optimization using an approximate Newton-type method[END_REF], has been developped further through the lens of relative regularity by Hendrikx et al. [2020]. The latter show that, on random datasets and with high probaiblity, f is smooth and strongly convex relative to h with a relative condition number  rel that satisfies

 rel = 1 + O ✓  eucl n prec ◆
where  eucl is the standard Euclidean condition number of f . Thus, this kernel allows to improve conditioning over the Euclidean setting. There is however a tradeoff, as computing the Bregman iterates becomes harder when n prec is large. Nevertheless, this method is advantageous in the setting of distributed optimization, where the Bregman subproblem is solved efficiently by a central server, and the goal is to reduce the total number of iterations performed, since computing the gradient of f requires costly communication between different machines.

Algorithmic reinterpretations: EM and Sinkhorn. Some classical iterative methods can also be interpreted as Bregman gradient descent on relatively-smooth functions. [START_REF] Kunstner | Homeomorphic-Invariance of EM: Non-Asymptotic Convergence in KL Divergence for Exponential Families via Mirror Descent[END_REF] show that the celebrated Expectation-Maximization (EM) for maximizing the likelihood of a statistical model with latent variables is an instance of BGD. For exponential families, the objective function of EM is smooth relative to the log-partition function of the distribution. This sheds light on a new theoretical analysis of EM algorithms, as well as improved local convergence rates.

In the same spirit, [START_REF] Léger | A Gradient Descent Perspective on Sinkhorn[END_REF] studies the Sinkhorn algorithm, which is widely used for solving matrix scaling and entropy-regularized optimal transport problems. He shows that this scheme can also be seen as an instance of BGD with a relatively-smooth objective1 . The theoretical analysis establishes that Sinkhorn converges at a subblinear rate, as in Proposition 1.4.

Extensions and variants of Bregman gradient descent

We now review the algorithmic variants of BGD and extensions of relative smoothness for various classes of optimization problems.

Bregman proximal gradient. An important class of problems in signal processing and machine learning is that of composite optimization problems of the form

min x2C f (x) + g(x)
where f is smooth relative to h and g is a (possibly nonsmooth) function that encodes a penalty on the solution, such as sparsity when g is an `1 norm. This problem can be tackled with the Bregman proximal gradient method

x k+1 = argmin u2C f (x k ) + g(u) + hrf (x k ), u x k i + 1 k D h (u, x k ), (BPG)
provided that g is simple enough so that the inner subproblem can be solved easily. This method generalizes the Euclidean proximal gradient method, also known as forward-backward splitting [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF]. In their work, Bauschke et al. [2017] prove that (BPG) enjoys a sublinear convergence rate of O(1/k) when f and g are convex; see also [START_REF] Bùi | Bregman Forward-Backward Operator Splitting[END_REF] for a more general analysis.

Nonconvex problems. [START_REF] Bolte | First Order Methods Beyond Convexity and Lipschitz Gradient Continuity with Applications to Quadratic Inverse Problems[END_REF] extend the Bregman proximal gradient method to relatively-smooth nonconvex objective functions, and establish the following convergence rate of the stationarity measure:

min i=1...k D h (x i , x i 1 )  O ✓ 1 k ◆ .
Moreover, they prove convergence to a critical point when the objective satisfies the Kurdyka-Lojasiewicz property [START_REF] Bolte | The Łojasiewicz Inequality for Nonsmooth Subanalytic Functions with Applications to Subgradient Dynamical Systems[END_REF] as well as some additional local regularity properties. In subsequent work, [START_REF] Heinz | On Linear Convergence of Non-Euclidean Gradient Methods without Strong Convexity and Lipschitz Gradient Continuity[END_REF] study different conditions for guaranteeing local linear convergence to a stationary point. [START_REF] Ahookhosh | A Bregman Forward-Backward Linesearch Algorithm for Nonconvex Composite Optimization: Superlinear Convergence to Nonisolated Local Minima[END_REF] describe a line-search method for determining the step size of BGD which enjoys favorable local convergence properties.

Inertial variants of BGD have been proposed for nonconvex problems with similar convergence guarantees and improved numerical performance [Mukkamala and Ochs, 2019[START_REF] Zhang | Bregman proximal gradient algorithm with extrapolation for a class of nonconvex nonsmooth minimization problems[END_REF][START_REF] Chandra Mukkamala | Convex-Concave Backtracking for Inertial Bregman Proximal Gradient Algorithms in Non-Convex Optimization[END_REF]. In these methods, the update of BGD is augmented with an additional momentum term of the form ↵ k (x k x k 1 ), where the coefficient ↵ k is determined by a line search procedure.

Stochastic gradients. [START_REF] Hanzely | Fastest Rates for Stochastic Mirror Descent Methods[END_REF], [START_REF] Davis | Stochastic model-based minimization under high-order growth[END_REF] analyze stochastic variants of the Bregman gradient method, under additional restrictive regularity assumptions. We study stochastic BGD for general relatively-smooth problems in Chapter 6, as well as variance reduction techniques for finite-sum objectives.

Bregman block coordinate methods. In some problems, it can be advantageous to leverage a block coordinate structure of the type

x = (x 1 , . . . , x p )
where x i 2 R d i and P p i=1 d i = d, and to assume that f is smooth relative to a different kernel h i on each coordinate block. This allows the application of Bregman block coordinate methods, which alternate updates of the form

x k+1 = argmin u2R d i hr i f (x k ), u x i k i + 1 k D h i (u, x i k ),
for i 2 {1 . . . p}, where r i f denotes the partial gradient of f on the block i. This approach has been considered in several works [START_REF] Teboulle | Novel proximal gradient methods for nonnegative matrix factorization with sparsity constraints[END_REF][START_REF] Ahookhosh | Multi-block Bregman proximal alternating linearized minimization and its application to orthogonal nonnegative matrix factorization[END_REF], Hendrikx et al., Gao et al., 2020]. It is particularly efficient in situations where Bregman block coordinate updates are much easier to compute, or if restricting on separate blocks allows to prove tighter relative smoothness inequalities than on the full coordinates, such as for certain matrix factorization problems.

Relative regularity beyond differentiable functions. The idea of choosing a Bregman geometry that is adapted to the singularities of the objective through relative regularity has also been extended to other settings than that of differentiable functions, such as nonsmooth convex functions [START_REF] Lu | Relative Continuity" for Non-Lipschitz Nonsmooth Convex Optimization Using Stochastic (or Deterministic) Mirror Descent[END_REF], online learning [START_REF] Kimon Antonakopoulos | Online and Stochastic Optimization beyond Lipschitz Continuity: A Riemannian Approach[END_REF], variational inequalities and monotone operators [START_REF] Antonakopoulos | An Adaptive Mirror-Prox Algorithm for Variational Inequalities with Singular Operators[END_REF][START_REF] Stonyakin | Inexact relative smoothness and strong convexity for optimization and variational inequalities by inexact model[END_REF][START_REF] Michael | Relative Lipschitzness in Extragradient Methods and a Direct Recipe for Acceleration[END_REF].

Variants of relative smoothness. [START_REF] Gutman | The condition number of a function relative to a set[END_REF] extend the notion to other type of distances than Bregman divergences, and to also take into account the geometry of the optimization set. [START_REF] Maddison | Dual Space Preconditioning for Gradient Descent[END_REF] propose a dual variant of relative smoothness, which is invariant to translation (unlike the standard version), with application to `p norm regression.

The issue of acceleration and non-homogeneity

Algorithm 2 Accelerated Bregman gradient method in generic form Input:

initial point x 0 2 int C Set z 0 = x 0 .
for k = 0,1,. . . do Determine interpolation coefficient ↵ k and step size k ,

y k = (1 ↵ k )x k + ↵ k z k , z k+1 = argmin {hrf (y k ), u y k i + 1 k D h (u, z k ) | u 2 R n }, x k+1 = (1 ↵ k )x k + ↵ k z k+1 ,

end for

Acceleration. For convex relatively-smooth functions, Proposition 1.4 states that BGD converges at a rate of O(1/k) in function values, which is a rather slow rate, both in theory and in practice. An important question is whether there exists a first-order Bregman algorithm that achieves a better convergence rate on the same class of problems.

In the Euclidean setting, the answer is positive and was provided in the seminal work of [START_REF] Nesterov | A Method of Solving A Convex Programming Problem With Convergence rate O(1/kˆ2)[END_REF] 2 . He showed that the accelerated gradient method, obtained by adding a simple and well-chosen linear extrapolation step to gradient descent, achieves a rate of O(1/k 2 ) on L-smooth convex functions, and that this worst-case rate is optimal among all first-order Euclidean methods [START_REF] Nesterov | Introductory Lectures on Convex Optimization: A Basic Course[END_REF]. Significant improvement in convergence speed is also observed in practice. This so-called acceleration technique has been adapted to improve the efficiency of several other algorithms in continuous optimization; see the monograph by D'Aspremont et al.

[2021] for a recent survey.

Naturally, the acceleration technique has been applied to Bregman gradient methods; most instances of accelerated BGD follow the generic form described in Algorithm 2. The choice of the interpolating coefficients ↵ k and the step size parameters k are crucial for guaranteeing fast convergence. Several choices have been proposed in different settings:

• Smooth objective and strongly convex kernel: the first instance of accelerated BGD was studied by [START_REF] Auslender | Interior Gradient and Proximal Methods for Convex and Conic Optimization[END_REF], under the assumption that f is

L k•k E - smooth and h is k•k E -strongly convex with respect to some norm k • k E .
In this setting, they show that Algorithm 2 with the same choice of coefficients as Nesterov's accelerated gradient method (that is, k / k and ↵ k / 1/k) yields the improved convergence rate of

O(1/k 2 ). However, the algorithm's performance is dependent on k•k E /L k•k E ,
which can be very small or even equal to 0 in the context of problems with unbounded curvature.

• Relatively-smooth objective: in this setting, the first tentative of acceleration was proposed by [START_REF] Hanzely | Accelerated Bregman Proximal Gradient Methods for Relatively Smooth Convex Optimization[END_REF]. Their analysis rely on the triangle scaling exponent of the Bregman divergence, that is, the largest exponent such that

D h (1 ✓) x + ✓z, (1 ✓) x + ✓y  ✓ D h (x, y), 8✓ 2 [0, 1], (TSE)
for every x, y, z 2 int C. Under such property, the authors show that an appropriate choice of coefficients in Algorithm 2 allows to achieve a O(1/k ) convergence rate. However, besides quadratic kernels for which = 2 and the entropy where = 1, the inequality (TSE) usually does not hold globally. It is nevertheless valid on small bounded subsets under sufficient regularity of h, as the Bregman divergence locally behaves as a quadratic function. Thus, the authors propose adaptive variants of accelerated BGD to take advantage of this property. Qualitative arguments can then be used to show that the improved rate of O(1/k 2 ) is reached asymptotically, and numerical experiments demonstrate improved performance over standard BGD.

Similar results are obtained by Hendrikx et al. [2020] with application to statistical preconditioning for distributed optimization. They propose an adaptive variant of Algorithm 2 which also leverages relative strong convexity. The analysis shows that an improved linear convergence rate of O ⇣ (1

q µ L ) k ⌘ is reached asymptotically.
Although numerical experiments are encouraging, there is still a gap in our theoretical understanding of accelerated Bregman gradient methods in the relatively-smooth setting. In Chapter 3, we show that in the general case, the O(1/k) rate of BGD is optimal among Bregman first-order methods and therefore that global acceleration is out of reach without additional regularity assumptions.

Non-homogeneity and other difficulties. The main difficulty when trying to transpose convergence results of Euclidean methods to the Bregman setting arises from the lack of homogeneity and translational invariance of the Bregman divergence. Indeed, for non-quadratic kernels, the quantity

D h (x + v, x)
depends on the point x, and is not proportional to 2 , whereas it reduces to 2 kvk 2 in the Euclidean setting. This issue is met when analyzing methods that combine gradients taken at different points and those that rely on step size scaling. This is the case for accelerated methods as well as variance reduction techniques for stochastic optimization (see Chapter 6).

Although it is a different setting, the difficulty of accelerating gradient methods also arises in Riemannian optimization for similar reasons [START_REF] Zhang | An Estimate Sequence for Geodesically Convex Optimization[END_REF]; see [START_REF] Hamilton | No-go Theorem for Acceleration in the Hyperbolic Plane[END_REF] for a negative result on the hyperbolic plane and an interesting intuitive explanation.

Besides acceleration, there is also difficulty in proving the convergence of the iterates of Bregman methods towards the optimal point. [START_REF] Bolte | Curiosities and counterexamples in smooth convex optimization[END_REF] provide an an instance of BGD where the iterate sequence {x k } k 0 does not converge and cycles indefinitely. In their example, the objective function is linear; the issue is created uniquely by the choice of a highly pathological Bregman kernel.

A brief history of Bregman methods before relative smoothness

Optimization methods with Bregman divergences have also been considered for a long time outside of the relatively-smooth setting. In this section, we provide a historical overview.

Bregman projections. Given a kernel h, the Bregman projection of a point x on a convex set

A ⇢ R d is proj h A (x) = argmin u2A D h (u, x). (1.6)
This tool was first introduced by [START_REF] Bregman | The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming[END_REF], along with the corresponding divergence, in order to find a point in the intersection of N convex sets A 1 . . . A N by alternating such projections. The main interest of choosing a specific kernel h is that it induces a bias that drives the algorithm towards the part of the intersection set that minimizes h3 . More precisely, with a proper initialization, the method of alternating Bregman projections converges to the solution of the problem

min x2dom h h(x) s.t. x 2 A 1 \ • • • \ A N .
See also [START_REF] Heinz | Legendre functions and the method of random Bregman projections[END_REF] for a more detailed analysis of Bregman projections.

Bregman proximal methods. For a convex proper function g : R d ! R [ {+1} and a kernel h, the Bregman proximal map of g is defined as

prox h g (x) = argmin u2R d g(u) + D h (u, x),
which generalizes the Moreau proximal map [START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF] to the non-Euclidean setting. Note also that taking g as the indicator function of a convex set recovers the Bregman projection (1.6). This operator has been introduced by [START_REF] Censor | Proximal Minimization Algorithm with D-functions[END_REF], who study the Bregman proximal point algorithm

x k+1 = argmin u2R d g(u) + 1 k D h (u, x k ) (BPP)
for minimizing the (possibly nonsmooth) function g on C. Simultaneously, [START_REF] Teboulle | Entropic Proximal Mappings with Applications to Nonlinear Programming[END_REF] studied the proximal map with Bregman divergence and also an alternative class of divergences called -divergences. He showed that these maps enjoyed a smoothing effect, as does the Moreau proximal map in the Euclidean case. The (BPP) scheme was then extended by [START_REF] Eckstein | Nonlinear Proximal Point Algorithms Using Bregman Functions, with Applications to Convex Programming[END_REF] to the more general setting of monotone operators, which allows to build nonquadratic versions of augmented Lagrangian methods. One of the main advantages of using a Bregman divergence adapted to C in proximal methods is that the kernel acts as a barrier, forcing the iterates to belong to the interior of C, thus ensuring favorable regularity properties. For instance, augmented Lagrangian methods with a well-chosen Bregman kernel have a twice continuously differentiable Lagrangian, unlike the Euclidean counterpart, allowing the use of efficient second-order methods for computing the iterates. Further analysis is provided in [START_REF] Chen | Convergence Analysis of a Proximal-Like Minimization Algorithm Using Bregman Functions[END_REF], [START_REF] Iusem | Entropy-Like Proximal Methods in Convex Programming[END_REF].

Nonsmooth kernels for ill-posed inverse problems. Although Bregman methods are generally used with differentiable kernels, nonsmooth convex kernels can also be used, provided that the subgradient defining D h is well-chosen. This extension was first studied by [START_REF] Krzysztof | Proximal minimization methods with generalized Bregman functions[END_REF] for the Bregman proximal point method.

In this situation, the Bregman kernel is chosen so as to exploit the idea of implicit bias: in the same way as for Bregman projections (see the paragraph above), the iterates of (BPP) converge towards the part of argmin f that also minimizes h. This can be useful for solving ill-posed inverse problems: for instance, taking h to be a `1 norm allows to induce sparsity in the iterates of the algorithm. This idea has had great success in compressed sensing and image denoising, where the Bregman proximal point scheme is also known as Bregman iteration, or Bregman iterative regularization [START_REF] Osher | An iterative regularization method for total variation-based image restoration[END_REF][START_REF] Yin | Bregman iterative algorithms for l1-minimization with applications to compressed sensing[END_REF][START_REF] Goldstein | The split Bregman method for L1-regularized problems[END_REF], see also [START_REF] Benning | Bregman Methods for Large-Scale Optimisation with Applications in Imaging[END_REF] for a recent survey.

Mirror descent subgradient method. The mirror descent algorithm, designed for minimizing a possibly nonsmooth convex function f over a set C, writes

x k+1 = argmin u2R d hf 0 (x k ), u x k i + 1 k D h (u, x k ), (MD)
where f 0 (x k ) 2 @f (x k ) is any subgradient of f at x k , and the sequence of step sizes { k } k 0 is decreasing to 0. Mirror descent was originally proposed by [START_REF] Nemirovski | Efficient methods for large-scale convex optimization problems[END_REF], [START_REF] Nemirovski | Problem Complexity and Method Efficiency in Optimization[END_REF] in the form of

x k+1 = rh ⇤ ⇥ rh(x k ) k f 0 (x k ) ⇤ (MD')
and applied by [START_REF] Ben-Tal | The Ordered Subsets Mirror Descent Optimization Method with Applications to Tomography[END_REF] to Poisson inverse problems (which were then classified as nonsmooth problems, because of the unbounded curvature). Later, [START_REF] Beck | Mirror Descent and Nonlinear Projected Subgradient Methods for Convex Optimization[END_REF] provided a simplified analysis and showed that, under appropriate assumptions on h, the update (MD') can be rewritten as the subgradient scheme with Bregman divergence desribed in (MD).

A central element in the analysis of Mirror Descent is the choice of a (possibly non-Euclidean) norm k • k E for the space E of iterates x k , along with the corresponding dual norm k • k E ⇤ for the space E ⇤ of gradients rh(x k ), f 0 (x k ). The kernel h is assumed to be strongly convex with respect to k • k E , and the subgradients of f are assumed to be bounded with respect to k

• k E ⇤ by a constant L k•k E on C.
The efficiency estimate of (MD) is proportional to

L k•k E (D h (x ⇤ , x 0 )) 1/2 .
Therefore, the goal is to choose the geometry that minimizes this estimate on the set C. The most classical example is the unit simplex C = {x 2 R d + :

P i x i = 1}
, for which choosing the `1 norm along with the entropy kernel allows to improve the efficiency by a factor of (d/ log d) 1/2 over the Euclidean geometry. This is a considerable theoretical and practical gain for largescale problems. As an additional advantage, the Bregman projection on the simplex can be computed in closed form, unlike the Euclidean one. A detailed analysis of the computational complexity of mirror descent in various situations is provided in the survey by [START_REF] Juditsky | First Order Methods for Nonsmooth Convex Large-Scale Optimization , I : General Purpose Methods[END_REF].

Beyond nonsmooth convex minimization, the mirror descent scheme has been studied in the context of stochastic convex minimization and saddle-point problems [START_REF] Nemirovski | Robust stochastic approximation approach to stochastic programming[END_REF][START_REF] Duchi | Composite objective mirror descent[END_REF] as well as online learning [START_REF] Hazan | The Convex Optimization Approach to Regret Minimization[END_REF][START_REF] Bubeck | Introduction to online optimization[END_REF]. As the simplex models discrete probability measures, entropic mirror descent has had great success in online learning, where it is also known as multiplicative weight updates. More recently, mirror descent schemes have also been used for solving reinforcement learning problems [START_REF] Tomar | Mirror Descent Policy Optimization[END_REF][START_REF] Lan | Policy Mirror Descent for Reinforcement Learning: Linear Convergence, New Sampling Complexity, and Generalized Problem Classes[END_REF].

From mirror descent to Bregman gradient method and relative smoothness. The mirror descent method (MD) has a similar form to the Bregman gradient scheme. However, mirror descent has mainly been applied in the context of nonsmooth, stochastic and online convex optimization. In that setting, the kernel h is usually chosen accordingly to the constraint set C and the step size sequence { k } k 0 needs to be decreasing towards 0 to account for nonsmoothness.

The idea of adapting the kernel to the curvature of the objective function through the simple condition r 2 f (x) Lr 2 h(x) has opened new possibilities for applying Bregman methods to differentiable objectives. Interestingly, the relative smoothness assumption was suggested for the first time by [START_REF] Birnbaum | Distributed algorithms via gradient descent for fisher markets[END_REF] in the context of algorithmic game theory, but remained unnoticed by the optimization community. Later, [START_REF] Heinz | A Descent Lemma Beyond Lipschitz Gradient Continuity: First-Order Methods Revisited and Applications[END_REF] discovered the same concept independently and popularized the subject among researchers in this field.

Chapter 2

Quartic First-Order Methods for Low-Rank Minimization

Introduction

Consider a low-rank semidefinite program, written

min F (Y ) subject to Y ⌫ 0, rank(Y )  r (SDP-r)
in the variable Y 2 R n⇥n , where F : R n⇥n ! R is typically a differentiable convex function with Lipschitz continuous gradient and r  n is the target rank. Fundamental applications of (SDP-r) arise in various areas including matrix completion [START_REF] Candès | Exact Matrix Completion Via Convex Optimization[END_REF][START_REF] Cai | A Singular Value Tresholding Algorithm for Matrix Completion[END_REF][START_REF] Jain | Low-rank Matrix Completion using Alternating Minimization[END_REF], matrix sensing [START_REF] Recht | Guaranteed Minimum-Rank Solutions of Linear Matrix Equations via Nuclear Norm Minimization[END_REF], Euclidean matrix completion [Mishra et al., 2011, Fang and[START_REF] Haw | Euclidean distance matrix completion problems[END_REF], phase retrieval [START_REF] Candès | Phase retrieval via wirtinger flow: Theory and algorithms[END_REF], robust principal component analysis [START_REF] Chen | Fast low-rank estimation by projected gradient descent: General statistical and algorithmic guarantees[END_REF], to name a few.

A popular approach to solving (SDP-r), known as the Burer-Monteiro formulation [START_REF] Burer | Local Minima and Convergence in Low-Rank Semidefinite Programming[END_REF], consists in representing Y as Y = XX T to solve min (X) , F (XX T ) + g(X)

(P) in the variable X 2 R n⇥r , where g is a lower semicontinuous simple function (in a sense that will be made clear later) that enforces additional constraints or penalties on the factor X. We will write f : R n⇥r ! R the factorized function with f (X) = F (XX T ). This reformulation has been thoroughly studied for various low rank minimization problems [START_REF] Mishra | Fixed-rank matrix factorizations and Riemannian low-rank optimization[END_REF][START_REF] Tu | Low-rank Solutions of Linear Matrix Equations via Procrustes Flow[END_REF][START_REF] Candès | Phase retrieval via wirtinger flow: Theory and algorithms[END_REF][START_REF] Zheng | A Convergent Gradient Descent Algorithm for Rank Minimization and Semidefinite Programming from Random Linear Measurements[END_REF][START_REF] Chen | Fast low-rank estimation by projected gradient descent: General statistical and algorithmic guarantees[END_REF][START_REF] Bhojanapalli | Dropping Convexity for Faster Semidefinite Optimization[END_REF][START_REF] Zhao | A Nonconvex Optimization Framework for Low Rank Matrix Estimation[END_REF][START_REF] Sun | Guaranteed Matrix Completion via Nonconvex Factorization[END_REF][START_REF] Park | Finding low-rank solutions to matrix problems , efficiently and provably[END_REF]. It has several benefits: first, the rank constraint is directly enforced by the factorization XX T . Second, it allows adding penalty functions g to impose additional structural properties on the factor, like nonnegativity or sparsity. More importantly, the problem size is reduced from n 2 to nr variables, which makes it far more efficient on large-scale problems, since r is usually much smaller than n.

This formulation comes however with an important drawback: the objective function of (P) becomes nonconvex, even when F is originally convex. Therefore, local optimization methods can generally only hope to find a stationary point, or at best a local minimum. Nevertheless, recent work shows convergence towards a global optimum for a close enough initialization [START_REF] Tu | Low-rank Solutions of Linear Matrix Equations via Procrustes Flow[END_REF][START_REF] Bhojanapalli | Dropping Convexity for Faster Semidefinite Optimization[END_REF][START_REF] Park | Finding low-rank solutions to matrix problems , efficiently and provably[END_REF], or under additional statistical assumptions about the problem [START_REF] Chen | Fast low-rank estimation by projected gradient descent: General statistical and algorithmic guarantees[END_REF][START_REF] Zheng | A Convergent Gradient Descent Algorithm for Rank Minimization and Semidefinite Programming from Random Linear Measurements[END_REF][START_REF] Ge | Matrix Completion has No Spurious Local Minimum[END_REF]. Although these global optimality results often impose restrictive assumptions that may not be satisfied in practice, they help to explain why the use of local algorithms for solving (P) usually leads to satisfactory solutions in practice.

The most commonly used algorithm to solve these problem formulations is some variant of the proximal gradient method. However, a critical issue with gradient schemes is the choice of step sizes, which significantly impacts performance. This step size choice is closely related to the smoothness of the objective. In particular, when it has a L-Lipschitz continuous gradient with respect to the Euclidean norm, standard gradient methods can be applied with a step size lying in ]0, 1/L]. This smoothness assumption is used in the broad majority of theoretical analyses of gradient algorithms, yet there are many cases where it is not satisfied [START_REF] Heinz | A Descent Lemma Beyond Lipschitz Gradient Continuity: First-Order Methods Revisited and Applications[END_REF][START_REF] Bolte | First Order Methods Beyond Convexity and Lipschitz Gradient Continuity with Applications to Quadratic Inverse Problems[END_REF]. In particular, it does not hold for the general Burer-Monteiro low-rank problem, as we will show in what follows.

Of course, there is a way to circumvent this issue in classical Euclidean methods, by using an Armijo line search [START_REF] Lin | Projected Gradient Methods for Nonnegative Matrix Factorization[END_REF]. However, in some cases, this naive line search strategy generates very small step sizes which in turn involve costly subroutines. Other approaches impose a step size that is only proven to be valid in a small neighborhood of the optimum [START_REF] Bhojanapalli | Dropping Convexity for Faster Semidefinite Optimization[END_REF][START_REF] Park | Finding low-rank solutions to matrix problems , efficiently and provably[END_REF].

Bregman gradient methods. We adopt an original approach based on a recent line of work on non-Euclidean gradient methods [START_REF] Heinz | A Descent Lemma Beyond Lipschitz Gradient Continuity: First-Order Methods Revisited and Applications[END_REF][START_REF] Bolte | First Order Methods Beyond Convexity and Lipschitz Gradient Continuity with Applications to Quadratic Inverse Problems[END_REF] and subsequent work [START_REF] Lu | Relatively-Smooth Convex Optimization by First-Order Methods, and Applications[END_REF]. Unlike standard gradient descent that uses the uniform Euclidean geometry, the NoLips method, also known as Bregman proximal gradient/Mirror descent, uses the Bregman divergence induced by a well-chosen convex kernel function. This allows the algorithm to take gradient steps that are more adapted to the geometry of the problem, advancing faster in directions where the gradient of the objective changes slowly, thus improving convergence speed. The kernel function is chosen so that the objective function satisfies a compatibility condition called relative smoothness [START_REF] Heinz | A Descent Lemma Beyond Lipschitz Gradient Continuity: First-Order Methods Revisited and Applications[END_REF][START_REF] Lu | Relatively-Smooth Convex Optimization by First-Order Methods, and Applications[END_REF], which is a generalization of the usual smoothness assumption mentioned earlier.

In our setting, the objective has a quartic growth, hence choosing the geometry induced by a quartic polynomial will prove to be efficient.

Contributions.

In this work, we focus on deriving efficient algorithms to find stationary points of nonconvex low-rank problems. Our main contribution is to identify favorable non-Euclidean geometries for these problems, induced by well-chosen quartic kernels. We first study a simple quartic norm kernel that is compatible with various regularization terms. We then introduce a novel family of quartic kernels that we call Gram kernels, which can be applied to unregularized problems. They provide richer geometries which greatly improve convergence speed with little impact on the iteration complexity. We also extend the NoLips/Bregman Gradient scheme to Dyn-NoLips, allowing for adaptive step size strategies.

To highlight the benefits of our approach, we study applications to symmetric nonnegative matrix factorization and Euclidean distance matrix completion and show competitive numerical performance compared to specialized algorithms for these problems.

Notation. For a square matrix M , we denote its trace Tr M = P n i=1 M ii . For two matrices X and Y of same size, we denote the standard Euclidean inner product and norm by hX, Y i = Tr(X T Y ) and kXk = p Tr(X T X). For a function f : R n⇥r ! R, we denote by rF (X) its gradient matrix rF (X) ij = @F (X) @x ij and by r 2 F (X)[U, V ] the second derivative at X in the directions U, V 2 R n⇥r . I r denotes the identity matrix of size r ⇥ r. For two square matrices X, Y , we write X Y if the matrix Y X is positive semidefinite. We write kAk op for the operator norm of a linear application A.

Quartic geometries for low-rank minimization

Problem Setup

We first state our standing assumptions for Problem (P).

Assumption 2.1. (a) F : R n⇥n ! R is a twice continuously differentiable function which is µ F -strongly convex and L F -smooth, i.e.,

hrF (X) rF (Y ), X Y i µ F kX Y k 2 , krF (X) rF (Y )k  L F kX Y k 8X, Y 2 R n⇥n , (b) g : R n⇥n ! R [ {+1} is a closed convex proper function, (c) min R n⇥r > 1.
Our analysis will involve the following lemma.

Lemma 2.1. Let F : R n⇥n ! R be a twice differentiable µ F -strongly convex and L F -smooth function. Then, the function

G := F µ F 2 k • k 2 is convex and (L F µ F )-smooth. Proof.
It suffices to use the second-order characterization [START_REF] Nesterov | Introductory Lectures on Convex Optimization: A Basic Course[END_REF] and notice that, for

Y, U 2 R n⇥n , we have r 2 G(Y )[U, U ] = r 2 F (Y )[U, U ] µ F kU k 2 and hence µ F kU k 2  r 2 F (Y )[U, U ]  L F kU k 2 =) 0  r 2 G(Y )[U, U ]  (L F µ F )kU k 2 .
We write f : R n⇥r ! R the factorized function defined by f (X) := F (XX T ).

Therefore Problem (P) can be written as a standard composite optimization problem min X2R n⇥r (X) = f (X) + g(X).

(2.1)

Relative Smoothness and the Bregman Proximal Gradient Map

In this section, we recall the framework of [START_REF] Heinz | A Descent Lemma Beyond Lipschitz Gradient Continuity: First-Order Methods Revisited and Applications[END_REF], [START_REF] Bolte | First Order Methods Beyond Convexity and Lipschitz Gradient Continuity with Applications to Quadratic Inverse Problems[END_REF] to derive non-Euclidean gradient methods for solving composite problems of the form (2.1).

The first essential step is the choice of a distance kernel. In our context we choose a differentiable strictly convex function h : R n⇥r ! R, with dom h = R n⇥r (although more general distance kernels can be used). The distance kernel h induces in turn a Bregman divergence

D h (X, Y ) = h(X) h(Y ) hrh(Y ), X Y i.
Note that D h is not a proper distance, as it is not symmetric in general. However D h enjoys a distance-like separation property: D h (X, X) = 0 and D h (X, Y ) > 0 for X 6 = Y . The choice of a distance kernel suited to the function f is guided by the following relative smoothness condition, also called generalized Lipschitz property.

Definition 2.2 (Relative smoothness). We say that a differentiable function f on R n⇥r is Lsmooth relative to the kernel function h if there exists L > 0 such that for every X, Y 2 R n⇥r ,

f (X)  f (Y ) + hrf (Y ), X Y i + L D h (X, Y ). (RelSmooth)
For twice differentiable functions, relative smoothness has an elementary characterization: f is L-smooth relative to h if and only if

r 2 f (X)[U, U ]  L r 2 h(X)[U, U ], 8X, U 2 R n⇥r .
(2.2)

Notice that if h(X) = 1 2 kXk 2 , then D h (X, Y ) = 1 2 kX Y k 2 and we recover the standard Euclidean descent lemma that would be implied by Lipschitz continuity of the gradient of f . Bregman proximal gradient map. Now that we are equipped with a non-Euclidean geometry generated by h, we define the Bregman proximal gradient map with step size as follows.

T (X) = argmin U 2R n⇥r ⇢ g(U ) + f (X) + hrf (X), U Xi + 1 D h (U, X) , (2.3)
which consists in minimizing a surrogate for where f has been replaced by the upper approximation given by (RelSmooth) and the nonsmooth part g is kept intact, generalizing thus the approach used in the proximal gradient method. The relative smoothness condition ensures that this operation decreases the objective when 2]0, 1/L]. This map is the basic brick for non-Euclidean methods à la Bregman. The simplest method is the Bregman proximal gradient method, also known as NoLips [START_REF] Heinz | A Descent Lemma Beyond Lipschitz Gradient Continuity: First-Order Methods Revisited and Applications[END_REF] and its extension Dyn-NoLips (Algorithm 3), which simply amounts to iterating X k+1 = T k (X k ), but other possibilities exist using momentum ideas [START_REF] Auslender | Interior Gradient and Proximal Methods for Convex and Conic Optimization[END_REF][START_REF] Hanzely | Accelerated Bregman Proximal Gradient Methods for Relatively Smooth Convex Optimization[END_REF][START_REF] Chandra Mukkamala | Convex-Concave Backtracking for Inertial Bregman Proximal Gradient Algorithms in Non-Convex Optimization[END_REF].

The Quartic Geometry

In order to provide some insight into the quartic geometry of our problem, let us consider the example where F is a quadratic function, i.e.,

F (Y ) = 1 2 hAY, Y i + hB, Y i 8Y 2 R n⇥n , (2.4)
where B 2 R n⇥n and A : R n⇥n ! R n⇥n is some linear map. Then, f writes

f (X) = F (XX T ) = 1 2 hA(XX T ), XX T i + hBX, Xi 8X 2 R n⇥r .
Clearly, f is a quartic function and its gradient is not Lipschitz continuous on R n⇥r , as the Hessian "grows" to infinity when kXk ! 1. In other words, (RelSmooth) does not hold with the Euclidean kernel h = 1 2 k • k 2 . We now show that relative smoothness holds with a family of well-chosen quartic kernels, which are more adapted to the geometry of f .

The Quartic Norm Kernel

We begin with a simple quartic kernel, which depends solely on the Frobenius norm of X.

Define the norm kernel h

N as h N (X) = ↵ 4 kXk 4 + 2 kXk 2 8 X 2 R n⇥r ,
where ↵, > 0 are fixed parameters. Note that this kernel is not new by itself, as it has been already studied in [START_REF] Bolte | First Order Methods Beyond Convexity and Lipschitz Gradient Continuity with Applications to Quadratic Inverse Problems[END_REF] for vectors in R n . Our first contribution is to show that it is adapted to every function of our class of problems.

Proposition 2.3 (Norm kernel). The function f is 1-smooth relative to the norm kernel h N for ↵ 6L F and 2krF (0)k.

Proof. As F is twice differentiable, then so is f and we can use the Hessian characterization (2.2). For X, U 2 R n⇥r , the second derivative of h N is written

r 2 h N (X)[U, U ] = ↵ kXk 2 kU k 2 + 2hX, U i 2 + kU k 2 ↵kXk 2 kU k 2 + kU k 2 .
(2.5)

On the other hand, the second derivative of f is

r 2 f (X)[U, U ] = r 2 F (XX T )[UX T + XU T , UX T + XU T ] + 2hrF (XX T ), UU T i. (2.6)
Since F has a Lipschitz continuous gradient, the standard second derivative inequality yields

r 2 F (XX T )[UX T + XU T , UX T + XU T ]  L F kUX T + XU T k 2 .
Now, the second term can be bounded by using the triangle inequality, the Cauchy-Schwarz inequality and the gradient Lipschitz property, to get

hrF (XX T ), UU T i = hrF (0), UU T i + hrF (XX T ) rF (0), UU T i  krF (0)k kU k 2 + krF (XX T ) rF (0)k kUU T k  ⇣ krF (0)k + L F kXX T k ⌘ kU k 2 hence r 2 f (X)[U, U ]  L F kUX T + XU T k 2 + 2 L F kXX T k + krF (0)k kU k 2  2L F kUX T k 2 + kXU T k 2 + 2 L F kXX T k + krF (0)k kU k 2  6L F kXk 2 kU k 2 + 2krF (0)kkU k 2  ↵kXk 2 kU k 2 + kU k 2 (2.7)
where we used the submultiplicative property of the Frobenius norm, and our choice of parameters ↵, . Combining (2.5) and (2.7) gives that et al., 2017].

r 2 f (X)[U, U ]  r 2 h N (X)[U, U ] for all X, U 2 R n⇥r , hence that f is 1-smooth relative to h [Bauschke
The Bregman proximal gradient map (2.3) associated with the kernel h N can be computed easily in closed form. We give its expression in the unconstrained case [START_REF] Bolte | First Order Methods Beyond Convexity and Lipschitz Gradient Continuity with Applications to Quadratic Inverse Problems[END_REF].

Proposition 2.4 (Bregman gradient map for h N , unconstrained case). Assume that there is no penalty term, i.e., that g ⌘ 0. The Bregman gradient map of the norm kernel h N with step size > 0 is given by Compared to a standard gradient iteration, the additional operations are elementary and have a minimal impact on the arithmetic complexity.

T (X) = 1 ⌧ (↵kU k 2 ) U where U = rh N (X) rf (X) = (↵kXk 2 + )X rf (
Constraints and regularization terms. Following the ideas in [START_REF] Bolte | First Order Methods Beyond Convexity and Lipschitz Gradient Continuity with Applications to Quadratic Inverse Problems[END_REF], the Bregman proximal gradient map of h N can also be easily computed in closed form when g is the `1 norm or the `0 pseudonorm. As we will show in Section 2.4.1, this is also elementary when g is the indicator function of the nonnegative orthant.

A More Refined Kernel for Unregularized Problems: the Gram Kernel

While the kernel h N is simple and compatible with many penalties g, a better kernel can be derived for unconstrained instances by considering a richer geometry involving the Gram matrix. Define the Gram kernel as

h G (X) = ↵ 4 kXk 4 + 4 kX T Xk 2 + 2 kXk 2 8 X 2 R n⇥r ,
where ↵, 0, > 0 are given parameters. The Gram kernel is more refined than the previous norm kernel since it incorporates some nonisotropic information with the kX T Xk 2 term. To show where this term stems from, observe that following Lemma 2.1, F can be decomposed as

F = µ F 2 k • k 2 + F where F is (L F µ F )-smooth. Hence f writes f (X) = F (XX T ) = µ F 2 kXX T k 2 + F (XX T ).
Since kXX T k 2 = kX T Xk 2 , the first term can be directly incorporated into the kernel, which allows to prove a tighter relative smoothness inequality.

Proposition 2.5 (Gram kernel). f is 1-smooth relative to the Gram kernel h G a choice of constants satisfying ↵ 2(L F µ F ), 2L F and 2krF (0)k.

Proof. This amounts to refine the analysis of the proof of Proposition 2.3. Let X, U 2 R n⇥r .

The second derivative of h G at X in the direction U writes

r 2 h G (X)[U, U ] = ↵ kXk 2 kU k 2 + 2hX, U i 2 + ✓ 1 2 kUX T + XU T k 2 + kU T Xk 2 ◆ + kU k 2 ↵kXk 2 kU k 2 + ✓ 1 2 kUX T + XU T k 2 + kU T Xk 2 ◆ + kU k 2 .
On the other hand, following (2.6) the second derivative of f satisfies

r 2 f (X)[U, U ]  L F kUX T + XU T k 2 + 2hrF (XX T ), UU T i.
To bound the second term, we use Lemma 2.1 which states that the function G(Y ) := F (Y ) µ F kY k 2 /2 is convex and smooth with constant L F µ F . Using the gradient Lipschitz property of G yields

hrF (XX T ), UU T i = hrF (0), UU T i + µ F hXX T , UU T i + hrF (XX T ) rF (0) µ F (XX T 0), UU T i = hrF (0), UU T i + µ F kU T Xk 2 + hrG(XX T ) rG(0), UU T i  krF (0)k kU k 2 + µ F kU T Xk 2 + (L F µ F )kXX T k kUU T k  krF (0)k kU k 2 + L F kU T Xk 2 + (L F µ F )kXk 2 kU k 2 ,
using that µ F  L F , and so we have

r 2 f (X)[U, U ]  2(L F µ F )kXk 2 kU 2 k + L F kUX T + XU T k 2 + 2L F kU T Xk 2 + 2krF (0)k kU k 2  ↵kXk 2 kU k 2 + 2 kUX T + XU T k 2 + kU T Xk 2 + kU k 2  r 2 h G (X)[U, U ]
which shows that, for the prescribed choice of ↵, , , the function f is 1-smooth relative to h G .

Approximation quality for well-conditioned F . Let us illustrate the advantage of the Gram kernel when F is well-conditioned. For simplicity, assume here that F is a quadratic function, as in (2.4), i.e., F (Y ) = 1 2 hAY, Y i + hB, Y i where A is a positive semidefinite linear operator on R n⇥r , and hence f has a quartic and a quadratic term

f (X) = 1 2 hA(XX T ), XX T i + hB, XX T i.
The gap between f and h G with the choice of coefficients prescribed by Proposition 2.5 writes, for

X 2 R n⇥r , h G (X) f (X) = (L F µ F ) 2 kXk 4 + L F 2 kX T Xk 2 + krF (0)k kXk 2 1 2 hA(XX T ), XX T i hBX, Xi = (L F µ F ) 2 kXk 4 + 1 2 h(L F I A)(XX T ), XX T i | {z } d 4 (X) + h(krF (0)kI B)X, Xi | {z } d 2 (X)
where we separated the gap into a quartic term d 4 and a quadratic term d 2 . It can be seen from ( 2.2) that the quality of approximation of the kernel is given by the difference of the Hessians.

Focusing on the quartic part, the Hessian difference is

r 2 d 4 (X)[U, U ] = 2(L F µ F ) kXk 2 kU k 2 + 2hX, U i 2 + 2h(L F I A)(XX T ), UU T i + h(L F I A)(UX T + XU T ), UX T + XU T i  6(L F µ F )kXk 2 kU k 2 + kL F I Ak op 2kXX T k kUU T k + kUX T + XU T k 2
for X, U 2 R n⇥r . Recalling that F is L F -smooth and µ F -strongly convex, we have that

kL F I Ak op  (L F µ F ), therefore r 2 d 4 (X)[U, U ]  (L F µ F ) 6kXk 2 kU k 2 + 2kXX T k kUU T k + kUX T + XU T k 2  12L F (1 µ F L F )kXk 2 kU k 2
which shows that the quality of approximation of the quartic part of f by the Gram kernel depends on the condition number  F := L F /µ F of F . Note that one could actually refine the analysis by replacing  F with the condition number of F restricted to the set of matrices of rank at most 2r, which can be much smaller. This is the case when the linear map A satisfies the restricted isometry property (RIP), which occurs with high probability in matrix sensing applications with a sufficiently large number n of samples [START_REF] Recht | Guaranteed Minimum-Rank Solutions of Linear Matrix Equations via Nuclear Norm Minimization[END_REF][START_REF] Meka | Guaranteed Rank Minimization via Singular Value Projection[END_REF][START_REF] Jain | Low-rank Matrix Completion using Alternating Minimization[END_REF].

Computing the Bregman gradient map. We now show that, when there is no penalty term g, the Bregman gradient map of h G can be computed efficiently, as it involves solving an easy quartic minimization subproblem of size r.

Proposition 2.6 (Bregman gradient map, Gram kernel). Assume that g ⌘ 0. For X 2 R n⇥r , the Bregman gradient map of f for the Gram kernel h G with step size > 0 is given by

T (X) = V [↵ Tr(Z)I r + Z + I r ] 1
where the matrices V, Z are computed through the routine:

• set V = rh G (X) rf (X), • diagonalize V T V as V T V = P T DP where P 2 O r and D = diag(⌘ 2 1 , . . . , ⌘ 2 r ),
• let µ = (µ 1 , . . . , µ r ) be the unique solution of the convex minimization problem

min x2R r (x) := ↵ 4 kxk 4 + 4 r X i=1 x 4 i + 2 kxk 2 r X i=1 ⌘ i x i , • finally set Z = P T diag ⇥ µ 2 1 , . . . , µ 2 r ⇤ P.
Proof. When g ⌘ 0, The Bregman gradient map of h G writes, for X 2 R n⇥r ,

T (X) = argmin U 2R n⇥r ⇢ hrf (X), U Xi + 1 D h G (U, X) = argmin U 2R n⇥r {h G (U ) hV, U i} (2.8)
where we remove constant terms and defined V := rh G (X) rf (X). Write for the sake of clarity U ? := T (X). The optimization problem (2.8) is strictly convex and the unique solution

U ? satisfies rh G (U ? ) = V , meaning that U ? ↵kU ? k 2 I r + U ?T U ? + I r = V.
(2.9)

Define Z := U ?T U ? 2 R r⇥r . Then, the knowledge of Z determines U ? , since kU ? k 2 = Tr(Z) and therefore U ? = V (↵ Tr(Z)I r + Z + I r ) 1 . Now, taking (2.9) and multiplying by its transpose implies that

(↵kU ? k 2 I r + Z + I r ) 2 Z = V T V. (2.10)
This shows that V T V is a polynomial in Z, and therefore that they admit the same eigenvectors. Write the diagonalization

V T V = P T diag(⌘ 2 1 , . . . , ⌘ 2 r )P, Z = P T diag(µ 2 1 , . . . , µ 2 r )P
where P 2 O r and µ i , ⌘ i 0 for i = 1 . . . r. It follows from diagonalizing (2.10) and taking the square root that 0

@ ↵ 0 @ r X j=1 µ 2 j 1 A + µ 2 i + 1 A µ i = ⌘ i 8i = 1, . . . , r
This is exactly the first-order optimality condition on µ = (µ 1 , . . . , µ r ) for the problem

µ = argmin x2R r ↵ 4 kxk 4 + 4 r X i=1 x 4 i + 2 kxk 2 r X i=1 ⌘ i x i .
(2.11)

Note that we do not need to enforce the nonnegativity constraint on x, since we chose ⌘ i 0 it follows that the optimal solution will be nonnegative. Hence, we can reconstruct Z from the diagonalization of V T V and the solution of Problem (2.11), and thus we get the procedure described in the theorem for computing U ⇤ = T (X).

Complexity. Note that the order of multiplication is important: we only need to compute the eigendecomposition of V T V , which is of size r ⇥ r. We additionally need to solve a small minimization problem of size r, which can be done efficiently using the quartic Bregman gradient algorithm with norm kernel (see Appendix 2.5 for implementation details). Due to this, the complexity of computing the Bregman gradient map of h G is O(nr 2 + r 3 + Kr), where K is the number of iterations needed to solve the subproblem. Since r is usually much smaller than n by several orders of magnitude, the main computational bottleneck remains in most applications computing the gradient rf (X).

How to choose the most appropriate kernel?

In order to devise efficient methods, one should search for the kernel h such that the upper approximation of f in (RelSmooth) is as tight as possible, or, equivalently, such that the Hessian of the residual Lh f is small. On the other hand, h has to be simple enough so that the Bregman gradient map (2.3) is easy to compute (which precludes choosing h = f , as the iteration would be as hard to solve as the initial problem). This trade-off is key in choosing the appropriate kernel. Let us review these two conflicting criteria in our situation.

Complexity of the Bregman gradient map. For the norm kernel h N , one iteration involves computing the gradient of f , then solving a simple scalar equation. The Gram kernel h G involves solving a subproblem which requires O(nr 2 + r 3 ) additional operations. This overhead is negligible for the typical regime where r ⌧ n; however, the iterate can be computed easily only for unconstrained problems.

Quality of Hessian approximation. We showed in Section 2.2.3 that the quality of the approximation of the quartic component of f by the Gram kernel is bounded by

O(1 µ F /L F ).
Therefore, it is expected to show good performance when F is sufficiently well-conditioned. The norm kernel, however, has no such property, as its approximation of f is much coarser. The difference stems from the supplementary kX T Xk 2 term, which can be much smaller than kXk 4 , especially when the columns of X are nearly orthogonal. Note that even if F is not globally strongly convex or µ F is unknown, the Gram kernel can take advantage of local strong convexity through adaptive step sizes, as we show in the sequel.

Algorithms for quartic low-rank minimization

Now that we are equipped with a non-Euclidean geometry induced by one of the kernels h N and h G , we are ready to define the minimization scheme Dyn-NoLips in Algorithm 3. It extends the Bregman proximal gradient, a.k.a NoLips, algorithm from [START_REF] Bolte | First Order Methods Beyond Convexity and Lipschitz Gradient Continuity with Applications to Quadratic Inverse Problems[END_REF] to allow step sizes larger than the theoretical value 1/L.

Algorithm 3 Dyn-NoLips

Input: A distance kernel h such that f is smooth relative to h and a maximal step size max

Initialize X 0 2 R n⇥r such that (X 0 ) < 1. for k = 1,2,. . . do
Choose a step size k  max such that the sufficient decrease condition (2.12) holds Set X k = T k (X k 1 ) end for

Step size choice The step size k is chosen so that the new iterate

X k = T k (X k 1 ) satisfies f (X k )  f (X k 1 ) + hrf (X k 1 ), X k X k 1 i + 1 k D h (X k , X k 1 ).
(2.12)

There are two ways to ensure this condition holds.

• Fixed step size. Since f is L-smooth relative to h, (2.12) holds as soon as 0 < k  1/L.

• Dynamical step size. In some cases, the relative Lipschitz constant might be too conservative, and better numerical performance can be achieved by taking larger steps.

We therefore can use a dynamical strategy for extending the step size, ensuring that (2.12) holds at each iteration. There are many strategies to efficiently adjust the step size; see, e.g., [START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF]. In our case, we choose a simple strategy similar in spirit to the Armijo line search: at iteration k, start with a tentative step size k , then find the smallest nonnegative integer j such that (2.12) is satisfied with step size 2 j k . Then, set k+1 = 2 j+1 k .

Convergence to a stationary point. We now extend the theoretical convergence results from [START_REF] Bolte | First Order Methods Beyond Convexity and Lipschitz Gradient Continuity with Applications to Quadratic Inverse Problems[END_REF] to handle the dynamical step size strategy.

Theorem 2.7 (Convergence results). Let {X k } k 0 be the sequence generated by Algorithm 3. Assume that 1. f is L-smooth relative to a distance kernel h such that h is strongly convex and twice continuously differentiable on R n⇥r , and the penalty function g is convex,

2. max 1/(2L),

3. the function = f + g is coercive (meaning that (X) ! +1 when kXk ! +1) and semialgebraic.

Then, the sequence { (X k )} k 0 is nonincreasing, and the sequence {X k } k 0 converges towards a critical point X ⇤ of problem .

Proof. First, the step size k can be bounded for k 0 as

1 2L  k  max .
Indeed, the upper bound holds by construction of the algorithm. The lower bound comes from the relative smoothness property: condition (2.12) is true for every 2 (0, 1 L ], so the inner loop will stop whenever gets below 1/L.

Let us now prove the result. Since Condition (2.12) holds at each iteration k, we can write

f (X k+1 )  f (X k ) + hrf (X k ), X k+1 X k i + 1 k D h (X k+1 , X k ).
(2.13)

On the other hand, the optimality condition characterizing

X k+1 = T k (X k ) writes 0 2 k ⇣ @g(X k+1 ) + rf (X k ) ⌘ + rh(X k+1 ) rh(X k ), (2.14)
where @g denotes the subdifferential of the convex function g. Combining (2.14) with the subgradient inequality for g yields

g(X k+1 )  g(X k ) + 1 k hrh(X k ) rh(X k+1 ), X k+1 X k i (2.15) hrf (X k ), X k+1 X k i.
Summing (2.13) and (2.15) gives

(X k+1 )  (X k ) + 1 k ⇥ D h (X k+1 , X k ) + hrh(X k ) rh(X k+1 ), X k+1 X k i ⇤ ,
which yields

(X k+1 )  (X k ) 1 k D h (X k , X k+1 ).
From this inequality, we can now prove the same convergence properties as for the standard BPG scheme. Indeed, the monotonicity of the sequence

{ (X k )} k 0 is a direct consequence of the above. Since k  max , it follows that at every iteration k 0, (X k ) (X k+1 ) 1 max D h (X k , X k+1 ).
Now, this inequality is the same as the one needed to prove convergence in the case of the fixed step size in [START_REF] Bolte | First Order Methods Beyond Convexity and Lipschitz Gradient Continuity with Applications to Quadratic Inverse Problems[END_REF]. Thus, global convergence towards a critical point is a consequence of [Bolte et al., 2018, Th. 4.1], since all the assumptions are met: the kernel h is defined over the entire space R n⇥r , it is strongly convex, and rh is Lipschitz continuous on bounded subsets of R n⇥r (because we assumed it is C 2 ). We also need the fact that the sequence {X k } k 0 is bounded, which is a consequence of the monotonicity of { (X k )} k 0 and the fact that the function is coercive.

The semialgebraicity assumption is needed to establish the crucial nonsmooth Lojasiewicz property [START_REF] Bolte | The Łojasiewicz Inequality for Nonsmooth Subanalytic Functions with Applications to Subgradient Dynamical Systems[END_REF], required to show convergence to a critical point. It holds for all the applications we cited, since the class of semialgebraic functions includes polynomial functions, `1 and `2 norms, the `0 seminorm and indicators of polynomial sets.

Applications

We now illustrate applications of our methodology to two different low-rank problems, symmetric nonnegative matrix factorization and Euclidean distance matrix completion. We show that good numerical performance can be reached using the dynamical step strategy, and that, for Euclidean matrix completion, it can be further improved by using the Gram kernel.

Symmetric Nonnegative Matrix Factorization

Symmetric Nonnegative Matrix Factorization (SymNMF) is the task of finding, given a symmetric nonnegative matrix M 2 R n⇥n , a nonnegative matrix X 2 R n⇥r such that M ⇡ XX T . This is done by solving

min 1 2 kM XX T k 2 F subject to X 0 (SymNMF)
in the variable X 2 R n⇥r , where the inequality constraint is meant componentwise and r  n is the target rank.

(SymNMF) is used as a probabilistic clustering or graph clustering technique [START_REF] Ding | On the Equivalence of Nonnegative Matrix Factorization and Spectral Clustering[END_REF][START_REF] He | Symmetric nonnegative matrix factorization: Algorithms and applications to probabilistic clustering[END_REF]. Numerical experiments by [START_REF] Kuang | SymNMF: nonnegative low-rank approximation of a similarity matrix for graph clustering[END_REF] have shown that it achieves state-of-the-art clustering accuracy on several text and image datasets.

Solving SymNMF

While (SymNMF) looks similar to the well-known asymmetric NMF problem min

X,Y 1 2 kM XY T k,
it is actually harder. This is because NMF has a favorable block structure that allows the application of efficient alternating algorithms [START_REF] Kim | Accelerated Proximal Point Method and Forward Method for Monotone Inclusions[END_REF]Park, 2013, Cichocki and[START_REF] Cichocki | Fast local algorithms for large scale nonnegative matrix and tensor factorizations[END_REF]. SymNMF, however, does not enjoy the same block structure. Current solvers fall into two categories: Direct solvers. There have been several attempts at solving the original problem, including multiplicative update rules [START_REF] He | Symmetric nonnegative matrix factorization: Algorithms and applications to probabilistic clustering[END_REF], projected gradient algorithm quasi-Newton schemes [START_REF] Kuang | SymNMF: nonnegative low-rank approximation of a similarity matrix for graph clustering[END_REF], and coordinate descent [START_REF] Vandaele | Efficient and non-convex coordinate descent for symmetric nonnegative matrix factorization[END_REF].

Nonsymmetric relaxations. Another idea is to use a mere penalty method [START_REF] Kuang | SymNMF: nonnegative low-rank approximation of a similarity matrix for graph clustering[END_REF][START_REF] Lu | A Nonconvex Splitting Method for Symmetric Nonnegative Matrix Factorization : Convergence Analysis and Optimality[END_REF][START_REF] Zhu | Dropping Symmetry for Fast Symmetric Nonnegative Matrix Factorization[END_REF], relaxing (SymNMF) to the following penalized nonsymmetric problem min

1 2 kM XY T k 2 F + µkX Y k 2 F subject to X, Y 0, (P-NMF)
in the variables X, Y 2 R n⇥r , with parameter µ 0. This formulation is very similar to asymmetric NMF and can be solved by the same fast alternating algorithms that exploit the block structure, such as Alternating Nonnegative Least Squares (ANLS) and Hierarchical Alternating Least Squares [START_REF] Zhu | Dropping Symmetry for Fast Symmetric Nonnegative Matrix Factorization[END_REF] (HALS), which are arguably the fastest SymNMF solvers.

Applying NoLips. We propose to apply the NoLips/Bregman proximal gradient algorithm for optimizing the original objective function. Problem (SymNMF) falls within our framework with F (Y ) = 1 2 kM Y k 2 , which has a Lipschitz gradient with constant 1, and g(X) = i{X 0} the indicator function of the nonnegative orthant. Therefore, Proposition 2.3 implies that f (X) := 1 2 kM XX T k 2 is 1-smooth relatively to the kernel h N with ↵ = 6 and = 2krF (0)k = 2kM k. Since, in addition, f is polynomial and g is the indicator of a polynomial set, f + g is semialgebraic, and it is also coercive, so Theorem 2.7 guarantees that NoLips will converge towards a stationary point of problem (SymNMF).

In this problem, the Bregman iteration map is solved by simply adding a projection step

T (X) = 1 ⌧ (↵k⇧ + (U )k 2 ) ⇧ + (U ),
where U = rh N (X) rf (X), ⌧ has been defined in Proposition 2.4 and ⇧ + is the projection on the nonnegative orthant: ⇧ + (U ) = max(U, 0) (entrywise).

Computational complexity for NoLips. The computational complexity of an iteration is dominated by gradient computations and objective function evaluations, as all other operations are linear in the size of the variable.

If M is a n ⇥ n dense matrix, each gradient and function evaluation uses O(n 2 r + nr 2 ) floating point operations. If M is represented as a sparse matrix with p << n 2 nonzero elements, then we can take advantage of this structure [Vandaele et al., 2016, Rmk. 2] by using

f (X) = 1 2 kXX T M k 2 = 1 2 kM k 2 + 1 2 kX T Xk 2 hMX, Xi rf (X) = 2X(X T X) 2MX
which yields a much improved O (r 2 + p)n complexity per iteration.

Numerical experiments

We implemented the following algorithms: Algorithm 3 with dynamical step size and the norm kernel (Dyn-NoLips), the -SNMF scheme from [START_REF] He | Symmetric nonnegative matrix factorization: Algorithms and applications to probabilistic clustering[END_REF], where we set = 0.99 as advised by the authors, the projected gradient algorithm (PG) with Armijo line search from [START_REF] Kuang | SymNMF: nonnegative low-rank approximation of a similarity matrix for graph clustering[END_REF], where we use the line search parameters = 0.1 and = 0.01, the coordinate descent scheme (CD) from [START_REF] Vandaele | Efficient and non-convex coordinate descent for symmetric nonnegative matrix factorization[END_REF], the ADMM algorithm [START_REF] Lu | A Nonconvex Splitting Method for Symmetric Nonnegative Matrix Factorization : Convergence Analysis and Optimality[END_REF], and the two fast algorithms from [START_REF] Zhu | Dropping Symmetry for Fast Symmetric Nonnegative Matrix Factorization[END_REF] for solving the penalized problem (P-NMF): SymANLS and SymHALS. For the last two, we tuned the µ penalization parameter for best performance. We left out the quasi-Newton algorithm from [START_REF] Kuang | SymNMF: nonnegative low-rank approximation of a similarity matrix for graph clustering[END_REF] because of its prohibitive O(n 3 ) complexity for large datasets.

All algorithms were implemented in Julia [Jeff Bezanson, Alan Edelman and Shah, 2017] which is a highly-optimized numerical computing language. Since our algorithms have different complexity per iteration, it is essential to compare them in terms of running time, and Julia provides a fairly accurate way to do so as there is little interpreter overhead in loops. Tests were run on a PC Intel CORE i7-4910MQ CPU @ 2.90 GHz x 8 with 32 Go RAM.

We used two image and two text datasets.

• Image.

-CBCL 1 : 2,429 images of faces of size 19 ⇥ 19 -Coil-20 2 : 1440 images of size 128⇥128 representing 20 objects under various angles.

• Text.

-TDT2 3 : dataset of 11,201 news articles classified in 96 semantic categories. -Reuters 3 : dataset of news articles, which we restricted to the largest 25 categories, leaving a total of 7,963 documents.

For all image and text datasets, we construct a sparse similarity matrix M following the procedure described in [Kuang et al., 2015, Section 7.1]. We begin by computing the similarity graph between data points, using cosine similarity on term frequency vectors for text, and a Gaussian kernel for image (with the self-tuning method for the scale). The graph obtained is sparsified by keeping only the edges connecting the k-nearest neighbors, with k = blog 2 nc + 1. Then, M is taken as a normalized version of the graph adjacency matrix.

We use the usual convergence criterion for constrained nonconvex problems where r P f (X) is the projected gradient defined as

kr P f (X k )k kr P f (X 0 )k  ✏ (2.16)
(r P f (X)) ij = ( rf (X) ij if X ij > 0, min (rf (X) ij , 0) if X ij = 0.
Table 2.1 reports the average time needed to reach a convergence criterion of ✏ = 10 3 , for 10 random initializations. For each dataset, we test several values for the rank parameter r. In addition, Figure 1 shows the average evolution of the normalized objective gap f (X k ) f min /kM k 2 , where f min is the minimal objective value encountered in all initializations.

Overall, the algorithm that shows the best convergence speed is SymHALS, but it has the disadvantage of needing to tune the penalization parameter µ. In the experiments we report, small values of µ yielded optimal performance, while the convergence theory of [START_REF] Zhu | Dropping Symmetry for Fast Symmetric Nonnegative Matrix Factorization[END_REF] only holds for large values for which the algorithm is much slower. By contrast, Dyn-NoLips is hyperparameter-free and has the second best overall performance. The gap with the other methods is particularly significant on the larger TDT2 and Reuters datasets, showing that the method scales well with problem dimension.

Euclidean Distance Matrix Completion

Euclidean distance matrix completion (EDMC) is the task of recovering the position of n points x ⇤ 1 , . . . , x ⇤ n 2 R r , given the knowledge of a partial set of pairwise distances

d ij = kx ⇤ i x ⇤ j k 2 for (i, j) 2 ⌦, where ⌦ ⇢ [1, n] ⇥ [1, n]. It is a fundamental problem with applications in
X (i,j)2⌦ kX i X j k 2 d ij 2 (EDMC)
in the variable X 2 R n⇥r . It can be rewritten f (X) = 1 2 kP ⌦ ((XX T ) D)k 2 where D is the matrix of known distances, P ⌦ denotes the projection operator such that P ⌦ (Y ) ij = Y ij if (i, j) 2 ⌦, and P ⌦ (Y ) ij = 0 elsewhere, and  is the linear operator defined for Y 2 R n⇥n by

(Y ) ij = Y ii + Y jj 2Y ij for 1  i, j  n.
Applying NoLips with the norm kernel. Problem (EDMC) falls within our framework with

F (Y ) = 1 2 kP ⌦ ((Y ) D)k 2
, which can be shown to have a Lipschitz gradient with constant L EDM := 9 max i=1...n |{j|(i, j) 2 ⌦}|. Therefore, as in the case of SymNMF, the norm kernel h N can be used with a initial step size 1 and parameters ↵ = 6L EDM and = 1 3 krF (0

)k = 2kP ⌦ (D)k.
Using the Gram kernel As the problem is unconstrained, we can also apply minimization using the Gram kernel h G . We use the parameters ↵ = 2L EDM , = L EDM and = 2kP ⌦ (D)k, which ensure that f is 1-smooth relatively to h G by Proposition 2.5.

Computational complexity for NoLips. As before, the main computational bottleneck for an iteration consists in computing the value and gradient of the objective function. If p = |⌦| denotes the number of known distances, then the computational complexity is O(pr). If the Gram kernel is used, each iteration requires an additional O(nr 2 + r 3 ) flops (See Section 2.2.3), which is negligible compared to the latter in the usual setting where p >> n and r is small.

Numerical experiments

We implement the following algorithms: NoLips with a dynamical step size and the norm kernel (Dyn-NoLips), NoLips with a dynamical step size and the Gram kernel (Dyn-NoLips-Gram), gradient descent with Armijo line search (GD), the Riemannian trust region algorithm from [START_REF] Mishra | Low-rank optimization for distance matrix completion[END_REF] (TR). We leave out semidefinite relaxations because of their memory requirement which is prohibitive on large data. As the implementation for TR is provided in Matlab, we run our experiments on Matlab as well, with the same setup as in Section 2.4.1.

We try the algorithms on a standard EDMC problem, the 3-dimensional Helix dataset [START_REF] Mishra | Low-rank optimization for distance matrix completion[END_REF] which is generated as X i = (cos(3t i ), sin(3t i ), 2t i ) where {t i } n i=1 are sampled uniformly in [0, 2⇡]. We randomly keep only 10 % on the pairwise distances, and test on two different problem sizes: n = 2000 and n = 5000. Figure 2.2 reports the normalized root mean squared error (RMSE) over all distances (known and unknown) averaged on 10 random initializations. All the algorithms manage to recover the ground truth; the Dyn-NoLips-Gram algorithm shows the best numerical performance, which demonstrates the advantage of using the Gram geometry.

Conclusion

We proposed a generic approach for solving Burer-Monteiro formulations of low rank minimization problems using the methodology of Bregman gradient methods and relative smoothness. We studied two quartic kernels, including a new Gram kernel, and demonstrated their benefits on numerical experiments. In future work, performance could be improved further by studying inertial variants [START_REF] Chandra Mukkamala | Convex-Concave Backtracking for Inertial Bregman Proximal Gradient Algorithms in Non-Convex Optimization[END_REF][START_REF] Hanzely | Accelerated Bregman Proximal Gradient Methods for Relatively Smooth Convex Optimization[END_REF]. New kernels could also be explored beyond the class of quartic functions to tackle other problems with inherent non-Euclidean geometries.

Appendix: solving the subproblem for computing the Bregman gradient map of the Gram kernel While it seems that computing the Bregman iteration map of the Gram kernel involves solving another difficult quartic subproblem, it is actually of small size (r is typically not larger than a few dozens) and can be solved efficiently with the NoLips scheme.

Indeed, the objective function of problem (2.11) is 1-smooth relatively to the norm kernel in R r h N (x) = ↵u 4 kxk 4 + u 2 kxk 2 with a choice of parameters ↵ u = ↵ + 3 and u = . Algorithm 4 details the procedure. We initialize µ with the values for the previous iteration of the outer procedure. This proves to be efficient as the values will not vary much from one iteration to another. For the stopping criterion, we use the scaled gradient norm kr (v)k/k⌘k and a tolerance value ✏ = 10 6 .

The subproblem being very well conditionned, it is minimized easily; in numerical experiments, it usually convergences in no more than 20 iterations.

Algorithm 4 Computing the Bregman iteration map of the Gram kernel

Input: Matrix X 2 R n⇥r , gradient of the objective rf (X), step size > 0, parameters ↵, , > 0, subproblem tolerance ✏, and (optionally), values µ of µ computed at the previous iteration.

Form V = rh G (X) rf (X) = ↵kXk 2 I r + X T X + I r X rf (X) Compute V T V Form the eigendecomposition of V T V = P T DP where P 2 O r and D = diag(⌘ 2 1 , . . . , ⌘ 2 r )
Initialize µ as µ if provided, and as (0, . . . , 0) otherwise. repeat Compute r (µ) where r (µ

) i = ↵kµk 2 µ i + µ 3 i + µ i ⌘ i Compute rh N (µ) where rh N (µ) i = (↵ + 3 )kµk 2 µ i + µ i Form v = rh N (µ) r (u) Set µ ⇥ ⌧ (↵ + 3 )kvk 2 ⇤ 1 v
where ⌧ has been defined in Proposition 2.4 until stopping criterion has been satisfied, i.e., kr

(v)k/k⌘k < ✏ Form Z = P T diag(µ 2 1 , . . . , µ 2 r )P Compute T (X) = V [↵ Tr(Z)I r + Z + I r ] 1 Output: Bregman gradient iterate T (X)
Chapter 3

A Lower Bound for Relatively-Smooth Convex Mimization

Chapter Abstract

We prove a lower bound showing that the O(1/k) convergence rate of the Bregman gradient descent method (a.k.a. NoLips or mirror descent) is optimal for the class of problems satisfying the relative smoothness assumption, among all first-order methods that use Bregman divergences and linear operations. As a consequence, no algorithm of this class can achieve a better rate than Bregman gradient descent for generic kernel functions.

Reference: this chapter is based on a publication in Mathematical Programming [Dragomir et al., 2021c]. Part of this work has been done in collaboration with Adrien Taylor.

Introduction

We consider the constrained minimization problem

min x2C f (x) (P)
where f is a convex continuously differentiable function and C is a nonempty closed convex subset of R d . In large-scale settings, first-order methods are particularly popular due to their simplicity and their low cost per iteration. The (projected) gradient descent (PG) is a classical method for solving (P), and consists in successively minimizing quadratic approximations of f , with

x k+1 = argmin u2C f (x k ) + hrf (x k ), u x k i + 1 2 ku x k k 2 , (PG)
where k • k is the Euclidean norm. Although standard, there is often no good reason for making such approximations, beyond our capability of solving this intermediate optimization problem.

In other words, this traditional approximation typically does not reflect neither the geometry of f nor that of C. A powerful generalization of PG consists in performing instead a Bregman gradient step

x k+1 = argmin u2C f (x k ) + hrf (x k ), u x k i + 1 D h (u, x k ), (BGD)
where the Euclidean distance has been replaced by the Bregman divergence

D h (x, y) := h(x) h(y) hrh(y), x yi
induced by some strictly convex and continuously differentiable kernel function h. A wellchosen h allows designing first-order algorithms adapted to the geometry of the constraint set and/or the objective function. Of course, a conflicting goal is to choose h such that each iteration (BGD) can be solved efficiently in practice, discarding choices such as h = f (for which performing an iteration would be as hard as solving the original problem).

Recently, Bauschke et al.

[2017] introduced a natural condition for analyzing this scheme, assuming that the inner objective in the iteration (BGD) is an upper bound on f . This ensures that performing an iteration decreases the function values f (x k ). This assumption, known as relative smoothness, generalizes the standard L-smoothness assumption implied by Lipschitz continuity of rf . The Bregman gradient algorithm, also called NoLips in the setting of [START_REF] Heinz | A Descent Lemma Beyond Lipschitz Gradient Continuity: First-Order Methods Revisited and Applications[END_REF], is thus a natural extension of gradient descent (PG) to objective functions whose geometry is better modeled by a non-quadratic kernel h. Practical examples of relative smoothness arise in Poisson inverse problems [START_REF] Heinz | A Descent Lemma Beyond Lipschitz Gradient Continuity: First-Order Methods Revisited and Applications[END_REF], quadratic inverse problems [START_REF] Bolte | First Order Methods Beyond Convexity and Lipschitz Gradient Continuity with Applications to Quadratic Inverse Problems[END_REF], rank minimization [Dragomir et al., 2021a] and regularized higherorder tensor methods [START_REF] Nesterov | Implementable tensor methods in unconstrained convex optimization[END_REF].

Can we accelerate Bregman gradient descent? In the Euclidean setting where h is the squared Euclidean norm 1 2 k • k 2 , accelerated projected gradient methods exhibit faster convergence than the vanilla projected gradient algorithm. These methods, which can be traced back to [START_REF] Nesterov | A Method of Solving A Convex Programming Problem With Convergence rate O(1/kˆ2)[END_REF], are proven to be optimal for L-smooth functions and have found a number of successful applications, in e.g., imaging [START_REF] Beck | A Fast Iterative Shrinkage-Thresholding Algorithm[END_REF]. A natural question is therefore to understand whether the (BGD) algorithm can be accelerated in the relatively-smooth setting.

This question has been raised in several works, including [Bauschke et al., 2017, Section 6], [Lu et al., 2018, Section 3.4], and the survey [Teboulle, 2018, Section 6]. Partial answers have already been provided under somewhat strict additional regularity assumptions (see e.g., Auslender and [START_REF] Auslender | Interior Gradient and Proximal Methods for Convex and Conic Optimization[END_REF], [START_REF] Walid | Accelerated Mirror Descent In Continuous and Discrete Time[END_REF], [START_REF] Hanzely | Accelerated Bregman Proximal Gradient Methods for Relatively Smooth Convex Optimization[END_REF] and discussions in the sequel), while the general case was apparently still open, and relevant in practical applications.

In this work, we establish a lower complexity bound proving that BGD is optimal for the general relatively-smooth setting, and therefore that generic acceleration is impossible.

In order to do so, we adopt the standard black-box model used for studying complexity of first-order methods [START_REF] Nemirovski | Problem Complexity and Method Efficiency in Optimization[END_REF]. We consider that both f and h are described by first-order oracles, so as to obtain generic complexity results, and we look for worst-case couples of functions (f, h) satisfying the relative smoothness assumption. A central idea in our approach is the fact that, when studying the worst-case behavior of Bregman methods in the relatively-smooth setting, f and h can get arbitrarily close to some limiting pathological nonsmooth functions.

Main inspiration: Performance Estimation Problems. The worst-case instance used for proving our lower bound was inferred from the solution to a Performance Estimation Problem (PEP). The PEP methodology allows to study worst-case behavior of first-order methods by solving appropriate semidefinite programs, and was pioneered by [START_REF] Drori | Performance of First-Order Methods for Smooth Convex Minimization: A Novel Approach[END_REF] in the context of smooth convex minimization. In Chapter 4, we show how the PEP technique can be extended to Bregman first-order methods in the relatively-smooth setting, and how it allowed us to discover the key elements for proving our lower bound. However, we emphasize that the proof we give in this chapter is self-contained and can be read independently.

Related work on lower complexity bounds. The first-order black-box model, developed initially in the works of [START_REF] Nemirovski | Problem Complexity and Method Efficiency in Optimization[END_REF] and later [START_REF] Nesterov | Introductory Lectures on Convex Optimization: A Basic Course[END_REF] has allowed to prove optimal complexity results for several classes of problems in first-order optimization [START_REF] Drori | The Exact Information-Based Complexity of Smooth Convex Minimization[END_REF]. These results usually rely on well-chosen worst-case functions whose structure makes them difficult to minimize for all methods within a given class.

Our worst-case instances are obtained from pointwise maxima of affine functions, reminiscent of lower bounds for nonsmooth convex minimization [START_REF] Nemirovski | Efficient methods for large-scale convex optimization problems[END_REF]Yudin, 1983, Woodworth and[START_REF] Woodworth | Lower Bound for Randomized First Order Convex Optimization[END_REF]. Our construction then involves smoothing those piecewise affine functions, making them differentiable. This technique is also used in the very related work of [START_REF] Guzmán | On Lower Complexity Bounds for Large-Scale Smooth Convex Optimization[END_REF], which studies lower bounds for minimization of convex functions that are smooth with respect to `p norms. To the best of our knowledge, the lower bound obtained in the sequel is not a particular case of those in [START_REF] Guzmán | On Lower Complexity Bounds for Large-Scale Smooth Convex Optimization[END_REF], as smoothness with respect to a certain norm is different from relative smoothness with respect to the same (squared) norm, beyond the `2-norm.

Notation, We use C to denote the closure of a set C, int C for its interior and @C for its boundary. We denote (e 1 , . . . , e n ) the canonical basis of R d , and for p 2 {1, . . . n} we write E p = Span(e 1 , . . . , e p ) the set of vectors supported by the first p coordinates. Subscripts on a vector denote the iteration counter, while a superscript such as x (i) denotes the i-th coordinate.

Algorithmic setup

In this section, we briefly recall the base ingredients and technical assumptions on f and h that are used within Bregman first-order methods. For a more detailed presentation and examples of applications, we refer the reader to Chapter 1.

Kernel functions

Let C be a nonempty closed convex subset of R d . The first step in defining Bregman methods is the choice of a kernel (or reference) function h on C.

Definition 3.1 (Kernel function). A function h : R d ! R [ {+1} is called a kernel function on C if (i) h is closed convex proper (c.c.p.),
(ii) h is continuously differentiable and strictly convex on int C, (iii) the Bregman gradient iterates are well-posed, i.e., for every p 2 R d , the problem min u2C hp, ui + h(u) has a unique minimizer, which belongs to int C.

A kernel function h induces a Bregman divergence D h defined as

D h (x, y) = h(x) h(y) hrh(y), x yi 8x 2 dom h, y 2 dom rh.
Note that D h is not a distance in the classical sense, however it enjoys a separation property; due to the strict convexity of h we have D h (x, y) 0 8x 2 dom h, y 2 dom rh, and

D h (x, y) = 0 iff x = y.
Examples. We list some of the most classical examples of kernel functions (see Chapter 1 for more):

• the Euclidean kernel h(x) = 1 2 kxk 2 with domain R d , and for which D h (x, y) = 1 2 kx yk 2 is the Euclidean distance,

• the Boltzmann-Shannon entropy h(x) = P i x (i) log x (i) extended to 0 by setting 0 log 0 = 0, whose domain is thus R d + ,

• the Burg entropy h(x) =

P i log x (i) with domain R d ++ ,
• the quartic kernel h(x) = 1 4 kxk 4 + 1 2 kxk 2 with domain R d [START_REF] Bolte | First Order Methods Beyond Convexity and Lipschitz Gradient Continuity with Applications to Quadratic Inverse Problems[END_REF].

Convex conjugate. If h is a kernel function, we define its convex conjugate h ⇤ as

h ⇤ (y) = sup u2R d hu, yi h(u)
If, for every y 2 R d , the supremum in the definition of h ⇤ (y) is attained, then h ⇤ is differentiable and its gradient satisfies for every u 2 dom rh ⇤ rh ⇤ (y) = argmax

u2R d hu, yi h(u).

Relatively-smooth optimization problems

We now recall the framework of relatively-smooth optimization [START_REF] Heinz | A Descent Lemma Beyond Lipschitz Gradient Continuity: First-Order Methods Revisited and Applications[END_REF][START_REF] Lu | Relatively-Smooth Convex Optimization by First-Order Methods, and Applications[END_REF] for solving the minimization problem

min x2C f (x)
For simplicity, we present the setting without nonsmooth regularization term; our lower bound is a fortiori valid for the Bregman proximal gradient algorithm designed for solving composite problems [Bauschke et al., 2017, Eq. (12)].

In addition to these assumptions, the central property we need in order to apply the Bregman gradient method is the so-called relative smoothness property [START_REF] Heinz | A Descent Lemma Beyond Lipschitz Gradient Continuity: First-Order Methods Revisited and Applications[END_REF][START_REF] Lu | Relatively-Smooth Convex Optimization by First-Order Methods, and Applications[END_REF]. Definition 3.2 (Relative smoothness). Let h be a kernel function on C, and f a function such that dom h ⇢ dom f . We say that f is smooth relative to h if it is differentiable on int C and if there exists a constant L > 0 such that Lh f is convex on dom h.

(LC)

Relative smoothness allows to build a simple global majorant of f ; indeed, (LC) is equivalent to the condition (see, e.g, Bauschke et al. [2017])

f (x)  f (y) + hrf (y), x yi + LD h (x, y) 8x 2 dom h, y 2 dom rh,
and the NoLips/Bregman gradient method consists in successively minimizing this upper approximation.

We use the following convenient notation to characterize the class of relatively-smooth problems.

Definition 3.3. We say that the couple of functions (f, h) is a relatively-smooth instance, and write

(f, h) 2 B L (C) if (i) h is a kernel function on C, (ii) f : R d ! R [ {+1} is a closed convex proper function, (iii) f is L-smooth relative to h on C.
Finally, let us denote by B L the union of B L (C) for all closed convex sets C:

B L = [ d 1 [ C⇢R d C closed convex B L (C)

The NoLips/Bregman Gradient algorithm

Previous assumptions allow defining the Bregman Gradient (BG)/NoLips algorithm for minimizing f . For simplicity, we consider here the constant step size method.

Algorithm 5 Bregman Gradient (BG) / NoLips [START_REF] Heinz | A Descent Lemma Beyond Lipschitz Gradient Continuity: First-Order Methods Revisited and Applications[END_REF] Input:

(f, h) 2 B L (C), x 0 2 int dom h, step size 2 (0, 1/L].
for k = 0,1,. . . do

x k+1 = argmin u2R d hrf (x k ), u x k i + 1 D h (u, x k ) (3.1)
end for

Using first-order optimality conditions, update (3.1) can alternatively be written as

x k+1 = rh ⇤ [rh(x k ) rf (x k )]
involving the gradient rh ⇤ which is usually referred to as the mirror map. The three operations rf, rh and rh ⇤ are the basic building blocks of Bregman-type methods, which we now define formally.

Defining a class of Bregman first-order methods

For proving a general lower bound for relatively-smooth optimization, we need to specify the oracle model and the class of methods under consideration. We adopt the first-order black-box model, where information about a function can be gained by calling an oracle returning the value and gradient of f at a given point. In the Bregman setting, we assume that we also have access to the first-order oracles of the kernel function h and its conjugate h ⇤ . Definition 3.4. An algorithm A is called a Bregman first-order algorithm if, for a given problem instance (f, h) 2 B L and number of iterations T 2 N, it generates at each time step t 2 {0, . . . , T }, a set of primal points X t and dual points Y t from the following process:

1. Set X 0 = {x 0 }, where x 0 2 int dom h is some initialization point, and

Y 0 = {rf (x 0 ), rh(x 0 )}.
2. For each t = 1, . . . T , perform one of the two following operations:

• either call the primal oracle (rf, rh) at some point x t chosen such as x t 2 Span(X t 1 ) \ dom rh and update the dual set as

Y t = Y t 1 [ {rf (x t ), rh(x t )}.
• Or call the mirror oracle rh ⇤ at some dual point y t chosen such as

y t 2 Span(Y t 1 ) with rh ⇤ (y t ) = argmin u2C h(u) hy t , ui
and update the primal set as

X t = X t 1 [ {rh ⇤ (y t )}.
3. Output some point x T 2 Span(X T ).

Such structural assumptions on the class of algorithms are classical from complexity analyses of Euclidean first-order methods and are used to prove e.g., optimality of accelerated first order methods [START_REF] Nesterov | Introductory Lectures on Convex Optimization: A Basic Course[END_REF]. Definition 3.4 is a natural extension to the Bregman setting, allowing additional uses of the oracles associated with the kernel function h. This model can often be relaxed through the use of more involved information theoretic arguments, see e.g., [START_REF] Nemirovski | Problem Complexity and Method Efficiency in Optimization[END_REF], [START_REF] Guzmán | On Lower Complexity Bounds for Large-Scale Smooth Convex Optimization[END_REF], [START_REF] Drori | The Exact Information-Based Complexity of Smooth Convex Minimization[END_REF], [START_REF] Woodworth | Lower Bound for Randomized First Order Convex Optimization[END_REF].

Here, we focus on Definition 3.4 as it is general enough to encompass all Bregman-type methods that only use oracles for rf, rh, which we call the primal oracles, the map rh ⇤ , which we call the mirror oracle, as well as linear operations. One can verify that all known Bregman gradient methods, including NoLips and inertial variants such as IGA [START_REF] Auslender | Interior Gradient and Proximal Methods for Convex and Conic Optimization[END_REF] or the recent algorithm in [START_REF] Hanzely | Accelerated Bregman Proximal Gradient Methods for Relatively Smooth Convex Optimization[END_REF], fit in this model.

Observe that, as BGD performs one primal oracle call and one mirror call per iteration, an iteration of BGD corresponds actually to two time steps of the formal procedure in Definition 3.4. This is why, in order to avoid ambiguity, we state our lower bound as a function of the number of oracle calls.

Convergence rate and optimality of Bregman gradient descent

In this section, we start by recalling the O(1/k) convergence rate bound for the Bregman gradient/NoLips algorithm in the setting where (f, h) 2 B L (C). We then proceed to prove that NoLips is an optimal algorithm for the class B L (C), by showing that this rate is also a lower bound for a generic class of Bregman gradient algorithms that we define below. The key elements for proving the lower bound were empirically discovered through the solution to a Performance Estimation Problem (PEP), which is detailed in Chapter 4.

Upper bound

We first recall the O(1/k) convergence rate for Bregman gradient descent/NoLips.

Theorem 3.5 (NoLips convergence rate). Let L > 0, C be a nonempty closed convex subset of R d and (f, h) 2 B L (C) be an relatively-smooth instance. Then the sequence {x k } k 0 generated by Algorithm 5 with constant step size 2 (0, 1/L] satisfies for all k 0

f (x k ) f (u)  D h (u, x 0 ) k (3.2)
for every u 2 dom h.

Proof. See [Lu et al., 2018, Thm 3.1.]; we also give an alternative simple proof in Chapter 4 (Section 4.4.1), whose analytical form has been inferred from the solution to a performance estimation problem.

Faster algorithms under additional assumptions. It is natural to ask whether an accelerated Bregman algorithm can be obtained, with a better convergence rate than O(1/k). This has already been achieved under additional regularity assumptions, as follows:

• in the Euclidean setting, when h(x) = 1 2 kxk 2 and f is L-smooth, the seminal accelerated gradient method of [START_REF] Nesterov | A Method of Solving A Convex Programming Problem With Convergence rate O(1/kˆ2)[END_REF] enjoys a O(1/k 2 ) convergence rate, which is optimal for this class of functions [START_REF] Nesterov | Introductory Lectures on Convex Optimization: A Basic Course[END_REF].

• When h is a strongly convex kernel with closed domain and f is L-smooth (which, as discussed in Section 1.2.1, is a particular case of relative smoothness), the Improved Interior Gradient Algorithm (IGA) of [START_REF] Auslender | Interior Gradient and Proximal Methods for Convex and Conic Optimization[END_REF] also admits a O(1/k 2 ) convergence rate using the same momentum technique as Nesterov-type methods.

• Recently, [START_REF] Hanzely | Accelerated Bregman Proximal Gradient Methods for Relatively Smooth Convex Optimization[END_REF] proposed an accelerated Bregman proximal gradient algorithm with rate O(1/k ), where 2 [1, 2] is determined by some crucial triangle scaling property of the Bregman divergence, whose genericity is unclear.

However, the existence of an accelerated algorithm for the general relatively-smooth setting was still an open question prior to this work. Indeed, many applications such as Poisson inverse problems [START_REF] Heinz | A Descent Lemma Beyond Lipschitz Gradient Continuity: First-Order Methods Revisited and Applications[END_REF] or D-optimal design [START_REF] Lu | Relatively-Smooth Convex Optimization by First-Order Methods, and Applications[END_REF] do not satisfy the supplementary assumptions made in the works mentioned above. In the next section, we prove that, up to a constant factor of 2, the bound (3.2) is not improvable in general for Bregman-type methods, making NoLips an optimal algorithm in the black box setting for (f, h) 2 B L .

A lower bound for relatively-smooth Bregman optimization

We show in Theorem 3.12 below that for any ✏ 2 (0, 1) and number of oracle calls N , there is a pair of functions (f, h) 2 B L (R 2N +1 ) and some x 0 2 R 2N +1 such that for any Bregman gradient algorithm initialized at x 0 , the output x N returned after performing at most N oracle calls satisfies

f (x N ) min R 2N +1 f (1 ✏) LD h (x 0 , x ⇤ ) 2N + 1 .
Proof intuition. For finding an instance (f, h) which is difficult for all Bregman methods, we use two main ideas. The first is the well-known technique used by [START_REF] Nesterov | Introductory Lectures on Convex Optimization: A Basic Course[END_REF] for proving that O(1/k 2 ) is the optimal complexity for L-smooth convex minimization. He defines a "worst function in the world" that allows any gradient method to discover only one dimension per iteration, hence hiding the minimizer from the algorithm in the remaining unexplored dimensions.

The second idea is more specific to our setting, and relies on the fact that the set of relativelysmooth problems B L (C) is not closed. In particular, a limit of differentiable functions need not be differentiable. Thence, we actually build a worst-case sequence of differentiable functions parameterized by some parameter µ, whose limit when µ ! 0 is a nonsmooth pathological function.

Choosing the objective function. Let us fix a dimension d 1 and a positive constant

⌘ > 0. Define the convex function f for x 2 R d by f (x) = max i=1,...,d |x (i) 1 ⌘ i | = kx x ⇤ k 1
which has an optimal value f⇤ = 0 attained at x ⇤ := (1 + ⌘, 1 + ⌘ 2 , . . . , 1 + ⌘ d ). The behavior of f as a pathological function comes from the fact that if at least one of the coordinates of x is zero, then f (x) f⇤ 1. Let us first prove a technical lemma about the subdifferential of f . Lemma 3.6. Let x 2 R d and v 2 @ f (x) be a subgradient of f at x.

Then (i) kvk 1  1. (ii) Let i 2 {1 . . . n}. If v (i) 6 = 0 then |x (i) x (i) ⇤ | = kx x ⇤ k 1 . Proof. Write f as f (x) = max 1id fi (x) with fi (x) = |x (i) x (i) ⇤ |.
Then, by [Nesterov, 2003, Lemma 3.1.10], we have

@ f (x) = Conv {@ fi (x)|i 2 I(x)}
where

I(x) = {i 2 {1 . . . d} | fi (x) = f (x)}.
Hence, (i) follows immediately from the well-known property that the subgradients of the absolute value lie in

[ 1, 1]. (ii) is a consequence of the fact that if v (i) 6 = 0, then i 2 I(x), which means that |x (i) x (i) ⇤ | = kx x ⇤ k 1 .
Note that f is nonsmooth hence does not meet our assumptions. We approach it with a differentiable function by considering its Moreau envelope f µ given by for n = 2, µ = 0.2 and ⌘ = 1/2. Lemma 3.8 states that if µ is small enough compared to ⌘, the behaviors of f and f µ at x 0 = 0 are the same. Indeed, the size of the smoothed region where the corners are "rounded" decreases when µ goes to 0.

f µ (x) = min u2R d f (u) + 1 2µ kx uk 2 (3.3) (a) f (b) fµ
where µ 2 (0, 1) is a small parameter. f µ is a smoothed version of f , which behaves similarly to f when we choose µ small enough. Figure 3.1 illustrates this phenomenon in two dimensions.

For general properties of the Moreau proximal envelope, we refer to [START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF]. Let us state some properties of f µ that we need for the analysis.

Lemma 3.7. f µ is a differentiable convex function, whose minimizers are the same as those of f . Its gradient at a point x 2 R d is given by rf

µ (x) = µ 1 ⇣ x prox µ f (x)
⌘ where

prox µ f (x) = argmin u2R d f (u) + 1 2µ kx uk 2
is the Moreau proximal map. Moreover, rf µ is Lipschitz continuous with constant 1/µ.

Let us now prove the central property of f µ , which states that when the last d p coordinates of x are small enough, the gradient rf µ (x) is supported on the first p + 1 coordinates. Recall that we denote (e 1 , . . . , e d ) the canonical basis of R d and write, for p 2 {1 . . . d}, E p = Span(e 1 , . . . , e p ) and E 0 = {(0, . . . , 0)}.

Lemma 3.8. Assume that µ 2 (0, 1) and ⌘ > 4µd 2 . Let p 2 {0 . . . d 1}. For any vector

x 2 R d such that max i=p+1,...,d |x (i) |  µ
we have that rf µ (x) 2 E p+1 . In addition, we have krf µ (x)k 1  1.

Proof. Take x 2 R d such that max i=p+1,...,n |x i |  µ. By Lemma 3.7, rf µ is given by

rf µ (x) = 1 µ (x prox µ f (x)) (3.4) Write y = prox µ f (x).
The optimality condition defining the proximal map yields

y x + µv = 0 (3.5)
where v 2 @ f (y), and therefore the combination of (3.4) and (3.5) implies

rf µ (x) = v 2 @ f (y). (3.6)
Now, let us assume by contradiction that rf µ (x) is not in E p+1 , meaning that there exists an index l 2 {p + 2 . . . d} such that v (l) 6 = 0. It follows from Lemma 3.6 that |(y x ⇤ ) (l) | = ky x ⇤ k 1 . Hence we have in particular that |y (l) x

(l) ⇤ | |y (p+1) x (p+1) ⇤ |. Using Condition (3.5) to replace y we get |x (l) ⇤ + µv (l) x (l) | |x (p+1) ⇤ + µv (p+1) x (p+1)
|, and recalling the definition of x ⇤ we have

|1 + ⌘ l + µv (l) x (l) | |1 + ⌘ p + 1 + µv (p+1) x (p+1) |.
By Lemma 3.6, kvk 1  1, so for all i we have 1 + µv (i) 1 µkvk 1 0. In addition, we assumed that max i=p+1,...,d |x (i) |  µ < ⌘ 4d 2 which implies ⌘ i x (i) 0 for all i p + 1. Therefore, both terms inside the absolute values are nonnegative, it follows that we can drop absolute values and

µ(v (l) v (p+1) ) ⌘ p + 1 ⌘ l + x (l) x (p+1) ⌘ • l (p + 1) l(p + 1) 2µ ⌘ l(p + 1) 2µ ⌘ d 2 2µ,
and therefore v (l) v (p+1) ⌘ µd 2
2 > 2, because we assumed ⌘ > 4µd 2 . This is a contradiction since (v (l) v (p+1) )  2kvk 1  2. Finally, the second part of the lemma is a consequence of (3.6) and kvk 1  1.

We also need the following lemma for relating the values of f and f µ .

Lemma 3.9. Let µ > 0 and x 2 R d . Then f µ (x) f (x) µ.
Proof. Write y = prox µ f (x). By definition of f µ and the proximal map we have

f µ (x) = f (y) + 1 2µ ky xk 2 f (y) = ky x ⇤ k 1 kx x ⇤ k 1 kx yk 1 .
Recall the optimality conditions defining the proximal map can be written as µ 1 (x y) 2 @f (y),and, since all subgradients of f have coordinates smaller than 1 (Lemma 3.6), we reach

kx yk 1  µ. It follows that f µ (x) kx x ⇤ k 1 kx yk 1 kx x ⇤ k 1 µ = f (x) µ,
which concludes the proof.

Choosing the kernel. As for the objective function f µ , let us pick a family of kernels h µ , whose behavior approach those of a nonsmooth function as µ ! 0.

Let us first define a unidimensional convex function µ : R ! R by

µ (t) = ⇢ t µ/2 if t µ, 1 2µ t 2 elsewhere.
Note that µ is sometimes known as the Huber function, which is a smooth approximation of the absolute value and also appears as a worst-case function for first-order methods in L-smooth minimization [START_REF] Taylor | Smooth Strongly Convex Interpolation and Exact Worst-Case Performance of First-Order Methods[END_REF].

Define

d µ : R d ! R through d µ (x) = µ 2 kxk 2 + n X i=1 µ (x (i) ), x 2 R d . (3.7)
d µ is a differentiable strictly convex function, whose gradient satisfies, for x 2 R d and i 2

{1 . . . n}, rd µ (x) (i) = µx (i) + min(1, x (i) /µ).
From the expression above, we can deduce two crucial properties that we need in the sequel: for x 2 R d and i 2 {1 . . . n}, we have rd µ (x) (i) = 0 if and only if x (i) = 0, (3.8)

|rd µ (x) (i) |  1 implies |x (i) |  µ.
(3.9)

Let L > 0. We define the kernel h µ for x 2 R d as

h µ (x) = 1 L (f µ (x) + d µ (x)) . (3.10)
By construction, Lh µ f µ is convex, so the relative smoothness property holds. It is easy to see that Definition 3.1 is satisfied as h µ is strongly convex, so we have

(f µ , h µ ) 2 B L (R d ).
Proving the zero-preserving property of the oracles. Now that the functions are defined, we are ready to prove that all oracles involved in the Bregman algorithm allow to discover only one dimension per oracle call.

Proposition 3.10 (Zero-preserving property of rf µ , rh µ , rh ⇤ µ ). Assume that µ 2 (0, 1) and ⌘ > 4µd 2 . Let p 2 {0 . . . d 1}, and x 2 R d \ E p a vector supported by the p first coordinates. Then

rf µ (x), rh µ (x), rh ⇤ µ (x) 2 E p+1 .
Proof. Let x 2 E p . Then x satisfies the assumption of Lemma 3.8 which proves that rf µ (x) 2 E p+1 . By Property (3.8) of d µ , we also have that rd µ (x) 2 E p , which allows us to conclude that

rh µ (x) = L 1 (rf µ (x) + rd µ (x)
) 2 E p+1 . It remains to prove the result for rh ⇤ µ (x). Write z = rh ⇤ µ (x), which amounts to say that rh µ (z) = x, that is rf µ (z) + rd µ (z) = Lx using (3.10). Let l 2 {p + 1 . . . d}. We have x 2 E p , hence the l th coordinate of x is zero and

rf µ (z) (l) + rd µ (z) (l) = 0.
Using the second part of Lemma 3.8, we have that krf µ (z)k 1  1; it follows that |rd µ (z) (l) |  1, which implies that |z (l) |  µ, by property (3.9) of d µ . Since this holds for any l p + 1, we have established max l=p+1,...,d

|z (l) |  µ.
Applying Lemma 3.8 to z, we obtain that rf µ (z) 2 E p+1 . Remembering that rh µ (z) = x 2 E p by construction, we get

rd µ (z) = Lrh µ (z) rf µ (z) 2 E p+1 . By Property (3.8) of d µ , it follows that z 2 E p+1 , which concludes the proof.
We can now use Proposition 3.10 inductively to state a lower bound on the performance of any Bregman gradient algorithm applied to (f µ , h µ ).

Proposition 3.11. Let N 1 and choose the dimension d = 2N + 1. Let µ 2 (0, 1) and ⌘ > 4µd 2 . Consider the functions f µ , h µ : R d ! R defined in (3.3) and (3.10) respectively. Then, for any Bregman gradient method satisfying Definition 3.4, applied to (f µ , h µ ) and initialized at x 0 = (0, . . . 0), the output x returned after performing at most N calls to each one of the primal and mirror oracles satisfies

f µ (x) min R d f µ LD hµ (x ⇤ , x 0 ) 2N + 1 • 1 µ 1 + µ + ⌘ + µ 2 (1 + ⌘) 2 .
Proof. The zero-preserving property and the structure of Bregman gradient algorithms described in Definition 3.4 implies that the set of primal points X t and dual points Y t at iteration t are supported by the t first coordinates, i.e.,

X t , Y t ⇢ E t .
Indeed, since we initialized X 0 = {x 0 } ⇢ E 0 , this follows by induction. Assume that at time t, we have X t , Y t ⇢ E t . If the primal oracle is chosen at iteration t + 1, since the query point x t+1 is taken as a linear combination of points in X t it also lies in E t , and thus Proposition 3.10 states that the new dual vectors rf µ (x t+1 ), rh µ (x t+1 ) belong to E t+1 . If, on the other hand, the mirror oracle is chosen, then with the same argument we have that y t+1 2 E t an by Proposition 3.10 that rh ⇤ µ (y t+1 ) 2 E t+1 . Now, because the algorithm has called at most N times each oracle, it has performed at most 2N steps and thus the output point satisfies x 2 E 2N , which means that x (2N +1) = 0.

We use Lemma 3.9 to relate f µ (x) and f (x). Recalling that min f µ = f⇤ = 0, we get

f µ (x) min R d f µ = f µ (x) f (x) µ |x (2N +1) x (2N +1) ⇤ | µ = 1 + ⌘ 2N + 1 µ 1 µ, (3.11)
where we used the definition of f and the fact that x (2N +1) = 0.

Let us now upper bound the initial diameter. Remembering that

Lh µ = f µ + d µ in (3.10), we have LD hµ (x ⇤ , x 0 ) = D fµ (x ⇤ , x 0 ) + D dµ (x ⇤ , x 0 ).
by definition of the Bregman divergence. To deal with the first term, we recall that f µ (x ⇤ ) = 0 and write

D fµ (x ⇤ , x 0 ) = f µ (x ⇤ ) f µ (x 0 ) hrf µ (x 0 ), x ⇤ x 0 i = f µ (x 0 ) hrf µ (x 0 ), x ⇤ x 0 i  f (x 0 ) + µ hrf µ (x 0 ), x ⇤ x 0 i = 1 ⌘ + µ hrf µ (x 0 ), x ⇤ x 0 i,
where we used again Lemma 3.9 at x 0 = (0, . . . , 0). Now, Lemma 3.8 applies to x 0 with p = 0 and allows to state that rf µ (x 0 ) 2 E 1 and that krf µ (x 0 )k 1  1. Therefore

|hrf µ (x 0 ), x ⇤ x 0 i| = |rf µ (x 0 ) (1) (x (1) ⇤ x (1) 0 )|  |x (1) ⇤ x (1) 0 | = 1 + ⌘.
Hence we have the following upper bound

D fµ (x ⇤ , x 0 )  1 ⌘ + µ + |hrf µ (x 0 ), x ⇤ x 0 i|  µ.
(3.12)

The second term can be directly computed from Definition (3.7) of d µ , recalling that x

(i) ⇤ 1 µ for i 2 {0 . . . d}, D dµ (x ⇤ , x 0 ) = d µ (x ⇤ ) d µ (x 0 ) hrd µ (x 0 ), x ⇤ x 0 i = d µ (x ⇤ ) = 2N +1 X k=1 h µ 2 (1 + ⌘ k ) 2 + 1 + ⌘ k µ 2 i  (2N + 1) h µ 2 (1 + ⌘) 2 + ⌘ + 1 i .
(3.13) Combining (3.12) and (3.13) gives

LD hµ (x ⇤ , x 0 ) = D fµ (x ⇤ , x 0 ) + D dµ (x ⇤ , x 0 )  µ + (2N + 1) h µ 2 (1 + ⌘) 2 + ⌘ + 1 i  (2N + 1) h µ + µ 2 (1 + ⌘) 2 + ⌘ + 1 i .
This bound, along with (3.11), yields

f µ (x) min R d f µ 1 µ LD hµ (x ⇤ , x 0 ) 2N + 1 • 1 µ 1 + µ + ⌘ + µ 2 (1 + ⌘) 2
whence the desired result.

Since constants µ, ⌘ can be taken arbitrarily small, we now use Proposition 3.10 to show that the bound can be approached to any precision and thus prove our main result.

Theorem 3.12 (Lower complexity bound for B L ). Let N 1, a precision ✏ 2 (0, 1) and let x 0 2 R 2N +1 be a starting point. Then, there exist functions (f, h) 2 B L (R 2N +1 ) such that for any Bregman gradient method A satisfying Definition 3.4 and initialized at x 0 , the output x returned after performing at most N calls to each one of the primal and mirror oracles satisfies

f (x) min R 2N +1 f LD h (x ⇤ , x 0 ) 2N + 1 • (1 ✏).
Proof. Consider a number N of oracle calls and a target precision ✏ 2 (0, 1). Choose the functions f µ , h µ defined respectively in Equations (3.3) and (3.10) on R d with d = 2N + 1. These functions satisfy the conditions in Definition 3.3, since their domain is R d , they are convex, differentiable, and h µ is strongly convex. Moreover, relative smoothness holds because

Lh µ f µ = d µ is convex by construction. Hence (f µ , h µ ) 2 B L (R d ).
Because the class of problems B L (R d ) is invariant by translation, we can assume without loss of generality that the algorithm is initialized at x 0 = (0, . . . 0). Recall that the only conditions our analysis imposed on the parameters ⌘, µ are that µ 2 (0, 1) and ⌘ > 4µd 2 .

Let us then choose ⌘ = ✏/4 and µ = ⌘/(5d 2 ) = ✏/(20d 2 ). Under these conditions, Proposition 3.11 applies and gives that for any point x returned by a Bregman gradient algorithm that is initialized at x 0 and which performs at most N calls to each oracle we have

f µ (x) min R 2N +1 f µ LD hµ (x ⇤ , x 0 ) 2N + 1 • 1 µ 1 + µ + ⌘ + µ 2 (1 + ⌘) 2 .
The last term can be bounded from below, using our choice of µ, ⌘, and the fact that ⌘ < 1, as

1 µ 1 + ⌘ + µ + µ 2 (1 + ⌘) 2 1 µ 1 + ⌘ + 3µ = 1 ✏ 20d 2 1 + ✏ 4 + 3✏ 20d 2 1 ✏
yielding the desired result.

Remark 3. One can refine the result above in the case where the primal and mirror oracles are not used the same number of times. Indeed, if the primal oracles are called N 1 times and the mirror oracle is called N 2 times, then the same reasoning shows that the lower bound remains true by replacing 2N with N 1 + N 2 . Our lower bound involves the relative smoothness constant L instead of the step size in (3.2), but it is equivalent (up to a factor 2) when choosing = 1/L, which is actually the best possible step size choice. This shows the optimality of NoLips within the class of Bregman first-order methods (up to a universal constant).

Connection with Conditional Gradient and the `1 setting. The worst-case function used for the lower bound involves the smoothing of an `1 norm. As pointed out by one of the referees, there might be a connection between the hardness of the relatively-smooth setting and the lower bound for smooth minimization on the `1 ball as done in [START_REF] Guzmán | On Lower Complexity Bounds for Large-Scale Smooth Convex Optimization[END_REF]. This lower bound, which is also O(1/k), is used by the authors to prove that the rate of the Conditional Gradient algorithm is near-optimal in this setting. It might be insightful to examine connections between these settings in future works, for example by exploiting duality between Bregman gradient methods and Conditional Gradient, as in [START_REF] Bach | Duality Between Subgradient and Conditional Gradient Methods[END_REF].

Conclusion

We proved that, in the general relatively-smooth setting, the Bregman gradient/NoLips algorithm is optimal within a large class of Bregman first-order methods. The fundamental idea used for proving our lower bound is that of limiting nonsmooth pathological behavior. Indeed, the class of relatively-smooth problem instances is not closed and the worst case functions are reached as (f, h) approach some non-differentiable functions. This idea, along with the corresponding worst-case functions, have been discovered from the solution to a Performance Estimation Problem, which we detail in the next chapter.

Our result shows that additional assumptions on functions f and h are needed in order to prove better bounds or devise faster algorithms than Bregman gradient descent. If the usual properties of L-smoothness and strong convexity are too restrictive and do not hold in many applications, the future challenge is to find weaker assumptions, that define a larger class of functions where improved rates can be obtained. Another possible approach would to enlarge the oracle model and to find algorithms that do not fit in Definition 3.4, for instance by including second-order oracles of h, in the case when h is simple enough.

Introduction

We consider the constrained minimization problem

min x2C f (x)
where f is a continuously differentiable convex function and C is a nonempty closed convex subset of R d . We are interested in the setting where f satisfies a relative smoothness condition with respect to some convex function h [START_REF] Heinz | A Descent Lemma Beyond Lipschitz Gradient Continuity: First-Order Methods Revisited and Applications[END_REF]. The standard method for solving such a problem is the Bregman gradient descent scheme (BGD), which writes

x k+1 = argmin u2C f (x k ) + hrf (x k ), u x k i + 1 D h (u, x k ), (BGD)
This simpler problem can then be solved numerically with semidefinite optimization packages.

We showcase this approach on several examples, including the proof that (4.1) is the exact worst-case complexity of BGD, a new result on the convergence rate of the stationarity measure D h (x k , x k+1 ) for BGD, and the analysis of the inertial Bregman method from [START_REF] Auslender | Interior Gradient and Proximal Methods for Convex and Conic Optimization[END_REF]. Finally, we show how numerically generated worst-case functions were used to infer the lower bound for general Bregman methods from Chapter 3.

Outline. This chapter is organized as follows. In Section 4.3, we describe the methodology and theoretical guarantees of performance estimation problems for Bregman methods. In Section 4.4, we provide numerical experiments and applications.

Notation. We use C to denote the closure of a set C, int C for its interior and @C for its boundary. S n denotes the set of symmetric matrices of size n ⇥ n. If (P) is an optimization problem, then val(P) stands for its (possibly infinite) value.

Subscripts on a vector denote the iteration counter, while a superscript such as x (i) denotes the i-th coordinate. The set I = {0, 1, . . . N, ⇤} is often used to index the first d iterates of an optimization algorithm as well as the optimal point:

{x i } i2I = {x 0 , x 1 , . . . , x N , x ⇤ }.
We use the standard notation h•, •i for the Euclidean inner product, and k • k for the corresponding Euclidean norm. For a vector x 2 R d , we write kxk 1 = max i=1...n |x (i) | for its `1 norm.

Problem setup

In this section, we briefly recall the blanket assumptions and definitions for relatively-smooth optimization. For a general introduction, we refer the reader to Chapter 1.

Definition 4.1 (Kernel function). A function

h : R d ! R [ {+1} is called a kernel function on C if (i) h is closed convex proper (c.c.p.),
(ii) h is continuously differentiable and strictly convex on int C, (iii) the Bregman gradient iterates are well-posed, i.e., for every p 2 R d , the problem min u2C hp, ui + h(u) has a unique minimizer, which belongs to int C. Definition 4.2 (Relative smoothness). Let h be a kernel function on C, and f a function such that dom h ⇢ dom f . We say that f is smooth relative to h if it is differentiable on int C and if there exists a constant L > 0 such that Lh f is convex on int C.

We recall the notation for a general relatively-smooth problem instance. Note that such an instance is composed of a couple of functions (f, h), as we seek guarantees that are independent of the chosen kernel. Definition 4.3. We say that the couple of functions (f, h) is a relatively-smooth instance, and write

(f, h) 2 B L (C) if (i) h is a kernel function on C, (ii) f : R d ! R [ {+1} is a closed convex proper function, (iii) f is L-smooth relative to h on C.

Worst-case scenarios of Bregman gradient methods through optimization

Formulation of the performance estimation problem

We formulate the task of finding the worst-case performance of (BGD) as an optimization problem. While we focus on the analysis of the vanilla BGD algorithm for ease of presentation, the same ideas are directly applicable to other Bregman-type algorithms like the inertial variant of [START_REF] Auslender | Interior Gradient and Proximal Methods for Convex and Conic Optimization[END_REF], as we illustrate in the sequel.

For simplicity, we first focus on the unconstrained setting where functions have full domain, i.e., C = R d for some d 1. In this setting, the set B L (R d ) can be rewritten as

B L (R d ) = 8 > > > > < > > > > : f is convex and differentiable, h is strictly convex and differentiable, (f, h) : R d ! R Lh f is convex, 8p 2 R d , the function u 7 ! hp, ui + h(u)
has a unique minimizer.

> > > > = > > > > ; ,

The general case when C is a convex subset of R d can be treated along the same approach, as we show in the appendix of this chapter. In fact, from the perspective of performance estimation, every problem in B L (C) can be reduced to some problem in B L (R d ) with equivalent convergence rate.

Performance estimation problem. Throughout this section, we fix a number of iterations N 1, a relative smoothness parameter L > 0, and a step size > 0. In the currently known analyses of Bregman gradient descent, worst-case guarantees have the following form

f (x N ) f (x ⇤ )  ✓(N, L, )D h (x ⇤ , x 0 ). (4.2)
For instance, Equation (4.1) states this result with ✓(N, L, 1/L) = L/N . We then naturally seek the smallest ✓(N, L, ) such that the bound (4.2) holds for any couple

(f, h) 2 B L (R d ), that is, solve the optimization problem maximize f (x N ) f (x ⇤ ) /D h (x ⇤ , x 0 ) subject to (f, h) 2 B L (R d ),
x ⇤ is a minimizer of f, x 1 , . . . , x N are generated from x 0 by BGD with step size , (PEP) in the variables f, h, x 0 , . . . , x N , x ⇤ , d. We refer to this problem as a performance estimation problem (PEP). We use the convention 0/0 = 0, so that the objective is well defined when x ⇤ = x 0 . Optimizing over the dimension d to get dimension-free bounds allows for the problem to admit efficient convex reformulations, as we see in the sequel. Importantly, we look for guarantees that are independent of the kernel h, therefore h is part of the optimization variables.

Let us start by simplifying the problem. First, due to the strict convexity of h, the BGD iteration can be equivalently formulated via the first-order optimality conditions

rh(x i+1 ) = rh(x i ) rf (x i ) 8i 2 {0 . . . N 1},
and, since the domain is R d , the condition that x ⇤ minimizes f reduces to requiring rf (x ⇤ ) = 0.

Second, the problem is homogeneous in (f, h) (i.e., from a feasible couple (f, h), take any constant c > 0 and observe that the couple (cf, ch) is also feasible with the same objective value), hence optimizing the objective function f (x N ) f (x ⇤ ) under the additional constraint D h (x ⇤ , x 0 ) = 1 produces the same optimal value as the problem above. Finally, we use the same argument as in [START_REF] Drori | Performance of First-Order Methods for Smooth Convex Minimization: A Novel Approach[END_REF], [START_REF] Taylor | Smooth Strongly Convex Interpolation and Exact Worst-Case Performance of First-Order Methods[END_REF] and observe that the objective of (PEP) and the algorithmic constraints mentioned above depend solely on the values of the first-order oracles of f and h at the points x 0 , . . . , x N , x ⇤ . Denoting I = {0, 1, . . . , N, ⇤} the indices associated with the points involved in the problem, we proceed to write these values as

{(f i , g i )} i2I = f (x i ), rf (x i ) i2I , {(h i , s i )} i2I = { h(x i ), rh(x i ) } i2I .
Using those elements, the iterations of BGD can be expressed as

s i+1 = s i g i
for i 2 {0 . . . N 1}, and the normalization constraint D h (x ⇤ , x 0 ) = 1 becomes

h ⇤ h 0 hs 0 , x ⇤ x 0 i = 1.
Using those discrete representations of f and h, we can reformulate (PEP) as

maximize f N f ⇤ subject to f i = f (x i ), g i = rf (x i ), h i = h(x i ), s i = rh(x i ), for all i 2 I and some (f, h) 2 B L (R d ), g ⇤ = 0, s i+1 = s i g i for i 2 {1 . . . N 1}, h ⇤ h 0 hs 0 , x ⇤ x 0 i = 1, in the variables d, {(x i , f i , g i , h i , s i )} i2I .
The equivalence with the initial problem is guaranteed by the first two constraints which are called the interpolation conditions.

It turns out that interpolation conditions for the class B L (R d ) are delicate to establish, due to assumptions on h. Fortunately, there exist two classes B L (R d ) and B L (R d ) for which they can be derived. The first class is a restriction of B L (R d ) where f and Lh f are both assumed to be strictly convex:

B L (R d ) = B L (R d ) \ {(f, h) : R d ! R | f and Lh f are strictly convex},
whereas the second class consists in considering a relaxation with possibly nonsmooth functions:

B L (R d ) = {(f, h) : R d ! R | f and Lh f are convex}.
The following inclusions then directly hold

B L (R d ) ⇢ B L (R d ) ⇢ B L (R d ).
With theses classes, we can now define two easier problems. The first one is a restriction of (PEP) defined on the class B L (R d ), under the additional constraint that all iterates are distinct:

maximize f N f ⇤ subject to f i = f (x i ), g i = rf (x i ), h i = h(x i ), s i = rh(x i ), for all i 2 I and some (f, h) 2 B L (R d ), g ⇤ = 0, s i+1 = s i g i for i 2 {1 . . . N 1}, h ⇤ h 0 hs 0 , x ⇤ x 0 i = 1, x i 6 = x j for i 6 = j 2 I, (PEP) in the variables d, {(x i , f i , g i , h i , s i )} i2I .
The second problem is a relaxation of (PEP), where (f, h) 2 B L (R d ) are possibly nonsmooth and g i , s i are thus subgradients:

maximize f N f ⇤ subject to f i = f (x i ), g i 2 @f (x i ), h i = h(x i ), s i 2 @h(x i ), Ls i g i 2 @(Lh f )(x i ) for all i 2 I and some (f, h) 2 B L (R d ), g ⇤ = 0, s i+1 = s i g i for i 2 {1 . . . N 1}, h ⇤ h 0 hs 0 , x ⇤ x 0 i = 1, (PEP)
in the variables d, {(x i , f i , g i , h i , s i )} i2I . We added the technical constraints

Ls i g i 2 @(Lh f )(x i ),
which are redundant for differentiable functions; but that are necessary in order to establish interpolation conditions in the nonsmooth case. Because of the inclusions between the feasible sets of these problems, we naturally have val(PEP)  val(PEP)  val(PEP).

We prove in the sequel that (PEP) can be solved via a semidefinite program and that val(PEP) = val(PEP) (Theorem 4.8), allowing to reach our claims.

Note that the relaxed problem (PEP) does not correspond to any practical algorithm, as BGD is not properly defined for nonsmooth functions h. However, we see in the sequel that feasible points of this problem correspond to accumulation points of (PEP). In other words, instances of BGD can get arbitrarily close to pathological nonsmooth functions whose behaviors are captured by (PEP).

In the following sections, we show that problems (PEP) and (PEP) can be cast as semidefinite programs (SDP) [START_REF] Vandenberghe | Semidefinite Programming[END_REF] and solved numerically using standard packages [START_REF] Lofberg | YALMIP : A Toolbox for Modeling and Optimization in MATLAB[END_REF][START_REF] Aps Mosek | The MOSEK optimization toolbox for MATLAB manual[END_REF]. The main ingredient consists in showing that interpolation constraints can actually be expressed using quadratic inequalities, as detailed in the next section.

Interpolation involving differentiability and strict convexity

In this section, we show how to reformulate interpolation constraints for (PEP) and (PEP) as quadratic inequalities. We start by recalling interpolation conditions for the class of L-smooth and µ-strongly convex functions.

Theorem 4.4 (Smooth strongly convex interpolation, [START_REF] Taylor | Smooth Strongly Convex Interpolation and Exact Worst-Case Performance of First-Order Methods[END_REF]). Let I be a finite index set,

{(x i , f i , g i )} i2I 2 (R d ⇥ R ⇥ R d ) |I| and 0  µ  L  +1.
The following statements are equivalent:

(i) There exists a proper closed convex function f : R d ! R [ {+1} such that f is µ-strongly convex, has a L-Lipschitz continuous gradient and

f i = f (x i ), g i 2 @f (x i ) 8i 2 I.
(ii) For every i, j 2 I we have

f i f j hg j , x i x j i 1 2L kg i g j k 2 + µ 2(1 µ/L) kx i x j 1 L (g i g j )k 2 .
In particular, when L = +1 (meaning that we require no smoothness) and µ = 0, those conditions reduce to the simpler convex interpolation conditions, reminiscent of subgradient inequalities:

f i f j hg j , x i x j i 0. (4.3)
In our setting, we want to avoid working with smoothness and strong convexity, so we provide interpolation conditions for the class of differentiable strictly convex functions. Proposition 4.5 (Differentiable and strictly convex interpolation). Let I be a finite index set and

{(x i , f i , g i )} i2I 2 (R d ⇥ R ⇥ R d ) |I| .
The following statements are equivalent:

(i) There exists a convex function f : R d ! R such that f is differentiable, strictly convex and

f i = f (x i ), g i = rf (x i ) 8i 2 I.
(ii) For every i, j 2 I we have ⇢ f i f j hg j , x i x j i > 0 if x i 6 = x j , f i = f j and g i = g j otherwise.

(4.4)

Proof. (i) =) (ii). Assume that (i) holds, and choose such a function f . The first inequality of (4.4) follows from strict convexity of f , and the second line is a consequence of the fact that a differentiable convex function has a unique subgradient at each point [Rockafellar, 1970, Thm 25.1].

(ii) =) (i). Assume (ii). If for all i, j 2 I, we have g i = g j and x i = x j , then there is only one point and one subgradient to be interpolated, and the result follows immediatly from considering a well-chosen definite quadratic function. In the other case, define

⌫ = min i,j2I x i 6 =x j f i f j hg j , x i x j i.
Because of (4.4) and the finiteness of I, we have that ⌫ > 0. Now, define r as r = max i,j2I kg i g j k 2 + kx i x j k 2 so that r > 0. Condition (4.4) implies that for all i, j 2 I we have

f i f j hg j , x i x j i ⌫ r kg i g j k 2 + kx i x j k 2 . (4.5)
Indeed, if x i 6 = x j , this follows from the definition of ⌫ and r. If x i = x j both sides of the inequality are 0 because of the second line in (4.4). Let us choose two constants

0 < µ < L < +1 such that 1 L µ  ⌫ r , µ 1 µ/L  ⌫ r ,
which is possible as it suffices to take L large enough and µ small enough. We now proceed to show that the interpolation conditions of Theorem 4.4 hold with the constants µ, L defined above. Using the inequality ku vk 2  2kuk 2 + 2kvk 2 and (4.5), we get that for all i, j,

1 2L kg i g j k 2 + µ 2(1 µ/L) kx i x j 1 L (g i g j )k 2  ✓ 1 2L + µ L(L µ) ◆ kg i g j k 2 + µ 1 µ/L kx i x j k 2  ✓ 1 L + µ L(L µ) ◆ kg i g j k 2 + µ 1 µ/L kx i x j k 2 = 1 L µ kg i g j k 2 + µ 1 µ/L kx i x j k 2  ⌫ r kg i g j k 2 + ⌫ r kx i x j k 2  f i f j hg j , x i x j i.
Under those conditions, Theorem 4.4 states that there exists a convex function f that interpolates {(x i , f i , g i )} i2I which is µ-strongly convex and has L-Lipschitz continuous gradients. A fortiori, since µ > 0 and L < 1, f is differentiable and strictly convex. Finally, f is finite on R d since it is L-smooth.

Using these results, we can now formulate interpolation conditions for the problems (PEP) and (PEP) involving the classes B L (R d ) and B L (R d ) that were defined above.

Corollary 4.6 (Interpolation conditions for (PEP)). Let I be a finite index set and

{(x i , f i , g i , h i , s i )} i2I 2 (R d ⇥ R ⇥ R d ⇥ R ⇥ R d ) |I| .
The following statements are equivalent.

(i) There exist functions (f, h) 2 B L (R d ) such that

f i = f (x i ), g i 2 @f (x i ), h i = h(x i ), s i 2 @h(x i ), Ls i g i 2 @(Lh f )(x i ).
(ii) For all i, j 2 I such that i 6 = j, we have

f i f j hg j , x i x j i 0, (Lh i f i ) (Lh j f j ) hLs j g j , x i x j i 0. (4.6)
Proof. (i) =) (ii) follows immediately from the definition of a subgradient applied to convex functions f and Lh f . Assume that (ii) holds. By the specialization of (4.3) in Theorem 4.4, conditions (ii) imply that there exist two convex functions f, d : R d ! R such that

f i = f (x i ), g i 2 @f (x i ), Lh i f i = d(x i ), Ls i g i 2 @d(x i ).
Defining the convex function h = (f + d)/L, we have that d = Lh f , hence

Ls i g i 2 @(Lh f )(x i )
for i 2 I. We also get

h i = h(x i ), s i 2 @h(x i ),
where we used the fact that Ls i 2 @f ( [Rockafellar, 1970, Thm 23.8] for the subdifferential of a sum of convex functions). Hence (i) holds.

x i ) + @d(x i ) ⇢ @(f + d)(x i ) = L@h(x i ) (see
Corollary 4.7 (Interpolation conditions for (PEP)). Let I be a finite index set and

{(x i , f i , g i , h i , s i )} i2I 2 (R d ⇥ R ⇥ R d ⇥ R ⇥ R d ) |I| .
Assume that x i 6 = x j for every i 6 = j 2 I. The following statements are equivalent.

(i) There exist functions (f, h) 2 B L (R d ) such that

f i = f (x i ), g i = rf (x i ), h i = h(x i ), s i = rh(x i ).
(ii) For all i, j 2 I such that i 6 = j we have

f i f j hg j , x i x j i > 0, (Lh i f i ) (Lh j f j ) hLs j g j , x i x j i > 0. (4.7)
Proof. Note that since we assumed x i 6 = x j for every i 6 = j, interpolation conditions of Proposition 4.5 reduce to requiring a strict inequality in (4.4) for every i 6 = j. As before, define d the function d = Lh f.

Then since (f, h) 2 B L (R d ) the functions f and d are differentiable strictly convex, hence the implication (i) =) (ii) follows simply from strict convexity of these functions. Conversely, assume (ii). By using Proposition 4.5 again, we can interpolate differentiable strictly convex functions f and d and recover h with h = (f + d)/L, thus we have naturally Lh f convex. The function h is thus also differentiable and strictly convex. Moreover, it can be seen from the proof of Proposition 4.5 that the interpolating functions can actually be chosen strongly convex, hence with this choice the well-posedness conditions in Definition 4.1 holds, and we can conclude that (f, h) 2 B L (R d ).

Semidefinite reformulations

Now that we established the interpolation conditions for (PEP) and (PEP), we may use them to obtain semidefinite performance estimation formulations as in [START_REF] Drori | Performance of First-Order Methods for Smooth Convex Minimization: A Novel Approach[END_REF], [START_REF] Taylor | Smooth Strongly Convex Interpolation and Exact Worst-Case Performance of First-Order Methods[END_REF]. This is made possible by observing that interpolation conditions (4.6)-(4.7) are quadratic inequalities in the problem variables.

Let {(x i , f i , g i , h i , s i )} i2I be a feasible point of one of the PEPs in dimension d. We write G 2 S 3(N +2) the Gram matrix that contains all dot products between x i , g i , s i for i 2 I, with

G = 0 @ G xx G gx G sx G gx> G gg G gs G sx> G gs> G ss 1 A ⌫ 0
whose size is independent of the dimension d, where the blocks are defined as

G xx ij = hx i , x j i, G gx ij = hg i , x j i, G gs ij = hg i , s j i, G gg ij = hg i , g j i, G sx ij = hs i , x j i, G ss ij = hs i , s j i,
for i, j 2 I. Denote by

F = (f 0 , . . . , f N , f ⇤ ) 2 R N +2 , H = (h 0 , . . . , h N , h ⇤ ) 2 R N +2 ,
the vectors representing the function values of f, h at the iterates. Finally observe that all the constraints of (PEP) and (PEP) can be expressed using only G, F and H. For instance, interpolation conditions (4.6) for B L (R d ) rewrite for all i, j 2 I as

f i f j G gx ji + G gx jj 0, (Lh i f i ) (Lh j f j ) L(G sx ji G sx jj ) + G gx ji G gx jj 0.
This allows us to reformulate the relaxation (PEP) as a semidefinite program, written

maximize f N f ⇤ subject to f i f j G gx ji + G gx jj 0, (Lh i f i ) (Lh j f j ) L(G sx ji G sx jj ) + G gx ji G gx jj 0 for i, j 2 I, G gg ⇤⇤ = 0, G sx i+1,j = G sx ij G gx ij for i 2 {0 . . . N 1}, j 2 I, h ⇤ h 0 G sx 0⇤ + G sx 00 = 1, G ⌫ 0, (sdp-PEP) in the variables G 2 S 3(N +2) and F, H 2 R N +2 .
Any feasible point of (PEP) can be cast into an admissible point of (sdp-PEP) by computing the semidefinite Gram matrix G. Conversely, if G, F, H is an admissible point of (sdp-PEP), then the vectors {(x i , g i , s i )} i2I can be recovered by performing, for instance, a Cholesky decomposition of G. Note that we expressed the algorithmic constraint s i+1 = s i g i only through scalar products with the x i 's in the SDP, since only the projection of the gradients on Span({x i } i2I ) is relevant in the PEPs. Because interpolation conditions from Corollary 4.6 are necessary and sufficient, we conclude that the problems are equivalent, that is val(sdp-PEP) = val(PEP).

The rank of G determines the dimension of the interpolated problem. If we look instead for a solution that has a given dimension d, this would mean imposing a nonconvex rank constraint on G. Our formulation, on the other hand, is convex and finds the best convergence bound that is dimension-independent, which is a usual requirement for large-scale settings.

In the same way, the value of (PEP) can be computed as

maximize f N f ⇤ subject to f i f j G gx ji + G gx jj > 0, (Lh i f i ) (Lh j f j ) L(G sx ji G sx jj ) + G gx ji G gx jj > 0 for i 6 = j 2 I, G gg ⇤⇤ = 0, G sx i+1,j = G sx ij G gx ij for i 2 {0 . . . N 1}, j 2 I, h ⇤ h 0 G sx 0⇤ + G sx 00 = 1, G xx ii + G xx jj 2G xx ij > 0 for i 6 = j 2 I, G ⌫ 0, (sdp-PEP)
in the variables G 2 S 3(N +2) and F, H 2 R N +2 , where we used interpolation conditions for B L (R d ) from Corollary 4.7, since all points {x i } i2I are constrained to be distinct. Therefore, as above we infer that val(sdp-PEP) = val(PEP).

Recalling the hierarchy between the problems, we thus have val(sdp-PEP)  val(PEP)  val(sdp-PEP).

By comparing the two semidefinite programs stated above, one can notice that the only difference is that (sdp-PEP) imposes some inequalities of (sdp-PEP) to be strict. In the next section, we use topological arguments to prove that the values of the two problems are actually equal. In fact, strict inequalities have little meaning in numerical optimization (the value of (sdp-PEP) is actually a supremum and not a maximum); in our experiments, we focus on (sdp-PEP) as solvers usually admit only closed feasible sets.

Tightness of the approach: nonsmooth limit behaviors

We are now ready to prove the main result of this section.

Theorem 4.8. The value of the performance estimation problem (PEP) for BGD is equal to the value of the nonsmooth relaxation (PEP), which can be computed by solving the semidefinite program (sdp-PEP).

Proof. We show that the closure of the feasible set of (sdp-PEP) is the feasible set of (sdp-PEP).

We first need to prove that the strengthened problem (PEP) is feasible, by finding an instance of BGD where f and Lh f are strictly convex and such that all iterates are distinct. It suffices for instance to consider two one-dimensional quadratic functions. Define f, h : R ! R with

f (x) = ↵ 2 x 2 , h(x) = 1 2 x 2 where ↵ = min ✓ 1 2 , L 2
◆ .
Then f is strictly convex and so is Lh

f = L ↵ 2 x 2 since L ↵ L 2 > 0. The optimum is x ⇤ = 0. Choose x 0 = p 2
for which we have D h (x ⇤ , x 0 ) = x 2 0 /2 = 1. Then, BGD is equivalent to gradient descent and the iterates satisfy

x N = (1 ↵) N x 0 .
Since ↵  1/2 < 1, all the iterates are distinct and therefore we constructed a feasible point of (PEP). Let us therefore write (G, F, H) a corresponding feasible point of (sdp-PEP), and (G, F , H) a feasible point of (sdp-PEP). Define the sequence

{(G k , F k , H k )} k 1 as G k = 1 k G + (1 1 k )G, F k = 1 k F + (1 1 k )F , H k = 1 k H + (1 1 k )H.
Then, for every k 1, (G k , F k , H k ) is still a feasible point of (sdp-PEP), because of convexity of the constraints and the fact that adding a strict inequality to a weak inequality gives a strict inequality. Moreover, the sequence converges to the point (G, f, h) when k ! +1.

Hence we proved that for any feasible point of (sdp-PEP), there is a sequence of admissible points of (sdp-PEP) that converge to it. Since the objective is linear in the vector F therefore continuous, we deduce that the two problems have the same value: val(sdp-PEP) = val(sdp-PEP), which means that val(PEP) = val(PEP). As val(PEP) lies in between these two values, we conclude that they are all equal.

Theorem 4.8 states that the value of the original problem (PEP) can be computed numerically with a semidefinite solver applied to (sdp-PEP). The result itself also helps us gain some theoretical insight: it tells us that the worst-case for BGD on B L (R d ) might be reached as (f, h) approach possibly pathological limiting nonsmooth functions in B L (R d ).

Numerical evidence and computer-assisted proofs

We now provide several applications of the performance estimation framework that we developed for Bregman methods.

Finding the exact worst-case convergence rate of BGD

We first start by the most direct application, that is finding the exact worst-case performance of BGD. Theorem 4.8 states that it can be computed by solving the semidefinite program (sdp-PEP). The link to the Matlab implementation is provided in Section 4.5.

To simplify our setting, note that we can assume without loss of generality that the relative smoothness constant L is 1, since we can replace h by a scaled version Lh. Recall that we know from Theorem 3.5, that val(PEP)  1 N .

Table 4.1 shows the result of solving (sdp-PEP) for several values of N up to 100, for a step size = 1/L. We observe that with high precision, val(sdp-PEP) is equal to the theoretical bound 1/( N ). Other values of . One can wonder how the numerical value evolves when one varies the step size . Our experimental observations are as follows:

• For any 2 (0, 1/L], val(PEP) is exactly equal to the theoretical bound 1/( N ).

• For any > 1/L, val(PEP) = +1, hence BGD does not converge in general with these step size values. This suggests that the maximal step size value allowed for BGD is indeed 1/L, unlike the Euclidean setting where gradient descent can be applied with a step size that goes up to 2/L [START_REF] Nesterov | Introductory Lectures on Convex Optimization: A Basic Course[END_REF].

While results above suggest that 1/( N ) is the exact worst-case rate of BGD, they provide only numerical evidence. We can however use them to deduce formal guarantees, both for proving an upper bound and a lower bound.

Upper bound guarantee through duality. As noticed in previous work on PEPs [START_REF] Drori | The Exact Information-Based Complexity of Smooth Convex Minimization[END_REF]Teboulle, 2014, Taylor et al., 2015], solving the dual of (sdp-PEP) can be used to deduce a proof. Indeed, the dual solution gives a combination of the constraints that, when transposed to analytical form, leads to a formal guarantee. This provides the following proof for the O(1/k) convergence rate of Theorem 3.5, which we recall here.

Theorem 4.9 (NoLips convergence rate, recall of Theorem 3.5). Let L > 0, C be a nonempty closed convex subset of R d and (f, h) 2 B L (C) be an relatively-smooth instance. Then the sequence {x k } k 0 generated by Algorithm 5 with constant step size 2 (0, 1/L] satisfies for all k 0

f (x k ) f (u)  D h (u, x 0 ) k for every u 2 dom h.
Proof. The proof relies on the fact that, since Lh f is convex we have that 1 h f is convex for any 2 (0, 1 L ], and only consists in performing the following weighted sum of inequalities:

• convexity of f , between u and x i (i = 0, . . . , k) with weights ⇤,i = 1 k :

f (u) f (x i ) + hrf (x i ), u x i i,
• convexity of f , between x i and x i+1 (i = 0, . . . , k 1) with weights i,i+1 = i k :

f (x i ) f (x i+1 ) + hrf (x i+1 ), x i x i+1 i,
• convexity of 1 h f , between u and x k with weight µ ⇤,k = 1 k :

1 h(u) f (u) 1 h(x k ) f (x k ) + h 1 rh(x k ) rf (x k ), u x k i,
• convexity of 1 h f , between x i+1 and x i (i = 0, . . . , k 1) with weight

µ i+1,i = i+1 k 1 h(x i+1 ) f (x i+1 ) 1 h(x i ) f (x i ) + h 1 rh(x i ) rf (x i ), x i+1 x i i,
• convexity of 1 h f , between x i and x i+1 (i = 0, . . . , k 1) with weight

µ i,i+1 = i k 1 h(x i ) f (x i ) 1 h(x i+1 ) f (x i+1 ) + h 1 rh(x i+1 ) rf (x i+1 ), x i x i+1 i.
The weighted sum is written as

0 k X i=0 ⇤,i [f (x i ) f (u) + hrf (x i ), u x i i] + k 1 X i=0 i,i+1 [f (x i+1 ) f (x i ) + hrf (x i+1 ), x i x i+1 i] + µ ⇤,k ⇥ 1 h(x k ) f (x k ) ( 1 h(u) f (u)) + h 1 rh(x k ) rf (x k ), u x k i ⇤ + k 1 X i=0 µ i+1,i ⇥ 1 h(x i ) f (x i ) ( 1 h(x i+1 ) f (x i+1 )) + h 1 rh(x i ) rf (x i ), x i+1 x i i ⇤ + k 1 X i=0 µ i,i+1 ⇥ 1 h(x i+1 ) f (x i+1 ) ( 1 h(x i ) f (x i )) + h 1 rh(x i+1 ) rf (x i+1 ), x i x i+1 i ⇤ .
By substitution of rh(x i+1 ) = rh(x i ) rf (x i ) (i = 0, . . . , k 1), one can reformulate the weighted sum exactly as (i.e., there is no residual):

0 f (x k ) f (u) h(u) h(x 0 ) hrh(x 0 ),u x 0 i k ,
yielding the desired result. As µ goes to 0, functions f µ , h µ tend to a pathological behavior where all iterates are equal and for which we have exactly f

x 0 = • • • = x 3 x ⇤ f (x) h(x) x 0 x 1 x 2 x 3 x ⇤ f µ (x) h µ (x)
(x N ) f ⇤ = D h (x ⇤ , x 0 )/N .
Lower bound through worst-case functions. As (PEP) computes the exact worst-case performance of BGD, experiments above suggest that 1/( N ) is also a lower bound, meaning that for every ✏ > 0, there exist functions (f, h) 2 B L such that the iterates of BGD satisfy

f (x N ) f ⇤ D h (x ⇤ , x 0 ) N ✏.
We detail here how such functions can be constructed from the solution of (sdp-PEP). The numerical solver allows us to find a maximizer G, F , H (recall that only the relaxed problem has a maximizer as the feasible set is closed), and by factorizing the matrix G as P T P , we can thus recover the corresponding discrete representation {x i , g i , f i , h i , s i } i2I . This discretization can in turn be interpolated to get the corresponding functions (f, h) 2 B L . There are multiple ways to perform this interpolation; see [Taylor et al., 2017, Thm. 1] for a constructive approach.

Recall that since functions (f, h) yield a solution to (PEP), they belong to B L and might thus form a pathological nonsmooth limiting worst-case. They can be approached by valid instances (f µ , h µ ) 2 B L by performing for instance smoothing through Moreau evelopes (as in Chapter 3) and adding a small quadratic to h to make it strictly convex.

There are however many possible maximizers of (sdp-PEP). If we seek a low-dimensional example that may be easily interpretable, we can search for a maximizer such that the Gram matrix G has minimal rank. Using rank minimization heuristics, we were able to find onedimensional worst-case functions. Fix a number of iterations N 1, assume = 1/L = 1 and define f, h : R ! R as

f (x) = |x 1|, h(x) = f (x) + max( Nx, 0),
and set x 0 = 0, x ⇤ = 1. Then clearly (f, h) 2 B L (R). Figure 4.1 shows the functions f, h as well as their smoothed versions (f µ , h µ ) 2 B L (R). Note that the pathological behavior also reflects in the iterates: in the limiting instance, all iterates x 0 , . . . , x N are equal. In the smoothed version, iterates are distinct (since h µ is strictly convex), but they get closer and closer as the smoothing parameter µ goes to 0.

The smoothed function f µ is a Huber function, which is also the worst-case instance for Euclidean gradient descent on L-smooth functions described in [START_REF] Taylor | Smooth Strongly Convex Interpolation and Exact Worst-Case Performance of First-Order Methods[END_REF]. This analysis could be formalized to prove the 1/k lower bound for BGD; however, this bound is just a particular case of the stronger result for general Bregman gradient methods derived in Chapter 3.

Extension to other criteria

In our performance estimation problem, we focused on studying bounds of the form

f (x N ) f ⇤  ✓(N, L, )D h (x ⇤ , x 0 ).
However, we are not limited to this criterion, and different convergence measures might be considered by changing the objective and constraints in (PEP). For instance, another popular criterion is the stationarity measure D h (x k , x k+1 ), which boils down to the squared gradient norm in the unconstrained Euclidean case. By adapting (PEP), we get the following new convergence result for BGD.

Proposition 4.10 (BGD convergence rate, take II). Let L > 0, C be a nonempty closed convex subset of R d and (f, h) 2 B L (C) a relatively-smooth problem instance. Then the sequence {x k } k 0 generated by Bregman gradient descent with constant step size 2 (0, 1/L] satisfies for k 2

min 1ik D h (x i 1 , x i )  2D h (x ⇤ , x 0 ) k(k 1)
for every x ⇤ 2 argmin C f \ dom h.

Proof. In the same way as before, the formal guarantee has been obtained by examining the dual of the corresponding PEP. The proof relies on the fact that 1 h f is convex for any 2 (0, 1 L ], and only consists in performing the following weighted sum of inequalities:

• convexity of f , between x ⇤ and x i (i = 0, . . . , k) with weights ⇤,i = 2 k(k 1) :

f (x ⇤ ) f (x i ) + hrf (x i ), x ⇤ x i i, • optimality of x ⇤ for each x k with weight k,⇤ = 2 k 1 : f (x k ) f (x ⇤ ),
• convexity of 1 h f , between x ⇤ and x k with weight µ ⇤,k = 2 k(k 1) :

1 h(x ⇤ ) f (x ⇤ ) 1 h(x k ) f (x k ) + h 1 rh(x k ) rf (x k ), x ⇤ x k i,
• convexity of 1 h f , between x i+1 and x i (i = 0, . . . , k 1) with weight µ i+1,i = 2 (i+1)

k(k 1) 1 h(x i+1 ) f (x i+1 ) 1 h(x i ) f (x i ) + h 1 rh(x i ) rf (x i ), x i+1 x i i,
• definition of smallest residual among the iterates (i = 1, . . . , k) with weights ⌧ i = 2(i 1) k(k 1) :

h(x i 1 ) h(x i ) hrh(x i ), x i 1 x i i min 1jk {D h (x j 1 , x j )}.
The weighted sum is written as

0 k X i=0 ⇤,i [f (x i ) f (x ⇤ ) + hrf (x i ), x ⇤ x i i] + k,⇤ [f (x ⇤ ) f (x k )] + µ ⇤,k [ 1 h(x k ) f (x k ) ( 1 h(x ⇤ ) f (x ⇤ )) + h 1 rh(x k ) rf (x k ), x ⇤ x k i] + k 1 X i=0 µ i+1,i [ 1 h(x i ) f (x i ) ( 1 h(x i+1 ) f (x i+1 )) + h 1 rh(x i ) rf (x i ), x i+1 x i i] + k X i=1 ⌧ i [ min 1jk {D h (x j 1 , x j )} (h(x i 1 ) h(x i ) hrh(x i ), x i 1 x i i)].
By substitution of rh(x i+1 ) = rh(x i ) rf (x i ) (i = 0, . . . , k 1), one can reformulate the weighted sum exactly as (i.e., there is no residual):

0 min 1jk {D h (x j 1 , x j )} 2 • h(x ⇤ ) h(x 0 ) hrh(x 0 ), x ⇤ x 0 i k(k 1) ,
yielding the desired result.

Beyond BGD: accelerated Bregman algorithms

Our approach is not limited to the vanilla BGD algorithm. For instance, we can also solve the performance estimation problem for the accelerated Bregman algorithm proposed by [START_REF] Auslender | Interior Gradient and Proximal Methods for Convex and Conic Optimization[END_REF], a.k.a. the Improved Interior Gradient Algorithm (IGA). We recall its simplified formulation in Algorithm 6, in the case where there are no affine constraints.

Algorithm 6 Improved Interior Gradient Algorithm (IGA) [START_REF] Auslender | Interior Gradient and Proximal Methods for Convex and Conic Optimization[END_REF] Input: Functions f, h, initial point x 0 2 int dom h, step size .

Set z 0 = x 0 and t 0 = 1.

for k = 0,1,. . . do

y k = (1 1 t k )x k + 1 t k z k z k+1 = argmin {hrf (y k ), u y k i + 1 t k D h (u, z k ) | u 2 R d } x k+1 = (1 1 t k )x k + 1 t k z k+1 t k+1 = (1 + q 1 + 4t 2 k )/2. end for
In the setting where f has L-Lipschitz continuous gradients and h is a -strongly convex kernel function, IGA with step size = / L enjoys the following convergence rate [Auslender Iteration k

(f (x k ) f⇤)/D h (x⇤, x 0 ) min 1ik {D h (x i 1 , x i )}/D h (x⇤, x 0 ) 10 0 10 1 10 2 10 4 10 3 10 2 10 1 10 0 10 1 Iteration k 4 L/( k 2 ) (theory) (f (x k ) f⇤)/(D h (x⇤, x 0 ) + f (x 0 ) f⇤)
Figure 4.2: Numerical worst-case guarantees obtained from PEPs as functions of the iteration counter k (shown in log scale as rates are sublinear). Left: guarantees for Bregman Gradient descent for two different convergence measures. Numerical values confirm exactly the theoretical rates of Theorem 3.5 and Proposition 4.10. Right: guarantees for IGA with no affine constraints (Algorithm 6) under the assumption that h is 1-strongly convex and f is 1-smooth, compared to the theoretical bound from [START_REF] Auslender | Interior Gradient and Proximal Methods for Convex and Conic Optimization[END_REF]. Notice that the theoretical bound is not tight in this case, as it is obtained by making some approximations in the proof.

and Teboulle, 2006, Thm. 5.2]:

f (x N ) f ⇤  4 L N 2 (D h (x ⇤ , x 0 ) + f (x 0 ) f ⇤ ) . (4.8)
Our PEP framework can also be applied to this algorithm, in order to find the smallest value of ✓(N, L, , ) which satisfies

f (x N ) f ⇤  ✓(N, L, , ) (D h (x ⇤ , x 0 ) + f (x 0 ) f ⇤ )
for every instance of IGA with the supplementary assumptions made above. In this case, we use the standard interpolation conditions of Theorem 4.4 for L-smooth and strongly convex functions. Results are shown in Figure 4.2. The exact numerical worst-case performance of IGA is slightly below the theoretical bound above, since the proof in [START_REF] Auslender | Interior Gradient and Proximal Methods for Convex and Conic Optimization[END_REF] makes some approximations.

IGA in the general relatively-smooth case: failure of acceleration. We pointed out in Section 1.2 that the setting in which f is L-smooth and h is -strongly convex is a particular case of relative smoothness with constant L = L/ . The natural question that was also raised in [Teboulle, 2018, Section 6] is therefore: does IGA converge for the general class B L ? Solving the corresponding PEP yields the following results. For Algorithm 6 with the setting that (f, h) 2 B L and several choices of step size in (0, 1/L], the solver states that the value of the corresponding performance estimation problem is unbounded, i.e., there does not exist any ✓ such that the bound (4.8) holds for every instance (f, h) 2 B L .

As suggested by the anonymous reviewers, one could legitimately wonder whether there exist other sequences {t k } k 0 with t k > 1, perhaps less aggressive than the one in Algorithm 6, such that the method converges (note that choosing t k = 1 for all k would yield the standard BGD scheme). After solving the PEP with several choices of such sequences and observing that it is unbounded, we formulate the following conjecture: for any sequence {t k } k 0 , in IGA, such that t k 0 > 1 for some k 0 , it is not possible to bound f (x N ) f ⇤ in general. Of course, this constitutes numerical evidence and not a formal proof. The conjecture could be proved by constructing worst-case functions in the same spirit as in Chapter 3, with some pathological lack of smoothness that would cause the iterates to diverge when taking a step size larger than 1/L.

These experiments lead us to believe that accelerated Bregman methods with non-adaptive coefficients fail to converge in the general relatively-smooth setting.

From worst-case functions for BGD to a lower bound for general Bregman methods

We briefly explain how, with the PEP methodology, the worst case functions from Chapter 3 were discovered.

We described in the previous section how a one-dimensional worst-case instance (f, h) for BGD was discovered from low-rank solutions of (sdp-PEP). However, this instance may not be difficult enough for a more generic Bregman algorithm that can use abritrary linear combinations of gradients, and thus cannot be used to prove a general lower bound.

Our objective now is to find worst-case instances that are difficult for any Bregman gradient algorithm. A desirable property would be that these instances allow to explore only one dimension per oracle call, so that the function hides information in the unexplored dimensions. This is similar in spirit to the so-called "worst function in the world" of [START_REF] Nesterov | Introductory Lectures on Convex Optimization: A Basic Course[END_REF]. In order to achieve this goal, we propose to search for functions f for which all gradients rf (x i ) are orthogonal, guaranteeing that one new dimension is explored at each step. Note that a similar approach has been used in some previous work on PEPs to find lower bounds or optimal methods e.g., in [START_REF] Drori | The Exact Information-Based Complexity of Smooth Convex Minimization[END_REF], [START_REF] Drori | Efficient First-Order Methods for Convex Minimization: a Constructive Approach[END_REF]. This amounts to adding some orthogonality constraints to (PEP) and solving

maximize f (x N ) f (x ⇤ ) /D h (x ⇤ , x 0 ) subject to (f, h) 2 B L (R d ),
x ⇤ is a minimizer of f, x 1 , . . . , x N are generated from x 0 by BGD with step size , hrf (x i ), rf (x j )i = 0 for i 6 = j 2 I,

(PEP-orth)
in the variables f, h, x 0 , . . . , x N , x ⇤ , d. In the same spirit as before, we were able to find a dimension-N solution of (PEP-orth). This allows us to interpolate the following worst-case pathological instance:

f (x) = kx (1, . . . , 1)k 1 , h(x) = f (x) + N X i=2 max( x (i) , 0).
Again, these are nonsmooth functions and, as such, they do not form valid instances for BGD. However, they can be approached by a sequence of such functions, for instance by applying smoothing with the Moreau evelope, and adding a small quadratic term to make h strictly convex. Along with a few tweaks, this is how we found the worst-case instance that was used to prove the general lower bound for B L in Chapter 3.

Conclusion

In this chapter, we showed how the worst-case complexity of Bregman gradient methods with generic kernel can be computed with the aid of performance estimation problems. To this end, we established interpolation condition for the class of differentiable and strictly convex functions. As we showcased on several examples, PEPs allow to conjecture and prove convergence rates in various situations.

In the setting of Bregman methods, a fundamental concept arising from this work is that of limiting nonsmooth pathological behavior. When looking for worst-case guarantees over a class of functions that is open such as the class of differentiable convex functions, the performance estimation problem is a supremum and the worst-case maximizing sequence might approach some function that is not in this class, e.g., one that is nonsmooth in our case. This idea, observed by analyzing the equivalence between (PEP) and the nonsmooth relaxation (PEP), was used in the proof of the lower bound in Chapter 3. Moreover, the worst-case sequence of functions was directly inspired by examining particular solutions of (PEP).

Code. Experiments have been run in MATLAB, using the semidefinite solver MOSEK Mosek [2019] as well as the modeling toolbox YALMIP [START_REF] Lofberg | YALMIP : A Toolbox for Modeling and Optimization in MATLAB[END_REF]. The support for Bregman methods has been added to the Performance Estimation Toolbox (PESTO, [START_REF] Taylor | Performance Estimation Toolbox (PESTO): Automated Worst-Case Analysis of First-Order Optimization Methods[END_REF]) for which we provide some examples. The code can be downloaded from the repository https://github.com/RaduAlexandruDragomir/BregmanPerformanceEstimation Appendix: extension of performance estimation to the case when

C is a general closed convex subset of R d
For simplicity of the presentation, we left out in Section 4.3 the case when the domain C is a proper subset of R d . We show in this section that it actually corresponds to the same minimization problem (sdp-PEP).

Let us formulate the performance estimation problem for Bregman gradient descent in the general case. Recall that we denote B L the union of B L (C) for all closed convex subsets of R d and for every d 1. The performance estimation problem writes

maximize f (x N ) f (x ⇤ ) /D h (x ⇤ , x 0 ) subject to (f, h) 2 B L ,
x ⇤ is a minimizer of f on dom h such that x 2 dom h, x 1 , . . . , x N are generated from x 0 by BGD with step size , (PEP-C) in the variables f, h, x 0 , . . . , x N , x ⇤ , d. Now, as (PEP-C) is a problem that includes (PEP) in the special case where C = R d , its value is larger:

val(PEP)  val(PEP-C)
Let us show that val(PEP-C) is upper bounded by the same relaxation val(PEP), which allows to conclude that the values are equal. We recall that the problem (PEP) can be written, using interpolation conditions of Corollary 4.6, as

maximize f N f ⇤ subject to f i f j hg j , x i x j i 0, (Lh i f i ) (Lh j f j ) hLs j g j , x i x j i 0 for i, j 2 I, g ⇤ = 0, s i+1 = s i g i for i 2 {1, . . . , N 1}, h ⇤ h 0 hs 0 , x ⇤ x 0 i = 1,
in the variables n, {(x i , f i , g i , h i , s i )} i2I . We show that every admissible point of (PEP-C) can be cast into an admissible point of (sdp-PEP). This actually amounts to show that, from the point of view of performance estimation, an instance (f, h) 2 B L (C) is actually equivalent to some instance in B L (R d ).

Let f, h, x 0 , . . . , x N , x ⇤ , d be a feasible point of (PEP-C). We distinguish two cases.

Case 1: x ⇤ 2 int dom h. This is the simplest case, as the necessary conditions are the same as in the situation where C = R d . Indeed, then we have x 0 , . . . , x N , x ⇤ 2 int dom h, since x 0 is constrained to be in the interior and the next iterates are in int dom h by definition of a kernel function. Since f and h are differentiable on int dom h, convexity of f and Lh f imply that the first two constraints of (PEP) hold for all i, j 2 I. Finally, g ⇤ = 0 follows from the fact that x ⇤ minimizes f and that it lies on the interior of the domain. Hence the discrete representation satisfies the constraints of (sdp-PEP).

Case 2: x ⇤ 2 @dom h. In this case, f and h are not necessarily differentiable at x ⇤ , but are still differentiable still at x 0 , . . . , x N for the same reasons. But we can still, with a small modification at x ⇤ , derive a discrete representation that fits the constraints of (PEP) and whose objective is the same. Indeed, define

(g i , f i , s i , h i ) = (rf (x i ), f(x i ), rh(x i ), h(x i )) for i = 0, . . . , N, (g ⇤ , f ⇤ , s ⇤ , h ⇤ ) = (0, f(x ⇤ ), v, h(x ⇤ )) ,
where v 2 R d is a vector that are specified later. Then, for i 2 I and j 2 {0 . . . N}, convexity of f and Lh f imply that the constraints

f i f j hg j , x i x j i 0 (Lh i f i ) (Lh j f j ) hLs j g j , x i x j i 0
hold. It remains to verify them for i 2 {0 . . . N} and j = ⇤. The first one holds because x ⇤ minimizes f on dom h, so with g ⇤ = 0 we have f i f ⇤ 0. We now show that the second one is satisfied, i.e., that we can choose v 2 R d so that

(Lh i f i ) (Lh ⇤ f ⇤ ) hLv, x i x ⇤ i 0 8i 2 {0 . . . N}.
To this extent, we use the fact that x ⇤ 2 @dom h and that x i 2 int dom h for i = 0 . . . N. This means that {x ⇤ } \ int dom h = ;, and therefore by the hyperplane separation theorem [Rockafellar, 1970, Thm 11.3], there exists a hyperplane that separates the convex sets {x ⇤ } and int dom h properly, meaning that there exists a vector u 2 R d such that

hx i x ⇤ , ui < 0 8i 2 {0, . . . , N}. Set ↵ = min i=0...N (Lh i f i ) (Lh ⇤ f ⇤ ), = min i=0,...,N hx i x ⇤ , ui > 0,
where > 0 because of the separation result. Choose s ⇤ = v as v = |↵| L u. Then we have

(Lh i f i ) (Lh ⇤ f ⇤ ) hLs ⇤ , x i x ⇤ i ↵ + L |↵| L ↵ + |↵| 0.
This eventually provides an instance {(x i , g i , f i , h i , s i )} i2I that is admissible for (PEP).

To conclude, we proved that in both cases, an admissible point of (PEP-C) can be turned into an admissible point of (sdp-PEP) with the same objective value. Hence we have val(PEP-C)  val(sdp-PEP).

Recalling that val(PEP)  val(PEP-C) and that val(sdp-PEP) = val(PEP) by Theorem 4.8, we get val(PEP-C) = val(PEP).

In other words, solving the performance estimation problem (PEP-C) for functions with any closed convex domain is equivalent to solving the performance estimation problem (PEP) restricted to functions that have full domain.

Chapter 5

Computer-Aided Analyses of Entropic-Smooth Minimization Methods

Chapter Abstract

In the previous chapter, we showed that the problem of finding the worstcase behavior of Bregman methods for generic kernel could be formulated as a semidefinite optimization problem. This lead to a lower bound proving the impossibility of acceleration, relying on a pathological worst-case kernel, and demonstrated the necessity of making additional regularity assumptions in order to devise faster algorithms. In this work, we propose to tackle a more restricted setting, by focusing on the entropic kernel. We show that finding the worst-case behavior of Bregman methods on functions that are convex and smooth relative to the entropy can be formulated as a finite dimensional convex program. This problem involves a new set called the Kullback-Leibler cone with log-linear constraints, for which no solver is currently available. However, we manage to solve small instances with heuristic methods and provide some conjectures based on preliminary numerical results.

Collaboration: the content presented in this chapter has not yet been published. Part of this work has been done in collaboration with Dmitrii Ostrovskii.

Introduction

As in the previous chapter, we consider the minimization problem

min x2C f (x)
where C is a convex subset of R d and f a continuously differentiable convex function which is smooth relative to some convex kernel h. The standard method for solving such problem is the Bregman gradient descent scheme (BGD), which writes

x k+1 = argmin u2C f (x k ) + hrf (x k ), u x k i + 1 D h (u, x k ), (BGD)
where

D h (x, y) = h(x) h(y) hrh(y), x yi
is the Bregman divergence induced by h. An important question, among others, is to determine if there exists an algorithm that achieves a better convergence rate than (BGD) on the same class of functions.

In the previous chapters, we gave a rather negative answer to this question for general kernel h. Indeed, we proved in Chapter 3 a lower bound implying that BGD is optimal among a large class of Bregman-type methods. The worst-case instance involved in this lower bound was discovered numerically by solving a performance estimation problem (PEP), which we described in Chapter 4. PEPs allow to formulate the problem of finding the worst-case complexity of gradient methods as numerical optimization problems. They form a great tool for conjecturing and proving new results in first-order optimization [START_REF] Drori | Performance of First-Order Methods for Smooth Convex Minimization: A Novel Approach[END_REF][START_REF] Taylor | Exact Worst-Case Performance of First-Order Methods for Composite Convex Optimization[END_REF][START_REF] Kim | Optimized First-Order Methods for Smooth Convex Minimization[END_REF][START_REF] Drori | Efficient First-Order Methods for Convex Minimization: a Constructive Approach[END_REF][START_REF] Taylor | Stochastic first-order methods: non-asymptotic and computer-aided analyses via potential functions[END_REF].

A critical limitation of the results from the previous chapters is that we consider generic kernels, that is, we search for the worst possible couple (f, h) such that f is convex and smooth relative to h. By doing so, the worst case kernels that we find are pathological and somewhat unrealistic. For instance, the proof of the lower bound in Chapter 3 involves a sequence of functions {h µ } µ 0 which approaches a nondifferentiable convex function as µ ! 0. In practical applications of relative smoothness, h is more regular, as is it usually a C 1 function. This is why, in this chapter, we propose to look at the particular case of the discrete entropy

h e (x) = d X i=1 x (i) log x (i) x (i) ,
on the set C = R d + , for which the corresponding Bregman divergence is the well-known Kullback-Leibler divergence, and the update of (BGD) writes

x k+1 = x k exp ( rf (x k )) .
Therefore, we consider the class of functions f that are convex and smooth relative to the entropy, which we call the class of entropic-smooth convex functions. We focus on the entropic kernel for two main reasons. First, it is one of the most common example of non-Euclidean kernels, with applications such as optimal transport and nonnegative entropy regression (see Chapter 1). Secondly, the entropy has some favorable properties which are crucial for our analysis, including the fact that the corresponding Bregman divergence is jointly convex [START_REF] Heinz | Joint and Separate Convexity of the Bregman Distance[END_REF]. Moreover, studying the entropic geometry might also help gain insight about other kernels with similar structure on the nonnegative orthant, such as the log-barrier.

In this chapter, we consider only the unnormalized case where C = R d + ; the setting where C is the unit simplex {x 2 R d + :

P d i=1 x (i) = 1}
is left for future work as it poses additional difficulties for performance estimation.

Contributions and outline. We propose to apply the performance estimation framework of [START_REF] Drori | Performance of First-Order Methods for Smooth Convex Minimization: A Novel Approach[END_REF] for finding the worst-case behavior of Bregman methods on entropicsmooth convex functions.

We start by recalling in Section 5.2 the basic setup and defining the performance estimation problem (PEP). In Section 5.3, we establish the interpolation conditions which are needed for formulating the PEP as an equivalent finite-dimensional optimization problem. To do so, we generalize the Euclidean interpolation conditions of [START_REF] Taylor | Smooth Strongly Convex Interpolation and Exact Worst-Case Performance of First-Order Methods[END_REF] to the entropic-smooth setting. We rely on the smoothing properties of the Bregman-Moreau envelope [START_REF] Kan | The Moreau envelope function and proximal mapping in the sense of the Bregman distance[END_REF][START_REF] Heinz | Regularizing with Bregman-Moreau envelopes[END_REF][START_REF] Laude | Bregman Proximal Mappings and Bregman-Moreau Envelopes under Relative Prox-Regularity[END_REF], which is the Bregman counterpart of the standard Moreau envelope [START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF]. Although part of our analysis can be done for a more general class of kernels, we need additional restrictive assumptions to reach our claims (such as joint convexity of the Bregman divergence), which are satisfied by the entropy.

Then, in Section 5.4, we show that the PEP can be written as a convex problem on the set of matrices of pairwise Kullback-Leibler divergences, which we call the Kullback-Leibler cone, with additional log-linear constraints. Unfortunately, we are not aware of an explicit description of such set, or any solvers for this type of problem. We are only able to describe the set in the simple case without log-linear constraints: in that case, we show that the set is trivial, as any nonnegative matrix with null diagonal can be approached by a matrix of pairwise Kullback-Leibler divergences.

Finally, we show in Section 5.5 some preliminary numerical experiments, made by solving the PEP with basic heuristics. This leads to the following pessimistic conjecture: among the class of linear-span Bregman methods with fixed coefficients, only vanilla BGD is guaranteed to converge in the worst-case.

Notation. Throughout this chapter, we use the notation h for a generic convex kernel function, while h e denotes the entropic kernel

h e (x) = d X i=1 x (i) log x (i) x (i)
defined for x 2 R d + , with the convention that 0 log 0 = 0. The corresponding Bregman divergence is the Kullback-Leibler divergence

D he (x, y) = d X i=1 x (i) log x (i) y (i) x (i) + y (i) defined for x 2 R d + , y 2 R d ++ .
We write F h L the class of convex functions that are L-smooth relative to h. 1 m denotes the vector (1, . . . , 1) 2 R m . Co(A) is the convex cone spanned by the elements of the set A. For x, y 2 R d , x y denotes the componentwise product of x and y. For a scalar function : R ! R and a vector x 2 R d , (x) denotes the result of applying componentwise to the vector x, i.e., (x) = (x (1) ), . . . , (x (d) ) . Let X 2 R d⇥m a matrix, and denote X = [x 1 ; . . . ; x m ] its column vectors. We define the pairwise divergence matrix D

h [X, X] 2 R m⇥m as (D h [X, X]) ij = D h (x i , x j ), i,j = 1 . . . m.

Problem setup

Relatively-smooth optimization

We briefly recall the basic ingredients and definitions for relatively-smooth optimization problems. For a general introduction, we refer the reader to Chapter 1.

Definition 5.1 (Kernel function). A function h : R d ! R [ {+1} is called a kernel function on C if (i) h is closed convex proper (c.c.p.),
(ii) h is continuously differentiable and strictly convex on int C, (iii) the Bregman gradient iterates are well-posed, i.e., for every p 2 R d , the problem

min u2C hp, ui + h(u)
has a unique minimizer, which belongs to int C.

The entropy h e is a particular case of kernel on

C = R d + .
We also recall the central notion of relative smoothness [START_REF] Heinz | A Descent Lemma Beyond Lipschitz Gradient Continuity: First-Order Methods Revisited and Applications[END_REF].

Definition 5.2 (Relative smoothness). Let h be a kernel function on C, and f a function such that dom h ⇢ dom f . We say that f is smooth relative to h if it is differentiable on int C and if there exists a constant L > 0 such that

Lh f is convex on int C, or equivalently, if f (x)  f (y) + hrf (y), x yi + LD h (x, y) 8x, y 2 int C.
(5.1)

For a given kernel h, we write f 2 F h L if the function f is convex and L-smooth relative to h.

Entropic smoothness. When a function f is smooth relative to the entropy h e , we say that f is entropic-smooth. The canonical example of convex entropic-smooth functions is given by functions of the form et al., 2017]. Such a function arises in nonnegative entropy regression and entropy-regularized optimal transport [START_REF] Chizat | Scaling algorithms for unbalanced optimal transport problems[END_REF][START_REF] Mishchenko | Sinkhorn Algorithm as a Special Case of Stochastic Mirror Descent[END_REF].

f A,b (x) = D he (Ax, b) for some matrix A 2 R p⇥d + and vector b 2 R p ++ . The function f A,b is entropic-smooth with constant L = max j P i A ij [Bauschke

Convex conjugate. Recall that if h is a kernel function, we define its convex conjugate h

⇤ as h ⇤ (y) = sup u2R d hu, yi h(u).
By Definition 5.1 (iii), h ⇤ is finite-valued on R d . Moreover, its gradient satisfies for every u 2 R d [Rockafellar, 1970, Sect. 26]

rh ⇤ (y) = argmax u2R d hu, yi h(u).
We recall the relationship between the Bregman divergence of h and its conjugate; see [Bauschke and Borwein, 1997, Thm. 3.7] for the proof.

Lemma 5.3 (Bregman divergence of conjugate). Let h be a kernel function on C. Then for every x, y 2 int C, we have

D h (x, y) = D h ⇤ (rh(y), rh(x)) .
The conjugate of the entropy function is the exponential kernel: for u 2 R d we have

h ⇤ e (u) = d X i=1 exp(u (i) ).
Let us check that Lemma 5.3 holds for this example. Let x, y 2 R d ++ and define

u = rh e (x) = log(x), v = rh e (y) = log(y).
Then we have

D he (x, y) = d X i=1 x (i) log x (i) y (i) x (i) + y (i) = d X i=1 exp(v (i) ) exp(u (i) ) ⇣ 1 + v (i) u (i) ⌘ = D h ⇤ e (v, u).
(5.2)

Performance estimation problems for Bregman methods with fixed kernel

In this section, we formulate the problem of finding the worst-case behavior of a Bregman firstorder method on the class F h L . The fundamental difference with the approach from Chapter 4 is that we now look for the worst case for a given kernel h. In the sequel, we will particularly focus on the entropy h e .

Let h be a kernel function, N 1 be a number of iterations, and A a Bregman firstorder method. We wish to compute the worst-case performance of A on the class F h L after N iterations. Since performance guarantees in the Bregman setting have the form

f (x N ) f (x ⇤ )  ✓(N, L)D h (x ⇤ , x 0 ),
for some function ✓, it is natural to seek the worst-case performance as the maximal value of the ratio

f (x N ) f ⇤ /D h (u, x 0 )
over the problem class F h L . To this extent, we propose to solve the minimization problem, for an initial radius R > 0,

maximize f (x N ) f (x ⇤ ) subject to D h (x ⇤ , x 0 ) = R f 2 F h L , x N is the output of algorithm A initialized at x 0 after N iterations, (PEP)
in the variables f, x 0 , x N , x ⇤ . We refer to this problem as a Performance Estimation Problem (PEP). Note that we do not need impose that x ⇤ minimizes f , this will occur naturally in the solutions of (PEP) as to maximize the gap f (x N ) f (x ⇤ ). Although it seems that (PEP) should be solved for every value of R, we will use homogeneity arguments in the sequel to show that it suffices to consider R = 1.

Solving (PEP) seems intractable, as it involves the infinite-dimensional variable f . Following [START_REF] Drori | Performance of First-Order Methods for Smooth Convex Minimization: A Novel Approach[END_REF], [START_REF] Taylor | Smooth Strongly Convex Interpolation and Exact Worst-Case Performance of First-Order Methods[END_REF], we proceed to reformulate it as a finite dimensional problem. Since we assume that the algorithm A is a first-order method and we adopt the black box model, it can gain information on f only through the first-order oracle (f, rf ) at query points x 0 , . . . , x N . Thus we formally write

x N = A ⇣ x 0 , {(f (x i ) , rf (x i ))} N 1 i=0 ⌘
and, introducing the variables

f i = f (x i ), g i = rf (x i ) for i in the index set I = {0, . . . , N, ⇤}, Problem (PEP) is equivalent to maximize f N f ⇤ subject to D h (x ⇤ , x 0 ) = R, f i = f (x i ), g i = rf (x i ) for all i 2 I and some f 2 F h L , x N = A ⇣ x 0 , {(f i , g i )} N 1 i=0 ⌘ , (PEP-disc)
in the variables {(x i , f i , g i )} i2I . We now deal with a finite-dimensional optimization problem, the equivalence with the previous PEP being guaranteed by the so-called interpolation constraints. We introduce the following convenient definition for interpolable sets [START_REF] Taylor | Smooth Strongly Convex Interpolation and Exact Worst-Case Performance of First-Order Methods[END_REF].

Definition 5.4. Let I be a finite index set and

S = {(x i , f i , g i )} i2I 2 (R d ⇥ R ⇥ R d ) |I|
for some dimension d. Let F be a class of convex functions. We say that the set S is F-interpolable if there exists a function f 2 F such that

f i = f (x i ), g i 2 @f (x i ) 8i 2 I.
Hence, the second line in (PEP-disc) amounts to requiring that the set {(x i , f i , g i )} i2I is F h L -interpolable. We now proceed to establish necessary and sufficient conditions in the case of the entropic kernel.

Entropic-smooth convex interpolation

In the Euclidean case, [START_REF] Taylor | Smooth Strongly Convex Interpolation and Exact Worst-Case Performance of First-Order Methods[END_REF] provide sufficient and necessary conditions for a finite set {(x i , f i , g i )} i2I to be interpolated by a L-smooth strongly convex function (Theorem 4.4). Our aim is to generalize this approach and to provide interpolation conditions for the class F h L of functions that are smooth relative to a given kernel h (recall that, contrary to the previous chapter where we considered interpolation conditions with generic kernels, h is now considered to be fixed).

It turns out that this task requires additional assumptions on h, as we show in the sequel. While necessary interpolation conditions can be derived for any kernel, sufficient conditions require joint convexity of the Bregman divergence D h (•, •), which is a rather restrictive assumption. Fortunately, it is valid for the important case of the entropy [START_REF] Heinz | Joint and Separate Convexity of the Bregman Distance[END_REF]. Finally, we show that the entropy satisfies a favorable algebraic property which allows to prove the equivalence between the necessary and the sufficient conditions.

Necessary conditions using Bregman co-coercivity

We start by establishing necessary conditions for interpolability. They are based on the following consequence of relative smoothness, which we call Bregman co-coercivity. It generalizes a known inequality for L-smooth convex functions [Nesterov, 2003, Eqn. (2.1.7)].

Proposition 5.5. Let h be a kernel function on C and f a convex functions such that f is L-smooth relative to h. Then for every x, y 2 int C we have

D f (x, y) LD h ⇤ ✓ rh(x) 1 L rf (x) rf (y) , rh (x)
◆ .

Proof. Let y 2 int C. Consider the function y defined for x 2 int C by y (x) = D f (x, y). Since y is built by adding a linear function to f , it is also convex and L-smooth relative to h. Thus, we can apply the descent inequality (5.1) which writes for u, x 2 C,

y (u)  y (x) + hr y (x), u xi + LD h (u, x) := Q x,y (u) (5.3)
The right-hand side Q x,y (u) is a strictly convex function which is minimized for û such that

1 L r y (x) + rh(û) rh(x) = 0, (5.4)
Thus, using (5.4) to substitute r y (x) in (5.3) yields

y (û)  Q x,y (û) = y (x) + Lhrh(x) rh(û), u xi + LD h (û, x) = y (x) LD h (x, û)
Recalling that 0  x by convexity of f and using Lemma 5.3 we get

0  y (x) LD h (x, û) = D f (x, y) LD h ⇤ (rh(û), rh(x)) = D f (x, y) LD h ⇤ ✓ rh(x) 1 L r y (x), rh (x)
◆ which yields the desired result as r y (x) = rf (x) rf (y).

As a consequence, we deduce necessary conditions for F h L -interpolability. Proposition 5.6 (Necessary interpolation conditions). Let h be a kernel on C. Let I be a finite index set and

S = {(x i , f i , g i )} i2I 2 (int C ⇥ R ⇥ R d ) |I| . If the set S is F h L -interpolable, then f i f j hg j , x i x j i LD h ⇤ ✓ rh(x i ) 1 L (g i g j ), rh(x i ) ◆ 8i, j 2 I.
Proof. This simply amounts to choosing an interpolating function f 2 F h L and applying the Bregman co-coercivity property (Proposition 5.5) at pairs of points (x i , x j ).

Sufficient conditions with Bregman-Moreau envelopes

To establish sufficient conditions, we need a constructive procedure for building a function in

F h L which interpolates the finite set {(x i , f i , g i )} i2I .
In the Euclidean setting, [START_REF] Taylor | Smooth Strongly Convex Interpolation and Exact Worst-Case Performance of First-Order Methods[END_REF] proceed by first building a nonsmooth convex function as a supremum of affine functions, and then smoothing by taking its Moreau envelope [START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF] 1 . Recall that the (Euclidean) Moreau envelope of a convex function f with parameter µ > 0 is given by

env µ f (x) = min u2R d f (u) + 1 2µ ku xk 2 .
(5.5)

If f is closed convex and proper, then env µ f is convex and µ 1 -smooth [Bauschke and Combettes, 2011, Sect. 12.4]. Thus, the Moreau envelope posses a smoothing property, and the level of smoothing is controlled by the coefficient µ.

Bregman-Moreau envelopes. It is natural to wonder if by replacing the squared Euclidean distance with a Bregman divergence in (5.5), relatively-smooth functions can be built based on the same principle. Note that there are two ways to perform this extension, as the Bregman divergence D h is not symmetric in general. In this work, we show that, in order to reach the desired property, we need to consider the right Bregman-Moreau envelope defined as

env h µ f (x) = min u2int C f (u) + 1 µ D h (x, u), (5.6)
Historically, the variant that has been considered first is the left Bregman-Moreau envelope which writes min u2R d f (u) + µ 1 D h (u, x), because it is the operation that naturally appears in Bregman proximal methods [START_REF] Kan | The Moreau envelope function and proximal mapping in the sense of the Bregman distance[END_REF]. The idea of considering the right Bregman-Moreau envelope appeared first in Bauschke [2006], and its smoothing properties were studied in Bauschke et al. [2018]; see also [START_REF] Laude | Bregman Proximal Mappings and Bregman-Moreau Envelopes under Relative Prox-Regularity[END_REF] for a recent extension to nonconvex functions f . The study of the right Bregman-Moreau envelope requires additional structural assumptions on h; in particular, we need to assume that D h (•, •) is jointly convex so that env h µ f is well defined.

1 Actually, the procedure described by [START_REF] Taylor | Smooth Strongly Convex Interpolation and Exact Worst-Case Performance of First-Order Methods[END_REF] does not mention explicitly using a Moreau envelope.

They rather propose to build a smooth convex function by adding a quadratic to the conjugate of a nonsmooth convex function, and taking the conjugate back. One can show that this process is equivalent to applying the Moreau envelope, by exploiting the relationship between conjugacy and inf-convolution [Bauschke and Combettes, 2011, Prop. 13.21].

(ii) The set

Ŝ := {(u i , fi , ĝi )} i2I
is (R d )-interpolable, where (R d ) denotes the set of convex closed proper functions on R d and the set Ŝ is defined as

u i = rh ⇤  rh(x i ) 1 L g i , fi = f i LD h (x i , u i ), ĝi = Lr 2 h [u i ] (x i u i ),
(5.11)

for i 2 I.

(iii) The set S is F h L -interpolable.

Then (i) =) (ii) =) (iii).
Proof. (i) =) (ii). Define {u i , fi , ĝi } as in (5.11). The definition of u i (5.10) implies that L(rh(x i ) rh(u i )) = g i , (5.12) and hence we have fi fj hĝ j , u i u j i

(5.11)

= f i LD h (x i , u i ) f j + LD h (x j , u j ) Lhr 2 h(u j )(x j u j ), u i u j i = f i LD h (x i , u i ) f j + LD h (x j , u j ) Lhr 2 h(u j )(x j u j ), u i u j i hg j , x i x j i + Lhrh(x j ) rh(u j ), x i x j i
(5.12)

= f i f j hg j , x i x j i L ⇣ D h (x i , u i ) D h (x j , u j )
hrh(x j ) rh(u j ), x i x j i hr 2 h(u j )(u j x j ), u i u j i

⌘

(5.9)

= f i f j hg j , x i x j i LD D h (x i , u i ), (x j , u j )

(i) 0
Hence, the set Ŝ = {u i , fi , ĝi } i2I satisfies the interpolation conditions for convex interpolation [Taylor et al., 2017, Thm. 1], and therefore there exists a convex function f : R d ! R that interpolates Ŝ, thus proving (ii).

(ii) =) (iii). Let f be a convex function that interpolates the set Ŝ. Moreover, let us choose f such that f is finite valued and bounded below on R d (this is always possible since the set S is finite). Let µ = 1/L and define the function f as the Bregman-Moreau envelope

f := env h µ f.
Since f is convex and bounded below, Proposition 5.7 ensures that f 2 F h L . Let i 2 I. As f interpolates Ŝ, we have ĝi 2 @ f (u i ) which means that

Lr 2 h(u i )(x i u i ) 2 @ f (u i )
by definition of ĝi (5.11). Hence, we have that u i satisfies the first-order optimality conditions defining the Bregman-Moreau proximal map of f at x i (Proposition 5.8), therefore

u i = prox h µ f (x i ).
As f interpolates Ŝ we have f (u i ) = fi and thus

f (x i ) (5.8) = f (u i ) + LD h (x i , u i ) = fi + LD h (x i , u i ) (5.11) = f i ,
and similarly

rf (x i ) (5.7) = 1 µ ⇣ rh(x i ) rh ⇣ prox h µ f (x i ) ⌘⌘ = L (rh(x i ) rh(u i ))
(5.12)

= g i , which achieves to prove that f interpolates the set S.

Equivalence for the entropic kernel

By examining the necessary conditions from Proposition 5.6 and the sufficient conditions form Proposition 5.9, one can notice that they do not seem to coincide in general. However, in the case of the entropic kernel, some favorable algebra allows to show that they are in fact equivalent.

Lemma 5.10. Let L > 0, h e be the entropic kernel, x, x 0 2 R d ++ , g, g 0 2 R d and define

u = rh ⇤ e (rh e (x) 1 L g), u 0 = rh ⇤ e (rh e (x 0 ) 1 L g 0 ).
Then D D he (x, u), (x 0 , u 0 ) = D h ⇤ e rh e (x) 1 L (g g 0 ), rh e (x)

Proof. For the entropy, we have

u = rh ⇤ e (rh e (x) 1 L g) = x exp( g/L).
We use the identity for the double Bregman divergence of the entropy [Bauschke, 2006, Example 2.9] which yields D D he (x, u), (x 0 , u 0 ) = D he (x, u x 0 /u 0 ), and thus

D D he (x, u), (x 0 , u 0 ) = D he (x, x exp( g 0 g L )) = D he ✓ x, rh ⇤ e (rh e (x) + 1 L (g 0 g)) ◆ = D h ⇤ e ✓ rh e (x) + 1 L (g 0 g), rh e (x)
◆ .

where we used the duality property of the Bregman divergence (Lemma 5.3).

Using this identity, we can show that the necessary and sufficient conditions are equivalent and formulate interpolability conditions for entropic-smooth functions.

Theorem 5.11 (Entropic-smooth interpolability). Let L > 0, d 1 a dimension, I be a finite index set and

S = {x i , f i , g i } i2I 2 (R d ++ ⇥ R ⇥ R d ) |I| .
Consider the entropic kernel h e . The set S is F he L -interpolable if and only if for every i, j 2 I,

f i f j hg j , x i x j i LD h ⇤ e  rh e (x i )
1 L (g i g j ), rh e (x i ) .

(5.13)

Proof. First, note that the entropy satisfies Assumption 5.1. The implication is a consequence from the necessary condition in Proposition 5.6. Conversely, assume that (5.13) holds. By the identity stated in Lemma 5.10, we have that

D h ⇤ e  rh e (x i ) 1 L (g i g j ), rh e (x i ) = D D he (x i , u i ), (x j , u j ) 8i, j 2 I with u i = x i exp( g i /L
) for i 2 I and therefore, the set S satisfies Condition (i) in Proposition 5.9, which implies that S is F he L -interpolable.

Remark 4. In fact, the "favorable algebra" from Lemma 5.10 also holds for the squared Euclidean norm; therefore, our analysis can be seen as a strict extension of the Euclidean interpolability conditions from [START_REF] Taylor | Smooth Strongly Convex Interpolation and Exact Worst-Case Performance of First-Order Methods[END_REF] to the kernels that satisfy this property as well as Assumption 5.1. However, it seems that the quadratic kernels and the entropy are the only kernels with such favorable properties. As stated by Bauschke and [START_REF] Heinz | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF], these two geometries seem to be "limiting cases in a profound sense".

Writing the PEP as a convex program on the Kullback-Leibler cone

We now use the interpolation conditions to formulate the problem of finding the worst-case scenario of Bregman methods in entropic-smooth settings as a finite-dimensional problem. Then, we show that by relaxing over the dimension d, this problem becomes a convex program on a new cone that we call the Kullback-Leibler cone.

Formulating the PEP on the set of pairwise Kullback-Leibler matrices

Recall that the performance estimation problem (PEP) for entropic-smooth convex minimization writes

maximize f N f ⇤ subject to D he (x ⇤ , x 0 ) = R, f i = f (x i ), g i = rf (x i ) for all i 2 I and some f 2 F he L , x N = A ⇣ x 0 , {(f i , g i )} N 1 i=0 ⌘ , in the variables d, {(x i , f i , g i )} i2I .
Using the interpolations conditions from Theorem 5.11, we can now write explicitly the constraints on the second line, and thus the PEP is equivalent to

maximize f N f ⇤ subject to D he (x ⇤ , x 0 ) = R, f i f j hg j , x i x j i LD h ⇤ e ⇥ rh e (x i ) 1 L (g i g j ), rh e (x i ) ⇤ , i,j 2 I, x N = A ⇣ x 0 , {(f i , g i )} N 1 i=0 ⌘ ,
in the variables d, {(x i , f i , g i )} i2I . We now need to express the algorithmic constraints. For simplicity, we first focus on the case when A is the vanilla Bregman gradient descent (BGD) algorithm with fixed step size > 0, which writes

rh e (x i+1 ) = rh e (x i ) rf (x i ), i = 0 . . . N 1.
Other Bregman-type algorithms which perform an arbitrary linear combination of past gradients can be treated with the same approach.

It is now convenient to introduce the dual variables {s i } i2I such that

s i = rh e (x i ) = log(x i )
for i 2 I. We use the three-point identity [Bauschke et al., 2017, Lemma 3] to express the inner product as

hg j , x i x j i = hrh e (x j ) (rh e (x j ) g j ), x i x j i = D he (x i , rh ⇤ e [rh e (x j ) g j ]) D he (x i , x j ) D he (x j , rh ⇤ e [rh e (x j ) g j ]) = D h ⇤ e (s j g j , s i ) D h ⇤ e (s j , s i ) D h ⇤ e (s j g j , s j )
where the last line uses the result on the Bregman divergence of the conjugate (Lemma 5.3). Hence, using this formulation on the dual space, the PEP rewrites (5.14) in the variables d, {(s i , f i , g i )} i2I . We now observe that the constraints depend linearly on the vector of function values F = (f 0 , . . . , f N , f ⇤ ), as well as the values of the Bregman divergence between some dual vectors which are linearly constrained. We make this idea clear in the following proposition, in which we write the PEP in a generic form.

maximize f N f ⇤ subject to D h ⇤ e (s 0 , s ⇤ ) = R, (f i f j ) D h ⇤ e (s j g j , s i ) + D h ⇤ e (s j , s i ) + D h ⇤ e (s j g j , s j ) LD h ⇤ e ⇥ s i 1 L (g i g j ), s i ⇤ , i,j 2 I, s i+1 = s i g i , i = 0 . . . N 1,
Proposition 5.12. There exist integers m, p, a constraint matrix A 2 R q⇥m , vectors c 0 , c 1 , . . . , c p in R N +2 and matrices B 0 , B 1 , . . . , B p 2 R m⇥m such that A1 m = 0 and the PEP (5.14) is equivalent to the following optimization problem

maximize hc 0 , F i subject to D 2 K m (A), Tr(B 0 D) = R, hc i , F i + Tr(B i D) 0, i = 1 . . . q, (PEP-K)
in the variables F 2 R N +2 and D 2 R m⇥m , where K m (A) is defined as the set of constrained pairwise KL-distance matrices as follows (5.16) in the variables d, {(s i , s 0 i , f i , g i )} i2I , {s 00 ij } i,j2I . Eliminating the redundant {g i } variables, the last three lines of the constraints rewrite

K m (A) := n D h ⇤ e [S,
f N f ⇤ subject to D h ⇤ e (s 0 , s ⇤ ) = R, (f i f j ) D h ⇤ e (s 0 j , s i ) + D h ⇤ e (s j , s i ) + D h ⇤ e (s 0 j , s j ) LD h ⇤ e h s 00 ij , s i i , i,j 2 I, s 00 ij = s i 1 L (g i g j ), i,j 2 I s 0 i = s i g i , i 2 I, s i+1 = s i g i , i = 0 . . . N 1,
s 00 ij = s i 1 L (s i s 0 i s j + s 0 j ), i,j 2 I, s i+1 = s 0 i , i = 0 . . . N 1.
(5.17)

Let us stack the variables {s i } i2I , {s 0 i } i2I , {s 00 i,j } i,j2I in a matrix S 2 R d⇥m with m = 3(N + 2)(N + 2). Define index reorderings , 0 : I ! {1 . . . m}, and 00 : I ⇥ I ! {1 . . . m} such that. S[:, (i)] = s i , i 2 I, S[:, 0 (i)] = s 0 i , i 2 I, S[:, 00 (i, j)] = s 00 ij , i,j 2 I,

The constraints (5.17) are linear and thus can be written

AS T = 0
for some matrix A 2 R q⇥m with q = 2N + 2. Notice that A1 = 0 as the constraints (5.17) are invariant by adding a constant term to each variable. Using the index reorderings and inroducing the distance matrix

D = D h ⇤ e [S, S] 2 R m⇥m , Problem (5.16) becomes equivalent to maximize F N F ⇤ subject to D (0), (⇤) = R,
(F i F j ) D 0 (j), (i) + D (j), (i) + D 0 (j), (j) LD 00 (i,j), (i) , i,j 2 I, D = D h ⇤ e [S, S] for some d 2 N and S 2 R d⇥m such that AS T = 0, in the variables F 2 R N +2 , D 2 R m⇥m . Due to the definition (5.15), the last constraint rewrites D 2 K m (A), and the others are linear, hence the problem is of the form (PEP-K).

The Kullback-Leibler cone

Problem (PEP-K) involves linear constraints, and the set K m (A) that we call the Kullback-Leibler cone with log-linear constraints. Indeed, the matrix A encodes linear constraints on the dual variables s i = log x i . We first show that it is indeed a convex cone.

Proposition 5.13. Let m, q 1, and a constraint matrix A 2 R q⇥m such that A1 m = 0. The the set K m (A) defined in (5.15) is a convex cone.

Additivity. Let D, D 0 2 K m (A) and ↵ > 0. Choose S 2 R n 0 ⇥m , S 0 2 R n⇥m such that D = D h ⇤ e [S, S], D 0 = D h ⇤ e [S 0 , S 0 ].
Consider the matrix S 00 2 R (n+n 0 )⇥m obtained by concatenating S and S 0 as S 00 = [S; S 0 ]. Then A(S 00 ) T = AS T + A(S 0 ) T = 0, and we have

D h ⇤ e [S 00 , S 00 ] = D h ⇤ e [S, S] + D h ⇤ e [S 0 , S 0 ] = D + D 0
which proves that D + D 0 2 K m (A) and concludes the proof.

The set K m (A) is the entropic counterpart of the Euclidean Distance Matrix cone, which can be expressed as a linear transformation of the cone of positive semidefinite matrices [START_REF] Krislock | Euclidean Distance Matrices and Applications[END_REF]. Moreover, in the Euclidean setting, linear constraints on the generating points are equivalent to a set of linear constraints on the distance matrix, which allows formulating the PEP as a semidefinite program [START_REF] Drori | Performance of First-Order Methods for Smooth Convex Minimization: A Novel Approach[END_REF].

To the best of our knowledge, the Kullback-Leibler cone is novel and no method is known for solving problems of type (PEP-K). However, as a first result, we are able to describe the cone in the case where there are no linear constraints, i.e. A = 0. In that case, we simply denote

K m := K m (0)
The unconstrained Kullback-Leibler cone is trivial. It is known that the Kullback-Leibler divergence is less regular than the Euclidean distance, for instance, it is not symmetric. Here, we analyze this structure a step further by showing that the closure of the KL cone K m is trivial. That is, any matrix with nonnegative entries and null diagonal can be approached by a sequence of pairwise KL divergence matrices.

Proposition 5.14 (The closure of K m is trivial). Let m 3. The closure of the cone K m is

K m = n D 2 R m⇥m + : diag(D) = (0 . . . 0) o .
We start by proving a lemma. Let us denote {E ij } 1i,jm the canonical basis of the set of matrices of size m ⇥ m.

Lemma 5.15. Let m 2. The matrix E 1,2 belongs to K m , the closure of K m .

Proof. Let us first show the result for m = 3. Consider the points x t , y t , z t 2 R defined as x t = t, y t = te 1/t , z t = te 1/s where t, s > 0 are to be chosen later. Then the pairwise KL divergences of x t , y t , z t are D he (x t , y t ) = 1 + t(e 1/t 1)

D he (y t , x t ) = t e 1/t (1 + t) D he (x t , z t ) = t( 1 s + e 1/
1 A = 0 @ 0 1 0 0 0 0 0 e 1/s 0 1 A which shows, by taking s ! 0, that E 12 belongs to the closure of K 3 (remember that D h [X t , X t ] = D h ⇤ [log X t , log X t ]
). For m = 2, it suffices to take the set of points {x t , y t }. For m 4, using the same reasoning and by considering the set of points {x t , y t , z t . . . z t }, we show that E 12 2 K m .

We are now able to prove the desired result. Proof of Proposition 5.14. By symmetry of the roles of x 1 , . . . , x m , Lemma 5.15 implies that for every i 6 = j 2 {1 . . . m}, E ij 2 K m . Hence, since it is a convex cone, the cone spanned by these elements also belongs to K m , i.e.,

Cone({E

ij } 1i6 =jm ) ⇢ K m , which proves the result since K m ⇢ D 2 R m⇥m + : diag(D) = (0 . . . 0) = Co({E ij } 1i6 =jm ).
This result demonstrates the lack of regularity of the Kullback-Leibler divergence, since any combination of nonnegative distances can be approached by a pairwise KL matrix, unlike the Euclidean setting where distance matrices have a nontrivial structure [START_REF] Krislock | Euclidean Distance Matrices and Applications[END_REF].

However, the general case of K m (A), where additional log-linear constraints are imposed, cannot be described as easily. Numerical exploration seems to suggest that adding such constraints make the cone non-trivial; but we were not able to provide an explicit description in such setting.

Preliminary numerical results and conjectures

Although we formulated the performance estimation problem for entropic-smooth minimization as a convex program, there is no known solver for such problem, as we are not able to easily describe the Kullback-Leibler cone K m (A) (beyond the setting A = 0, where it has a trivial structure).

In this section, we propose a basic heuristic method for solving the PEP, and provide preliminary numerical results. 2. Solve the restricted linear problem on Co(D 1 , . . . , D ns ), writing maximize hc 0 , F i subject to

D = ↵ 1 D 1 + • • • + ↵ n S D n S , ↵ 0, Tr(B 0 D) = R, hc i , F i + Tr(B i D) 0, i = 1 . . . q, in the variables D, F, ↵ 2 R n S .

A simple class of linear-span Bregman methods

We propose to analyze the following class of Bregman first-order methods:

rh(x k+1 ) = rh(x k ) k X i=0 k,i rf (x k ), k = 0, 1, . . . (BFOM)
where = { i,k } 0ik is a sequence of fixed coefficients (i.e., the method is not adaptive). While choosing k,k = and i,k = 0 when k 6 = i yields the vanilla Bregman Gradient descent method, one can wonder if there exist some other choice of with better worst-case guarantees. To this extend, we propose to use the PEP methodology.

For a given number of iterations N , the worst-case behavior of (BFOM) on the class F h L of entropic-smooth convex functions can be computed by solving the performance estimation

problem maximize f (x N ) f (x ⇤ ) subject to D h (x ⇤ , x 0 ) = R, f 2 F he L , rh(x k+1 ) = rh(x k ) P k i=0 k,i rf (x k ), k = 0, . . . , N 1 (PEP-BFOM)
in the variables f, x 0 , . . . , x N , x ⇤ , and we denote by E N ( ) = val(PEP-BFOM) its optimal value. In the Euclidean setting with h = k • k 2 /2, the class of algorithms (BFOM) covers a large number of first-order methods, including gradient descent and inertial algorithms such as Nesterov's accelerated gradient method [START_REF] Nesterov | A Method of Solving A Convex Programming Problem With Convergence rate O(1/kˆ2)[END_REF]. The approach of finding the "best" coefficients through worst-case performance estimation has been explored in [START_REF] Drori | Performance of First-Order Methods for Smooth Convex Minimization: A Novel Approach[END_REF], [START_REF] Kim | Optimized First-Order Methods for Smooth Convex Minimization[END_REF], [START_REF] Drori | Efficient First-Order Methods for Convex Minimization: a Constructive Approach[END_REF], [START_REF] Kim | Accelerated Proximal Point Method and Forward Method for Monotone Inclusions[END_REF].

However, for Bregman gradient methods on relatively-smooth functions, there is currently no known method with better guarantee than BGD. While, in Chapter 3, we proved that BGD is optimal for generic kernel h, we propose to study if results can be improved with the specific entropic geometry by a method of the class (BFOM).

Resolution for N = 2. As a simple problem, we propose to solve (PEP-BFOM) for N = 2 iterations and several values of the three parameters 0,0 , 1,0 , 1,1 . Hence we study the scheme

rh e (x 1 ) = rh e (x 0 ) 0,0 rf (x 0 ) rh e (x 2 ) = rh e (x 1 ) 1,1 rf (x 1 ) 1,0 rf (x 0 ) (BFOM-2-it)
For homogeneity reasons, we can assume w.l.o.g that L = 1 and R = 1. Following the approach detailed in Section 5.4, we write (PEP-BFOM) as a convex program on the Kullback-Leibler cone, and solve it with the heuristic described above.

In Figure 5.1, we show the numerical results for a fixed value of 0,0 = 1/L, as a function of 1,1 , 1,0 . We compare with the worst-case performance of the same algorithm in the Euclidean setting with h = 1 2 k • k 2 , obtained with the standard PEP methodology [START_REF] Drori | The Exact Information-Based Complexity of Smooth Convex Minimization[END_REF]Teboulle, 2014, Taylor et al., 2017]. Numerical evidence shows the following:

• in the Euclidean case, there is a non-trivial set of values for 1,0 that yields an improved rate, which include well-known momentum-based algorithms.

• In the entropic setting, however, there is no value of 1,0 6 = 0 that allows the functions value to be bounded, i.e., the value of (PEP-BFOM) is finite only for 1,1 , 0,0 2 (0, 1/L] and 1,0 = 0.

In other words, among the class of methods described by (BFOM-2-it), only Bregman Gradient Descent with step size bounded by (0, 1/L] has guaranteed worst-case convergence. ). Black denotes an infinite value. In the Euclidean case (left), the dependence is continous and a large set of values for the parameters is allowed. For the entropic kernel (right), only the choice 11 2 (0, 1/L] and 10 = 0 (corresponding to vanilla BGD) guarantees a finite worst-case objective gap.

From numerical evidence to formal counter-examples. Numerical results shown above are not formal proofs. However, for negative results, they can lead to rigorous arguments by examining the corresponding worst-case functions, which we deduce from the solutions to the PEP. For instance; let us show that the worst-case performance of (BFOM-2-it) is not bounded for 0,0 = 1,1 = 1 and 1,0 > 0 in the entropic setting. Let us choose the one-dimensional function

f f (x) = D h (x, x ⇤ )
for some x ⇤ > 0 to be fixed later. f is convex and entropic-smooth with constant 1. Applying (BFOM-2-it) to f yields

x 1 = x 0 exp ( 0,0 rf (x 0 )) = x ⇤ , x 2 = x 1 exp ( 1,1 rf (x 1 ) 1,0 rf (x 0 )) = x ⇤ ✓ x ⇤ x 0 ◆ 1,0
Then, the PEP objective writes

f (x 2 ) f ⇤ D h (x ⇤ , x 0 ) = x ⇤ ⇣ x⇤ x 0 ⌘ 1,0 ⇣ 1,0 log x⇤ x 0 1 ⌘ + x ⇤ x ⇤ log x⇤ x 0 x ⇤ + x 0 Denoting u := x 0 /x ⇤ , the quantity rewrites f (x 2 ) f ⇤ D h (x ⇤ , x 0 ) = u 1,0 ( 1,0 log u 1) + 1 log u 1 + u ! which shows, for 1,0 > 0, that lim x 0 x⇤ !0 + f (x 2 ) f ⇤ D h (x ⇤ , x 0 ) = +1
and hence that the worst-case performance is unbounded. Similar arguments can be used to show that the value is unbounded in the other situations described before.

case, stochastic gradient methods have proven to be particularly effective thanks to their low cost per iteration. The simplest one, stochastic gradient descent (SGD), consists in updating x t as

x t+1 = arg min u2C hg t , u xi + 1 2⌘ t ku x t k 2
where g t is a gradient estimate such that E [g t ] = rf (x t ). In our case, a natural choice would be g t = rf ⇠t (x t ) for some ⇠ t . The choice of the step size ⌘ t is crucial for obtaining good performances and is typically related to the smoothness of f with respect to the Euclidean norm.

Beyond simply adapting the step size, a powerful generalization of SGD consists in refining the geometry and performing instead Bregman gradient (a.k.a mirror) steps as

x t+1 = arg min u2C hg t , u xi + 1 2⌘ t D h (u, x t ), (BSGD)
where the Euclidean distance has been replaced by the Bregman divergence with respect to a kernel function h, which writes:

D h (x, y) = h(x) h(y) hrh(y), x yi, (6.2)
for all x, y 2 int C.The standard SGD algorithm corresponds to the case where h = 1 2 k • k 2 . However, a different choice of h might better fit the geometry of the set C and the curvature of the function, allowing the algorithm to take larger steps in directions where the objective gradient changes slowly. This choice is guided by the notion of relative smoothness and strong convexity, introduced in Bauschke et al. [2017], [START_REF] Lu | Relatively-Smooth Convex Optimization by First-Order Methods, and Applications[END_REF]. Instead of the squared Euclidean norm for standard smoothness, relative regularity is measured with respect to the reference function h. In most situations, the update of BSGD can be computed in closed form. In some other cases, it might require running an optimization subroutine. However, it can be worth doing so to reduce the overall complexity, if the kernel h is chosen so that the inner objective (BSGD) is a good approximation of the objective function.

Contributions. In this work, we study Bregman stochastic gradient methods for relativelysmooth problems. First, we analyze the convergence of vanilla BSGD with fixed step size, relying on a variance condition that depends on the stochastic gradients at the optimum. We show that this condition leads to a tighter result than that of the previous work from [START_REF] Hanzely | Fastest Rates for Stochastic Mirror Descent Methods[END_REF], especially for the interpolation setting where the stochastic gradients are null at the optimum.

Then, we show that variance-reduction techniques, which are widely used to accelerate traditional Euclidean stochastic methods when the objective has a finite-sum structure, can be adapted to the Bregman setting. Although this generally requires stronger regularity assumptions (such as global smoothness of h and Lipschitz continuity of r 2 h ⇤ ) and an accurate step size schedule, we show that the asymptotical rate of convergence solely depends on relative regularity constants. The same type of results (asymptotic speedup under additional smoothness assumptions) is observed when applying Nesterov-type acceleration to Bregman gradient methods [Hendrikx et al., 2020[START_REF] Hanzely | Accelerated Bregman Proximal Gradient Methods for Relatively Smooth Convex Optimization[END_REF].

Finally, we illustrate the efficiency of our approach on inverse problems with Poisson noise, which are a typical application of relatively-smooth optimization.

Outline. We start by discussing related work in Section 6.2. Section 6.3 introduces the setup and useful lemmas. We study Bregman stochastic gradient descent in Section 6.4, variance reduction techniques in Section 6.5, and describe numerical experiments in Section 6.6.

Related work

Euclidean stochastic gradient methods and variance reduction. Stochastic optimization methods, and in particular SGD, are very efficient when the number of samples is high [START_REF] Bottou | Stochastic gradient descent tricks[END_REF] and are often referred to as "the workhorse of machine learning". The problem with SGD is that, in general, it only converges to a neighbourhood of the optimum unless a diminishing step-size is used. Variance reduction can be used to counter this problem, and many variance-reduced methods have been developed, such as SAG [START_REF] Schmidt | Minimizing finite sums with the stochastic average gradient[END_REF], SDCA [START_REF] Shalev-Shwartz | SDCA without duality, regularization, and individual convexity[END_REF]Zhang, 2013, Shalev-Shwartz, 2016], SVRG [START_REF] Johnson | Accelerating stochastic gradient descent using predictive variance reduction[END_REF] or SAGA [START_REF] Defazio | SAGA: A Fast Incremental Gradient Method With Support for Non-Strongly Convex Composite Objectives[END_REF].

Mirror descent for nonsmooth functions. In the deterministic setting, Bregman gradient descent was originally proposed as the mirror descent scheme for nonsmooth convex optimization [START_REF] Nemirovski | Efficient methods for large-scale convex optimization problems[END_REF]Yudin, 1983, Beck and[START_REF] Beck | Mirror Descent and Nonlinear Projected Subgradient Methods for Convex Optimization[END_REF] 1 . Mirror descent was also extended to the stochastic setting [START_REF] Nemirovski | Robust stochastic approximation approach to stochastic programming[END_REF][START_REF] Lan | An optimal method for stochastic composite optimization[END_REF][START_REF] Zhou | Stochas-tic mirror descent in variationally coherent optimization problems[END_REF][START_REF] Zhang | On the convergence rate of stochastic mirror descent for nonsmooth nonconvex optimization[END_REF][START_REF] Kimon Antonakopoulos | Online and Stochastic Optimization beyond Lipschitz Continuity: A Riemannian Approach[END_REF], but remained focused on nonsmooth objectives.

To the best of our knowledge, variance reduction for Bregman stochastic methods was only studied in [START_REF] Shi | Bregman divergence for stochastic variance reduction: Saddle-point and adversarial prediction[END_REF] in the context of stochastic saddle-point optimization, but without leveraging relative regularity assumptions like we do in this work.

Stochastic Bregman methods for relatively-smooth problems. The introduction of relative smoothness by [START_REF] Heinz | A Descent Lemma Beyond Lipschitz Gradient Continuity: First-Order Methods Revisited and Applications[END_REF] has opened new perspectives for applying Bregman methods to differentiable objectives. In this more recent context, we are only aware of a few works that study the stochastic variant of Bregman gradient descent. [START_REF] Hanzely | Fastest Rates for Stochastic Mirror Descent Methods[END_REF] consider the same setting as us and obtain comparable convergence rates for Bregman SGD, but with a much looser notion of variance, which we discuss more in details in the next section. This is problematic since their bound on the variance is thus proportional to the magnitude of the gradients along the trajectory, and may thus be very large when far from the optimum. In contrast, our definition of variance leverages the stochastic gradients at the optimum, which allows us to obtain significant results without bounded gradients and in the interpolation regime (zero gradients at the optimum). [START_REF] Davis | Stochastic model-based minimization under high-order growth[END_REF] also analyze a similar setting for Bregman SGD, but again with more restrictive assumptions on the noise and boundedness of the gradients.

When preparing the final version of this work, we became aware of the recent paper by [START_REF] Latafat | Bregman finito/miso for nonconvex regularized finite sum minimization without lipschitz gradient continuity[END_REF], who also study variance-reduced Bregman stochastic algorithms for finitesum minimization. Unlike us, they focus more on nonconvex objective functions. They also prove a convergence rate on strongly convex objectives under additional regularity assumptions on f and h, but with a weaker dependence on relatively smoothness constants than ours.

Problem setup and preliminary lemmas

We begin by recalling the base ingredients for relatively-smooth minimization, and prove preliminary lemmas which will be useful for our analysis. For a more general introduction to Bregman methods, we refer the reader to Chapter 1. Let us begin with the blanket assumptions on the kernel function h. Moreover, we additionally assume that h is twice continuously differentiable on int C.

We also recall the notions of relative smoothness and strong convexity [START_REF] Heinz | A Descent Lemma Beyond Lipschitz Gradient Continuity: First-Order Methods Revisited and Applications[END_REF][START_REF] Lu | Relatively-Smooth Convex Optimization by First-Order Methods, and Applications[END_REF]. Definition 6.1. The function f is said to be L-relatively smooth and µ-relatively strongly convex with respect to h if it is differentiable and for all x, y 2 int C,

µD h (x, y)  D f (x, y)  LD h (x, y), (6.3)
where D f is defined similarly to (6.2).

Note that if µ = 0, the left-hand side inequality reduces to assuming convexity of f . Similarly, if h = 1 2 k • k 2 , then D h (x, y) = 1 2 kx yk 2 , and the usual notions of smoothness and strong convexity are recovered. If both functions are two times differentiable, Equation (6.3) can be turned into an equivalent condition on the Hessians:

µr 2 h(x) r 2 f (x) Lr 2 h(x).
Throughout the chapter, we will generally write µ f/h and L f/h to insist on the relative aspect.

Recall that the conjugate h ⇤ of h is defined for y 2 R d as

h ⇤ (y) = sup x2R d hx, yi h(x).
Under Assumption 6.1, h ⇤ is convex and differentiable on R d [Bauschke and Combettes, 2011, Cor. 18.12], and rh ⇤ (rh(y)) = y for y 2 int C. We recall the relation between the Bregman divergence of h and that of its conjugate h ⇤ . Lemma 6.5. Let x 2 int C and ⌘ > 0. If x + satisfies rh(x + ) = rh(x) ⌘g with E [g] = rf (x), rf (x ? ) = 0, then:

E ⇥ D h (x ? , x + ) ⇤ = D h (x ? , x) ⌘D f (x ? , x) ⌘D f (x, x ? ) + E ⇥ D h (x, x + ) ⇤ ,
where the expectation is taken with respect to the choice of g.

Proof. By using the simple algebraic three-point identity which follows form the definition of Bregman divergence [Bauschke et al., 2017, Lemma 3], we get

D h (x ? , x) D h (x ? , x + ) D h (x + , x) = hrh(x + ) rh(x), x ? x + i = ⌘hg, x ? x + i = ⌘hg, x ? xi ⌘hg, x x + i (6.4)
For the first term of the right-hand side, we use the fact that E [g] = rf (x) and rf (x ? ) = 0 to obtain

E [ ⌘hg, x ? xi] = ⌘hrf (x) rf (x ? ), x ? xi = ⌘D f (x ? , x) + ⌘D f (x, x ? ). (6.5)
For the second term, we write

E ⇥ ⌘hg, x x + i ⇤ = E ⇥ hrh(x) rh(x + ), x x + i ⇤ = E ⇥ D h (x, x + ) D h (x + , x) ⇤ . (6.6)
By replacing Equations (6.5) and (6.6) back in the expectation of (6.4), we get

D h (x ? , x) E ⇥ D h (x ? , x + ) + D h (x + , x) ⇤ = ⌘D f (x ? , x)+⌘D f (x, x ? ) E ⇥ D h (x, x + ) + D h (x + , x) ⇤ ,
which yields the desired result after re-arranging terms.

Lemma 6.5 holds for any instance of BSGD that uses an unbiased gradient estimate. The methods that we study in the sequel essentially differ in the way the variance of this estimate is handled. In the case of vanilla Bregman SGD, we now specify a natural condition for bounding this variance. Assumption 6.2. The stochastic gradients g t are such that g t = rf ⇠t (x t ), with E ⇠t [f ⇠t ] = f and f ⇠t is convex and L f/h -relatively smooth with respect to h for all ⇠ t . Besides, there exists a constant 2 0 such that for every t 0,

1 2⌘ 2 E ⇠t h D h ⇤ rh(x t ) 2⌘rf ⇠t (x ? ), rh(x t ) i  2 .
The assumption that the stochastic gradients are actual gradients of stochastic functions which are themselves smooth with respect to h is rather natural, as already discussed in the introduction. It is in particular verified when solving (empirical) risk minimization problems. Assumption 6.2 is a Bregman adaptation of the usual variance at the optimum definition used for instance in [START_REF] Bach | Non-asymptotic analysis of stochastic approximation algorithms for machine learning[END_REF], [START_REF] Mansel Gower | Sgd: General analysis and improved rates[END_REF]. Note that if h ⇤ is µ h -strongly convex with respect to the Euclidean norm, then the assumption is verified for instance when the variance is bounded in `2 norm, that is if

E ⇠t ⇥ krf ⇠t (x ⇤ )k ⇤  µ h 2 ,
since in that case we have

1 2⌘ 2 D h ⇤ rh(x t ) 2⌘rf ⇠t (x ? ), rh(x t )  1 µ h krf ⇠t (x ⇤ )k 2 .
We used the fact that if h is µ h -strongly convex, then h ⇤ is 1/µ h -smooth and its Bregman divergence is upper bounded by the squared Euclidean distance, see, e.g., Bauschke and Combettes [2010].

Remark 5. [START_REF] Hanzely | Fastest Rates for Stochastic Mirror Descent Methods[END_REF] also consider Bregman SGD for relatively-smooth problems, but with a different assumption on the variance.

1 ⌘ t E ⇠t [hrf (x t ) rf ⇠t (x t ), x t+1 xt+1 i]  2 , (6.7)
for t 0, where g t is the stochastic gradient estimate and xt+1 is the output of the (theoretical) Bregman gradient step taken with the true gradient, that is, rh(x t+1 ) = rh(x t ) ⌘ t rf (x t ). Thus, their condition can be written:

1 ⌘ 2 t E ⇠t [D h (x t+1 , xt+1 ) + D h (x t+1 , x t+1 )]  2 ,
so that 2 bounds at each step the distance (in the Bregman sense) between x t+1 and xt+1 . To illustrate why our assumption is weaker, let us consider the case where h is µ h -strongly convex.

In this setting, a sufficient condition for (6.7) to hold is that

1 µ h E ⇠t ⇥ krf (x t ) rf ⇠t (x t )k 2 ⇤  2 ,
while a sufficient condition for our variance definition to hold is (using that rf (x ? ) = 0):

1 µ h E ⇠t ⇥ krf (x ? ) rf ⇠t (x ? )k 2 ⇤  2 ,
which only depends on the magnitude of the gradients at the optimum instead of the variance along the full trajectory since x t is replaced by x ? . In particular, in the interpolation setting where rf ⇠ (x ? ) = 0 for every ⇠, 2 = 0 with our condition.

We now state our convergence result for Bregman SGD. To avoid notation clutter, we generally omit with respect to which variable expectations are taken when clear from the context. Theorem 6.6. If f is L f/h -smooth and µ f/h -strongly convex relative to h with µ f/h > 0, and Assumptions 6.1 and 6.2 hold, then for a step size ⌘  1/(2L f/h ), the iterates produced by Bregman stochastic gradient (BSGD) satisfy

E [D h (x ? , x t )]  (1 ⌘µ f/h ) t D h (x ? , x 0 ) + ⌘ 2 µ f/h .
Proof. By using Lemma 6.5 with x = x t , we obtain:

E ⇠t [D h (x ? , x t+1 )] = D h (x ? , x t ) ⌘D f (x ? , x t ) ⌘D f (x t , x ? ) + E ⇠t [D h (x t , x t+1 )]  D h (x ? , x t ) ⌘D f (x ? , x t ) + E ⇠t [D h (x t , x t+1 )] ,
(6.8) since D f is nonnegative by convexity of f . Using Lemma 6.3 with

x = x t , g 1 = 2⌘ [rf ⇠t (x t ) rf ⇠t (x ? )] , g 2 = 2⌘rf ⇠t (x ? ),
we have that rh(x t+1 ) = rh(x t ) (g 1 + g 2 )/2 and hence the last term can be bounded as

D h (x t , x t+1 )  1 2 D h ⇤ (rh(x t ) g 1 , rh(x t )) + 1 2 D h ⇤ (rh(x t ) g 2 , rh(x t ))
We use Lemma 6.4 (Bregman co-coercivity) to bound the first term as :

D h ⇤ (rh(x t ) g 1 , rh(x t )) = D h ⇤ (rh(x t ) 2⌘ [rf ⇠t (x t ) rf ⇠t (x ? )] , rh(x t ))  2⌘D f ⇠ t (x t , x ? ), so that 1 2 E [D h ⇤ (rh(x t ) g 1 , rh(x t ))]  ⌘D f (x t , x ? ),
by linearity of D f in f and the fact that E ⇠t [f ⇠t ] = f . To bound the second term, we use the variance condition from Assumption 6.2:

1 2 E [D h ⇤ (rh(x t ) g 2 , rh(x t ))] = 1 2 E [D h ⇤ (rh(x t ) 2⌘rf ⇠t (x ? ), rh(x t ))]  ⌘ 2 2 .
By relative strong convexity of f we have D f (x ? , x t )  µ f/h D h (x ? , x t ) and thus inserting these bounds in (6.8) yields

E ⇠t ⇥ D h (x ? , x t+1 ) ⇤  (1 ⌘µ f/h )D h (x ? , x t ) + ⌘ 2 2 .
Remark 6 (Interpolation). In the interpolation setting (when rf ⇠t (x ? ) = 0 for all ⇠ t ), we have that 2 = 0. Theorem 6.6 thus proves linear convergence in this case. For instance, when solving objectives of the form D KL (Ax, b) (which has applications in optimal transport [START_REF] Mishchenko | Sinkhorn Algorithm as a Special Case of Stochastic Mirror Descent[END_REF]) or D KL (b, Ax) (which has application in deblurring or tomographic reconstruction), then the variance as defined in [START_REF] Hanzely | Fastest Rates for Stochastic Mirror Descent Methods[END_REF] may be unbounded, whereas the variance as we define it is equal to 0 if there exists z such that Az = b.

When f is convex but not relatively-strongly convex (µ f/h = 0), the analysis of Theorem 6.6 can be adapted to obtain a 1/T decrease of the error up to a noise region. Theorem 6.7 (Convex case). Under the same assumptions as Theorem 6.6, if

µ f/h = 0, then E " 1 T T X t=0 D f (x ? , x t ) #  D h (x ? , x 0 ) ⌘T + ⌘ 2 (6.9)
Proof. We follow the same analysis as the beginning of the proof of Theorem 6.6. Bounding the terms of (6.8) in same way but by keeping D f (x ? , x t ), we get

⌘D f (x ? , x t ) = D h (x ? , x t ) E ⇠t [D h (x ? , x t+1 )] + ⌘ 2 2 .
Averaging over t and dividing by ⌘ leads to (6.9).

Contrary to the Euclidean case, we do not obtain a guarantee on the average iterate in general. This is because the bound is on the average of D f (x ? , x t ) instead of D f (x t , x ? ), and Bregman divergences are not necessarily convex in their second argument (except for the particular cases of the Euclidean distance and the Kullback-Leibler divergence).

Variance reduction

We have shown in the previous section that BSGD enjoys guarantees that are similar to that of its Euclidean counterpart, although the notion of variance needs to be adapted. We show in this section that it is also possible to apply variance reduction techniques to accelerate convergence, in the case where the objective is a finite sum of the form

min x2C f (x) := 1 n n X i=1 f i (x), (6.10)
As previously, we also assume that the minimizer x ? belongs to int C, so that rf (x ? ) = 0. To solve (6.10), we propose in Algorithm 7 a Bregman adaptation of the SAGA algorithm [START_REF] Defazio | SAGA: A Fast Incremental Gradient Method With Support for Non-Strongly Convex Composite Objectives[END_REF], which is one of the most standard variance-reduced stochastic methods. Following its Euclidean counterpart, Bregman-SAGA uses previously computed gradients rf j ( t j ) in order to reduce the variance of the estimate g t at the current point x t . Note that one can check that the gradient estimate is ubiased as E it [g t ] = rf (x t ).

Algorithm 7 Bregman-SAGA((⌘ t ) t 0 , x 0 )

1: i = x 0 for i = 1, ..., n 2: for t = 0, 1, 2, . . . do

3:

Pick i t 2 {1, ..., n} uniformly at random 4:

g t = rf it (x t ) rf it ( t it ) + 1 n P n j=1 rf j ( t j ) 5: x t+1 = arg min x n hg t , x x t i + 1 ⌘t D h (x, x t ) o 6: t+1 it = x t
, and store rf it ( t+1 it ).

7:

t+1 j = t j for j 6 = i t . 8: end for We now introduce our standing assumptions for analyzing Bregman-SAGA. Compared to the previous section, we need an additional crucial regularity condition on D h . Assumption 6.3. For all i 2 {1, • • • , n}, f i is L f/h -smooth relative to h, and f is µ f/h -strongly convex relative to h. Moreover, there exists a gain function G such that for any x, y, v 2 R d and 2 [ 1, 1], D h ⇤ (x + v, x)  G(x, y, v) 2 D h ⇤ (y + v, y) .

Such structural assumptions appear to be essential for analyzing Bregman-type methods that use information provided by gradients of past iterates. The function G models the fact that the function D h ⇤ (x + v, x) is not generally homogeneous in v nor invariant to translation in x (except for the Euclidean case where it is equal to kvk 2 /2). Note that such difficulties are also encountered for obtaining accelerated rates with inertial variants of Bregman descent, where similar assumptions are needed [START_REF] Hanzely | Accelerated Bregman Proximal Gradient Methods for Relatively Smooth Convex Optimization[END_REF]. This seems unavoidable, as suggested by the lower bound from Chapter 3.

Although the gain function G is relatively abstract at this point, it plays a key role in defining the step-size, and convergence guarantees similar those of Euclidean SAGA can be obtained provided G can be chosen small enough. We first state the general Theorem 6.8 (convergence proof for Algorithm 7), and then detail how G can be bounded in several interesting cases.

For t 0 and step-sizes ⌘ t > 0, define

H t = 1 n n X i=1 D f i ( t i , x ? ),
and the potential t as follows:

t = 1 ⌘ t D h (x ? , x t ) + n 2 H t .
First note that by convexity of h and of the f i , t 0 for all t. Our goal in this section is to show that { t } t 0 converges to 0 at a given speed. Indeed, since D h (x ? , x t )  t , this implies (as in Section 6.4) that x t converges to x ? at the same rate. To ease notations, we define ↵t = 1 n n X j=1 rf j ( t j ), and ↵t i = rf i ( t i ) ↵t .

Expectations are taken either on the full past, or on the index i t conditionally on the past, but we generally omit this dependence to avoid notation clutter. Similarly, we write i instead of i t when clear from the context.

Theorem 6.8. Assume that Algorithm 7 is run with a step size sequence {⌘ t } t 0 satisfying ⌘ t = 1/(8L f/h G t ) for every t 0, with G t decreasing in t and such that for all j 2 {1, • Using the gain function with the fact that ⌘ t  1/L f/h , we have

E it [D 1 ] = E i [D h ⇤ (rh(x t ) 2⌘ t (rf i (x t ) rf i (x ? )) , rh(x t ))]  4L 2 f/h ⌘ 2 t E i h G ✓ x t , x t , 1 L f/h (rf i (x t ) rf i (x ? )) ◆ ⇥ D h ⇤ rh(x t ) 1 L f/h (rf i (x t ) rf i (x ? )) , rh(x t ) i  4L f/h ⌘ 2 t E i  G ✓ x t , x t , 1 L f/h (rf i (x t ) rf i (x ? )) ◆ D f (x t , x ? )  4L f/h ⌘ 2 t G t D f (x t , x ? ).
where we used Lemma 6.4, and the last line is implied by the choice of G t in Equation (6.11). Note that we can pull the G t term out of the expectation over the choice of i since G t holds for all i.

For bounding D 2 , we use a Bregman counterpart of the bias-variance inequality [START_REF] Pfau | A generalized bias-variance decomposition for bregman divergences[END_REF], which we prove in Appendix 6.A for completeness (Lemma 6.14). Let us denote with V the (random) vector V = 2⌘ t (rf i (x ? ) rf i ( t i )). Then, notice that, the expectation of V with respect to the choice of i is

E i [V ] = 2⌘ t n n X j=1 rf j ( t t ) = 2⌘ t ↵ t ,
since rf (x ? ) = 0, and that

D 2 = D h ⇤ (rh(x t ) + V E i [V ] , rh(x t )) .
Therefore, Lemma 6.14 leads to

E i [D 2 ] = E i [D h ⇤ (rh(x t ) + V E i [V ] , rh(x t ))]  E i [D h ⇤ (rh(x t ) + V E i [V ] , rh(x t ) E i [V ])]  4⌘ 2 t L 2 f/h E i " G ✓ rh(x t ) 2⌘ t ↵ t , rh( t j ), 1 L f/h rf i ( t j ) rf i (x ? ) ◆ ⇥ D h ⇤  rh( t i ) 1 L f/h rf i ( t i ) rf i (x ? ) , rh( t i ) #  4L f/h ⌘ 2 t G t E i ⇥ D f i ( t i , x ? )
⇤ .

Where we used the gain function for translating D h ⇤ , and the choice of G t in Equation (6.12).

Recall that H t = 1 n P n j=1 D f j ( t j , x ? ). Plugging the upper bounds of D 1 and D 2 into Equation (6.14), we obtain:

E it [D h (x ? , x t+1 )] D h (x ? , x t )  ⌘ t D f (x ? , x t ) ⌘ t D f (x t , x ? )+2L f/h ⌘ 2 t G t [D f (x t , x ? ) + H t ] .

Following [START_REF] Hofmann | Variance reduced stochastic gradient descent with neighbors[END_REF], we write:

E it [H t+1 ] = ✓ 1 1 n ◆ H t + 1 n D f (x t , x ? ),
Indeed, t+1 j = t j with probability 1 1/n, and t+1 i = x t with probability 1/n. Therefore, we can use the H t /n term to control the excess term from bounding D h (x t , x t+1 ). In the end, using that G t is decreasing and so ⌘ t is increasing, we obtain the following recursion: ).

E it [ t+1 ] t = 1 ⌘ t+1 D h (x
If we choose ⌘ t  1/(8L f/h G t ) then the last term is positive and 1 4⌘ t L f/h G t 1/2, so that using the relative strong convexity of f leads to:

E it [ t+1 ]  (⌘ 1 t µ f/h )D h (x ? , x t ) + ✓ 1 1 2n ◆ n 2 H t  ✓ 1 min ✓ ⌘ t µ f/h , 1 2n ◆◆ t .
The result can then be obtained by chaining this inequality. If µ f/h = 0 then we start back from Equation (6.5), use that D f (x ? , x t ) 0 and the same fact that 1 4⌘ t L f/h G t 1/2 to obtain:

1 4 [D f (x t , x ? ) + H t ]  t E it [ t+1 ] .
The result is obtained by averaging over T , since the right hand side yields a telescopic sum, leading to the 1/T rate of Equation (6.13).

After proving the general convergence theorem, we show how the step size ⌘ t can be chosen in several settings. We start with the most favorable situation, the case where h is quadratic. Corollary 6.9. If h is a quadratic function, Assumption 6.3 is satisfied with G = 1, so that

E [ t ]  ✓ 1 min ✓ 1 8 f/h , 1 2n
◆◆ t 0 ,
where  f/h = L f/h /µ f/h is the relative condition number.

Thus, we recover the result of [START_REF] Defazio | SAGA: A Fast Incremental Gradient Method With Support for Non-Strongly Convex Composite Objectives[END_REF] in the Euclidean setting. Another case where a simple gain function G can be found is if h is smooth and strongly convex with respect to some norm. Corollary 6.10. If h is L h -smooth and µ h -strongly convex with respect to the some norm k•k E , then the stepsize can be chosen constant as

⌘ t = µ h 8L h L f/h
, and

E [ t ]  ✓ 1 min ✓ 1 8 h  f/h , 1 2n
◆◆ t 0 .
Proof. In this setting, the gain function can be chosen constant equal to G = L h /µ h , since we have for x, y, v 2 R d and 2 [ 1, 1]

D h ⇤ (x + v, x)  1 2µ h k vk 2 E  2 L h µ h D h ⇤ (y + v, y)
where we used the fact that h ⇤ is 1/µ h -smooth and 1/L h -strongly convex [START_REF] Heinz | The Baillon-Haddad Theorem Revisited[END_REF]. Then, Theorem 6.8 yields the result.

This convergence result comes with an important drawback: the rate depends on the condition number  h = L h /µ h of h, which can be excessively large or even infinite in the context of relatively-smooth problems.

Asymptotical rates. In the general case, bounding the gain function is a more tedious task. As stated at the beginning of this section, one of the problems is that Bregman divergences lack translation invariance and homogeneity. However, as the algorithm converges, one can expect these conditions to hold locally, as D h ⇤ (x + v, x) is approximated by 1 2 kvk 2 r 2 h ⇤ (x ⇤ ) for small enough v, and x close enough to x ⇤ . This is indeed what happens under enough regularity assumptions on h. Here, we consider the restricted setting where h is globally smooth and its conjugate h ⇤ has Lipschitz Hessian. Proposition 6.11. If h is L h -smooth and the Hessian r 2 h ⇤ is M -smooth, then the gain function can be chosen as:

G(x, y, v) = 1 + 2ML h (ky xk + kvk) .

Proof. Writing the divergence in integral form, we have for x, y, v 2 R d and 2 [ 1, 1]

D h ⇤ (x + v, x) = 2 Z 1 0 Z t 0 v > r 2 h ⇤ (x + s v)v ds dt  2 Z 1 0 Z t 0 ⇣ v > r 2 h ⇤ (y + sv)v + M ky + sv x svkkvk 2 ⌘ ds dt  2 Z 1 0 Z t 0 ⇣ v > r 2 h ⇤ (y + sv)v + M (ky xk + 2skvk) kvk 2 ⌘ ds dt = 2 D h ⇤ (y + v, y) + M (ky xk + kvk)kvk 2 .
Using the fact that is h is L h -smooth, h ⇤ is 1/L h -strongly convex and hence kvk 2  2L h D h ⇤ (y + v, y), leading to

D h ⇤ (x + v, x)  2 [1 + 2ML h (ky xk + kvk)] D h ⇤ (y + v, y).
Note that, even if the regularity conditions of Proposition 6.11 do not hold globally (such as for problems with unbounded curvature), they are at least valid on every bounded subset of int C, as soon as h is C 3 on int C, and we show in the sequel that the dependence on these constants at least disappears asymptotically. We now explicit a possible explicit choice for G t in this setting. Corollary 6.12. Assume that h is L h -smooth, µ h -strongly convex and that the Hessian r 2 h ⇤ is M -smooth. Then, there exists an explicit constant C such that if Algorithm 7 is run with a step size ⌘ t = 1/(8L f/h G t ) with G t decreasing and satisfying

G t min L f/h L h µ h , 1 + C ⇣ n X j=1 kx t t j k + k n X j=1 rf j ( t j )k ⌘ !
, (6.15)

then we have the convergence rate

E [ t+1 ]  ✓ 1 min ✓ 1 8G t  f/h , 1 2n ◆◆ t ,
where lim t!1 G t = 1, or, more precisely,

E [G t ]  1 + O ✓ 1 min ✓ 1 8 h  f/h , 1 2n
◆◆ t .
(6.16)

The explicit expression for the constant C is provided in Appendix 6.A.2 along with the proof. Although the result involves smoothness constants of h which can be large in the relatively-smooth setting, this dependence vanishes asymptotically. Hence, after some time t, which we can roughly estimate using Equation (6.16), we obtain that G t = O(1). Thus, we reach the same kind of convergence rate as in the ideal quadratic case, which depends only on the relative condition number  f/h , but with more general functions h, and thus possibly much better conditioning.

We note however that the choice of G t in (6.15) is rather conservative and has a more theoretical value for now; in our experiments, we will simply choose a constant G t and show that the algorithm behaves well.

Application to Poisson inverse problems

We consider the minimization problem

min x2R d + f (x) = 1 n D KL (b, Ax) (6.17)
where D KL (u, v) = P n i=1 u i log(u i /v i ) u i +v i is the Kullback-Leibler divergence, and A 2 R n⇥d is a typically sparse matrix that models the measurement process. Problem (6.17) models the maximum likelihood estimation problem when assuming the statistical model b ⇠ Poisson(Ax ⇤ ) where x ⇤ is the true unknown signal. Inverse problems with Poisson noise arise in various signal processing applications such as astronomy or computerized tomography, see [START_REF] Bertero | Image Deblurring with Poisson Data: From Cells to Galaxies[END_REF] and references therein. Synthetic problem in the interpolation setting. In Figure 6.1(a), we simulate a synthetic problem the gradients at optimum are zero, by choosing b = Ax ? for some random A 2 R n⇥d and x ⇤ 2 R d (indices sampled uniformly between 0 and 1), with n = 10000 and d = 1000. This corresponds to the interpolation setting, for which Theorem 6.7 predicts fast convergence towards the optimum of Bregman SGD. We observe that BGD is by far the slowest algorithm, but that BSGD is faster than MU thanks to the stochastic speedup. We also observe that BSGD does not plateau in a noise region and converges to the true solution, which is consistent with Theorem 6.7. The step-size for BGD and BSGD is chosen as 1/L f/h , whereas MU is parameter-free.

Tomographic reconstruction problem. Computerized tomography [START_REF] Avinash | Principles of computerized tomographic imaging[END_REF] is the task of reconstructing an object from cross-sectional projections, with fundamental applications to medical imaging. We study a classical synthetic toy problem for this task: the Shepp-Logan phantom (Figure 6.2(a)). In this setting, the observation matrix A corresponds to the discrete Radon transform, which is the cross-sectional projection of the original image x along different projection angles ✓ 1 , . . . , ✓ n (Figure 6.2(b)). That is, the objective writes

f (x) = 1 n D KL (b, Ax) = 1 n n X i=1 D KL (b ✓ i , A ✓ i x) (6.18)
where b ✓ i , A ✓ i correspond to the observation and projection matrix along the angle ✓ i . For stochastic algorithms, the formulation (6.18) naturally yields a finite-sum structure: we thus take f i (x) = D KL (b ✓ i , A ✓ i x) for i = 1 . . . n.

We corrupt the sinogram with Poisson noise, and apply our algorithms. We use n = 360 projection angles, and the image dimension is d = 100 2 . As the matrix A has a sparse structure, we use the relative smoothness constant provided by Proposition 6.13 for a better estimate. The step-size given by theory was rather conservative in this case, so we increased it by a factor of 5 for all Bregman algorithms (and even 10 for BGD).

Figure 6.1(b) shows again that stochastic algorithms drastically outperform BGD. Yet, BSGD quickly reaches a plateau because of the noise. On the other hand, BSAGA enjoys variance reduction and fast convergence to the optimum. In this case, BSAGA is on par with MU, the state-of-the-art algorithm for this problem. One possible explanation is that the log barrier allows relative smoothness to hold, but heavily slows down Bregman algorithms when coordinates are close to 0. Yet, these results are encouraging and one may hope for even faster convergence of BSAGA for tomographic reconstruction with a tighter reference function.

Conclusion

In this work, we have (i) given tight convergence guarantees for Bregman SGD that allow to accurately describe its behaviour in the interpolation setting, and (ii) introduced and analyzed Bregman analogs to the standard variance-reduced algorithm SAGA.

The convergence results on variance reduction require stronger assumptions on the objective than relative smoothness and strong convexity. We show that fast rates can be obtained asymptotically when h is nicely behaved, provided that the step sizes are well chosen. However, in practical experiments, the step size does not need to be taken as conservatively as the theory predicts, and there does not seem to be a slow transient regime. Therefore, a promising extension of our work is to bridge the gap between theory and practice by analyzing this regime and providing more accurate worst-case guarantees. Appendices 6.A Missing proofs for variance reduction

6.A.1 Bregman variance decomposition

First, we use the following Bregman counterpart of a standard variance identity [START_REF] Pfau | A generalized bias-variance decomposition for bregman divergences[END_REF], which we prove for completeness. Lemma 6.14 (Bregman variance decomposition). Let X be a random variable on R d . Then for any u 2 R d ,

E [D h ⇤ (X, u)] = D h ⇤ (E [X] , u) + E [D h ⇤ (X, E [X])] (6.19)
As a consequence, for any random variable V on R d and point y 2 R d we have

E [D h ⇤ (y + V E [V ] , y E [V ])] E [D h ⇤ (y + V E [V ] , y)] .
(6.20)

Proof. Denoting x := E [X], We have for u 2 R d D h ⇤ (x, u) + E [D h ⇤ (X, x)] = h ⇤ (x) h ⇤ (u) hrh ⇤ (u), (x u)i + E h h ⇤ (X) h ⇤ (x) rh ⇤ (x) > (X x) i = h ⇤ (u) rh ⇤ (u) > (x u) + E [h ⇤ (X)] = E h h ⇤ (X) h ⇤ (u) rh ⇤ (u) > (X u) i = E [D h ⇤ (X, u)]
which proves (6.19). Then, (6.20) follows from applying it to the point u = y E [V ] and the random variable X = y + V E [V ], along with using the nonnegativity of the Bregman divergence D h ⇤ (E [X] , u). 6.A.2 Lipschitz-Hessian setting: proof of Corollary 6.12 Corollary 6.15 (Recall of Corollary 6.12). Assume that h is L h -smooth, µ h -strongly convex and the Hessian r 2 h ⇤ is M -smooth. Then, there exists an explicit constant C such that if Algorithm 7 is run with a step size ⌘ t = 1/(8L f/h G t ) with G t decreasing in t and satisfying

G t min L f/h L h µ h , 1 + C ⇣ n X j=1 kx t t j k + k n X j=1 rf j ( t j )k ⌘ !
, then we have the convergence rate (6.21) where lim t!1 G t = 1, or, more precisely, As the quantities involving rf it (x ? ) are unknown, we provide an upper estimate. We can proceed in the following way, using the fact that, due to relative regularity, f i is also smooth with constant L h L f/h , and f is strongly convex with constant µ h µ f/h :

E it [ t+1 ]  ✓ 1 min ✓ 1 8G t  f/h , 1 2n ◆◆ t ,
E [G t ]  1 + O ✓ 1 min ✓ 1 8 h  f/h , 1
krf it (x t ) rf it (x ? )k 2  2L h L f/h D f i t (x t , x ? )  2L h L f/h n D f (x t , x ? )  L h L f/h µ h µ f/h n k 1 n n X j=1 rf j (x t )k 2   f  f/h n 0 @ n X j=1 krf j (x t ) rf j ( t j )k + k n X j=1 rf j ( t j )k 1 A 2   f  f/h n 0 @ n X j=1 L h L f/h kx t t j k + k n X j=1 rf j ( t j )k 1 A 2 .
And similarly, we can estimate the second term from krf it ( t it ) rf it (x ? )k  krf it (x t ) rf it (x ? )k + L h L f/h k t it x t k, which leads to the following upper estimate of the RHS of Condition (6.22):

1+2ML h ⇣ 1 L f/h krf it (x t ) rf it (x ? )k + 1 L f/h krf i ( t it ) rf i (x ? )k + krh(x t ) rh( t it ) 1 4nL f/h n X j=1 rf j ( t j )k ⌘  1 + 2ML h 2 L f/h krf it (x t ) rf it (x ? )k + L h k t it x t k + krh(x t ) rh( t it )k + 1 4nL f/h k n X j=1 rf j ( t j )k !  1 + 2ML h 2 r  h  f/h n L h n X j=1 kx t t j k + 2L h k t it x t k + 1 4nL f/h + 2 L f/h r  h  f/h n ! k n X j=1 rf j ( t j )k !  1 + C 0 @ n X j=1 kx t t j k + k n X j=1 rf j ( t j )k 1 A
where C is defined as

C = 2ML h max 4L h 1 + r  h  f/h n ! , 1 L f/h 1 4n + 2 r  h  f/h n !!
. Now, with such choice of G t , Theorem 6.8 applies and the convergence rate (6.21) holds. It remains to prove the estimate for the convergence rate of G t towards 1. To this end, we show that it is upper bounded by O(1

+ 1/2 t ) since 1 + C n X j=1 kx t t j k + k n X j=1 rf j ( t j )k !  1 + C 0 @ n X j=1 kx t t j k + k n X j=1 rf j (x t )k + n X j=1 krf j ( t j ) rf j (x t )k 1 A  1 + C 0 @ n X j=1 (1 + L h L f/h )kx t t j k + nkrf (x t )k 1 A  1 + C 0 @ n X j=1 (1 + L h L f/h ) kx t x ⇤ k + kx ⇤ t j k + nL h L f/h kx t x ⇤ k 1 A  1 + C 0 @ n(1 + 2L h L f/h kx t x ⇤ k + n X j=1 (1 + L h L f/h )kx ⇤ t j k 1 A  1 + C n(1 + 2L h L f/h ) r 2 µ h D h (x ⇤ , x t ) + n X j=1 (1 + L h L f/h ) s 2 µ h µ f/h D f j ( t j , x ⇤ ) ! = 1 + O 0 @ p D h (x ⇤ , x t ) + n X j=1 q D f j ( t j , x ⇤ ) 1 A = 1 + O ⇣ p t ⌘ (6.23)
Since we imposed a safeguard such that G t L f/h L h µ h , the convergence rate of t is bounded by

E [ t ] = O ✓ 1 min ✓ 1 8 h  f/h , 1 2n
◆◆ t as stated by Corollary 6.10. Indeed, the assumptions are verified as h is L h -smooth and µ hstrongly convex. This worst-case estimate for t , along with the majorization (6.23), gives the resulting rate for G t .

  iii 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Related work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Problem setup and preliminary lemmas . . . . . . . . . . . . . . . . . . . . . . . 6.4 Bregman stochastic gradient descent . . . . . . . . . . . . . . . . . . . . . . . . 6.5 Variance reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.6 Application to Poisson inverse problems . . . . . . . . . . . . . . . . . . . . . . 6.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.A Missing proofs for variance reduction . . . . . . . . . . . . . . . . . . . . . . . .

  Definition 1.1 (Kernel function). A functionh : R d ! R [ {+1} is called a kernel function on C if (i) h is closed convex proper,(ii) h is continuously differentiable and strictly convex on int C, (iii) the Bregman gradient iterates are well-posed, i.e., for every p 2 R d , the problem min u2C hp, ui + h(u) has a unique minimizer, which belongs to int C.

  Figure 1.1: Level lines of the objective function of a linear inverse problem on R 2 + , with the same data A 2 R 2⇥2 , b 2 R 2 and different loss functions. Left: Euclidean loss (Gaussian noise). Right: Kullback-Leibler loss (Poisson noise). Unlike the Euclidean loss, the Poisson loss is not globally smooth as its Hessian diverges near some boundary points of R 2 + .

  X)and ⌧ (c) denotes the unique real solution z to the cubic equation z 2 (z ) = c.Note that ⌧ (c) can be computed in closed form using Cardano's method:

  Figure 2.1: SymNMF normalized objective gap f (X k ) f min /kM k 2 averaged over 10 random initializations, for various sparse similarity matrices M 2 R n⇥n . Hyperparameters for SymHALS, SymANLS were tuned for best performance, while Dyn-NoLips is parameter-free.

  Figure 2.2: Euclidean matrix completion problems on the Helix dataset, with 10% known distances and two different problem sizes. We present the normalized RMSE over the full distance matrix versus CPU time. The results are averaged over 10 random initializations.

Figure 3

 3 Figure 3.1: Level curves of function f (left) and of its smoothed Moreau evelope f µ (right)for n = 2, µ = 0.2 and ⌘ = 1/2. Lemma 3.8 states that if µ is small enough compared to ⌘, the behaviors of f and f µ at x 0 = 0 are the same. Indeed, the size of the smoothed region where the corners are "rounded" decreases when µ goes to 0.

Figure 4 . 1 :

 41 Figure 4.1: Worst-case functions for BGD in dimension 1 with N = 3 iterations. The left figure shows the limiting instance (f, h) 2 B L (R), while the right plot represents the smooth approximation by a valid instance (f µ , h µ ) 2 B L (R), with smoothing parameter µ = 0.1. As µ goes to 0, functions f µ , h µ tend to a pathological behavior where all iterates are equal and for which we have exactly f (x N ) f ⇤ = D h (x ⇤ , x 0 )/N .

  Proof. Scalar multiplication. Let D 2 K m (A) and ↵ > 0. Let us pick S 2 R n⇥m such that D = D h ⇤ [S, S]. Consider the matrix S 0 = S + (log ↵)1 n⇥m . Then, since we assumed A1 m = 0, we have AS 0 = AS. Moreover, since the dual KL divergence D h ⇤ e satisfiesD h ⇤ e (u + log ↵, v + log ↵) = ↵D h ⇤ e (u, v) 8u, v 2 R d ,(see (5.2) for the expression of D h ⇤ e ) we have D h ⇤ e [S 0 , S 0 ] = ↵D h ⇤ e [S, S] = ↵D and hence ↵D 2 K m (A).

  5.5.1 Heuristic solver for Kullback-Leibler cone programsRecall the generic form of a K-L cone programmaximize hc 0 , F i subject to D 2 K m (A), Tr(B 0 D) = R, hc i , F i + Tr(B i D) 0, i = 1 .. . q, in the variables D, F . We propose here the following basic heuristic:1. Randomly sample n S elements D 1 , . . . , D n S of K m (A), and consider the cone spanned by these elements Co(D 1 , . . . , D n S ) ⇢ K m (A).

Figure 5

 5 Figure5.1: Worst-case performance of generic scheme (BFOM-2-it) for N = 2 iterations, with 0,0 = 1/L and for different values of ( 11 , 10 ). Black denotes an infinite value. In the Euclidean case (left), the dependence is continous and a large set of values for the parameters is allowed. For the entropic kernel (right), only the choice 11 2 (0, 1/L] and 10 = 0 (corresponding to vanilla BGD) guarantees a finite worst-case objective gap.

  (a) Toy problem (interpolation setting).

Figure 6 . 1 :

 61 Figure 6.1: Experiments on Poisson inverse problems. SAGA designates Bregman-SAGA (Algorithm 7).

Figure 6 . 2 :

 62 Figure 6.2: Illustration of the Radon transform on the Shepp-Logan phantom. On the sinogram, each column corresponds to the line integral of the image under different projection angles.

  by definition of smoothness and strong convexity with respect to the norm k • k E . However, this is a restrictive setting as many problem of interest are not globally L-smooth. Or, even if they are, tighter relative smoothness constants can be established by eliminating the need for the intermediate bound involving k • k E in the inequality (1.4). We provide such examples in the sequel.

Table 2 .

 2 1: CPU time (in seconds) needed to reach a decrease of ✏ = 10 3 in projected gradient norm (see (2.16) for definition). Results have been averaged over 10 random initializations. Hyperparameters for SymHALS, SymANLS and ADMM have been tuned for best performance. Missing values indicate failure of convergence.

	Dataset r	NoLips	PG	Beta	CD	SymHALS SymANLS ADMM
		10	24.7	51.4	-	26.2	7.0	32.3	-
	Coil-20	20 30	23.7 20.7	36.8 40.8	--	21.3 35.4	4.0 6.5	18.2 20.2	--
		40	21.7	49.5	-	57.6	7.5	28.4	-
		10	38.2	42.7	44.0	35.6	13.6	35.2	42.8
	CBCL	20 30	57.7 60.9	88.4 134.3	--	93.9 135.0	17.8 15.1	47.8 43.4	--
		40	50.8	126.4	-	90.0	23.7	52.5	-
		10	35.2	54.2	-	97.5	11.0	-	-
	TDT2	20 30	52.4 29.4	76.1 45.1	--	109.9 -	20.1 12.1	--	--
		40	28.0	49.8	-	-	17.7	-	-
		10	6.5	10.0	-	33.0	3.0	54.2	-
	Reuters	20 30	28.7 24.3	32.8 45.5	--	71.7 69.4	9.5 6.5	74.7 91.0	--
		40	40.2	68.5	-	83.2	10.6	108.3	-

Table 4 .

 4 1: Numerical value of the performance estimation problem (PEP) with = 1, L = 1. Rel. error denotes the relative error between val(PEP) and the theoretical bound of 1/N given by Theorem 3.5. Primal feasibility corresponds to the maximal absolute value of constraint violation returned by the MOSEK solver.

	N	val(PEP) Rel. error Primal feasibility
	1	1 . 0 0 0	1 . 8 e -1 1	4 . 3 e -1 0
	2	0 . 5 0 0	1 . 8 e -8	2 . 8 e -9
	3	0 . 3 3 3	1 . 8 e -8	2 . 8 e -9
	4	0 . 2 5 0	4 . 9 e -8	2 . 3 e -8
	5	0 . 2 0 0	1 . 8 e -1 0	6 . 4 e -1 1
	10	0.100	6.4e-11	1.3e-11
	20	0.050	1.1e-8	1.9e-10
	50	0.020	6.5e-6	5.0e-7
	100	0.01	7.2e-5	1.6e-6

  S] : there exists d 2 N and S 2 R d⇥m such that AS T = 0

	o
	.	(5.15)
	Proof. Introducing some auxiliary variables {s 0 i } i2I 2 (R d ) |I| and {s 00 ij } i,j2I 2 (R d ) |I| 2 , we reformulate (5.14) as
	maximize

  When t ! 0, the pairwise KL matrix D he [X t , X t ] with X = [x t ; y t ; z t ] satisfies D he [x t , y t , z t ] = lim D he (x t , x t ) D he (x t , y t ) D he (x t , z t ) D he (y t , x t ) D he (y t , y t ) D he (y t , z t ) D he (z t , x t ) D he (z t , y t ) D he (z t , z t )

				s 1)
		D he (z t , x t ) = t(1 e 1/s 1 s	e 1/s )
		D he (y t , z t ) = e 1/t (	t s	1) + t(e 1/s e 1/t )
		D he (z t , y t ) = e 1/s (1	t s	) + t(e 1/t e 1/s )
			0
	t!0 lim	t!0	@

  Assumption 6.1. The function h is a kernel function on C, that is,(i) h : R d ! R [ {+1} is a convex function,(ii) h is finite-valued, strictly convex and differentiable on int C, (iii) the Bregman gradient iterates are well-posed, i.e., for every p 2 R d , the problem min

u2C

hp, ui + h(u) has a unique minimizer, which belongs to int C.

  Proof. Since Bregman-SAGA also uses an update with un unbiased gradient estimate, we apply Lemma 6.5 which yieldsE it [D h (x ? , x t+1 )] = D h (x ? , x t ) ⌘ t D f (x ? , x t ) ⌘ t D f (x t , x ? ) + E it [D h (x t , x t+1 )] .(6.14)In the same way as before, 6.3 impliesD h (x t , x t+1 )  (D 1 + D 2 )/2, with D 1 = D h ⇤ (rh(x t ) 2⌘ t [rf i (x t ) rf i (x ? )] , rh(x t )), D 2 = D h ⇤ (rh(x t ) 2⌘ t (rf i (x ?) ↵t i ), rh(x t )).

											◆
											,	(6.11)
	G t G	⇣	rh(x t ) 2⌘ t	↵t , rh( t j ),	1 L f/h	(rf j ( t j ) rf j (x ? ))	⌘	.	(6.12)
	Then, under Assumptions 6.1 and 6.3, the potential t satisfies
			E it [ t+1 ] 	✓	1 min	✓	⌘ t µ f/h ,	1 2n	◆◆	t ,
	In the convex case (µ f/h = 0), we obtain that
			E	"	1 4T	T X t=1	D f (x t , x ? ) + H t	#		0 T	.	(6.13)

• • , n}: G t G ✓ rh(x t ), rh(x t ), 1 L f/h (rf j (x t ) rf j (x ? ))

  (x ? , x t+1 ) D h (x ? , x t )) +

	? , x t+1 ) +	n 2	H t+1	1 ⌘ t	D h (x ? , x t )	n 2	H t
	  D f (x ? , x t ) 1 ⌘ t (D h n 2 1 2 1 4⌘ t L f/h G t H t (H t+1 H t ) ✓ 1 2⌘ t L f/h G t	1 2	◆	D

f (x t , x ?

  2n◆◆ t .Proof. Using the gain function from Proposition 6.11, to satisfy the assumptions of Theorem 6.8 it is sufficient to choose G t such thatG t 1 + 2ML h ⇣ 1 L f/h krf it (x t ) rf it (x ? )k + 1 L f/hkrf it ( t it ) rf it (x ? )k

	+ krh(x t ) rh( t it )	1 4nL f/h	n X j=1	j )k rf j ( t	⌘	.	(6.22)

Interestingly, the author does not make the connection with other work on relative smoothness and seems to have discovered the same notion independently.

There were actually earlier versions of accelerated algorithms proposed by A. Nemirovski which were slightly more complex as they required to perform a two-dimensional minimization subroutine, see e.g.,[START_REF] Bubeck | Nemirovski's acceleration[END_REF].

This idea of implicit bias of gradient methods has been re-discovered recently in the context of mirror descent[START_REF] Gunasekar | Characterizing implicit bias in terms of optimization geometry[END_REF].

Note that mirror descent and Bregman gradient refer to the same algorithm, but that mirror descent is typically used when f is non-smooth, or in the online optimization community, whereas Bregman gradient is generally preferred when using the relative smoothness assumption. However, there is no consensus in the litterature; for instance,[START_REF] Hanzely | Fastest Rates for Stochastic Mirror Descent Methods[END_REF] use the mirror descent terminology although they consider the relatively-smooth setting.

Remerciements

Code

The code for reproducing experiments for SymNMF and Euclidean Distance Matrix Completion can be downloaded from the public repository https://github.com/RaduAlexandruDragomir/QuarticLowRankOptimization

Chapter 4

Computer-Aided Analyses of Bregman Gradient Methods with Generic Kernels

Chapter Abstract

We show how worst-case scenarios of Bregman gradient methods with general kernel functions can be computed by solving appropriate semidefinite programs. Additionally, the corresponding convergence proof can be inferred from the solution to the dual program. This technique, called performance estimation, has been pioneered by [START_REF] Drori | Performance of First-Order Methods for Smooth Convex Minimization: A Novel Approach[END_REF] in the context of Euclidean smooth convex minimization. We extend the performance estimation framework to Bregman methods for relatively-smooth problems, and use it to provide several complexity results in this setting. In particular, numerically generated worst-case examples were used as a basis for obtaining the general lower bound presented in Chapter 3. where D h (x, y) = h(x) h(y) hrh(y), x yi is the Bregman divergence induced by the kernel function h. If f is L-smooth relative to h, that is, if Lh f is convex, then the iterates of the (BGD) scheme with step size = 1/L satisfy the convergence rate

for every u 2 dom h (see Chapter 1). Despite this result, there are still some open questions regarding the complexity of Bregman methods for relatively-smooth minimization. In particular, is the bound (4.1) optimal? Is there another algorithm that achieves a better rate in the same setting, for generic kernels h?

Performance estimation problems.

In this work, we adopt a computer-aided technique for analyzing the worst-case behavior of first-order methods, using performance estimation problems (PEPs). PEPs were first introduced by Drori and [START_REF] Drori | Performance of First-Order Methods for Smooth Convex Minimization: A Novel Approach[END_REF] for analyzing the exact convergence rate of Euclidean gradient descent on L-smooth functions. By solving appropriate semidefinite programs, PEPs allow to:

1. compute (numerically) the exact worst-case complexity of an algorithm on a given class of problems after a fixed number of iterations, 2. study the corresponding worst-case functions, 3. infer an analytical worst-case guarantee by obtaining a feasible point to the dual PEP. Such dual feasible points correspond to combinations of inequalities that certify the convergence bound.

The PEP technique was used to analyze the worst-case complexity of gradient methods in several settings, such as smooth convex minimization [START_REF] Drori | Performance of First-Order Methods for Smooth Convex Minimization: A Novel Approach[END_REF][START_REF] Taylor | Exact Worst-Case Performance of First-Order Methods for Composite Convex Optimization[END_REF][START_REF] Kim | Optimized First-Order Methods for Smooth Convex Minimization[END_REF][START_REF] Drori | The Exact Information-Based Complexity of Smooth Convex Minimization[END_REF][START_REF] Drori | Efficient First-Order Methods for Convex Minimization: a Constructive Approach[END_REF][START_REF] Barré | Complexity Guarantees for Polyak Steps with Momentum[END_REF], nonsmooth convex minimization [START_REF] Drori | The Exact Information-Based Complexity of Smooth Convex Minimization[END_REF][START_REF] Drori | An Optimal Variant of Kelley's Cutting-Plane Method[END_REF]Taylor, 2019], stochastic optimization [Taylor andBach, 2019, Drori and[START_REF] Drori | The Complexity of Finding Stationary Points with Stochastic Gradient Descent[END_REF] and monotone operators [START_REF] Ryu | Operator Splitting Performance Estimation: Tight Contraction Factors and Optimal parameter Selection[END_REF][START_REF] Kim | Accelerated Proximal Point Method and Forward Method for Monotone Inclusions[END_REF].

Contributions. We propose to use the PEP methodology for computing the worst-case behavior of Bregman gradient methods on the set B L (C) of relatively-smooth problem instances with general kernels:

To this end, we adapt the interpolation conditions of [START_REF] Taylor | Smooth Strongly Convex Interpolation and Exact Worst-Case Performance of First-Order Methods[END_REF] to handle the class of differentiable and strictly convex functions, which appears in the relatively-smooth setting. We then use topological arguments to show that the resulting performance estimation problem is equivalent to a simpler limiting problem on a larger class B L (C), which can be seen as the closure of B L (C) and involves possibly nonsmooth functions:

Assumption (iv) the Bregman divergence

•) is strictly convex on C and coercive.

Among these conditions, item (iv) is particularly restrictive; fortunately, besides the squared Euclidean norm, it is also verified for the entropic kernel h e [START_REF] Heinz | Joint and Separate Convexity of the Bregman Distance[END_REF]. All other conditions from Assumption 5.1 also hold for h e .

Under these assumptions, we can now state the properties of the right Bregman-Moreau envelope studied by [START_REF] Heinz | Regularizing with Bregman-Moreau envelopes[END_REF].

Proposition 5.7. (Relative smoothness of the right Bregman-Moreau envelope) Assume that f is a finite-valued convex function on R d which is bounded from below. Let h be a kernel function on C which satisfies Assumption 5.1. Then, for µ > 0, the right Bregman-Moreau envelope env h µ f defined in (5.6) satisfies the following properties:

In other words, we have et al., 2018, Prop. 2.7 (ii)], and differentiability on int C follows from [Bauschke et al., 2018, Prop. 2.19 (ii)]. Finally, by [Bauschke et al., 2018, Prop. 2.4 (ii)] we have Combettes, 2011, Prop. 13.11], and thus that relative smoothness holds.

Let also define the right Bregman-Moreau proximal map for x 2 R d :

which we will often write u(x) for shortness when the context is clear. The next proposition caracterizes the Breman-Moreau proximal map and establishes the link with the gradient of the envelope.

Proposition 5.8. Let x 2 int C. Under the same assumptions as Proposition 5.7, the proximal map u

(5.7)

Proof. The first equation follows from the first-order optimality conditions defining u(x), and the second is a consequence of [Bauschke et al., 2018, Prop. 2.19 (ii)].

Using the definition of env h µ f and the proposition above, we can note that for u(x) = prox h µ f (x) we have

where r 1 and r 2 denote the partial gradients in each variable of D h (•, •).

Double Bregman divergence. As the function D

⇥R d , we can consider its Bregman divergence D D h , which we call the double Bregman divergence. For x, y, u, v 2 int C, it is defined by

(5.9)

The double divergence, which was already studied by [START_REF] Heinz | Joint minimization with alternating Bregman proximity operators[END_REF], will play a central role in the sufficient conditions for F h L -interpolability, which we are now ready to state. These conditions rely on a two-stage procedure: first, we reduce the problem of interpolating a relativelysmooth convex function to the problem of interpolating a convex function, by the means of the Bregman-Moreau envelope. Then, we use the interpolation result for convex functions from [START_REF] Taylor | Smooth Strongly Convex Interpolation and Exact Worst-Case Performance of First-Order Methods[END_REF] (which we recalled in Theorem 4.4).

Proposition 5.9 (Sufficient interpolation conditions). Consider a kernel h satisfying Asssumption 5.1, with zone C ⇢ R d . Let L > 0, I be a finite index set and

Consider the following statements.

(i) For all i, j 2 I we have

with

(5.10)

Conclusion

In this chapter, we proposed to study computer-aided analyses of Bregman methods for entropicsmooth convex functions. First, we established interpolation conditions that allow to write the PEP as a finite-dimensional problem. These conditions rely on the smoothing property of the Bregman-Moreau envelope, as well as favorable algebraic properties of the entropic kernel (such as joint convexity of the Bregman divergence). Then, we formulated the PEP as a convex problem on the cone of pairwise Kullback-Leibler divergences with log-linear constraints. Although we showed that in the setting with no constraints, this cone is trivial, we are not aware of a simple description in the general case. Therefore, we had to use heuristics for solving simple versions of the PEP.

In particular, we demonstrated on a small example of two iterations that among all general Bregman first-order methods with fixed coefficients, only Bregman gradient descent with step size smaller than 1/L is guaranteed to converge in the worst case. This shows that, even if we consider the specific entropic setting instead of generic kernels like in the previous chapters, building accelerated Bregman methods is a tedious task.

There are several important questions that arise for future work:

1. Can we find an algorithm that provably solves convex problems on the Kullback-Leibler cone with log-linear constraints? Note that, since PEPs are usually low-dimensional problems, we do not need efficient solvers, but rather ones with strong guarantees.

2. How to design accelerated Bregman methods in the entropic-smooth setting? An idea would be to consider methods like (BFOM) but with adaptive coefficients, in the spirit of the accelerated Bregman methods proposed in [START_REF] Hanzely | Accelerated Bregman Proximal Gradient Methods for Relatively Smooth Convex Optimization[END_REF], Hendrikx et al. [2020]. An additional difficulty in that case would be that the PEP methodology is not directly applicable to adaptive methods.

Chapter 6

Bregman Stochastic Gradient Descent and Variance Reduction

Chapter Abstract

We study the problem of minimizing an expectation of relatively-smooth functions using Bregman stochastic gradient methods. We first study Bregman stochastic gradient descent, and show convergence towards a region that depends on the magnitude of the gradients at the optimum. Then, we consider variance reduction methods for problems whose objective has a finite-sum structure. In this setting, we show that improved convergence rates can be obtained under additional regularity assumptions on the Bregman kernel. We provide numerical experiments on Poisson inverse problems.

Reference: this chapter is based on a publication in International Conference on Machine Learning [Dragomir et al., 2021b]. This work has been done in collaboration with Hadrien Hendrikx and Mathieu Even, and each student contributed equally.

Introduction

We consider the minimization problem

C is a convex subset of R d and f ⇠ are differentiable convex functions. These problems typically arise in machine learning when performing (empirical) risk minimization, in which case f ⇠ is for instance a loss function for some random sample ⇠. Problem (6.1) is also encountered in signal processing applications such as image deblurring or tomographic reconstruction inverse problems, in which the goal is to recover an unknown signal from a large number of noisy observations. First-order methods are often very efficient for solving such problems, but computing a gradient rf might be very expensive for large-scale problems (large number of components f ⇠ ), and even impossible in the case of true risk minimization (infinite number of f ⇠ ). In this Lemma 6.2 (Duality). For x, y 2 int C, we have D h (x, y) = D h ⇤ (rh(y), rh(x)).

See, e.g., Bauschke and Borwein [1997, Thm 3.7.] for the proof. Using duality, we prove the following key lemma. Lemma 6.3. Let x 2 int C, and g 1 , g 2 2 R d . Define the points x + , x + 1 , x + 2 as the unique points satisfying

Then we have

Proof. Using Lemma 6.2 (duality), we have:

where the inequality step is obtained by convexity of the Bregman divergence in its first argument. The final result is obtained by using Lemma 6.2 back.

In the Euclidean case h = k • k 2 , we recover k g 1 +g 2 2 k 2  1 2 kg 1 k 2 + kg 2 k 2 . We also recall the Bregman counterpart of an inequality linked to co-coercivity of the gradients [Nesterov, 2003, Eq. 2.1.7], which we proved in Chapter 5 (Proposition 5.5). Lemma 6.4. If a convex function f is L-smooth relative to the kernel h, then for any ⌘  1 L and x, y 2 int C,

Bregman stochastic gradient descent

In this section, we study the (BSGD) scheme with fixed step size ⌘ > 0. Recall that, by writing the first-order optimality conditions, and by strict convexity of h, the iterations can be equivalently defined as rh(x t+1 ) = rh(x t ) ⌘g t with E [g t ] = rf (x t ). We denote by x ? the minimizer of f , and we assume x ?

2 int C so that rf (x ? ) = 0; the case where the solution lies on the boundary of C poses additional technical difficulties and is left for future work. We first establish a lemma satisfied by the updates of Bregman stochastic gradient methods with unbiased gradient estimate.

As a motivating application of relative smoothness, Bauschke et al. [2017] prove that the Poisson objective f is relatively smooth with respect to the log-barrier reference function

P n j=1 b j /n. This constant can be quite conservative when A is a sparse matrix, and so we prove a better estimate by leveraging this structure. For j 2 {1 . . . n}, we denote S j the support of the j-th column of A, that is

Proposition 6.13. The Poisson objective function defined in (6.17) is L-smooth relative to the log-barrier for

Proof. Let us denote A 1 , . . . A n the row vectors of A. We refine the analysis from Bauschke et al. [2017, Lemma 7] and start by writing for

Applying the Jensen inequality to the function t 7 ! t 2 and weights

where we used the fact that w ij 2 [0, 1] if i 2 S j , and w ij = 0 otherwise.

The relative Lipschitz constant provided by Proposition 6.13 can be considerably smaller than P n j=1 b j /n when A is sparse, which is the case in practical applications. For our numerical experiments, we compare full-batch Bregman gradient descent (BGD), Bregman stochastic gradient descent (BSGD), and the Bregman SAGA scheme described in Algorithm 7. We also implement the Multiplicative Update (MU), also known as Lucy-Richardson or Expectation-Maximization [START_REF] Shepp | Maximum likelihood reconstruction for emission tomography[END_REF], which is the standard baseline for Poisson inverse problems.

Conclusion and Perspectives

In this thesis, we studied several aspects of Bregman methods for relatively-smooth optimization. We focused on practical applications to low-rank problems, extensions to stochastic variants for large-scale problems from machine learning and signal processing, as well as questions of theoretical complexity and acceleration. Yet, there are many open problems and research directions left to explore.

Homogeneity and additional regularity assumptions. As mentioned throughout this thesis, a central issue in the analyses and extensions of Bregman methods is the lack of homogeneity and translational invariance of the Bregman divergence D h . This problem appears for acceleration [START_REF] Hanzely | Accelerated Bregman Proximal Gradient Methods for Relatively Smooth Convex Optimization[END_REF], as well as for variance reduced stochastic methods (Chapter 6). In fact, it seems to be an obstacle for every algorithm that has memory and combines gradients taken at different points.

The lower bound from Chapter 3 demonstrates that that, in the worst case, it is indeed hopeless to prove accelerated rates for Bregman methods. The corresponding worst-case instance involves pathological functions that a nearly nonsmooth (for which the Bregman divergence is highly non-homogeneous) and that are not quite representative of usual applications. This shows that, in order to study methods with memory, additional regularity assumptions are needed. Several possibilities are available, such as Lipschitz regularity of r 2 h ⇤ , or selfconcordance ideas [START_REF] Sun | Generalized self-concordant functions: a recipe for Newton-type methods[END_REF]. The study on the specific setting of entropy (Chapter 5) can also be pursued.

Despite these theoretical difficulties, numerical experiments show good performance of both accelerated variants [Hendrikx et al., 2020] and variance reduction (Chapter 6). While the current analyses predict that these methods only work in a neighborhood of the optimum, or with asymptotical learning rates, a slow transient regime is not observed in experiments. These results should encourage future work towards bridging the gap between theory and practice.

Convergence on the boundary. A less known problem is that the behavior of Bregman methods is not well understood when the optimum x ⇤ lies on the boundary of the set C. In this setting, the O(1/k) complexity bound does not hold if h is infinite on the boundary (such as for the log kernel). Additionally, the iterate sequence {x k } k 0 is not guaranteed to converge, as demonstrated by the counter-example in [START_REF] Bolte | Curiosities and counterexamples in smooth convex optimization[END_REF]. This is a key theoretical problem and is relevant in many applications, such as Poisson inverse problems on the nonnegative orthant.

Adaptive variants for improving numerical performance. In this thesis, we mostly focused on theoretical analyses of complexity and studied the convergence rate of the algorithms in the worst-case. This approach is fruitful, as it has lead to many important advances in firstorder optimization [START_REF] Nemirovski | Efficient methods for large-scale convex optimization problems[END_REF]Yudin, 1983, Nesterov, 1983].

A complementary approach is to improve practical performance of these methods through adaptiveness. Indeed, in some situations, the relative smoothness inequality is too conservative and Bregman methods exhibit slow convergence speed (see the comparison between BGD and Lucy-Richardson in Chapter 6). In the future, it should be interesting to study methods for automatically adapting the Bregman geometry to the objective, in the spirit of adaptive steps sizes [START_REF] Duchi | Adaptive subgradient methods for online learning and stochastic optimization[END_REF] or Quasi-Newton methods [START_REF] Dennis | Quasi-Newton Methods, Motivation and Theory[END_REF].