Hongzhou Lin
email: hongzhou@mit.edu

Julien Mairal
email: julien.mairal@inria.fr

Zaid Harchaoui

An Inexact Variable Metric Proximal Point Algorithm

Introduction

Convex composite optimization arises in many scientific fields, such as image and signal processing or machine learning. It consists of minimizing a real-valued function composed of two convex terms:

min x∈R d f (x) f 0 (x) + ψ(x) , ( 1
)
where f 0 is smooth with Lipschitz continuous derivatives, and ψ is a regularization function which is not necessarily differentiable. A typical example from the signal and image processing literature is the 1 -norm ψ(x) = x 1 , which encourages sparse solutions [START_REF] Elad | Sparse and Redundant Representations[END_REF][START_REF] Mairal | Sparse modeling for image and vision processing[END_REF]; composite minimization also encompasses constrained minimization when considering extended-valued indicator functions ψ that may take the value +∞ outside of a convex set C and 0 inside (see [START_REF] Hiriart-Urruty | Convex Analysis and Minimization Algorithms I[END_REF]). In general, algorithms that are dedicated to composite optimization only require to be able to compute efficiently the proximal operator of ψ: p ψ (y) arg min

x∈R d ψ(x) + 1 2 x -y 2 ,
where • denotes the Euclidean norm. Note that when ψ is an indicator function, the proximal operator corresponds to the simple Euclidean projection.

To solve (1), significant efforts have been devoted to (i) extending techniques for smooth optimization to deal with composite terms; (ii) exploiting the underlying structure of the problem-is f a finite sum of independent terms? Is ψ separable in different blocks of coordinates? (iii) exploiting the local curvature of the smooth term f to achieve faster convergence than gradient-based approaches when the dimension d is large. Typically, the first point is well understood in the context of optimal first-order methods, see [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF][START_REF] Nesterov | Gradient methods for minimizing composite functions[END_REF], and the third point is tackled with effective heuristics such as L-BFGS when the problem is smooth [START_REF] Liu | On the limited memory BFGS method for large scale optimization[END_REF][START_REF]Updating quasi-Newton matrices with limited storage[END_REF]. Yet, addressing all these challenges at the same time is difficult, which is precisely the focus of this paper.

In particular, a problem of interest that initially motivated our work is that of empirical risk minimization (ERM); the problem arises in machine learning and can be formulated as the minimization of a composite function f : R d → R:

min x∈R d f (x) 1 n n i=1 f i (x) + ψ(x) , (2)
where the functions f i are convex and smooth with Lipschitz continuous derivatives, and ψ is a composite term, possibly non-smooth. The function f i measures the fit of some model parameters x to a specific data point indexed by i, and ψ is a regularization penalty to prevent over-fitting. To exploit the sum structure of f , a large number of randomized incremental gradient-based techniques have been proposed, such as SAG [START_REF] Schmidt | Minimizing finite sums with the stochastic average gradient[END_REF], SAGA [START_REF] Defazio | Saga: A fast incremental gradient method with support for non-strongly convex composite objectives[END_REF], SDCA [START_REF] Shalev-Shwartz | Accelerated proximal stochastic dual coordinate ascent for regularized loss minimization[END_REF], SVRG [START_REF] Xiao | A proximal stochastic gradient method with progressive variance reduction[END_REF], Finito [START_REF] Defazio | Finito: A faster, permutable incremental gradient method for big data problems[END_REF], or MISO [START_REF]Incremental majorization-minimization optimization with application to large-scale machine learning[END_REF]. These approaches access a single gradient ∇f i (x) at every iteration instead of the full gradient (1/n) n i=1 ∇f i (x) and achieve lower computational complexity in expectation than optimal first-order methods [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF][START_REF] Nesterov | Gradient methods for minimizing composite functions[END_REF] under a few assumptions. Yet, these methods are unable to exploit the curvature of the objective function; this is indeed also the case for variants that are accelerated in the sense of Nesterov [START_REF] Frostig | Un-regularizing: approximate proximal point and faster stochastic algorithms for empirical risk minimization[END_REF][START_REF] Lin | A universal catalyst for first-order optimization[END_REF][START_REF] Shalev-Shwartz | Accelerated proximal stochastic dual coordinate ascent for regularized loss minimization[END_REF].

To tackle [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF], dedicated first-order methods are often the default choice in machine learning, but it is also known that standard Quasi-Newton approaches can sometimes be surprisingly effective in the smooth case-that is when ψ = 0, see, e.g., [START_REF] Schmidt | Minimizing finite sums with the stochastic average gradient[END_REF] for extensive benchmarks. Since the dimension of the problem d is typically very large (d ≥ 10 000), "limited memory" variants of these algorithms, such as L-BFGS, are necessary to achieve the desired scalability [START_REF] Liu | On the limited memory BFGS method for large scale optimization[END_REF][START_REF]Updating quasi-Newton matrices with limited storage[END_REF]. The theoretical guarantees offered by L-BFGS are somewhat limited, meaning that it does not outperform accelerated first-order methods in terms of worst-case convergence rate, and also it is not guaranteed to correctly approximate the Hessian of the objective. Yet, L-BFGS remains one of the greatest practical success of smooth optimization. Adapting L-BFGS to composite and structured problems, such as the finite sum of functions [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF], is of utmost importance nowadays.

For instance, there have been several attempts to develop a proximal Quasi-Newton method [START_REF] Byrd | An inexact successive quadratic approximation method for L-1 regularized optimization[END_REF][START_REF] Lee | Proximal Newton-type methods for convex optimization[END_REF][START_REF] Scheinberg | Practical inexact proximal quasi-Newton method with global complexity analysis[END_REF][START_REF] Yu | A quasi-Newton approach to non-smooth convex optimization[END_REF]. These algorithms typically require computing many times the proximal operator of ψ with respect to a variable metric. Quasi-Newton steps were also incorporated as local search steps into accelerated first-order methods to further enhance their numerical performance [START_REF] Ghadimi | Generalized Uniformly Optimal Methods for Nonlinear Programming[END_REF]. More related to our work, L-BFGS is combined with SVRG for minimizing smooth finite sums in [START_REF] Gower | Stochastic block BFGS: Squeezing more curvature out of data[END_REF]. The scope of our approach is broader beyond the case of SVRG. We present a generic Quasi-Newton scheme, applicable to a large-class of first-order methods for composite optimization, including other incremental algorithms [START_REF] Defazio | Saga: A fast incremental gradient method with support for non-strongly convex composite objectives[END_REF][START_REF] Defazio | Finito: A faster, permutable incremental gradient method for big data problems[END_REF][START_REF]Incremental majorization-minimization optimization with application to large-scale machine learning[END_REF][START_REF] Schmidt | Minimizing finite sums with the stochastic average gradient[END_REF][START_REF] Shalev-Shwartz | Accelerated proximal stochastic dual coordinate ascent for regularized loss minimization[END_REF]] and block coordinate descent methods [START_REF] Razaviyayn | A unified convergence analysis of block successive minimization methods for nonsmooth optimization[END_REF][START_REF] Richtárik | Iteration complexity of randomized block-coordinate descent methods for minimizing a composite function[END_REF] More precisely, the main contribution of this paper is a generic meta-algorithm, called QNing (the letters "Q" and "N" stand for Quasi-Newton), which uses a given optimization method to solve a sequence of auxiliary problems up to some appropriate accuracy, resulting in faster global convergence in practice. QNing falls into the class of inexact proximal point algorithms with variable metric and may be seen as applying a Quasi-Newton algorithm with inexact (but accurate enough) gradients to the Moreau envelope of the objective. As a result, our approach is (i) generic, as stated previously; (ii) despite the smoothing of the objective, the sub-problems that we solve are composite ones, which may lead to exactly sparse iterates when a sparsity-inducing regularization is involved, e.g., the 1 -norm; (iii) when used with L-BFGS rules, it admits a worst-case linear convergence rate for strongly convex problems similar to that of gradient descent, which is typically the best guarantees obtained for L-BFGS schemes in the literature.

The idea of combining second-order or Quasi-Newton methods with Moreau envelope is in fact relatively old. It may be traced back to variable metric proximal bundle methods [START_REF] Chen | Proximal quasi-Newton methods for nondifferentiable convex optimization[END_REF][START_REF] Fukushima | A globally and superlinearly convergent algorithm for nonsmooth convex minimization[END_REF][START_REF] Mifflin | A quasi-second-order proximal bundle algorithm[END_REF], which aims to incorporate curvature information into the bundle methods. Our approach revisits this principle with a limited-memory variant (to deal with large dimension d), with a simple line search scheme, with several warm start strategies for the sub-problems and with a global complexity analysis that is more relevant than convergence rates that do not take into account the cost per iteration.

To demonstrate the effectiveness of our scheme in practice, we evaluate QNing on regularized logistic regression and regularized least-squares, with smooth and nonsmooth regularization penalities such as the Elastic-Net [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF]. We use large-scale machine learning datasets and show that QNing performs at least as well as the recently proposed accelerated incremental algorithm Catalyst [START_REF] Lin | A universal catalyst for first-order optimization[END_REF], and other Quasi-Newton baselines such as proximal Quasi-Newton methods [START_REF] Lee | Proximal Newton-type methods for convex optimization[END_REF] and Stochastic L-BFGS [START_REF] Moritz | A linearly-convergent stochastic L-BFGS algorithm[END_REF] in all numerical experiments, and significantly outperforms them in many cases.

The paper is organized as follows: Section 2 presents related work on Quasi-Newton methods such as L-BFGS; we introduce QNing in Section 3 and its convergence analysis in Section 4; Section 5 is devoted to numerical experiments and Section 6 concludes the paper.

Related work and preliminaries

The history of Quasi-Newton methods can be traced back to the 1950's [START_REF] Bonnans | Numerical Optimization: Theoretical and Practical Aspects[END_REF][START_REF] Hiriart-Urruty | Convex analysis and minimization algorithms[END_REF][START_REF] Nocedal | Numerical optimization[END_REF]. Quasi-Newton methods often lead to significantly faster convergence in practice compared to simpler gradient-based methods for solving smooth optimization problems [START_REF] Schmidt | Projected Newton-type methods in machine learning[END_REF]. Yet, a theoretical analysis of Quasi-Newton methods that explains their impressive empirical behavior is still an open topic. Here, we briefly review the well-known BFGS algorithm in Section 2.1, its limited memory variant [START_REF]Updating quasi-Newton matrices with limited storage[END_REF], and a few recent extensions in Section 2.2. Then, we present earlier works that combine proximal point algorithm and Quasi-Newton methods in Section 2.3.

Quasi-Newton methods for smooth optimization

The most popular Quasi-Newton method is probably BFGS, named after its inventors (Broyden-Fletcher-Goldfarb-Shanno), and its limited variant L-BFGS [START_REF] Nocedal | Numerical optimization[END_REF]. These approaches will be the workhorses of the QNing meta-algorithm in practice. Consider a smooth convex objective f to be minimized, the BFGS method constructs at iteration k a couple (x k , B k ) with the following update:

x k+1 = x k -α k B -1 k ∇f (x k ) and B k+1 = B k - B k s k s k B k s k B k s k + y k y k y k s k , ( 3
)
where α k is a suitable stepsize and

