OPTAMI: Global Superlinear Convergence of High-order Methods

## Abstract

Abstract Second-order methods for convex optimization outperform first-order methods in terms of theoretical iteration convergence, achieving rates up to O ⁢ ( k − 5 ) 𝑂 superscript 𝑘 5 O(k^{-5}) italic_O ( italic_k start_POSTSUPERSCRIPT - 5 end_POSTSUPERSCRIPT ) for highly-smooth functions. However, their practical performance and applications are limited due to their multi-level structure and implementation complexity. In this paper, we present new results on high-order optimization methods, supported by their practical performance. First, we show that the basic high-order methods, such as the Cubic Regularized Newton Method, exhibit global superlinear convergence for μ 𝜇 \mu italic_μ -strongly star-convex functions, a class that includes μ 𝜇 \mu italic_μ -strongly convex functions and some non-convex functions. Theoretical convergence results are both inspired and supported by the practical performance of these methods. Secondly, we propose a practical version of the Nesterov Accelerated Tensor method, called NATA. It significantly outperforms the classical variant and other high-order acceleration techniques in practice. The convergence of NATA is also supported by theoretical results. Finally, we introduce an open-source computational library for high-order methods, called OPTAMI. This library includes various methods, acceleration techniques, and subproblem solvers, all implemented as PyTorch optimizers, thereby facilitating the practical application of high-order methods to a wide range of optimization problems. We hope this library will simplify research and practical comparison of methods beyond first-order.

## 1 Introduction

In this paper, we consider the following unconstrained optimization problem:

$$ $\min\limits_{x\in\mathbb{E}}f(x),\vspace{-.1cm}$ (1) $$

where $\mathbb{E}$ is a $d$-dimensional real value space and $f(x)$ is a highly-smooth function:

###### Definition 1.1

Function $f$ has $L_{p}$ - Lipschitz-continuous $p$-th derivative, if

$$ $\|D^{p}f(x)-D^{p}f(y)\|_{op}\leq L_{p}\|x-y\|\qquad\forall x,y\in\mathbb{E},$ (2) $$

where $D^{p}f(x)$ is a $p$-th order derivative, and $\|\cdot\|_{op}$ is an operator norm.

In the paper, we primarily focus on three main cases: $p=1$, $p=2$, and $p=3$. We assume that function $f(x)$ has at least one minimum $x^{\ast}$. For convexity-type assumptions, we use star-convexity and uniform star-convexity for non-accelerated methods and convexity for accelerated methods.
Note, that the classes of star-convex or uniformly star-convex functions include some non-convex functions.

###### Definition 1.2

Let $x^{\ast}$ be a minimizer of the function $f$. For $q\geq 2$ and $\mu_{q}\geq 0$, the function $f$ is $\mu_{q}$-uniformly star-convex of degree $q$ with respect to $x^{\ast}$ if for all
$x\in\mathbb{R}^{d}$
and $\forall\alpha\in[0,1]$

$$ $f\left(\alpha x+(1-\alpha)x^{\ast}\right)\leq\alpha f(x)+(1-\alpha)f(x^{\ast}) \\ -\tfrac{\alpha(1-\alpha)\mu_{q}}{q}\|x-x^{\ast}\|^{q}.$ (3) $$

If $q=2$ then the function $f$ is $\mu$-strongly star-convex with respect to $x^{\ast}$.
If $\mu_{q}=0$ then the function $f$ is star-convex with respect to $x^{\ast}$. From this definition, we can additionally get the next useful inequality sometimes called $q$-order growth condition

$$ $\tfrac{\mu_{q}}{q}\|x-x^{\ast}\|^{q}\leq f(x)-f(x^{\ast}).$ (4) $$

We introduce two types of distance measures between the starting point and the solution: for non-accelerated methods, we consider the diameter of the level set $\mathcal{L}~{}=~{}\left\{x\in\mathbb{E}:f(x)\leq f(x_{0})\right\}$

$$ $D~{}=~{}\max\limits_{x\in\mathcal{L}}\|x-x^{\ast}\|;$ (5) $$

and for accelerated methods, we use the Euclidean distance given by

$$ $R=\|x_{0}-x^{\ast}\|.$ (6) $$

Second-order methods are widely used in optimization, finding applications in diverse fields such as machine learning, statistics, control, and economics (Polyak, ; Boyd and Vandenberghe, ; Nocedal and Wright, ; Nesterov, ). One of the main and well-known advantages of the second-order methods is their local quadratic convergence rate.
The classical (Damped) Newton method (Newton, ; Raphson, ; Kantorovich, ) can be written as:

$$ $x_{t+1}=x_{t}-\gamma_{t}\left[\nabla^{2}f(x_{t})\right]^{-1}\nabla f(x_{t}), \vspace{-.1cm}$ (7) $$

where $\gamma_{t}\in\mathbb{R}_{+}$ is a step-size or damping coefficient. The Newton step originates from the second-order Taylor expansion $\Phi_{2}(x,x_{t})$:

$$ $x_{t+1}=\operatornamewithlimits{argmin}_{x\in\mathbb{E}}\left\{\Phi_{2}(x,x_{t })=f(x_{t})+\left\langle\nabla f(x_{t}),x-x_{t}\right\rangle+\left\langle \nabla^{2}f(x_{t})(x-x_{t}),x-x_{t}\right\rangle\right\}.\vspace{-.2cm}$ (8) $$

The solution of this problem corresponds to (7) with $\gamma_{t}=1$.
Despite its fast local convergence, the Damped Newton method has a slow global convergence rate $O(T^{-1/3})$ (Berahas et al., ), which is even slower than the gradient method’s convergence rate $O(T^{-1})$. One reason for this slow convergence is that the approximation $\Phi_{2}(x,x_{t})$ is not an upper-bound for $f(x)$, and it is possible that $f(x)>\Phi_{2}(x,x_{t})$. Consequently, the step could be too aggressive and without a damping step-size $\gamma_{t}$, the Newton method could diverge (Nesterov, , Example 1.2.3).
To address this issue, the Cubic Regularized Newton (CRN) method was proposed by Nesterov and Polyak ().

$$ $x_{t+1}=\operatornamewithlimits{argmin}_{y\in\mathbb{E}}\left\{\Omega_{M_{2}}( x,x_{t})=\Phi_{2}(x,x_{t})+\tfrac{M_{2}}{6}\|x-x_{t}\|^{3}\right\}.\vspace{-.2cm}$ (9) $$

For the function $f(x)$ with $L_{2}$-Lipschitz Hessian, the model $\Omega_{M_{2}}(y,x_{t})$ is an upperbound of the function $f(x)$ for $M_{2}\geq L_{2}$; hence $\Omega_{M_{2}}(x,x_{t})\geq f(x)$.
This method is the first second-order method with a global convergence rate of $O(T^{-2})$, which is faster than the Gradient Method (GM). It also has a quadratic local convergence rate for strongly convex functions similar to the classical Newton method (Doikov and Nesterov, ).

In large-scale optimization problems, the computation of the inverse Hessian can be rather costly.
This limitation significantly restricts the practical application of second-order methods for large-scale problems.
One of the possible solutions to this problem is to approximate the Hessian using gradient information or Hessian-vector products (HVP). These methods are known as Quasi-Newton (QN) methods, which include BFGS (Broyden, ; Fletcher, ; Goldfarb, ; Shanno, ), L-BFGS (Liu and Nocedal, ), SR-1 (Conn et al., ; Khalfan et al., ; Byrd et al., ), among others (Nocedal and Wright, ; Gupta et al., ). The approximate Hessian can also have a form of diagonal preconditioning (Yao et al., ; Jahani et al., ; Sadiev et al., ). The classical form of the QN update looks as follows:

$$ $x_{t+1}=x_{t}-\gamma_{t}B_{t}^{-1}\nabla f(x_{t})=x_{t}-\gamma_{t}H_{t}\nabla f (x_{t}),\vspace{-.1cm}$ (10) $$

where $B_{t}\approx\nabla^{2}f(x_{t})$ is a Hessian approximation and $H_{t}\approx[\nabla^{2}f(x_{t})]^{-1}$ is an approximation of inverse Hessian, in a lot of cases $H_{t}=B_{t}^{-1}$. These methods have local superlinear convergence and could perform well globally in practice. However, they have poor global convergence $O(T^{-1/3})$ (Berahas et al., ), similar to the Newton method. In (Kamzolov et al., ), the Cubic Regularized Quasi-Newton method was proposed with a global convergence rate $O\left(\tfrac{\delta D^{2}}{T}+\tfrac{L_{2}R^{3}}{T^{2}}\right)$, where $\delta=\max_{x_{t}}\|B_{t}-\nabla^{2}f(x_{t})\|_{2}$
is a measure of inexactness of the Hessian approximation. Depending on the approximation, one can theoretically guarantee that $\delta\leq L_{1}$, which makes the method intermediate between the gradient method and the Cubic Regularized Newton method. This observation indicates that second-order methods can surpass first-order methods; even when relying solely on first-order information, the Cubic Quasi-Newton method can perform at least as well as, if not better than, the gradient method.

Such potential motivates us to study second-order methods. The Cubic Regularized Newton is the basic method in the line-up of second-order methods. There are two main directions for its improvement: accelerated second-order methods, including Nesterov-type acceleration (Nesterov, ; ), A-NPE near-optimal acceleration (Monteiro and Svaiter, ; Gasnikov et al., ), and optimal acceleration (Kovalev and Gasnikov, ; Carmon et al., ); and third-order methods with superfast subsolvers, which allows making a third-order step without computation of third-order derivative (Nesterov, ; ; ; Kamzolov, ).
In this paper, we explore both of these directions from theoretical and practical perspectives, which may provide valuable insights into the practical performance of various acceleration techniques and high-order methods.

We are particularly interested in $\mu$-strongly convex function $f(x)$ with $L_{2}$-Lipschitz Hessian. Existing results on second-order methods primarily focus on the local quadratic convergence. Specifically, these method are capable of achieving an $\varepsilon$-solution $x_{T}$, such that $f(x_{T})-f^{\ast}\leq\varepsilon$, within at most $\tilde{O}(\log\log\tfrac{1}{\varepsilon})$ iterations, where $\tilde{O}(\cdot)$ denotes complexity up to a logarithmic factors. To get this rate, the methods have to be inside of the quadratic convergence region $\mathcal{Q}=\left\{x\in\mathbb{E}:f(x)-f^{\ast}\leq r\right\}$, where the radius $r$ is given by

$$ $r=\tfrac{c_{1}\mu^{3}}{L_{2}^{2}}\vspace{-.1cm}$ (11) $$

for some numerical constant $c_{1}$. Now, let us illustrate this idea using the example of Cubic Regularized Newton method (CRN) (Doikov and Nesterov, ). Initially, the methods performs at most

$$ $T=O\left(\left(\tfrac{L_{2}D}{\mu}\right)^{1/2}\log\left(\tfrac{f(x_{0})-f^{ \ast}}{\varepsilon}\right)\right)$ (12) $$

iterations to reach $f(x_{T})-f^{\ast}\leq\varepsilon$ for $\varepsilon>r$.
If $\varepsilon\leq r$, the CRN method enters the quadratic convergence region after at most $T_{1}=O\left(\left(\tfrac{L_{2}D}{\mu}\right)^{1/2}\log\left(\tfrac{f(x_{0})-f
^{\ast}}{r}\right)\right)$ iterations. Once in this region, the method requires at most $T_{2}=O\left(\log\log\left(\tfrac{r}{\varepsilon}\right)\right)$ additional iterations. Therefore, the total number of CRN iterations to reach $f(x_{T})-f^{\ast}\leq\varepsilon$ for $\varepsilon\leq r$ is upper bounded by

$$ $T=O\left(\left(\tfrac{L_{2}D}{\mu}\right)^{1/2}\log\left(\tfrac{(f(x_{0})-f^{ \ast})L_{2}^{2}}{\mu^{3}}\right)+\log\log\left(\tfrac{\mu^{3}}{L_{2}^{2} \varepsilon}\right)\right).$ $$

For the class of $\mu$-strongly convex functions with $L_{2}$-Lipschitz Hessian, the lower bound is

$$ $\Omega\left(\left(\tfrac{L_{2}D}{\mu}\right)^{2/7}+\log\log\left(\tfrac{\mu^{3 }}{L_{2}^{2}\varepsilon}\right)\right)$ $$

for $\varepsilon\leq c_{3}\tfrac{\mu^{3}}{L_{2}^{2}}=c_{3}r$ as established by Arjevani et al. (). The power $2/7$ corresponds to the optimal accelerated method. However, this lower bound applies only when $\varepsilon\leq c_{3}r$, which corresponds to small values of $\varepsilon$. In the case when $\varepsilon\leq c_{3}r$, meaning the desired accuracy exceeds the radius of the quadratic convergence region, it may be possible to achieve faster global rates of Cubic Regularized Newton method than the rate given by (12).
Another problem, we are interested in, is the practical performance of Nesterov Accelerated Tensor Method (Nesterov, ). In contrast to first-order methods, classical acceleration is slower or similar to the non-accelerated version in practice, which leads to avoidance of the method (Scieur, ; Carmon et al., ; Antonakopoulos et al., ).
Lastly, there is a wide variety of acceleration techniques, which include Nesterov acceleration (Nesterov, ) with a rate $O(T^{-(p+1)})$; Near-Optimal Monteiro-Svaiter Acceleration (Monteiro and Svaiter, ; Bubeck et al., ; Gasnikov et al., ; Kamzolov, ) with a rate $\tilde{O}(T^{-(3p+1)/2})$; Near-Optimal Proximal-Point Acceleration (Nesterov, ) with the rate $\tilde{O}(T^{-(3p+1)/2})$;
Optimal Acceleration (Kovalev and Gasnikov, ; Carmon et al., ) with a rate $O(T^{-(3p+1)/2})$ and more (Nesterov, ). However, the current literature is lacking their practical comparison, especially for high-order methods with $p=3$. All these considerations lead us to the next questions about the performance of high-order methods:

Can second-order methods for strongly convex functions achieve global rates faster than linear?
Why may Nesterov Accelerated Tensor method underperform in practice, how to overcome this issue?
Which acceleration technique for high-order methods is faster in practice?

In our work, we attempt to answer these questions in a systematic manner.

Contributions. We summarize our key contributions as follows:

