# Advancing the lower bounds:
An accelerated, stochastic, second-order method with optimal adaptation to
inexactness

Artem Agafonov  Moscow Institute of Physics and Technology, Dolgoprudny,
Russia  Mohamed bin Zayed University of Artificial Intelligence, Abu Dhabi,
UAE  Dmitry Kamzolov  Mohamed bin Zayed University of Artificial Intelligence,
Abu Dhabi, UAE  Alexander Gasnikov  Moscow Institute of Physics and
Technology, Dolgoprudny, Russia  Institute for Information Transmission
Problems, Moscow, Russia  HSE University, Moscow, Russia  Ali Kavis
Laboratory for Information and Inference Systems, IEM STI EPFL, Lausanne,
Switzerland  Kimon Antonakopoulos  Laboratory for Information and Inference
Systems, IEM STI EPFL, Lausanne, Switzerland  Volkan Cevher  Laboratory for
Information and Inference Systems, IEM STI EPFL, Lausanne, Switzerland  Martin
Takáč  Mohamed bin Zayed University of Artificial Intelligence, Abu Dhabi, UAE

###### Abstract

We present a new accelerated stochastic second-order method that is robust to
both gradient and Hessian inexactness, which occurs typically in machine
learning. We establish theoretical lower bounds and prove that our algorithm
achieves optimal convergence in both gradient and Hessian inexactness in this
key setting. We further introduce a tensor generalization for stochastic
higher-order derivatives. When the oracles are non-stochastic, the proposed
tensor algorithm matches the global convergence of Nesterov Accelerated Tensor
method. Both algorithms allow for approximate solutions of their auxiliary
subproblems with verifiable conditions on the accuracy of the solution.

##  1 Introduction

In this paper, we consider the following general convex optimization problem:

$\min\limits_{x\in\mathbb{R}^{d}}f(x),$ (1)

where $f(x)$ is a convex and sufficiently smooth function. We assume that a
solution $x^{\ast}\in\mathbb{R}^{d}$ exists and we denote
$f^{\ast}:=f(x^{\ast})$.

######  Assumption 1.1.

The function $f(x)\in C^{2}$ has $L_{2}$-Lipschitz-continuous Hessian if for
any $x,y\in\mathbb{R}^{d}$

$\|\nabla^{2}f(x)-\nabla^{2}f(y)\|\leq L_{2}\|x-y\|.$

Since the calculation of an exact gradient is very expensive or impossible in
many applications in several domains, including machine learning, statistics,
and signal processing, efficient methods that can work with inexact stochastic
gradients are of great interest.

######  Assumption 1.2.

For all $x\in\mathbb{R}^{d}$, we assume that stochastic gradients
$g(x,\xi)\in\mathbb{R}^{d}$ satisfy

$\mathbb{E}[g(x,\xi)\mid x]=\nabla f(x),\quad\mathbb{E}\left[\|g(x,\xi)-\nabla
f(x)\|^{2}\mid x\right]\leq\sigma_{1}^{2}.$ (2)

Extensive research has been conducted on first-order methods, both from a
theoretical and practical perspective. For $L_{1}$-smooth functions with
stochastic gradients characterized by a variance of $\sigma_{1}^{2}$, lower
bound
$\Omega\left(\tfrac{\sigma_{1}R}{\sqrt{T}}+\tfrac{L_{1}R^{2}}{T^{2}}\right)$
has been established [47]. For $f(x)=\mathbb{E}[F(x,\xi)]$, the stochastic
approximation (SA) was developed, starting from the pioneering paper by
Robbins and Monro (1951) [61]. Important improvements of the SA method were
developed by Polyak [56], Polyak and Juditsky [57], and Nemirovski et al.
[46], where longer stepsizes with iterate averaging and proper step-size
modifications were proposed, obtaining the rate
$O\left(\tfrac{\sigma_{1}R}{\sqrt{T}}+\tfrac{L_{1}R^{2}}{T}\right)$. The
optimal method matching the lower bounds have been developed by Lan [41] with
a convergence rate
$O\left(\tfrac{\sigma_{1}R}{\sqrt{T}}+\tfrac{L_{1}R^{2}}{T^{2}}\right)$.
However, the literature on second-order methods is significantly limited for
the study of provable, globally-convergent stochastic second-order methods for
convex minimization.
Second-order methods. Although second-order methods have been studied for
centuries [55, 60, 64, 38, 45, 29], the most of the results are connected with
local quadratic convergence. The significant breakthroughs regarding the
global convergence have been achieved only recently starting from the paper on
Cubic Regularized Newton (CRN) method by Nesterov and Polyak [54], the first
second-order method with a global convergence rate
$O\left(\tfrac{L_{2}R^{3}}{T^{2}}\right)$. Following this work, Nesterov [49]
proposes an acceleration mechanism on top of CRN and achieve the convergence
rate of $O\left(\tfrac{L_{2}R^{3}}{T^{3}}\right)$, going beyond the
$O(1/T^{2})$ lower bound for first-order methods. Another cornerstone in this
field is the work by Monteiro and Svaiter [44], which achieves lower
complexity bound $\Omega\left(\tfrac{L_{2}R^{3}}{T^{7/2}}\right)$ [3, 5] up to
a logarithmic factor, for the first time in the literature. The gap between
upper and lower bounds was closed only in 2022 in subsequent works of Kovalev
and Gasnikov [40]; Carmon et al. [10].
One of the main limitations of the second-order methods is a high per-
iteration cost as they require computation of exact Hessian. Therefore, it is
natural to use approximations of derivatives instead of their exact values. In
[25], CRN method with $\delta_{2}$-inexact Hessian information and its
accelerated version were proposed, achieving convergence rate
$O\left(\frac{\delta_{2}R^{2}}{T^{2}}+\frac{L_{2}R^{3}}{T^{3}}\right)$. This
algorithm was later extended in [2] to handle $\delta_{1}$-inexact gradients
(and high-order derivatives) with a resulting convergence rate of
$O\left(\delta_{1}R+\frac{\delta_{2}R^{2}}{T^{2}}+\frac{L_{2}R^{3}}{T^{3}}\right)$.
A recent paper [4] proposes a stochastic adaptive second-order method based on
the extragradient method without line-search and with the convergence rate
$O\left(\frac{\sigma_{1}R}{\sqrt{T}}+\frac{\sigma_{2}R^{2}}{T^{3/2}}+\frac{L_{2}R^{3}}{T^{3}}\right)$
when gradients and Hessians are noisy with variances $\sigma_{1}^{2}$ and
$\sigma_{2}^{2}$.

In the light of these results, we identify several shortcomings and open
questions:

What are the lower bounds for inexact second-order methods?
What is the optimal trade-off between inexactness in the gradients and the
Hessian?

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

Table 1:  Comparison of existing results for second-order methods under
inexact feedback. $T$ denotes the number of iterations, and $L_{2}$ represents
the Lipschitz constant of the Hessian.

Algorithm | Inexactness |  Gradient convergence |  Hessian convergence |  Exact convergence
---|---|---|---|---
Accelerated Inexact Cubic Newton [25] |  exact gradient $\delta_{2}$-inexact Hessian 111 $\delta_{2}$-inexact Hessian: $-\delta_{2}I\preceq\nabla^{2}f(x)-B_{x}\preceq-\delta_{2}I/2$ | ✗ | $O\left(\frac{\delta_{2}R^{2}}{T^{2}}\right)$ | $O\left(\frac{L_{2}R^{3}}{T^{3}}\right)$
Accelerated Inexact Tensor Method 222 It is worth noting that the Accelerated Inexact Tensor Method can also be applied to the case of stochastic derivatives. Specifically, when $p=2$, the total number of stochastic gradient computations is on the order of $O(\varepsilon^{-7/3})$, while the total number of stochastic Hessian computations is on the order of $O(\varepsilon^{-2/3})$ [2]. In our work, we propose an algorithm that achieves the same number of stochastic Hessian computations, but significantly improves the number of stochastic gradient computations to $O(\varepsilon^{-2})$.  [2] |  $\delta_{1}$-inexact gradient 333 $\delta_{1}$-inexact gradient: $\|g_{x}-\nabla f(x)\|\leq\delta_{1}$ $\delta_{2}$-inexact Hessian 444 $\delta_{2}$-inexact Hessian: $\|\left(H_{x}-\nabla^{2}f(x)\right)(y-x)\|\leq\delta_{2}\|y-x\|$ | $O(\delta_{1}R)$ | $O\left(\frac{\delta_{2}R^{2}}{T^{2}}\right)$ | $O\left(\frac{L_{2}R^{3}}{T^{3}}\right)$
Extra-Newton [4] |  stochastic gradient (2) unbiased stochastic Hessian 555 Unbiased stochasic Hessian: $\mathbb{E}[H(x,\xi)\mid x]=\nabla^{2}f(x)$, $\mathbb{E}\left[\left\|H(x,\xi)-\nabla^{2}f(x)\right\|^{2}\mid x\right]\leq\sigma_{2}^{2}$ | $O\left(\frac{\sigma_{1}R}{\sqrt{T}}\right)$ | $O\left(\frac{\sigma_{2}R^{2}}{T^{3/2}}\right)$ | $O\left(\frac{L_{2}R^{3}}{T^{3}}\right)$
Accelerated Stochastic Second-order method [This Paper] |  stochastic gradient (2) stochastic Hessian (4) | $O\left(\frac{\sigma_{1}R}{\sqrt{T}}\right)$ |  $O\left(\frac{\sigma_{2}R^{2}}{T^{2}}\right)$ 666 Under assumption of $\delta_{2}$-inexact Hessian the convergence is $O\left(\frac{\delta_{2}R^{2}}{T^{2}}\right)$ | $O\left(\frac{L_{2}R^{3}}{T^{3}}\right)$
Lower bound [This Paper] |  stochastic gradient (2) stochastic Hessian (4) | $\Omega\left(\frac{\sigma_{1}R}{\sqrt{T}}\right)$ |  $\Omega\left(\frac{\sigma_{2}R^{2}}{T^{2}}\right)$ |  $\Omega\left(\frac{L_{2}R^{3}}{T^{7/2}}\right)$

Contributions. We summarize our contributions as follows:

  1. 1.

We propose an accelerated second-order algorithm that achieves the convergence
rate of
$O\left(\frac{\sigma_{1}R}{\sqrt{T}}+\frac{\sigma_{2}R^{2}}{T^{2}}+\frac{L_{2}R^{3}}{T^{3}}\right)$
for stochastic Hessian with variance $\sigma_{2}^{2}$ and
$O\left(\frac{\sigma_{1}R}{\sqrt{T}}+\frac{\delta_{2}R^{2}}{T^{2}}+\frac{L_{2}R^{3}}{T^{3}}\right)$
for $\delta_{2}$-inexact Hessian, improving the existing results [2, 4] (see
Table 1).

  2. 2.

We prove that the above bounds are tight with respect to the variance of the
gradient and the Hessian by developing a matching lower bound (see Table 1).

  3. 3.

Our algorithm involves solving a cubic subproblem that arises in several
globally convergent second-order methods [54, 49]. To address this, we propose
a criterion based on the accuracy of the subproblem’s gradient, along with a
dynamic strategy for selecting the appropriate level of inexactness. This
ensures efficient solution of the subproblems without sacrificing the fast
convergence of the initial method.

  4. 4.

We extend our method for higher-order minimization with stochastic/inexact
oracles. We achieve the
$O\left(\frac{\sigma_{1}R}{\sqrt{T}}+\sum\limits_{i=2}^{p}\frac{\delta_{i}R^{i}}{T^{i}}+\frac{L_{p}R^{p+1}}{T^{p+1}}\right)$
rate with $\delta_{i}$-inexact $i$-th derivative.

  5. 5.

We propose a restarted version of our algorithm for strongly-convex
minimization, which exhibits a linear rate. Via a mini-batch strategy, we
demonstrate that the total number of Hessian computations scales linearly with
the desired accuracy $\varepsilon$.

Related work. The idea of using high-order derivatives in optimization has
been known for a long time [31]. In 2009, M. Baes extended the cubic
regularization approach with second-order derivatives ($p=2$) from [54] to
high-order derivatives ($p>2$) in [6]. However, the subproblem of these
methods was non-convex, making them impractical. In 2018, Yu. Nesterov
proposed the implementable (Accelerated) Tensor Method [51], wherein the
convexity of the subproblem was reached by increasing a regularization
parameter. Hence, the convex subproblem could be efficiently solved by
appropriate subsolvers, making the algorithm practically applicable. In the
same work, a lower complexity bound for tensor methods under higher-order
smoothness assumption was proposed. Shortly after, near-optimal [23, 22, 9,
32] and optimal [40, 10] high-order methods were introduced. Furthermore,
under higher smoothness assumptions, _second-order_ methods [53, 52, 35, 18]
can surpass the corresponding lower complexity bound for functions with
Lipschitz-continuous Hessians. For more comprehensive information on high-
order methods, one can refer to the review [36].
In general, second and higher-order methods are known for their faster
convergence compared to first-order methods. However, their computational cost
per iteration is significantly higher due to the computation of high-order
derivatives. To alleviate this computational burden, it is common to employ
approximations of derivatives instead of exact values. While there is a wide
range of second-order and tensor methods available for the non-convex case,
assuming stochastic or inexact derivatives [12, 13, 14, 39, 67, 65, 42, 7, 8,
17], the same cannot be said for the convex case. In the context of convex
problems, there have been studies on high-order methods such as second-order
methods with inexact Hessian information [25], tensor methods with inexact and
stochastic derivatives [2], and Extra-Newton algorithm with stochastic
gradients and Hessians [4]. As a possible application of methods with inexact
Hessians, we highlight Quasi-Newton(QN) methods. Such methods approximate
second-order derivatives using the history of gradient feedback. Quasi-Newton
methods are known for their impressive practical performance and local
superlinear convergence. However, for the long period of time, the main
drawback of such methods was a slow theoretical global convergence, slower
than gradient descent. First steps to improve the global convergence of such
methods were done in [62, 26] but the methods could be still slower than
gradient descent. The first global Quasi-Newton methods that provably matches
the gradient descent were reached by cubic regularization in two consecutive
papers [37, 33]. It also opened a possibility for accelerated QN [37] that
theoretically matches fast gradient method and first-order lower-bounds. This
direction were further explored for different methods in [63, 34]. Another
possible application for high-order methods with inexact or stochastic
derivatives is distributed optimization [68, 15, 1, 21].
One of the main challenges of tensor and regularized second-order methods is
solving the auxiliary subproblem to compute the iterate update. In both the
second-order and higher-order cases, it usually requires running a subsolver
algorithm. The impact of the accuracy, up to which we solve the auxiliary
problem, on the convergence of the algorithm has been studied in several works
[27, 28, 20]. One actively developing direction relies on the constructions of
CRN with explicit step [59, 58, 43, 19, 18, 30].
Organization. The rest of the paper is organized as follows. Section 2
introduces the main notations and assumptions. In Section 3, we present the
Accelerated Stochastic Cubic Regularized Newton method and provide its
convergence analysis. The lower bound is discussed in Section 4. The tensor
generalization of the proposed algorithm is presented in Section 5. Lastly,
Section 6 introduces a restarted version of the proposed method.

##  2 Problem statement and preliminaries

Taylor approximation and oracle feedback. Our starting point for constructing
second-order method is based primarily on the second-order Taylor
approximation of the function $f(x)$

$\Phi_{x}(y)\stackrel{{\scriptstyle\text{ def }}}{{=}}f(x)+\left\langle\nabla
f(x),y-x\right\rangle+\tfrac{1}{2}\left\langle
y-x,\nabla^{2}f(x)(y-x)\right\rangle,\quad y\in\mathbb{R}^{d}.$

In particular, since the exact computation of the Hessians can be a quite
tiresome task, we attempt to employ more tractable inexact estimators $g(x)$
and $H(x)$ for the gradient and Hessian. These estimators are going to be the
main building blocks for the construction of the "inexact" second-order Taylor
approximation. Formally, this is given by:

$\phi_{x}(y)=f\left(x\right)+\left\langle
g(x),y-x\right\rangle+\tfrac{1}{2}\left\langle
y-x,H(x)(y-x)\right\rangle,\quad y\in\mathbb{R}^{d}.$ (3)

Therefore, by combining Assumption 1.1 with the aforementioned estimators, we
readily get the following estimation:

######  Lemma 2.1.

Let Assumption 1.1 hold. Then, for any $x,y\in\mathbb{R}^{d}$, we have

$\begin{split}|f(y)-\phi_{{x}}(y)|\leq\left(\|g(x)-\nabla
f(x)\|+\tfrac{1}{2}\|\left(H(x)-\nabla^{2}f(x)\right)(y-x)\|\right)\|y-x\|+\tfrac{L_{2}}{6}\|y-x\|^{3}.\end{split}$

$\begin{split}\|\nabla f(y)-\nabla\phi_{x}(y)\|\leq\|g(x)-\nabla
f(x)\|+\|\left(H(x)-\nabla^{2}f(x)\right)(y-x)\|+\tfrac{L_{2}}{3}\|y-x\|^{2}.\end{split}$

Now, having established the main toolkit concerning the approximation of $f$
in the rest of this section, we introduce the blanket assumptions regarding
the inexact gradients and Hessians (for a complete overview, we refer to Table
1). In particular, we assume that our estimators satisfy the following
statistical conditions.

######  Assumption 2.2 (Unbiased stochastic gradient with bounded variance and
bounded variance stochastic Hessian).

For all $x\in\mathbb{R}^{d}$, stochastic gradient $g(x,\xi)$ satisfies (2) and
stochastic Hessian $H(x,\xi)$ satisfies

$\mathbb{E}\left[\|H(x,\xi)-\nabla^{2}f(x)\|_{2}^{2}\mid
x\right]\leq\sigma_{2}^{2}.$ (4)

######  Assumption 2.3 (Unbiased stochastic gradient with bounded variance and
inexact Hessian).

For all $x\in\mathbb{R}^{d}$ stochastic gradient $g(x,\xi)$ satisfies (2). For
given $x,~{}y\in\mathbb{R}^{d}$ inexact Hessian $H(x)$ satisfies

$\displaystyle\|(H(x)-\nabla^{2}f({x}))[y-x]\|\leq\delta_{2}^{x,y}\|y-x\|.$
(5)

Assumptions 2.2 and 2.3 differ from Condition 1 in [2] by the unbiasedness of
the gradient. An unbiased gradient allows us to attain optimal convergence in
the corresponding term $O(1/\sqrt{T})$, while an inexact gradient slows down
the convergence to $O(1)$ since a constant error can misalign the gradient.
Note, that we do not assume the unbiasedness of the Hessian in all
assumptions. Finally, note that Assumption 2.3 does not require (5) to be met
for all $x,~{}y\in\mathbb{R}^{d}$. Instead, we only consider inexactness along
the direction $y-x$, which may be significantly less than the norm of the
difference between Hessian and its approximation $H(x)$.
Auxiliary problem. Most second-order methods with global convergence require
solving an auxiliary subproblem at each iteration. However, to the best of our
knowledge, existing works that consider second-order methods under inexact
derivatives do not account for inexactness in the solution of the subproblem.
To address this gap, we propose incorporating a gradient criteria for the
subproblem solution, given by