s k = x k+1 -x k , y k = ∇f (x k+1 ) -∇f (x k ).
The matrix B k aims to approximate the Hessian matrix at the iterate x k . When f is strongly convex, the positive definiteness of B k is guaranteed, as well as the condition y k s k > 0, which ensures that (3) is well defined. The stepsize α k is usually determined by a line search strategy. For instance, applying Wolfe's line-search strategy provides linear convergence rate for strongly convex objectives. Moreover, under stronger conditions that the objective f is twice differentiable and its Hessian is Lipschitz continuous, the algorithm can asymptotically achieve superlinear convergence rate [START_REF] Nocedal | Numerical optimization[END_REF]. However, when the dimension d is large, storing the d-by-d matrix B k is infeasible. The limited memory variant L-BFGS [START_REF]Updating quasi-Newton matrices with limited storage[END_REF] overcomes this issue by restricting the matrix B k to be low rank. More precisely, instead of storing the full matrix, a "generating list" of at most l pairs of vectors {(s k i , y k i )} i=0...j is kept in memory. The low rank matrix B k can then be recovered by performing the matrix update recursion in (3) involving all pairs of the generating list. Between iteration k and k + 1, the generating list is incrementally updated, by removing the oldest pair in the list (when j = l) and adding a new one. What makes the approach appealing is the ability of computing the matrix-vector product H k z = B -1 k z with only O(ld) floating point operations for any vector z. This procedure entirely relies on vector-vector product which does not explicitly construct the d-by-d matrix B k or H k . The price to pay is that superlinear convergence becomes out of reach.

L-BFGS is thus appropriate for high-dimensional problems (when d is large), but it still requires computing the full gradient at each iteration, which may be cumbersome in the large sum setting (2). This motivated stochastic counterparts of the Quasi-Newton methods (SQN) [START_REF] Byrd | A stochastic quasi-Newton method for large-scale optimization[END_REF][START_REF] Mokhtari | Global convergence of online limited memory bfgs[END_REF][START_REF] Schraudolph | A stochastic quasi-newton method for online convex optimization[END_REF]. Unfortunately, directly substituting the full gradient ∇f (x k ) by its stochastic counterpart does not lead to a convergent scheme. Instead, the SQN method [START_REF] Byrd | A stochastic quasi-Newton method for large-scale optimization[END_REF] uses updates with sub-sampled Hessian-vector products, which leads to a sublinear convergence rate. Later, a linearly-convergent SQN algorithm is proposed by exploiting a variance reduction scheme [START_REF] Gower | Stochastic block BFGS: Squeezing more curvature out of data[END_REF][START_REF] Moritz | A linearly-convergent stochastic L-BFGS algorithm[END_REF]. However, it is unclear how to extend these techniques into the composite setting.

Quasi-Newton methods for composite optimization

Different approaches have been proposed to extend Quasi-Newton methods to composite optimization problems. A first approach consists in minimizing successive quadratic approximations, also called proximal Quasi-Newton methods [START_REF] Byrd | An inexact successive quadratic approximation method for L-1 regularized optimization[END_REF][START_REF] Ghanbari | Proximal quasi-newton methods for convex optimization[END_REF][START_REF] Lee | Inexact successive quadratic approximation for regularized optimization[END_REF][START_REF] Lee | Proximal Newton-type methods for convex optimization[END_REF][START_REF] Liu | An inexact subsampled proximal newton-type method for large-scale machine learning[END_REF][START_REF] Scheinberg | Practical inexact proximal quasi-Newton method with global complexity analysis[END_REF]. More concretely, a local quadratic approximation q k is minimized at each iteration:

q k (x) f 0 (x k ) + ∇f 0 (x k ), x -x k + 1 2 (x -x k ) T B k (x -x k ) + ψ(x), ( 4
)
where B k is a Hessian approximation based on Quasi-Newton methods. The minimizer of q k provides a descent direction, which is subsequently used to build the next iterate. However, a closed form solution of ( 4) is usually not available since B k changes over the iterations. Thus, one needs to apply an optimization algorithm to approximately solve (4). The composite structure of the subproblem naturally leads to choosing a first-order optimization algorithm, such as randomized coordinate descent algorithms. Then, superlinear complexity becomes out of reach since it requires the subproblems (4) to be solved with "high accuracy" [START_REF] Lee | Proximal Newton-type methods for convex optimization[END_REF].

The global convergence rate of this inexact variant has been for instance analyzed in [START_REF] Scheinberg | Practical inexact proximal quasi-Newton method with global complexity analysis[END_REF], where a sublinear convergence rate is obtained for convex problems; later, the analysis has been extended to strongly convex problems in [START_REF] Liu | An inexact subsampled proximal newton-type method for large-scale machine learning[END_REF], where linear convergence rate is achieved. A second approach of extending Quasi-Newton methods to composite optimization problems is based on smoothing techniques. More precisely, a Quasi-Newton method is applied to a smoothed version of the objective. For instance, one may use the forward-backward envelope [START_REF] Becker | A quasi-newton proximal splitting method[END_REF][START_REF] Stella | Forward-backward quasi-newton methods for nonsmooth optimization problems[END_REF]. The idea is to mimic forward-backward splitting methods and apply Quasi-Newton instead of gradient steps on top of the envelope. Another well known smoothing technique is to apply the Moreau-Yosida regularization [START_REF] Moreau | Fonctions convexes duales et points proximaux dans un espace hilbertien[END_REF][START_REF] Yosida | Functional analysis[END_REF] which gives the smoothed function called Moreau envelope. Then apply Quasi-Newton methods on it leads to the family of variable metric proximal point algorithms [START_REF] Burke | On the superlinear convergence of the variable metric proximal point algorithm using Broyden and BFGS matrix secant updating[END_REF][START_REF] Chen | Proximal quasi-Newton methods for nondifferentiable convex optimization[END_REF][START_REF] Fuentes | Descentwise inexact proximal algorithms for smooth optimization[END_REF][START_REF] Fukushima | A globally and superlinearly convergent algorithm for nonsmooth convex minimization[END_REF]. Our method pursues this line of work by developing a practical inexact variant with global complexity guarantees.

Combining the proximal point algorithm and Quasi-Newton methods

We briefly recall the definition of the Moreau envelope and its basic properties. Definition 1. Given an objective function f and a smoothing parameter κ > 0, the Moreau envelope of f is the function F obtained by performing the infimal convolution

F (x) min z∈R d f (z) + κ 2 z -x 2 . ( 5
)
When f is convex, the sub-problem defined in ( 5) is strongly convex which provides an unique minimizer, called the proximal point of x, which we denote by p(x).

Proposition 1 (Basic properties of the Moreau Envelope

). If f is convex, the Moreau envelope F defined in (5) satisfies 1. F has the same minimum as f , i.e.

min x∈R d F (x) = min x∈R d f (x),
and the solution set of the two above problems coincide with each other.

F is continuously differentiable even when f is not and

∇F (x) = κ(x -p(x)). ( 6
)
Moreover the gradient ∇F is Lipschitz continuous with constant L F = κ.

3. F is convex; moreover, when f is µ-strongly convex with respect to the Euclidean norm, F is µ F -strongly convex with µ F = µκ µ+κ .

4.

F is upper-bounded by f . More precisely, for any x ∈ R d ,

F (x) + 1 2κ ∇F (x) 2 ≤ f (x). (7)
Interestingly, F inherits all the convex properties of f and more importantly it is always continuously differentiable, see [START_REF] Lemaréchal | Practical aspects of the Moreau-Yosida regularization: Theoretical preliminaries[END_REF] for elementary proofs. Moreover, the condition number of F is given by

q = L F µ F = µ + κ µ , ( 8
)
which may be adjusted by the regularization parameter κ. Then, a naive approach to overcome the nonsmoothness of the function f is to transfer the optimization problem to its Moreau envelope F . More concretely, we may apply an optimization algorithm to minimize F and use the obtained solution as a solution to the original problem, since both functions share the same minimizers. This yields the following well-known algorithms.

Proximal point algorithm. Consider the gradient descent method with constant step size 1/L F :

x k+1 = x k - 1 κ ∇F (x k ).
By rewriting the gradient ∇F (x k ) as κ(x k -p(x k )), we obtain the proximal point algorithm [START_REF] Rockafellar | Monotone operators and the proximal point algorithm[END_REF]:

x k+1 = p(x k ) = arg min z∈R d f (z) + κ 2 z -x k 2 . ( 9
)
Accelerated proximal point algorithm. Since gradient descent on F yields the proximal point algorithm, it is natural to apply an accelerated first-order method to get faster convergence. To that effect, Nesterov's algorithm [START_REF] Nesterov | A method of solving a convex programming problem with convergence rate O(1/k 2 )[END_REF] uses a two-stage update, along with a specific extrapolation parameter β k+1 :

x k+1 = y k - 1 κ ∇F (y k ) and y k+1 = x k+1 + β k+1 (x k+1 -x k ),
and, given [START_REF] Bonnans | Numerical Optimization: Theoretical and Practical Aspects[END_REF], we obtain that

x k+1 = p(y k ) and y k+1 = x k+1 + β k+1 (x k+1 -x k ),
This is known as the accelerated proximal point algorithm introduced by Güler [START_REF] Güler | New proximal point algorithms for convex minimization[END_REF], which was recently extended in [START_REF] Lin | A universal catalyst for first-order optimization[END_REF][START_REF] Lin | Catalyst acceleration for first-order convex optimization: from theory to practice[END_REF].

Variable metric proximal point algorithm. One can also apply Quasi-Newton methods on F , which yields

x k+1 = x k -α k B -1 k ∇F (x k ), (10)
where B k is the Hessian approximation of F based on Quasi-Newton methods. This is known as the variable metric proximal point algorithm [START_REF] Burke | On the superlinear convergence of the variable metric proximal point algorithm using Broyden and BFGS matrix secant updating[END_REF][START_REF] Chen | Proximal quasi-Newton methods for nondifferentiable convex optimization[END_REF][START_REF] Fuentes | Descentwise inexact proximal algorithms for smooth optimization[END_REF][START_REF] Fukushima | A globally and superlinearly convergent algorithm for nonsmooth convex minimization[END_REF].

Towards an inexact variable metric proximal point algorithm. Quasi-Newton approaches have been applied after inexact Moreau envelope in various ways [START_REF] Burke | On the superlinear convergence of the variable metric proximal point algorithm using Broyden and BFGS matrix secant updating[END_REF][START_REF] Chen | Proximal quasi-Newton methods for nondifferentiable convex optimization[END_REF][START_REF] Fuentes | Descentwise inexact proximal algorithms for smooth optimization[END_REF][START_REF] Fukushima | A globally and superlinearly convergent algorithm for nonsmooth convex minimization[END_REF]. In particular, it is shown in [START_REF] Chen | Proximal quasi-Newton methods for nondifferentiable convex optimization[END_REF] that if the sub-problems (5) are solved up to high enough accuracy, then the inexact variable metric proximal point algorithm preserves the superlinear convergence rate. However, the complexity for solving the sub-problems with high accuracy is typically not taken into account in these previous works.

In the unrealistic case where p(x k ) can be obtained at no cost, proximal point algorithm can afford much larger step sizes than classical gradient methods, thus are more effective. For instance, when f is strongly convex, the Moreau envelope F can be made arbitrarily well-conditioned by making κ arbitrarily small, according to [START_REF] Byrd | A stochastic quasi-Newton method for large-scale optimization[END_REF]. Then, a single gradient step on F is enough to be arbitrarily close to the optimum. In practice, however, sub-problems are solved only approximately and the complexity of solving the sub-problems is directly related to the smoothing parameter κ. This leaves an important question: how to choose the smoothing parameter κ. A small κ makes the smoothed function F better conditioned, while a large κ is needed to improve the conditioning of the sub-problem [START_REF] Bertsekas | Convex Optimization Algorithms[END_REF].