- 1.
Global Superlinear Convergence of Second and High-order methods. Our main contribution is providing theoretical guarantees for global superlinear convergence of the Cubic Regularized Newton Method and the Basic Tensor Methods for $\mu$-strongly star-convex functions. These theoretical results are validated by practical performance. These results are a significant improvement over the current state-of-the-art in second-order methods.
- 2.
Nesterov Accelerated Tensor Method with $A_{t}$-Adaptation (NATA).
We propose a new practical variant of the Nesterov Accelerated Tensor Method, called NATA. This method addresses the practical limitations of the classical version of acceleration for high-order methods. We demonstrate the superior performance of NATA compared to both the classical Nesterov Accelerated Tensor Method and Basic Tensor Method for $p=2$ and $p=3$. We also prove a convergence theorem for NATA that matches the classical convergence rates.
- 3.
Comparative Analysis of High-Order Acceleration Methods.
We provide a practical comparison of state-of-the-art (SOTA) acceleration techniques for high-order methods, with a focus on the cases $p=2$ and $p=3$. Our experiments show that the proposed NATA method consistently outperforms all SOTA acceleration techniques, including both optimal and near-optimal methods.
- 4.
Open-Source Computational Library for Optimization Methods (OPTAMI(^1^11[https://github.com/OPTAMI/OPTAMI](https://github.com/OPTAMI/OPTAMI))).
We introduce OPTAMI, an open-source library for high-order optimization methods. It facilitates both practical research and applications in this field. Its modular architecture supports various combinations of acceleration techniques with basic methods and their subsolvers. All methods are implemented as PyTorch optimizers. This allows for seamless application of high-order methods to a wide range of optimization problems, including neural networks.

###### Definition 1.1

###### Definition 1.2

## 2 Basic Methods

Notation.
In the paper, we consider a $d$-dimensional real value space $\mathbb{E}$. $\mathbb{E}^{\ast}$ is a dual space, composed of all linear functionals on $\mathbb{E}$. For a functional $g\in\mathbb{E}^{\ast}$, we denote by $\left\langle g,x\right\rangle$ its value at $x\in\mathbb{E}$. For $p\geq 1$, we define $D^{p}f(x)[h_{1},\ldots,h_{p}]$ as a directional $p$-th order derivative of $f(x)$ along $h_{i}\in\mathbb{E},\,i=1,\ldots,p$. If all $h_{i}=h$ we simplify $D^{p}f(x)[h_{1},\ldots,h_{p}]$ as $D^{p}f(x)[h]^{p}$.
So, for example, $D^{1}f(x)[h]=\langle\nabla f(x),h\rangle$ and
$D^{2}f(x)[h]^{2}=\left\langle\nabla^{2}f(x)h,h\right\rangle$. Note, that
$\nabla f(x)\in\mathbb{E}^{\ast}$, $\nabla^{2}f(x)h\in\mathbb{E}^{\ast}$. Now, we introduce different norms for spaces $\mathbb{E}$ and $\mathbb{E}^{\ast}$. For a self-adjoint positive-definite operator $B:\mathbb{E}\rightarrow\mathbb{E}^{\ast}$, we can endow these spaces with conjugate Euclidian norms:

$$ $\|x\|=\left\langle Bx,x\right\rangle^{1/2},\quad x\in\mathbb{E},\qquad\|g\|_{ \ast}=\left\langle g,B^{-1}g\right\rangle^{1/2},\quad g\in\mathbb{E}^{\ast}. \vspace{-.1cm}$ $$

So, for an identity matrix $B=I$, we get the classical $2$-norm $\|x\|_{2}=\|x\|_{I}=\left\langle x,x\right\rangle^{1/2}$.

High-order model.
High-order Taylor approximation of a function $f(x)$ can be written as follows:

$$ $\Phi_{x,p}(y)=f(x)+\sum_{k=1}^{p}\tfrac{1}{k!}D^{k}f(x)\left[y-x\right]^{k}, \quad x,y\in\mathbb{E},\vspace{-.1cm}$ (13) $$

where, for $p=1$, we simplify notation to $\Phi_{x}(y)$.
From (2), we can get the next upper-bound of the function $f(x)$ (Nesterov, ; )

$$ $|f(y)-\Phi_{x,p}(y)|\leq\tfrac{L_{p}}{(p+1)!}\|y-x\|^{p+1},\vspace{-.1cm}$ (14) $$

which leads us to the high-order model

$$ $\Omega_{x,M_{p}}(y)=\Phi_{x,p}(y)+\tfrac{M_{p}}{(p+1)!}\|y-x\|^{p+1}.\vspace{- .1cm}$ (15) $$

Now, we can formulate the Basic Tensor method

$$ $x_{t+1}=\operatornamewithlimits{argmin}_{y\in\mathbb{E}}\left\{\Omega_{x_{t},M _{p}}(y)\right\},\vspace{-.1cm}$ (16) $$

where $M_{p}\geq pL_{p}$. For $p=1$ and $M_{1}\geq L_{1}$, it is the gradient descent step

$$ $x_{t+1}=x_{t}-\tfrac{1}{M_{1}}\nabla f(x_{t})\vspace{-.1cm}$ (17) $$

with the convergence rate $O\left(\tfrac{M_{1}R^{2}}{T}\right)$ for convex functions.
For $p=2$ and $M_{2}\geq L_{2}$, it is a Cubic Regularized Newton Method (Nesterov and Polyak, ):

$$ $x_{t+1}=x_{t}+\operatornamewithlimits{argmin}_{h\in\mathbb{E}}\left\{f(x_{t})+ \nabla f(x_{t})\left[h\right]+\tfrac{1}{2}\nabla^{2}f(x_{t})\left[h\right]^{2} +\tfrac{M_{2}}{6}\|h\|^{3}\right\},\vspace{-.1cm}$ (18) $$

with the convergence rate $O\left(\tfrac{M_{2}D^{3}}{T^{2}}\right)$.
For $p=3$ and $M_{3}\geq 3L_{3}$, it is a Basic Third-order Method (Nesterov, ):

$$ $x_{t+1}=x_{t}+\operatornamewithlimits{argmin}_{h\in\mathbb{E}}\left\{f(x_{t})+ \nabla f(x_{t})\left[h\right]+\tfrac{1}{2}\nabla^{2}f(x_{t})\left[h\right]^{2} +\tfrac{1}{6}D^{3}f(x_{t})\left[h\right]^{3}+\tfrac{M_{3}}{24}\|h\|^{4}\right \},\vspace{-.1cm}$ (19) $$

with the convergence rate $O\left(\tfrac{M_{3}D^{4}}{T^{3}}\right)$. The step (19) can be performed with almost the same computational complexity (up to a logarithmic factor) by using the Bregman Distance Gradient Method as a subsolver (Nesterov, ; ). The details are written in the Appendix C.1.

## 3 Global Superlinear Convergence of High-order Methods for Strongly Star-Convex Functions

In this section, we establish the global superlinear convergence of high-order methods for uniformly star-convex functions. We begin by defining global superlinear convergence.

###### Definition 3.1

A method is said to exhibit a global superlinear convergence rate with respect to the functional gap if there exists a sequence $\zeta_{t}$ for all $t\in\left\{0,\ldots,T\right\}$ such that

$$ $\tfrac{f(x_{t+1})-f^{\ast}}{f(x_{t})-f^{\ast}}\leq\zeta_{t},\qquad 1>\zeta_{t} >\zeta_{t+1}\quad\forall t\in\left\{0,\ldots,T\right\},\qquad\text{and}\qquad \zeta_{t}\rightarrow 0\,\,\,\textit{for}\,\,\,t\rightarrow 0.\vspace{-.4cm}$ (20) $$

The essence of this definition lies in the fact that the scaling coefficient $\zeta_{t}$ decreases with each iteration. If $\zeta_{t}$ remains constant, the method achieves linear convergence. Conversely, if $\zeta_{t}$ increases over time (i.e., $\zeta_{t}<\zeta_{t+1}$), the convergence becomes sublinear. For example, gradient descent exhibits global linear convergence for strongly convex functions, with $\zeta_{t}=1-\tfrac{\mu}{L_{1}+\mu}$.
Now, to get some intuition on the performance of the methods, we begin with two simple and classical examples: the $l_{2}$-regularized logistic regression problem and the $l_{2}$-regularized Nesterov’s lower-bound function. The logistic regression problem can be formulated as

$$ $f(x)=\tfrac{1}{n}\textstyle{\sum}_{i=1}^{n}\log\left(1+e^{-b_{i}\left\langle a _{i},x\right\rangle}\right)+\tfrac{\mu}{2}\|x\|_{2}^{2},\vspace{-.1cm}$ (21) $$

where $a_{i}\in\mathbb{R}^{d}$ are data features and $b_{i}\in\left\{-1;1\right\}$ are data labels for $i=1,\ldots,n$. The $l_{2}$-regularized third order Nesterov’s lower-bound function from (Nesterov, ) has the next form

$$ $f(x)=\tfrac{1}{4}\textstyle{\sum}_{i=1}^{d-1}(x_{i}-x_{i+1})^{4}-x_{1}+\tfrac{ \mu}{2}\|x\|_{2}^{2}.\vspace{-.1cm}$ (22) $$

Next, we present two figures illustrating the convergence of the methods for both functions.

Figure: (a) Logistic Regression for a9a dataset from the starting point $x_{0}=3e$, where the regularizer $\mu=1e-4$ and $e$ is a vector of all ones.
Refer to caption: x1.png

Figure 1 illustrates that both the Cubic Newton method and Basic Tensor method have areas of superlinear convergence where the graphics are going down faster with each iteration (concave downward). In contrast, gradient descent demonstrates linear convergence.
To verify the behavior of these methods, we plot the values $\alpha_{t}=1-\tfrac{f(x_{t+1})-f^{\ast}}{f(x_{t})-f^{\ast}}\leq 1$, which typically represent the per-iteration decrease rate from $f(x_{t+1})-f^{\ast}\leq(1-\alpha_{t})\left(f(x_{t})-f^{\ast}\right)$. The bigger $\alpha_{t}$ means faster convergence. In the classical proofs of linear convergence, first, the first step is to show that there exists some $\alpha\leq\alpha_{t}$ such that $f(x_{t+1})-f^{\ast}\leq(1-\alpha)\left(f(x_{t})-f^{\ast}\right)$. Next, we derive that $f(x_{T+1})-f^{\ast}\leq(1-\alpha)^{T}\left(f(x_{0})-f^{\ast}\right)\leq e^{-
\alpha\cdot T}\left(f(x_{0})-f^{\ast}\right)$. From the last inequality, we can estimate the number of iterations $T=O\left(\alpha\log\left(\tfrac{f(x_{0})-f^{\ast}}{\varepsilon}\right)\right)$ to reach $\varepsilon$-solution, where $f(x_{T+1})-f^{\ast}\leq\varepsilon$. For gradient descent with $M_{1}=L_{1}$ in (17), $\alpha=\kappa=\tfrac{\mu}{\mu+L_{1}}$.

Figure: (a) Logistic Regression for a9a dataset from the starting point $x_{0}=3e$, where the regularizer $\mu=1e-4$ and $e$ is a vector of all ones.
Refer to caption: x3.png

In Figure 2, we observe that at the beginning all methods slow down for both cases because a smaller $\alpha_{t}=1-\tfrac{f(x_{t+1})-f^{\ast}}{f(x_{t})-f^{\ast}}$ indicates slower convergence. This phase corresponds to the region where the sublinear rate outperforms the linear rate. For example, in the case of gradient descent, this occurs when the guarantee $f(x_{t+1})\leq f(x_{t})-\tfrac{1}{2L_{1}}\|\nabla f(x_{t+1})\|^{2}$ is better than $f(x_{t+1})-f^{\ast}\leq\left(1-\tfrac{\mu}{\mu+L_{1}}\right)\left(f(x_{t})-f^{
\ast}\right)$. As iterations proceed, gradient descent stabilizes around $\alpha_{t}=10^{-3}$, which corresponds to the theoretical convergence rate $\kappa$. The Cubic Newton method and the Basic Tensor method, however, start to accelerate and switch to a superlinear convergence rate. This practical performance gives the intuition for the global superlinear convergence of high-order methods.
Now, we present the theoretical results demonstrating that basic high-order methods indeed have a global superlinear convergence for $\mu$-strongly star-convex functions. We start with a simplified version of the theorem which includes the linear convergence and then we present the full version.

###### Theorem 3.2

For $\mu$-strongly star-convex (3) function $f$ with $L_{2}$-Lipschitz-continuous Hessian (2), Cubic Regularized Newton Method from (18) with $M_{2}\geq L_{2}$ converges with the rate

$$ $f(x_{t+1})-f^{\ast}\leq(1-\alpha_{t})\left(f(x_{t})-f^{\ast}\right),$ (23) $$

for all

$$ $\alpha_{t}\in\left[0;\alpha_{t}^{\ast}\right],\text{ where }\alpha_{t}^{\ast}= \tfrac{-1+\sqrt{1+4\kappa_{t}}}{2\kappa_{t}}\quad\text{and}\quad\kappa_{t}= \tfrac{(M_{2}+L_{2})\|x_{t}-x^{\ast}\|}{3\mu}.$ (24) $$

This range includes the classical linear rate

$$ $f(x_{t})-f(x^{\ast})\leq(1-\alpha^{low})^{t}\left(f(x_{0})-f(x^{\ast})\right) \quad\text{for}\quad\alpha^{low}=\min\left\{\tfrac{1}{2};\sqrt{\tfrac{3\mu}{4( M_{2}+L_{2})D}}\right\}\vspace{-.1cm}$ (25) $$

Proof.
We start the proof by using an upper-bound (14)

$$ $\displaystyle f(x_{t+1})\stackrel{{\scriptstyle\eqref{eq:upper-bound}}}{{\leq} }\Phi_{x_{t},2}(x_{t+1})+\tfrac{L_{2}}{6}\|x_{t+1}-x_{t}\|^{3}\stackrel{{ \scriptstyle\eqref{eq:cubic}}}{{\leq}}\min\limits_{y\in\mathbb{R}^{n}}\left\{ \Phi_{x_{t},2}(y)+\tfrac{M_{2}}{6}\|y-x_{t}\|^{3}\right\}$ $\displaystyle\stackrel{{\scriptstyle\eqref{eq:upper-bound}}}{{\leq}}\min \limits_{y\in\mathbb{R}^{n}}\left\{f(y)+\tfrac{M_{2}+L_{2}}{6}\|y-x_{t}\|^{3} \right\}\stackrel{{\scriptstyle y=x_{t}+\alpha_{t}(x^{\ast}-x_{t})}}{{\leq}}f( (1-\alpha_{t})x_{t}+\alpha_{t}x^{\ast})+\alpha_{t}^{3}\tfrac{M_{2}+L_{2}}{6}\| x^{\ast}-x_{t}\|^{3}$ $\displaystyle\stackrel{{\scriptstyle\eqref{eq:star-convex-strong}}}{{\leq}}(1- \alpha_{t})f(x_{t})+\alpha_{t}f(x^{\ast})-\tfrac{\alpha_{t}(1-\alpha_{t})\mu}{ 2}\|x_{t}-x^{\ast}\|^{2}+\alpha_{t}^{3}\tfrac{M_{2}+L_{2}}{6}\|x_{t}-x^{\ast} \|^{3}.$ $$

From the second inequality, we get that the method is monotone and $f(x_{t+1})\leq f(x_{t})$. Now, by subbing $f(x^{\ast})$ from both sides, we get

$$ $f(x_{t+1})-f(x^{\ast})\leq(1-\alpha_{t})\left(f(x_{t})-f(x^{\ast})\right)- \tfrac{\alpha_{t}}{2}\|x_{t}-x^{\ast}\|^{2}\left((1-\alpha_{t})\mu-\alpha_{t}^ {2}\tfrac{M_{2}+L_{2}}{3}\|x_{t}-x^{\ast}\|\right),\vspace{-.1cm}$ $$

By choosing $\alpha_{t}$ such that

$$ $\alpha_{t}^{2}\tfrac{M_{2}+L_{2}}{3}\|x_{t}-x^{\ast}\|+\mu\alpha_{t}-\mu\leq 0,$ (26) $$

we get (23).
By solving the quadratic inequality (26), we get that the method (9) converges with the rate (23) for all (24). Next, we present Lemma 3.3 with the useful properties of $\alpha_{t}^{\ast}$ from (24). The more general Lemma B.2 with the detailed proof is in Appendix B.

###### Lemma 3.3

For $z>0$, the function

$$ $\alpha^{\ast}(z)=\tfrac{-1+\sqrt{1+4z}}{2z}\vspace{-.1cm}$ (27) $$

has next constant lower and upper-bound

$$ $\min\left\{1,\tfrac{1}{\sqrt{z}}\right\}>\alpha^{\ast}(z)>\min\left\{\tfrac{1} {2};\tfrac{1}{2\sqrt{z}}\right\},\vspace{-.1cm}$ (28) $$

and it is monotonically decreasing

$$ $\forall z,y>0:\quad z<y\quad\Rightarrow\quad\alpha^{\ast}(z)>\alpha^{\ast}(y). \vspace{-.1cm}$ (29) $$

The convergence rate is well-defined as
$0<\alpha_{t}^{\ast}\leq 1$ from (28). As $\|x_{t}-x^{\ast}\|\leq D$ from (5) and $\alpha^{low}\leq\alpha^{\ast}$ by (28), we get the linear convergence rate (25).
$\Box$

Now, we move to the second Theorem and prove the global superlinear convergence. The main idea of the proof is to notice that $\|x_{t}-x^{\ast}\|$ in (24) is decreasing for $\mu$-strongly star-convex functions and hence we can exploit it and show that $\kappa_{t}$ is decreasing and hence $\alpha_{t}^{\ast}$ is increasing from (29), which gives us a superlinear convergence.

###### Theorem 3.4

For $\mu$-strongly star-convex (3) function $f$ with $L_{2}$-Lipschitz-continuous Hessian (2), Cubic Regularized Newton Method from (18) with $M_{2}\geq L_{2}$ converges globally superlinearly as defined in (20) with $\zeta_{t}=1-\alpha_{t}^{sl}$

$$ $f(x_{t+1})-f^{\ast}\leq(1-\alpha_{t}^{sl})\left(f(x_{t})-f^{\ast}\right),$ (30) $$

where

$$ $\alpha_{t}^{sl}=\tfrac{-1+\sqrt{1+4\kappa_{t}^{sl}}}{2\kappa_{t}^{sl}}\quad \text{for}\quad\kappa_{t}^{sl}=\tfrac{(M_{2}+L_{2})\sqrt{2}}{3\mu^{3/2}}(1- \alpha^{low})^{t/2}\left(f(x_{0})-f(x^{\ast})\right)^{1/2}.$ (31) $$

The aggregated convergence rate for $T\geq 1$ equals to

$$ $f(x_{T})-f(x^{\ast})\leq(f(x_{0})-f(x^{\ast}))\textstyle{\prod}_{t=1}^{T}(1- \alpha^{sl}_{t}).$ (32) $$

Proof.
From $\mu$-strongly star-convexity (3), we can upper-bound $\|x_{t}-x^{\ast}\|$ in (24) by

$$ $\displaystyle\|x_{t}-x^{\ast}\|\leq\left(\tfrac{2}{\mu}\left(f(x_{t})-f(x^{ \ast})\right)\right)^{1/2}\stackrel{{\scriptstyle\eqref{eq:al_low}}}{{\leq}} \left(\tfrac{2}{\mu}\left((1-\alpha^{low})^{t}(f(x_{0})-f(x^{\ast}))\right) \right)^{1/2}.$ $$

So, we got that $\|x_{t}-x^{\ast}\|$ is linearly decreasing to zero.
From that, we get a new superlinear $\alpha_{t}^{sl}\leq\alpha_{t}^{\ast}$ from (31).
As $\kappa_{t}^{sl}$ is getting smaller within each iteration $\kappa_{t}^{sl}>\kappa_{t+1}^{sl}$, we get that $\alpha(\kappa_{t}^{sl})<\alpha(\kappa_{t+1}^{sl})$ from (29). Finally, for $\zeta_{t}=1-\alpha(\kappa_{t}^{sl})$, we get $\zeta_{t}>\zeta_{t+1}$ in (30). This finishes the proof of global superlinear convergence. The aggregated convergence rate is equal to (32).
$\Box$

Similar results hold for Basic Tensor methods from (16) in general for $p\geq 2$. Next, we present the theorem for global superlinear convergence of Basic Tensor methods.

###### Theorem 3.5

For $\mu_{q}$-uniformly star-convex (3) function $f$ of degree $q\geq 2$ with $L_{p}$ - Lipschitz-continuous $p$-th derivative ($p\geq q\geq 2$) (2), Basic Tensor Method from (16) with $M_{p}\geq pL_{p}$ converges converges globally superlinearly as defined in (20) with $\zeta_{t,p}=1-\alpha_{t,p}^{sl}$

$$ $f(x_{t+1})-f^{\ast}\leq(1-\alpha_{t,p}^{sl})\left(f(x_{t})-f^{\ast}\right),$ (33) $$

where $\alpha_{t,p}^{sl}$ is such that

$$ $\begin{gathered}h_{\kappa_{t,p}^{sl}}(\alpha_{t,p}^{sl})=0,\quad\text{where} \quad h_{\kappa}(\alpha)=\alpha^{p}\kappa+\alpha-1,\quad\alpha^{low}_{p}=\min \left\{\tfrac{1}{2};\tfrac{1}{2}\left(\tfrac{(p+1)!\mu}{q(M_{p}+L_{p})D^{p-q+1 }}\right)^{1/p}\right\}\\ \text{and}\quad\kappa_{t,p}^{sl}=\tfrac{(M_{p}+L_{p})q^{(q+1)/q}}{(p+1)!\mu^{( q+1)/q}}(1-\alpha^{low}_{p})^{t/q}\left(f(x_{0})-f(x^{\ast})\right)^{1/q}.\end {gathered}$ (34) $$

The aggregated convergence rate for $T\geq 1$ equals to

$$ $f(x_{T})-f(x^{\ast})\leq(f(x_{0})-f(x^{\ast}))\textstyle{\prod}_{t=1}^{T}(1- \alpha^{sl}_{t,p}).$ (35) $$

To sum up, we showed the global superlinear convergence of Cubic Regularized Newton Method for $\mu$-strongly star-convex functions and Basic Tensor methods for $\mu_{q}$-uniformly star-convex functions. More details and proofs are presented in the Appendix B.

###### Definition 3.1

###### Theorem 3.2

###### Lemma 3.3

###### Theorem 3.4

###### Theorem 3.5

## 4 Nesterov Accelerated Tensor Method with A t subscript 𝐴 𝑡 A_{t} italic_A start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT -Adaptation (NATA)

In this section, we focus on Nesterov Acceleration for tensor methods proposed in (Nesterov, ; ). Theoretically, the Nesterov Accelerated Tensor Method (NATM), with a covergence rate $O\left(T^{-(p+1)}\right)$, dominates Basic Tensor Method with the rate $O\left(T^{-p}\right)$. However, in practice, it is often observed that the NATM method is slower, particularly in the initial stages (Scieur, ; Carmon et al., ; Antonakopoulos et al., ). To illustrate this, we present a practical example using the logistic regression problem (21) with $\mu=0$.
In Figure 3(a), the accelerated versions appear slower, which contradicts the theoretical expectations. In this section, we investigate the causes of this underperformance and propose a solution. We begin by revisiting the classical Nesterov Accelerated method, with further details provided in Appendix C.2.

Figure: Algorithm 1 Nesterov Accelerated Tensor Method

According to the theoretical convergence result $f(x_{t})-f(x^{\ast})\leq\tfrac{\|x^{\ast}-x_{0}\|^{p+1}}{(p+1)A_{t}}$ from (Nesterov, , Theorem 2.3), the sequence $A_{t}$ is directly connected with the method’s performance - the larger the $A_{t}$, the faster the convergence. Therefore, our goal is to maximize $A_{t}$. Theoretically, $A_{t}$ should be defined as $A_{t}=\tfrac{\nu_{p}}{L_{p}}t^{p+1},$ where $\nu_{2}=\tfrac{1}{24}$ for $M_{2}=L_{2}$ and $\nu_{3}=\tfrac{5}{3024}$ for $M_{3}=6L_{3}$. However, the values of $\nu_{p}$ appear to be quite small, which limits the speed of convergence. Can these values be increased? The answer is yes.
We propose the Nesterov Accelerated Tensor Method with $A_{t}$-Adaptation, which selects these parameters more aggressively, leading to faster convergence.

Figure: Algorithm 2 Nesterov Accelerated Tensor Method with $A_{t}$-Adaptation (NATA)

###### Theorem 4.1

For convex function $f$ with $L_{p}$-Lipschitz-continuous $p$-th derivative, to find $x_{T}$ such that $f(x_{T})-f(x^{\ast})\leq\varepsilon$, it suffices to perform no more than $T\geq 1$ iterations of the Nesterov Accelerated Tensor Method with $A_{t}$-Adaptation (NATA) with $M_{p}\geq pL_{p}$ (Algorithm 2), where

$$ $T=O\left(\left(\tfrac{M_{p}R^{p+1}}{\varepsilon}\right)^{\tfrac{1}{p+1}}+\log_ {\theta}\left(\tfrac{\nu^{\max}}{\nu^{\min}}\right)\right).$ (36) $$

The proof is presented in the Appendix C.3. Next, we demonstrate the practical improvements of NATA compared to the classical methods. As shown in Figure 3(b), one can see that the Cubic and Tensor variants of NATA significantly outperform the classical Basic and Nesterov Accelerated Methods. We also included a version of Cubic NATA with fixed $\nu_{t}=10$, where $\nu_{t}$ is a tunable hyperparameter. This more aggressive variant of NATA can exhibit even faster practical performance, though it may diverge if $\nu_{t}$ is not chosen carefully.

Figure: (a) Comparison of Basic Methods vs Nesterov Accelerated Methods
Refer to caption: x5.png

###### Theorem 4.1

## 5 Computational Comparison of Acceleration Methods

In this section, we present the practical comparison of different acceleration techniques for tensor methods for convex optimization. We choose five main variants: Nesterov Accelerated Tensor Method (Algorithm 1) with a rate $O(T^{-(p+1)})$ (Nesterov, ); Near-Optimal Tensor Acceleration (Algorithm 5) with a rate $\tilde{O}(T^{-(3p+1)/2})$ (Bubeck et al., ; Gasnikov et al., ; Kamzolov, ); Near-Optimal Proximal-Point Acceleration Method with Segment Search (Algorithm 6) with the rate $\tilde{O}(T^{-(3p+1)/2})$ (Nesterov, );
and Optimal Acceleration (Algorithm 7) with a rate $O(T^{-(3p+1)/2})$ (Kovalev and Gasnikov, ), where $\tilde{O}(\cdot)$ means up to a logarithmic factor. Next, we present experiments for logistic regression (21) on the a9a dataset in Figure 4 and on the w8a dataset in Figure 5.

Figure: Figure 4: Comparison of different cubic and tensor acceleration methods on Logistic Regression for a9a dataset from the starting point $x_{0}=3e$, where $e$ is a vector of all ones.
Refer to caption: x7.png

Let us now discuss the performance of the methods. The new NATA acceleration outperforms all other methods. We attribute this to NATA’s strategy of maximizing $\tilde{a}_{t}$ and $\tilde{A}_{t}$, which enables even faster convergence in the later stages. The second-best performer is the Near-Optimal Tensor Acceleration method. Although it struggles initially due to a large number of line-search iterations per step, it gradually requires fewer line-search iterations—less than two per step on average—as parameters from previous line-search steps become well-suited for the current iteration. With fewer line-search iterations, the method accelerates and outpaces the remaining competitors. A promising direction for improving this method would be to refine the line-search process through an advanced line-search strategy. Next, the classical Accelerated method starts off slower than the basic method without acceleration. Eventually, the method accelerates and overtakes the basic version. Near-Optimal Proximal-Point Acceleration Method with Segment Search performs very similarly to Basic Methods. It has much fewer iterations, but it does a safe segment search with an average of 3 Basic steps per search. Lastly, the Optimal Acceleration method performs the worst in practice. We believe the main issue lies in the internal parameters, which need tuning and adaptation, as we used the theoretical parameters in our implementation. This leads to many inner iterations without significant global progress. Improving these parameters presents an open question for future research.
More details can be found in the Appendix D.

Figure: Figure 5: Comparison of different cubic and tensor acceleration methods on Logistic Regression for w8a dataset from the starting point $x_{0}=3e$, where $e$ is a vector of all ones.
Refer to caption: x9.png

## 6 Conclusion

Limitations. This paper primarily focuses on high-order methods which come with certain limitations. First of all, they have computational and memory limitations in high-dimensional spaces, due to the need for Hessian calculations. There are, however, approaches to overcome this, such as using first-order subsolvers,
or inexact Hessian approximations like Quasi-Newton approximations (BFGS, L-SR1).
In this paper, we focus on the exact Hessian to analyze methods’ peak performance.
Another limitation arises from the specific function classes and the theoretical results considered. Nonetheless, many of the proposed methods can be practically applied to a broader set of problems. For instance, the CRN performs competitively from general non-convex to strongly convex functions.

Conclusion and Future work. In the paper, we demonstrated that the basic high-order methods exhibit global superlinear convergence for $\mu$-strongly star-convex functions. This result is significant because it shows that high-order methods accelerate with each iteration, in stark contrast to first-order methods, which typically have a steady linear convergence rate. This raises intriguing questions: Can superlinear convergence rates be established for accelerated high-order methods as well? What is the best possible global per-iteration decrease that we can theoretically guarantee?
In the second part of the paper, we proposed NATA, a practical acceleration technique. NATA employs a more aggressive schedule adaptation for $A_{t}$, enabling faster convergence. Our experimental results show that NATA significantly outperforms both basic and accelerated methods, including near-optimal and optimal methods. This opens up another interesting question: Could other high-order methods be optimized by addressing practical issues that arise due to overly conservative theoretical guarantees?
Finally, we introduced OPTAMI, an open-source library designed to make high-order optimization methods more accessible and easier to experiment with. We plan to expand this library to cover a wider range of settings and incorporate more optimization methods in the future.

## Appendix A Related Works

The origins of Newton method trace back to the foundational works on root-finding algorithms Newton [], Raphson [], Simpson [], and Bennett []. The next breakthrough in applying Newton method to optimization and proving its local quadratic convergence rates was done by Kantorovich [, , , , , , ]. Over the following decades, Newton’s method have been studied in depth, modified and improved in works of Moré [] Griewank []; Nesterov and Polyak []. Today, Newton’s method is widely used in industrial and scientific computing. For a more detailed history of Newton method, see Boris T. Polyak’s paper [Polyak, ].

Recently, second-order methods have taken a new direction in development with the introduction of globally convergent methods achieving convergence rates of $O(T^{-2})$ [Nesterov and Polyak, ] and $O(T^{3})$ [Nesterov, ] convergence rate, surpassing the performance of first-order methods [Nesterov, ]. These advancements were later extended to higher-order (tensor) methods by Baes []. However, the tensor subproblem in these methods is nonconvex, leading to implementation challenges. This issue was addressed by the introduction of the (Accelerated) Tensor Method in [Nesterov, ], which resolved the nonconvexity by increasing the scaling coefficient of the regularization term, making the subproblem convex. The basic $p$-th order Tensor Method achieves a rate of $O\left(T^{-p}\right)$ , while the accelerated version improves this to $O\left(T^{-(p+1)}\right)$. Earlier work by Monteiro and Svaiter [] demonstrated that even faster convergence for second-order methods is possible with the Accelerated Proximal Extragradient method (A-HPE), achieving a rate of $\tilde{O}\left(T^{-7/2}\right)$. Lower bounds for second-order and higher-order methods of $\Omega\left(T^{-(3p+1)/2}\right)$ were established in [Arjevani et al., , Nesterov, ], demonstrating that the A-HPE method is nearly optimal for second-order convex optimization. Subsequently, three independent research groups [Gasnikov et al., , Bubeck et al., , Jiang et al., ] extended the A-HPE framework to develop tensor methods with a convergence rate of $\tilde{O}\left(T^{-(3p+1)/2}\right)$, achieving near-optimal complexity for these higher-order methods. Truly optimal methods with a rate of $O\left(T^{-(3p+1)/2}\right)$ were later proposed in [Kovalev and Gasnikov, , Carmon et al., ]. Moreover, when assuming higher levels of smoothness, *second-order* methods [Nesterov, , , Kamzolov, , Doikov et al., ] have been shown to exceed the established lower complexity bounds for problems with Lipschitz-continuous Hessians. For an in-depth exploration of higher-order methods, see the review in [Kamzolov et al., ].

Since second-order and higher-order methods generally incur greater computational costs due to the need for calculating higher-order derivatives, it is natural to consider inexact or stochastic algorithms to reduce these overheads. In convex optimization, several studies have explored globally convergent second-order methods with inexact Hessians [Ghadimi et al., ], higher-order methods with inexact and stochastic derivatives [Agafonov et al., , Kamzolov et al., ], and adaptive stochastic methods [Antonakopoulos et al., ]. In [Agafonov et al., ], a lower bound of
$\Omega\left(\frac{\sigma_{1}}{\sqrt{T}}+\frac{\sigma_{2}}{T^{2}}+\frac{1}{T^{7
/2}}\right)$ was established for stochastic globally convergent second-order methods, where $\sigma_{1}$ and $\sigma_{2}$ represent the variances of the stochastic gradients and Hessians, respectively. Additionally, the Accelerated Stochastic Cubic Newton method was introduced, achieving a convergence rate of $O\left(\frac{\sigma_{1}}{\sqrt{T}}+\frac{\sigma_{2}}{T^{2}}+\frac{1}{T^{3}}\right)$, which, to the best of our knowledge, represents the state-of-the-art result. Inexact second-order methods enable the use of Quasi-Newton Hessian approximations, which are well-regarded for their strong practical performance. Although classical Quasi-Newton (QN) methods are known for local superlinear convergence but lack global convergence, their integration with cubic regularization has led to globally convergent methods that also feature relatively inexpensive subproblem solutions [Kamzolov et al., , Scieur, , Jiang et al., ]. Second-order methods with inexact or stochastic derivatives also hold promise for distributed optimization [Shamir et al., , Reddi et al., , Zhang and Lin, , Daneshmand et al., , Agafonov et al., , Dvurechensky et al., , Agafonov et al., ], offering an effective way to manage the computational demands typically encountered in distributed settings. One actively developing direction relies on the constructions of Cubic Newton with explicit steps in order to reduce the complexity of solving methods’ subproblems [Polyak, , , Mishchenko, , Doikov and Nesterov, , Doikov et al., , Hanzely et al., ].

## Appendix B Global Superlinear convergence

In this section, we show the theoretical global superlinear convergence of high-order methods ($p\geq 2$) for $\mu$-strongly star-convex functions.

###### Theorem B.1

For $\mu_{q}$-uniformly star-convex (3) function $f$ of degree $q\geq 2$ with $L_{p}$ - Lipschitz-continuous $p$-th derivative ($p\geq q\geq 2$) (2), Basic Tensor Method from (16) with $M_{p}\geq(p-1)L_{p}$ converges with the rate

$$ $f(x_{t+1})-f^{\ast}\leq(1-\alpha_{t,p})\left(f(x_{t})-f^{\ast}\right),$ (37) $$

for all $\alpha_{t,p}^{\ast}\geq\alpha_{t,p}\geq 0$ such that

$$ $h_{\kappa_{t,p}}(\alpha_{t,p})\leq 0\,\,\text{and}\,\,h_{\kappa_{t,p}}(\alpha_ {t,p}^{\ast})=0,\,\,\text{where}\,\,h_{\kappa}(\alpha)=\alpha^{p}\kappa+\alpha -1\,\,\text{and}\,\,\kappa_{t,p}=\tfrac{q(M_{p}+L_{p})\|x_{t}-x^{\ast}\|^{p-q+ 1}}{(p+1)!\mu}.$ (38) $$

This range includes the classical linear rate

$$ $f(x_{t})-f(x^{\ast})\leq(1-\alpha^{low}_{p})^{t}\left(f(x_{0})-f(x^{\ast}) \right)\,\,\text{for}\,\,\alpha^{low}_{p}=\min\left\{\tfrac{1}{2};\tfrac{1}{2} \left(\tfrac{(p+1)!\mu}{q(M_{p}+L_{p})D^{p-q+1}}\right)^{1/p}\right\}$ (39) $$

Proof.
We start the proof from an upper-bound (14)

$$ $\displaystyle f(x_{t+1})\stackrel{{\scriptstyle\eqref{eq:upper-bound}}}{{\leq} }\Phi_{x_{t},p}(x_{t+1})+\tfrac{L_{p}}{(p+1)!}\|x_{t+1}-x_{t}\|^{p+1}\stackrel {{\scriptstyle\eqref{eq:model_step}}}{{\leq}}\min\limits_{y\in\mathbb{R}^{n}} \left\{\Phi_{x_{t},p}(y)+\tfrac{M_{p}}{(p+1)!}\|y-x_{t}\|^{p+1}\right\}$ $\displaystyle\stackrel{{\scriptstyle\eqref{eq:upper-bound}}}{{\leq}}\min \limits_{y\in\mathbb{R}^{n}}\left\{f(y)+\tfrac{M_{p}+L_{p}}{(p+1)!}\|y-x_{t}\| ^{p+1}\right\}\stackrel{{\scriptstyle y=x_{t}+\alpha_{t,p}(x^{\ast}-x_{t})}}{{ \leq}}f((1-\alpha_{t,p})x_{t}+\alpha_{t,p}x^{\ast})+\alpha_{t,p}^{p+1}\tfrac{M _{p}+L_{p}}{(p+1)!}\|x^{\ast}-x_{t}\|^{p+1}$ $\displaystyle\stackrel{{\scriptstyle\eqref{eq:star-convex-strong}}}{{\leq}}(1- \alpha_{t,p})f(x_{t})+\alpha_{t,p}f(x^{\ast})-\frac{\alpha_{t,p}(1-\alpha_{t,p })\mu}{q}\|x_{t}-x^{\ast}\|^{q}+\alpha_{t,p}^{p+1}\tfrac{M_{p}+L_{p}}{(p+1)!} \|x^{\ast}-x_{t}\|^{p+1}.$ $$

From the third inequality, we get that the method is monotone and $f(x_{t+1})\leq f(x_{t})$. Next, we subtract $f(x^{\ast})$ from the both sides and get

$$ $f(x_{t+1})-f(x^{\ast})\leq(1-\alpha_{t,p})\left(f(x_{t})-f(x^{\ast})\right)- \tfrac{\alpha_{t,p}}{q}\|x_{t}-x^{\ast}\|^{q}\left((1-\alpha_{t})\mu-\alpha_{t ,p}^{p}\tfrac{q(M_{p}+L_{p})}{(p+1)!}\|x_{t}-x^{\ast}\|^{p+1-q}\right),$ $$

If we choose $\alpha_{t}$ such that

$$ $\alpha_{t,p}^{p}\tfrac{q(M_{p}+L_{p})}{(p+1)!}\|x_{t}-x^{\ast}\|^{p-q+1}+\mu \alpha_{t,p}-\mu\leq 0,$ $$

or equivalent version

$$ $\alpha_{t,p}^{p}\tfrac{q(M_{p}+L_{p})}{(p+1)!\mu_{q}}\|x_{t}-x^{\ast}\|^{p-q+1 }+\alpha_{t,p}-1\leq 0,$ $$

we get (37).
To understand the solutions of such inequality, we present Lemma B.2 with the useful properties. From this Lemma, the convergence rate is well-defined as
$0<\alpha_{t,p}^{\ast}\leq 1$ from (41). As $\|x_{t}-x^{\ast}\|\leq D$ from (5) and $\alpha^{low}_{p}\leq\alpha_{t,p}^{\ast}$ by (41), we get the linear convergence rate (39).
$\Box$

###### Lemma B.2

For $z>0$, the solution $\alpha^{\ast}(z)$ of

$$ $h_{z}(\alpha^{\ast}(z))=0,\qquad\text{where}\qquad h_{z}(\alpha)=\alpha^{p}z+ \alpha-1,$ (40) $$

has the next constant lower and upper-bound

$$ $\min\left\{1,\tfrac{1}{z^{1/p}}\right\}>\alpha^{\ast}(z)>\min\left\{\tfrac{1}{ 2};\tfrac{1}{2z^{1/p}}\right\}.$ (41) $$

This bounds show that, for $z\geq 1$, the solution $\alpha^{\ast}(z)$ is similar to $z^{-1/p}$ up to a constant factor as $z^{-1/p}>\alpha^{\ast}(z)>0.5z^{-1/p}$.

For $z\leq 1$, we get the next improved upper and lower-bound

$$ $1-\tfrac{z}{p+1}\geq\alpha^{\ast}(z)\geq 1-z,$ (42) $$

which means that for $z\rightarrow+0$ we have $\alpha^{\ast}(z)\rightarrow 1$.

The solution $\alpha^{\ast}(z)$ is monotonically decreasing

$$ $\forall z,y>0:\quad z<y\quad\Rightarrow\quad\alpha^{\ast}(z)>\alpha^{\ast}(y).$ (43) $$

Proof.
We start the proof from the upper-bound inequality. Note, the function $h_{z}(\alpha)=\alpha^{p}z+\alpha-1$ is monotonically increasing for $\alpha\geq 0$, as $h_{z}(\alpha)^{\prime}=p\alpha^{p-1}z+1>0$. As $h_{z}(0)=-1$ and $h_{z}(1)=z>0$, we get that the solution is unique and $\alpha^{\ast}(z)\in[0;1]$. Next, for $\alpha=\tfrac{1}{z^{1/p}}$, we have

$$ $h_{z}(\tfrac{1}{z^{1/p}})=\tfrac{1}{z}z+\tfrac{1}{z^{1/p}}-1=\tfrac{1}{z^{1/p} }>0,$ $$

which means that $\min\left\{1,\tfrac{1}{z^{1/p}}\right\}>\alpha^{\ast}(z)$ and we proved the upper-bound.

Next, we move to the lower-bound inequality. For $p\geq 2$, $z\geq 1$ and $\alpha=\tfrac{1}{2z^{1/p}}$, we have

$$ $h_{z}(\tfrac{1}{2z^{1/p}})=\tfrac{1}{2^{p}z}z+\tfrac{1}{2z^{1/p}}-1=\tfrac{1}{ 2^{p}}+\tfrac{1}{2z^{1/p}}-1\leq\tfrac{1}{2^{p}}+\tfrac{1}{2}-1<0,$ $$

where the first inequality is coming from $z\geq 1$. The second part of lower-bound holds for $0<z<1$ because

$$ $h_{z}(\tfrac{1}{2})=\tfrac{1}{2^{p}}z+\tfrac{1}{2}-1<\tfrac{1}{2^{p}}-\tfrac{1 }{2}<0.$ $$

We proved (40). Now, to understand the behavior of $\alpha^{\ast}(z)$ for $0<z\leq 1$, we improve the upper and lower-bound for $0<z\leq 1$. For $0<z\leq 1$ and $\alpha=1-z$, we get the improved lower-bound

$$ $h_{z}(1-z)=(1-z)^{p}z+1-z-1=(1-z)^{p}z-z=((1-z)^{p}-1)z<0.$ $$

For $p\geq 2$, $0<z\leq 1$ and $\alpha=1-\tfrac{z}{p+1}$, we get the improved upper-bound

$$ $\displaystyle h_{z}\left(1-\tfrac{z}{p+1}\right)=\left(1-\tfrac{z}{p+1}\right) ^{p}z+1-\tfrac{z}{p+1}-1=\left(1-\tfrac{z}{p+1}\right)^{p}z-\tfrac{z}{p+1}$ $\displaystyle\geq\left(\left(1-\tfrac{1}{p+1}\right)^{p}-\tfrac{1}{p+1}\right) z=\left(\tfrac{p^{p}-(p+1)^{p-1}}{(p+1)^{p}}\right)z>0,$ (44) $$

where to use the last inequity or $p\geq 2$ we need to use some additional analysis. We introduce an additional function and its derivatives

$$ $\displaystyle s(x)=x\log(x)-(x-1)\log(x+1),$ $\displaystyle s(x)^{\prime}=\tfrac{2}{1+x}+\log\left(1-\tfrac{1}{1+x}\right),$ $\displaystyle s(x)^{\prime\prime}=\tfrac{1-x}{x(1+x)^{2}}.$ $$

It is clear that $s(x)^{\prime\prime}<0$ for $x>1$. It means that $s(x)^{\prime}$ is monotonically decreasing. $s(1)^{\prime}=1-log(2)>0$ and the limit $\lim_{x\rightarrow+\infty}s(x)^{\prime}=0$, hence $s(x)^{\prime}\geq 0$ and $s(x)$ is a monotonically increasing function. $s(1)=0$, hence $s(x)>0$ for $x>1$ and finally $x^{x}>(x+1)^{x-1}$ for $x>1$, which proves (44) and finishes the proof of the improved upper-bound (42).

Finally, we show that the solution $\alpha^{\ast}(z)$ is monotonically decreasing with $z$. Let $0<z<y$ and $\alpha^{\ast}(z)$ and $\alpha^{\ast}(y)$ are such that $h_{z}(\alpha^{\ast}(z))=0$ and $h_{y}(\alpha^{\ast}(y))=0$, then

$$ $\alpha^{\ast}(z)^{p}y+\alpha^{\ast}(z)-1\overset{\eqref{eq:al_z_tensor}}{=} \alpha^{\ast}(z)^{p}y+1-\alpha^{\ast}(z)^{p}z-1=\alpha^{\ast}(z)^{p}(y-z)>0,$ $$

which proves that $\alpha^{\ast}(z)>\alpha^{\ast}(y)$ and hence the solution $\alpha^{\ast}(z)$ is monotonically decreasing.
$\Box$

Now, we proceed to the second theorem to establish the global superlinear convergence of high-order methods. The key idea behind the proof is to observe that $|x_{t}-x^{\ast}|$ in (38) decreases for $\mu_{q}$-uniformly star-convex functions. This allows us to notice the fact that $\kappa_{t,p}$ is also decreasing, hence $\alpha_{t,p}$ increases according to (43), ultimately leading to superlinear convergence.

###### Theorem B.3

(Copy of Theorem 3.5)
For $\mu_{q}$-uniformly star-convex (3) function $f$ of degree $q\geq 2$ with $L_{p}$ - Lipschitz-continuous $p$-th derivative ($p\geq q\geq 2$) (2), Basic Tensor Method from (16) with $M_{p}\geq(p-1)L_{p}$ converges converges globally superlinearly as defined in (20) with $\zeta_{t,p}=1-\alpha_{t,p}^{sl}$

$$ $f(x_{t+1})-f^{\ast}\leq(1-\alpha_{t,p}^{sl})\left(f(x_{t})-f^{\ast}\right),$ (45) $$

where $\alpha_{t,p}^{sl}$ is such that

$$ $\begin{gathered}h_{\kappa_{t,p}^{sl}}(\alpha_{t,p}^{sl})=0,\quad\text{where} \quad h_{\kappa}(\alpha)=\alpha^{p}\kappa+\alpha-1,\quad\alpha^{low}_{p}=\min \left\{\tfrac{1}{2};\tfrac{1}{2}\left(\tfrac{(p+1)!\mu}{q(M_{p}+L_{p})D^{p-q+1 }}\right)^{1/p}\right\}\\ \text{and}\quad\kappa_{t,p}^{sl}=\tfrac{(M_{p}+L_{p})q^{(q+1)/q}}{(p+1)!\mu^{( q+1)/q}}(1-\alpha^{low}_{p})^{t/q}\left(f(x_{0})-f(x^{\ast})\right)^{1/q}.\end {gathered}$ (46) $$

The aggregated convergence rate for $T\geq 1$ equals to

$$ $f(x_{T})-f(x^{\ast})\leq(f(x_{0})-f(x^{\ast}))\textstyle{\prod}_{t=1}^{T}(1- \alpha^{sl}_{t,p}).$ (47) $$

Proof.
From $\mu$-uniform star-convexity (3), we can upper-bound $\|x_{t}-x^{\ast}\|$ in (38) by

$$ $\displaystyle\|x_{t}-x^{\ast}\|\leq\left(\tfrac{q}{\mu}\left(f(x_{t})-f(x^{ \ast})\right)\right)^{1/q}\stackrel{{\scriptstyle\eqref{eq:al_low_tensor}}}{{ \leq}}\left(\tfrac{q}{\mu}\left((1-\alpha^{low})^{t}(f(x_{0})-f(x^{\ast})) \right)\right)^{1/q}.$ $$

So, we got that $\|x_{t}-x^{\ast}\|$ is linearly decreasing to zero.
From that, we get a new superlinear $\alpha_{t,p}^{sl}\leq\alpha_{t,p}^{\ast}$ from (46).
As $\kappa_{t}^{sl}$ is getting smaller within each iteration $\kappa_{t,p}^{sl}>\kappa_{t+1,p}^{sl}$, we get that $\alpha(\kappa_{t,p}^{sl})<\alpha(\kappa_{t+1,p}^{sl})$ from (43). Finally, for $\zeta_{t,p}=1-\alpha(\kappa_{t,p}^{sl})$, we get $\zeta_{t,p}>\zeta_{t+1,p}$ in (45). This finishes the proof of global superlinear convergence. The aggregated convergence rate equals to (47).
$\Box$

###### Theorem B.1

###### Lemma B.2

###### Theorem B.3

## Appendix C Methods

### C.1 Subsolover for Basic Tensor Method

In this section, we introduce the subsolver, called Bregman Distance Gradient Method (BDGM), for the Basic Tensor Method of order $p=3$ (19):

$$ $x_{t+1}=x_{t}+\operatornamewithlimits{argmin}_{h\in\mathbb{E}}\left\{f(x_{t})+ \nabla f(x_{t})\left[h\right]+\tfrac{1}{2}\nabla^{2}f(x_{t})\left[h\right]^{2} +\tfrac{1}{6}D^{3}f(x_{t})\left[h\right]^{3}+\tfrac{M_{3}}{24}\|h\|^{4}\right\}.$ (48) $$

The first effective subsolver was introduced by Nesterov in [Nesterov, , Section 5] and later improved in [Nesterov, ]. Next, we describe the BDGM subsolver by following the [Nesterov, ].

#### Relatively inexact p 𝑝 p italic_p -th order solution.

First, we introduce the relatively inexact $p$-th order solution of (16)

$$ $\mathcal{N}^{\gamma}_{M_{p}}(x)=\left\{y\in\mathbb{E}:\|\nabla\Omega_{x,M_{p}} (y)\|_{\ast}\leq\gamma\|\nabla f(y)\|_{\ast}\right\},$ (49) $$

where $\gamma\in[0,1)$ is an accuracy parameter. Then from [Nesterov, , Theorem 2.1], for $\gamma$ and $M_{p}$ such that $\gamma+\tfrac{L_{p}}{M_{p}}\leq\tfrac{1}{p}$ any point $y\in\mathcal{N}^{\gamma}_{M_{p}}(x)$ satisfies

$$ $f(x)-f(y)\geq c_{\gamma,M_{p}}\|\nabla f(y)\|_{\ast}^{\frac{p+1}{p}},\quad \text{where}\quad c_{\gamma,M_{p}}=\left[\tfrac{(1-\gamma)p!}{L_{p}+M_{p}} \right]^{\tfrac{1}{p}}.$ $$

Note, that for the exact solution, we get the same improvement guarantee with $\gamma=0$.
For $p=3$ and (48), we choose $\gamma=1/6$ and $M_{3}=6L_{3}$, then $\mathcal{N}_{L_{3}}(x)=\mathcal{N}^{1/6}_{L_{3}}(x)$ and the method

$$ $x_{t+1}\in\mathcal{N}_{L_{3}}(x_{t})$ $$

converge with the same rate up to a constant as an exact version [Nesterov, , Theorem 2.2]. Note, $M_{3}\geq 3L_{3}$ is also required for the convexity of the subproblem (48). In our implementation, all third-order basic steps are solved with this relative inexactness and $M_{3}=6L_{3}$. This approach creates practical and parameter-free stopping criteria for the subproblem solvers.

#### Relative smoothness and relative strong convexity.

Now, we move on to the concept of relative smoothness and relative strong convexity proposed in [Lu et al., ]. Similarly to classical smoothness and strong convexity, we say that function $\phi(h)$ is relatively $L_{\rho}$-smooth and relatively $\mu_{\rho}$ strongly convex with respect to scaling function $\rho(h)$ if

$$ $\mu_{\rho}\nabla^{2}\rho(h)\preceq\nabla^{2}\phi(h)\preceq L_{\rho}\nabla^{2} \rho(h).$ $$

In classical regime, $\rho(h)=\tfrac{1}{2}\|h\|^{2}$ and $\nabla^{2}\rho(h)$ is an identity matrix. For the scaling function $\rho(h)$, we introduce its Bregman distance

$$ $\beta_{\rho}(h,y)=\rho(y)-\rho(h)-\left\langle\nabla\rho(h),y-h\right\rangle.$ $$

Now the gradient method with respect to this Bregman distance is called Bregman Distance Gradient Method (BDGM) and has the next form

$$ $h_{k+1}=\operatornamewithlimits{argmin}_{y\in\mathbb{E}}\left\{\left\langle \nabla\phi(h_{k}),y-h_{k}\right\rangle+2L_{\rho}\beta_{\rho}(h_{k},y)\right\}.$ $$

The convergence rate of such method is $O\left(\tfrac{L_{\rho}}{\mu_{\rho}}\log\left(\tfrac{\phi(h_{0})-\phi(h^{\ast})
}{\varepsilon}\right)\right)$.

#### Bregman Distance Gradient Method (BDGM) for ( 48 ).

Let’s apply this approach to the solution of subproblem (48) with $M_{3}=6L_{3}$. In [Nesterov, , Section 4], it was shown that the subproblem function $\phi(h)=\nabla f(x_{t})\left[h\right]+\tfrac{1}{2}\nabla^{2}f(x_{t})\left[h
\right]^{2}+\tfrac{1}{6}D^{3}f(x_{t})\left[h\right]^{3}+\tfrac{L_{3}}{4}\|h\|^
{4}$ is relatively smooth and relatively strongly convex with respect to

$$ $\rho(h)=\tfrac{1}{2}\nabla^{2}f(x_{t})\left[h\right]^{2}+\tfrac{L_{3}}{4}\|h\| ^{4}$ $$

with constants $L_{\rho}=1+\tfrac{1}{\sqrt{2}}$ and $\mu_{\rho}=1-\tfrac{1}{\sqrt{2}}$. It means that the method has an incredibly fast convergence rate $O\left(\tfrac{\sqrt{2}+1}{\sqrt{2}-1}\log\left(\tfrac{\phi(h_{0})-\phi(h^{\ast
})}{\varepsilon}\right)\right)$. The details and more formal convergence results are presented in [Nesterov, ].

Now, we present the explicit formulation of the BDGM for (48). First, we have the general form

$$ $h_{k+1}=\operatornamewithlimits{argmin}_{y\in\mathbb{E}}\left\{\left\langle \nabla\phi(h_{k}),y-h_{k}\right\rangle+2L_{\rho}\beta_{\rho}(h_{k},y)\right\}.$ (50) $$

Let us calculate $\nabla\phi(h_{k})$ first. It equals to

$$ $\nabla\phi(h_{k})=\nabla f(x_{t})+\nabla^{2}f(x_{t})h_{k}+\tfrac{1}{2}D^{3}f(x _{t})\left[h_{k}\right]^{2}+L_{3}\|h_{k}\|^{2}h_{k}.$ $$

In [Nesterov, ], the universal approximation for $D^{3}f(x_{t})\left[h_{k}\right]^{2}$ is presented by using the finite differences approach. However, in practice, we recommend using autogradient computation of $D^{3}f(x_{t})\left[h_{k}\right]^{2}$ if it is possible. The computation by autogradient is much more precise while having the same computational complexity. The computation complexity of $D^{3}f(x_{t})\left[h_{k}\right]^{2}$ by autogradient is similar to calculating three gradients as $D^{3}f(x)\left[h\right]^{2}=\nabla_{x}(\nabla^{2}f(x)[h]^{2})=\nabla(\nabla
\left\{\nabla f(x)[h]\right\}[h])$. Also, autogradient computations are commonly used in modern frameworks such as PyTorch, Jax, and others. So, essentially we still have access to third-order information but with the complexity of a gradient computation.

Now, let us calculate explicit $\beta_{\rho}(h_{k},y)$

$$ $\displaystyle\beta_{\rho}(h_{k},y)=\rho(y)-\rho(h)-\left\langle\nabla\rho(h),y -h\right\rangle$ $\displaystyle=\tfrac{1}{2}\nabla^{2}f(x_{t})\left[y\right]^{2}+\tfrac{L_{3}}{4 }\|y\|^{4}-\tfrac{1}{2}\nabla^{2}f(x_{t})\left[h_{k}\right]^{2}-\tfrac{L_{3}}{ 4}\|h_{k}\|^{4}$ $\displaystyle-\left\langle\nabla^{2}f(x_{t})\left[h_{k}\right]+L_{3}\|h_{k}\|^ {2}h_{k},y-h_{k}\right\rangle.$ $$

Note, that the constant terms are useless for finding the argminimum in (50), hence we can remove them. We also can divide all parts of (50) by $2L_{\rho}=2+\sqrt{2}$ for simplicity and unite the linear parts together

$$ $\begin{gathered}g_{k}=\tfrac{1}{2+\sqrt{2}}\nabla\phi(h_{k})-\nabla^{2}f(x_{t} )\left[h_{k}\right]-L_{3}\|h_{k}\|^{2}h_{k}\\ =\tfrac{2-\sqrt{2}}{2}\left(\nabla f(x_{t})+\nabla^{2}f(x_{t})h_{k}+\tfrac{1}{ 2}D^{3}f(x_{t})\left[h_{k}\right]^{2}+L_{3}\|h_{k}\|^{2}h_{k}\right)-\nabla^{2 }f(x_{t})\left[h_{k}\right]-L_{3}\|h_{k}\|^{2}h_{k}\\ =\tfrac{2-\sqrt{2}}{2}\left(\nabla f(x_{t})+\tfrac{1}{2}D^{3}f(x_{t})\left[h_{ k}\right]^{2}\right)-\tfrac{\sqrt{2}}{2}\left(\nabla^{2}f(x_{t})\left[h_{k} \right]+L_{3}\|h_{k}\|^{2}h_{k}\right)\end{gathered}$ $$

So, we finally get the next explicit BDGM step

$$ $h_{k+1}=\operatornamewithlimits{argmin}_{y\in\mathbb{E}}\left\{\left\langle g_ {k},y\right\rangle+\tfrac{1}{2}\nabla^{2}f(x_{t})\left[y\right]^{2}+\tfrac{L_{ 3}}{4}\|y\|^{4}\right\}.$ (51) $$

This step doesn’t require the computation of a full third-order derivative and is similar to the Cubic Regularized Newton step. Hence, we count it as a second-order method. So, the total complexity of Basic Tensor Method for convex functions is $\tilde{O}\left(\tfrac{L_{3}D^{4}}{T^{3}}\right)$ steps of (51), where $\tilde{O}(\cdot)$ means number of iterations up to a logarithmic factor.

#### Inner subsolver for ( 51 ).

The last part is to solve (51). We solve it similarly to the Cubic Regularized Step by ray-search with eigenvalue decomposition (EVD). First, we apply eigenvalue decomposition to $\nabla^{2}f(x_{t})$

$$ $\nabla^{2}f(x_{t})=USU^{\top},$ (52) $$

where $S\in\mathbb{R}^{d\times d}$ is a diagonal matrix with eigenvalues and $U\in\mathbb{R}^{d\times d}$ is an orthoganal matrix such that $UU^{\top}=I$. Then, we denote $v=U^{\top}y$ and $\hat{g}=U^{\top}g_{k}$. Now we can formulate a dual one-dimensional problem.

$$ $\displaystyle\min_{y\in\mathbb{E}}\left\{\left\langle g_{k},y\right\rangle+ \tfrac{1}{2}\left\langle\nabla^{2}f(x_{k})y,y\right\rangle+\tfrac{L_{3}}{4}\|y \|^{4}\right\}$ $\displaystyle=\min_{y\in\mathbb{E}}\left\{\left\langle U^{\top}g_{k},U^{\top}y \right\rangle+\tfrac{1}{2}\left\langle USU^{\top}y,y\right\rangle+\tfrac{L_{3} }{4}\|U^{\top}y\|^{4}\right\}$ $\displaystyle=\min_{v\in\mathbb{E}}\left\{\left\langle\tilde{g},v\right\rangle +\tfrac{1}{2}\left\langle Sv,v\right\rangle+\tfrac{L_{3}}{4}\|v\|^{4}\right\}$ $\displaystyle=\min_{v\in\mathbb{E}}\max_{\tau\geq 0}\left\{\left\langle\tilde{ g},v\right\rangle+\tfrac{1}{2}\left\langle Sv,v\right\rangle+\tfrac{\sqrt{2L_{ 3}}}{2}\|v\|^{2}\tau-\tfrac{1}{2}\tau^{2}\right\}$ $\displaystyle=\max_{\tau\geq 0}\min_{v\in\mathbb{E}}\left\{\left\langle\tilde{ g},v\right\rangle+\tfrac{1}{2}\left\langle Sv,v\right\rangle+\tfrac{\sqrt{2L_{ 3}}}{2}\|v\|^{2}\tau-\tfrac{1}{2}\tau^{2}\right\}$ $\displaystyle=\max_{\tau\geq 0}\left\{-\tfrac{1}{2}\left\langle\left(S+\tau \sqrt{2L_{3}}\right)^{-1}\tilde{g},\tilde{g}\right\rangle-\tfrac{1}{2}\tau^{2} \right\},$ (53) $$

where $\tau^{\ast}=\tfrac{\sqrt{2L_{3}}}{2}\|v\|^{2}$ for the third equality and $v=-\left(S+\tau\sqrt{2L_{3}}\right)^{-1}\tilde{g}$ in the last equality. By solving (53) with one-dimensional ray-search, we find optimal $\tau^{\ast}$ then we can calculate $v$ and $y$, which we found the solution for subproblem (51). In our code, we use eigenvalue decomposition for efficiency of the ray-search, but it is also possible to just inverse the regularized matrix multiple times in (53) or apply some efficient first-order method for quadratic problems such as conjugate gradient.

To finalize, in this section we presented the subsolver which allows us to efficiently implement the Basic Tensor Method for $p=3$ with the complexity same up to a logarithmic factor as Cubic Regularized Newton Method.

### C.2 Nesterov Accelerated Tensor Methods

In this section, we present Nesterov Acceleration for tensor methods proposed in [Nesterov, , ].
First, let us introduce the main parts of the method. The key part of such acceleration is the estimated sequences technique. It is based on linear approximations of function $f(x)$ in a sequence of points $x_{t}$, which allows to construct the estimating function $\psi_{t}(x)$ for a scaling sequence $a_{t}\in\mathbb{R}_{+}$:

$$ $\psi_{t+1}(z)=\psi_{t}(z)+a_{t+1}\left(f(x)+\left\langle\nabla f(x),z-x\right \rangle\right),\quad\text{where}\quad\psi_{0}(z)=\tfrac{1}{p+1}\|z-x_{0}\|^{p+ 1}.$ (54) $$

Additionally, we introduce the sequence

$$ $A_{t+1}=A_{t}+a_{t}.$ (55) $$

Now, we are ready to present the accelerated method.

Figure: Algorithm 3 Nesterov Accelerated Tensor Method

For the convergence results, the sequence $A_{t}$ should be defined in the following way

$$ $A_{t}=\tfrac{\nu_{p}}{L_{p}}t^{p+1},\quad{where}\quad\nu_{p}=\tfrac{2p-1}{(p+1 )(2p+1)}\cdot\tfrac{(p-1)!}{(2p)^{p}}.$ (56) $$

Then, $a_{t+1}=\tfrac{\nu_{p}}{L_{p}}\left((t+1)^{p+1}-t^{p+1}\right)$. With such parameters, we can present the convergence theorem from [Nesterov, , Theorem 2.3]

###### Theorem C.1

Let sequence $\left\{x_{t}\right\}_{t\geq 0}$ be generated by method 3. Then, for any $T\geq 1$, we have

$$ $f(x_{T})-f(x^{\ast})\leq O\left(\tfrac{L_{p}R^{p+1}}{T^{p+1}}\right).$ $$

###### Theorem C.1

### C.3 Nesterov Accelerated Tensor Method with A t subscript 𝐴 𝑡 A_{t} italic_A start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT -Adaptation (NATA)

In this subsection, we present the proof of Theorem 4.1

Figure: Algorithm 4 Nesterov Accelerated Tensor Method with $A_{t}$-Adaptation (NATA)

###### Theorem C.2

(Copy of Theorem 4.1)
For convex function $f$ with $L_{p}$-Lipschitz-continuous $p$-th derivative, to find $x_{T}$ such that $f(x_{T})-f(x^{\ast})\leq\varepsilon$, it suffices to perform no more than $T\geq 1$ iterations of the Nesterov Accelerated Tensor Method with $A_{t}$-Adaptation (NATA) with $M_{p}\geq pL_{p}$ (Algorithm 2), where

$$ $T=O\left(\left(\tfrac{L_{p}R^{p+1}}{\varepsilon}\right)^{\tfrac{1}{p+1}}+\log_ {\theta}\left(\tfrac{\nu^{\max}}{\nu^{\min}}\right)\right).$ (57) $$

Proof.

Let us present the convergence analysis of Algorithm 2. The proof is based on the proof from [Nesterov, ].

First of all, by convexity and definition of $\psi_{t}(x)$, it is easy to show that

$$ $\psi_{t}(x^{\ast})\leq\tilde{A}_{t}f(x^{\ast})+\tfrac{1}{p+1}\|x^{\ast}-x_{0} \|^{p+1}.$ (58) $$

Now, let us assume that the condition on Line 10 is satisfied for every step. Then, we get

$$ $\tilde{A}_{t}f(x_{t})\leq\psi_{t}(v_{t})\leq\psi_{t}(x^{\ast})\leq\tilde{A}_{t }f(x^{\ast})+\tfrac{1}{p+1}\|x^{\ast}-x_{0}\|^{p+1},$ (59) $$

where in the second inequality we use the definition of $v_{t}$. Next, by simple calculations, we get the convergence result

$$ $f(x_{t})-f(x^{\ast})\leq\tfrac{\|x^{\ast}-x_{0}\|^{p+1}}{(p+1)\tilde{A}_{t}}.$ (60) $$

From that inequality, one can see that the bigger $\tilde{A}_{t}$ means the faster convergence. That is the reason, we want to have a more aggressive $\tilde{a}_{t}$ and start the search of $\nu$ from the maximal value. Now, we need to show that the condition in Line 10 is always can be satisfied.

Let us prove it by induction of the following relation:

$$ $\psi^{*}_{t}=\psi_{t}(v_{t})\geq\tilde{A}_{t}f(x_{t}),\quad t\geq 0.$ (61) $$

For $t=0$, we have $\psi^{*}_{0}=0$ and $A_{0}=0$. Hence, (61) is valid.

Assume it is valid for some $t\geq 0$. Then,

$$ $\displaystyle\psi^{*}_{t+1}=\psi_{t}(v_{t+1})+\tilde{a}_{t+1}\left(f(x_{t+1})+ \left\langle\nabla f(x_{t+1}),v_{t+1}-x_{t+1}\right\rangle\right)$ $\displaystyle\geq\psi^{*}_{t}+\tfrac{1}{(p+1)2^{p-1}}\|v_{t+1}-v_{t}\|^{p+1}+ \tilde{a}_{t+1}\left(f(x_{t+1})+\left\langle\nabla f(x_{t+1}),v_{t+1}-x_{t+1} \right\rangle\right),$ $$

where the last inequality is coming from uniform convexity of $\|\cdot\|^{p+1}$.
Now, we can use the structure of the method in previous inequality and get

$$ $\displaystyle\psi^{*}_{t+1}-\tfrac{1}{(p+1)2^{p-1}}\|v_{t+1}-v_{t}\|^{p+1} \overset{\eqref{superfast:induction_tesis}}{\geq}\tilde{A}_{t}f(x_{t})+\tilde{ a}_{t+1}\left(f(x_{t+1})+\left\langle\nabla f(x_{t+1}),v_{t+1}-x_{t+1}\right \rangle\right)$ $\displaystyle\geq\tilde{A}_{t+1}f(x_{t+1})+\langle\nabla f(x_{t+1}),\tilde{a}_ {t+1}(v_{t+1}-x_{t+1})+\tilde{A}_{t}(x_{t}-x_{t+1})\rangle$ $\displaystyle=\tilde{A}_{t+1}f(x_{t+1})+\langle\nabla f(x_{t+1}),\tilde{a}_{t+ 1}(v_{t+1}-v_{t})+\tilde{A}_{t+1}(y_{t}-x_{t+1})\rangle,$ $$

where, for the second inequality, we use convexity and, for the last equality, we use the definition of $y_{t}$ from Line 6 of the Algorithm 2.

Further, we use inequality $\frac{\alpha}{p+1}\tau^{p+1}-\beta\tau\geq-\frac{p}{p+1}\alpha^{-1/p}\beta^{(p
+1)/p},\quad\tau\geq 0$, for all
$x\in\mathbb{E}$ and we have

$$ $\tfrac{1}{(p+1)2^{p-1}}\|v_{t+1}-v_{t}\|^{p+1}+\tilde{a}_{t+1}\left\langle \nabla f(x_{t+1}),v_{t+1}-v_{t}\right\rangle\geq-\tfrac{p}{p+1}2^{\tfrac{p-1}{ p}}\left(\tilde{a}_{t+1}\|\nabla f(x_{t+1})\|_{*}\right)^{\frac{p+1}{p}}.$ (62) $$

Next, for $x_{t+1}\in\mathcal{N}_{L_{p}}(y_{k})$, from [Nesterov, , Theorem 2.1],
we get

$$ $\langle\nabla f(x_{t+1}),y_{t}-x_{t+1}\rangle\geq c_{p}\|\nabla f(x_{t+1})\|^{ \frac{p+1}{p}}_{*},$ $$

where
$c_{p}=\left[\tfrac{2p-1}{2p(2p+1)}\tfrac{p!}{L_{p}}\right]^{1/p}$ for relative inexact $p$-th order solution.

Putting all these inequalities together, we obtain

$$ $\displaystyle\psi^{*}_{t+1}\geq\tilde{A}_{t+1}f(x_{t+1})-\tfrac{p}{p+1}2^{ \tfrac{p-1}{p}}\left(\tilde{a}_{t+1}\|\nabla f(x_{t+1})\|_{*}\right)^{\tfrac{p +1}{p}}+\tilde{A}_{t+1}c_{p}\|\nabla f(x_{t+1})\|^{\tfrac{p+1}{p}}_{*}$ $\displaystyle=\tilde{A}_{t+1}f(x_{t+1})+\|\nabla f(x_{t+1})\|^{\tfrac{p+1}{p}} _{*}\left(\tilde{A}_{t+1}c_{p}-\tfrac{p}{p+1}2^{\tfrac{p-1}{p}}\tilde{a}^{ \tfrac{p+1}{p}}_{t+1}\right).$ $$

Finally, by the choice of $\nu_{t}$ in Algorithm 2, $\nu_{t}\geq\nu_{p}$ and $\tilde{a}_{t+1}\geq a_{t+1}$, where $a_{t+1}=\tfrac{\nu_{p}}{L_{p}}((t+1)^{p+1}-t^{p+1})$ is the theoretical value of $a_{t+1}$. Hence, $\tilde{A}_{t+1}\geq A_{t+1}$, where $A_{t+1}=\tfrac{\nu_{p}}{L_{p}}(t+1)^{p+1}$ is the theoretical value of $A_{t+1}$. So, in the final inequality, we prove that there exists $\nu_{t}=\nu_{p}$ such that

$$ $\tilde{A}_{t+1}c_{p}\geq A_{t+1}c_{p}\geq\tfrac{p}{p+1}2^{\tfrac{p-1}{p}}a^{ \tfrac{p+1}{p}}_{t+1},$ $$

where the last inequality holds from [Nesterov, , Equation 25].
Thus, we have proved the induction step.

The search of $\nu_{t}$ takes maximal total of $\log_{\theta}\left(\tfrac{\nu^{\max}_{t}}{\nu^{\min}_{t}}\right)+T$ additional steps, where $\nu^{\max}_{t}=\max_{t\in[0;T]}\nu_{t}\leq\nu^{\max}$ and $\nu^{\min}_{t}=\min_{t\in[0;T]}\nu_{t}\geq\nu^{\min}=\nu_{p}$. The $T$ term in the sum is coming from Line 11 in Algorithm 2. If we want to make the Algorithm less aggressive, we can remove this Line then $\nu_{t}$ will only decrease.

The total number of iterations hence is equal to $T=O\left(\left(\tfrac{L_{p}R^{p+1}}{\varepsilon}\right)^{\tfrac{1}{p+1}}+\log_
{\theta}\left(\tfrac{\nu^{\max}_{t}}{\nu^{\min}_{t}}\right)\right)$, which finishes the proof.
$\Box$

###### Theorem C.2

### C.4 Near-optimal Tensor methods and Hyperfast Second-order method

#### Near-optimal Tensor methods.

Monteiro and Svaiter [] demonstrated that the global convergence rate of second-order methods can be further improved from $O\left(\varepsilon^{-1/3}\right)$ to $O\left(\varepsilon^{-2/7}\log\left(1/\varepsilon\right)\right)$ . This improvement was achieved through the development of the Accelerated Hybrid Proximal Extragradient (A-HPE) framework, which, when combined with a trust-region Newton-type method, resulted in the Accelerated Newton Proximal Extragradient (A-NPE) method that achieves the improved rate. A lower bound of $O\left(\varepsilon^{-2/7}\right)$ was established by Arjevani et al. [], rendering that the A-NPE method is nearly optimal.

Near-optimal tensor methods [Gasnikov et al., , Bubeck et al., , Jiang et al., ], with a convergence rate of $O\left(\varepsilon^{-2/(3p+1)}\log\left(1/\varepsilon\right)\right)$, are based on the A-HPE framework. Similar to A-HPE, these tensor methods require an additional binary search procedure at each iteration. The cost of these procedures introduces an extra $O(\log(1/\varepsilon))$ factor in the overall convergence rate.

Figure: Algorithm 5 Inexact $p$-th order Near-optimal Accelerated Tensor Method [Kamzolov, , Algorithm 1]

One version of the near-optimal tensor methods is presented in Algorithm 5. This version was initially proposed by Bubeck et al. [] and later improved by Kamzolov [], who introduced the handling of inexact solution to subproblem (64). Note that line (4) of Algorithm 5 requires finding the pair $\left(x_{t+1},\lambda_{t+1}\right)$, which cannot be done explicitly. Specifically, $\lambda_{t+1}$ depends on $x_{t+1}$ via (63), which in turn depends on $y_{t}$ through (64). Furthermore, $y_{t}$ depends on $a_{t+1}$, which itself depends on $\lambda_{t+1}$ as per (65). This recursive dependence implies that $\lambda_{t+1}$ relies on itself, making it impossible to solve in closed form.

To find the pair $\left(x_{t+1},\lambda_{t+1}\right)$, a binary search procedure is employed. Below, we provide the approach used by Bubeck et al. []. Let us denote $\theta=\tfrac{A_{t}}{A_{t+1}}\in[0,1]$. Thus, both $y_{t}$ and $x_{t+1}$ depend on $\theta$,

$$ $y_{t}(\theta):=y_{t}\stackrel{{\scriptstyle\eqref{eq:hyper_a}}}{{=}}\theta x_{ t}+(1-\theta)v_{t},\quad x_{t+1}(\theta):=x_{t+1}\stackrel{{\scriptstyle\eqref {prox_step}}}{{=}}\mathcal{N}_{p,M_{p}}^{\gamma}(y_{t}(\theta)).$ $$

Since $\lambda_{t+1}=\tfrac{a_{t+1}^{2}}{A_{t+1}}$, we have that $\lambda_{t+1}=\tfrac{(1-\theta)^{2}}{\theta}A_{t}$. Thus, in terms of $\theta$, (63) can be rewritten as

$$ $\frac{1}{2}\leq\zeta(\theta)\leq\frac{p}{p+1},\quad\text{where}\quad\zeta( \theta)=\frac{(1-\theta)^{2}}{\theta}\frac{A_{t}M_{p}\cdot\|x_{t+1}(\theta)-y_ {t}(\theta)\|^{p-1}}{(p-1)!}.$ (66) $$

Note that $\zeta(0)\to+\infty$ and $\zeta(1)=0$. Hence, one can use binary search to find $\theta$ such that (66) holds true. The complexity of this procedure is $O\left(\log(1/\varepsilon)\right)$, and a theoretical analysis of binary search procedure can be found in [Bubeck et al., ]. Below we present the total complexity of Algorithm 5.

###### Theorem C.3 ( [Kamzolov, 2020 , Theorem 1] )

For convex function $f$ with $L_{p}$-Lipschitz-continuous $p$-th derivative, to find $x_{T}$ such that $f(x_{T})-f^{*}\leq\epsilon$, it suffices to perform no more than $T\geq 1$ iterations of Algorithm 5 with $H_{p}=\xi L_{p}$, where $\xi$ and $\gamma$ satisfy $1\geq 2\gamma+\tfrac{1}{\xi(p+1)}$, and

$$ $T=\tilde{O}\left(\frac{H_{p}R^{p+1}}{\varepsilon}\right).$ $$

###### Theorem C.3 ( [Kamzolov, 2020 , Theorem 1] )

#### Hyperfast Second-order method.

Interestingly, the lower bound for second-order convex optimization, $O\left(\epsilon^{-2/7}\right)$, can be surpassed under higher smoothness assumptions on the objective. Nesterov [] showed that, under the assumption of an $L_{3}$-Lipschitz third derivative, Algorithm 1 can be implemented using only a second-order oracle, with the third-order derivative approximated via finite gradient differences. This results in a second-order method with $O\left(\epsilon^{-1/4}\right)$ calls to the second-order oracle. The same idea can be applied to Algorithm 1, improving the convergence rate of the second-order method to $\tilde{O}\left(\epsilon^{-1/5}\right)$ [Kamzolov, ].

###### Theorem C.4 ( [Kamzolov, 2020 , Theorem 2] )

For a convex function $f$ with an $L_{3}$-Lipschitz-continuous third derivative, to find $x_{T}$ such that $f(x_{T})-f^{*}\leq\epsilon$, it suffices to perform no more than $N_{1}\geq 1$ gradient calculations and $N_{2}\geq 1$ Hessian calculations in Algorithm 5 with BGDM as the subsolver for the subproblem (64), $H_{p}=3L_{p}/2$, $\gamma=1/6$, and

$$ $N_{1}=\tilde{O}\left(\left(\frac{L_{3}R^{4}}{\epsilon}\right)^{\frac{1}{5}} \log\left(\frac{G+H}{\epsilon}\right)\right),$ $$

$$ $N_{2}=\tilde{O}\left(\left(\frac{L_{3}R^{4}}{\epsilon}\right)^{\frac{1}{5}} \right),$ $$

where $G$ and $H$ are the uniform upper bounds for the norms of the gradients and Hessians computed at the points generated by the main algorithm.

###### Theorem C.4 ( [Kamzolov, 2020 , Theorem 2] )

### C.5 Proximal Point Method with Segment Search

Another approach for constructing near-optimal tensor methods involves high-order proximal-point type methods [Nesterov, , ], which are based on the $p$-th-order proximal-point operator:

$$ $\text{prox}_{p,H}(y)=\operatornamewithlimits{argmin}\limits_{x\in\mathbb{E}} \left\{f_{y,p,H}(x):=f(x)+\tfrac{H}{p+1}\|x-y\|^{p+1}\right\}.$ (67) $$

Nesterov [] demonstrated that using a single step of a $p$-th-order tensor method to solve (67) results in a convergence rate of $O(\epsilon^{-1/p})$, and moreover, this approach can be accelerated to achieve a rate of $O(\epsilon^{-1/(p+1)})$. Another significant contribution of Nesterov [] is the introduction of a proximal-point operator with segment search:

$$ $\text{Sprox}_{p,H}(y,u)=\operatornamewithlimits{argmin}\limits_{x\in\mathbb{E} ,~{}\tau\in[0,1]}\left\{f(x)+\tfrac{H}{p+1}\|x-y-\tau u\|^{p+1}\right\}.$ (68) $$

Assuming that (68) can be solved exactly, Nesterov [] showed that convergence rate of $O(\varepsilon^{-2/(3p+1)})$ can be achieved via different acceleration scheme.

A more practical algorithm was introduced in [Nesterov, ]. Following Nesterov [], the authors assumed that the problem (67) can be solved under the following approximate condition:

$$ $\mathcal{A}_{p,H}^{\gamma}(y)=\left\{x\in\mathbb{E}:\|\nabla f_{y,p,H}(x)\|_{* }\leq\beta\|\nabla f(x)\|_{*}\right\},$ $$

where $\gamma\in[0,1)$ is a tolerance parameter. Furthermore, a specific approach for approximating the solution to subproblem (68) was proposed. The resulting method, called the Inexact $p$-th-order Proximal Point Method with Segment Search, is presented in Algorithm 6. Lines 5-14 of Algorithm 6 detail the steps for the approximate solution of (68).

Figure: Algorithm 6 Inexact $p$-th-order Proximal Point Method with Segment Search [Nesterov, , Method (3.6)]

###### Theorem C.5 ( [Nesterov, 2021a , Theorem 2] )

For smooth convex function $f$ to find $x_{T}$ such that $f(x_{T})-f^{*}\leq\epsilon$, it suffices to perform no more than $T\geq 1$ iterations of Algorithm 6, where

$$ $T=O\left(\left[\frac{HR^{p+1}}{\varepsilon}\right]^{\tfrac{2}{3p+1}}\right).$ $$

Line 12 requires additional bisection search with complexity of $O\left(\frac{HD^{p+1}}{\varepsilon}\right)$ [Nesterov, , Theorem 4]. This results in the following upper bound for the number of evaluations of $w\in\mathcal{A}_{p,H}^{\gamma}(x)$ during the execution of Algorithm 6 $O\left(\left[\frac{HD^{p+1}}{\varepsilon}\right]^{\frac{2}{3p+1}}\log\frac{HD^
{p+1}}{\varepsilon}\right)$.

Under the additional assumption of an $L_{p}$-Lipschitz continuous $p$-th derivative of $f$, the inclusion $w\in\mathcal{A}_{p,H}^{\gamma}(x)$ can be achieved by performing one inexact tensor step with specific choice of parameters $\beta$ and $M_{p}$: $w\in\mathcal{N}_{p,M_{p}}^{\beta}(x)$ [Nesterov, , Section 3] [Nesterov, , Section 5.1]. This makes Algorithm 6 a near-optimal tensor method, comparable to [Gasnikov et al., , Bubeck et al., , Jiang et al., ]. However, it differs in nature: while the latter methods are based on A-NPE-type approaches, Algorithm 6 follows an interior-point-type framework.

For the case when $p=3$, the tensor step can be efficiently performed using BDGM in $O\left(\log 1/\varepsilon\right)$ iterations. As demonstrated in [Nesterov, , Kamzolov, ], a second-order implementation of a third-order tensor method can be achieved by approximating the third-order derivative using finite gradient differences. However, in practice, this approximation may suffer from numerical instability. For Algorithm 6 another approach is available: the interior-point subproblem (67) can be solved using a second-order method [Nesterov, ], which provides a more reliable alternative to finite gradient differences. Under the assumption of an $L_{3}$-Lipschitz continuous third derivative of $f$, Algorithm 6 achieves convergence $\tilde{O}\left(\varepsilon^{-1/5}\right)$.

###### Theorem C.5 ( [Nesterov, 2021a , Theorem 2] )

### C.6 Optimal

An Optimal Tensor Method was recently proposed by Kovalev and Gasnikov [], Carmon et al. [], improving upon the convergence of near-optimal tensor methods [Gasnikov et al., , Bubeck et al., , Jiang et al., ]. The convergence rate was enhanced from $O\left(\varepsilon^{-2/(3p+1)}\log\left(1/\varepsilon\right)\right)$ to $O\left(\varepsilon^{-2/(3p+1)}\right)$, matching the lower bound $\Omega\left(\varepsilon^{-2/(3p+1)}\right)$ [Arjevani et al., ]. Similar to near-optimal methods, the Optimal Tensor Method is based on the A-HPE framework proposed by Monteiro and Svaiter [].

Before describing the Optimal Tensor Method, we introduce some necessary notations. Let $\Phi_{p}^{g}$ denote the $p$-th order Taylor approximation of the function $g$:

$$ $\Phi_{p}^{g}(x,y)=g(y)+\sum\limits_{k=1}^{p}\tfrac{1}{k!}D^{k}g(y)[x-y]^{k}.$ (69) $$

Additionally, note that $\Phi_{p}^{f}(x,y)=\Phi_{p}(x,y)$ as defined in (13). We also define the function $g_{\lambda}(x,y)=f(x)+\frac{1}{2\lambda}\|x-y\|^{2}$.

The main distinction from near-optimal methods lies in the procedure used to find the pair $(x_{t+1},\lambda_{t+1})$. Instead of first computing $x_{t+1}$ and then using a binary search to determine $\lambda_{t+1}$, as done in previous approaches, Kovalev and Gasnikov [] first select the parameter $\lambda_{t+1}$ and then compute $x_{t+1}$. This procedure, known as the Tensor Extragradient Method, is shown in lines 6- 10 of Algorithm 7. This method converges in a constant number of iterations, leading to the optimal convergence rate of $O\left(\varepsilon^{-2/(3p+1)}\right)$ for Algorithm 7.

Figure: Algorithm 7 Optimal Tensor Method [Kovalev and Gasnikov, , Algorithm 4]

###### Theorem C.6 ( [Kovalev and Gasnikov, 2022 , Theorem 5] )

Let $M_{p}=L_{p}$ and $\sigma=1/2$. Let

$$ $\displaystyle\nu=\left(\frac{(3p+1)^{p}C_{p}(M_{p},\sigma)R^{p-1}}{2^{p}\sqrt{ p}}\cdot\left(\frac{1+\sigma}{1-\sigma}\right)^{\frac{p-1}{2}}\right)^{-1},$ $\displaystyle\text{where}\quad C_{p}(M_{p},\sigma)=\frac{p^{p}M_{p}^{p}(1+ \sigma^{-1})}{p!(pM-L_{p})^{p/2}(pM_{p}+L_{p})^{p/2-1}}.$ $$

Then, for convex function $f$ with $L_{p}$-Lipschitz-continuous $p$-th derivative, to find $x_{T}$ such that $f(x_{T})-f^{*}\leq\epsilon$, it suffices to perform no more than $T\geq 1$ iterations of Algorithm 7, where

$$ $T=5D_{p}\cdot\left({L_{p}R^{p+1}}/{\epsilon}\right)^{\frac{2}{3p+1}}+7,$ $$

with $D_{p}$ is defined as follows:

$$ $D_{p}=\left(\frac{3^{\frac{p+1}{2}}(3p+1)^{p+1}p^{p}(p+1)}{2^{p+2}\sqrt{p}p!(p ^{2}-1)^{\frac{p}{2}}}\right)^{\frac{2}{3p+1}}.$ $$

###### Theorem C.6 ( [Kovalev and Gasnikov, 2022 , Theorem 5] )

## Appendix D Experimental details

#### Setup.

All methods and experiments were performed using Python 3.11, PyTorch 2.2.2, on a 13-inch MacBook Pro 2019 with 1,4 GHz Quad-Core Intel Core i5 and 8GB memory. All computations are done in torch.double.
All methods are implemented as PyTorch 2 optimizers.

#### Logistic Regression.

The logistic regression problem can be formulated as

$$ $f(x)=\tfrac{1}{n}\textstyle{\sum}_{i=1}^{n}\log\left(1+e^{-b_{i}\left\langle a _{i},x\right\rangle}\right)+\tfrac{\mu}{2}\|x\|_{2}^{2},\vspace{-.1cm}$ (70) $$

where $a_{i}\in\mathbb{R}^{d}$ are data features and $b_{i}\in\left\{-1;1\right\}$ are data labels for $i=1,\ldots,n$.

We present results on the a9a dataset $(d=123,n=32561)$ and w8a $(d=300,n=49749)$ from LibSVM by Chang and Lin []. We choose the starting point $x_{0}=3e$, where $e$ is a vector of all ones. This choice of $x_{0}$ allows us to show the convergence of the methods from a far point.
For Figures 1(a) and 2(a), we choose the regularizer $\mu=10^{-4}$ to get strongly-convex function $f$. For Figures 3(a),4, and 5, we choose the regularizer $\mu=0$ to get a convex function $f$.
For the better conditioning, we normalize data features $\|a_{i}\|=1$. For the normalized case, we choose theoretical $L_{2}=0.1$. We set $L_{3}=L_{2}=0.1$ to demonstrate the convergence rates for the same constants $L$. Note, that actual $L_{3}$ is smaller than $0.1$.

#### Third-order Nesterov’s lower-bound function.

The $l_{2}$-regularized third order Nesterov’s lower-bound function from Nesterov [] has the next form

$$ $f(x)=\tfrac{1}{4}\textstyle{\sum}_{i=1}^{d-1}(x_{i}-x_{i+1})^{4}-x_{1}+\tfrac{ \mu}{2}\|x\|_{2}^{2}.$ (71) $$

For Figures 1(b) and 2(b), we set $d=20$, $\mu=10^{-3}$, we’ve tuned $L_{3}=L_{2}=10.$.

#### Poisson regression.

Poisson regression is a type of generalized linear model used for analyzing count data and contingency tables. It assumes that the response variable $b_{i}$ follows a Poisson distribution, and the logarithm of its expected value can be expressed as a linear combination of unknown parameters.
The Poisson regression function has the next form

$$ $f(x)=\textstyle{\sum}_{i=1}^{n}e^{\left\langle a_{i},x\right\rangle}-b_{i} \left\langle a_{i},x\right\rangle,\vspace{-.1cm}$ (72) $$

where $a_{i}\in\mathbb{R}^{d}$ are data features and $b_{i}\in\left\{0,1,\ldots,k,\ldots\right\}$ are countable targets.

We present results for synthetic data: $d=21$, $n=6000$. We set $L_{1}=L_{2}=L_{3}=1$ and $x_{0}=e$ is all ones.

Figure: Figure 6: Comparison of different cubic and tensor accelerated methods on Poisson Regression.
Refer to caption: x11.png

The Cubic Regularized Newton (CRN) method and NATA with a tuned parameter $\nu$ demonstrate the best performance in Figure 6 (Left). Notably, CRN exhibits rapid superlinear convergence, likely due to the strong convexity properties of the loss function. Interestingly, NATA with the tuned $\nu$ manages to match CRN’s convergence rate. While Optimal Acceleration is slower than both CRN and NATA, it also achieves global superlinear convergence. In Figure 6 (Right) for $p=3$, the Tensor Nata method is the fastest, followed by the Basic Tensor Method, with the Optimal Tensor method ranking third. All three methods exhibit global superlinear convergence. The classical Nesterov Tensor Acceleration method is the slowest, likely due to its small default $\nu$. Notably, the tensor-based methods outperform their cubic counterparts.

Figure: Figure 7: Comparison of the methods by the relative value $1-\tfrac{f(x_{t+1})-f^{\ast}}{f(x_{t})-f^{\ast}}$.
Refer to caption: x13.png

The global superlinear performance of these accelerated second-order methods in Figure 7 raises the hope of establishing theoretical results on global superlinear convergence for accelerated second-order methods.

## Peer Reviews

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: How do the theoretical results proved in the paper compare to the literature? Was this setting already studied?

* Minor comments/typos:
    - p.3, l.140: $c_3$ is never introduced before.
    - p.4, l.188,196,201: "the function $f(x)$" should be " the function $f$"
    - p.5, l.263: $e$ is both the vector of all ones and $1e-4=10^{-4}$
    - p.6, l.283: idem
    - p.6, l.287: the bigger $\rightarrow$ the larger
    - p.6, l.288: repetition "first, the first..."
    - p.6, l.293: I think the paragraph is a bit too long especially the explanation on gradient descent.
    - p.6, l.321: I understand but I think it is not very clear.
    - p.7, l.331: " has next constant" is a bit odd to me
    - p.7, l.338: "second Theorem" $\rightarrow$ "second theorem"
    - p.7, l.338: the sentence is too long and contains two times "hence"
    - p.18, l.930: "subsolover"
    - p.24, l.1274: "Optimal" seems to be a typo.
soundness: 3
strengths: Overall, I think that the contributions of the paper are valuable and the technical content would be sufficient for publishing the paper in ICLR. Moreover, I did not find any major mistake in the proofs. The library seems to be qualitative and the numerical experiments are satisfactory to me.
summary: The paper focuses on the analysis of high order methods for nearly convex functions (i.e. star-convex or convex) with additional growth properties. The authors leverage a strong star-convexity and a uniform star-convexity assumption to prove the global superlinear convergence of the Cubic Regularized Newton Method and the Basic Tensor Method. In addition, they introduce an adaptive variant of Nesterov Accelerated Tensor Method called NATA (for Nesterov Accelerated Tensor Method with $A_t$-Adaptation) which solves the problem of having too conservative parameters. Theoretical convergence guarantees are given as well as numerical experiments highlighting good performance. The authors also provide OPTAMI, a python library for high order methods.
weaknesses: I am not convinced by the structure of the paper and the way it is written. It seems to me that the authors try to answer too many questions for a 10 pages paper. Also, I think that the theoretical claims are not discussed enough. Detailed comments can be found below. The main problem I have with the current version of the paper is that it lacks a clear unified story and it seems  an aggregation of results. In my opinion, this could be improved by modifying the structure of the paper.

The OPTAMI library is only mentioned but never properly introduced. I think that at least a paragraph should be dedicated to it or else, that it should be removed from the title.
About the structure:

   1) I believe that the introduction is already too technical although it is well explained. It also states too many problems and, due to a lack of structure, it can confuse the reader. This is highlighted by the fact that there are three questions at the end of the introduction. My suggestion is to keep the detailed discussions for later sections and do a way shorter and high-level introduction. Also, I would try to use as few equations as possible.

    2) It is related to the previous comment but I think that the second section "Basic methods" could be enriched with some comments from the introduction.

    3) I think that the third section is interesting with intuitive proofs. However, I regret that there is no discussion on the stated theorems: is there any similar result in the literature, was it expected, is it tight?

    4) The fourth section seems to come bit out of nowhere and lacks a proper introduction. As a reader, it can seem odd to come from superlinear convegence rates (which is in the name of the paper) to adaptive techniques and a new method without a paragraph that bridges both sections. As said before, Theorem 4.1 should be commented and I would have expected a comparison with the results ensured by vanilla NATM.

    5) I believe that the fifth section should be merged with the fourth one or defined as a subsection.