$\displaystyle\textstyle{\min}_{y\in\mathbb{R}^{d}}\omega_{x}(y)~{}~{}\text{such
that}~{}~{}\|\nabla\omega_{x}(y)\|\leq\tau,$ (6)

where $\omega_{x}(y)$ is the objective of subproblem and $\tau\geq 0$ is a
tolerance parameter. We highlight that this criterion is verifiable at each
step of the algorithm, which facilitates determining when to stop. By setting
a constant tolerance parameter $\tau$, we get the following relationship
between the absolute accuracy $\epsilon$ required for the initial problem and
$\tau$: $\tau=O\left(\epsilon^{\frac{5}{6}}\right)$. In practice, it may not
be necessary to use a very small accuracy in the beginning. Later, we will
discuss strategies for choosing the sequence of $\delta_{t}$ based on the
number of iterations $t$.

##  3 The method

In this section, we present our proposed method, dubbed as Accelerated
Stochastic Cubic Regularized Newton’s method. In particular, extending on
recent accelerated second-order algorithms [51, 25, 2], we propose a new
variant of the accelerated cubic regularization method with stochastic
gradients that achieves optimal convergence in terms corresponding to gradient
and Hessian inexactness. Moreover, the proposed scheme allows for the
approximate solution of the auxiliary subproblem, enabling a precise
determination of the required level of subproblem accuracy.
We begin the algorithm description by introducing the main step. Given
constants $\bar{\delta}>0$ and $M\geq L_{2}$, we define a model of the
objective

$\omega^{M,\bar{\delta}}_{x}(y):=\phi_{x}(y)+\tfrac{\bar{\delta}}{2}\|x-y\|^{2}+\tfrac{M}{6}\|x-y\|^{3}.$

At each step of the algorithm, we aim to find
$u\in\operatornamewithlimits{argmin}_{y\in\mathbb{R}^{d}}\omega_{x}^{M,\bar{\delta}}(y)$.
However, finding the respective minimizer is a separate challenge. Instead of
computing the exact minimum, we aim to find a point $s\in\mathbb{R}^{d}$ with
a small norm of the gradient.

######  Definition 3.1.

Denote by $s^{M,\bar{\delta},\tau}(x)$ a $\tau$-inexact solution of
subproblem, i.e. a point $s:=s^{M,\bar{\delta},\tau}(x)$ such that

$\|\nabla\omega^{M,\bar{\delta}}_{x}(s)\|\leq\tau.$

Next, we employ the technique of estimating sequences to propose the
Accelerated Stochastic Cubic Newton method. Such acceleration is based on
aggregating stochastic linear models given by

$l(x,y)=f(y)+\left\langle g(x,\xi),x-y\right\rangle$

in function $\psi_{t}(x)$. The method is presented in detail in Algorithm 1.

Throughout the paper we define $R=\|x_{0}-x^{*}\|.$

Algorithm 1 Accelerated Stochastic Cubic Newton

1: Input: $y_{0}=x_{0}$ is starting point; constants $M\geq 2L_{2}$; non-
negative non-decreasing sequences $\\{\bar{\delta}_{t}\\}_{t\geq
0},\\{\lambda_{t}\\}_{t\geq 0},\\{\kappa_{2}^{t}\\}_{t\geq 0}$,
$\\{\kappa_{3}^{t}\\}_{t\geq 0}$, and
$\alpha_{t}=\tfrac{3}{t+3},~{}~{}~{}A_{t}=\textstyle{\prod\limits_{j=1}^{t}}(1-\alpha_{j}),~{}~{}~{}A_{0}=1,$
$\psi_{0}(x):=\tfrac{\bar{\kappa}_{2}^{0}+\lambda_{0}}{2}\|x-x_{0}\|^{2}+\tfrac{\bar{\kappa}_{3}^{0}}{3}\|x-x_{0}\|^{3}.$

2: for $t\geq 0$ do

3: $v_{t}=(1-\alpha_{t})x_{t}+\alpha_{t}y_{t},\quad
x_{t+1}=s^{M,\bar{\delta}_{t},\tau}(v_{t})$

4: Compute  $\displaystyle y_{t+1}=\arg\min_{x\in\mathbb{R}^{n}}$
$\displaystyle\left\\{\psi_{t+1}(x):=\psi_{t}(x)+\tfrac{\lambda_{t+1}-\lambda_{t}}{2}\|x-x_{0}\|^{2}\right.$
$\displaystyle\left.+\textstyle{\sum\limits_{i=2}^{3}}\tfrac{\bar{\kappa}^{t+1}_{i}-\bar{\kappa}^{t}_{i}}{i}\|x-x_{0}\|^{i}+\tfrac{\alpha_{t}}{A_{t}}l(x,x_{t+1})\right\\}.$

######  Theorem 3.2.

Let Assumption 1.1 hold and $M\geq 2L_{2}$.

  * •

Let Assumption 2.2 hold. After $T\geq 1$ with parameters

$\displaystyle\bar{\kappa}_{2}^{t+1}=\tfrac{2\bar{\delta}_{t}\alpha_{t}^{2}}{A_{t}},~{}\bar{\kappa}_{3}^{t+1}=\frac{8M}{3}\tfrac{\alpha_{t+1}^{3}}{A_{t+1}},~{}\lambda_{t}=\tfrac{\sigma_{1}}{R}(t+3)^{\frac{5}{2}},~{}\bar{\delta}_{t}=2\sigma_{2}+\tfrac{\sigma_{1}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}{R}(t+3)^{\frac{3}{2}},$
(7)

we get the following bound

$\displaystyle\mathbb{E}\left[f(x_{T})-f(x^{\ast})\right]$ $\displaystyle\leq
O\left(\tfrac{\tau
R}{\sqrt{T}}+\tfrac{\sigma_{1}R}{\sqrt{T}}+\tfrac{\sigma_{2}R^{2}}{T^{2}}+\tfrac{MR^{3}}{T^{3}}\right).$
(8)

  * •

Let Assumption 2.3 hold. After $T\geq 1$ with parameters defined in (7) and
$\sigma_{2}=\delta_{2}=\max\limits_{t=1,\ldots,T}\delta_{t}^{v_{t-1},x_{t}}$,
we get the following bound

$\displaystyle\mathbb{E}[f(x_{T})-f(x^{\ast})]$ $\displaystyle\leq
O\left(\tfrac{\tau
R}{\sqrt{T}}+\tfrac{\sigma_{1}R}{\sqrt{T}}+\tfrac{\delta_{2}R^{2}}{T^{2}}+\tfrac{MR^{3}}{T^{3}}\right).$
(9)

This result provides an upper bound for the objective residual after $T$
iterations of Algorithm 1. The last term in the RHS of (8) and (9) corresponds
to the case of exact Accelerated Cubic Newton method [49]. The remaining terms
reveal how the convergence rate is affected by the imprecise calculation of
each derivative and by inexact solution of subproblem. We provide sufficient
conditions for the inexactness in the derivatives to ensure that the method
can still obtain an objective residual smaller than $\varepsilon$.
Specifically, this result addresses the following question: given that the
errors are controllable and can be made arbitrarily small, how small should
each derivative’s error be to achieve an $\varepsilon$-solution?

######  Corollary 3.3.

Let assumptions of Theorem 3.2 hold and let $\varepsilon>0$ be the desired
solution accuracy.

  * •

Let the levels of inexactness in Assumption 2.2 be:

$\displaystyle\tau=O\left(\varepsilon^{\frac{5}{6}}\left(\tfrac{M}{R^{3}}\right)^{\frac{1}{6}}\right),\quad\sigma_{1}=O\left(\varepsilon^{\frac{5}{6}}\left(\tfrac{M}{R^{3}}\right)^{\frac{1}{6}}\right),\quad\sigma_{2}=O\left(\varepsilon^{\frac{1}{3}}M^{\frac{2}{3}}\right)$

  * •

Let the levels of inexactness in Assumption 2.3 be:

$\displaystyle\tau=O\left(\varepsilon^{\frac{5}{6}}\left(\tfrac{M}{R^{3}}\right)^{\frac{1}{6}}\right),\quad\sigma_{1}=O\left(\varepsilon^{\frac{5}{6}}\left(\tfrac{M}{R^{3}}\right)^{\frac{1}{6}}\right),\quad\delta_{2}=O\left(\varepsilon^{\frac{1}{3}}M^{\frac{2}{3}}\right)$

And let the number of iterations of Algorithm 1 satisfy
$T=O\left(\frac{MR^{3}}{\varepsilon}\right)^{\frac{1}{3}}$. Then $x_{T}$ is an
$\varepsilon$-solution of problem (1), i.e.
$f(x_{T})-f(x^{\ast})\leq\varepsilon$.

In practice, achieving an excessively accurate solution for the subproblem on
the initial iterations is not essential. Instead, a dynamic strategy can be
employed to determine the level of accuracy required for the subproblem.
Specifically, we can choose a dynamic precision level according to
$\tau_{t}=\frac{c}{t^{3/2}}$, where $c>0$. As a result, the convergence rate
term associated with the inexactness of the subproblem becomes
$O\left(\frac{c}{T^{3}}\right)$, which matches the convergence rate of the
Accelerated Cubic Newton method.

##  4 Lower bound

In this section, we present a novel lower bound for inexact second-order
methods with stochastic gradient and inexact (stochastic) Hessian. The proof
technique draws inspiration from the works [16, 51, 50]. For this section, we
assume that the function $f(x)$ is convex and has $L_{1}$-Lipschitz-continuous
gradient and $L_{2}$-Lipschitz-continuous Hessian.
To begin, we describe the information and structure of stochastic second-order
methods. At each point $x_{t}$, the oracle provide us with an unbiased
stochastic gradient $g_{t}=g(x_{t},\xi)$ and an inexact (stochastic) Hessian
$H_{t}=H(x_{t},\xi)$. The method can compute the minimum of the following
models:

$h_{t+1}=\textstyle{\operatornamewithlimits{argmin}}_{h}\left\\{\phi_{x_{t}}(h)=a_{1}\left\langle
g_{t},h\right\rangle+a_{2}\left\langle
H_{t}h,h\right\rangle+b_{1}\|h\|^{2}+b_{2}\|h\|^{3}\right\\}.$

Now, we formulate the main assumption regarding the method’s ability to
generate new points.

######  Assumption 4.1.

The method generates a recursive sequence of test points ${x_{t}}$ that
satisfies the following condition

$x_{t+1}\in x_{0}+\operatorname{Span}\left\\{h_{1},\ldots,h_{t+1}\right\\}$

Most first-order and second-order methods, including accelerated versions,
typically satisfy this assumption. However, we highlight that randomized
methods are not covered by this lower bound. Randomized lower bound even for
exact high-order methods is still an open problem. More details on randomized
lower bounds for first-order methods are presented in [66, 47]. Finally, we
present the main lower bound theorem for stochastic second-order methods.

######  Theorem 4.2.

Let some second-order method $\mathcal{M}$ with exactly solved subproblem
satisfies Assumption 4.1 and has access only to unbiased stochastic gradient
and inexact Hessian satisfying Assumption 2.2 or Assumption 2.3 with
$\sigma_{2}=\delta_{2}=\max\limits_{t=1,\ldots,T}\delta_{t}^{x_{t-1},x_{t}}$.
Assume the method $\mathcal{M}$ ensures for any function $f$ with
$L_{1}$-Lipschitz-continuous gradient and $L_{2}$-Lipschitz-continuous Hessian
the following convergence rate

$\min\limits_{0\leq t\leq T}\mathbb{E}\left[f(x_{t})-f(x^{\ast})\right]\leq
O(1)\max\left\\{\tfrac{\sigma_{1}R}{\Xi_{1}(T)};\tfrac{\sigma_{2}R^{2}}{\Xi_{2}(T)};\tfrac{L_{2}R^{3}}{\Xi_{3}(T)}\right\\}.$
(10)

Then for all $T\geq 1$ we haven

$\Xi_{1}(T)\leq\sqrt{T},\qquad\Xi_{2}(T)\leq T^{2},\qquad\Xi_{3}(T)\leq
T^{\frac{7}{2}}.$ (11)

###### Proof.

We prove this Theorem from contradiction. Let assume that there exist the
method $\mathcal{M}$ that satisfies conditions of the Theorem 4.2 and it is
faster in one of the bounds from (11).

The first case, $\Xi_{1}(T)>\sqrt{T}$ or $\Xi_{2}(T)>T^{2}$. Let us apply this
method for the first-order lower bound function. It is well-known, that for
the first-order methods, the lower bound is
$\Omega\left(\frac{\sigma_{1}R}{\sqrt{T}}+\frac{L_{1}R^{2}}{T^{2}}\right)$
[48]. Also, the first-order lower bound function has $0$-Lipschitz-continuous
Hessian. It means, that the method $\mathcal{M}$ can be applied for the first-
order lower-bound function. We fix stochastic Hessian oracle as
$H(x,\xi)=2L_{1}I$. It means that $\sigma_{2}=2L$ for such inexact Hessian.
With such matrix $H(x,\xi)=2L_{1}I$, the method $\mathcal{M}$ has only the
first-order information and lies in the class of first-order methods. Hence,
we apply the method $\mathcal{M}$ to the first-order lower bound function and
get the rate $\min\limits_{0\leq t\leq
T}\mathbb{E}\left[f(x_{t})-f(x^{\ast})\right]\leq
O(1)\max\left\\{\frac{\sigma_{1}R}{\Xi_{1}(T)};\frac{\sigma_{2}R^{2}}{\Xi_{2}(T)}\right\\}$,
where $\Xi_{1}(T)>\sqrt{T}$ or $\Xi_{2}(T)>T^{2}$. It means that we’ve got a
faster method than a lower bound. It is a contradiction, hence the rates for
the method $\mathcal{M}$ are bounded as $\Xi_{1}(T)\leq\sqrt{T},\Xi_{2}(T)\leq
T^{2}$. The second case, $\Xi_{3}(T)>T^{\frac{7}{2}}$. It is well-known, that
the deterministic second-order lower bound is
$\Omega\left(\frac{L_{2}R^{3}}{T^{\frac{7}{2}}}\right)$. Let us apply the
method $\mathcal{M}$ for the second-order lower bound function, where the
oracle give us exact gradients and exact Hessians, then $\sigma_{1}=0$,
$\sigma_{2}=0$ and the method $\mathcal{M}$ is in class of exact second-order
methods but converges faster than the lower bound. It is a contradiction,
hence the rate for the method $\mathcal{M}$ is bounded as $\Xi_{3}(T)\leq
T^{\frac{7}{2}}$. ∎

##  5 Tensor generalization

In this section we propose a tensor generalization of Algorithm 1. We start
with introducing the standard assumption on the objective $f$ for tensor
methods.

######  Assumption 5.1.

Function $f$ is convex, $p$ times differentiable on $\mathbb{R}^{d}$, and its
$p$-th derivative is Lipschitz continuous, i.e. for all $x,y\in\mathbb{R}^{d}$

$\|\nabla^{p}f(x)-\nabla^{p}f(y)\|\leq L_{2}\|x-y\|.$

We denote the $i$-th directional derivative of function $f$ at $x$ along
directions $s_{1},\ldots,s_{i}\in\mathbb{R}^{n}$ as
$\nabla^{i}f(x)[s_{1},\ldots,s_{i}].$ If all directions are the same we write
$\nabla^{i}f(x)[s]^{i}.$ For a $p$-th order tensor $U$, we denote by $\|U\|$
its tensor norm recursively induced [11] by the Euclidean norm on the space of
$p$-th order tensors:

$\|U\|=\max\limits_{\|s_{1}\|=\ldots=\|s_{p}\|=1}\\{|U[s_{1},\ldots,s_{p}]|\\},$

where $\|\cdot\|$ is the standard Euclidean norm.
We construct tensor methods based on the $p$-th order Taylor approximation of
the function $f(x)$, which can be written as follows:

$\Phi_{x,p}(y)\stackrel{{\scriptstyle\text{ def
}}}{{=}}f(x)+\textstyle{\sum}_{i=1}^{p}\tfrac{1}{i!}\nabla^{i}f(x)[y-x]^{i},\quad
y\in\mathbb{R}^{d}.$

Using approximations $G_{i}(x)$ for the derivatives
$\nabla^{i}f\left(x\right)$ we create an inexact $p$-th order Taylor series
expansion of the objective

$\phi_{x,p}(y)=f\left(x\right)+\textstyle{\sum}_{i=1}^{p}\tfrac{1}{i!}G_{i}(x)[y-x]^{i}.$

Next, we introduce a counterpart of Lemma 2.1 for high-order methods.

######  Lemma 5.2.

Let Assumption 5.1 hold. Then, for any $x,y\in\mathbb{R}^{d}$, we have

$\begin{split}|f(y)-\phi_{{x},p}(y)|\leq\textstyle{\sum\limits_{i=1}^{p}}\frac{1}{i!}\|(G_{i}(x)-\nabla^{i}f({x}))[y-x]^{i-1}\|\|y-x\|+\frac{L_{p}}{(p+1)!}\|y-x\|^{p+1},\end{split}$
$\begin{split}\|\nabla
f(y)-\nabla\phi_{x,p}(y)\|\leq\textstyle{\sum\limits_{i=1}^{p}}\frac{1}{(i-1)!}\|(G_{i}(x)-\nabla^{i}f({x}))[y-x]^{i-1}\|+\frac{L_{p}}{p!}\|y-x\|^{p},\end{split}$

where we use the standard convention $0!=1$.

Following the assumptions for the second-order method we introduce analogical
assumptions for high-order method.

######  Assumption 5.3 (Unbiased stochastic gradient with bounded variance and
bounded variance stochastic high-order derivatives).

For any $x\in\mathbb{R}^{d}$ stochastic gradient $G_{1}(x,\xi)$ and stochastic
high-order derivatives $G_{i}(x,\xi),~{}i=2,\ldots,p$ satisfy

$\displaystyle\mathbb{E}[G_{1}(x,\xi)\mid x]=\nabla
f(x),\quad\mathbb{E}\left[\|G_{1}(x,\xi)-\nabla f(x)\|^{2}\mid
x\right]\leq\sigma_{1}^{2},$ (12)
$\displaystyle\mathbb{E}\left[\|G_{i}(x,\xi)-\nabla^{i}f(x)\|^{2}\mid
x\right]\leq\sigma_{i}^{2},~{}i=2,\ldots,p.$