The main contribution of our paper is to close this gap by providing a global complexity analysis which takes into account the complexity of solving the subproblems. More concretely, in the proposed QNing algorithm, we provide i) a practical stopping criterion for the sub-problems; ii) several warm-start strategies; iii) a simple line-search strategy which guarantees sufficient descent in terms of function value. These three components together yield the global convergence analysis, which allows us to use first-order method as a subproblem solver. Moreover, the global complexity we developed depends on the smoothing parameter κ, which provides some insight about how to practically choose this parameter.

Solving the subproblems with first-order algorithms. In the composite setting, both proximal Quasi-Newton methods and the variable metric proximal point algorithm require solving sub-problems, namely (4) and [START_REF] Bertsekas | Convex Optimization Algorithms[END_REF], respectively. In the general case, when a generic first-order method-e.g., proximal gradient descentis used, our worst-case complexity analysis does not provide a clear winner, and our experiments in Section 5.4 confirm that both approaches perform similarly. However, when it is possible to exploit the specific structure of the sub-problems in one case, but not in the other one, the conclusion may differ.

For instance, when the problem has a finite sum (2) structure, the proximal point algorithm approach leads to sub-problems that can be solved in O(n log(1/ε)) iterations with first-order incremental methods such as SVRG [START_REF] Xiao | A proximal stochastic gradient method with progressive variance reduction[END_REF], SAGA [START_REF] Defazio | Saga: A fast incremental gradient method with support for non-strongly convex composite objectives[END_REF] or MISO [START_REF]Incremental majorization-minimization optimization with application to large-scale machine learning[END_REF], by using the same choice of smoothing parameter κ = 2L/n as Catalyst [START_REF] Lin | Catalyst acceleration for first-order convex optimization: from theory to practice[END_REF]. Assuming that computing a gradient of a function f i and computing the proximal operator of ψ are both feasible in O(d) floating point operations, our approach solves each sub-problem with enough accuracy in Õ(nd) operations. 1 On the other hand, we cannot naively apply SVRG to solve the proximal Quasi-Newton update (4) at the same cost due to the following reasons. First, the variable metric matrix B k does not admit a natural finite sum decomposition. The naive way of writing it into n copies results an increase of computational complexity for evaluating the incremental gradients. More precisely, when B k has rank l, computing a single gradient now requires to compute a matrix-vector product with cost at least O(dl), resulting in l-fold increase per iteration. Second, the previous iteration-complexity O(n log(1/ε)) for solving the sub-problems would require the subproblems to be well-conditioned, i.e. B k (L/n)I, forcing the Quasi-Newton metric to be potentially more isotropic. For these reasons, existing attempts to combine SVRG with Quasi-Newton principles have adopted other directions [START_REF] Gower | Stochastic block BFGS: Squeezing more curvature out of data[END_REF][START_REF] Moritz | A linearly-convergent stochastic L-BFGS algorithm[END_REF].

Algorithm 1 QNing: a Quasi-Newton meta-algorithm input Initial point x 0 in R d ; number of iterations K; smoothing parameter κ > 0; optimization algorithm M; optionally, budget T M for solving the sub-problems. 1: Initialization: (g 0 , F 0 , z 0 ) = ApproxGradient (x 0 , M); H 0 = 1 κ I. 2: for k = 0, . . . , K -1 do

3:

Initialize η k = 1.

4:

Perform the Quasi-Newton step

x test = x k -(η k H k + (1 -η k )H 0 ) g k . 5:
Estimate the gradient and function value of the Approximate Moreau envelope at x test

(g test , F test , z test ) = ApproxGradient (x test , M) . 6: while F test > F k -1 4κ g k 2 do 7:
Decrease the value of the line search parameter η k in [0, 1] and re-evaluate x test .

8:

Re-evaluate (g test , F test , z test ) = ApproxGradient (x test , M) .

9:

end while 10:

Accept the new iterate: (x k+1 , g k+1 , F k+1 , z k+1 ) = (x test , g test , F test , z test ).

11:

Update H k+1 (for example, use L-BFGS update with s k = x k+1 -x k , and y k = g k+1 -g k ). 12: end for output inexact proximal point z K (solution).

Algorithm 2 Generic procedure ApproxGradient input Current point x in R d ; smoothing parameter κ > 0; optionally, budget T M .

1: Compute the approximate proximal mapping using an optimization method M:

z ≈ arg min w∈R d h(w) f (w) + κ 2 w -x 2 , ( 11
)
using one of the following stopping criteria:

• Stop when the approximate solution z satisfies

h(z) -h * ≤ κ 36 z -x 2 . ( 12
)
• Stop when we reach the pre-defined constant budget T M (for instance one pass over the data).

2: Estimate the gradient ∇F (x) of the Moreau envelope using

g = κ(x -z).
output gradient estimate g, objective value estimate F a h(z), proximal mapping estimate z.

The main algorithm

We now discuss the main algorithm components and its main features.

Outer-loop: inexact variable metric proximal point algorithm. We apply variable metric algorithms with a simple line search strategy similar to [START_REF] Scheinberg | Practical inexact proximal quasi-Newton method with global complexity analysis[END_REF] on the Moreau envelope F . Given a positive definite matrix H k and a step size η k in [0, 1], the algorithm computes the update

x k+1 = x k -(η k H k + (1 -η k )H 0 )g k , (LS)
where H 0 = 1 κ I. When η k = 1, the update uses the metric H k , and when η k = 0, it uses an inexact proximal point update

x k+1 = x k -(1/κ)g k .
In other words, when the quality of the metric H k is not good enough, due to the inexactness of the gradients used in its construction, the update is corrected towards a simple proximal point update, whose convergence is well understood when the gradients are inexact.

In order to determine the stepsize η k , we introduce the following descent condition,

F k+1 ≤ F k - 1 4κ g k 2 . ( 13
)
We show that the descent condition ( 13) is always satisfied when η k = 0, thus the finite termination of the line search follows, see Section 4.3 for more details. In our experiments, we observed empirically that the stepsize η k = 1 was almost always selected. In practice, we try the values η k in {1, 1/2, 1/4, 1/8, 0} starting from the largest one and stops whenever condition ( 13) is satisfied.

Example of variable metric algorithm: inexact L-BFGS method.

The L-BFGS rule we consider is the standard one and consists in updating incrementally a generating list of vectors {(s i , y i )} i=1...j , which implicitly defines the L-BFGS matrix. We use here the two-loop recursion detailed in [START_REF] Nocedal | Numerical optimization[END_REF]Algorithm 7.4] and use skipping steps when the condition s i y i > 0 is not satisfied, in order to ensure the positive-definiteness of the L-BFGS matrix H k (see [START_REF] Friedlander | Hybrid deterministic-stochastic methods for data fitting[END_REF]).

Inner-loop: approximate Moreau envelope. The inexactness of our scheme comes from the approximation of the Moreau envelope F and its gradient. The procedure ApproxGradient () calls an minimization algorithm M and apply M to minimize the sub-problem [START_REF] Byrd | Global convergence of a case of quasi-Newton methods on convex problems[END_REF]. When the problem is solved exactly, the function returns the exact values g = ∇F (x), F a = F (x), and z = p(x). However, this is infeasible in practice and we can only expect approximate solutions. In particular, a stopping criterion should be specified. We consider the following variants:

(a) we define an adaptive stopping criterion based on function values and stop M when the approximate solution satisfies the inequality [START_REF] Cartis | Global convergence rate analysis of unconstrained optimization methods based on probabilistic models[END_REF]. In contrast to standard stopping criterion where the accuracy is an absolute constant, our stopping criterion is adaptive since the righthand side of ( 12) also depends on the current iterate z. More detailed theoretical insights will be given in Section 4. Typically, checking whether or not the criterion is satisfied requires computing a duality gap, as in Catalyst [START_REF] Lin | Catalyst acceleration for first-order convex optimization: from theory to practice[END_REF].

(b) using a pre-defined budget T M in terms of number of iterations of the method M, where T M is a constant independent of k.

Note that such an adaptive stopping criterion is relatively classical in the literature of inexact gradient-based methods [START_REF] Byrd | Sample size selection in optimization methods for machine learning[END_REF]. As we will see later in Section 4, when T M is large enough, criterion ( 12) is guaranteed.

Requirements on M. To apply QNing, the optimization method M needs to have linear convergence rates for strongly-convex problems. More precisely, for any strongly-convex objective h, the method M should be able to generate a sequence of iterates (w t ) t≥0 such that

h(w t ) -h * ≤ C M (1 -τ M ) t (h(w 0 ) -h * ) for some constants C M > 0 and 1 > τ M > 0, ( 14
)
where w 0 is the initial point given to M. The notion of linearly-convergent methods extends naturally to non-deterministic methods where ( 14) is satisfied in expectation:

E[h(w t ) -h * ] ≤ C M (1 -τ M ) t (h(w 0 ) -h * ). ( 15
)
The linear convergence condition typically holds for many primal gradient-based optimization techniques, including classical full gradient descent methods, block-coordinate descent algorithms [START_REF] Nesterov | Efficiency of coordinate descent methods on huge-scale optimization problems[END_REF][START_REF] Richtárik | Iteration complexity of randomized block-coordinate descent methods for minimizing a composite function[END_REF], or variance reduced incremental algorithms [START_REF] Defazio | Saga: A fast incremental gradient method with support for non-strongly convex composite objectives[END_REF][START_REF] Schmidt | Minimizing finite sums with the stochastic average gradient[END_REF][START_REF] Xiao | A proximal stochastic gradient method with progressive variance reduction[END_REF]. In particular, our method provides a generic way to combine incremental algorithms with Quasi-Newton methods which are suitable for large-scale optimization problems.

For the simplicity of the presentation, we only consider the deterministic variant ( 14) in the analysis. However, it is possible to show that the same complexity results still hold for non-deterministic methods in expectation, as discussed in Section 4.5. We emphasize that we do not assume any convergence guarantee of M on non-strongly convex problems since our sub-problems are always strongly convex.

Warm starts for the sub-problems. Employing an adequate initialization for solving each sub-problem plays an important rule in our analysis. The warm start strategy we proposed here ensures that the stopping criterion in each subproblem can be achieved in a constant number of iterations: Consider the minimization of a sub-problem min

w∈R d h(w) f (w) + κ 2 w -x 2 .
Then, our warm start strategy depends on the nature of f :

• when f is smooth, we initialize with w 0 = x;

• when f = f 0 + ψ is composite, we initialize with

w 0 = arg min w∈R d f 0 (x) + ∇f 0 (x), w -x + L + κ 2 w -x 2 + ψ(w) ,
which performs an additional proximal step comparing to the smooth case.