---

confidence: 4
contribution: 3
initial_rating: 8
presentation: 4
questions: 1. How do the authors interpret this trade-off between accuracy and rate of convergence? If the requirement on accuracy is relaxed and the conditions are different, how can the superlinear convergence rate be seen as an improvement compared to previous results?
2. In addition to the first question, Song et. al. 2021 proposed an acceleration framework that matches the lower bound established by Arjevani et al. 2019. Shouldn't that be seen as the optimal rate for this setting of optimizing strongly convex functions with second-order methods?

Chaobing Song, Yong Jiang, and Yi Ma. Unified acceleration of high-order algorithms under general holder
continuity. SIAM Journal on Optimization, 31(3):1797–1826, 2021.
soundness: 4
strengths: 1. The proposed acceleration variant of tensor method achieves better empirical performance than the existing (near)-optimal accelerated second-order methods.
2. All methods are systematically implemented and released as a library.
summary: This paper first establishes the global superlinear convergence rate for second-order methods in optimizing strongly star-convex functions, and its higher-order extension for optimizing uniformly star-convex functions, then proposes a variant of Nesterov accelerated tensor method which demonstrates superior performance in numerical experiments consistent with the theoretical faster convergence rate, and finally presents a systematical numerical comparison across mainstream second-order methods.
weaknesses: 1. The global linear rate is established by relaxing the required accuracy to exceed the radius of the quadratic convergence region.
2. Some typos: 
- line 143 "$\epsilon \leq c_3 r$" --> $\epsilon > c_3 r$ 
- Eq (20) $t \rightarrow 0$ --> $t \rightarrow \infty$?

