Contracting Proximal Methods for Smooth Convex Optimization

In this paper, we propose new accelerated methods for smooth convex optimization, called contracting proximal methods. At every step of these methods, we need to minimize a contracted version of the objective function augmented by a regularization term in the form of Bregman divergence. We provide global convergence analysis for a general scheme admitting inexactness in solving the auxiliary subproblem. In the case of using for this purpose high-order tensor methods, we demonstrate an acceleration effect for both convex and uniformly convex composite objective functions. Thus, our construction explains acceleration for methods of any order starting from one. The augmentation of the number of calls of oracle due to computing the contracted proximal steps is limited by the logarithmic factor in the worst-case complexity bound.

1. Introduction. One of the classical iterative methods in theoretical optimization is the proximal point algorithm [24]. This method, as applied to minimizing a convex function f : dom f → R, consists of solving at each iteration the following subproblem: where · is the standard Euclidean norm, and {a k } k≥1 is a sequence of positive coefficients. In general, we can hope only to use an inexact solution of the subproblem (1.1) (see [10,27,26] for the convergence analysis). An important observation is that the regularized objective in (1.1) is strongly convex. Therefore, we can hope that computing an (inexact) proximal step is usually simpler than solving the initial problem.
For a function f ∈ F 1,1 L (convex differentiable functions with Lipschitz continuous gradients), we can set all values of the coefficients a k equal to a positive constant. This gives a global sublinear rate of convergence of the iterations (1.1) in functional residual of the order O(1/k). This rate is the same rate as that of the gradient method [21].
For the same class of functions, we can get a faster rate of convergence of the order O(1/k 2 ) using the accelerated gradient method [18]. This is the best possible rate achievable for the first-order black-box optimization on F 1,1 variant of the proximal point algorithm with the optimal rate of convergence was proposed in [11] (see also [25,13,14,12] for extensions and some applications).
In this paper, we present a new family of proximal-type algorithms for smooth convex optimization called contracting proximal methods, which includes an accelerated algorithm from [11] as a particular case and provides a systematic way for constructing faster proximal accelerated methods for high-order optimization. Thus, for the class of convex functions, whose p-th derivative is Lipschitz continuous (p ≥ 1), our new methods achieve the O(1/k p+1 )-rate of convergence for the outer proximal iterations, while the inner subproblems can be efficiently solved up to desired accuracy by the high-order tensor methods [22]. Note that this rate can also be achieved by a direct acceleration scheme, utilizing the notion of estimating sequences [2,22]. It can be improved up to the level O(1/k 3p+1 2 ) by using a special line-search on each iteration [16,7]. The latter rate was shown to be the optimal one [1].
The main difference between contracting proximal methods and the classical approach (1.1) consists in employing the contracted objective function (which provides the methods with their name) and the Bregman divergence (notation β d (x; y)) instead of the usual Euclidean norm. The exact form of our method is very simple: Thus, we use a sequence of auxiliary points {v k } k≥0 and the scaling coefficients Let us illustrate the basic idea behind this construction by the simplest Euclidean setting, when β d (x; y) ≡ 1 2 x − y 2 . We are going to ensure at each iteration k ≥ 0 the following condition: We can propagate inequality (1.3) to the next iteration by a trivial observation: where the last inequality is due to convexity of the objective. Note that the first step of contracting proximal method (1.2) is defined exactly as follows: Hence, by strong convexity of h k+1 (·), we finally justify that Thus, for the Euclidean setting, iteration (1.2) immediately results in the convergence guarantee (1.4). However, we are still free in the choice of coefficients {a k } k≥1 . The only reason for bounding their growth consists in keeping the complexity of the optimization problem (1.5) at an acceptable level. 1 For f ∈ F 1,1 L , the recommended choice of a k+1 corresponds to the quadratic equation [18]: It is easy to see that this choice results in the optimal O(1/k 2 )-rate of convergence for the method. On the other hand, it makes the condition number of the problem (1.5) equal to an absolute constant. Let us assume for simplicity that f is two times continuously differentiable. Then, in view of the presence of the regularization term, On the other hand, 2I.
Hence, we are able to solve the problem (1.5) very efficiently by a usual gradient method (see the details in section 4). It is remarkable that exactly the same reasoning justifies the accelerated versions of all high-order tensor methods (p ≥ 2). The only difference consists in the degree of the proximal term, which must be compatible with the order of optimization scheme used for solving the problem (1.5).
Our first-order contracting proximal method for the Euclidean setting (described above) produces the same sequence of points as the accelerated Proximal Point Algorithm from [11]. However, now we can also employ the Bregman divergence, which sometimes is more suitable to the topology of our function and ensures faster convergence.
The rest of the paper is organized as follows. Section 2 introduces notation used throughout the paper and describes our problem of interest in the composite form. We also give a definition of Bregman divergence and mention some of its properties.
In section 3, we introduce a general contracting proximal method (formulated as Algorithm 1). We present its convergence analysis for a problem in composite form and arbitrary Bregman divergence. We study both convex and strongly convex cases under inexactness in proximal steps. Theorem 3.2 specifies how the parameters of the algorithm and inner accuracy affect the convergence rate.
In section 4, we discuss implementation of one iteration of our method, under the assumption that the p-th derivative (p ≥ 1) of the smooth part of the objective is Lipschitz continuous. We present a fully defined optimization scheme (Algorithm 2), with incorporated steps of the tensor method of a certain degree. The resulting algorithm achieves the accelerated rate of convergence, with an additional logarithmic factor for the number of total oracle calls. The final complexity estimate for this scheme is given by Theorems 4.6 and 4.7.
Section 5 contains numerical experiments. Section 6 has some final remarks.

