Exploiting higher-order derivatives in convex optimization methods

It is well known since the works of I. Newton [64] and L. Kantorovich [45] that the second-order derivative of the objective function can be used in numerical algorithms for solving optimization problems and nonlinear equations and that such algorithms have better convergence guarantees. Higher-order derivatives can be efficiently used for solving nonlinear equations as was shown by P. Chebyshev [17, 24]. For optimization problems, the basic idea known at least since 1970’s [39] and developed in later works [69, 68, 9] is to approximate, at each iteration of an algorithm, the objective by its Taylor polynomial at the current iterate, optionally add a regularization, and minimize this Taylor approximation to obtain the next iterate. In this way, the first Taylor polynomial leads to first-order methods that are very well understood, see, e.g., [54, 58], with the optimal methods existing since 1980’s [54, 56].The second-order methods are using the second Taylor polynomial. The most famous representative is the Newton’s method that minimizes at each iteration the second-order quadratic approximation of the objective function. If the second derivative is Lipschitz continuous, the objective is strongly convex, and the starting point is sufficiently close to the solution, then this algorithm has very fast quadratic convergence and requires $\log\log\varepsilon^{-1}$ iterations to reach an $\varepsilon$-solution in terms of the objective value [45]. This bound is optimal [54] even for univariate optimization problems with the possibility of using in algorithms derivatives of any order. Different modifications of the basic algorithm, such as the Damped Newton’s method or the Levenberg–Marquardt algorithm achieve global convergence, but have a slower, i.e., linear, convergence [67, 65]. Second-order methods played also the central role in the development by Yu. Nesterov and A. Nemirovskii of interior point methods that have global linear convergence rate and allow proving polinomial solvability of a large class of convex optimization problems [62]. This theory was based on the analysis of the Damped Newton’s method for the class of self-concordant functions that, in particular, include functions without Lipschitz derivatives.

An important idea that eventually led to the current developments of tensor methods was the cubic regularization of the Newton’s method, which dates back to the work of A. Griewank [37]. The global performance of the Cubic regularized Newton’s method was analysed in 2006 by Yu. Nesterov and B. Polyak in [63] for a large list of settings with the main assumption that the second derivative is Lipschitz continuous. The main idea is based on the fact that, for functions with Lipshitz second derivative, the objective’s model consisting of the second Taylor polynomial regularized by the cube of the norm with sufficiently large regularization parameter is an upper bound for the objective function. This model is minimized in each iteration of the algorithm, which, in particular, allowed to obtain global convergence rate $1/k^{2}$ for minimizing convex functions with Lipschitz second derivative. Moreover, the authors showed that the complexity of minimizing the third-order polynomial model of the objective is of the same order as the complexity of the standard Newton’s method step. The Cubic regularized Newton’s method was further accelerated in [55] to reach $1/k^{3}$ convergence rate for convex problems with Lipschitz second derivative. In 2012 R. Monteiro and B. Svaiter proposed in [53] a very perspective accelerated-proximal envelope (see also the work [16]) that allowed them to develop even faster second-order method with the convergence rate $1/k^{3.5}$ for minimizing convex objectives with Lipschitz second derivative.

Another important step in the development of tensor methods was made by M. Baes in 2009 [8], where the Newton–Nesterov–Polyak algorithm was generalized to the setting of convex minimization with the objective having Lipschitz $p$-th derivative for $p\geq 2$. The idea was to construct a $(p+1)$-th order polynomial model that upper bounds the objective function by taking the $p$-th Taylor polynomial of the objective and regularizing it by the $(p+1)$-th power of the Euclidean norm with sufficiently large regularization parameter. The author showed $1/k^{p}$ global convergence rate for methods that minimize such model in each iteration and proposed an accelerated version with the rate $1/k^{p+1}$, all under the assumption of Lipschitz $p$-th derivative. As in the world of first-order methods, where the interest in optimal methods was one of the central driving forces in 2000–2020, a natural question arose on what are the lower bounds for second- and higher-order methods and which algorithms are optimal. The results in [7, 5, 29] gave the answer that the lower bound is $1/k^{(3p+1)/2}$, revealing the gap when $p\geq 3$ between the rates of existing methods an lower bounds. One of the drawbacks of tensor methods at this stage was that each iteration required to minimize a higher-order polynomial that may not be convex, leading to the high cost of each iteration and impracticality of such methods.