---

confidence: 3
contribution: 2
initial_rating: 6
presentation: 3
questions: I have a few other minor comments for the authors to consider:

- Line 143: I think it should be $\epsilon > c_3 r$?

- Line 295: "where the sublinear rate outperforms the linear rate" is a little confusing. Maybe the authors meant the high-order methods have not entered the linear convergence region and the convergence rate is sublinear at the beginning?

- For the CRN method, did the authors prove the superlinear convergence for $\mu_q$-uniformly star-convex functions? If so, it may be better to state the results in Theorem 3.2. Currently, I cannot find the results for $\mu_q$-uniformly star-convex functions and the CRN method.
soundness: 3
strengths: The paper is well-organized and easy to follow. The results are novel and should be interesting to the audience from machine learning and optimization fields. The problems studied in this work are important and applicable to certain practical problems where the computational time is not a critical constraint. The proof in the main manuscript should be correct, while I do not have time to check the proofs in the appendix due to the time limit.
summary: In this work, the authors studied the convergence rate of high-order methods, e.g., the cubic regularized Newton method and the basic tensor method. The authors proved the global superlinear convergence of both methods for $\mu$-strongly star-convex functions and $\mu_q$-uniformly star-convex functions, respectively. In addition, a variant of accelerated high-order method, named NATA, was proposed and compared with other accelerated high-order methods.
weaknesses: The main problem with the paper is that the first part (superlinear convergence rate of high-order methods) and the second part (NATA algorithm) seem to be independent and can be separated into two papers. I feel that these two parts considered two different topics. The first one is mostly theoretical and is about non-accelerated methods, while the second one is about accelerated methods and their empirical performance. I would suggest the authors split the paper into two and include more details to the content. For example, the intuition behind the design of NATA.