Notation.
In what follows, we denote by E a finite-dimensional real vector space and by E * its dual space, which is a space of linear functions on E. The value of function s ∈ E * at point x ∈ E is denoted by s, x .
Let us fix some arbitrary (possibly non-Euclidean) norm · on space E and define the dual norm · * on E * in the standard way: For a smooth function f , its gradient at point x is denoted by ∇f (x), and its Hessian is Higher derivatives are denoted as D p f (x) [·], which are p-linear symmetric forms on E, and the norm is induced: For convex but not necessary differentiable function ψ, we denote by ∂ψ(x) ⊂ E * its subdifferential at point x ∈ dom ψ.
Our goal is to solve the following composite minimization problem: where f is several times differentiable on its open domain convex function, with some reasonable assumptions on the growth of its derivatives (for example, that its p-th derivative is Lipschitz continuous for some p ≥ 1), and ψ : E → R ∪ {+∞} is a proper closed convex function, which we assume to be simple, but possibly nondifferentiable, with dom ψ ⊂ dom f . We also assume that solution x * ∈ dom F of problem (2.1) does exist, denoting F * = F (x * ). Let us fix arbitrary differentiable strictly convex function d : dom ψ → R, which we call theprox function. Then, we denote by β d (x; y) the corresponding Bregman divergence, centered at x: We say that function d is uniformly convex of degree p + 1 (with respect to the norm · ) with constant σ p+1 (d) > 0 if it holds for all x, y ∈ dom d: x − y p+1 .
The main example, which naturally appears in tensor methods (see [22]) and which we use in section 4, is the following prox function.
for some p ≥ 1. For the Euclidean norm (when x ≡ Bx, x 1/2 for a fixed postivedefinite linear operator B = B * 0) this prox function is uniformly convex of degree p + 1 with constant 2 1−p (see Lemma 5 in [5]), so it holds that For more examples of available prox functions see [3,15]. The definition of Bregman divergence can be extended onto nondifferentiable function ψ by specifying a particular subgradient ψ (x) ∈ ∂ψ(x): However, we will use simpler notation β ψ (x; y) if no ambiguity arises.
We say that function ψ is strongly convex with respect to d (see [28,3,15]) with constant σ d (ψ) > 0 if it holds for all x, y ∈ dom ψ and for all ψ (x) ∈ ∂ψ(x) that Inequality (2.4) always holds with σ d (ψ) = 0 just by convexity. An interesting illustration of this concept is given by a regularized Taylor polynomial of degree 3 for convex function (see [22]).
Example 2.2. Let f : dom f → R be convex, with Lipschitz continuous third derivative: Denote by Ω 3 (f, x; y) its Taylor approximation of degree 3 around some fixed point x, and consider its regularization of degree 4, with some τ > 1: Then, for the Euclidean norm, the function g(·) is strongly convex with respect to the following prox function (see Lemma 4 in [22]): Let us summarize some basic properties of Bregman divergence, which follow directly from its definition. For any pair f 1 , f 2 of convex functions and all x, y ∈ dom (f 1 + f 2 ) we have For any linear function (x) = a + g, x we have Therefore, from (2.5) and (2.6) we conclude that when f (y) = β d (z; y) for some fixed z. Now, consider the following simple but general construction, which we use in a core of our analysis. Let h be a regularized composite objective: h(y) = g(y) + aψ(y) + γβ d (z; y), a, γ ≥ 0, where g and ψ are arbitrary closed convex functions, and ψ is strongly convex with respect to d with some constant σ d (ψ) ≥ 0. Then we have, for every x, y ∈ dom h and every h (x) ∈ ∂h(x) In particular, for the exact minimum T = argmin y∈E h(y), we have 3. Contracting proximal method. In our general scheme, we are going to maintain the following inequality, for every x ∈ dom ψ and k ≥ 0: and {γ k } k≥0 is a sequences of nondecreasing proximal coefficients, We would prefer to have functions C k (x) as big as possible. Thus, if it happens to be C k (x * ) ≥ 0 for all k ≥ 1, then from (3.1) we have a convergence guarantee, and the rate of convergence is determined by the growth of coefficients A k toward infinity. However, generally C k (x) may have arbitrary sign.
Let us discus a simple possibility for propagating relation (3.1) to the next iteration.
where the last inequality is due to the convexity of f . Let us consider a contracted objective with a regularizer from the last step: This function is strongly convex with respect to d(·) with parameter If we are able to compute the exact minimum then by (2.9) we see that And it is natural to set v k+1 = T and Thus we would obtain guarantee (3.1) for the next step, with Now, instead of computing the exact minimum (3.5), let us relax v k+1 ∈ dom ψ to be a point with a small norm of subgradient: Note that condition (3.7) can be easily verified algorithmically since in composite setting we are able to compute points with a small subgradient of h k+1 (see [22]). Thus, we come to the following general scheme.