######  Assumption 5.4 (Unbiased stochastic gradient with bounded variance and
inexact high-order derivatives).

For any $x\in\mathbb{R}^{d}$ stochastic gradient $G_{1}(x,\xi)$ satisfy (12).
For given $x,~{}y\in\mathbb{R}^{d}$ inexact high-order derivatives
$G_{i}(x),~{}i=2,\ldots,p$ satisfy

$\displaystyle\|(G_{i}(x)-\nabla^{i}f({x}))[y-x]^{i-1}\|\leq\delta_{i}^{x,y}\|y-x\|^{i-1}.$

To extend Algorithm 1 to tensor methods, we introduce a $p$-th order model of
the function:

$\omega^{M,\bar{\delta}}_{x,p}(y):=\phi_{x,p}(y)+\tfrac{\bar{\delta}}{2}\|x-y\|^{2}+\textstyle{\sum\limits_{i=3}^{p}}\tfrac{\eta_{i}\delta_{i}}{i!}\|x-y\|^{i}+\tfrac{pM}{(p+1)!}\|x-y\|^{p+1},$

where $\eta_{i}>0,~{}3\leq i\leq p$. Next, we modify Definition 3.1 for the
high order derivatives case

######  Definition 5.5.

Denote by $S_{p}^{M,\bar{\delta},\tau}(x)$ a point
$S:=S_{p}^{M,\bar{\delta},\tau}(x)$ such that
$\|\nabla\omega_{x,p}^{M,\bar{\delta}}(S)\|\leq\tau.$

Now, we are prepared to introduce the method and state the convergence
theorem.

Algorithm 2 Accelerated Stochastic Tensor Method

1: Input: $y_{0}=x_{0}$ is starting point; constants $M\geq\frac{2}{p}L_{p}$;
$\eta_{i}\geq 4$, $3\leq i\leq p$; starting inexactness $\bar{\delta}_{0}\geq
0$; nonnegative nondecreasing sequences $\\{\bar{\kappa}^{t}_{i}\\}_{t\geq 0}$
for $i=2,\ldots,p+1$, and
$\alpha_{t}=\tfrac{p+1}{t+p+1},~{}~{}~{}A_{t}=\textstyle{\prod\limits_{j=1}^{t}(1-\alpha_{j})},~{}~{}~{}A_{0}=1.$
(13)
$\psi_{0}(x):=\tfrac{\bar{\kappa}_{2}^{0}+\lambda_{0}}{2}\|x-x_{0}\|^{2}+\textstyle{\sum\limits_{i=3}^{p}}\tfrac{\bar{\kappa}_{i}^{0}}{i!}\|x-x_{0}\|^{i}.$

2: for $t\geq 0$ do

3: $v_{t}=(1-\alpha_{t})x_{t}+\alpha_{t}y_{t},\quad
x_{t+1}=S_{p}^{M,\bar{\delta}_{t},\tau}(v_{t})$

4: Compute  $\displaystyle y_{t+1}=\arg\min_{x\in\mathbb{R}^{n}}$
$\displaystyle\left\\{\psi_{t+1}(x):=\psi_{t}(x)+\tfrac{\lambda_{t+1}-\lambda_{t}}{2}\|x-x_{0}\|^{2}\right.$
$\displaystyle\left.+\textstyle{\sum\limits_{i=2}^{p}}\tfrac{\bar{\kappa}^{t+1}_{i}-\bar{\kappa}^{t}_{i}}{i!}\|x-x_{0}\|^{i}+\tfrac{\alpha_{t}}{A_{t}}l(x,x_{t+1})\right\\}.$

######  Theorem 5.6.

Let Assumption 5.1 hold and $M\geq\frac{2}{p}L_{p}$.

  * •

Let Assumption 5.3 hold. After $T\geq 1$ with parameters

$\begin{gathered}\bar{\kappa}_{2}^{t}=O\left(\tfrac{\bar{\delta}_{t}\alpha_{t}^{2}}{A_{t}}\right),~{}\bar{\kappa}_{i}^{t+1}=O\left(\tfrac{\alpha_{t+1}^{i}\delta_{i}}{A_{t+1}}\right),~{}\bar{\kappa}_{p+1}^{t+1}=O\left(\tfrac{\alpha_{t+1}^{p+1}M}{A_{t+1}}\right),\\\
\lambda_{t}=O\left(\tfrac{\sigma_{1}}{R}t^{p+1/2}\right),~{}\delta_{t}=O\left(\sigma_{2}+\tfrac{\sigma_{1}+\tau}{R}t^{\frac{3}{2}}\right)\end{gathered}$
(14)

we get the following bound

$\displaystyle\mathbb{E}\left[f(x_{T})-f(x^{\ast})\right]$ $\displaystyle\leq
O\left(\tfrac{\tau
R}{\sqrt{T}}+\tfrac{\sigma_{1}R}{\sqrt{T}}+\textstyle\sum_{i=2}^{p}\tfrac{\sigma_{i}R^{i}}{T^{i}}+\tfrac{MR^{p+1}}{T^{p+1}}\right).$

  * •

Let Assumption 5.4 hold. After $T\geq 1$ with parameters defined in (7) and
$\sigma_{i}=\delta_{i}=\max\limits_{t=1,\ldots,T}\delta_{i,t}^{v_{t-1},x_{t}}$
we get the following bound

$\displaystyle\mathbb{E}\left[f(x_{T})-f(x^{\ast})\right]$ $\displaystyle\leq
O\left(\tfrac{\tau
R}{\sqrt{T}}+\tfrac{\sigma_{1}R}{\sqrt{T}}+\textstyle{\sum}_{i=2}^{p}\tfrac{\delta_{i}R^{i}}{T^{i}}+\tfrac{MR^{p+1}}{T^{p+1}}\right).$

##  6 Strongly convex case

######  Assumption 6.1.

Function $f$ is $\mu$-strongly convex, $p$ times differentiable on
$\mathbb{R}^{d}$, and its $p$-th derivative is Lipschitz continuous, i.e. for
all $x,y\in\mathbb{R}^{d}$

$\|\nabla^{p}f(x)-\nabla^{p}f(y)\|\leq L_{p}\|x-y\|.$

To exploit the strong convexity of the objective function and attain a linear
convergence rate, we introduce a restarted version of Restarted Accelerated
Stochastic Tensor Method (Algorithm 2). In each iteration of Restarted
Accelerated Stochastic Tensor Method (Algorithm 3), we execute Algorithm 2 for
a predetermined number of iterations as specified in equation
$\eqref{eq:restarts_iter}$. The output of this run is then used as the initial
point for the subsequent iteration of Algorithm 1, which resets the
parameters, and this process repeats iteratively.

Algorithm 3 Restarted Accelerated Stochastic Tensor Method

Input: $z_{0}\in\mathbb{R}^{d}$, strong convexity parameter $\mu>0$, $M\geq
L_{p}$, and $R_{0}>0$ such that $\|z_{0}-x^{*}\|\leq R_{0}$. For
$s=1,2,\ldots$:

  1. 1.

Set $x_{0}=z_{s-1}$, $r_{s-1}=\frac{R_{0}}{2^{s-1}}$, and
$R_{s-1}=\|z_{s-1}-x^{*}\|$.

  2. 2.

Run Algorithm 2 for $t_{s}$ iterations, where

$t_{s}=O(1)\textstyle{\max}\left\\{1,\left(\tfrac{\tau}{\mu
r_{s-1}}\right)^{2},\left(\tfrac{\sigma_{1}}{\mu
r_{s-1}}\right)^{2},\textstyle{\max\limits_{i=2,\ldots,p}}\left(\tfrac{\delta_{i}R_{s-1}^{i-2}}{\mu}\right)^{\frac{1}{i}},\left(\tfrac{L_{p}R^{p-1}_{s-1}}{\mu}\right)^{\frac{1}{p+1}}\right\\}.$
(15)

  3. 3.

Set $z_{s}=x^{t_{s}}$.

######  Theorem 6.2.

Let Assumption 6.1 hold and let parameters of Algorithm 1 be chosen as in
(14). Let $\\{z_{s}\\}_{s\geq 0}$ be generated by Algorithm 3 and $R>0$ be
such that $\|z_{0}-x^{*}\|\leq R$. Then for any $s\geq 0$ we have

$\displaystyle\mathbb{E}\|z_{s}-x^{*}\|^{2}$ $\displaystyle\leq 4^{-s}R^{2},$
$\displaystyle\mathbb{E}f(z_{s})-f(x^{*})$ $\displaystyle\leq 2^{-2s-1}\mu
R^{2}.$ (16)

Moreover, the total number of iterations to reach desired accuracy
$\varepsilon:~{}f(z_{s})-f(x^{*})\leq\varepsilon$ in expectation is

$O\left(\tfrac{(\tau+\sigma_{1})^{2}}{\mu\varepsilon}+\left(\sqrt{\tfrac{\sigma_{2}}{\mu}}+1\right)\log\tfrac{f(z_{0})-f(x^{*})}{\varepsilon}+\textstyle{\sum}_{i=3}^{p}\left(\tfrac{\sigma_{i}R^{i-2}}{\mu}\right)^{\frac{1}{i}}+\left(\tfrac{L_{p}R^{p-1}}{\mu}\right)^{\frac{1}{p+1}}\right).$

Now, let us make a few observations regarding the results obtained in Theorem
3.2. For simplicity let solution of the subproblem be exact and $p=2$, i.e. we
do the restarts of the Accelerated Stochastic Cubic Newton, so the total
number of iterations is

$O\left(\tfrac{\sigma_{1}^{2}}{\mu\varepsilon}+\left(\sqrt{\tfrac{\sigma_{2}}{\mu}}+1\right)\log\tfrac{f(z_{0})-f(x^{*})}{\varepsilon}+\left(\tfrac{L_{2}R}{\mu}\right)^{\frac{1}{3}}\right).$
(17)

Next, let’s consider solving the stochastic optimization problem
$\min\limits_{x\in\mathbb{R}^{d}}F(x)=\mathbb{E}[f(x,\xi)]$ using the mini-
batch Restarted Accelerated Stochastic Cubic Newton method (Algorithm 3) with
$p=2$. In this approach, the mini-batched stochastic gradient is computed as
$\tfrac{1}{r_{1}}\textstyle\sum_{i=1}^{r_{1}}\nabla f(x,\xi_{i})$ and the
mini-batched stochastic Hessian is computed as
$\tfrac{1}{r_{2}}\textstyle\sum_{i=1}^{r_{2}}\nabla^{2}f(x,\xi_{i})$, where
$r_{1}$ and $r_{2}$ represent the batch sizes for gradients and Hessians,
respectively.
From the convergence estimates in (16) and (17), we can determine the required
sample sizes for computing the batched gradients and batched Hessians.
Specifically, we have
$r_{1}=\tilde{O}\left(\tfrac{\sigma_{1}^{2}}{\varepsilon\mu^{2/3}}\right)$ and
$r_{2}=O\left(\tfrac{\sigma_{2}}{\mu^{1/3}}\right)$. Consequently, the overall
number of stochastic gradient computations is
$O\left(\tfrac{\sigma_{1}^{2}}{\varepsilon\mu^{2/3}}\right)$, which is similar
to the accelerated SGD method [24]. Interestingly, the number of stochastic
Hessian computations scales linearly with the desired accuracy $\varepsilon$,
i.e., $O\left(\tfrac{\sigma_{2}}{\mu^{1/3}}\log\frac{1}{\varepsilon}\right)$.
This result highlights the practical importance of second-order methods. Since
the batch size of the Hessian is constant, there is no need to adjust it as
the desired solution accuracy increases. This is particularly useful in
distributed optimization problems under the assumption of beta similarity
[68]. In methods with such assumption [68, 15, 2], the server stores a Hessian
sample that provides a "good" approximation of the exact Hessian of the
objective function. Algorithms utilize this approximation instead of
exchanging curvature information with the workers. The constant batch size
allows for accurately determining the necessary sample size to achieve fast
convergence to any desired accuracy.