In addition, I think the current introduction section is too lengthy. Considering the page limit of the conference, the background knowledge can be simplified and moved to the appendix, since it can be easily found in textbooks and literature. 

Finally, I think the authors could include more details to the experiments. For example, in Figure 3, I wonder if the Tensor NATA converges faster with a carefully chosen $v_t$? This is not discussed. Also, it would be better if the authors could provide more intuition behind the current design of searching $v_t$ in sub-iterations instead of fixing $v_t$ to be a constant. Since the performance of a fixed $v_t$ is better than that of a adaptive $v_t$, I wonder if this design is unnecessary. Furthermore, I think it will be helpful if the authors could provide the running time comparison. This is because the solving time of the sub-problem in each iteration may be non-negligible and comparable to the computation time of the Hessian matrix. Especially, for the NATA method, the sub-problem needs to be solved several times in each iteration.

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 2
questions: 1) What are the total computational complexities for the Cubic Regularized Newton Method and Basic Tensor Method if the costs of solving the subproblems in each iteration can be evaluated? How attractive are they compared to the total computational complexities of other algorithms? 

2) If the subproblems of the Cubic Regularized Newton Method and Basic Tensor Method are to be solved iteratively, can't we allow them to be solved inexactly and include the errors in total computational costs?

3) About the proposed algorithm, Newterov Accelerated Tensor Method with At-Adaptation (NATA). Can the authors show anything about the iteration complexities for the strongly (star-)convex case? If it is difficult, show us the best iteration complexity for the classical Nesterov accelerated tensor method when the function $f$ is strongly convex, and tell us why it is difficult to derive the complexity for NATA. 