Handling composite objective functions. In machine learning or signal processing, convex composite objectives (1) with a non-smooth penalty ψ are typically formulated to encourage solutions with specific characteristics; in particular, the 1 -norm is known to provide sparsity. Smoothing techniques [START_REF] Nesterov | Smooth minimization of non-smooth functions[END_REF] may allow us to solve the optimization problem up to some chosen accuracy, but they provide solutions that do not inherit the properties induced by the non-smoothness of the objective. To illustrate what we mean by this statement, we may consider smoothing the 1 -norm, leading to a solution vector with small coefficients, but not with exact zeroes. When the goal is to perform model selection-that is, understanding which variables are important to explain a phenomenon, exact sparsity is seen as an asset, and optimization techniques dedicated to composite problems such as FISTA [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF] are often preferred (see [START_REF] Mairal | Sparse modeling for image and vision processing[END_REF]).

Then, one might be concerned that our scheme operates on the smoothed objective F , leading to iterates (x k ) k≥0 that may suffer from the above "non-sparse" issue, assuming that ψ is the 1 -norm. Yet, our approach does not directly output the iterates (x k ) k≥0 but their proximal mappings (z k ) k≥0 . In particular, the 1 -regularization is encoded in the proximal mapping [START_REF] Byrd | Global convergence of a case of quasi-Newton methods on convex problems[END_REF], thus the approximate proximal point z k may be sparse. For this reason, our theoretical analysis presented in Section 4 studies the convergence of the sequence (f (z k )) k≥0 to the solution f * .

Convergence and complexity analysis

In this section, we study the convergence of the QNing algorithm-that is, the rate of convergence of the quantities (F (x k ) -F * ) k≥0 and (f (z k ) -f * ) k≥0 , and also the computational complexity due to solving the sub-problems [START_REF] Byrd | Global convergence of a case of quasi-Newton methods on convex problems[END_REF]. We start by stating the main properties of the gradient approximation in Section 4.1. Then, we analyze the convergence of the outer loop algorithm in Section 4.2, and Section 4.3 is devoted to the properties of the line search strategy. After that, we provide the cost of solving the sub-problems in Section 4.4 and derive the global complexity analysis in Section 4.5.

Properties of the gradient approximation

The next lemma is classical and provides approximation guarantees about the quantities returned by the ApproxGradient procedure (Algorithm 2); see [START_REF] Bertsekas | Convex Optimization Algorithms[END_REF][START_REF] Fukushima | A globally and superlinearly convergent algorithm for nonsmooth convex minimization[END_REF]. We recall here the proof for completeness.

Lemma 1 (Approximation quality of the gradient approximation). Consider a vector x in R d , a positive scalar ε and an approximate proximal point z ≈ arg min

w∈R d h(w) f (w) + κ 2 w -x 2 , such that h(z) -h * ≤ ε,
where h * = min w∈R d h(w). As in Algorithm 2, we define the gradient estimate g = κ(x -z) and the function value estimate F a = h(z). Then, the following inequalities hold

F (x) ≤ F a ≤ F (x) + ε, ( 16
)
z -p(x) ≤ 2ε κ , ( 17
)
g -∇F (x) ≤ √ 2κε. ( 18
)
Moreover, F a is related to f by the following relationship

f (z) = F a - 1 2κ g 2 . ( 19
)
Proof. [START_REF] Defazio | Finito: A faster, permutable incremental gradient method for big data problems[END_REF] and [START_REF] Elad | Sparse and Redundant Representations[END_REF] are straightforward by definition of h(z). Since f is convex, the function h is κ-strongly convex, and (17) follows from

κ 2 z -p(x) 2 ≤ h(z) -h(p(x)) = h(z) -h * ≤ ε,
where we recall that p(x) minimizes h. Finally, we obtain (18) from

g -∇F (x) = κ(x -z) -κ(x -p(x)) = κ(p(x) -z),
by using the definitions of g and the property (6).

This lemma allows us to quantify the quality of the gradient and function value approximations, which is crucial to control the error accumulation of inexact proximal point methods. Moreover, the relation [START_REF] Elad | Sparse and Redundant Representations[END_REF] establishes a link between the approximate function value of F and the function value of the original objective f ; as a consequence, it is possible to relate the convergence rate of f from the convergence rate of F . Finally, the following result is a direct consequence of Lemma 1: Lemma 2 (Bounding the exact gradient by its approximation). Consider the same quantities introduced in Lemma 1. Then,

1 2 g 2 -2κε ≤ ∇F (x) 2 ≤ 2( g 2 + 2κε). ( 20
)
Proof. The right-hand side of Eq. ( 20) follows from 18)).