## References

  * [1] A. Agafonov, P. Dvurechensky, G. Scutari, A. Gasnikov, D. Kamzolov, A. Lukashevich, and A. Daneshmand.  An accelerated second-order method for distributed stochastic optimization.  In 2021 60th IEEE Conference on Decision and Control (CDC), pages 2407–2413. IEEE, 2021.
  * [2] A. Agafonov, D. Kamzolov, P. Dvurechensky, and A. Gasnikov.  Inexact tensor methods and their application to stochastic convex optimization.  arXiv preprint arXiv:2012.15636, 2020.
  * [3] N. Agarwal and E. Hazan.  Lower bounds for higher-order convex optimization.  In S. Bubeck, V. Perchet, and P. Rigollet, editors, Proceedings of the 31st Conference On Learning Theory, volume 75 of Proceedings of Machine Learning Research, pages 774–792. PMLR, 06–09 Jul 2018.
  * [4] K. Antonakopoulos, A. Kavis, and V. Cevher.  Extra-newton: A first approach to noise-adaptive accelerated second-order methods.  arXiv preprint arXiv:2211.01832, 2022.
  * [5] Y. Arjevani, O. Shamir, and R. Shiff.  Oracle complexity of second-order methods for smooth convex optimization, 2017.
  * [6] M. Baes.  Estimate sequence methods: extensions and approximations.  Institute for Operations Research, ETH, Zürich, Switzerland, 2(1), 2009.
  * [7] S. Bellavia and G. Gurioli.  Stochastic analysis of an adaptive cubic regularization method under inexact gradient evaluations and dynamic hessian accuracy.  Optimization, 71(1):227–261, 2022.
  * [8] S. Bellavia, G. Gurioli, B. Morini, and P. L. Toint.  Adaptive regularization for nonconvex optimization using inexact function values and randomly perturbed derivatives.  Journal of Complexity, 68:101591, 2022.
  * [9] S. Bubeck, Q. Jiang, Y. T. Lee, Y. Li, and A. Sidford.  Near-optimal method for highly smooth convex optimization.  In Conference on Learning Theory, pages 492–507. PMLR, 2019.
  * [10] Y. Carmon, D. Hausler, A. Jambulapati, Y. Jin, and A. Sidford.  Optimal and adaptive monteiro-svaiter acceleration.  In A. H. Oh, A. Agarwal, D. Belgrave, and K. Cho, editors, Advances in Neural Information Processing Systems, 2022.
  * [11] C. Cartis, N. Gould, and P. Toint.  Improved second-order evaluation complexity for unconstrained nonlinear optimization using high-order regularized models.  08 2017.
  * [12] C. Cartis, N. I. Gould, and P. L. Toint.  Adaptive cubic regularisation methods for unconstrained optimization. part i: motivation, convergence and numerical results.  Mathematical Programming, 127(2):245–295, 2011.
  * [13] C. Cartis, N. I. Gould, and P. L. Toint.  Adaptive cubic regularisation methods for unconstrained optimization. part ii: worst-case function-and derivative-evaluation complexity.  Mathematical programming, 130(2):295–319, 2011.
  * [14] C. Cartis and K. Scheinberg.  Global convergence rate analysis of unconstrained optimization methods based on probabilistic models.  Mathematical Programming, 169:337–375, 2018.
  * [15] A. Daneshmand, G. Scutari, P. Dvurechensky, and A. Gasnikov.  Newton method over networks is fast up to the statistical precision.  In International Conference on Machine Learning, pages 2398–2409. PMLR, 2021.
  * [16] O. Devolder, F. Glineur, and Y. Nesterov.  First-order methods of smooth convex optimization with inexact oracle.  Mathematical Programming, 146:37–75, 2014.
  * [17] N. Doikov, E. M. Chayti, and M. Jaggi.  Second-order optimization with lazy hessians.  arXiv preprint arXiv:2212.00781, 2022.
  * [18] N. Doikov, K. Mishchenko, and Y. Nesterov.  Super-universal regularized newton method.  arXiv preprint arXiv:2208.05888, 2022.
  * [19] N. Doikov and Y. Nesterov.  Gradient regularization of newton method with bregman distances.  arXiv preprint arXiv:2112.02952, 2021.
  * [20] N. Doikov and Y. E. Nesterov.  Inexact tensor methods with dynamic accuracies.  In ICML, pages 2577–2586, 2020.
  * [21] P. Dvurechensky, D. Kamzolov, A. Lukashevich, S. Lee, E. Ordentlich, C. A. Uribe, and A. Gasnikov.  Hyperfast second-order local solvers for efficient statistically preconditioned distributed optimization.  EURO Journal on Computational Optimization, 10:100045, 2022.
  * [22] A. Gasnikov, P. Dvurechensky, E. Gorbunov, E. Vorontsova, D. Selikhanovych, and C. A. Uribe.  Optimal tensor methods in smooth convex and uniformly convexoptimization.  In Conference on Learning Theory, pages 1374–1391. PMLR, 2019.
  * [23] A. Gasnikov, P. Dvurechensky, E. Gorbunov, E. Vorontsova, D. Selikhanovych, C. A. Uribe, B. Jiang, H. Wang, S. Zhang, S. Bubeck, Q. Jiang, Y. T. Lee, Y. Li, and A. Sidford.  Near optimal methods for minimizing convex functions with lipschitz $p$-th derivatives.  In A. Beygelzimer and D. Hsu, editors, Proceedings of the Thirty-Second Conference on Learning Theory, volume 99 of Proceedings of Machine Learning Research, pages 1392–1393. PMLR, 25–28 Jun 2019.
  * [24] S. Ghadimi and G. Lan.  Optimal stochastic approximation algorithms for strongly convex stochastic composite optimization, ii: shrinking procedures and optimal algorithms.  SIAM Journal on Optimization, 23(4):2061–2089, 2013.
  * [25] S. Ghadimi, H. Liu, and T. Zhang.  Second-order methods with cubic regularization under inexact information.  arXiv preprint arXiv:1710.05782, 2017.
  * [26] H. Ghanbari and K. Scheinberg.  Proximal quasi-newton methods for regularized convex optimization with linear and accelerated sublinear convergence rates.  Computational Optimization and Applications, 69:597–627, 2018.
  * [27] G. Grapiglia and Y. Nesterov.  On inexact solution of auxiliary problems in tensor methods for convex optimization.  Optimization Methods and Software, 36(1):145–170, 2021.
  * [28] G. N. Grapiglia and Y. Nesterov.  Tensor methods for minimizing convex functions with hölder continuous higher-order derivatives.  SIAM Journal on Optimization, 30(4):2750–2779, 2020.
  * [29] A. Griewank.  The modification of Newton’s method for unconstrained optimization by bounding cubic terms.  Technical report, Department of Applied Mathematics and Theoretical Physics, University of Cambridge, 1981.  Technical Report NA/12.
  * [30] S. Hanzely, D. Kamzolov, D. Pasechnyuk, A. Gasnikov, P. Richtárik, and M. Takáč.  A damped newton method achieves global $\mathcal{O}\left(\frac{1}{k^{2}}\right)$ and local quadratic convergence rate.  In A. H. Oh, A. Agarwal, D. Belgrave, and K. Cho, editors, Advances in Neural Information Processing Systems, 2022.
  * [31] K. H. Hoffmann and H. J. Kornstaedt.  Higher-order necessary conditions in abstract mathematical programming.  Journal of Optimization Theory and Applications, 26(4):533–568, Dec 1978.
  * [32] B. Jiang, H. Wang, and S. Zhang.  An optimal high-order tensor method for convex optimization.  In Conference on Learning Theory, pages 1799–1801. PMLR, 2019.
  * [33] R. Jiang, Q. Jin, and A. Mokhtari.  Online learning guided curvature approximation: A quasi-newton method with global non-asymptotic superlinear convergence.  CoRR, abs/2302.08580, 2023.
  * [34] R. Jiang and A. Mokhtari.  Accelerated quasi-newton proximal extragradient: Faster rate for smooth convex optimization.  arXiv preprint arXiv:2306.02212, 2023.
  * [35] D. Kamzolov.  Near-optimal hyperfast second-order method for convex optimization.  In Y. Kochetov, I. Bykadorov, and T. Gruzdeva, editors, Mathematical Optimization Theory and Operations Research, pages 167–178, Cham, 2020. Springer International Publishing.
  * [36] D. Kamzolov, A. Gasnikov, P. Dvurechensky, A. Agafonov, and M. Takáč.  Exploiting higher-order derivatives in convex optimization methods.  arXiv preprint arXiv:2208.13190, 2022.
  * [37] D. Kamzolov, K. Ziu, A. Agafonov, and M. Takáč.  Accelerated adaptive cubic regularized quasi-newton methods.  arXiv preprint arXiv:2302.04987, 2023.
  * [38] L. V. Kantorovich.  On newton’s method.  Trudy Matematicheskogo Instituta imeni VA Steklova, 28:104–144, 1949.
  * [39] J. M. Kohler and A. Lucchi.  Sub-sampled cubic regularization for non-convex optimization.  In International Conference on Machine Learning, pages 1895–1904. PMLR, 2017.
  * [40] D. Kovalev and A. Gasnikov.  The first optimal acceleration of high-order methods in smooth convex optimization.  In A. H. Oh, A. Agarwal, D. Belgrave, and K. Cho, editors, Advances in Neural Information Processing Systems, 2022.
  * [41] G. Lan.  An optimal method for stochastic composite optimization.  Mathematical Programming, 133(1):365–397, 2012.
  * [42] A. Lucchi and J. Kohler.  A sub-sampled tensor method for nonconvex optimization.  IMA Journal of Numerical Analysis, 2022.
  * [43] K. Mishchenko.  Regularized newton method with global ${\boldsymbol{\mathcal{o}(1/{k}^{2})}}$ convergence.  SIAM Journal on Optimization, 33(3):1440–1462, 2023.
  * [44] R. D. C. Monteiro and B. F. Svaiter.  An accelerated hybrid proximal extragradient method for convex optimization and its implications to second-order methods.  SIAM Journal on Optimization, 23(2):1092–1125, 2013.
  * [45] J. J. Moré.  The Levenberg-Marquardt algorithm: implementation and theory.  In Numerical analysis, pages 105–116. Springer, 1978.
  * [46] A. Nemirovski, A. Juditsky, G. Lan, and A. Shapiro.  Robust stochastic approximation approach to stochastic programming.  SIAM Journal on optimization, 19(4):1574–1609, 2009.
  * [47] A. S. Nemirovski and D. B. Yudin.  Problem complexity and method efficiency in optimization.  1983\.
  * [48] A. Nemirovskii, D. B. Yudin, and E. Dawson.  Problem complexity and method efficiency in optimization.  1983\.
  * [49] Y. Nesterov.  Accelerating the cubic regularization of newton’s method on convex problems.  Mathematical Programming, 112(1):159–181, 2008.
  * [50] Y. Nesterov.  Lectures on convex optimization, volume 137.  Springer, 2018.
  * [51] Y. Nesterov.  Implementable tensor methods in unconstrained convex optimization.  Mathematical Programming, 186(1):157–183, 2021.
  * [52] Y. Nesterov.  Inexact high-order proximal-point methods with auxiliary search procedure.  SIAM Journal on Optimization, 31(4):2807–2828, 2021.
  * [53] Y. Nesterov.  Superfast second-order methods for unconstrained convex optimization.  Journal of Optimization Theory and Applications, 191(1):1–30, 2021\.
  * [54] Y. Nesterov and B. T. Polyak.  Cubic regularization of newton method and its global performance.  Mathematical Programming, 108(1):177–205, 2006.
  * [55] I. Newton.  Philosophiae naturalis principia mathematica.  Jussu Societatis Regiae ac Typis Josephi Streater, 1687.
  * [56] B. T. Polyak.  A new method of stochastic approximation type.  Avtomatika i telemekhanika, (7):98–107, 1990.
  * [57] B. T. Polyak and A. B. Juditsky.  Acceleration of stochastic approximation by averaging.  SIAM journal on control and optimization, 30(4):838–855, 1992.
  * [58] R. Polyak.  Complexity of the regularized newton method.  arXiv preprint arXiv:1706.08483, 2017.
  * [59] R. A. Polyak.  Regularized Newton method for unconstrained Convex optimization.  Mathematical Programming, 120(1):125–145, 2009.
  * [60] J. Raphson.  Analysis Aequationum Universalis Seu Ad Aequationes Algebraicas Resolvendas Methodus Generalis & Expedita, Ex Nova Infinitarum Serierum Methodo, Deducta Ac Demonstrata.  Th. Braddyll, 1697.
  * [61] H. Robbins and S. Monro.  A stochastic approximation method.  The annals of mathematical statistics, pages 400–407, 1951.
  * [62] K. Scheinberg and X. Tang.  Practical inexact proximal quasi-newton method with global complexity analysis.  Mathematical Programming, 160:495–529, 2016.
  * [63] D. Scieur.  Adaptive quasi-newton and anderson acceleration framework with explicit global (accelerated) convergence rates.  arXiv preprint arXiv:2305.19179, 2023.
  * [64] T. Simpson.  Essays on several curious and useful subjects, in speculative and mix’d mathematicks. Illustrated by a variety of examples. Printed by H. Woodfall, jun. for J. Nourse, at the Lamb without Temple-Bar, 1740.
  * [65] N. Tripuraneni, M. Stern, C. Jin, J. Regier, and M. I. Jordan.  Stochastic cubic regularization for fast nonconvex optimization.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 31, pages 2899–2908. Curran Associates, Inc., 2018.
  * [66] B. Woodworth and N. Srebro.  Lower bound for randomized first order convex optimization.  arXiv preprint arXiv:1709.03594, 2017.
  * [67] P. Xu, F. Roosta, and M. W. Mahoney.  Newton-type methods for non-convex optimization under inexact hessian information.  Mathematical Programming, 184(1):35–70, 2020.
  * [68] Y. Zhang and X. Lin.  Disco: Distributed optimization for self-concordant empirical loss.  In International conference on machine learning, pages 362–370. PMLR, 2015.

##  Appendix A Proof of Lemmas 2.1 and 5.2

######  Lemma 5.2.

Let Assumption 5.1 hold. Then, for any $x,y\in\mathbb{R}^{d}$, we have

$\begin{split}|f(y)-\phi_{{x},p}(y)|\leq\textstyle{\sum\limits_{i=1}^{p}}\frac{1}{i!}\|(G_{i}(x)-\nabla^{i}f({x}))[y-x]^{i-1}\|\|y-x\|+\frac{L_{p}}{(p+1)!}\|y-x\|^{p+1},\end{split}$
$\begin{split}\|\nabla
f(y)-\nabla\phi_{x,p}(y)\|\leq\textstyle{\sum\limits_{i=1}^{p}}\frac{1}{(i-1)!}\|(G_{i}(x)-\nabla^{i}f({x}))[y-x]^{i-1}\|+\frac{L_{p}}{p!}\|y-x\|^{p},\end{split}$

The proof of that lemma is the same as proofs of Lemmas 1, 2 from [2]. Lemma
2.1 is a special case of the Lemma 5.2 for $p=2$.

##  Appendix B Proof of Theorem 3.2

The full proof is organized as follows:

  * •

Lemma B.1 provides an upper bound for the estimating sequence $\psi_{t}(x)$;

  * •

Lemmas B.2, B.6 present the efficiency of Inexact Cubic Newton step
$x_{t+1}=S_{M,\delta_{t}}(v_{t})$.

  * •

Lemma B.6 provides a lower bound on $\psi_{t}(x)$ based on results of
technical Lemmas B.5\- B.4;

  * •

Everything is combined together in Theorem B.7 in order to prove convergence
and obtain convergence rate.

The following lemma shows that the sequence of functions $\bar{\psi}_{t}(x)$
can be upper bounded by the properly regularized objective function.

######  Lemma B.1.

For convex function $f(x)$ and $\psi_{t}(x)$, we have

$\psi_{t}(x^{\ast})\leq\frac{f(x^{\ast})}{A_{t-1}}+\frac{\bar{\kappa}_{2}^{t}+\lambda_{t}}{2}\|x^{\ast}-x_{0}\|^{2}+\frac{\bar{\kappa}^{t}_{3}}{6}\|x^{\ast}-x_{0}\|^{3}+err_{t}^{up},$
(18)

where

$err^{up}_{t}=\sum\limits_{j=0}^{t-1}\frac{\alpha_{j}}{A_{j}}\left\langle
g(x_{j+1})-\nabla f(x_{j+1}),x^{\ast}-x_{j+1}\right\rangle$ (19)

###### Proof.

For $t=0$, let us define $A_{-1}$ such that $\tfrac{1}{A_{-1}}=0$ then
$\tfrac{f(x^{\ast})}{A_{-1}}=0$ and

$\psi_{0}(x^{\ast})\leq\frac{\bar{\kappa}^{0}_{2}+\lambda_{0}}{2}\|x^{\ast}-x_{0}\|^{2}+\frac{\bar{\kappa}^{0}_{3}}{3}\|x^{\ast}-x_{0}\|^{3}$

From

$\psi_{t+1}(x):=\psi_{t}(x)+\tfrac{\lambda_{t+1}-\lambda_{t}}{2}\|x-x_{0}\|^{2}+\textstyle{\sum\limits_{i=2}^{3}}\tfrac{\bar{\kappa}^{t+1}_{i}-\bar{\kappa}^{t}_{i}}{i}\|x-x_{0}\|^{i}+\tfrac{\alpha_{t}}{A_{t}}l(x,x_{t+1}),$
(20)

we have

$\displaystyle\psi_{t}(x^{\ast})=\frac{\bar{\kappa}^{t}_{2}+\lambda_{t}}{2}\|x^{\ast}-x_{0}\|^{2}+\frac{\bar{\kappa}^{t}_{3}}{3}\|x^{\ast}-x_{0}\|^{3}+\sum\limits_{j=0}^{t-1}\frac{\alpha_{j}}{A_{j}}l(x^{\ast},x_{j+1}).$
(21)

From (1), we have that, for all $j\geq 1$, $A_{j}=A_{j-1}(1-\alpha_{j})$,
which leads to $\frac{\alpha_{j}}{A_{j}}=\frac{1}{A_{j}}-\frac{1}{A_{j-1}}$.
Hence, we have
$\sum\limits_{j=0}^{t-1}\frac{\alpha_{j}}{A_{j}}=\frac{1}{A_{t-1}}-\frac{1}{A_{-1}}=\frac{1}{A_{t-1}}$
and, using the convexity of the objective $f$, we get

$\displaystyle\sum\limits_{j=0}^{t-1}\frac{\alpha_{j}}{A_{j}}l(x^{\ast},x_{j+1})\leq\sum\limits_{j=0}^{t-1}\frac{\alpha_{j}}{A_{j}}\bar{l}(x^{\ast},x_{j+1})+\sum\limits_{j=0}^{t-1}\frac{\alpha_{j}}{A_{j}}\left\langle
g(x_{j+1})-\nabla f(x_{j+1}),x^{\ast}-x_{j+1}\right\rangle$ (22)
$\displaystyle\leq
f(x^{\ast})\sum\limits_{j=0}^{t-1}\frac{\alpha_{j}}{A_{j}}+err_{t}^{up}=\frac{f(x^{\ast})}{A_{t-1}}+err_{t}^{up},$
(23)

where $\bar{l}(x,y)=f(y)+\langle\nabla f(x),y-x\rangle$. Finally, combining
all the inequalities from above, we obtain

$\displaystyle\psi_{t}(x^{\ast})\stackrel{{\scriptstyle\eqref{eq:lem_upper_seq_pr2},\eqref{eq:lem_upper_seq_pr3}}}{{\leq}}\frac{f(x^{\ast})}{A_{t-1}}+\frac{\bar{\kappa}_{2}^{t}+\lambda_{t}}{2}\|x^{\ast}-x_{0}\|^{2}+\frac{\bar{\kappa}^{t}_{3}}{6}\|x^{\ast}-x_{0}\|^{3}+err_{t}^{up}.$
(24)

∎

######  Lemma B.2.

For the function $f(x)$ with $L_{2}$-Lipschitz-continuous Hessian and
$H(x_{t})$ is $\delta_{2}^{t}$-inexact Hessian for
$v_{t},x_{t+1}\in\mathbb{R}^{d}$ we have

$\|\nabla\phi_{v_{t}}(x_{t+1})-\nabla
f(x_{t+1})\|\leq\delta_{2}^{t}\|x_{t+1}-v_{t}\|+\frac{L_{2}}{2}\|x_{t+1}-v_{t}\|^{2}+\|g(v_{t})-\nabla
f(v_{t})\|,$ (25)

where we denote $\delta_{2}^{t}=\delta_{t}^{v_{t},x_{t+1}}$ to simplify the
notation.

###### Proof.

$\displaystyle\|\nabla\phi_{v_{t}}(x_{t+1})-\nabla f(x_{t+1})\|$
$\displaystyle=\|\nabla\phi_{v_{t}}(x_{t+1})-\Phi_{v_{t}}(x_{t+1})+\Phi_{v_{t}}(x_{t+1})-\nabla
f(x_{t+1})\|$
$\displaystyle\leq\|\nabla\phi_{v_{t}}(x_{t+1})-\Phi_{v_{t}}(x_{t+1})\|+\|\Phi_{v_{t}}(x_{t+1})-\nabla
f(x_{t+1})\|$
$\displaystyle=\|(\nabla^{2}f(v_{t})-B_{v_{t}})(x_{t+1}-v_{t})\|+\|\Phi_{v_{t}}(x_{t+1})-\nabla
f(x_{t+1})\|$ $\displaystyle+\|g(v_{t})-\nabla f(v_{t})\|$
$\displaystyle{\leq}\delta_{2}^{t}\|x_{t+1}-v_{t}\|+\frac{L_{2}}{2}\|x_{t+1}-v_{t}\|^{2}+\|g(v_{t})-\nabla
f(v_{t})\|$

∎

The next Lemma characterizes the progress of the inexact cubic step
$S_{M,\delta_{t}}(v_{t})$ in Algorithm 1.

######  Lemma B.3.

Let $\\{x_{t},v_{t}\\}_{t\geq 1}$ be generated by Algorithm 1. Then, for any
$\bar{\delta}_{t}\geq 4{\delta}_{2}^{t}$, $M\geq 4L_{2}$ and
$x_{t+1}=S_{M,\bar{\delta}_{t}}(v_{t})$, the following holds

$\begin{gathered}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{\tau^{2}}{\bar{\delta}_{t}}}+\frac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}2}}{\bar{\delta}_{t}}\|g(v_{t})-\nabla
f(v_{t})\|^{2}+\langle\nabla f(x_{t+1}),v_{t}-x_{t+1}\rangle\\\
\geq\min\left\\{\|\nabla
f(x_{t+1})\|^{2}\left(\tfrac{1}{4\delta_{t}}\right),\right.\left.\|\nabla
f(x_{t+1})\|^{\frac{3}{2}}\left(\tfrac{1}{3M}\right)^{\frac{1}{2}}\right\\}.\end{gathered}$
(26)

###### Proof.

For simplicity, we denote $r_{t+1}=\|x_{t+1}-v_{t}\|$ and

$\zeta_{t+1}=\bar{\delta}_{t}+\frac{M}{2}\|x_{t+1}-v_{t}\|.$ (27)

By Definition 3.1 for $x_{t+1}=S^{M,\bar{\delta}_{t},\tau}(v_{t})$

$\begin{gathered}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau\geq}\|\nabla\phi_{v_{t}}(x_{t+1})+\bar{\delta}_{t}(x_{t+1}-v_{t})+\frac{M}{2}\|x_{t+1}-v_{t}\|(x_{t+1}-v_{t})\|\\\
\stackrel{{\scriptstyle\eqref{big_sum_lemma6}}}{{=}}\|\nabla\phi_{v_{t}}(x_{t+1})+\zeta_{t+1}(x_{t+1}-v_{t})\|.\end{gathered}$
(28)

We start with getting an upper bound for $\|\nabla\phi_{v_{t}}(x_{t+1})-\nabla
f(x_{t+1})\|$.

$\|\nabla\phi_{v_{t}}(x_{t+1})-\nabla
f(x_{t+1})\|\stackrel{{\scriptstyle\eqref{eq:grad_model_bound}}}{{\leq}}\delta_{2}^{t}r_{t+1}+\frac{L_{2}}{2}r_{t+1}^{2}+\|g(v_{t})-\nabla
f(v_{t})\|.$ (29)

From inexact solution of the subproblem we get

$\begin{gathered}\|\nabla f(x_{k+1})+\zeta_{t+1}(x_{t+1}-v_{t})\|^{2}\\\
\leq\|\nabla\phi_{v_{t}}(x_{t+1})-\nabla
f(x_{t+1})-\nabla\phi_{v_{t}}(x_{t+1})-\zeta_{t+1}(x_{t+1}-v_{t})\|^{2}\\\
\leq 2\|\nabla\phi_{v_{t}}(x_{t+1})-\nabla
f(x_{t+1})\|^{2}+2\|\nabla\phi_{v_{t}}(x_{t+1})+\zeta_{t+1}(x_{t+1}-v_{t})\|^{2}\\\
\stackrel{{\scriptstyle\text{Def.~{}\ref{def:inexact_sub}}}}{{\leq}}2\|\nabla\phi_{v_{t}}(x_{t+1})-\nabla
f(x_{t+1})\|^{2}+2\tau^{2}\end{gathered}$ (30)

Next, from the previous inequality and (29), we get