4) For the strongly (star-)convex case, how does NATA perform in numerical experiments? More precisely, what would Figure 3(b) look like in the case of strongly convex (21) with positive $\mu$ (perhaps the horizontal axis might need to be the iteration number)?
soundness: 3
strengths: It is great to show that cubic Newton and Basic Tensor methods can achieve superlinear convergence for strongly convex functions.
This paper also includes the necessary equations and propositions for following the proof of superlinear convergence in the main text, and explains the proof strategy leading to the most important theorem, Theorem 3.4, in an easy-to-understand manner. 

 I also think that providing an open-source computational library for high-order methods would be very helpful to subsequent researchers.
summary: The paper presents two main results: the first shows that the basic high-order methods (the Cubic Regularized Newton Method and Basic Tensor Method) achieve global superlinear convergence for μ-strongly star-convex functions. 
Second, it proposes a technique for choosing the parameters of an existing method (Nesterov Accelerated Tensor Method), evaluates its convergence speed, and shows its usefulness in numerical experiments. 
Furthermore, the paper introduces an open-source computational library called OPTAMI for high-order methods.
weaknesses: [Structure of the paper]
Section 3 and Section 4 are not well connected, and the paper seems to contain two disparate contents; it would be better to divide the paper into two papers, and strengthen the contents of each, for example, in the following ways. The algorithms and classes of functions are different between Sections 3 and 4.