In 2018, a breakthrough was made by Yu. Nesterov [59] in understanding the significance of higher-order methods in modern convex optimization when the $p$-th derivative satisfies the Lipschitz condition. The breakthrough involved increasing the regularization parameter in M. Baes’s approach for the regularized $p$-th Taylor polynomial of the objective. It was demonstrated that this modified Taylor model is convex. Additionally, it was shown that the convex Taylor model can be efficiently minimized for $p=3$, with similar complexity to the standard Newton’s method step. Furthermore, it was proven that minimizing the convex Taylor model doesn’t require the calculation of the entire tensor of third derivatives. Instead, directional third derivatives can be calculated, for example, through automatic differentiation. Lastly, a lower bound of $1/k^{(3p+1)/2}$ was established for convex problems with Lipschitz $p$-th derivative. For each $p\geq 2$, worst-case functions within the class of convex functions with Lipschitz $p$-th derivatives were constructed, and a research program was outlined to develop optimal tensor methods [58]. It followed from the obtained lower bounds for convex optimization problems that the Monteiro–Svaiter second-order method is optimal up to a logarithmic factor since it used a complicated line-search procedure. Based on the Monteiro–Svaiter method and the convex Taylor model of the objective, in [30, 26, 27] and independently in [11, 42], the authors proposed near-optimal tensor methods for convex problems with Lipschitz $p$-th derivatives, which have the rate $1/k^{(3p+1)/2}$ up to logarithmic factors. In Section 3, the Monteiro–Svaiter approach with Nesterov’s implementable tensor steps is described. As it was mentioned in [58, 21], the auxiliary line-search procedures that lead to logarithmic factors in the convergence rates of the near-optimal tensor methods may also significantly slow down the convergence in practice. In 2022, two independent works [48] and [14] proposed optimal tensor methods with convergence rates without additional logarithmic factors. At the same time, such methods were proposed also for monotone variational inequalities [1, 51] under higher-order regularity of the operator.

The developments in the theory of tensor methods had also an important impact on the theory of second-order methods. In 2020, Yu. Nesterov proposed an implementation of the third-order method using inexact third derivative constructed via finite differences of gradients [61]. In other words, it appeared to be possible to implement the third-order method using only first and second derivatives, which lead to <<superfast>> second-order method with the convergence rate $1/k^{4}$ violating the lower bound thanks to additional assumption that the third derivative is Lipschitz continuous. The key observation was that in the class of functions with Lipschitz second derivative the worst-case function [59] for second-order methods does not have Lipschitz third derivative. These results were further improved in [43, 60], where second-order algorithms with the rate $1/k^{5}$ corresponding to optimal third-order methods were proposed for minimizing convex functions with Lipschitz third derivative. These developments are in large contrast to first-order methods, for which the worst-case function is quadratic [57] and has Lipschitz derivatives of any order, thus, preventing improvements under additional assumptions. This line of research was continued in [6], where such reductions were shown to be possible for higher-order methods. The ideas used in superfast methods are described in Section 4.

Tensor methods remain a very active area of research with many extensions of the existing methods. In particular, there are adaptive variants for the setting when the Lipschitz constant is not known [41, 36], universal generalizations for the setting when the $p$-th derivative is Hölder continuous [33, 70, 19], versions with inexact minimization of the Taylor model [21, 35], tensor methods for finding approximate stationary points of convex functions [34, 22]. The ideas described above, especially the Monteiro–Svaiter accelerated proximal point method, turned out to be productive also in the areas not directly related to tensor methods, see non-trivial examples in [10, 12, 15]. Modern second-order and third-order methods demonstrate also their efficiency in Data Science and Machine Learning applications [18, 2, 23, 13]. More details of such applications are given in Section 5.

The following optimization problem is considered

where $f$ and $g$ are convex functions. Denote $x^{*}$ – the solution of the problem (1). If the solution is not unique, let $x^{*}$ be a solution that is the closest to starting point $x_{0}$ in $2$-norm.

Denote $\|\,\cdot\,\|$ the Euclidean $2$-norm in $\mathbb{R}^{d}$,

Assume that $f$ has Lipschitz derivatives of order $p$ ($p\geq 1$):

Here and below (see e.g. (5)) it is considered that $x,y\in\mathbb{R}^{d}$ belongs to the ball in $2$-norm centred at $x^{*}$ with the radius $O(\|x_{0}-x^{*}\|)$ [60].

The $p$-th Taylor polynomial of $f$ is defined as

Note that from (2) it follows [59] that

$F$ satisfies $r$-growth condition ($p+1\geq r\geq 2$) (see [66]) with constant $\sigma_{r}>0$ iff