$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}4}\|g(v_{t})-\nabla
f(v_{t})\|^{2}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}4}(\delta_{2}^{t}+\tfrac{L_{2}}{2}r_{t+1})^{2}r_{t+1}^{2}\stackrel{{\scriptstyle\eqref{eq:progress_step_error_upper_bound}}}{{\geq}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}2}\|\nabla\phi_{v_{t},p}(x_{t+1})-\nabla
f(x_{t+1})\|^{2}$
$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\stackrel{{\scriptstyle\eqref{eq:progress_step_inexact_upper_bound}}}{{\geq}}}\left\|\nabla
f(x_{t+1})+\zeta_{t+1}(x_{t+1}-v_{t})\right\|^{2}-{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}2\tau^{2}}$
$\displaystyle=2\langle\nabla
f(x_{t+1}),x_{t+1}-v_{t}\rangle\zeta_{t+1}+\|\nabla
f(x_{t+1})\|^{2}+\zeta_{t+1}^{2}\|x_{t+1}-v_{t}\|^{2}-{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color<EMAIL_ADDRESS>
Hence,

$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{\tau^{2}}{\zeta_{t+1}}}+\frac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}2}}{\zeta_{t+1}}\|g(v_{t})-\nabla
f(v_{t})\|^{2}+\langle\nabla
f(x_{t+1}),v_{t}-x_{t+1}\rangle\geq\frac{1}{2\zeta_{t+1}}\|\nabla
f(x_{t+1})\|^{2}+\frac{1}{2\zeta_{t+1}}\left(\zeta_{t+1}^{2}-4(\delta_{2}^{t}+\tfrac{L_{2}}{2}r_{t+1})^{2}\right)r_{t+1}^{2},$

and finally by using defenetion of $\zeta_{t+1}$, we get

$\displaystyle\frac{\tau^{2}}{\bar{\delta}_{t}}+\frac{2}{\bar{\delta}_{t}}\|g(v_{t})-\nabla
f(v_{t})\|^{2}+\langle\nabla
f(x_{t+1}),v_{t}-x_{t+1}\rangle\geq\frac{1}{2\zeta_{t+1}}\|\nabla
f(x_{t+1})\|^{2}+\frac{1}{2\zeta_{t+1}}\left(\zeta_{t+1}^{2}-4(\delta_{2}^{t}+\tfrac{L_{2}}{2}r_{t+1})^{2}\right)r_{t+1}^{2}$

Next, we consider $2$ cases depending on which term dominates in the
$\zeta_{t+1}$.

  * •

If $\bar{\delta}_{t}\geq\frac{M}{2}\|x_{t+1}-v_{t}\|$, then we get the
following bound

$\displaystyle\frac{\tau^{2}}{\bar{\delta}_{t}}+\frac{2}{\bar{\delta}_{t}}\|g(v_{t})-\nabla
f(v_{t})\|^{2}+\langle\nabla
f(x_{t+1}),v_{t}-x_{t+1}\rangle\geq\frac{1}{2\zeta_{t+1}}\|\nabla
f(x_{t+1})\|^{2}+\frac{1}{2\zeta_{t+1}}\left(\zeta_{t+1}^{2}-4({\delta_{2}^{t}}+\tfrac{L_{2}}{2}r_{t+1})^{2}\right)r_{t+1}^{2}$
$\displaystyle\geq\tfrac{1}{4\bar{\delta}_{t}}\|\nabla f(x_{t+1})\|^{2}$

  * •

If $\bar{\delta}_{t}<\frac{M}{2}\|x_{t+1}-v_{t}\|$ , then similarly to
previous case, we get

$\displaystyle\frac{\tau^{2}}{\bar{\delta}_{t}}+\frac{2}{\bar{\delta}_{t}}\|g(v_{t})-\nabla
f(v_{t})\|^{2}+\langle\nabla f(x_{t+1}),v_{t}-x_{t+1}\rangle$
$\displaystyle\geq\frac{\|\nabla
f(x_{t+1})\|^{2}}{2\zeta_{t+1}}+\left(\left(\bar{\delta}_{t}+\tfrac{M}{2}r_{t+1}\right)^{2}-4\left({\delta_{2}^{t}}+\tfrac{L_{2}}{2}r_{t+1}\right)^{2}\right)\frac{r_{t+1}^{2}}{2\zeta_{t+1}}$
$\displaystyle=\frac{\|\nabla
f(x_{t+1})\|^{2}}{2\zeta_{t+1}}+\left(\left(\bar{\delta}_{t}-2\delta_{2}^{t}+\tfrac{M-{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}2}L_{2}}{2}r_{t+1}\right)\left(\bar{\delta_{t}}+2\delta_{2}^{t}+\tfrac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}2}L_{2}+M}{2}r_{t+1}\right)\right)\frac{r_{t+1}^{2}}{2\zeta_{t+1}}$
$\displaystyle\geq\frac{\|\nabla
f(x_{t+1})\|^{2}}{2Mr_{t+1}}+\frac{M^{2}-{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}4}L_{2}^{2}}{4}\frac{r_{t+1}^{3}}{2M}\geq\frac{\|\nabla
f(x_{t+1})\|^{2}}{2Mr_{t+1}}+\frac{2M}{32}r_{t+1}^{3}$
$\displaystyle\geq\left(\frac{1}{3M}\right)^{\frac{1}{2}}\|\nabla
f(x_{t+1})\|^{\frac{3}{2}},$

where for the last inequality, we use $\tfrac{\alpha}{r}+\tfrac{\beta
r^{3}}{3}\geq\frac{4}{3}\beta^{1/4}\alpha^{3/4}$.

∎

We use the following lemma [2, Lemma 7].

######  Lemma B.4.

Let $h(x)$ be a convex function, $x_{0}\in\mathbb{R}^{n}$, $\theta_{i}\geq 0$
for $i=2,\ldots,p+1$ and

$\bar{x}=\arg\min\limits_{x\in\mathbb{R}^{n}}\\{\bar{h}(x)=h(x)+\sum\limits_{i=2}^{p+1}\theta_{i}d_{i}(x-x_{0})\\},$

where $d_{i}(x)=\tfrac{1}{i}\|x\|^{i}$ is a power-prox function. Then, for all
$x\in\mathbb{R}^{n}$,

$\bar{h}(x)\geq\bar{h}(\bar{x})+\sum\limits_{i=2}^{p+1}\left(\frac{1}{2}\right)^{i-2}\theta_{i}d_{i}(x-\bar{x}).$

We will also use the next technical lemma [49, 25] on Fenchel conjugate for
the $p$-th power of the norm.

######  Lemma B.5.

Let $g(z)=\frac{\theta}{p}\|z\|^{p}$ for $p\geq 2$ and $g^{*}$ be its
conjugate function i.e., $g^{*}(v)=\sup_{z}\\{\langle v,z\rangle-$ $g(z)\\}.$
Then, we have

$g^{*}(v)=\frac{p-1}{p}\left(\frac{\|v\|^{p}}{\theta}\right)^{\frac{1}{p-1}}$

Moreover, for any $v,z\in\mathbb{R}^{n}$, we have $g(z)+g^{*}(v)-\langle
z,v\rangle\geq 0.$

Finally, the last step is the next Lemma which prove that
$\frac{f(x_{t})}{A_{t-1}}\leq\min\limits_{x}\bar{\psi}_{t}(x)=\bar{\psi}_{t}^{\ast}$.

######  Lemma B.6.

Let $\\{x_{t},y_{t}\\}_{t\geq 1}$ be generated by Algorithm 1. Then

$\psi_{t}^{\ast}=\min\limits_{x}\psi_{t}(x)\geq\frac{f(x_{t})}{A_{t-1}}-err_{t}^{v}-err_{t}^{x}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-err_{t}^{\tau}},$
(31)

where

$err_{t}^{v}=\sum\limits_{j=0}^{t-1}\frac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}2}}{A_{j}\bar{\delta}_{j}}\|g(v_{j})-\nabla
f(v_{j})\|^{2},$ (32)
$err_{t}^{\tau}=\sum\limits_{j=0}^{t-1}\frac{\tau^{2}}{A_{j}\bar{\delta}_{j}},$
(33)

and

$err_{t}^{x}=\sum\limits_{j=0}^{t-1}\frac{\alpha_{j}^{2}}{2A_{j}^{2}\lambda_{j}}\|g(x_{j+1})-\nabla
f(x_{j+1})\|^{2}+\sum\limits_{j=0}^{t-1}\frac{\alpha_{j}}{A_{j}}\left\langle
g(x_{j+1})-\nabla f(x_{j+1}),y_{j}-x_{j+1}\right\rangle.$ (34)

###### Proof.

We prove Lemma by induction. Let us start with $t=0$, we define $A_{-1}$ such
that $\tfrac{1}{A_{-1}}=0$. Then $\tfrac{f(x_{0})}{A_{-1}}=0$ and
$\psi_{0}^{\ast}=0$, hence, $\psi_{0}^{\ast}\geq\tfrac{f(x_{0})}{A_{-1}}$. Let
us assume that (31) is true for $t$ and show that (31) is true for $t+1$. By
definition,

$\displaystyle\psi_{t}(x)=\frac{\lambda_{t}+\bar{\kappa}^{t}_{2}}{2}\|x-x_{0}\|^{2}+\frac{\bar{\kappa}^{t}_{3}}{3}\|x-x_{0}\|^{3}+\sum\limits_{j=0}^{t-1}\frac{\alpha_{j}}{A_{j}}l(x,x_{j+1}).$

Next, we apply Lemma B.4 with the following choice of parameters:
$h(x)=\sum\limits_{j=0}^{t-1}\frac{\alpha_{j}}{A_{j}}l(x,x_{j+1})$,
$\theta_{2}=\lambda_{t}+\bar{\kappa}^{t}_{2}$, and
$\theta_{3}=\bar{\kappa}^{t}_{3}$.
By (20),
$y_{t}=\operatornamewithlimits{argmin}\limits_{x\in\mathbb{R}^{d}}h(x)$, and
we have

$\displaystyle\psi_{t}(y_{t+1})\geq\psi_{t}^{\ast}+\frac{\bar{\kappa}^{t}_{2}+\lambda_{t}}{2}\|x-y_{t}\|^{2}+\frac{\bar{\kappa}^{t}_{3}}{6}\|x-y_{t}\|^{3}$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:lemma_ass}}}{{\geq}}\frac{f(x_{t})}{A_{t-1}}+\frac{\bar{\kappa}^{t}_{2}+\lambda_{t}}{2}\|x-y_{t}\|^{2}+\frac{\bar{\kappa}^{t}_{3}}{6}\|x-y_{t}\|^{3}-err_{t}^{v}-err_{t}^{x}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-err_{t}^{\tau}},$

where the last inequality follows from the assumption of the lemma.

By the definition of $\psi_{t+1}(x)$, the above inequality, and convexity of
$f$, we obtain

$\displaystyle\psi_{t+1}(y_{t+1})=\psi_{t}(y_{t+1})+\frac{\lambda_{t+1}-\lambda_{t}}{2}\|y_{t+1}-x_{0}\|^{2}+\frac{\bar{\kappa}^{t+1}_{2}-\bar{\kappa}_{2}^{t}}{2}\|y_{t+1}-x_{0}\|^{2}$
$\displaystyle+\frac{\bar{\kappa}^{t+1}_{3}-\bar{\kappa}_{3}^{t}}{3}\|y_{t+1}-x_{0}\|^{3}+\frac{\alpha_{t}}{A_{t}}l(y_{t+1},x_{t+1})$
$\displaystyle\geq\frac{f(x_{t})}{A_{t-1}}+\frac{\bar{\kappa}^{t}_{2}+\lambda_{t}}{2}\|y_{t+1}-y_{t}\|^{2}+\frac{\bar{\kappa}^{t}_{3}}{6}\|y_{t+1}-y_{t}\|^{3}+\frac{\alpha_{t}}{A_{t}}l(y_{t+1},x_{t+1})-err_{t}^{v}-err_{t}^{x}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-err_{t}^{\tau}}$
$\displaystyle=\frac{f(x_{t})}{A_{t-1}}+\frac{\bar{\kappa}^{t}_{2}+\lambda_{t}}{2}\|y_{t+1}-y_{t}\|^{2}+\frac{\bar{\kappa}^{t}_{3}}{6}\|y_{t+1}-y_{t}\|^{3}-err_{t}^{v}-err_{t}^{x}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-err_{t}^{\tau}}$
$\displaystyle\frac{\alpha_{t}}{A_{t}}(f(x_{t+1})+\left\langle\nabla
f(x_{t+1}),y_{t+1}-x_{t+1}\right\rangle+\left\langle g(x_{t+1})-\nabla
f(x_{t+1}),y_{t+1}-y_{t}\right\rangle+\left\langle g(x_{t+1})-\nabla
f(x_{t+1}),y_{t}-x_{t+1}\right\rangle)$
$\displaystyle\geq\frac{f(x_{t})}{A_{t-1}}+\frac{\bar{\kappa}^{t}_{2}}{2}\|y_{t+1}-y_{t}\|^{2}+\frac{\bar{\kappa}^{t}_{3}}{6}\|y_{t+1}-y_{t}\|^{3}-err_{t}^{v}-err_{t}^{x}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-err_{t}^{\tau}}$
$\displaystyle+\frac{\alpha_{t}}{A_{t}}\left(f(x_{t+1})+\left\langle\nabla
f(x_{t+1}),y_{t+1}-x_{t+1}\right\rangle\right)-\frac{\alpha_{t}^{2}}{2A_{t}^{2}\lambda_{t}}\|g(x_{t+1})-\nabla
f(x_{t+1})\|^{2}+\frac{\alpha_{t}}{A_{t}}\left\langle g(x_{t+1})-\nabla
f(x_{t+1}),y_{t}-x_{t+1}\right\rangle$
$\displaystyle=\frac{f(x_{t})}{A_{t-1}}+\frac{\bar{\kappa}^{t}_{2}}{2}\|y_{t+1}-y_{t}\|^{2}+\frac{\bar{\kappa}^{t}_{3}}{6}\|y_{t+1}-y_{t}\|^{3}-err_{t}^{v}-err_{t+1}^{x}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-err_{t}^{\tau}}+\frac{\alpha_{t}}{A_{t}}\left(f(x_{t+1})+\left\langle\nabla
f(x_{t+1}),y_{t+1}-x_{t+1}\right\rangle\right)$
$\displaystyle\geq\frac{\bar{\kappa}^{t}_{2}}{2}\|y_{t+1}-y_{t}\|^{2}+\frac{\bar{\kappa}^{t}_{3}}{6}\|y_{t+1}-y_{t}\|^{3}-err_{t}^{v}-err_{t+1}^{x}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-err_{t}^{\tau}}$
$\displaystyle+\frac{1}{A_{t-1}}(f(x_{t+1})+\langle\nabla
f(x_{t+1}),x_{t}-x_{t+1}\rangle)+\frac{\alpha_{t}}{A_{t}}\left(f(x_{t+1})+\left\langle\nabla
f(x_{t+1}),y_{t+1}-x_{t+1}\right\rangle\right).$

Next, we consider the sum of two linear models from the last inequality:

$\displaystyle\frac{1}{A_{t-1}}(f(x_{t+1})+\langle\nabla
f(x_{t+1}),x_{t}-x_{t+1}\rangle)+\frac{\alpha_{t}}{A_{t}}(f(x_{t+1})+\langle\nabla
f(x_{t+1}),y_{t+1}-x_{t+1}\rangle)$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:alphas}}}{{=}}\frac{1-\alpha_{t}}{A_{t}}f(x_{t+1})+\frac{1-\alpha_{t}}{A_{t}}\langle\nabla
f(x_{t+1}),x_{t}-x_{t+1}\rangle+\frac{\alpha_{t}}{A_{t}}f(x_{t+1})+\frac{\alpha_{t}}{A_{t}}\langle\nabla
f(x_{t+1}),y_{t+1}-x_{t+1}\rangle$
$\displaystyle{=}\frac{f(x_{t+1})}{A_{t}}+\frac{1-\alpha_{t}}{A_{t}}\langle\nabla
f(x_{t+1}),\frac{v_{t}-\alpha_{t}y_{t}}{1-\alpha_{t}}-x_{t+1}\rangle+\frac{\alpha_{t}}{A_{t}}\langle\nabla
f(x_{t+1}),y_{t+1}-x_{t+1}\rangle$
$\displaystyle=\frac{f(x_{t+1})}{A_{t}}+\frac{1}{A_{t}}\langle\nabla
f(x_{t+1}),v_{t}-x_{t+1}\rangle+\frac{\alpha_{t}}{A_{t}}\langle\nabla
f(x_{t+1}),y_{t+1}-y_{t}\rangle.$

As a result, by (20), we get

$\begin{gathered}\psi^{\ast}_{t+1}=\psi_{t+1}(y_{t+1})\geq\frac{f(x_{t+1})}{A_{t}}+\frac{1}{A_{t}}\langle\nabla
f(x_{t+1}),v_{t}-x_{t+1}\rangle+\frac{\bar{\kappa}^{t}_{2}}{2}\|y_{t+1}-y_{t}\|^{2}\\\
+\frac{\bar{\kappa}^{t}_{3}}{6}\|y_{t+1}-y_{t}\|^{3}+\frac{\alpha_{t}}{A_{t}}\langle\nabla
f(x_{t+1}),y_{t+1}-y_{t}\rangle-
err_{t}^{v}-err_{t+1}^{x}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color<EMAIL_ADDRESS>(35)

To complete the induction step, we show, that the sum of all terms in the RHS
except $\frac{f(x_{t+1})}{A_{t}}$ is non-negative (except err).

Lemma B.3 provides the lower bound for $\langle\nabla
f(x_{t+1}),v_{t}-x_{t+1}\rangle$. Let us consider the case when the minimum in
the RHS of (26) is attained at the first term. By Lemma B.5 with the following
choice of the parameters

$z=y_{t}-y_{t+1},~{}~{}v=\frac{\alpha_{t}}{A_{t}}\nabla
f(x_{t+1}),~{}~{}\theta=\bar{\kappa}_{i}^{t},$

we have

$\frac{\bar{\kappa}_{2}^{t}}{2}\|y_{t}-y_{t+1}\|^{2}+\frac{\alpha_{t}}{A_{t}}\langle\nabla
f(x_{t+1}),y_{t+1}-y_{t}\rangle\geq-\frac{1}{2}\left(\frac{\|\frac{\alpha_{t}}{A_{t}}\nabla
f(x_{t+1})\|^{2}}{\bar{\kappa}_{2}^{t}}\right).$ (36)

Hence,