For example, for the contents of Section 3, if the cost of solving the subproblem in each iteration can be evaluated, what is the total computational complexity? How attractive is it compared to the total computational complexities of other algorithms? Even if the number of iterations (i.e., the worst-case iteration complexity) is reduced, on the other hand, if the computational complexity at each iteration explodes, it will not be attractive as an algorithm, so it would be better to provide the total computational complexity. If the subproblem is to be solved iteratively, can't we allow it to be solved inexactly and include the error in the iteration complexity and total computational cost?

For example, for the contents of Section 4, if we consider strongly star-convex, how can the theoretical guarantee given in Theorem 4.1 be changed? As future work, the authors wrote this type of question, but what trends do you see, at least in numerical experiments? More precisely, what would Figure 3(b) look like in the case of strongly convex (21) with positive $\mu$ (perhaps the horizontal axis might need to be the iteration number)?

----------

[Various unclear descriptions]
There are unclear descriptions and various minor errors, giving the impression that this paper was written in haste.

For example, lines 62-74 introduce the existing studies and describe their global convergence rates, but the assumption of the function $f$ is not clearly stated. Before that, there is a definition of star-convex, etc., but I do not think that those existing studies assume star-convex. This is because the Hessian matrix in (7) is not necessarily positive definite when the function is nonconvex. 
The others are line 80: upperbound --> upper bound, NPE in line 109 should be written as Newton Proximal Extragradient (NPE), in line 123 the methods performs --> the method performs, etc. Line 193 defines the norm $||x||$ using a matrix B, but which kind of matrices are specifically used in your claims? If we always assume $B=I$, why use $B$ to define the norm? I won't point out any more, but I would like the authors to review the paper again and and reduce this type of typo.
 
----------

[Insufficient reference]
The survey of existing studies does not seem to be sufficient. For example, in line 97, the authors write “poor global convergence” for the quasi-Newton method, but I do not think they are aware of the following paper. The paper also shows “Global Non-Asymptotic Superlinear Convergence” as well as this paper.

Qiujiang Jin, Ruichen Jiang, Aryan Mokhtari, "Non-asymptotic Global Convergence Rates of BFGS with Exact Line Search",
arXiv:2404.01267, 2024.

 Although the authors mention the Jiang, Jin, and Mokhtari paper (COLT 2023) in Appendix A, I think the global convergence rates of quasi-Newton methods should be mentioned around lines 97-98, including the above-mentioned papers. Can it no longer be called “poor global convergence”?

Meta-review: The paper focuses on the analysis of high order methods for mu-strongly star-convex functions. They propose an adaptive variant of Nesterov Accelerated Tensor Method called NATA. The authors also provide OPTAMI, a python library for high order methods. The reviewers found the paper to be solid, though with some concerns that the presentation of the paper is trying to cover too much ground. Based on my own reading, the contribution of the paper is solid and worthy of acceptance.

Decision: Accept (Poster)