The following algorithm is taken from [28], for $g\equiv 0$ see [11].

Note that Theorem 1 is true also when $H\geq 2L_{p,f}$ (independently of $p\geq 1$). It can be derived from (4). The condition $H\geq(p+1)L_{p,f}$ was used only because it guarantees the convexity of auxiliary problem (6), see [59]. Under the conditions $H\geq(p+1)L_{p,f}$, $g\equiv 0$ for $p=1,2,3$ there exist efficient ways to solve auxiliary problem (6), see [59]. For $p=1$ there exists an explicit formula for the solution of (6). For $p=2,3$ the complexity of (6) is almost the same (up to a logarithmic factor in $\varepsilon$) as a complexity of Newton’s method iteration [59], see also Section 4. It is important that there is no need to solve (6) accurately. It is sufficient to find such $\tilde{y}_{k+1}$ that satisfied

where

Such a practical modification slow down the theoretical convergence rate in a factor $12/5$ in the right hand side of (7) [43, 35, 60].

In May 2022, D. Kovalev et al. [48] proposed for $g\equiv 0$ an explicit policy for choosing a pair $(\lambda_{k+1},y_{k+1})$ in Algorithm 1 and modify the stopping criteria (8) according to [40]. It allows to solve (6) on each iteration at average only two times rather $O\left(\ln\left(\varepsilon^{-1}\right)\right)$ times as it was before. The final complexity bound for $p$-order oracle calls in [48] matched the lower oracle complexity bound from [59] (see also [25]) obtained for the worst-case function

In concurrent and independent paper [14], the authors also proposed a way of reducing additional logarithmic factors. Note, that approach of [14] does not require any a priory knowledge of smoothness parameters (including Holder continuity assumption instead of Lipschitz one).

If additionally $F$ satisfies $r$-growth condition (5) then optimal method can be developed based on restarts procedure [63, 27] – see Algorithm 2.

Everything that was noted after Theorem 1 will also take place in this case.

In particular, when $g\equiv 0$, $p\geq 2$ and $r=2$ MSN algorithm in Theorem 2 can be replaced by Kovalev’s variant of MSN [48] to obtain

$p$-order oracle complexity bound. This upper bound corresponds to the $\sigma_{2}$-strongly convex case lower bound from [46] and improves (10) on a logarithmic factor. The dependence on $\varepsilon$ is $\log\log\varepsilon^{-1}$ as it should be locally for tensor methods [26, 20] (see also [54, 63, 5, 19]). But (11) describe two regimes: the first term (complexity independent on $\varepsilon$) describe the complexity to reach the vicinity of quadratic convergence, the second term (complexity is $\log\log\varepsilon^{-1}$) describe Newton’s type local convergence.

Note that due to the presence of composite term $g$ the described above MSN algorithm and its variations can be used for splitting oracle complexities [44, 47]. Namely, assume that $f$ and $g$ have Lipshitz $p$-th order derivatives with different Lipschitz constants. Based on MSN algorithm one can propose a general framework to accelerate tensor methods by splitting computational complexities. As a result, one can get near-optimal oracle complexity for each function in the sum separately for any $p\geq 1$, including the first-order methods. To be more precise, if the near optimal complexity to minimize $f$ is $N_{f}(\varepsilon)$ and to minimize $g$ is $N_{g}(\varepsilon)$, then MSN-based sliding algorithm [44, 47] requires no more than $\tilde{O}\left(N_{f}(\varepsilon)\right)$ oracle calls for $f$ and $\tilde{O}\left(N_{g}(\varepsilon)\right)$ oracle calls for $g$ to minimze $F=f+g$.

Not also that for $p=1$ and $f\equiv 0$ one can take $H>0$ in MSN and obtain Monteiro–Svaiter accelerated proximal envelop. In the cycle of recent papers [40, 16, 47, 49, 48] it was shown that this type of proximal (Catalyst-type [50]) envelop is logarithmic-free due to the well developed stopping rule criteria for the inner (auxiliary) problem. So it opens up different perspectives for improving existing Catalyst-based algorithms, see e.g. [72, 71].

In this section, the computationally efficient solution of the tensor subproblem for $p=3$ is discussed. With the proposed procedure, it is possible to implement third-order methods without the computation and storing of the third-order derivative. At the beginning of the section, it is shown that this subproblem is convex. Then the Bregman-distance gradient method is used as a subsolver for the tensor step. (6).