$\displaystyle\frac{1}{A_{t}}\langle\nabla
f(x_{t+1}),v_{t}-x_{t+1}\rangle+\frac{\bar{\kappa}_{2}^{t}}{2}\|y_{t+1}-y_{t}\|^{2}+\frac{\alpha_{t}}{A_{t}}\langle\nabla
f(x_{t+1}),y_{t+1}-{y_{t}}\rangle$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:1case}}}{{\geq}}\frac{1}{A_{t}}\langle\nabla
f(x_{t+1}),v_{t}-x_{t+1}\rangle-\frac{\|\frac{\alpha_{t}}{A_{t}}\nabla
f(x_{t+1})\|^{i}}{2\bar{\kappa}_{2}^{t}}$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:scalar_lb_cases}}}{{\geq}}\frac{1}{A_{t}}\|\nabla
f(x_{t+1})\|^{2}\left(\frac{1}{4\bar{\delta}_{t}}\right)-\frac{\|\frac{\alpha_{t}}{A_{t}}\nabla
f(x_{t+1})\|^{2}}{2\bar{\kappa}_{2}^{t+1}}-\frac{2}{A_{t}\bar{\delta}_{t}}\|g(v_{t})-\nabla
f(v_{t})\|^{2}-\frac{\tau^{2}}{A_{t}\bar{\delta}_{t}}$
$\displaystyle\geq-\frac{2}{A_{t}\bar{\delta}_{t}}\|g(v_{t})-\nabla
f(v_{t})\|^{2}-\frac{\tau^{2}}{A_{t}\bar{\delta}_{t}},$

where the last inequality holds by our choice of the parameters

$\bar{\kappa}_{2}^{{t+1}}\geq\frac{2\bar{\delta}_{t}\alpha_{t}^{2}}{A_{t}}.$
(37)

Next, we consider the case when the minimum in the RHS of (26) is achieved on
the second term. Again, by Lemma B.5 with the same choice of $z,v$ and with
$\theta=\frac{\bar{\kappa}_{3}^{t}}{2}$, we have

$\frac{\bar{\kappa}_{3}^{t}}{6}\|y_{t}-y_{t+1}\|^{3}+\frac{\alpha_{t}}{A_{t}}\langle\nabla
f(x_{t+1}),y_{t+1}-{y_{t}}\rangle\geq-\frac{2}{3}\left(\frac{2\|\frac{\alpha_{t}}{A_{t}}\nabla
f(x_{t+1})\|^{3}}{\bar{\kappa}_{3}^{t}}\right)^{\frac{1}{2}}.$ (38)

Hence, we get

$\displaystyle\frac{1}{A_{t}}\langle\nabla
f(x_{t+1}),v_{t}-x_{t+1}\rangle+\frac{\bar{\kappa}_{3}^{t}}{6}\|y_{t+1}-y_{t}\|^{3}+\frac{\alpha_{t}}{A_{t}}\langle\nabla
f(x_{t+1}),y_{t+1}-{y_{t}}\rangle$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:2case}}}{{\geq}}\frac{1}{A_{t}}\langle\nabla
f(x_{t+1}),v_{t}-x_{t+1}\rangle-\frac{2}{3}\left(\frac{2\|\frac{\alpha_{t}}{A_{t}}\nabla
f(x_{t+1})\|^{3}}{\bar{\kappa}_{3}^{t}}\right)^{\frac{1}{2}}$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:scalar_lb_cases}}}{{\geq}}\frac{1}{A_{t}}\|\nabla
f(x_{t+1})\|^{\frac{3}{2}}\left(\frac{1}{3M}\right)^{\frac{1}{2}}-\frac{2}{3}\left(\frac{2\|\frac{\alpha_{t}}{A_{t}}\nabla
f(x_{t+1})\|^{3}}{\bar{\kappa}_{3}^{t}}\right)^{\frac{1}{2}}-\frac{2}{A_{t}\bar{\delta}_{t}}\|g(v_{t})-\nabla
f(v_{t})\|^{2}-\frac{\tau^{2}}{A_{t}\bar{\delta}_{t}}$
$\displaystyle\geq-\frac{2}{A_{t}\bar{\delta}_{t}}\|g(v_{t})-\nabla
f(v_{t})\|^{2}-\frac{\tau^{2}}{A_{t}\bar{\delta}_{t}},$

where the last inequality holds by our choice of $\bar{\kappa}_{3}^{t}$:

$\bar{\kappa}_{3}^{t}\geq\frac{8M}{3}\frac{\alpha_{t}^{3}}{A_{t}}.$ (39)

As a result, we unite both cases and get

$\displaystyle\psi^{\ast}_{t+1}\geq\frac{f(x_{t+1})}{A_{t}}+\frac{1}{A_{t}}\langle\nabla
f(x_{t+1}),v_{t}-x_{t+1}\rangle+\frac{\bar{\kappa}^{t}_{2}}{2}\|y_{t+1}-y_{t}\|^{2}$
$\displaystyle+\frac{\bar{\kappa}^{t}_{3}}{6}\|y_{t+1}-y_{t}\|^{3}+\frac{\alpha_{t}}{A_{t}}\langle\nabla
f(x_{t+1}),y_{t+1}-y_{t}\rangle-
err_{t}^{v}-err_{t+1}^{x}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-err_{t}^{\tau}}$
$\displaystyle\geq\frac{f(x_{t+1})}{A_{t}}-\frac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}2}}{A_{t}\bar{\delta}_{t}}\|g(v_{t})-\nabla
f(v_{t})\|^{2}-{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{\tau^{2}}{A_{t}\bar{\delta}_{t}}}-err_{t}^{v}-err_{t+1}^{x}$
$\displaystyle=\frac{f(x_{t+1})}{A_{t}}-err_{t+1}^{v}-err_{t+1}^{x}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-err_{t+1}^{\tau}}$

To sum up, by our choice of the parameters $\bar{\kappa}_{i}^{t}$, $i=2,3$, we
prove the induction step.

∎

Finally, we are in a position to prove the convergence rate theorem The proof
uses the following technical assumption. Let $R$ be such that

$\|x_{0}-x^{\ast}\|\leq R.$ (40)

######  Theorem B.7.

Let Assumption 1.1 hold and $M\geq 4L_{2}$. Let Assumption 2.3 hold. After
$T\geq 1$ with parameters defined in (7) and
$\sigma_{2}=\delta_{2}=\max\limits_{t=1,\ldots,T}\delta_{t}^{v_{t-1},x_{t}}$
we get the following bound for the objective residual

$\displaystyle\mathbb{E}\left[f(x_{T})-f(x^{\ast})\right]$
$\displaystyle\leq\frac{10\tau
R}{\sqrt{T+2}}+\frac{19\sigma_{1}R}{\sqrt{T+1}}+\frac{18\delta_{2}R^{2}}{(T+3)^{2}}+\frac{20L_{2}R^{3}}{(T+1)^{3}}.$

###### Proof.

First of all, let us bound $A_{T}$.

$\alpha_{t}=\frac{3}{t+3},~{}t\geq 1.$ (41)

Then, we have

$A_{T}=\prod_{t=1}^{T}\left(1-\alpha_{t}\right)=\prod_{t=1}^{T}\frac{t}{t+3}=\frac{T!3!}{(T+3)!}=\frac{6}{(T+1)(T+2)(T+3)}.$
(42)

And from [2] we get

$\sum_{t=0}^{T}\frac{A_{T}\alpha_{t}^{i}}{A_{t}}\leq\frac{3^{i}}{(T+3)^{i-1}}$
(43)

From Lemmas B.6 and B.1, we obtain that, for all $t\geq 1$,

$\displaystyle\frac{f(x_{t+1})}{A_{t}}-err_{t+1}^{v}-err_{t+1}^{x}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-err_{t+1}^{\tau}}\stackrel{{\scriptstyle\eqref{eq:lemma_ass}}}{{\leq}}\psi_{t+1}^{\ast}\leq\psi_{t+1}(x^{\ast})$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:acc_upper_bound}}}{{\leq}}\frac{f(x^{\ast})}{A_{t}}+\frac{\bar{\kappa}^{t}_{2}+\lambda_{t}}{2}\|x^{\ast}-x_{0}\|^{2}+\frac{\bar{\kappa}^{t}_{3}}{6}\|x^{\ast}-x_{0}\|^{3}+err_{t+1}^{up}.$

Next, we apply expectation

$\mathbb{E}\left[f(x_{T+1})-f(x^{\ast})\right]\leq
A_{T}\mathbb{E}\left[\frac{\bar{\kappa}_{2}^{T}+\lambda_{T}}{2}R^{2}+\frac{\bar{\kappa}_{3}^{T}}{6}R^{3}+err_{T+1}^{up}+err_{T+1}^{v}+err_{T+1}^{x}+err_{T+1}^{\tau}\right].$
(44)

Let us choose

$\displaystyle\bar{\delta}_{t}=\delta_{2}+\frac{\tau+\sigma}{R}(t+3)^{3/2},$
(45) $\displaystyle\lambda_{t}=\frac{\sigma}{R}(t+3)^{5/2}.$ (46)

Then, we bound terms in (68) step by step. We start from deterministic terms.

$\displaystyle
A_{T}\mathbb{E}\left[\frac{\bar{\kappa}_{2}^{T}+\lambda_{T}}{2}R^{2}+\frac{\bar{\kappa}_{3}^{T}}{6}R^{3}\right]\stackrel{{\scriptstyle\eqref{eq:params}}}{{=}}\frac{18\bar{\delta}_{T}R^{2}}{(T+3)^{2}}+\frac{6\lambda_{T}R^{2}}{(T+3)^{3}}+\frac{72MR^{3}}{(T+3)^{3}}$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:delta},\eqref{eq:lambda}}}{{=}}\frac{18\tau
R}{(T+3)^{1/2}}+\frac{18\sigma
R}{(T+3)^{1/2}}+\frac{18\delta_{2}R^{2}}{(T+3)^{2}}+\frac{6\sigma
R}{(T+3)^{1/2}}+\frac{72MR^{3}}{(T+3)^{3}}.$

Now, we bound expectation of all error terms. Firstly, we consider
$err_{T+1}^{up}$

$\displaystyle
A_{T}\mathbb{E}\left[err_{T+1}^{up}\right]\stackrel{{\scriptstyle\eqref{eq:err_up}}}{{=}}A_{T}\mathbb{E}\left[\sum\limits_{j=0}^{T}\frac{\alpha_{j}}{A_{j}}\left\langle
g(x_{j+1})-\nabla f(x_{j+1}),x^{\ast}-x_{j+1}\right\rangle\right]=0.$

Next, we bound $A_{T}\mathbb{E}\left[err_{T+1}^{v}\right]$

$\displaystyle
A_{T}\mathbb{E}\left[err_{T+1}^{v}\right]\stackrel{{\scriptstyle\eqref{eq:err_v}}}{{=}}A_{T}\mathbb{E}\left[\sum\limits_{j=0}^{T}\frac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}2}}{A_{j}\bar{\delta}_{j}}\|g(v_{j})-\nabla
f(v_{j})\|^{2}\right]\leq
2\sigma^{2}\sum\limits_{j=0}^{T}\frac{A_{T}}{A_{j}\bar{\delta}_{j}}$
$\displaystyle\stackrel{{\scriptstyle\eqref{alpha_t},\eqref{eq:delta}}}{{=}}\frac{2\sigma
R}{3^{3/2}}\sum\limits_{j=0}^{T}\frac{A_{T}\alpha_{j}^{3/2}}{A_{j}}\stackrel{{\scriptstyle\eqref{eq:alpha_sum_bound}}}{{\leq}}\frac{2\sigma
R}{(T+3)^{1/2}}$

Now we calculate $A_{T}\mathbb{E}\left[err_{T+1}^{x}\right]$

$\displaystyle
A_{T}\mathbb{E}\left[err_{T+1}^{x}\right]\stackrel{{\scriptstyle\eqref{eq:err_x}}}{{=}}A_{T}\mathbb{E}\left[\sum\limits_{j=0}^{T}\frac{\alpha_{j}^{2}}{2A_{j}^{2}\lambda_{j}}\|g(x_{j+1})-\nabla
f(x_{j+1})\|^{2}+\sum\limits_{j=0}^{T}\frac{\alpha_{j}}{A_{j}}\left\langle
g(x_{j+1})-\nabla f(x_{j+1}),y_{j}-x_{j+1}\right\rangle\right]$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:lambda}}}{{=}}\frac{\sigma
R}{2}\sum\limits_{j=0}^{T}\frac{A_{T}\alpha_{j}^{2}}{A_{j}^{2}(j+3)^{5/2}}\stackrel{{\scriptstyle\eqref{alpha_t}}}{{=}}\frac{\sigma
R}{3^{5/2}2}\sum\limits_{j=0}^{T}\frac{A_{T}\alpha_{j}^{9/2}}{A_{j}^{2}}\leq\frac{\sigma
R}{3^{5/2}2A_{T}}\sum\limits_{j=0}^{T}\frac{A_{T}\alpha_{j}^{9/2}}{A_{j}}$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:alpha_sum_bound}}}{{\leq}}\frac{3\sigma
R}{2A_{T}(T+3)^{7/2}}\stackrel{{\scriptstyle\eqref{eq:A_t_bound}}}{{\leq}}\frac{\sigma
R}{4(T+3)^{1/2}}.$

Finally, we consider $err_{T+1}^{\tau}$

$\displaystyle
A_{T}\mathbb{E}\left[err_{T+1}^{\tau}\right]\stackrel{{\scriptstyle\eqref{eq:err_tau}}}{{=}}\sum\limits_{j=0}^{T}\frac{A_{T}\tau^{2}}{A_{j}\bar{\delta}_{j}}\stackrel{{\scriptstyle\eqref{eq:delta}}}{{=}}\frac{\tau
R}{3^{3/2}}\sum\limits\frac{A_{T}\alpha_{t}^{3/2}}{A_{j}}\stackrel{{\scriptstyle\eqref{eq:alpha_sum_bound}}}{{\leq}}\frac{\tau
R}{(T+3)^{1/2}}.$

Combining all bounds from above we achieve convergence rate

$\displaystyle\mathbb{E}\left[f(x_{T+1})-f(x^{\ast})\right]\leq\frac{19\tau
R}{(T+3)^{1/2}}+\frac{27\sigma
R}{(T+3)^{1/2}}+\frac{18\delta_{2}R^{2}}{(T+3)^{2}}+\frac{72MR^{3}}{(T+3)^{3}}.$

∎

The case of stochastic Hessian (Theorem 3.2 under Assumption 2.3) can be
obtained in the same way by taking expectation in Lemma B.3.

##  Appendix C Proof of Theorem 5.6

Algorithm 2 is the tensor generalization of Algorithm 1. So, we will follow
the same steps as in Appendix B.

First of all, we provide tensor counterpart of Lemmas B.1, B.2, B.3. The
proofs directly follow the proofs of Lemmas for $2$-nd order case.

######  Lemma C.1.

For convex function $f(x)$ and $\psi_{t}(x)$, we have

$\psi_{t}(x^{\ast})\leq\frac{f(x^{\ast})}{A_{t-1}}+\frac{\bar{\kappa}_{2}^{t}+\lambda_{t}}{2}\|x^{\ast}-x_{0}\|^{2}+\sum\limits_{i=3}^{p+1}\frac{\bar{\kappa}^{t}_{i}}{i!}\|x^{\ast}-x_{0}\|^{i}+err_{t}^{up},$
(47)

where

$err^{up}_{t}=\sum\limits_{j=0}^{t-1}\frac{\alpha_{j}}{A_{j}}\left\langle
g(x_{j+1})-\nabla f(x_{j+1}),x^{\ast}-x_{j+1}\right\rangle$ (48)

######  Lemma C.2.

For the function $f(x)$ with $L_{p}$-Lipschitz-continuous Hessian and
$G_{i}(x_{t})$ is $\delta_{i}^{t}$-inexact $i$-th order derivative for
$v_{t},x_{t+1}\in\mathbb{R}^{d}$ we have

$\|\nabla\phi_{v_{t}}(x_{t+1})-\nabla
f(x_{t+1})\|\leq\sum\limits_{i=2}^{p}\frac{\delta_{i}^{t}}{(i-1)!}\|x_{t+1}-v_{t}\|^{i}+\frac{L_{p}}{p!}\|x_{t+1}-v_{t}\|^{p+1}+\|g(v_{t})-\nabla
f(v_{t})\|,$ (49)

where we denote $\delta_{2}^{t}=\delta_{t}^{v_{t},x_{t+1}}$ to simplify the
notation.

######  Lemma C.3.

Let $\\{x_{t},v_{t}\\}_{t\geq 1}$ be generated by Algorithm 1. Then, for any
$\bar{\delta}_{t}\geq 4{\delta}_{2}^{t}$, $\eta_{i}\geq 4$,
$M\geq\tfrac{2}{p}L_{p}$ and $x_{t+1}=S_{M,\bar{\delta}_{t}}(v_{t})$, the
following holds

$\begin{gathered}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{\tau^{2}}{\delta_{t}}}+\frac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}2}}{\delta_{t}}\|g(v_{t})-\nabla
f(v_{t})\|^{2}+\langle\nabla f(x_{t+1}),v_{t}-x_{t+1}\rangle\\\
\geq\min\left\\{\frac{2}{p}\|\nabla
f(x_{t+1})\|^{\frac{p+1}{p}}\left(\frac{(p-1)!}{M}\right)^{\frac{1}{p}};~{}\tfrac{1}{4\bar{\delta}_{t}}\|\nabla
f(x_{t+1})\|^{2};\right.\\\
\left.\min\limits_{i=3,\ldots,p}\left(\frac{2}{p}\|\nabla
f(x_{t+1})\|^{\frac{i}{i-1}}\left(\frac{(i-1)!}{\eta_{i}\delta_{i}^{t}}\right)^{\frac{1}{i-1}}\right)\right\\}.\end{gathered}$
(50)

###### Proof.

For simplicity, we denote $r_{t+1}=\|x_{t+1}-v_{t}\|$ and

$\zeta_{t+1}=\bar{\delta}_{t}+\sum\limits_{i=3}^{p}\frac{\eta_{i}\delta_{i}}{(i-1)!}\|x_{t+1}-v_{t}\|^{i-2}+\frac{M}{(p-1)!}\|x_{t+1}-v_{t}\|^{p-1}.$
(51)

By Definition 5.5 for $x_{t+1}=S_{p}^{M,\bar{\delta}_{t},\tau}(v_{t})$

$\begin{gathered}\tau^{2}\geq\left\|\nabla\phi_{v_{t},p}(x_{t+1})+\bar{\delta}_{t}(x_{t+1}-v_{t})+\sum\limits_{i=3}^{p}\frac{\eta_{i}\delta_{i}^{t}}{(i-1)!}\|x_{t+1}-v_{t}\|^{i-2}(x_{t+1}-v_{t})\right.\\\
\left.+\frac{M}{(p-1)!}\|x_{t+1}-v_{t}\|^{p-1}(x_{t+1}-v_{t})\right\|^{2}\\\
\stackrel{{\scriptstyle\eqref{eq:zeta_p}}}{{=}}\|\nabla\phi_{v_{t},p}+\zeta_{t+1}(x_{t+1}-v_{t})\|^{2}.\end{gathered}$
(52)