∇F (x) 2 ≤ 2( ∇F (x) -g 2 + g 2 ) ≤ 2(2κε + g 2 ) (from (
Interchanging ∇F (x) and g gives the left-hand side inequality.

Corollary 1. If ε ≤ c κ g 2 with c < 1 4 , then 1 -4c 2 ≤ ∇F (x) 2 g 2 ≤ 2(1 + 2c). ( 21
)
This corollary is important since it allows to replace the unknown exact gradient ∇F (x) by its approximation g , at the cost of a constant factor, as long as the condition ε ≤ c κ g 2 is satisfied.

Convergence analysis of the outer loop

We are now in shape to establish the convergence of the QNing meta-algorithm, without considering yet the cost of solving the sub-problems [START_REF] Byrd | Global convergence of a case of quasi-Newton methods on convex problems[END_REF]. At iteration k, an approximate proximal point is evaluated:

(g k , F k , z k ) = ApproxGradient(x k , M) . ( 22
)
The following lemma characterizes the expected descent in terms of objective function value.

Lemma 3 (Approximate descent property). At iteration k, if the sub-problem ( 22) is solved up to accuracy ε k in the sense of Lemma 1 and the next iterate x k+1 satisfies the descent condition [START_REF] Chang | Libsvm: a library for support vector machines[END_REF], then,

F (x k+1 ) ≤ F (x k ) - 1 8κ ∇F (x k ) 2 + 3 2 ε k . ( 23
)
Proof. From ( 16) and [START_REF] Chang | Libsvm: a library for support vector machines[END_REF], 16) and ( 20)

F (x k+1 ) ≤ F k+1 ≤ F k - 1 4κ g k 2 ≤ F (x k ) + ε k - 1 8κ ∇F (x k ) 2 - ε k 2 (from (
) = F (x k ) - 1 8κ ∇F (x k ) 2 + 3 2 ε k .
This lemma gives us a first intuition about the natural choice of the accuracy ε k , which should be in the same order as ∇F (x k ) 2 . In particular, if

ε k ≤ 1 16κ ∇F (x k ) 2 , ( 24
)
then we have

F (x k+1 ) ≤ F (x k ) - 1 32κ ∇F (x k ) 2 , ( 25
)
which is a typical inequality used for analyzing gradient descent methods. Before presenting the convergence result, we remark that condition (24) cannot be used directly since it requires knowing the exact gradient ∇F (x k ) . A more practical choice consists of replacing it by the approximate gradient.

Lemma 4 (Practical choice of ε k ). The following condition implies inequality [START_REF] Ghadimi | Generalized Uniformly Optimal Methods for Nonlinear Programming[END_REF]:

ε k ≤ 1 36κ g k 2 . ( 26
)
Proof. From Corollary 1, Equation ( 26) implies

g k 2 ≤ 2 1 -4 36 ∇F (x k ) 2 = 9 4 ∇F (x k ) 2 and thus ε k ≤ 1 36κ g k 2 ≤ 1 16κ ∇F (x k ) 2 .
This is the first stopping criterion [START_REF] Cartis | Global convergence rate analysis of unconstrained optimization methods based on probabilistic models[END_REF] in Algorithm 2. Finally, we obtain the following convergence result for strongly convex problems, which is classical in the literature of inexact gradient methods (see Section 4.1 of [START_REF] Byrd | Sample size selection in optimization methods for machine learning[END_REF] for a similar result).

Proposition 2 (Convergence of Algorithm 1, strongly-convex objectives). Assume that f is µstrongly convex. Let (x k ) k≥0 be the sequences generated by Algorithm 1 where the stopping criterion ( 12) is used. Then,

F (x k ) -F * ≤ 1 - 1 16q k (F (x 0 ) -F * ), with q = µ + κ µ .
Proof. The proof follows directly from [START_REF] Ghanbari | Proximal quasi-newton methods for convex optimization[END_REF] and the standard analysis of the gradient descent algorithm for the µ F -strongly convex and L F -smooth function F by remarking that L F = κ and µ F = µκ µ+κ . Corollary 2. Under the conditions of Proposition 2, we have

f (z k ) -f * ≤ 1 - 1 16q k (f (x 0 ) -f * ). ( 27
)
Proof. From ( 19) and ( 26), we have

f (z k ) = F k - 1 2κ g k 2 ≤ F (x k ) + ε k - 1 2κ g k 2 ≤ F (x k ).
Moreover, F (x 0 ) is upper-bounded by f (x 0 ) following [START_REF] Burke | On the superlinear convergence of the variable metric proximal point algorithm using Broyden and BFGS matrix secant updating[END_REF].

It is worth pointing out that our analysis establishes a linear convergence rate whereas one would expect a superlinear convergence rate as for classical variable metric methods. The tradeoff lies in the choice of the accuracy ε k . In order to achieve superlinear convergence, the approximation error ε k needs to decrease superlinearly, as shown in [START_REF] Chen | Proximal quasi-Newton methods for nondifferentiable convex optimization[END_REF]. However, a fast decreasing sequence ε k requires an increasing effort in solving the sub-problems, which will dominate the global complexity. In other words, the global complexity may become worse even though we achieve faster convergence in the outer-loop. This will become clearer when we discuss the inner loop complexity in Section 4.4.

Next, we show that under a bounded level set condition, QNing enjoys the classical sublinear O(1/k) convergence rate when the objective is convex but not strongly convex.

Proposition 3 (Convergence of Algorithm 1 for convex, but not strongly-convex objectives).

Let f be a convex function with bounded level sets. Then, there exists a constant R > 0, which depends on the initialization point x 0 , such that the sequences (x k ) k≥0 and (z k ) k≥0 generated by Algorithm 1 with stopping criterion [START_REF] Cartis | Global convergence rate analysis of unconstrained optimization methods based on probabilistic models[END_REF], satisfies

F (x k ) -F * ≤ 32κR 2 k and f (z k ) -f * ≤ 32κR 2 k .
Proof. We defer the proof and the proper definition of the bounded level set assumption to Appendix A.

So far, we have assumed in our analysis that the iterates satisfy the descent condition [START_REF] Chang | Libsvm: a library for support vector machines[END_REF], which means the line search strategy will always terminate. We prove this is indeed the case in the next section and provide some additional conditions under which a non-zero step size will be selected.

Conditions for non-zero step sizes η k and termination of the line search

At iteration k, a line search is performed on the stepsize η k to find the next iterate

x k+1 = x k -(η k H k + (1 -η k )H 0 )g k ,
such that x k+1 satisfies the descent condition [START_REF] Chang | Libsvm: a library for support vector machines[END_REF]. We first show that the descent condition holds when η k = 0 before giving a more general result. Lemma 5. If the sub-problems are solved up to accuracy ε k ≤ 1 36κ g k 2 , then the descent condition ( 13) holds when η k = 0.

Proof. When η k = 0, x k+1 = x k -1 κ g k = z k . Then, F k+1 ≤ F (x k+1 ) + 1 36κ g k+1 2
(from Eq. ( 16))

≤ F (x k+1 ) + 1 36κ 2 1 -4 36
∇F (x k+1 ) 2 (from Eq. ( 21) )

< F (x k+1 ) + 1 2κ ∇F (z k+1 ) 2 ≤ f (x k+1 ) = f (z k ) (from Eq. ( 7
)) = F k - 1 2κ g k 2 (from Eq. ( 19
)).
Therefore, it is theoretically sound to take the trivial step size η k = 0, which implies the termination of our line search strategy. In other words, the descent condition always holds by taking an inexact gradient step on the Moreau envelope F , which corresponds to the update of the proximal point algorithm. However, the purpose of using variable metric method is to exploit the curvature of the function, which is not the case when η k = 0. Thus, the trivial step size should be only considered as a backup plan and we show in the following some sufficient conditions for taking non-zero step sizes and even stronger, unit step sizes. Lemma 6 (A sufficient condition for satisfying the descent condition (13)). If the sub-problems are solved up to accuracy ε k ≤ 1 36κ g k 2 , then the sufficient condition [START_REF] Chang | Libsvm: a library for support vector machines[END_REF] holds for any

x k+1 = x k -A k g k where A k is a positive definite matrix satisfying 1-α κ I A k 1+α κ I with α ≤ 1 3 .
As a consequence, a line search strategy consisting of finding the largest η k of the form γ i , with i = 1, . . . , +∞ and γ in (0, 1) always terminates in a bounded number of iterations if the sequence of variable metric (H k ) k≥0 is bounded, i.e. there exists 0 < m < M such that for any k, mI H k M I. This is the case for L-BFGS update: Proof of Lemma 6. First, we recall that z k = x k -1 κ g k and we rewrite

F k+1 = F k+1 -F (x k+1 ) E1 + F (x k+1 ) -F (z k ) E2 +F (z k ).
We are going to bound the two error terms E 1 and E 2 by some factors of g k 2 . Noting that the subproblems are solved up to ε k ≤ c κ g k 2 with c = 1 36 , we obtain by construction

E 1 = F k+1 -F (x k+1 ) ≤ ε k+1 ≤ c κ g k+1 2 ≤ 2c (1 -4c)κ ∇F (x k+1 ) 2 , ( 28
)
where the last inequality comes from Corollary 1. Moreover,

∇F (x k+1 ) ≤ ∇F (z k ) + ∇F (x k+1 ) -∇F (z k ) ≤ ∇F (z k ) + κ x k+1 -z k . Since x k+1 -z k = ( 1 κ -A k )g k , we have x k+1 -z k ≤ α κ g k . This implies that ∇F (x k+1 ) ≤ ∇F (z k ) + α g k , ( 29
)
and thus

E 1 ≤ 4c (1 -4c)κ ∇F (z k ) 2 + α 2 g k 2 . ( 30
)
Second, by the κ-smoothness of F , we have

E 2 = F (x k+1 ) -F (z k ) ≤ ∇F (z k ), x k+1 -z k + κ 2 x k+1 -z k 2 ≤ 1 4κ ∇F (z k ) 2 + κ x k+1 -z k 2 + κ 2 x k+1 -z k 2 ≤ 1 4κ ∇F (z k ) 2 + 3α 2 2κ g k 2 , ( 31
)
where the last inequality follows from x k+1 -z k ≤ α κ g k . Combining ( 30) and ( 31) yields

E 1 + E 2 ≤ 4c 1 -4c + 1 4 1 κ ∇F (z k ) 2 + 4c 1 -4c + 3 2 α 2 κ g k 2 . ( 32
)
When c ≤ 1 36 and α ≤ 1 3 , we have

E 1 + E 2 ≤ 1 2κ ∇F (z k ) 2 + 1 4κ g k 2 .
Therefore,

F k+1 ≤ F (z k ) + E 1 + E 2 ≤ F (z k ) + 1 2κ ∇F (z k ) 2 + 1 4κ g k 2 ≤ f (z k ) + 1 4κ g k 2 = F k - 1 4κ g k 2 , ( 33
)
where the last equality follows from [START_REF] Elad | Sparse and Redundant Representations[END_REF], this completes the proof.

Note that in practice, we consider a set of step sizes η k = γ i for i ≤ i max or η k = 0, which naturally upper-bounds the number of line search iterations to i max . More precisely, all experiments performed in this paper use γ = 1/2 and i max = 3. Moreover, we observe that the unit stepsize is very often sufficient for the descent condition to hold, as empirically studied in Appendix C.2.

The following result shows that under a specific assumption on the Moreau envelope F , the unit stepsize is indeed selected when the iterate are close to the optimum. The condition, called Dennis-Moré criterion [START_REF] Dennis | A characterization of superlinear convergence and its application to quasi-newton methods[END_REF], is classical in the literature of Quasi-Newton methods. Even though we cannot formally show that the criterion holds for the Moreau envelope F , since it requires F to be twice continuously differentiable, which is not true in general, see [START_REF] Lemaréchal | Practical aspects of the Moreau-Yosida regularization: Theoretical preliminaries[END_REF], it provides a sufficient condition for the unit step size. Therefore, the lemma below should not be seen as an formal explanation for the choice of step size η k = 1, but simply as a reasonable condition that leads to this choice. Lemma 8 (A sufficient condition for unit stepsize). Assume that f is strongly convex and F is twicecontinuously differentiable with Lipschitz continuous Hessian ∇ 2 F . If the sub-problems are solved up to accuracy ε k ≤ µ 2 128κ(µ+κ) 2 g k 2 and the Dennis-Moré criterion [START_REF] Dennis | A characterization of superlinear convergence and its application to quasi-newton methods[END_REF] is satisfied, i.e.

lim k→∞ (B -1 k -∇ 2 F (x * ) -1 )g k g k = 0, (DM)
where x * is the minimizer of the problem and

B k = H -1 k
is the variable metric matrix, then the descent condition ( 13) is satisfied with η k = 1 when k is large enough.

We remark that the Dennis-Moré criterion we use here is slightly different from the standard one since our criterion is based on approximate gradients g k . If the g k 's are indeed the exact gradients and the variable metric B k are bounded, then our criterion is equivalent to the standard Dennis-Moré criterion. The proof of the lemma is close to that of similar lemmas appearing in the proximal Quasi-Newton literature [START_REF] Lee | Proximal Newton-type methods for convex optimization[END_REF], and is relegated to the appendix. Interestingly, this proof also suggests that a stronger stopping criterion ε k such that ε k = o( g k 2 ) could lead to superlinear convergence. However, such a choice of ε k would significantly increase the complexity for solving the sub-problems, and overall degrade the global complexity.

Complexity analysis of the inner loop

In this section, we evaluate the complexity of solving the sub-problems [START_REF] Byrd | Global convergence of a case of quasi-Newton methods on convex problems[END_REF] up to the desired accuracy using a linearly convergent method M. Our main result is that all sub-problems can be solved in a constant number T M of iterations (in expectation if the method is non-deterministic) using the proposed warm start strategy. Let us consider the sub-problem with an arbitrary prox center x,

min w∈R d h(w) = f (w) + κ 2 w -x 2 . ( 34
)
The number of iterations needed is determined by the ratio between the initial gap h(w 0 ) -h * and the desired accuracy. We are going to bound this ratio by a constant factor.

Lemma 9 (Warm start for primal methods -smooth case). If f is differentiable with L-Lipschitz continuous gradients, we initialize the method M with w 0 = x. Then, we have the guarantee that

h(w 0 ) -h * ≤ L + κ 2κ 2 ∇F (x) 2 . ( 35
)
Proof. Denote by w * the minimizer of h. Then, we have the optimality condition ∇f (w * ) + κ(w * -x) = 0. As a result,

h(w 0 ) -h * = f (x) -f (w * ) + κ 2 w * -x 2 ≤ f (w * ) + ∇f (w * ), x -w * + L 2 x -w * 2 -f (w * ) + κ 2 w * -x 2 = L + κ 2 w * -x 2 = L + κ 2κ 2 ∇F (x) 2 .
The inequality in the above proof relies on the smoothness of f , which does not hold for composite problems. The next lemma addresses this issue.

Lemma 10 (Warm start for primal methods -composite case). Consider the composite optimization problem

f = f 0 + ψ, where f 0 is L-smooth. By initializing with w 0 = arg min w∈R d f 0 (x) + ∇f 0 (x), w -x + L + κ 2 w -x 2 + ψ(w) , ( 36
)
we have,

h(w 0 ) -h * ≤ L + κ 2κ 2 ∇F (x) 2 .
Proof. We use the inequality corresponding to Lemma 2.3 in [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF]: for any w,

h(w) -h(w 0 ) ≥ L 2 w 0 -x 2 + L w 0 -x, x -w , ( 37
)
with L = L + κ. Then, we apply this inequality to w = w * , and

h(w 0 ) -h * ≤ - L 2 w 0 -x 2 -L w 0 -x, x -w * ≤ L 2 x -w * 2 = L + κ 2κ 2 ∇F (x) 2 .
We get an initialization of the same quality in the composite case as in the smooth case, by performing an additional proximal step. It is important to remark that the above analysis do not require strong convexity of f , which allows us to derive the desired inner-loop complexity. Lemma 10. Assume that the optimization method M applied in the inner loop produces a sequence (w t ) t≥0 for each sub-problem [START_REF] Lin | Catalyst acceleration for first-order convex optimization: from theory to practice[END_REF] such that

Proposition 4 (Inner-loop complexity for Algorithm 1). Consider Algorithm 1 with the warm start strategy described in Lemma 9 or in

h(w t ) -h * ≤ C M (1 -τ M ) t (h(w 0 ) -h * ) for some constants C M , τ M > 0. ( 38
)
Then, the stopping criterion ε ≤ 1 72κ g 2 is achieved in at most T M iterations with

T M = 1 τ M log 38C M L + κ κ .
Proof. Consider at iteration k, we apply M to approximate the proximal mapping according to x. With the given T M (which we abbreviate by T ), we have 2 (By Lemma 9 and Lemma 10)

h(w T ) -h * ≤ C M (1 -τ M ) T (h(w 0 ) -h * ) ≤ C M e -τ M T (h(w 0 ) -h * ) ≤ C M e -τ M T L + κ 2κ 2 ∇F (x)
= 1 76κ ∇F (x) 2 ≤ 1 36κ g 2 ,
where the last inequality follows from Lemma 2.

Next, we extend the previous result obtained with deterministic methods M to randomized ones, where linear convergence is only achieved in expectation. The proof is a simple application of Lemma C.1 in [START_REF] Lin | A universal catalyst for first-order optimization[END_REF] (see also [START_REF] Cartis | Global convergence rate analysis of unconstrained optimization methods based on probabilistic models[END_REF] for related results on the expected complexity of randomized algorithms).

Remark 1 (When M is non-deterministic). Assume that the optimization method M applied to each sub-problem (34) produces a sequence (w

t ) t≥0 such that E[h(w t ) -h * ] ≤ C M (1 -τ M ) t (h(w 0 ) -h * ) for some constants C M , τ M > 0.
We define the stopping time T M by

T M = inf t ≥ 1 | h(w t ) -h * ≤ 1 36κ g t 2 , where g t = κ(x -w t ), ( 39
)
which is the random variable corresponding to the minimum number of iterations to guarantee the stopping condition [START_REF] Cartis | Global convergence rate analysis of unconstrained optimization methods based on probabilistic models[END_REF]. Then, when the warm start strategy described in Lemma 9 or in Lemma 10 is applied, the expected number of iterations satisfies

E[T M ] ≤ 1 τ M log 76C M L + κ τ M κ + 1. ( 40
)
Remark 2 (Checking the stopping criterium). It is worth to notice that the stopping criterium [START_REF] Cartis | Global convergence rate analysis of unconstrained optimization methods based on probabilistic models[END_REF], i.e. h(w) -h * ≤ κ 36 w -x 2 , can not be directly checked since h * is unknown. Nevertheless, an upper bound on the optimality gap h(w) -h * is usually available. In particular,

• When f is smooth, which implies h is smooth, we have

h(w) -h * ≤ 1 2(µ + κ) ∇h(w) 2 . ( 41
)
• Otherwise, we can evaluate the Fenchel conjugate function, which is a natural lower bound of h * , see Section D.2.3 in [START_REF]Sparse coding for machine learning, image processing and computer vision[END_REF].

Global complexity of QNing

Finally, we can use the previous results to upper-bound the complexity of the QNing algorithm in terms of iterations of the method M for minimizing f up to ε.

Proposition 5 (Worst-case global complexity for Algorithm 1). Given a linearly-convergent method M satisfying [START_REF] Chen | Proximal quasi-Newton methods for nondifferentiable convex optimization[END_REF], we apply M to solve the sub-problems of Algorithm 1 with the warm start strategy given in Lemma 9 or Lemma 10 up to accuracy

ε k ≤ 1 36κ g k 2 .
Then, the number of iterations of the method M to guarantee the optimality condition

f (z k ) -f * ≤ ε is • for µ-strongly-convex problems: O T M × µ + κ µ log f (x 0 ) -f * ε = O µ + κ τ M µ log f (x 0 ) -f * ε log 38C M L + κ κ .
• for convex problems with bounded level sets:

O T M × 2κR 2 ε = O 2κR 2 τ M ε log 38C M L + κ κ .
Proof. The total number of calls of method M is simply T M times the number of outer-loop iterations times the potential number of line search steps at each iteration (which is hidden in the O(.) notation since this number can be made arbitrarily small).

Remark 3. For non-deterministic methods, applying [START_REF] Mairal | Sparse modeling for image and vision processing[END_REF] yields a global complexity in expectation similar to the previous result with additional constant 2/τ M in the last log factor.

As we shall see, the global complexity of our algorithm is mainly controlled by the smoothing parameter κ. Unfortunately, under the current analysis, our algorithm QNing does not lead to an improved convergence rate in terms of the worst-case complexity bounds. It is worthwhile to underline, though, that this result is not surprising since it is often the case for L-BFGS-type methods, for which an important gap remains between theory and practice. Indeed, L-BFGS often outperforms the vanilla gradient descent method in many practical cases, but never in theory, which turns out to be the bottleneck in our analysis.

We give below the worst-case global complexity of QNing when applied to two optimization methods M of interest. Proposition 5 and its application to the two examples show that, in terms of worse-case complexity, the QNing scheme leaves the convergence rate almost unchanged. Example 1. Consider gradient descent with fixed constant step-size 1/L as the optimization method M.

Directly applying gradient descent (GD) to minimize f requires

O(L/µ log(1/ε))
iterations to achieve ε accuracy. The complexity to achieve the same accuracy with QNing-GD is in the worst case Õ((L + κ)/µ log(1/ε)).

Example 2. Consider the stochastic variance-reduced gradient (SVRG) as the optimization method M. SVRG minimizes f to ε accuracy in

O max n, L µ log 1 ε
iterations in expectation. QNing-SVRG achieves the same result with the worst-case expected complexity

Õ max µ + κ µ n, L + κ µ log 1 ε .
Choice of κ. Minimizing the above worst-case complexity respect to κ suggests that κ should be chosen as small as possible. However, such statement is based on the pessimistic theoretical analysis of L-BFGS-type method, which is not better than standard gradient descent methods. Noting that for smooth functions, L-BFGS method often outperforms Nesterov's accelerated gradient method, it is reasonable to expect they achieve a similar complexity bound. In other words, the choice of κ may be substantially different if one is able to show that L-BFGS-type method enjoys an accelerated convergence rate.

In order to illustrate the difference, we heuristically assume that L-BFGS method enjoys a similar convergence rarte as Nesterov's accelerated gradient method. Then, the global complexity of our algorithm QNing matches the complexity of the related Catalyst acceleration scheme [START_REF] Lin | A universal catalyst for first-order optimization[END_REF], which will be

Õ 1 τ M µ + κ µ log 1 ε ,
for µ-strongly-convex problems. In such case, the complexity of QNing-GD and QNing-SVRG will be

Õ L + κ (µ + κ)µ log(1/ε) and Õ max µ + κ µ n, L + κ (µ + κ)µ log 1 ε ,
which do enjoy acceleration by taking κ = O(L) and κ = O(L/n) respectively. In the following section, we will experiment with this heuristic, as if L-BFGS method enjoys an accelerated convergence rate. More precisely, we will choose the smoothing parameter κ as in the related Catalyst acceleration scheme [START_REF] Lin | A universal catalyst for first-order optimization[END_REF], and we present empirical evidence in support of this heuristic.

Experiments and practical details

In this section, we present the experimental results obtained by applying QNing to several first-order optimization algorithms. We start the section by presenting various benchmarks and practical parametertuning choices. Then, we study the performance of QNing applied to SVRG (Section 5.3) and to the proximal gradient algorithm ISTA (Section 5.4), which reduces to gradient descent (GD) in the smooth case. We demonstrate that QNing can be viewed as an acceleration scheme: by applying QNing to an optimization algorithm M, we achieve better performance than when applying M directly to the problem. Besides, we also compare QNing to existing stochastic variants of L-BFGS algorithm in Section 5.3. Finally, we study the behavior of QNing under different choice of parameters in Section 5.5. The code used for all the experiments is available at https://github.com/hongzhoulin89/Catalyst-QNing/.

Formulations and datasets

We consider three common optimization problems in machine learning and signal processing, including logistic regression, Lasso and linear regression with Elastic-Net regularization. These three formulations all admit the composite finite-sum structure but differ in terms of smoothness and strong-convexness. Specifically, the three formulations are listed below.

• 2 2 -regularized Logistic Regression:

min x∈R d 1 n n i=1 log 1 + exp(-b i a T i x) + µ 2 x 2 ,
which leads to a µ-strongly convex smooth optimization problem.

• 1 -regularized Linear Regression (LASSO):

min x∈R d 1 2n n i=1 (b i -a T i x) 2 + λ x 1 ,
which is convex and non-smooth, but not strongly convex.

• 1 -2 2 -regularized Linear Regression (Elastic-Net):

min x∈R d 1 2n n i=1 (b i -a T i x) 2 + λ x 1 + µ 2 x 2 ,
which is based on the Elastic-Net regularization [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF] leading to non-smooth strongly-convex problems.

For each formulation, we consider a training set (a i , b i ) n i=1 of n data points, where the b i 's are scalars in {-1, +1} and the a i 's are feature vectors in R d . Then, the goal is to fit a linear model x in R d such that the scalar b i can be well predicted by the inner-product a i x, or by its sign. Since we normalize the feature vectors a i , a natural upper-bound on the Lipschitz constant L of the unregularized objective can be easily obtained with L logistic = 1/4, L elastic-net = 1 and L lasso = 1.

In the experiments, we consider relatively ill-conditioned problems with the regularization parameter µ = 1/(100n). The 1 -regularization parameter is set to λ = 1/n for the Elastic-Net formulation; for the Lasso problem, we consider a logarithmic grid 10 i /n, with i = -3, -2, . . . , 3, and we select the parameter λ that provides a sparse optimal solution closest to 10% non-zero coefficients.

Datasets. We consider five standard machine learning datasets with different characteristics in terms of size and dimension, which are described below: The first three data sets are standard machine learning data sets from LIBSVM [START_REF] Chang | Libsvm: a library for support vector machines[END_REF]. We normalize the features, which provides a natural estimate of the Lipschitz constant as mentioned previously. The last two data sets are coming from computer vision applications. MNIST and CIFAR-10 are two image classification data sets involving 10 classes. The feature representation of each image was computed using an unsupervised convolutional kernel network [START_REF]End-to-end kernel learning with supervised convolutional kernel networks[END_REF]. We focus here on the task of classifying class #1 vs. other classes.

Choice of hyper-parameters and variants

We now discuss the choice of default parameters used in the experiments as well as the different variants.

First, to deal with the high-dimensional nature of the data, we systematically use the L-BFGS metric H k and maintain the positive definiteness by skipping updates when necessary (see [START_REF] Friedlander | Hybrid deterministic-stochastic methods for data fitting[END_REF]).

Choice of method M. We apply QNing to proximal SVRG algorithm [START_REF] Xiao | A proximal stochastic gradient method with progressive variance reduction[END_REF] and proximal gradient algorithm. The proximal SVRG algorithm is an incremental algorithm that is able to exploit the finite-sum structure of the objective and can deal with the composite regularization. We also consider the gradient descent algorithm and its proximal variant ISTA, which allows us to perform a comparison with the natural baselines FISTA [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF] and L-BFGS.

Stopping criterion for the inner loop. The default stopping criterion consists of solving each subproblem with accuracy ε

k ≤ 1 36 g k 2 .
Although we have shown that such accuracy is attainable in some constant T = Õ(n) number of iterations for SVRG with the choice κ = L/2n, a natural heuristic proposed in Catalyst [START_REF] Lin | Catalyst acceleration for first-order convex optimization: from theory to practice[END_REF] consists of performing exactly one pass over the data T = n in the inner loop without checking any stopping criterion. In particular, for gradient descent or ISTA, one pass over the data means a single gradient step, because the evaluation of the full gradient requires passing through the entire dataset. When applying QNing to SVRG and ISTA, we call the default algorithm using stopping criterion [START_REF] Ghadimi | Generalized Uniformly Optimal Methods for Nonlinear Programming[END_REF] QNing-SVRG, QNing-ISTA and the one-pass variant QNing-SVRG1, QNing-ISTA1, respectively.

Choice of regularization parameter κ. We choose κ as in the Catalyst algorithm [START_REF] Lin | Catalyst acceleration for first-order convex optimization: from theory to practice[END_REF], which is κ = L for gradient descent/ISTA and κ = L/2n for SVRG. Indeed, convergence of L-BFGS is hard to characterize and its theoretical rate of convergence can be pessimistic as shown in our theoretical analysis. Noting that for smooth functions, L-BFGS often outperforms Nesterov's accelerated gradient method, it is reasonable to expect QNing achieves a similar complexity bound as Catalyst. Later in Section 5.5, we make a comparison between different values of κ to demonstrate the effectiveness of this strategy.

Choice of limited memory parameter l. The default setting is l = 100. We show later in Section 5.5 a comparison with different values to study the influence of this parameter.

Implementation of the line search.

As mentioned earlier, we consider the stepsizes η k in the set {1, 1/2, 1/4, 1/8, 0} and select the largest one that satisfies the descent condition.

Evaluation metric.

For all experiments, we use the number of gradient evaluations as a measure of complexity, assuming this is the computational bottleneck of all methods considered. This is indeed the case here since the L-BFGS step cost O(ld) floating-point operations [START_REF] Nocedal | Numerical optimization[END_REF], whereas evaluating the gradient of the full objective costs O(nd), with l n.

QNing-SVRG for minimizing large sums of functions

We now apply QNing to SVRG and compare different variants.

• SVRG: the Prox-SVRG algorithm of [START_REF] Xiao | A proximal stochastic gradient method with progressive variance reduction[END_REF] with default parameters m = 1 and η = 1/L, where L is the upper-bound on Lipschitz constant of the gradient, as described in the Section 5.1.

• Catalyst-SVRG: The Catalyst meta-algorithm of [START_REF] Lin | Catalyst acceleration for first-order convex optimization: from theory to practice[END_REF] applied to Prox-SVRG, using the variant (C3) that performs best among the different variants of Catalyst.

• L-BFGS/Orthant: Since implementing effectively L-BFGS with a line-search algorithm is a bit involved, we use the implementation of Mark Schmidt2 , which has been widely used in other comparisons [START_REF] Schmidt | Minimizing finite sums with the stochastic average gradient[END_REF]. In particular, the Orthant-wise method follows the algorithm developed in [START_REF] Andrew | Scalable training of L 1 -regularized log-linear models[END_REF]. We use L-BFGS for the logistic regression experiment and the Orthant-wise method [START_REF] Andrew | Scalable training of L 1 -regularized log-linear models[END_REF] for elastic-net and lasso experiments. The limited memory parameter l is set 100.

• QNing-SVRG: the algorithm according to the theory by solving the sub-problems until

ε k ≤ 1 36 g k 2 .
• QNing-SVRG1: the one-pass heuristic.

The result of the comparison is presented in Figure 1 and leads to the conclusions below, showing that QNing-SVRG1 is a safe heuristic, which never decreases the speed of the method SVRG:

• L-BFGS/Orthant is less competitive than other approaches that exploit the sum structure of the objective, except on the dataset real-sim; the difference in performance with the SVRG-based approaches can be important (see dataset alpha).

• QNing-SVRG1 is significantly faster than or on par with SVRG and QNing-SVRG.

• QNing-SVRG is significantly faster than, or on par with, or only slightly slower than SVRG.

• QNing-SVRG1 is significantly faster, or on par with Catalyst-SVRG. This justifies our choice of κ which assumes "a priori" that L-BFGS performs as well as Nesterov's method.

So far, we have shown that applying QNing with SVRG provides a significant speedup compared to the original SVRG algorithm or other acceleration scheme such as Catalyst. Now we compare our algorithm to other variable metric approaches including Proximal L-BFGS [START_REF] Lee | Proximal Newton-type methods for convex optimization[END_REF] and Stochastic L-BFGS [START_REF] Moritz | A linearly-convergent stochastic L-BFGS algorithm[END_REF]:

• Proximal L-BFGS: We apply the Matlab package PNOPT3 implemented by [START_REF] Lee | Proximal Newton-type methods for convex optimization[END_REF]. The sub-problems are solved by the default algorithm up to desired accuracy. We consider one sub-problem as one gradient evaluation in our plot, even though it often requires multiple passes.

• Stochastic L-BFGS (for smooth objectives): We apply the Matlab package StochBFGS4 implemented by [START_REF] Moritz | A linearly-convergent stochastic L-BFGS algorithm[END_REF]. We consider the 'prev' variant which has the best practical performance.

The result of the comparison is presented in Figure 2 and we observe that QNing-SVRG1 is significantly faster than Proximal L-BFGS and Stochastic L-BFGS:

• Proximal L-BFGS often outperforms Orthant-based methods but it is less competitive than QNing.

• Stochastic L-BFGS is sensitive to parameters and data since the variable metric is based on stochastic information which may have high variance. It performs well on dataset covtype but becomes less competitive on other datasets. Moreover, it only applies to smooth problems.

The previous results are complemented by Appendix C.1, which also presents some comparison in terms of outer-loop iterations, regardless of the cost of the inner-loop.

Relative function value (log scale)

Figure 1: Experimental study of the performance of QNing-SVRG for minimizing large sums of functions. We plot the value F (x k )/F * -1 as a function of the number of gradient evaluations, on a logarithmic scale; the optimal value F * is estimated with a duality gap. We plot the value F (x k )/F * -1 as a function of the number of gradient evaluations, on a logarithmic scale; the optimal value F * is estimated with a duality gap.

QNing-ISTA and comparison with L-BFGS

The previous experiments have included a comparison between L-BFGS and approaches that are able to exploit the sum structure of the objective. It is then interesting to study the behavior of QNing when applied to a basic proximal gradient descent algorithm such as ISTA. Specifically, we now consider

• GD/ISTA: the classical proximal gradient descent algorithm ISTA [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF] with back-tracking line-search to automatically adjust the Lipschitz constant of the gradient objective;

• Acc-GD/FISTA: the accelerated variant of ISTA from [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF].

• QNing-ISTA, and QNing-ISTA1, as in the previous section replacing SVRG by GD/ISTA. Figure 3: Experimental study of the performance of QNing-ISTA. We plot the value F (x k )/F * -1 as a function of the number of gradient evaluations, on a logarithmic scale; the optimal value F * is estimated with a duality gap.

The results are reported in Figure 3 and lead to the following conclusions

• L-BFGS is slightly better on average than QNing-ISTA1 for smooth problems, which is not surprising since we use a state-of-the-art implementation with a well-calibrated line search.

• QNing-ISTA1 is always significantly faster than ISTA and QNing-ISTA.

• The QNing-ISTA approaches are significantly faster than FISTA in 12 cases out of 15.

• There is no clear conclusion regarding the performance of the Orthant-wise method vs other approaches.

For three datasets, covtype, alpha and mnist, QNing-ISTA is significantly better than Orthant-wise. However, on the other two datasets, the behavior is different Orthant-wise method outperforms QNing-ISTA.

Experimental study of hyper-parameters l and κ

In this section, we study the influence of the limited memory parameter l and of the regularization parameter κ in QNing. More precisely, we start with the parameter l and try the method QNing-SVRG1 with the values l = 1, 2, 5, 10, 20, 100. Note that all previous experiments were conducted with l = 100, which is the most expensive in terms of memory and computational cost for the L-BFGS step. The results are presented in Figure 4. Interestingly, the experiment suggests that having a large value for l is not necessarily the best choice, especially for composite problems where the solution is sparse, where l = 10 seems to perform reasonably well.

The next experiment consists of studying the robustness of QNing to the smoothing parameter κ. We present in Figure 5 an experiment by trying the values κ = 10 i κ 0 , for i = -3, -2, . . . , 2, 3, where κ 0 is the default parameter that we used in the previous experiments. The conclusion is clear: QNing clearly slows down when using a larger smoothing parameter than κ 0 , but it is robust to small values of κ (and in fact it even performs better for smaller values than κ 0 ).

Discussions and concluding remarks

A few questions naturally arise regarding the QNing scheme: one may wonder whether or not our convergence rates may be improved, or if the Moreau envelope could be replaced by another smoothing technique. In this section, we discuss these two points and present concluding remarks.

Discussion of convergence rates

In this paper, we have established the linear convergence of QNing for strongly convex objectives when sub-problems are solved with enough accuracy. Since QNing uses Quasi-Newton steps, one might have expected a superlinear convergence rate as several Quasi-Newton algorithms often enjoy [START_REF] Byrd | Global convergence of a case of quasi-Newton methods on convex problems[END_REF]. The situation is as follows. Consider the BFGS algorithm (without limited memory), as shown in [START_REF] Chen | Proximal quasi-Newton methods for nondifferentiable convex optimization[END_REF], if the sequence (ε k ) k≥0 decreases super-linearly, then, it is possible to design a scheme similar to QNing that indeed enjoys a super-linear convergence rate. There is a major downside though: a super-linearly decreasing sequence (ε k ) k≥0 implies an exponentially growing number of iterations in the inner-loops, which will degrade the global complexity of the algorithm. This issue makes the approach proposed in [START_REF] Chen | Proximal quasi-Newton methods for nondifferentiable convex optimization[END_REF] impractical.

Another potential strategy for obtaining a faster convergence rate consists in interleaving a Nesterov-type extrapolation step in the QNing algorithm. Indeed, the convergence rate of QNing scales linearly in the condition number µ F /L F , which suggests that a faster convergence rate could be obtained using a Nesterovtype acceleration scheme. Empirically, we did not observe any benefit of such a strategy, probably because of the pessimistic nature of the convergence rates that are typically obtained for Quasi-Newton approaches based on L-BFGS. Obtaining a linear convergence rate for an L-BFGS algorithm is still an important sanity check, but to the best of our knowledge, the gap in performance between these worst-case rates and practice has always been huge for this class of algorithms. Figure 5: Experimental study of influence of the smoothing parameter κ for QNing-SVRG1. κ 0 denotes the default choice used in the previous experiments. We plot the value F (x k )/F * -1 as a function of the number of gradient evaluations, on a logarithmic scale; the optimal value F * is estimated with a duality gap.

Other types of smoothing

The Moreau envelope we considered is a particular instance of infimal convolution smoothing [START_REF] Beck | Smoothing and first order methods: A unified framework[END_REF], whose family also encompasses the so-called Nesterov smoothing [START_REF] Beck | Smoothing and first order methods: A unified framework[END_REF]. Other ways to smooth a function include randomization techniques [START_REF] Duchi | Randomized smoothing for stochastic optimization[END_REF] or specific strategies tailored for the objective at hand.

One of the main purposes of applying the Moreau envelope is to provide a better conditioning. As recalled in Proposition 1, the gradient of the smoothed function F is κ-Lipschitz continuous regardless of whether the original function is continuously differentiable or not. Furthermore, the conditioning of F is improved with respect to the original function, with a condition number depending on the amount of smoothing. As highlighted in [START_REF] Beck | Smoothing and first order methods: A unified framework[END_REF], this property is also shared by other types of infimal convolutions. Therefore, QNing could potentially be extended to such types of smoothing in place of the Moreau envelope. A major advantage of our approach, though, is its outstanding simplicity.

Concluding remarks

To conclude, we have proposed a generic mechanism, QNing, to accelerate existing first-order optimization algorithms with Quasi-Newton-type methods. QNing's main features are the compatibility with variable metric update rule and composite optimization. Its ability of combining with incremental approaches makes it a promising tool for solving large-scale machine learning problems. A few questions remain however open regarding the use of the method in a pure stochastic optimization setting, and the gap in performance between worst-case convergence analysis and practice is significant. We are planning to address the first question about stochastic optimization in future work; the second question is unfortunately difficult and is probably one of the main open questions in the literature about L-BFGS methods.

Lemma 11. If f has bounded level sets, then its Moreau envelope F has bounded level sets as well.

Proof. First, from Proposition 1, the minimum of F is attained at x * . Next, we reformulate the bounded level set property by contraposition: for any x, there exists R x > 0 such that ∀y ∈ R d s.t. y -x * > R x then f (y) > f (x). This completes the proof.

We are now in shape to prove the proposition.

Proof. From (25), we have

F (x k+1 ) ≤ F (x k ) - 1 32κ ∇F (x k ) 2 .
Thus F (x k ) is decreasing. From the bounded level set property of F , there exists R > 0 such that x k -x * ≤ R for any k. By the convexity of F , we have

F (x k ) -F * ≤ ∇F (x k ), x k -x * ≤ ∇F (x k ) x k -x * ≤ R ∇F (x k ) .
Therefore,

F (x k+1 ) -F * ≤ F (x k ) -F * - 1 32κ ∇F (x k ) 2 ≤ F (x k ) -F * - (F (x k ) -F * ) 2 32κR 2 .
Let us define r k f (x k ) -f * . Thus,

1 r k+1 ≥ 1 r k (1 -r k 32κR 2 ) ≥ 1 r k 1 + r k 32κR 2 = 1 r k + 1 32κR 2 .
Then, after exploiting the telescoping sum,

1 r k+1 ≥ 1 r 0 + k + 1 32κR 2 ≥ k + 1 32κR 2 .

B Proof of Lemma 8

Proof. Let us denote δ k = -B -1 k g k and let the subproblems solved to accuracy ε k ≤ c κ g k 2 . We show that when c ≤ µ 2 128(µ+κ) 2 , the following two inequalities hold:

F (x k + δ k ) ≤ F (x k ) - 3 8κ g k 2 + o( g k 2 ), (42)
and

F k+1 ≤ F (x k + δ k ) + 1 16κ g k 2 + o( g k 2 ). ( 43
)
Then summing up the above inequalities yields

F k+1 ≤ F (x k ) - 5 16κ g k 2 + o( g k 2 ) ≤ F k - 1 4κ g k 2 ,
where the last inequality holds since F (x k ) ≤ F k and o( g k 2 ) ≤ 1 4κ g k 2 when k is large enough. This is the desired descent condition [START_REF] Chang | Libsvm: a library for support vector machines[END_REF].

We first prove [START_REF] Mokhtari | Global convergence of online limited memory bfgs[END_REF] which relies on the somoothness and Lipschitz Hessian assumption of F . More concretely,

F (x k + δ k ) -F (x k ) ≤ ∇F (x k ) T δ k + 1 2 δ T k ∇ 2 F (x k )δ k + L 2 6 δ k 3 = (∇F (x k ) -g k ) T δ k + g T k δ k + 1 2 δ T k (∇ 2 F (x k ) -B k )δ k + 1 2 δ T k B k δ k =-1 2 g T k δ k + L 2 6 δ k 3 = 1 2 g T k δ k E1 + (∇F (x k ) -g k ) T δ k E2 + 1 2 δ T k (∇ 2 F (x k ) -B k )δ k E3 + L 2 6 δ k 3 E4
.

We are going upper bound each term one by one. First,

E 1 = 1 2 g T k δ k = - 1 2 g k B -1 k g k = - 1 2 g k ∇ 2 F (x * ) -1 g k - 1 2 g k (B -1 k -∇ 2 F (x * ) -1 )g k ≤ - 1 2κ g k 2 + o( g k 2 ),
where the last inequality uses (DM) and the κ-smoothness of F which implies ∇ 2 F (x * ) κI. Second,

E 2 = (∇F (x k ) -g k ) T δ k ≤ ∇F (x k ) -g k δ k ≤ √ 2c g k B -1 k g k (from (18)) ≤ √ 2c g k ∇ 2 F (x * ) -1 g k + B -1 k -∇ 2 F (x * ) -1 g k ≤ √ 2c 1 µ F g k 2 + o( g k 2 ) = 1 8κ g k 2 + o( g k 2 ).
Third,

E 3 = 1 2 δ T k (∇ 2 F (x k ) -B k )δ k ≤ 1 2 δ k (∇ 2 F (x k ) -B k )δ k ≤ 1 2 δ k (∇ 2 F (x k ) -∇ 2 F (x * ))δ k + (∇ 2 F (x * ) -B k )δ k ≤ L 2 2 x k -x * δ k 2 + ∇ 2 F (x * ) (B -1 k -∇ 2 F (x * ) -1 )g k = o( g k 2 )
, where the last line comes from (DM) and the fact that x k -x * → 0. Last, since

δ k = -∇ 2 F (x * ) -1 g k =O( g k ) + (∇ 2 F (x * ) -1 -B -1 k )g k =o( g k ) by Dennis-Moré condition
and g k → 0, we have

E 4 = L 2 6 δ k 3 = o( g k 2 ).
Summing up above four inequalities yields [START_REF] Mokhtari | Global convergence of online limited memory bfgs[END_REF]. Next, we prove the other desired inequality [START_REF] Moreau | Fonctions convexes duales et points proximaux dans un espace hilbertien[END_REF]. The main effort is to bound g k+1 by a constant factor times g k . From the inexactness of the subproblem, we have

F k+1 ≤ F (x k+1 ) + c κ g k+1 2 ≤ F (x k+1 ) + 2c (1 -4c)κ ∇F (x k+1 ) 2
Moreover,

∇F (x k+1 ) -∇F (x k ) -∇ 2 F (x * )(x k+1 -x k ) = 1 0 ∇ 2 F (x k + τ (x k+1 -x k )) -∇ 2 F (x * ) (x k+1 -x k )dτ Therefore, ∇F (x k+1 ) -∇F (x k ) -∇ 2 F (x * )(x k+1 -x k ) ≤ max { x k -x * , x k+1 -x * } x k+1 -x k = o( g k ). ∇F (x k+1 ) ≤ ∇F (x k ) + ∇ 2 F (x * )(x k+1 -x k ) + o( g k ) ≤ ∇F (x k ) -g k + g k + ∇ 2 F (x * )(x k+1 -x k ) =o( g k ) by Dennis-Moré condition +o( g k ) ≤ √ 2c g k + o( g k ).
As a result,

F k+1 ≤ F (x k+1 ) + 4c 2 (1 -4c)κ ∇g k 2 + o( g k 2 ) ≤ F (x k+1 ) + 1 16κ ∇g k 2 + o( g k 2 ) when c ≤ 1 16 .
This completes the proof.

complexity of the sub-problems into account, QNing-SVRG never outperforms QNing-SVRG1. This suggests that it is not beneficial to solve the sub-problem up to high accuracy as long as the algorithm converge.

C.2 Empirical frequency of choosing the unit stepsize

In this section, we evaluate how often the unit stepsize is taken in the line search. When the unit stepsize is taken, the variable metric step provides sufficient decrease, which is the key for acceleration. The statistics of QNing-SVRG1 (one-pass variant) and QNing-SVRG (the sub-problems are solved until the stopping criteria ( 12) is satisfed ) are given in Table 1 and Table 2, respectively. As we can see, for most of the iterations (> 90%), the unit stepsize is taken. The first column is in the form N/D, where N is the number of times over the iterations the unit stepsize was picked and D is the total number of iterations. The total number of iterations D varies a lot since we stop our algorithm as soon as the relative function gap is smaller than 10 -10 or the maximum number of iterations 100 is reached. It implicitly indicates how easy the problem is. The setting are the same as in Table 1.

Lemma 7 (

 7 Boundedness of L-BFGS metric matrix, see Chap 8,9 of [51]). The variable metric matrices (B k ) k constructed by the L-BFGS rule are positive definite and bounded.

1 Relative

 1

1 Relative

 1

1 Relative

 1

1 Relative

 1

1 Relative

 1

1 Relative

 1

Figure 2 :

 2 Figure 2: Comparison to proximal L-BFGS and Stochastic L-BFGS. We plot the value F (x k )/F * -1 as a function of the number of gradient evaluations, on a logarithmic scale; the optimal value F * is estimated with a duality gap.

Figure 4 :

 4 Figure4: Experimental study of influence of the limited-memory parameter l for QNing-SVRG1. We plot the value F (x k )/F * -1 as a function of the number of gradient evaluations, on a logarithmic scale; the optimal value F * is estimated with a duality gap.

2 .•

 2 Given x in R d , we show that∀y ∈ R d s.t. y -x * > 2(f (x) -f * ) κ + R x then F (y) > F (x).Let y satisfies the above inequality, by definition,F (y) = f (p(y)) + κ 2 p(y) -yFrom the triangle inequality,y -p(y) + p(y) -x * ≥ y -x * > 2(f (x) -f * ) κ + R x .Then either y -p(y) > 2(f (x)-f * )κ or p(y) -x * > R x . • If y -p(y) > 2(f (x)-f * ) κ ,thenF (y) = f (p(y)) + κ 2 p(y) -y 2 > f (p(y)) + f (x) -f * ≥ f (x) ≥ F (x). If p(y) -x * > R x , then F (y) = f (p(y)) + κ 2 p(y) -y 2 ≥ f (p(y)) > f (x) ≥ F (x).

Table 1 :

 1 Relative frequency of picking the unit stepsize η k = 1 of QNing-SVRG1

		Logistic	Elastic-net	Lasso
	covtype 24/27 89% 54/56 96%	19/21	90%
	alpha	8/8	100%	6/6	100%	6/6	100%
	real-sim 60/60 100% 71/76 93%	14/14	100%
	mnist	53/53 100% 80/80 100% 100/100 100%
	cifar-10 58/58 100% 75/75 100%	42/44	95%

Table 2 :

 2 Relative frequency of picking the unit stepsize η k = 1 of QNing-SVRG

		Logistic	Elastic-net	Lasso
	covtype 18/20 90% 23/25 92% 16/16 100%
	alpha	6/6	100%	4/4	100%	3/3	100%
	real-sim 27/27 100% 20/23 87%	8/8	100%
	mnist	27/27 100% 28/28 100% 28/28 100%
	cifar-10 25/25 100% 29/29 100% 31/31 100%

QNing: a Quasi-Newton meta-algorithmWe now present the QNing method in Algorithm 1, which consists of applying variable metric algorithms on the smoothed objective F with inexact gradients. Each gradient approximation is the result of a minimization problem tackled with the algorithm M, used as a sub-routine. The outer loop of the algorithm performs Quasi-Newton updates. The method M can be any algorithm of the user's choice, as long as it enjoys linear convergence rate for strongly convex problems. More technical details are given in Section 3.1.[START_REF] Andrew | Scalable training of L 1 -regularized log-linear models[END_REF] The notation Õ hides logarithmic quantities.

available here http://www.cs.ubc.ca/~schmidtm/Software/minFunc.html

available here https://web.stanford.edu/group/SOL/software/pnopt

available here https://perso.telecom-paristech.fr/rgower/software.html

Acknowledgements

The authors would like to thank the editor and the reviewers for their constructive and detailed comments. This work was supported by the ERC grant SOLARIS (number 714381), a grant from ANR (MACARON project ANR-14-CE23-0003-01), and the program "Learning in Machines and Brains" (CIFAR).

A Proof of Proposition 3

First, we show that the Moreau envelope F inherits the bounded level set property from f . Definition 2. We say that a convex function f has bounded level sets if f attains its minimum at x * in R d and for any x, there exists R x > 0 such that

C Additional experiments

In this section, we provide additional experimental results including experimental comparisons in terms of outer loop iterations, and an empirical study regarding the choice of the unit step size η k = 1.

C.1 Comparisons in terms of outer-loop iterations

In the main paper, we have used the number of gradient evaluations as a natural measure of complexity.

Here, we also provide a comparison in terms of outer-loop iterations, which does not take into account the complexity of solving the sub-problems. While interesting, the comparison artificially gives an advantage to the stopping criteria [START_REF] Cartis | Global convergence rate analysis of unconstrained optimization methods based on probabilistic models[END_REF] since achieving it usually requires multiple passes. The result of the comparison is presented in Figure 6. We observe that the theoretical grounded variant QNing-SVRG always outperform the one-pass heuristic QNing-SVRG1. This is not surprising since the sub-problems are solved more accurately in the theoretical grounded variant. However, once we take the