For this section, it is assumed that $g=0$ and $p=3$. Then $L_{3}=L_{3,f}$, the third-order Taylor’s polynomial is

The regularized third-order model is

The basic step for every tensor method is formulated as

This step is a major part of almost every third-order method. Next, it is described how to effectively solve this subproblem without the computation of the third-order derivative and only using second-order information.

First, it is clarified that the function (13) is convex for $H\geq 3L_{3}$. Note, that if $H\geq L_{3}$ the model (13) is possibly non-convex because of third-order derivative. It means that the minimizing non-convex subproblem is harder than minimizing original convex problem. In 2018 Yu. Nesterov made the breakthrough in [59]. It was shown that if $H\geq 3L_{3}$ then the model (13) is convex and high-order methods are implementable.

To give some intuition, a sketch of the convexity proof is presented. The following inequality can be derived from the Lipschitz-continuous condition and the convexity of function $f$.

Now, to finish the proof, one need to choose $H$ such that

This lead to the crucial detail, that was misleading before

So, from the fact that $x$ is a vector and $xx^{T}$ is a singular matrix, the factor $3$ in the last inequality is losing. Finally, if $H\geq 3L_{3}$ then the function (13) is convex. For more details one can check Theorem 1 in [59].

Thus, one can move to the effective subsolver of the problem (14) by Bregman-distance gradient method for relatively smooth functions from [52]. The main idea is to show that the optimized function $\phi(y)$ is $L_{\rho}$-smooth and $\mu_{\rho}$-strongly convex with respect to some convex function $\rho(y)$ or

Then Bregman-distance gradient method make next steps

where

is a Bregman-distance generated by $\rho(y)$. Bregman-distance gradient method converges linearly with condition number $\kappa=\tfrac{L_{\rho}}{\mu_{\rho}}$ and convergence rate

One can show that the model (13) with $H=6L_{3}$ can be optimized as $\phi(y)=\tilde{\Omega}(f,x_{k};y)$ by gradient method with Bregman-distance generated by

$L_{\rho}=1+\tfrac{1}{\sqrt{2}}$, $\mu_{\rho}=1-\tfrac{1}{\sqrt{2}}$, and $\kappa=\tfrac{1}{\left(1+\sqrt{2}\right)^{2}}=O(1)$. This means that this method is very fast and converges with a fixed number of iterations. It is worth noting that for each step, the full hessian for $\beta$ is computed, but the full third-order derivative is unnecessary because one only need the derivative-vector-vector product $D^{3}f(x_{k})\left[y-x\right]^{2}$ to compute $\nabla\tilde{\Omega}(f,x_{k};y)$. Derivative-vector-vector product can be efficiently and precisely computed by autogradient or approximated by finite difference. To summarize, the tensor method with a convergence rate of $1/k^{5}$ is obtained for third-order methods but only the first and second-order derivatives are computed. This approach was firstly proposed in [59] and then it was improved in [61]. Section 5 from [59] is good for a general understanding of this approach. For details and precise proofs, one can read [61].

As well as first-order methods, modifications of tensor methods are proposed for problems of the finite sum type, stochastic and distributed optimization.

In each iteration higher-order (also called tensor) methods minimize a regularized Taylor expansion of the objective function, which leads to faster convergence rates if the corresponding higher-order derivative is Lipschitz-continuous. Recently a series of lower iteration complexity bounds for such methods were proved, and a gap between upper an lower complexity bounds was revealed. Moreover, it was shown that such methods can be implementable since the appropriately regularized Taylor expansion of a convex function is also convex and, thus, can be minimized in polynomial time. Only very recently an algorithm with optimal convergence rate $1/k^{((3p+1)/2)}$ was proposed for minimizing convex functions with Lipschitz $p$-th derivative. For convex functions with Lipschitz third derivative, these developments allowed to propose a second-order method with convergence rate $1/k^{5}$, which is faster than the rate $1/k^{3.5}$ of existing second-order methods.

90C25,65K10,90C15

• •

  Automatic Differentiation: Calculation of Newton Steps; Dixon, L.

• •

  Computational Optimal Transport ; Dvurechensky P., Dvinskikh D., Tupitsa N., Gasnikov A.

• •

  Stochastic gradient descent; Lan G.

• •

  Unified analysis of SGD-type methods; Gorbunov E.

• •

  Decentralized Convex Optimization over Time-Varying Graphs; Rogozin A., Gasnikov A., Beznosikov A., Kovalev D.

• •

  Stochastic Quasi-Newton Scheme; Jalilzadeh, A.