From inexact solution of subproblem we get

$\begin{gathered}\|\nabla f(x_{t+1})+\zeta_{t+1}(x_{t+1}-v_{t})\|^{2}\\\
\leq\|\nabla\phi_{v_{t},p}(x_{t+1})-\nabla
f(x_{t+1})-\nabla\phi_{v_{t},p}(x_{t+1})-\zeta_{t+1}(x_{t+1}-v_{t})\|^{2}\\\
\leq 2\|\nabla\phi_{v_{t},p}(x_{t+1})-\nabla
f(x_{t+1})\|^{2}+2\|\nabla\phi_{v_{t},p}(x_{t+1})+\zeta_{t+1}(x_{t+1}-v_{t})\|^{2}\\\
\leq 2\|\nabla\phi_{v_{t},p}(x_{t+1})-\nabla
f(x_{t+1})\|^{2}+2\tau^{2}.\end{gathered}$ (53)

Next, from previous inequality and Lemma C.2

$\begin{gathered}4\|g(v_{t})-\nabla
f(v_{t})\|^{2}+4\left(\sum\limits_{i=2}^{p}\frac{\delta_{i}^{t}}{(i-1)!}r_{t+1}^{i}+\frac{L_{p}}{p!}r_{t+1}^{p+1}\right)^{2}+2\tau^{2}\\\
\geq 2\|\nabla\phi_{v_{t},p}(x_{t+1}-v_{t})\|+2\tau^{2}\\\ \geq\|\nabla
f(x_{t+1})+\zeta_{t+1}(x_{t+1}-v_{t})\|\\\ \geq\|\nabla
f(x_{t+1})\|^{2}+2\zeta_{t+1}\langle\nabla
f(x_{t+1}),x_{t+1}-v_{t}\rangle+\zeta_{t+1}^{2}\|x_{t+1}-v_{t}\|^{2}.\end{gathered}$
(54)

Hence,

$\displaystyle\frac{\tau^{2}}{\zeta_{t+1}}+\frac{2}{\zeta_{t+1}}\|g(v_{t})-\nabla
f(v_{t})\|^{2}+\langle\nabla f(x_{t+1}),v_{t}-x_{t+1}\rangle\geq$
$\displaystyle\frac{1}{2\zeta_{t+1}}\|\nabla
f(x_{t+1})\|^{2}+\frac{1}{2\zeta_{t+1}}\left(\zeta_{t+1}^{2}-4\left(\sum\limits_{i=2}^{p}\frac{\delta_{i}^{t}}{(i-1)!}r_{t+1}^{i-2}+\frac{L_{p}}{p!}r_{t+1}^{p-1}\right)^{2}\right)r_{t+1}^{2},$

and finally by using definition of $\zeta_{t+1}$, we get

$\displaystyle\frac{\tau^{2}}{\bar{\delta}_{t}}+\frac{2}{\bar{\delta}_{t}}\|g(v_{t})-\nabla
f(v_{t})\|^{2}+\langle\nabla f(x_{t+1}),v_{t}-x_{t+1}\rangle\geq$
$\displaystyle\frac{1}{2\zeta_{t+1}}\|\nabla
f(x_{t+1})\|^{2}+\frac{1}{2\zeta_{t+1}}\left(\zeta_{t+1}^{2}-4\left(\sum\limits_{i=2}^{p}\frac{\delta_{i}^{t}}{(i-1)!}r_{t+1}^{i-2}+\frac{L_{p}}{p!}r_{t+1}^{p-1}\right)^{2}\right)r_{t+1}^{2}.$

Next, we consider $p$ cases depending on which term dominates in
$\zeta_{t+1}$.

  * •

If
$\bar{\delta}_{t}\geq\sum\limits_{i=3}^{p}\frac{\eta_{i}\delta_{i}^{t}}{(i-1)!}r_{t+1}^{i-2}+\frac{M}{(p-1)!}r_{t+1}^{p-1}$,
then we get the following bound

$\displaystyle\frac{\tau^{2}}{\bar{\delta}_{t}}+\frac{2}{\bar{\delta}_{t}}\|g(v_{t})-\nabla
f(v_{t})\|^{2}+\langle\nabla f(x_{t+1}),v_{t}-x_{t+1}\rangle$
$\displaystyle\geq\frac{1}{2\zeta_{t+1}}\|\nabla
f(x_{t+1})\|^{2}+\frac{1}{2\zeta_{t+1}}\left(\zeta_{t+1}^{2}-4\left(\sum\limits_{i=2}^{p}\frac{\delta_{i}^{t}}{(i-1)!}r_{t+1}^{i-2}+\frac{L_{p}}{p!}r_{t+1}^{p-1}\right)^{2}\right)r_{t+1}^{2}$
$\displaystyle\geq\tfrac{1}{2p\bar{\delta}_{t}}\|\nabla f(x_{t+1})\|^{2}.$

  * •

If
$\frac{\eta_{i}\delta_{i}^{t}}{(i-1)!}r_{t+1}^{i-1}\geq\bar{\delta}_{t}+\sum\limits_{j=3,j\neq
i}^{p}\frac{\eta_{j}\delta_{j}}{(j-1)!}r_{t+1}^{j-1}+\frac{M}{(p-1)!}r_{t+1}^{p-1}$
for $i,~{}3\leq i\leq p$, we get

$\displaystyle\frac{\tau^{2}}{\bar{\delta}_{t}}+\frac{2}{\bar{\delta}_{t}}\|g(v_{t})-\nabla
f(v_{t})\|^{2}+\langle\nabla f(x_{t+1}),v_{t}-x_{t+1}\rangle$
$\displaystyle\geq\frac{1}{2\zeta_{t+1}}\|\nabla
f(x_{t+1})\|^{2}+\frac{1}{2\zeta_{t+1}}\left(\zeta_{t+1}^{2}-4\left(\sum\limits_{i=2}^{p}\frac{\delta_{i}^{t}}{(i-1)!}r_{t+1}^{i}+\frac{L_{p}}{p!}r_{t+1}^{p+1}\right)^{2}\right)r_{t+1}^{2}$
$\displaystyle\geq\frac{(i-1)!\|\nabla
f(x_{t+1})\|^{2}}{p\eta_{i}\delta_{i}^{t}r_{t+1}^{i-1}}+\left(\bar{\delta}_{t}-2\delta_{2}^{t}+\sum\limits_{i=3}^{p}\frac{\eta_{i}\delta_{i}^{t}-2\delta_{i}^{t}}{(i-1)!}r_{t+1}^{i-2}+\frac{pM-2L_{p}}{p!}r_{t+1}^{p-1}\right)$
$\displaystyle\times\left(\bar{\delta}_{t}+2\delta_{2}^{t}+\sum\limits_{i=3}^{p}\frac{\eta_{i}\delta_{i}^{t}+2\delta_{i}^{t}}{(i-1)!}r_{t+1}^{i-2}+\frac{pM+2L_{p}}{p!}r_{t+1}^{p-1}\right)\frac{r_{t+1}^{2}}{\zeta_{t+1}}$
$\displaystyle\geq\frac{(i-1)!\|\nabla
f(x_{t+1})\|^{2}}{p\eta_{i}\delta_{i}^{t}r_{t+1}^{i-2}}+\frac{\eta_{i}\delta_{i}^{t}-2\delta_{i}^{t}}{(i-1)!}\frac{\eta_{i}\delta_{i}^{t}+2\delta_{i}^{t}}{(i-1)!}\frac{r_{t+1}^{i}(i-1)!}{p\eta_{i}\delta_{i}^{t}}$
$\displaystyle\geq\frac{(i-1)!\|\nabla
f(x_{t+1})\|^{2}}{p\eta_{i}\delta_{i}^{t}r_{t+1}^{i-2}}+\frac{\eta_{i}\delta_{i}^{t}}{(i-1)!p}r_{t+1}^{i}.$
$\displaystyle\geq\frac{2}{i}\left(\frac{(i-2)(i-1)!\|\nabla
f(x_{t+1})\|^{2}}{p\eta_{i}\delta_{i}^{t}}\right)^{\frac{i}{2(i-1)}}\left(\frac{i\eta_{i}\delta_{i}^{t}}{(i-1)!p}\right)^{\frac{i-2}{2(i-1)}}$
$\displaystyle\geq\frac{2}{p}\|\nabla
f(x_{t+1})\|^{\frac{i}{i-1}}\left(\frac{(i-1)!}{\eta_{i}\delta_{i}^{t}}\right)^{\frac{1}{i-1}},$

where we used $\frac{\alpha}{(i-2)r^{i-2}}+\frac{\beta
r^{i}}{i}\geq\frac{2}{i}\alpha^{\frac{i}{2(i-1)}}\beta^{\frac{i-2}{2(i-1)}}$.

  * •

If
$\frac{M}{(p-1)!}r_{t+1}^{p-1}\geq\bar{\delta}_{t}+\sum\limits_{i=3}^{p}\frac{\eta_{i}\delta_{i}^{t}}{(i-1)!}r_{t+1}^{i-1}$,
then similarly to previous case, we get

$\displaystyle\frac{\tau^{2}}{\bar{\delta}_{t}}+\frac{2}{\bar{\delta}_{t}}\|g(v_{t})-\nabla
f(v_{t})\|^{2}+\langle\nabla f(x_{t+1}),v_{t}-x_{t+1}\rangle$
$\displaystyle\geq\frac{1}{2\zeta_{t+1}}\|\nabla
f(x_{t+1})\|^{2}+\frac{1}{2\zeta_{t+1}}\left(\zeta_{t+1}^{2}-4\left(\sum\limits_{i=2}^{p}\frac{\delta_{i}^{t}}{(i-1)!}r_{t+1}^{i}+\frac{L_{p}}{p!}r_{t+1}^{p+1}\right)^{2}\right)r_{t+1}^{2}$
$\displaystyle\frac{(p-1)!}{pMr_{t+1}^{p-1}}\|\nabla
f(x_{t+1})\|^{2}+\frac{pM}{2p(p!)}r_{t+1}^{p+1}$
$\displaystyle\geq\frac{2}{p}\|\nabla
f(x_{t+1})\|^{\frac{p+1}{p}}\left(\frac{(p-1)!}{M}\right)^{\frac{1}{p}}.$

∎

Next we will need technical Lemmas B.4, B.5.

######  Lemma C.4.

Let $\\{x_{t},y_{t}\\}_{t\geq 1}$ be generated by Algorithm 2. Then

$\psi_{t}^{\ast}=\min\limits_{x}\psi_{t}(x)\geq\frac{f(x_{t})}{A_{t-1}}-err_{t}^{v}-err_{t}^{x}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-err_{t}^{\tau}},$
(55)

where

$err_{t}^{v}=\sum\limits_{j=0}^{t-1}\frac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}2}}{A_{j}\bar{\delta}_{j}}\|g(v_{j})-\nabla
f(v_{j})\|^{2},$ (56)
$err_{t}^{\tau}=\sum\limits_{j=0}^{t-1}\frac{\tau^{2}}{A_{j}\bar{\delta}_{j}},$
(57)

and

$err_{t}^{x}=\sum\limits_{j=0}^{t-1}\frac{\alpha_{j}^{2}}{2A_{j}^{2}\lambda_{j}}\|g(x_{j+1})-\nabla
f(x_{j+1})\|^{2}+\sum\limits_{j=0}^{t-1}\frac{\alpha_{j}}{A_{j}}\left\langle
g(x_{j+1})-\nabla f(x_{j+1}),y_{j}-x_{j+1}\right\rangle.$ (58)

###### Proof.

We prove Lemma by induction. Let us start with $t=0$, we define $A_{-1}$ such
that $\tfrac{1}{A_{-1}}=0$. Then $\tfrac{f(x_{0})}{A_{-1}}=0$ and
$\psi_{0}^{\ast}=0$, hence, $\psi_{0}^{\ast}\geq\tfrac{f(x_{0})}{A_{-1}}$. Let
us assume that (55) is true for $t$ and show that (55) is true for $t+1$.

Following the steps of Lemma B.6 for $p$-th order case, we get

$\begin{gathered}\psi^{\ast}_{t+1}=\psi_{t+1}(y_{t+1})\geq\frac{f(x_{t+1})}{A_{t}}+\frac{1}{A_{t}}\langle\nabla
f(x_{t+1}),v_{t}-x_{t+1}\rangle\\\
+\sum\limits_{i=2}^{p+1}\left(\frac{1}{2}\right)^{i-2}\frac{\bar{\kappa}_{i}^{t}}{(i-1)!}d_{i}(y_{t+1}-y_{t})+\frac{\alpha_{t}}{A_{t}}\langle\nabla
f(x_{t+1}),y_{t+1}-y_{t}\rangle-
err_{t}^{v}-err_{t+1}^{x}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color<EMAIL_ADDRESS>(59)

To complete the induction step we will show, that the sum of all terms in the
RHS except $\frac{f(x_{t+1})}{A_{t}}$ and error terms is non-negative.

Lemma C.3 provides the lower bound for $\langle\nabla
f(x_{t+1}),v_{t}-x_{t+1}\rangle$. Let us consider the case when the minimum in
the RHS of (50) is attained at the term with particular $i=3,\ldots,p$. By
Lemma B.5 with the following choice of the parameters

$z=y_{t}-y_{t+1},~{}~{}v_{t}=\frac{\alpha_{t}}{A_{t}}\nabla
f(x_{t+1}),~{}~{}\theta=\left(\frac{1}{2}\right)^{i-2}\frac{\bar{\kappa}_{i}^{t}}{(i-1)!},$

we have

$\frac{1}{i}\left(\frac{1}{2}\right)^{i-2}\frac{\bar{\kappa}_{i}^{t}}{(i-1)!}\|y_{t}-y_{t+1}\|^{i}+\frac{\alpha_{t}}{A_{t}}\langle\nabla
f(x_{t+1}),y_{t+1}-y_{t}\rangle\geq-\frac{i-1}{i}\left(\frac{\|\frac{\alpha_{t}}{A_{t}}\nabla
f(x_{t+1})\|^{i}}{\left(\frac{1}{2}\right)^{i-2}\frac{\bar{\kappa}_{i}^{{t}}}{(i-1)!}}\right)^{\frac{1}{i-1}}.$
(60)

Hence,

$\displaystyle\frac{f(x_{t+1})}{A_{t}}+\frac{1}{A_{t}}\langle\nabla
f(x_{t+1}),v_{t}-x_{t+1}\rangle+\left(\frac{1}{{2}}\right)^{i-2}\frac{\bar{\kappa}_{i}^{t}}{(i-1)!}d_{i}(y_{t+1}-y_{t})+\frac{\alpha_{t}}{A_{t}}\langle\nabla
f(x_{t+1}),y_{t+1}-{y_{t}}\rangle$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:1case}}}{{\geq}}\frac{f(x_{t+1})}{A_{t}}+\frac{1}{A_{t}}\langle\nabla
f(x_{t+1}),v_{t}-x_{t+1}\rangle-\frac{i-1}{i}\left(\frac{\|\frac{\alpha_{t}}{A_{t}}\nabla
f(x_{t+1})\|^{i}}{\left(\frac{1}{{2}}\right)^{i-{2}}\frac{\bar{\kappa}_{i}^{{t}}}{{(i-1)}!}}\right)^{\frac{1}{i-1}}$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:scalar_lb_cases_p}}}{{\geq}}\frac{f(x_{t+1})}{A_{t}}+\frac{2}{p}\|\nabla
f(x_{t+1})\|^{\frac{i}{i-1}}\left(\frac{(i-1)!}{\eta_{i}\delta_{i}^{t}}\right)^{\frac{1}{i-1}}-\frac{i-1}{i}\left(\frac{\|\frac{\alpha_{t}}{A_{t}}\nabla
f(x_{t+1})\|^{i}}{\left(\frac{1}{{2}}\right)^{i-{2}}\frac{\bar{\kappa}_{i}^{{t}}}{{(i-1)}!}}\right)^{\frac{1}{i-1}}$
$\displaystyle\geq\frac{f(x_{t+1})}{A_{t}},$

where the last inequality holds by our choice of the parameters

$\bar{\kappa}_{i}^{{t}}\geq\frac{p^{i-1}}{2}\frac{\alpha_{t}^{i}}{A_{t}}\eta_{i}\delta_{i}^{t}.$
(61)

Let us consider the case when the minimum in the RHS of (50) is achieved on
the second term. Following similar steps, we get

$\bar{\kappa}_{2}^{t}\geq\frac{2p\alpha_{t}^{2}}{A_{t}}\bar{\delta}_{t}.$ (62)

Next, we consider the case when the minimum in the RHS of (50) is achieved on
the first term. Again, by Lemma B.5 with the same choice of $z,v$ and with
$\theta={\left(\frac{1}{2}\right)^{p-1}\frac{\bar{\kappa}_{p+1}^{t-1}}{(p-1)!}}$,
we have

${\frac{1}{p}}\left(\frac{1}{{2}}\right)^{p-1}\frac{\bar{\kappa}_{p+1}^{t-1}}{{{p}!}}\|y_{t}-y_{t+1}\|^{p+1}+\frac{\alpha_{t}}{A_{t}}\langle\nabla
f(x_{t+1}),y_{t+1}-{y_{t}}\rangle\geq-\frac{p}{p+1}\left(\frac{\|\frac{\alpha_{t}}{A_{t}}\nabla
f(x_{t+1})\|^{p+1}}{{\left(\frac{1}{2}\right)^{p-1}\frac{\bar{\kappa}_{p+1}^{t}}{(p-1)!}}}\right)^{\frac{1}{p}}.$
(63)

Hence, we get

$\displaystyle\frac{f(x_{t+1})}{A_{t}}+\frac{1}{A_{t}}\langle\nabla
f(x_{t+1}),v_{t}-x_{t+1}\rangle+\frac{1}{2^{p-1}}\frac{\bar{\kappa}_{p+1}^{t-1}}{{p}!}d_{p+1}(y_{t+1}-y_{t})+\frac{\alpha_{t}}{A_{t}}\langle\nabla
f(x_{t+1}),y_{t+1}-{y_{t}}\rangle$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:2case_p}}}{{\geq}}\frac{f(x_{t+1})}{A_{t}}+\frac{1}{A_{t}}\langle\nabla
f(x_{t+1}),v_{t}-x_{t+1}\rangle-\frac{p}{p+1}\left(\frac{\|\frac{\alpha_{t}}{A_{t}}\nabla
f(x_{t+1})\|^{p+1}}{{\left(\frac{1}{2}\right)^{p-1}\frac{\bar{\kappa}_{p+1}^{t}}{{p}!}}}\right)^{\frac{1}{p}}$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:scalar_lb_cases_p}}}{{\geq}}\frac{f(x_{t+1})}{A_{t}}+\frac{1}{A_{t}}\frac{2}{p}\|\nabla
f(x_{t+1})\|^{\frac{p+1}{p}}\left(\frac{(p-1)!}{M}\right)^{\frac{1}{p}}-\frac{p}{p+1}\left(\frac{\|\frac{\alpha_{t}}{A_{t}}\nabla
f(x_{t+1})\|^{p+1}}{{\left(\frac{1}{2}\right)^{p-1}\frac{\bar{\kappa}_{p+1}^{t}}{{p}!}}}\right)^{\frac{1}{p}}$
$\displaystyle\geq\frac{f(x_{t+1})}{A_{t}},$