Algorithm 1: Contracting Proximal Method
Require: 2: Denote contracted objective with regularizer: At this moment, we need one additional assumption. It relates the dual norm · * (used at step 4) with the Bregman divergence β d (v; x).
Let us write down the convergence guarantees of the method.
Let Assumption 3.1 hold. Then for Algorithm 1 at all iterations k ≥ 0 we have Proof. First, let us ensure by induction in k ≥ 0 the following inequality: where s i ∈ ∂h i (v i ). It is obviously true for k = 0. Let it hold for some k ≥ 0 and consider the next. Note that (3.10) is exactly (3.1) with Therefore, we have This is (3.10) for the next step. Now, plugging x ≡ x * into (3.10) and taking into account the nonnegativity of all terms in the left-hand side, we get Now, we need to estimate the right-hand side from above. Using the uniform convexity (2.2), we conclude that for every k ≥ 0 In order to finish the proof, it is enough to bound from above the value α k , for which we have the following recurrence: Dividing both sides by α 1 p+1 k and using the monotonicity of this sequence, we get Finally, from the last inequality we obtain which is the right-hand side of (3.8).
We see that accuracies δ k for subgradients of the subproblems appear in (3.9) in an additive form, weighted by the coefficients γ − 1 p+1 k . They should be chosen in a way making the right-hand side of (3.8) small enough. Let us consider the simplest case, when all δ k are the same. Corollary 3.3. Let δ k = δ > 0 for all k ≥ 1. Assume that the coefficients A k grow sublinearly: with some constant c > 0. Then for every Proof. Indeed, Summing up these two inequalities we obtain (3.14).
Proof. Indeed, Now, note that the following inequality holds for all x ≥ 0: And we obtain Estimates (3.13) and (3.17) show that the bound for the inner accuracy δ has a reasonable dependency on the absolute accuracy ε required for the initial problem (2.1). Thus, in both cases, in step 4 of the algorithm we need to find a point v k+1 with subgradient s ∈ ∂h k+1 (v k+1 ): This is a reachable goal, especially for methods minimizing h k+1 (·) with a linear rate of convergence.
In practice, it may be reasonable not to use very small inner accuracy on a first stage but to decrease it over the iterations. Then, the following simple choice of {δ k } k≥0 can work.
Corollary 3.5. Let us define δ k ≡ c k s with fixed absolute constants c > 0 and s > 1. Then, Therefore, we have