where the last inequality holds by our choice of $\bar{\kappa}_{p+1}^{t}$:

$\bar{\kappa}_{p+1}^{t}\geq\frac{(p+1)^{p+1}}{2}\frac{\alpha_{t}^{p+1}}{A_{t}}M.$
(64)

To sum up, by our choice of the parameters $\bar{\kappa}_{i}^{t}$,
$i=2,...,p$, we obtain

$\psi_{t+1}^{\ast}\geq\frac{f(x_{t+1})}{A_{t}}-err_{t+1}^{v}-err_{t+1}^{x}-err_{t+1}^{\tau}.$

∎

######  Theorem C.5.

Let Assumption 5.1 hold and $M\geq\frac{2}{p}L_{p}$. Let Assumption 5.4 hold.
After $T\geq 1$ with parameters defined in (14) and
$\sigma_{2}=\delta_{2}=\max\limits_{t=1,\ldots,T}\delta_{t}^{v_{t-1},x_{t}}$
we get the following bound for the objective residual

$\displaystyle\mathbb{E}\left[f(x_{T})-f(x^{\ast})\right]$
$\displaystyle\leq\frac{2(p+1)^{3}\tau
R}{(T+p+1)^{1/2}}+\frac{3(p+1)^{3}\sigma
R}{(T+p+1)^{1/2}}+\sum\limits_{i=3}^{p}\frac{2(p+1)^{2i-1}\delta_{i}R^{i}}{i!(T+p+1)^{i}}+\frac{(p+1)^{2(p+1)}}{(p+1)!}\frac{MR^{p+1}}{(T+p+1)^{p+1}}$

###### Proof.

First of all, let us bound $A_{T}$.

$\alpha_{t}=\frac{p+1}{t+p+1},~{}t\geq 1.$ (65)

Then, we have

$\frac{1}{(T+p+1)^{p+1}}\leq A_{T}\leq\frac{(p+1)!}{(T+1)^{p+1}}.$ (66)

And from [2] we get

$\sum_{t=0}^{T}\frac{A_{T}\alpha_{t}^{i}}{A_{t}}\leq\frac{(p+1)^{i}}{(T+p+1)^{i-1}}$
(67)

From Lemmas C.4 and C.1, we obtain that, for all $t\geq 1$,

$\displaystyle\frac{f(x_{t+1})}{A_{t}}-err_{t+1}^{v}-err_{t+1}^{x}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-err_{t+1}^{\tau}}\stackrel{{\scriptstyle\eqref{eq:lemma_ass_p}}}{{\leq}}\psi_{t+1}^{\ast}\leq\psi_{t+1}(x^{\ast})$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:acc_upper_bound_p}}}{{\leq}}\frac{f(x^{\ast})}{A_{t}}+\frac{\bar{\kappa}^{t}_{2}+\lambda_{t}}{2}\|x^{\ast}-x_{0}\|^{2}+\frac{\bar{\kappa}^{t}_{3}}{6}\|x^{\ast}-x_{0}\|^{3}+err_{t+1}^{up}.$

Next, we apply expectation

$\mathbb{E}\left[f(x_{T+1})-f(x^{\ast})\right]\leq
A_{T}\mathbb{E}\left[\frac{\bar{\kappa}_{2}^{T}+\lambda_{T}}{2}R^{2}+\sum\limits_{i=3}^{p+1}\frac{\bar{\kappa}_{i}^{T}}{i!}R^{i}+err_{T+1}^{up}+err_{T+1}^{v}+err_{T+1}^{x}+err_{T+1}^{\tau}\right].$
(68)

Let us choose

$\displaystyle\bar{\delta}_{t}=\delta_{2}+\frac{\tau+\sigma}{R}(t+p+1)^{3/2},$
(69) $\displaystyle\lambda_{t}=\frac{\sigma}{R}(t+p+1)^{p+1/2}.$ (70)

Then, we bound terms in (68) step by step.

We start from deterministic terms.

$\displaystyle
A_{T}\mathbb{E}\left[\frac{\bar{\kappa}_{2}^{T}+\lambda_{T}}{2}R^{2}+\sum\limits_{i=3}^{p}\frac{\bar{\kappa}_{i}^{T}}{i!}R^{i}+\frac{\bar{\kappa}_{p+1}^{T}}{(p+1)!}R^{p+1}\right]$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:p_ord_params}}}{{=}}p\alpha_{T}^{2}\bar{\delta}_{T}R^{2}+\sum\limits_{i=3}^{p}\frac{4p^{i-1}}{i!2}\alpha_{T}^{i}\delta_{i}R^{i}+\frac{(p+1)^{p+1}}{(p+1)!}\alpha_{T}^{p+1}MR^{p+1}$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:delta_p},\eqref{eq:lambda_p}}}{{=}}\frac{(p+1)^{3}(\tau+\sigma)R}{(T+p+1)^{1/2}}+\frac{(p+1)^{3}\delta_{2}R^{2}}{(T+p+1)^{2}}+\sum\limits_{i=3}^{p}\frac{2(p+1)^{2i-1}\delta_{i}R^{i}}{i!(T+p+1)^{i}}+\frac{(p+1)^{2(p+1)}}{(p+1)!}\frac{MR^{p+1}}{(T+p+1)^{p+1}}.$

Now, we bound expectation of all error terms. Firstly, we consider
$err_{T+1}^{up}$

$\displaystyle
A_{T}\mathbb{E}\left[err_{T+1}^{up}\right]=A_{T}\mathbb{E}\left[\sum\limits_{j=0}^{T}\frac{\alpha_{j}}{A_{j}}\left\langle
g(x_{j+1})-\nabla f(x_{j+1}),x^{\ast}-x_{j+1}\right\rangle\right]=0.$

Next, we bound $A_{T}\mathbb{E}\left[err_{T+1}^{v}\right]$

$\displaystyle
A_{T}\mathbb{E}\left[err_{T+1}^{v}\right]=A_{T}\mathbb{E}\left[\sum\limits_{j=0}^{T}\frac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}2}}{A_{j}\bar{\delta}_{j}}\|g(v_{j})-\nabla
f(v_{j})\|^{2}\right]\leq
2\sigma^{2}\sum\limits_{j=0}^{T}\frac{A_{T}}{A_{j}\bar{\delta}_{j}}$
$\displaystyle\stackrel{{\scriptstyle\eqref{alpha_t_p},\eqref{eq:delta_p}}}{{=}}\frac{2\sigma
R}{(p+1)^{3/2}}\sum\limits_{j=0}^{T}\frac{A_{T}\alpha_{j}^{3/2}}{A_{j}}\stackrel{{\scriptstyle\eqref{eq:alpha_sum_bound_p}}}{{\leq}}\frac{2\sigma
R}{(T+p+1)^{1/2}}$

Now we calculate $A_{T}\mathbb{E}\left[err_{T+1}^{x}\right]$

$\displaystyle
A_{T}\mathbb{E}\left[err_{T+1}^{x}\right]\stackrel{{\scriptstyle\eqref{eq:err_x}}}{{=}}A_{T}\mathbb{E}\left[\sum\limits_{j=0}^{T}\frac{\alpha_{j}^{2}}{2A_{j}^{2}\lambda_{j}}\|g(x_{j+1})-\nabla
f(x_{j+1})\|^{2}+\sum\limits_{j=0}^{T}\frac{\alpha_{j}}{A_{j}}\left\langle
g(x_{j+1})-\nabla f(x_{j+1}),y_{j}-x_{j+1}\right\rangle\right]$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:lambda_p}}}{{=}}\frac{\sigma
R}{2}\sum\limits_{j=0}^{T}\frac{A_{T}\alpha_{j}^{2}}{A_{j}^{2}(j+3)^{p+1/2}}\stackrel{{\scriptstyle\eqref{alpha_t}}}{{=}}\frac{\sigma
R}{(p+1)^{5/2}2}\sum\limits_{j=0}^{T}\frac{A_{T}\alpha_{j}^{[p+5/2}}{A_{j}^{2}}\leq\frac{\sigma
R}{(p+1)^{5/2}2A_{T}}\sum\limits_{j=0}^{T}\frac{A_{T}\alpha_{j}^{p+5/2}}{A_{j}}$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:alpha_sum_bound}}}{{\leq}}\frac{(p+1)^{2}\sigma
R}{2A_{T}(T+p+1)^{p+3/2}}\stackrel{{\scriptstyle\eqref{eq:A_t_bound}}}{{\leq}}\frac{(p+1)^{2}\sigma
R}{2(T+p+1)^{1/2}}.$

Finally, we consider $err_{T+1}^{\tau}$

$\displaystyle
A_{T}\mathbb{E}\left[err_{T+1}^{\tau}\right]=\sum\limits_{j=0}^{T}\frac{A_{T}\tau^{2}}{A_{j}\bar{\delta}_{j}}\stackrel{{\scriptstyle\eqref{eq:delta_p}}}{{=}}\frac{\tau
R}{(p+1)^{3/2}}\sum\limits\frac{A_{T}\alpha_{t}^{3/2}}{A_{j}}\stackrel{{\scriptstyle\eqref{eq:alpha_sum_bound}}}{{\leq}}\frac{\tau
R}{(T+p+1)^{1/2}}.$

Combining all bounds from above we achieve convergence rate

$\displaystyle\mathbb{E}\left[f(x_{T+1})-f(x^{\ast})\right]\leq\frac{2(p+1)^{3}\tau
R}{(T+p+1)^{1/2}}+\frac{3(p+1)^{3}\sigma
R}{(T+p+1)^{1/2}}+\sum\limits_{i=3}^{p}\frac{2(p+1)^{2i-1}\delta_{i}R^{i}}{i!(T+p+1)^{i}}+\frac{(p+1)^{2(p+1)}}{(p+1)!}\frac{MR^{p+1}}{(T+p+1)^{p+1}}.$

∎

Again, the case of stochastic Hessian (Theorem 5.6 under Assumption 5.3) can
be obtained in the same way by taking expectation in Lemma C.3.

##  Appendix D Restarts for strongly-convex function

######  Theorem 6.2.

Let Assumption 6.1 hold and let parameters of Algorithm 1 be chosen as in
(14). Let $\\{z_{s}\\}_{s\geq 0}$ be generated by Algorithm 3 and $R>0$ be
such that $\|z_{0}-x^{*}\|\leq R$. Then for any $s\geq 0$ we have

$\displaystyle\mathbb{E}\|z_{s}-x^{*}\|^{2}$ $\displaystyle\leq 4^{-s}R^{2},$
(71) $\displaystyle\mathbb{E}f(z_{s})-f(x^{*})$ $\displaystyle\leq
2^{-2s-1}\mu R^{2}.$ (72)

Moreover, the total number of iterations to reach desired accuracy
$\varepsilon:~{}f(z_{s})-f(x^{*})\leq\varepsilon$ in expectation is

$O\left(\tfrac{(\tau+\sigma_{1})^{2}}{\mu\varepsilon}+\left(\sqrt{\tfrac{\sigma_{2}}{\mu}}+1\right)\log\tfrac{f(z_{0})-f(x^{*})}{\varepsilon}+\sum\limits_{i=3}^{p}\left(\tfrac{\sigma_{i}R^{i-2}}{\mu}\right)^{\frac{1}{i}}+\left(\tfrac{L_{p}R^{p-1}}{\mu}\right)^{\frac{1}{p+1}}\right).$
(73)

###### Proof.

We prove by induction that $\mathbb{E}\|z_{s}-x^{*}\|^{2}\leq
4^{-s}\|z_{0}-x^{*}\|^{2}=4^{-s}R^{2}_{0}$. For $s=0$ this obviously holds. By
strong convexity and convergence of Algorithm 2

$\mathbb{E}\left[f(x_{T})-f(x^{\ast})\right]\leq\frac{C_{\tau}\tau
R}{\sqrt{T}}+\frac{C_{1}\sigma_{1}R}{\sqrt{T}}+\sum\limits_{i=2}^{p}\frac{C_{i}\sigma_{i}R^{i}}{T^{i}}+\frac{C_{p+1}L_{p}R^{p+1}}{T^{p+1}}$

we get

$\displaystyle\mathbb{E}_{[z_{s+1}|z_{s},z_{s-1},\ldots,z_{0}]}{\|z_{s+1}-x^{*}\|^{2}}\leq\frac{2}{\mu}\mathbb{E}_{[z_{s+1}|z_{s},\ldots,z_{0}]}{(f(x_{t_{s}+1})-f(x^{*}))}$
$\displaystyle\leq\frac{2}{\mu}\left(\frac{C_{\tau}\tau
R}{\sqrt{T}}+\frac{C_{1}\sigma_{1}R}{\sqrt{T}}+\sum\limits_{i=2}^{p}\frac{C_{i}\sigma_{i}R^{i}}{T^{i}}+\frac{C_{p+1}L_{p}R^{p+1}}{T^{p+1}}.\right)$
$\displaystyle\leq\frac{2}{\mu}\left(\frac{\mu
C_{\tau}\sigma\|{z_{s}}-x^{*}\|}{8(p+2)C_{\tau}\sigma\|z_{s}-x^{*}\|}+\frac{\mu
C_{1}\sigma\|{z_{s}}-x^{*}\|}{8(p+2)C_{1}\sigma\|z_{s}-x^{*}\|}\right.$
$\displaystyle+\sum\limits_{i=2}^{p}\frac{\mu
C_{i}\delta_{i}\|{z_{s}}-x^{*}\|^{i}}{8(p+2)C_{i}\delta_{i}\|z_{s}-x^{*}\|^{i}}\left.+\frac{\mu
C_{p+1}L_{p}\|{z_{s}}-x^{*}\|^{p+1}}{8(p+2)C_{p+1}L_{p}\|z_{s}-x^{*}\|^{p+1}}\right)$
$\displaystyle\leq\frac{p}{4(p+2)}\|z_{s}-x^{*}\|^{2}+\frac{2}{4(p+2)}r_{s}\|{z_{s}}-x^{*}\|.$

Then by taking full expectation we obtain

$\displaystyle\mathbb{E}{\|z_{s+1}-x^{*}\|^{2}}\leq\frac{p}{4(p+1)}\mathbb{E}\|z_{s}-x^{*}\|^{2}+\frac{2}{4(p+2)}r_{s}\mathbb{E}\|z_{s}-x^{*}\|\leq\frac{1}{4^{s+1}}R_{0}^{2}$

Thus, by induction, we obtain that (71), (72) hold.

Next we provide the corresponding complexity bounds. From the above induction
bounds, we obtain that after $S$ restarts the total number of iterations of
Algorithm 2

$\displaystyle\mathbb{E}T=\mathbb{E}\sum_{s=1}^{S}t_{s}\leq\sum_{s=1}^{S}\max\left\\{1,\left(\tfrac{8(p+2)C_{\tau}\tau}{\mu
r_{s-1}}\right)^{2},\left(\tfrac{8(p+2)C_{1}\sigma_{1}}{\mu
r_{s-1}}\right)^{2},\max\limits_{i=1,\ldots,p}\left(\tfrac{8(p+2)C_{i}\delta_{i}R_{s-1}^{i-2}}{\mu}\right)^{\frac{1}{i}},\left(\tfrac{8(p+2)C_{p+1}L_{p}R^{p-1}_{s-1}}{\mu}\right)^{\frac{1}{p+1}}\right\\}$
$\displaystyle\leq\sum\limits_{s=1}^{S}\left(1+\left(\tfrac{8(p+2)C_{\tau}\tau}{\mu
r_{s-1}}\right)^{2}+\left(\tfrac{8(p+2)C_{1}\sigma_{1}}{\mu
r_{s-1}}\right)^{2}+\textstyle{\sum\limits_{i=2}^{p}}\left(\tfrac{8(p+2)C_{i}\delta_{i}R_{s-1}^{i-2}}{\mu}\right)^{\frac{1}{i}}+\left(\tfrac{8(p+2)C_{p+1}L_{p}R^{p-1}_{s-1}}{\mu}\right)^{\frac{1}{p+1}}\right)$
$\displaystyle\leq S+\left(\tfrac{8(p+2)C_{\tau}\tau}{\mu
R_{0}}\right)^{2}4^{S}+\left(\tfrac{8(p+2)C_{1}\sigma_{1}}{\mu
R_{0}}\right)^{2}4^{S}+\left(\tfrac{8(p+2)C_{2}\delta_{2}}{\mu}\right)^{\frac{1}{2}}S+\textstyle{\sum\limits_{i=3}^{p}}\left(\tfrac{8(p+2)C_{i}\delta_{i}R_{0}^{i-2}}{\mu}\right)^{\frac{1}{i}}+\left(\tfrac{8(p+2)C_{p+1}L_{p}R^{p-1}_{0}}{\mu}\right)^{\frac{1}{p+1}}$
$\displaystyle\leq\log\frac{f(z_{0})-f(x^{*})}{\varepsilon}+\left(\tfrac{8(p+2)\tau}{\mu
R_{0}}\right)^{2}\tfrac{f(z_{0})-f(x^{*})}{\varepsilon}+\left(\tfrac{8(p+2)C_{1}\sigma_{1}}{\mu
R_{0}}\right)^{2}\tfrac{f(z_{0})-f(x^{*})}{\varepsilon}$
$\displaystyle+\left(\tfrac{8(p+2)C_{2}\delta_{2}}{\mu}\right)^{\frac{1}{2}}\log\tfrac{f(z_{0})-f(x^{*})}{\varepsilon}+\textstyle{\sum\limits_{i=3}^{p}}\left(\tfrac{8(p+2)C_{i}\delta_{i}R_{0}^{i-2}}{\mu}\right)^{\frac{1}{i}}+\left(\tfrac{8(p+2)L_{p}C_{p+1}R^{p-1}_{0}}{\mu}\right)^{\frac{1}{p+1}}.$

Therefore, the total oracle complexities are given by (73). ∎

For the case of stochastic high-order derivative the proof remains same w.r.t.
change $\delta_{i}$ to $\sigma_{i}$.