Application of tensor methods.
In this section, let us incorporate the high-order tensor methods [22] into Algorithm 1 for solving the corresponding inner subproblem (3.5). From now on, we restrict our attention to Euclidean norms. Let us fix symmetric positive-definite linear operator B : E → E * (notation B = B * 0) and use the following norm for the primal space: x ≡ Bx, x 1/2 , x ∈ E. The norms for multilinear forms on E are induced in the standard way (see section 2).
For this setup, we use the following simple prox function: Thus, the choice of prox function (4.2) is strictly related to the preferable degree p ≥ 1 of smoothness of function f . Let us define the Taylor approximation Ω p (f, x; y) of function f around the point x ∈ dom f : By Assumption 4.1, we are able to bound its accuracy in the following way: for all x, y ∈ dom f it holds that Let us look at our regularized objective h k+1 (·), which needs to be minimized at every step k ≥ 0: . This is a sum of two convex functions: smooth component g k+1 and possibly nonsmooth but simple component φ k+1 , which is strongly convex with respect to d.
Let us drop unnecessary indices and consider the subproblem in a general form: with g having bounded Lipschitz constant for some p ≥ 1: 0 < L p (g) < +∞. Since we assume the objective to be strongly convex with respect to d from (4.2) with parameter σ d (h) > 0, for every x, y ∈ dom h and all h (x) ∈ ∂h(x) we have Bound (4.3) motivates us to define the following point, and consider the following iteration process, For p = 1, the point (4.8) is used in the composite gradient method [20]. For p = 2, this is a step of composite cubic Newton [6,9]. It can be shown that for M ≥ pL p (g) the auxiliary optimization problem in (4.8) is convex for all p ≥ 1 (see Theorem 1 in [22]). Therefore it can be efficiently solved by different techniques of convex optimization and linear algebra (see also [23,22]).
Let us mention some properties of point T ≡ T M (h; x). Its characteristic condition is as follows: Therefore,

This inclusion justifies notation
In order to work with these objects, we use the following result (see Lemma 2 in [4]).
In particular, if β = p, then The next lemma describes the global behavior of method (4.9).
Lemma 4.3. Let β ≥ 1 and M = βL p (g). Then for any x, y ∈ dom h we have Now, we are ready to prove a convergence result on the iteration process (4.9).
Theorem 4.4 (convergence of tensor method). Let M = pL p (g). Then, for every t ≥ 0 and y ∈ dom h we have Proof. Let us consider the point z t+1 = T M (z t ). By (4.12), we have for any y ∈ dom h.
Denote x * h def = argmin y∈E h(y) and consider y = z t +α(x * h −z t ) for α ∈ [0, 1]. Then we have The minimum of the right-hand side is attained at Plugging it into (4.15) gives Therefore, for every t ≥ 0 we have for every y ∈ dom h. It remains to use Lemma 4.2 and finish the proof: Thus, we can see that, applying tensor method (4.9) of degree p ≥ 1 on step 4 of the general contracting proximal method (Algorithm 1), we obtain fast linear convergence for the norms of subgradients. Hence, we can estimate the total number of inner steps t k at iteration k ≥ 0 as follows.
using M := pL p (g k+1 ) and z 0 := v k . Then we have Proof. By definition, for all x ∈ dom ψ, we have Therefore, Then for y ≡ x * def = argmin y∈E F (y) we obtain It remains to use this bound together with (4.13) and the following estimation of the strong convexity parameter: By representation (4.5), we have a simple relations between Lipschitz constants of the derivatives for function g k+1 (·) and f (·): Therefore, we can control the condition number of our objective. Indeed, by (4.17), the main complexity factor in the minimization process for h k+1 (·) is the ratio We are able to keep this ratio small by applying an appropriate growth strategy for coefficients A k . Let us consider two cases: σ d (ψ) = 0 and σ d (ψ) > 0.
and we get Thus we obtain 2. σ d (ψ) > 0. For k = 0 we pick a 1 ≡ c(p + 1) as in the previous case. Now consider k ≥ 1. Denote and choose a k+1 from the equation Thus, in both cases, at every upper-level step we need to perform a logarithmic number of iterations of the inner method, multiplied by a small constant.
We are ready to specify the whole optimization procedure.
Let us present global complexity bounds for this method in convex and strongly convex cases.
Theorem 4.6 (convex case). Let for a given ε > 0, the inner accuracy δ be fixed as follows: Then, in order to achieve F (x K ) − F * ≤ ε it is enough to perform iterations of Algorithm 2. The total number of oracle calls N K def = K k=1 t k is bounded as (4.27) N K ≤ K · 3 + p+1 p log 4 1 + 1 γ0 (p + 1) Proof. Estimate (4.26) follows directly from (3.13), by substituting the value c = p!γ0 . Now, let us prove (4.27). By (4.17), we have In order to finish the proof, we need to bound the value under the logarithm. By the choice of a k , we have an upper bound for A k : Therefore, for every 0 ≤ k ≤ K: This completes the proof. Now, let us discuss the overall dependence of δ and K on p, given by the claim of Theorem 4.6. For simplicity, we fix x * ) = 1, and γ 0 = 1. Thus, we observe the functions One can see that log 2 δ(p) ≤ −p. Therefore, increasing the order of the method, it requires at least to double the precision of solving the subproblem. At the same time, we have (using Stirling's formula) lim p→+∞ K(p) = 1 + 2 exp lim p→+∞ (p+2) log(p+1)−log p! p+1 = 1 + 2 exp(1).
Hence, the value of K(p) is bounded from above by an absolute constant. The graphs of the dependence (4.29) are shown in Figure 1. Note that in practice, we are interested rather in small values of p.   .24). Let for a given ε > 0 the inner accuracy δ be fixed as follows: Then, in order to achieve F (x K ) − F * ≤ ε, it is enough to perform (4.31) The total number of oracle calls N K is bounded as follows: (4.32) Therefore it holds that λ 1 = 1/(1 + exp(α)) and λ n = 1/(1 + exp(−α)), so parameter α is related to the condition number of the problem.
In Table 1 we demonstrate the number of iterations and the total number of matrix-vector products, which are required for the methods to solve the problem up to ε = 10 −7 accuracy in functional residual.  Table 1: Minimization of quadratic function, q = λ min (A)/λ max (A).
We see that the contracting proximal method is always better than the usual proximal algorithm. It requires about the same number of iterations as the accelerated gradient methods, but it needs to spend more oracle calls per iteration, which confirms the theory.

Log-Sum-Exp.
In the next example we compare a performance of secondorder methods for unconstrained minimization of the following objective: where µ > 0 is a parameter, while coefficients of the vectors {a i } m i=1 and b are randomly generated, and we set m = 6n.
We compare the cubically regularized Newton method [23] and its accelerated variant from [19] with the contracting proximal cubic Newton (Algorithm 2 with p = 2) for minimizing the objective up to ε = 10 −8 accuracy in functional residual. In these algorithms we use the following Euclidean norm for the primal space, x = Bx, x 1/2 , with matrix B = m i=1 a i a T i , and fix the regularization parameter equal to 1. The results are shown in Table 2  We see that the contracting proximal method outperforms the direct methods in the number of iterations, but usually requires additional oracle calls for solving the subproblem.
6. Conclusion. In this work, we propose a general acceleration scheme, based on proximal iterations. There are two distinguishing features of our methods: employing the contraction of the smooth component of the objective (this provides the acceleration) and flexibility of prox-function (its choice should take into account both the geometry of the problem and the order of the smoothness).
One of the recent important applications of the accelerated proximal point methods in machine learning is the universal framework Catalyst, applicable to the firstorder methods [13,14]. This is a powerful approach for accelerating many specific optimization methods in a common way. We believe that the results of this paper can help in advancing in this direction, resulting in the faster high-order methods for many practical applications.