# A Multilevel Low-Rank Newton Method with Super-linear Convergence Rate and
its Application to Non-convex Problems

Nick Tsipinakis
Department of Mathematics
UniDistance Suisse
Brig, Switzerland
<EMAIL_ADDRESS>
&Panagiotis Tigkas
Department of Computer Science
Oxford University
Oxford, UK
<EMAIL_ADDRESS>
&Panos Parpas
Department of Computing
Imperial College London
London, UK
<EMAIL_ADDRESS>

###### Abstract

Second-order methods can address the shortcomings of first-order methods for
the optimization of large-scale machine learning models. However, second-order
methods have significantly higher computational costs associated with the
computation of second-order information. Subspace methods that are based on
randomization have addressed some of these computational costs as they compute
search directions in lower dimensions. Even though super-linear convergence
rates have been empirically observed, it has not been possible to rigorously
show that these variants of second-order methods can indeed achieve such fast
rates. Also, it is not clear whether subspace methods can be applied to non-
convex cases. To address these shortcomings, we develop a link between
multigrid optimization methods and low-rank Newton methods that enables us to
prove the super-linear rates of stochastic low-rank Newton methods rigorously.
Our method does not require any computations in the original model dimension.
We further propose a truncated version of the method that is capable of
solving high-dimensional non-convex problems. Preliminary numerical
experiments show that our method has a better escape rate from saddle points
compared to accelerated gradient descent and Adam and thus returns lower
training errors.

_Keywords_ multilevel methods $\cdot$ super-linear convergence $\cdot$ saddle-
free optimization $\cdot$ self-concordant functions

##  1 Introduction

When it comes to applying second-order optimization methods in machine
learning, there are two open questions: 1. Can second-order methods be
implemented efficiently? 2. Can second-order methods outperform standard
first-order methods for traditional ML metrics such as generalization error?
Regarding the latter question, several recent articles have argued the
efficacy of second-order methods in deep learning Singh and Alistarh (2020);
Pascanu and Bengio (2013); Dauphin et al. (2014), reinforcement learning Wu et
al. (2017) and variational inference Regier et al. (2017) (to name but a few).
Given the recent promising experimental results obtained via second-order
methods, the first question regarding the efficiency of these methods has
become even more urgent. In this paper we provide some answers to the first
question.

The development of randomized methods such as sketching and sub-sampling has
enabled the application of second-order methods to large scale ML models.
Super-linear convergence for sub-sampled methods has been shown in various
recent works e.g., Berahas et al. (2017); Bollapragada et al. (2019); Erdogdu
and Montanari (2015); Pilanci and Wainwright (2017a); Roosta-Khorasani and
Mahoney (2019). But all the existing methods require computations with the
full Hessian, which is of order $O(n^{2})$ to compute/store, and require
$O(n^{3})$ operations to compute the Newton direction. Stochastic multilevel
or sub-space methods, such as the one proposed in this paper, have $O(nN)$
cost for computing and storing the Hessian, and $O(n^{2}N)$ cost for computing
the Newton direction (where $N$ is the dimension of the sub-space). However,
it is unclear if stochastic multilevel methods can still converge super-
linearly when the hessian is approximated via a random method. Although it is
clear from numerical experiments that randomized Newton methods can obtain
super-linear convergence, rigorous proof under general conditions has not
appeared before Hanzely et al. (2020); Gower et al. (2019). In addition, to
the best of our knowledge there do not exist stochastic multilevel methods
which are capable of solving non-convex problems.

In terms of the theoretical contributions of this paper, the algorithms of
Pilanci and Wainwright (2017b) and Tsipinakis and Parpas (2021) are
particularly relevant since they are analyzed using self-concordant functions.
In Pilanci and Wainwright (2017b), restrictive assumptions are made regarding
the Hessian. In particular, the authors assume that the square root of the
Hessian is available or easily computable, and the method requires access to
the full Hessian and both the proof and hence the algorithm cannot be extended
to the non-convex case. In this paper we relax all these assumptions without
compromising the convergence rate. In Tsipinakis and Parpas (2021), super-
linear convergence was established, but the algorithm cannot be applied to the
non-convex case without non-trivial modifications. This paper links multilevel
optimization methods and randomized Newton methods. We exploit this link and
establish the global convergence of the algorithm and a local super-linear
convergence rate for self-concordant functions. Because of the way the low-
rank Hessian approximation is performed, our method can easily be extended to
the non-convex case and to very high dimensions (even when the Hessian is
dense).

Regarding the stochastic subspace methods, the work of Gower et al. (2019) is
quite relevant, however the method only achieves a linear rate, whereas we
achieve a super-linear rate for the same iteration complexity. In addition,
the method in Gower et al. (2019) cannot be extended to the non-convex case
without non-trivial modifications. In Hanzely et al. (2020), the authors
develop a related method that is based on the cubic Newton method. But the
method has a linear rate, and it is unclear if the method can be efficiently
applied to non-convex problems. Here, in order to apply the Newton method to
non-convex problems, we take a truncated SVD of the Hessian to keep the $N+1$
most informative eigenvalues and replace the rest with the $(N+1)^{\text{th}}$
eigenvalue while negative eigenvalues are replaced by their absolute value.
Therefore, as in the convex case, enforcing the Hessian matrix to be positive
definite and premultiplying the negative gradient by its inverse, we perform a
local change of coordinates and thus we should expect an accelerated
convergence behavior compared to the first-order methods. In addition, when
dealing with non-convex functions, the fact that we do not allow for
sufficiently small eigenvalues means that slow and flat manifolds around
saddle points will turn into saddles whose eigenvalues are large and we
therefore achieve a faster escape rate from the unstable region of the saddle
points. Similar approaches when computing the Newton direction have already
been explored O’Leary-Roseberry et al. (2019); Paternain et al. (2019);
Erdogdu and Montanari (2015). However, they all require the computation of the
full Hessian matrix and hence they are inefficient for high dimensional
problems. A work that scales to high dimensional problems which also
approximates the Hessian matrix was explored by Marteau-Ferey et al. (2019).
The difference to our approach is that it regularizes the Hessian matrix
instead of performing a Truncated SVD. However, selecting the regularization
parameter such that the Hessian matrix becomes positive definite is a
difficult task and thus the method is suitable to convex problems. In
addition, their method has limited applicability since it relies on the
existence of an $\ell_{2}$ regularizer to work. Such an assumption is not
present in this paper.

To this end, all the aforementioned advantages of our approach are
demonstrated in our numerical experiments. We show that our method is
extremely efficient in terms of computational requirements and it can be
applied to models with millions of parameters even when the Hessian is dense.
Our numerical experiments suggest that the proposed method behaves similarly
to the cubic Newton method (in terms of efficiency and ability to escape
saddle points) but with significantly less computational requirements.

##  2 Background and method

In this section, we briefly discuss the main components of multilevel
optimization methods, and introduce the coarse-grained model, which we later
modify to compute search directions. In the optimization literature,
multilevel methods are also referred to as multigrid methods. Since there is
no notion of ”grid" in the proposed algorithm we use the term multilevel, but
in the context of this paper the two terms are equivalent. Multilevel
optimization methods have two components: (1) A hierarchy of coarse models,
and (2) Linear prolongation and restriction operators used to “transfer"
information up and down the hierarchy of models. In this section, we describe
these two components; however, for simplicity, we consider a two-level
hierarchy. Our results can be extended to more than two levels without
substantial changes.

###  2.1 The Coarse-Grained Model and Main Assumptions

We first develop the proposed methodology for a convex optimization model. The
non-convex case is discussed in Section 3.1.1. Consider the following model,

$\mathbf{x}^{*}=\underset{\mathbf{x}\in{\mathbb{R}}^{n}}{\operatorname{arg\
min}}\ f(\mathbf{x}).$ (1)

where we assume that the function $f$ is a strictly convex and self-
concordant. Moreover, it has a closed sublevel set and is bounded below so
that $\mathbf{x}^{\star}$ exists. The definition and properties of self-
concordant functions can be found in Appendix A, and we refer the interested
reader to Boyd and Vandenberghe (2004); Nesterov (2004); Nesterov et al.
(2018); Nesterov and Nemirovskii (1994) for additional background information.
We will adopt the terminology used in the multilevel literature and call the
model in (1) the fine model. We will also assume that we can construct a lower
dimensional model called the coarse model. In traditional multilevel methods
the coarse model is typically derived by varying a discretization parameter.
In machine learning applications a coarse model can be derived by varying the
number of pixels in image processing applicationsGalun et al. (2015), or by
varying the dictionary size and fidelity in statistical pattern recognition
(see e.g. Hovhannisyan et al. (2016) for examples in face recognition, and
background extraction from video). We denote the coarse model as,

$\mathbf{y}^{*}=\underset{\mathbf{y}\in{\mathbb{R}}^{N}}{\operatorname{arg\
min}}\ F(\mathbf{y}).$ (2)

We assume that $N<n$ and typically $N\ll n$ and that $F$ is also a strictly
convex, self-concordant function. The property that $N<n$ justifies the use of
the terms coarse-grained or reduced order model. The role of the coarse model
is to help generate high quality search directions from the current incumbent
point $\mathbf{x}_{k}$ but at a reduced computational cost. Let
$\mathbf{R}_{k}\in{\mathbb{R}}^{N\times n}$ be the restriction operator with
which one may transfer information from the fine to coarse model and thus we
define the initial point in the coarse model as
$\mathbf{y_{0}}:=\mathbf{R}_{k}\mathbf{x}_{k}$ at iteration $k$. For
simplification purposes from now onwards we will omit the subscript $k$ from
the restriction operator. In order for the search direction to be useful (e.g.
a descent direction), we assume that first-order coherency condition below
holds,

$\mathbf{R}\nabla f(\mathbf{x_{k}})=\nabla F(\mathbf{y_{0}}),$ (3)

See Wen and Goldfarb (2009) for a discussion of the first-order coherency
condition in multilevel optimization. In addition to the restriction operator
we also assume that a prolongation operator,
$\mathbf{P}\in{\mathbb{R}}^{n\times N}$ is available. We make the following
assumption about these two operators.

######  Assumption 2.1.

The restriction and prolongation operators $\mathbf{R}$ and $\mathbf{P}$
satisfy, $\mathbf{P}=\sigma\mathbf{R}^{T},$ where $\sigma>0$, and $\mathbf{P}$
has full column rank, i.e., $\operatorname{rank}(\mathbf{P})=N$.

The assumptions regarding the restriction and prolongation operators are
standard Briggs et al. (2000) and, without loss of generality, we assume that
$\sigma=1$. A simple way to construct the prolongation and restriction
operators, which satisfies Assumption 2.1, arises from the naive
Nystr$\ddot{\text{o}}$m method (see Drineas and Mahoney (2005) for an
introduction). In particular, we construct $\mathbf{R}$ and $\mathbf{P}$ as
described in the definition below.

######  Definition 2.2.

Let $S_{n}=\left\\{1,2,\ldots,n\right\\}$ and denote $S_{N}\subset S_{n}$,
with the property that the elements $N<n$ are uniformly selected by the set
$S_{n}$ without replacement. Furthermore, assume that $s_{i}$ is the
$i^{\text{th}}$ element of $S_{N}$. Then the prolongation operator
$\mathbf{P}$ is generated as follows: the $i^{\text{th}}$ column of
$\mathbf{P}$ is the $s_{i}$ column of $\mathbf{I}_{n\times n}$ and,
furthermore, we set $\mathbf{R}$ as $\mathbf{P}^{T}$.

In addition to the first-order coherency condition, we will also assume that
the coarse model is also second-order coherent,

$\mathbf{R}\nabla
f^{2}(\mathbf{x_{k}})\mathbf{P}=\nabla^{2}F(\mathbf{y_{0}}).$ (4)

A simple and surprisingly effective method to construct the coarse model in
(2) while satisfying both the first and second order coherency conditions in
(3) and (4) is the so-called Galerkin model,

$\displaystyle F(\mathbf{y}):=\langle\mathbf{R}$ $\displaystyle\nabla
f(\mathbf{x_{k}}),\mathbf{y}-\mathbf{y}_{0}\rangle+$ (5)
$\displaystyle+\frac{1}{2}\langle\mathbf{R}\nabla^{2}f(\mathbf{x_{k}})\mathbf{P}(\mathbf{y}-\mathbf{y}_{0}),\mathbf{y}-\mathbf{y}_{0}\rangle.$

In the context of multilevel optimization methods, the Galerkin model was
experimentally tested in Gratton et al. (2010) where the authors found that it
compared favorably with other methods to construct coarse-grained models.
Motivated by the excellent numerical results in Gratton et al. (2010), and its
simplicity, we will adopt the model in (5) as the coarse-grained model in the
proposed algorithm. The model in (5) also has close links with the recently
proposed randomized Newton methods which we briefly discuss below (for more
details see Tsipinakis and Parpas (2021); Ho et al. (2019)).

We compute the coarse direction by
$\hat{\mathbf{d}}_{H,k}:=\mathbf{y}^{\star}-\mathbf{y}_{0}$. Furthermore, in
the case of the Galerkin model defined in (5) we have the closed form
solution,

$\begin{split}\hat{\mathbf{d}}_{H,k}&=\underset{\mathbf{d}_{H}\in\mathbb{R}^{N}}{\operatorname{arg\
min}}\left\\{\frac{1}{2}\|[\nabla^{2}F(\mathbf{y_{0}})]^{1/2}\mathbf{d}_{H}\|_{2}^{2}+\langle\nabla
F(\mathbf{y}_{0}),\mathbf{d}_{H}\rangle\right\\}\\\
&=-[\nabla^{2}F(\mathbf{y_{0}})]^{-1}\nabla F(\mathbf{y_{0}}).\end{split}$ (6)

Note that the coarse direction is a vector in ${\mathbb{R}}^{N}$. To correct
the incumbent solution $\mathbf{x}_{k}$ we must prolongate it to
${\mathbb{R}}^{n}$ using the prolongation operator,

$\hat{\mathbf{d}}_{h,k}:=\mathbf{P}\hat{\mathbf{d}}_{H,k}=-\mathbf{P}[\nabla^{2}F(\mathbf{y_{0}})]^{-1}\nabla
F(\mathbf{y}_{0}),$ (7)

where the $h$ and $H$ subscripts denote directions related to the fine and
coarse model respectively. If we naively set $\mathbf{P}$ as the identity
matrix, the search direction in (7) becomes the Newton direction. When
$\mathbf{P}$ is selected as in Definition 2.2 then we obtain the randomized
Newton method which is based on uniform sampling. In this case, the cost to
compute the reduced Hessian is $\mathcal{O}(nN)$ which is much cheaper than
$\mathcal{O}(n^{2})$ of the Newton method (for details on how to compute the
reduced Hessian matrix see Remark 2.1 in Tsipinakis and Parpas (2021)). The
general multilevel method of this section achieves a local super-linear
convergence for both self-concordant and strongly convex functions Tsipinakis
and Parpas (2021); Ho et al. (2019). However, when assuming self-concordant
function, one can additionally show that the general multilevel method enjoys
a global and scale invariant convergence analysis Tsipinakis and Parpas
(2021). In the next section, we develop a variant of the general multilevel
method which achieves a global analysis with a local super-linear convergence
rate, which is also applicable to non-convex functions.

##  3 Low-rank multilevel Newton methods

We saw above that the general multilevel method can be viewed as a randomized
Newton method. In this section we discuss connections of the multilevel scheme
with the low-rank Newton method through the naive Nyström method and propose
constructing the coarse model using a Truncated Singular Value Decomposition
(T-SVD). The version of our method which is suitable for non-convex
optimization problems is given in Algorithm 1.

Let $\mathbf{A}\in{\mathbb{R}}^{n\times n}$ be a positive definite matrix and
$\mathbf{Y}\in\mathbb{R}^{n\times N}$ with rank$(\mathbf{Y})=N<n$. The Nyström
method builds a $\operatorname{rank}$-$N$ approximation of $\mathbf{A}$,
namely $\mathbf{A}_{N}$, as follows,

$\mathbf{A}_{N}:=\mathbf{A}\mathbf{Y}(\mathbf{Y^{T}}\mathbf{A}\mathbf{Y})^{-1}(\mathbf{A}\mathbf{Y})^{T}\approx\mathbf{A}.$
(8)

The Nyström method has been extensively studied and has been shown to be an
efficient low-rank approximation method for different sampling techniques
which gives practitioners the freedom to choose from a wide range of random
matrices (for more details see Drineas and Mahoney (2005); Gittens (2011);
Smola and Schölkopf (2000); Williams and Seeger (2000)). Then, one may obtain
the naive Nyström method when $\mathbf{Y}$ is selected as in Definition 2.2.
Similarly, we may obtain a $\operatorname{rank}$-$N$ approximation of
$\mathbf{A}$ by T-SVD,

$\mathbf{A}_{N}:=\mathbf{U}_{N}\Sigma_{N}\mathbf{U}^{T}_{N}.$

For $\mathbf{A}$ being a symmetric positive definite matrix,
$\mathbf{\Sigma}_{N}\in{\mathbb{R}}^{N\times N}$ is a diagonal matrix
containing the $N$-largest eigenvalues of $\mathbf{A}$, with
$\sigma_{1}\geq\sigma_{2}\geq\cdots\geq\sigma_{N}$, and
$\mathbf{U}_{N}\in{\mathbb{R}}^{n\times N}$ contains the corresponding
eigenvectors. Although the naive Nyström method is more efficient compared to
the T-SVD method, for the algorithms we propose below we employ the T-SVD as
it gives us direct access to the eigenvalues of the Hessian matrix.

###  3.1 Low-rank Newton method with Convergence Analysis for Self-Concordant
Functions

Based on the Nyström method in (8), substituting $\mathbf{A}$ with
$\nabla^{2}f(\mathbf{x}_{k})$, $\mathbf{Y}$ with $\mathbf{P}$ and multiplying
right and left with $[\nabla^{2}f(\mathbf{x}_{k})]^{-1}$ we obtain,

$\displaystyle[\nabla^{2}f(\mathbf{x}_{k})]^{-1}\approx\mathbf{Q}_{h,k}^{-1}=\mathbf{P}(\mathbf{R}\nabla^{2}f(\mathbf{x}_{k})\mathbf{P})^{-1}\mathbf{R},$
(9)

which implies that the coarse direction in (7) can be viewed as an
approximation of the Newton direction that is based on a low-rank
approximation of the Hessian matrix. The fact that $\mathbf{Q}_{h,k}$ is a low
rank approximation of the Hessian matrix gives us the freedom to design search
directions using different low-rank approximation approaches. Here, similar to
the approach proposed in Erdogdu and Montanari (2015), we perform a T-SVD
approximation method (note that the work of Erdogdu and Montanari (2015) is
based on sub-sampling and thus it is different to ours and not directly
comparable). In particular, we construct $\mathbf{Q}_{h,k}^{-1}$ by computing
the $(N+1)^{\text{th}}$ T-SVD of $\nabla^{2}f(\mathbf{x}_{k})$ and then
replace all the eigenvalues after the $(N+1)^{\text{th}}$ eigenvalue with
$\sigma_{N+1}$. Formally, $\mathbf{Q}_{h,k}^{-1}$ is constructed as follows:

$\mathbf{Q}_{h,k}^{-1}:=\sigma_{N+1}^{-1}\mathbf{I}_{n\times
n}+\mathbf{U}_{N}(\mathbf{\Sigma}_{N}^{-1}-\sigma_{N+1}^{-1}\mathbf{I}_{N\times
N})\mathbf{U}_{N}^{T}.$ (10)

In other words, the matrix of eigenvalues of $\mathbf{Q}_{h,k}^{-1}$ has the
form
$\mathbf{\Sigma}^{-1}_{n}:=\operatorname{diag}\left(1/\sigma_{1},1/\sigma_{2},\ldots,1/\sigma_{N},1/\sigma_{N+1},\ldots,1/\sigma_{N+1}\right).$
Then, $\hat{\mathbf{d}}_{h,k}=-\mathbf{Q}_{h,k}^{-1}\nabla f(\mathbf{x}_{k})$
is a descent direction. As a result, for convex (stronlgy or self-concordant)
problems, this approach is efficient when the valuable second-order
information is concentrated on the first, say $N$, eigenvalues. In particular,
eigenvalues of large magnitude correspond to more informative directions since
in these directions the objective function has a large curvature, whereas for
eigenvalues that are close to zero the curvature becomes almost flat.
Therefore, by employing (10) in the computation of the search direction (7) we
aim to determine the subspace spanned by the dominant eigenvectors associated
with the largest eigenvalues, which are those that yield the fast convergence
rate of the Newton method. Such problem structures are typical in machine
learning problem Berahas et al. (2017), e.g., the Hessian matrix is (nearly)
low-rank and/or there is a big gap between the $\sigma_{N}$ and $\sigma_{N+1}$
eigenvalues. In the convergence analysis, we make use of the Newton decrement,
which is given by,

$\lambda(\mathbf{x}_{k}):=\left[\nabla
f(\mathbf{x}_{k})^{T}\nabla^{2}f(\mathbf{x}_{k})^{-1}\nabla
f(\mathbf{x}_{k})\right]^{1/2}.$

We are now in position to present the super-linear convergence of the above
scheme with its region depending on $\eta:=(3-\sqrt{9-4\varepsilon})/2$, where
$\varepsilon:=\min_{k\in\mathbb{N}}\\{\frac{\sigma_{k,n}}{\sigma_{k,N+1}}\\}$
and $\sigma_{k,j}$ is the $j^{\text{th}}$-eigenvalue of the Hessian matrix at
iteration $k$.

######  Theorem 3.1.

Suppose that the sequence $(\mathbf{x}_{k})_{k\in\mathbb{N}}$ is generated by
$\mathbf{x}_{k+1}=\mathbf{x}_{k}-t_{k}\mathbf{Q}_{h,k}^{-1}\nabla
f(\mathbf{x}_{k})$, with $\mathbf{Q}_{h,k}$ as in (10). Then, there exist
$\gamma>0$ and $\eta\in(0,\frac{3-\sqrt{5}}{2})$ such that

  1. (i)

if ${\lambda}(\mathbf{x}_{k})>\eta$, then
$f(\mathbf{x}_{k+1})-f(\mathbf{x}_{k})\leq-\gamma$

  2. (ii)

if ${\lambda}(\mathbf{x}_{k})\leq\eta$, then the line search selects the unit
step and

$\displaystyle\lambda(\mathbf{x}_{k+1})$
$\displaystyle<\lambda(\mathbf{x}_{k}),$ (11)

where in particular this process converges with a super-linear rate.

The proof of the theorem appears in Appendix A.2. As opposed to the classical
theory for strongly convex functions, the advantage of using self-concordance
as our main assumption results in obtaining a global analysis that has an
intuitive local super-linear convergence rate which only depends on the
difference between $\sigma_{n}$ and $\sigma_{N+1}$. Then, if the error
$\varepsilon$ is small, the low-rank Newton method is expected to converge
with a fast super-linear rate. Of course, the maximum range of the super-
linear convergence region is attained when $\sigma_{N+1}=\sigma_{n}$. In this
case, Theorem 3.1 is identical to that of the Newton method in Nesterov
(2004), and the method convergence quadratically. The cost per iteration of
the proposed method is $\mathcal{O}(Nn^{2})$. This is much better than
$\mathcal{O}(n^{3})$ of the full Newton method.

####  3.1.1 Extension to Non-Convex Problems

For self-concordant and other convex functions finding $\mathbf{x}^{*}$ can in
many instances be considered as an easy task (e.g., when descent methods are
applicable) since all local minima are also global. In this case, the unique
global minimum can be attained at an $\mathbf{x}$ for which $\|\nabla
f(\mathbf{x})\|=0$. On the other hand, for non-convex problems, finding
$\mathbf{x}^{*}$ is in general an NP hard problem. For the latter class of
problems, a point $\mathbf{x}$ for which $\|\nabla f(\mathbf{x})\|=0$ can be
one of the many local minima or a critical point such as a saddle, for which
the Hessian matrix is indefinite. Here, we are concerned with the task of
finding a local minimum of a general possible non-convex function $f$. To
achieve this, we propose a variant of the low-rank Newton method which we
conjecture will have a better escape rate from saddles compared to
(stochastic) first-order methods (our conjecture is validated by numerical
experiments).

Constructing the coarse direction using the definition in (10) is particularly
suitable for convex problems since all eigenvalues are positive. Then, for any
$N\in\\{1,2,\ldots,n-1\\}$ we have $\sigma_{N}>0$ which ensures the descent
property of $\hat{\mathbf{d}}_{h,k}$. However, when dealing with non-convex
functions such guarantee may not be true around saddle points. Since for the
Newton method the negative gradient is pre-multiplied by the inverse Hessian
matrix, a negative eigenvalue will result in changing the sign of the
corresponding gradient entry which may yield the Newton method converging to a
saddle point or a local maximum. On the other hand, the Newton method breaks
down when there are zero eigenvalues. To address these shortcomings, we
compute the $(N+1)^{\text{th}}$ T-SVD of the Hessian as above, but here we
replace all the negative eigenvalues by their absolute value. Further,
sufficiently small eigenvalues are replaced by a positive scalar to ensure the
non-singularity of the approximated Hessian matrix. Formally, we define
$g_{i}:{\mathbb{R}}^{N\times N}\rightarrow{\mathbb{R}},i=1,\ldots,N$,
$\text{for all}\ \ i\in\\{1,2,\dots,N\\}$, such that,

$g_{i}([\mathbf{\Sigma}_{N}]):=\max\\{|[\mathbf{\Sigma}_{N}]_{i}|,\nu\\},\
\nu>0,$ (12)

then we obtain the truncated low-rank approximation of the Hessian matrix as
follows,

$\displaystyle|{\mathbf{Q}}_{h,k}^{-1}|:=g(\sigma_{N+1})^{-1}$
$\displaystyle\mathbf{I}_{n\times
n}+\mathbf{U}_{N}(g(\mathbf{\Sigma}_{N})^{-1}-$ (13)
$\displaystyle-g(\sigma_{N+1})^{-1}\mathbf{I}_{N\times N})\mathbf{U}_{N}^{T},$

where
$g(\mathbf{\Sigma}_{N}):=\operatorname{diag}(g_{1}([\mathbf{\Sigma}_{N}]),\ldots,g_{N}([\mathbf{\Sigma}_{N}]))$
here is an element-wise mapping. Defining the approximated inverse Hessian
matrix as above is necessary in order to obtain a descent direction, which
together with an appropriately selected step-size parameter we can guarantee
the descent nature of the Low-Rank Newton method for non-convex problems.

Despite the lower iteration cost compared to the Newton method, there still
may be cases that forming the Hessian matrix and computing the T-SVD is too
expensive for some applications (when $n$ is extremely large). In the next
section, we address this issue.

###  3.2  Coarse-Grained Low-Rank Newton Method with Convergence Analysis for
Self-Concordant Functions

Algorithm 1 SigmaSVD

1: Input: $p,\ \nu,\ \epsilon\in(0,0.68^{2}),\alpha\in(0,0.5),\
\beta\in(0,1),\ \mathbf{P}_{k}\in\mathbb{R}^{n\times N}$,
$\mathbf{x}_{h,0}\in\mathbb{R}^{n}$

2: Compute $|{\mathbf{Q}}_{H,k}^{-1}|$ by (16) using the randomized T-SVD
Halko et al. (2011)

3: Form $|\hat{\mathbf{d}}_{h,k}|$ by (17)

4: Quit if $-\langle\nabla
f_{h,k}(\mathbf{x}_{h,k}),|\hat{\mathbf{d}}_{h,k}|\rangle\leq\epsilon\ \ $

5: Armijo search: while
$f_{h}(\mathbf{x}_{h,k}+t_{k}|\hat{\mathbf{d}}_{h,k}|)>f_{h}(\mathbf{x}_{h,k})+\alpha
t_{h,k}\nabla f^{T}_{h,k}(\mathbf{x}_{h,k})|\hat{\mathbf{d}}_{h,k}|,\
t_{h,k}\leftarrow\beta t_{h,k}$

6: Update:
$\mathbf{x}_{h,k+1}:=\mathbf{x}_{h,k}+t_{h,k}|\hat{\mathbf{d}}_{h,k}|$, go to
2

7: Return $\mathbf{x}_{h,k}$

The computational bottleneck of the procedure described in the previous
section arises from the fact that the computations are performed over the full
Hessian. To address this the we will perform the T-SVD on the Hessian matrix
of the reduced order model. This idea, along with the novel convergence
theory, is one of the main contributions of this paper. To keep the notation
simple, we do not specify connections with low-rank Newton methods formally as
before. However, one can show that the analysis in section 3.1 can be applied
to the reduced-order Hessian matrix and thus, as in equation (9), we may
obtain similar intuitive results for the reduced Hessian matrix. We begin by
describing a low-rank multilevel method which is suitable for convex
optimization and has a well-established theory for self-concordant functions.
Next, we present the extension for non-convex functions.

Let $\mathbf{Q}_{H,k}$ be a T-SVD-based $\operatorname{rank}$-$p$
approximation of the reduced Hessian matrix in (4), where $p<N<n$. Then, by
(8) we take

$\displaystyle\hat{\mathbf{Q}}_{h,k}^{-1}:=\mathbf{P}\mathbf{Q}_{H,k}^{-1}\mathbf{R}$
$\displaystyle\approx\mathbf{P}[\mathbf{R}\nabla^{2}f(\mathbf{x}_{k})\mathbf{P}]^{-1}\mathbf{R}$
$\displaystyle=\mathbf{Q}_{h,k}^{-1}\approx[\nabla^{2}f(\mathbf{x}_{k})]^{-1},$

which in fact is a $\operatorname{rank}$-$p$ approximation of the inverse
Hessian matrix with the difference that the computational cost of forming
$\hat{\mathbf{Q}}_{h,k}^{-1}$ is significantly reduced since we compute
$[\mathbf{U}_{p+1},\mathbf{\Sigma}_{p+1}]$ (by T-SVD) over the reduced Hessian
matrix (4). Thus, we formally define,

$\mathbf{Q}_{H,k}^{-1}:=\sigma_{p+1}^{-1}\mathbf{I}_{N\times
N}+\mathbf{U}_{p}(\mathbf{\Sigma}_{p}^{-1}-\sigma_{p+1}^{-1}\mathbf{I}_{p\times
p})\mathbf{U}_{p}^{T},$ (14)

where now $\mathbf{U}_{p}\in\mathbb{R}^{N\times p}$. Therefore, we can claim
that $\mathbf{\hat{d}}_{h,k}=-\hat{\mathbf{Q}}_{h,k}^{-1}\nabla
f(\mathbf{x}_{k})$ is an approximation of the Newton direction
$\mathbf{d}_{h,k}$. We call this method SigmaSVD as it is an extension of
Sigma in Tsipinakis and Parpas (2021) which constructs the reduced Hessian
matrix by T-SVD. Let us now define the approximate decrements for both
SigmaSVD and the general multilevel method of section 2,

$\begin{split}&\hat{\lambda}(\mathbf{x}_{k}):=\left[\nabla
f(\mathbf{x}_{k})^{T}\hat{\mathbf{Q}}_{h,k}^{-1}\nabla
f(\mathbf{x}_{k})\right]^{1/2},\\\
&\tilde{\lambda}(\mathbf{x}_{k}):=\left[(\mathbf{R}\nabla
f(\mathbf{x}_{k}))^{T}[\nabla^{2}F(\mathbf{y_{0}})]^{-1}\mathbf{R}\nabla
f(\mathbf{x}_{k})\right]^{1/2},\\\ \end{split}$ (15)

respectively. Both quantities are analogous to the Newton decrement
$\lambda(\mathbf{x}_{k})$ and they serve the same purpose. When assuming self-
concordant functions, the convergence analysis of the proposed method is
similar to that in Tsipinakis and Parpas (2021), with the difference that the
term that controls its convergence rate is given by
$\hat{e}_{k}:=(\lambda^{2}(\mathbf{x}_{k})-\hat{\lambda}^{2}(\mathbf{x}_{k}))^{1/2}$
instead of
$\tilde{e}_{k}:=(\lambda^{2}(\mathbf{x}_{k})-\tilde{\lambda}^{2}(\mathbf{x}_{k}))^{1/2}$.
It has been shown that
$\tilde{e}_{k}\leq(1-\mu_{k}^{2})^{1/2}\lambda(\mathbf{x}_{k})$,
$0<\mu_{k}\leq 1$ Tsipinakis and Parpas (2021). Intuitively, it should be
expected $\hat{e}_{k}$ to further incorporate the information carried out by
the T-SVD. Self-concordant functions allow us to prove such an informative
upper bound.

######  Lemma 3.2.

For all $k\in\mathbb{N}$ we have that

  1. 1.

$\mathbf{d}_{h,k}^{T}\nabla^{2}f(\mathbf{x}_{k})\hat{\mathbf{d}}_{h,k}=\hat{\lambda}^{2}(\mathbf{x}_{k})$

  2. 2.

$\hat{\mathbf{d}}_{h,k}^{T}\nabla^{2}f(\mathbf{x}_{k})\hat{\mathbf{d}}_{h,k}\leq\hat{\lambda}^{2}(\mathbf{x}_{k})$

  3. 3.

$\|[\nabla^{2}f(\mathbf{x}_{k})]^{\frac{1}{2}}(\mathbf{d}_{h,k}-\hat{\mathbf{d}}_{h,k})\|\leq\hat{e}_{k}$

  4. 4.

$\sqrt{\frac{\sigma_{k,N}}{\sigma_{k,p+1}}}\tilde{\lambda}(\mathbf{x}_{k})\leq\hat{\lambda}(\mathbf{x}_{k})\leq\tilde{\lambda}(\mathbf{x}_{k})\leq\lambda(\mathbf{x}_{k})$

Further, suppose that for all $k\in\mathbb{N}$ such that
$\mathbf{x}_{k}\neq\mathbf{x}^{*}$ it holds $\hat{\lambda}(\mathbf{x}_{k})>0$.
Then, there exists $\mu_{k}\in(0,1]$ such that
$\hat{e}_{k}\leq(1-\frac{\sigma_{k,N}}{\sigma_{k,p+1}}\mu_{k}^{2})^{1/2}\lambda(\mathbf{x}_{k})$

######  Definition 3.3.

Given $\mathbf{P}_{k}$ as in Definition 2.2, consider the random sequence
$M_{k}(\omega):=(\mathbf{x}_{k}(\omega),t_{k}(\omega),\hat{\mathbf{d}}_{k}(\omega),\mathbf{P}_{k}(\omega))$
that is generated by Algorithm 1 up to iteration $k$. We define
$\mathcal{F}_{k}:=\sigma(M_{0},M_{1},\ldots,M_{k})$ be the $\sigma$-algebra
generated by Algorithm 1.

Further, we define
$\hat{M}_{k}:=\mathbb{E}[\mu_{k}\sqrt{\frac{\sigma_{k,N}}{\sigma_{k,p+1}}}|\mathcal{F}_{k}]$,
$\hat{M}:=\min_{k\in\mathbb{N}}\\{\hat{M}_{k}\\}$ and
$\hat{\varepsilon}:=\sqrt{1-\hat{M}^{2}}$. Then, $0\leq\hat{\varepsilon}<1$,
and it can be shown that the proposed method converges super-linearly with its
region controlled by $\hat{\eta}:=\frac{3-\sqrt{5+4\hat{\varepsilon}}}{2}$.

######  Theorem 3.4.

Suppose that $\mathbf{P}$ is selected as in Definition 2.2. Suppose also that
the sequence $\\{\mathbf{x}_{k}\\}_{k\in\mathbb{N}}$ is generated by
$\mathbf{x}_{k+1}=\mathbf{x}_{k}-t_{k}\hat{\mathbf{Q}}_{h,k}^{-1}\nabla
f(\mathbf{x}_{k})$ and that the conditions of Lemma 3.2 are fulfilled. Then,
there exist constants $\hat{\gamma}>0$ and
$\hat{\eta}\in(0,\frac{3-\sqrt{5}}{2})$ such that

  1. (i)

if $\lambda(\mathbf{x}_{h,k})>\hat{\eta}$, then

$\mathbb{E}\left[f_{h}(\mathbf{x}_{h,k+1})-f_{h}(\mathbf{x}_{h,k})\right]\leq-\hat{\gamma},$

where $\gamma:=\mathbb{E}[\omega(\hat{\eta}\hat{M})]$ and
$\omega(x):=x-\log(1+x)$.

  2. (ii)

if $\lambda(\mathbf{x}_{h,k})\leq\hat{\eta}$, then

$\displaystyle\mathbb{E}[\lambda(\mathbf{x}_{h,k+1})]$
$\displaystyle<\mathbb{E}[\lambda(\mathbf{x}_{h,k})],$

and therefore
$\underset{k\rightarrow\infty}{\lim}\mathbb{E}[f_{h}(\mathbf{x}_{h,k})]=f_{h}(\mathbf{x}_{h}^{*})$
with probability one, where, in particular, this process converges with a
super-linear rate almost surely.

The proof of the theorem appears in Appendix A.3 and is analogue to that of
Tsipinakis and Parpas (2021) but now incorporates the results of Lemma 3.2.
Since now the analysis depends on $\hat{e}_{k}$ (instead of $\tilde{e}_{k}$),
one should expect the proposed method to converge to $\mathbf{x}^{*}$ with a
slower super-linear rate compared to that in Tsipinakis and Parpas (2021). In
particular, the convergence rate is expected to interpolate between a linear
and a super-linear rate of the general multilevel method. The conditions where
the general multilevel method converges with very fast super-linear or
quadratic rate to the global minimum were already discussed in Tsipinakis and
Parpas (2021). In our case, to achieve this, we additionally require that the
difference between $\sigma_{p}+1$ and $\sigma_{N}$ is sufficiently small. Even
in the case where randomness is induced by Definition 2.2 we are still able to
prove a global convergence analysis with a local super-linear convergence rate
that is independent of unknown constants.

We end this section by describing the convergence of SigmaSVD under different
assumptions. We provide two simple convergence results: (a) for strongly
convex functions with Lipschitz continuous gradients and (b) for non-convex
functions for which the Polyak-Lojasiewicz (PL) inequality and Lipschitz
continuity of the gradients hold (see Karimi et al. (2016) for a discussion on
the PL inequality) .

Based on the discussion in section 3.1.1 the truncated reduced Hessian matrix
is defined as follows,

$\displaystyle|{\mathbf{Q}}_{H,k}^{-1}|:=g(\sigma_{p+1})^{-1}$
$\displaystyle\mathbf{I}_{N\times
N}+\mathbf{U}_{p}(g(\mathbf{\Sigma}_{p})^{-1}$ (16)
$\displaystyle-g(\sigma_{p+1})^{-1}\mathbf{I}_{p\times p})\mathbf{U}_{p}^{T},$

then the truncated coarse direction is given by

$|\hat{\mathbf{d}}_{h,k}|:=-|\hat{\mathbf{Q}}_{h,k}^{-1}|\nabla
f(\mathbf{x}_{k})=-\mathbf{P}|{\mathbf{Q}}_{H,k}^{-1}|\mathbf{R}\nabla
f(\mathbf{x}_{k}).$ (17)

The full algorithm together with a step-size strategy is specified in
Algorithm 1. When the objective function is convex or self-concordant
Algorithm 1 coincides with SigmaSVD.

If the objective function is strongly convex with Lipschitz continuous
gradients then there exist $m,M>0$ (strong convexity and Lipschitz constants,
respectively) such that $m\mathbf{I}_{n}\leq\nabla^{2}f(\mathbf{x})\leq
M\mathbf{I}_{n}$ Boyd and Vandenberghe (2004).

######  Theorem 3.5.

Let $f(\mathbf{x})$ be a strongly convex function with Lipschitz continuous
gradients and that $\mathbf{x}_{k+1}$ generated by Algorithm 1. Then,

$\mathbb{E}[f(\mathbf{x}_{k+1})-f(\mathbf{x}_{k})]\leq-\frac{m}{2M}\mathbb{E}[\hat{M}_{k}\lambda(\mathbf{x}_{k})^{2}].$

See Appendix A.4 for a proof. Next, we drop the strong convexity assumption
and assume that the objective satisfies the PL inequality, i.e., there exist
$\xi>0$ such that

$\frac{1}{2}\|\nabla
f(\mathbf{x}_{k})\|^{2}\geq\xi(f(\mathbf{x}_{k})-f(\mathbf{x}^{*})).$ (18)

The PL condition is weaker than convexity and has been (empirically) shown
that it is often satisfied for over-parameterized neural networks (see Belkin
(2021) and references therein). We think that this result goes some way to
explain the behaviour of the algorithm (at least close enough to a minimum).
Using the PL inequality we obtain the following global result.

######  Theorem 3.6.

Assume that $f(\mathbf{x})$ has Lipschitz continuous gradients and that the PL
inequality holds. Assume further that $\mathbf{x}_{k+1}$ is generated by
Algorithm 1. Then,

$\mathbb{E}[f(\mathbf{x}_{k+1})-f(\mathbf{x}_{k})]\leq-\frac{\xi}{M}\mathbb{E}\left[\frac{\sigma_{k,p+1}^{2}}{\sigma_{k,1}^{2}}(f(\mathbf{x}_{k})-f(\mathbf{x}^{*}))\right],$

where $\sigma_{k,p+1}$ and $\sigma_{k,1}$ are the $(p+1)^{\text{th}}$ and
maximum eigenvalues of the truncated Hessian matrix.

See Appendix A.5 for a proof. Note that the above result is pessimistic
compared to the one in Theorem 3.5 because it does not involve the approximate
and the Newton decrements. Nevertheless, it shows that with a weaker
assumption than convexity we are still able to obtain global convergence
results. The advantages of using a truncated version of the descent direction
when dealing with non-convex functions have been already discussed in Section
3.1.1. A similar behavior should be also expected by the truncated multilevel
method. Nevertheless, the method of this section requires
$\mathcal{O}(nN+pN^{2})$ to form the reduced Hessian and compute its truncated
inverse which is a significant improvement compared to the low-rank Newton
method. This number is also better than the general multilevel method, and
thus SigmaSVD should be faster when the difference between the eigenvalues is
sufficiently small. Moreover, this number can even be reduced to
$\mathcal{O}(N+pN^{2})$ when performing computations in parallel. In non-
convex problems, it is expected that only the first few eigenvalues
concentrate the important second-order information; hence, in practice, both
$p$ and $N$ should be kept small Alain et al. (2019); Ghorbani et al. (2019).
This results in a method with low per-iteration cost, which is expected to
have a much faster escape rate from saddle points compared to first-order
methods. This is illustrated in the next section.

##  4 Numerical Results

(a) Gisette

(b) Gisette

(c) Gissete with $\mathbf{x}_{0}=0$

Figure 1: Non-convex minimization. All the methods in plot (a) are initialized
at the origin, while in plot (b) the initializer is selected randomly by
$\mathcal{N}(0,1)$. Plot (c) shows the convergence behavior of SigmaSVD for
different values of $p$.

We are now ready to validate the efficiency of the proposed method on
different machine learning models. We illustrate that the method is efficient
and compares favorably with other state-of-the-art algorithms on a broad class
of problems even when the assumptions imposed in this paper do not hold. In
particular, we consider minimizing a self-concordant function, strongly convex
functions with or without Lipschitz continuity and a non-convex function. We
also illustrate that the method is efficient on a very large problem with more
than a million dimensions. Due to space limitations in this section we only
report the results for the non-convex case. The additional numerical results
together with a detailed description of the algorithms, training datasets, and
the setup used to obtain the results appear in the Appendix B.

Suppose that we are given a training set
$\\{\mathbf{a}_{i},b_{i}\\}_{i=1}^{m}$ with
$\mathbf{a}_{i}\in{\mathbb{R}}^{n}$ and $b_{i}\in[0,1]$ and consider solving
the non-linear least-squares problem, i.e.,
$\min_{\mathbf{x}\in{\mathbb{R}}^{n}}\frac{1}{m}\sum_{i=1}^{m}(b_{i}-\phi(\mathbf{a}_{i}^{T}\mathbf{x}))^{2}$,
where $\phi(w):=1/(1+\exp(w)^{-1})$ is the sigmoid function. Since we apply
the mean squared error over the sigmoid function the non-linear least-squares
is a non-convex problem. We compare the performance of algorithm 1 against
gradient descent (GD), accelerated gradient descent (AGD), cubic Newton and
Adam. Although Adam, as a batch algorithm, is not directly comparable to
SigmaSVD, we include it in our comparisons to demonstrate the advantages of
second-order methods compared to the stochastic first-order method on problems
with several saddle points and flat areas. For algorithm 1, GD and AGD the
Armijo rule is applied to determine the step size parameter with constants
$\alpha=0.001,\beta=0.7$. A line search strategy is applied to select the
regularization parameter of the Cubic Newton method (for details see Nesterov
et al. (2018)). In all experiments we select the fixed eigenvalue threshold,
i.e., $\nu=10^{-10}$ in (12). The momentum parameters for Adam are selected as
suggested in Kingma and Ba (2014) while for AGD the momentum is selected
$0.5$.

In Figure 1 we demonstrate the performance of the different optimization
methods on the non-linear least-squares problem for the Gisette dataset (see
also Figure 2 in Appendix B for simulations on different data regimes). In
Figure 1(a), 1(b) we illustrate the reduction of the value function over GPU
time. Observe that SigmaSVD is able to return a better solution compared to
the first-order methods. Clearly, first-order methods are trapped in a saddle
point or flat area and thus they are stuck far from the local or global
minimum. This is not a surprise since in such areas the gradient is almost
zero and thus first-order methods are unable to progress and although in
theory they have been shown to always escape even ill-conditioned saddles,
here we see that they require an extremely large amount of iterations which
makes them inefficient for practical applications. On the other hand, Cubic
Newton is able to escape saddle points in one iterations however we observed
that in many experiments it converges rapidly to early local minima probably
due to the lack of randomness. Moreover, as suggested in Theorem 3.4, Figure 1
illustrates that SigmaSVD enjoys a fast local convergence rate and it is
faster than its counterparts. It is also faster and has the ability to achieve
lower training errors for different initialization points.

In Table 1 and Figure 1(c) we demonstrate how the choice of $N$ and $p$
affects the behavior of SigmaSVD around a saddle point or flat area. We
revisit the experiment of Figure 1(b) (i.e., $\mathbf{x}_{0}=0$) in which due
to the presence of a saddle point first order methods are trapped far from the
global minimum while SigmaSVD escapes from such a region in one iteration
(similar to Cubic Newton). Figure 1(c) shows the behavior of SigmaSVD for
different values of $p$ with fixed $N=0.5n$. Observe that for very small
values of $p$ SigmaSVD is trapped in the same saddle as the first-order
however for $p\geq 350$ it converges to the global minimum. Next, in Table 1
we fix $p=450$ and vary $N$ to show the escape rate probability from the
saddle, or otherwise the probability of convergence to the global minimum over
$50$ trials. Similarly to $p$, we see that the escape rate probability is
proportional to $N$. Both experiments verify Theorem 3.4 which indicates that
for very small values in $p$ and $N$, one should expect SigmaSVD to behave as
a first-order method. Furthermore, Table 1 and Figure 1 show that SigmaSVD is
able to reach the behavior of the Cubic Newton method (that is, escape from
the saddle point in one iteration) using only half of the dimensions of the
problem and about $20\%$ of the eigenvalues of the reduced Hessian, which
constitutes a significant improvement in total complexity as also depicted in
Figure 1(a).

Table 1: Probability of successfully escaping from saddle points and
convergence to the global minimum for various values of $N$ and fixed
$p=450,\mathbf{x}_{0}=0$. Each row shows the probability over $50$ trials.
Escape Rate Probability - Gisette
---
N | Probability
$0.1n$ | $18\%$
$0.13n$ | $46\%$
$0.26n$ | $52\%$
$0.36n$ | $66\%$
$0.42n$ | $80\%$
$0.46n$ | $92\%$

##  5 Conclusion

We develop a stochastic multilevel low-rank Newton type method that enjoys a
super-linear convergence rate for self-concordant functions with low per-
iteration cost. We further propose a variant that is applicable to non-convex
problems. Preliminary numerical experiments show that our method is efficient
for large-scale optimization problems with millions of dimensions. It is also
faster and has an improved escape rate from saddles and flat-areas compared to
first-order methods in non-convex cases. As a future direction, we aim to
develop a batch variant of our method that is efficient for training deep
neural networks. We also plan to provide a convergence analysis for non-convex
functions.

## References

  * Singh and Alistarh [2020] Sidak Pal Singh and Dan Alistarh.  Woodfisher: Efficient second-order approximation for neural network compression.  _Advances in Neural Information Processing Systems_ , 33:18098–18109, 2020.
  * Pascanu and Bengio [2013] Razvan Pascanu and Yoshua Bengio.  Revisiting natural gradient for deep networks.  _arXiv preprint arXiv:1301.3584_ , 2013.
  * Dauphin et al. [2014] Yann N Dauphin, Razvan Pascanu, Caglar Gulcehre, Kyunghyun Cho, Surya Ganguli, and Yoshua Bengio.  Identifying and attacking the saddle point problem in high-dimensional non-convex optimization.  _Advances in neural information processing systems_ , 27, 2014.
  * Wu et al. [2017] Yuhuai Wu, Elman Mansimov, Roger B Grosse, Shun Liao, and Jimmy Ba.  Second-order optimization for deep reinforcement learning using kronecker-factored approximation.  In _NIPS_ , 2017.
  * Regier et al. [2017] Jeffrey Regier, Michael I Jordan, and Jon McAuliffe.  Fast black-box variational inference through stochastic trust-region optimization.  _Advances in Neural Information Processing Systems_ , 30, 2017.
  * Berahas et al. [2017] Albert S Berahas, Raghu Bollapragada, and Jorge Nocedal.  An investigation of newton-sketch and subsampled newton methods.  _arXiv preprint arXiv:1705.06211_ , 2017.
  * Bollapragada et al. [2019] Raghu Bollapragada, Richard H Byrd, and Jorge Nocedal.  Exact and inexact subsampled newton methods for optimization.  _IMA Journal of Numerical Analysis_ , 39(2):545–578, 2019.
  * Erdogdu and Montanari [2015] Murat A Erdogdu and Andrea Montanari.  Convergence rates of sub-sampled newton methods.  In _Proceedings of the 28th International Conference on Neural Information Processing Systems-Volume 2_ , pages 3052–3060. MIT Press, 2015.
  * Pilanci and Wainwright [2017a] Mert Pilanci and Martin J. Wainwright.  Newton sketch: a near linear-time optimization algorithm with linear-quadratic convergence.  _SIAM J. Optim._ , 27(1):205–245, 2017a.  ISSN 1052-6234.  doi:10.1137/15M1021106.  URL https://doi.org/10.1137/15M1021106.
  * Roosta-Khorasani and Mahoney [2019] Farbod Roosta-Khorasani and Michael W Mahoney.  Sub-sampled newton methods.  _Mathematical Programming_ , 174(1):293–326, 2019\.
  * Hanzely et al. [2020] Filip Hanzely, Nikita Doikov, Yurii Nesterov, and Peter Richtarik.  Stochastic subspace cubic newton method.  In _International Conference on Machine Learning_ , pages 4027–4038. PMLR, 2020.
  * Gower et al. [2019] Robert Gower, Dmitry Kovalev, Felix Lieder, and Peter Richtárik.  Rsn: Randomized subspace newton.  _Advances in Neural Information Processing Systems_ , 32, 2019.
  * Pilanci and Wainwright [2017b] Mert Pilanci and Martin J Wainwright.  Newton sketch: A near linear-time optimization algorithm with linear-quadratic convergence.  _SIAM Journal on Optimization_ , 27(1):205–245, 2017b.
  * Tsipinakis and Parpas [2021] Nick Tsipinakis and Panos Parpas.  A multilevel method for self-concordant minimization.  _arXiv preprint arXiv:2106.13690_ , 2021.
  * O’Leary-Roseberry et al. [2019] Thomas O’Leary-Roseberry, Nick Alger, and Omar Ghattas.  Inexact newton methods for stochastic nonconvex optimization with applications to neural network training.  _arXiv preprint arXiv:1905.06738_ , 2019.
  * Paternain et al. [2019] Santiago Paternain, Aryan Mokhtari, and Alejandro Ribeiro.  A newton-based method for nonconvex optimization with fast evasion of saddle points.  _SIAM Journal on Optimization_ , 29(1):343–368, 2019.
  * Marteau-Ferey et al. [2019] Ulysse Marteau-Ferey, Francis Bach, and Alessandro Rudi.  Globally convergent newton methods for ill-conditioned generalized self-concordant losses.  _Advances in Neural Information Processing Systems_ , 32, 2019.
  * Boyd and Vandenberghe [2004] Stephen Boyd and Lieven Vandenberghe.  _Convex optimization_.  Cambridge University Press, Cambridge, 2004.  ISBN 0-521-83378-7.  doi:10.1017/CBO9780511804441.  URL https://doi.org/10.1017/CBO9780511804441.
  * Nesterov [2004] Yurii Nesterov.  _Introductory lectures on convex optimization_ , volume 87 of _Applied Optimization_.  Kluwer Academic Publishers, Boston, MA, 2004.  ISBN 1-4020-7553-7.  doi:10.1007/978-1-4419-8853-9.  URL https://doi.org/10.1007/978-1-4419-8853-9.  A basic course.
  * Nesterov et al. [2018] Yurii Nesterov et al.  _Lectures on convex optimization_ , volume 137.  Springer, 2018.
  * Nesterov and Nemirovskii [1994] Yurii Nesterov and Arkadii Nemirovskii.  _Interior-point polynomial algorithms in convex programming_ , volume 13 of _SIAM Studies in Applied Mathematics_.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1994.  ISBN 0-89871-319-6.  doi:10.1137/1.9781611970791.  URL https://doi.org/10.1137/1.9781611970791.
  * Galun et al. [2015] Meirav Galun, Ronen Basri, and Irad Yavneh.  Review of methods inspired by algebraic-multigrid for data and image analysis applications.  _Numer. Math. Theory Methods Appl._ , 8(2):283–312, 2015.  ISSN 1004-8979.  doi:10.4208/nmtma.2015.w14si.  URL https://doi.org/10.4208/nmtma.2015.w14si.
  * Hovhannisyan et al. [2016] Vahan Hovhannisyan, Panos Parpas, and Stefanos Zafeiriou.  MAGMA: multilevel accelerated gradient mirror descent algorithm for large-scale convex composite minimization.  _SIAM J. Imaging Sci._ , 9(4):1829–1857, 2016\.  ISSN 1936-4954.  doi:10.1137/15M104013X.  URL https://doi.org/10.1137/15M104013X.
  * Wen and Goldfarb [2009] Zaiwen Wen and Donald Goldfarb.  A line search multigrid method for large-scale nonlinear optimization.  _SIAM J. Optim._ , 20(3):1478–1503, 2009.  ISSN 1052-6234.  doi:10.1137/08071524X.  URL https://doi.org/10.1137/08071524X.
  * Briggs et al. [2000] William L. Briggs, Van Emden Henson, and Steve F. McCormick.  _A multigrid tutorial_.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, second edition, 2000.  ISBN 0-89871-462-1.  doi:10.1137/1.9780898719505.  URL https://doi.org/10.1137/1.9780898719505.
  * Drineas and Mahoney [2005] Petros Drineas and Michael W. Mahoney.  On the Nyström method for approximating a Gram matrix for improved kernel-based learning.  _J. Mach. Learn. Res._ , 6:2153–2175, 2005.  ISSN 1532-4435.
  * Gratton et al. [2010] Serge Gratton, Mélodie Mouffe, Annick Sartenaer, Philippe L. Toint, and Dimitri Tomanos.  Numerical experience with a recursive trust-region method for multilevel nonlinear bound-constrained optimization.  _Optim. Methods Softw._ , 25(3):359–386, 2010\.  ISSN 1055-6788.  doi:10.1080/10556780903239295.  URL https://doi.org/10.1080/10556780903239295.
  * Ho et al. [2019] Chin Pang Ho, Michal Kočvara, and Panos Parpas.  Newton-type multilevel optimization method.  _Optimization Methods and Software_ , pages 1–34, 2019.
  * Gittens [2011] Alex Gittens.  The spectral norm error of the naive nystrom extension.  _arXiv preprint arXiv:1110.5305_ , 2011.
  * Smola and Schölkopf [2000] AJ Smola and B Schölkopf.  Sparse greedy matrix approximation for machine learning.  In _Seventeenth International Conference on Machine Learning (ICML 2000)_ , pages 911–918. Morgan Kaufmann, 2000.
  * Williams and Seeger [2000] Christopher Williams and Matthias Seeger.  Using the nyström method to speed up kernel machines.  _Advances in neural information processing systems_ , 13, 2000.
  * Halko et al. [2011] N. Halko, P. G. Martinsson, and J. A. Tropp.  Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions.  _SIAM Rev._ , 53(2):217–288, 2011.  ISSN 0036-1445.  doi:10.1137/090771806.  URL https://doi.org/10.1137/090771806.
  * Karimi et al. [2016] Hamed Karimi, Julie Nutini, and Mark Schmidt.  Linear convergence of gradient and proximal-gradient methods under the polyak-łojasiewicz condition.  In _Machine Learning and Knowledge Discovery in Databases: European Conference, ECML PKDD 2016, Riva del Garda, Italy, September 19-23, 2016, Proceedings, Part I 16_ , pages 795–811. Springer, 2016.
  * Belkin [2021] Mikhail Belkin.  Fit without fear: remarkable mathematical phenomena of deep learning through the prism of interpolation.  _Acta Numerica_ , 30:203–248, 2021.
  * Alain et al. [2019] Guillaume Alain, Nicolas Le Roux, and Pierre-Antoine Manzagol.  Negative eigenvalues of the hessian in deep neural networks.  _arXiv preprint arXiv:1902.02366_ , 2019.
  * Ghorbani et al. [2019] Behrooz Ghorbani, Shankar Krishnan, and Ying Xiao.  An investigation into neural net optimization via hessian eigenvalue density.  In _International Conference on Machine Learning_ , pages 2232–2241. PMLR, 2019.
  * Kingma and Ba [2014] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Horn and Johnson [2013] Roger A. Horn and Charles R. Johnson.  _Matrix analysis_.  Cambridge University Press, Cambridge, second edition, 2013.  ISBN 978-0-521-54823-6.

##  Appendix A Proofs of theorems and background

###  A.1 Background

Throughout this paper, all vectors are denoted with bold lowercase letters,
i.e., $\mathbf{x}\in\mathbb{R}^{n}$ and all matrices with bold uppercase
letters, i.e., $\mathbf{A}\in\mathbb{R}^{m\times n}$. The function
$\|\mathbf{x}\|_{2}=\langle\mathbf{x}^{T},\mathbf{x}\rangle^{1/2}$ is the
$\ell_{2}$\- or Euclidean norm of $\mathbf{x}$. The spectral norm of
$\mathbf{A}$ is the norm induced by the $\ell_{2}$-norm on $\mathbb{R}^{n}$
and it is defined as
$\|\mathbf{A}\|_{2}:=\underset{\|\mathbf{x}\|_{2}=1}{\max}\|\mathbf{A}\mathbf{x}\|_{2}$.
It can be shown that $\|\mathbf{A}\|_{2}=\sigma_{1}(\mathbf{A})$, where
$\sigma_{1}(\mathbf{A})$ (or simply $\sigma_{1}$) is the largest singular
value of $\mathbf{A}$, see [Horn and Johnson, 2013, Section 5.6]. For two
symmetric matrices $\mathbf{A}$ and $\mathbf{B}$ we write
$\mathbf{A}\succeq\mathbf{B}$ when
$\mathbf{x}^{T}(\mathbf{A}-\mathbf{B})\mathbf{x}\geq 0$ for all
$\mathbf{x}\in\mathbb{R}^{n}$, or otherwise when the matrix
$\mathbf{A}-\mathbf{B}$ is positive semi-definite. Below we present the main
properties and inequalities for self-concordant functions. For a complete
analysis see Nesterov [2004], Boyd and Vandenberghe [2004]. A univariate
convex function $\phi:\mathbb{R}\rightarrow\mathbb{R}$ is called self-
concordant with constant $M_{\phi}\geq 0$ if,

$|\phi^{\prime\prime\prime}(x)|\leq M_{\phi}\phi^{\prime\prime}(x)^{3/2}.$
(19)

Examples of self-concordant functions include but not limited to linear,
convex quadratic, negative logarithmic and negative log-determinant function.
Based on the above definition, a multivariate convex function
$f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is called self-concordant if
$\phi(t):=f(\mathbf{x}+t\mathbf{u})$ satisfies (19) for all
$\mathbf{x}\in\operatorname{dom}f,\mathbf{u}\in\mathbb{R}^{n}$ and
$t\in\mathbb{R}$ such that $\mathbf{x}+t\mathbf{u}\in\operatorname{dom}f$.
Further, self-concordance is preserved under composition with any affine
function. In addition, for any convex self-concordant function $\tilde{\phi}$
with constant $M_{\tilde{\phi}}\geq 0$ it can be shown that
$\phi(x):=\frac{M_{\tilde{\phi}}^{2}}{4}\tilde{\phi}(x)$ is self-concordant
with constant $M_{\phi}=2$. Next, given $\mathbf{x}\in\operatorname{dom}f$ and
assuming that $\nabla^{2}f(\mathbf{x})$ is positive-definite we can define the
following norms

$\|\mathbf{u}\|_{\mathbf{x}}:=\langle\nabla^{2}f(\mathbf{x})\mathbf{u},\mathbf{u}\rangle^{1/2}\
\ \text{and}\ \
\|\mathbf{v}\|_{\mathbf{x}}^{*}:=\langle[\nabla^{2}f(\mathbf{x})]^{-1}\mathbf{v},\mathbf{v}\rangle^{1/2},$
(20)

for which it holds that
$|\langle\mathbf{u},\mathbf{v}\rangle|\leq\|\mathbf{u}\|_{\mathbf{x}}^{*}\|\mathbf{v}\|_{\mathbf{x}}$.
Therefore, the Newton decrement can be written as

$\lambda_{f}(\mathbf{x}):=\|\nabla
f(\mathbf{x})\|_{\mathbf{x}}^{*}=\|[\nabla^{2}f(\mathbf{x})]^{-1/2}\nabla
f(\mathbf{x})\|_{2}.$ (21)

In addition, we take into consideration two auxiliary functions, both
introduced in Nesterov [2004]. Define the functions $\omega$ and $\omega_{*}$
such that

$\omega(x):=x-\log(1+x)\quad\text{and}\quad\omega_{*}(x):=-x-\log(1-x),$ (22)

where $\operatorname{dom}\omega=\\{x\in\mathbb{R}:x\geq 0\\}$ and
$\operatorname{dom}\omega_{*}=\\{x\in\mathbb{R}:0\leq x<1\\}$, respectively,
and $\log(x)$ denotes the natural logarithm of $x$. Moreover, note that both
functions are convex and their range is the set of positive real numbers.
Further, from the definition (19) for $M_{\phi}=2$, we have that

$\left|\frac{d}{dt}\left(\phi^{\prime\prime}(t)^{-1/2}\right)\right|\leq 1,$

from which, after integration, we obtain the following bounds

$\frac{\phi^{\prime\prime}(0)}{(1+t\phi^{\prime\prime}(0)^{1/2})^{2}}\leq\phi^{\prime\prime}(t)\leq\frac{\phi^{\prime\prime}(0)}{(1-t\phi^{\prime\prime}(0)^{1/2})^{2}}$
(23)

where the lower and the uppers bounds hold for $t\geq 0$ and
$t\in[0,\phi^{\prime\prime}(0)^{-1/2})$, with $t\in\operatorname{dom}\phi$,
respectively (see also Boyd and Vandenberghe [2004]). Consider now self-
concordant functions on $\mathbb{R}^{n}$ and let
$S(\mathbf{x})=\left\\{\mathbf{y}\in\mathbb{R}^{n}:\|\mathbf{y}-\mathbf{x}\|_{\mathbf{x}}<1\right\\}$.
For any $\mathbf{x}\in\operatorname{dom}f$ and $\mathbf{y}\in S(\mathbf{x})$,
we have that (Nesterov [2004])

$(1-\|\mathbf{y}-\mathbf{x}\|_{\mathbf{x}})^{2}\nabla^{2}f(\mathbf{x})\preceq\nabla^{2}f(\mathbf{y})\preceq\frac{1}{(1-\|\mathbf{y}-\mathbf{x}\|_{\mathbf{x}})^{2}}\nabla^{2}f(\mathbf{x}).$
(24)

The following results will be used later for the convergence analysis of the
algorithms. For their proofs see Nesterov [2004].

######  Lemma A.1 (Nesterov [2004]).

Let $\mathbf{x},\mathbf{y}\in\operatorname{dom}f$. If
$\|\mathbf{y}-\mathbf{x}\|_{x}<1$, then,

$f(\mathbf{y})\leq f(\mathbf{x})+\langle\nabla
f(\mathbf{x}),\mathbf{y}-\mathbf{x}\rangle+\omega_{*}(\|\mathbf{y}-\mathbf{x}\|_{\mathbf{x}})$

######  Lemma A.2 (Nesterov [2004]).

Let $\mathbf{x},\mathbf{y}\in\operatorname{dom}f$. Then,

$f(\mathbf{y})\geq f(\mathbf{x})+\langle\nabla
f(\mathbf{x}),\mathbf{y}-\mathbf{x}\rangle+\omega(\|\nabla
f(\mathbf{y})-\nabla f(\mathbf{x})\|_{\mathbf{y}}^{*}).$

If in addition $\|\nabla f(\mathbf{y})-\nabla
f(\mathbf{x})\|_{\mathbf{y}}^{*}<1$, then,

$f(\mathbf{y})\leq f(\mathbf{x})+\langle\nabla
f(\mathbf{x}),\mathbf{y}-\mathbf{x}\rangle+\omega_{*}(\|\nabla
f(\mathbf{y})-\nabla f(\mathbf{x})\|_{\mathbf{y}}^{*}).$

Throughout this paper, we refer to notions such as super-linear and quadratic
convergence rates. Denote $\mathbf{x}_{k}$ the iterate generated by an
iterative process at the $k^{\text{th}}$ iteration. The sub-optimality gap of
the Newton method for self-concordant function satisfies the bound
$f(\mathbf{x}_{k})-f(\mathbf{x}^{*})\leq\lambda_{f}(\mathbf{x}_{k})^{2}$ which
holds for $\lambda_{f}(\mathbf{x}_{k})\leq 0.68$ (Boyd and Vandenberghe
[2004]) and thus one can estimate the convergence rate in terms of the local
norm of the gradient. It is known that the Newton method achieves a local
quadratic convergence rate. In this setting, we say that a process converges
quadratically if $\lambda_{f}(\mathbf{x}_{k+1})\leq
R\lambda_{f}(\mathbf{x}_{k})^{2}$, for $R>0$. In addition to quadratic
converge rate, we say that a process converges with super-linear rate if
$\lambda_{f}(\mathbf{x}_{k+1})/\lambda_{f}(\mathbf{x}_{k})\leq R(k)$, where
$R(k)\downarrow 0$.

###  A.2 Proof of Theorem 3.1

Recall the update rule of the Low-Rank Newton method at iteration $k$,

$\mathbf{x}_{k+1}=\mathbf{x}_{k}+t_{k}\hat{\mathbf{d}}_{h,k},$ (25)

where for the purposes of this section we use
$\hat{\mathbf{d}}_{h,k}:=-\mathbf{Q}_{h,k}^{-1}\nabla f(\mathbf{x}_{k})$,
where is defined in 10. Define also the corresponding approximate Newton
decrement $\bar{\lambda}(\mathbf{x}):=(\nabla
f(\mathbf{x}_{k})^{T}\mathbf{Q}_{h,k}^{-1}\nabla
f(\mathbf{x}_{k}))^{\frac{1}{2}}$.

######  Lemma A.3.

For any $k\in\mathbb{N}$ it holds that

$\displaystyle\bar{\lambda}(\mathbf{x}_{k})^{2}=-\nabla
f(\mathbf{x}_{k})^{T}\hat{\mathbf{d}}_{h,k}$ (26)
$\displaystyle\sqrt{\frac{\sigma_{k,n}}{\sigma_{k,N+1}}}\lambda(\mathbf{x}_{k})\leq\bar{\lambda}(\mathbf{x}_{k})\leq\lambda(\mathbf{x}_{k})$
(27)
$\displaystyle\|\hat{\mathbf{d}}_{h,k}\|_{\mathbf{x}}\leq\bar{\lambda}(\mathbf{x}_{k})$
(28)
$\displaystyle\left\|[\nabla^{2}f({\mathbf{x}_{k}})]^{1/2}\left(\mathbf{\hat{d}}_{h,k}-\mathbf{d}_{h,k}\right)\right\|_{2}\leq\left(1-\frac{\sigma_{k,n}}{\sigma_{k,N+1}}\right)\lambda(\mathbf{x}_{k}).$
(29)

where $\|\cdot\|_{\mathbf{x}}$ is defined in (20).

###### Proof.

Equality (26) follows directly by the definition of $\hat{\mathbf{d}}_{h,k}$.
For the right-hand side of (27) we get

$\displaystyle\bar{\lambda}(\mathbf{x}_{k})^{2}-\lambda(\mathbf{x}_{k})^{2}$
$\displaystyle=\nabla f(\mathbf{x}_{k})^{T}(\mathbf{Q}_{h,k}^{-1}-\nabla
f^{2}(\mathbf{x}_{k})^{-1})\nabla f(\mathbf{x}_{k})^{T}$ $\displaystyle=\nabla
f(\mathbf{x}_{k})^{T}[\nabla f^{2}(\mathbf{x}_{k})]^{-\frac{1}{2}}([\nabla
f^{2}(\mathbf{x}_{k})]^{\frac{1}{2}}\mathbf{Q}_{h,k}^{-1}[\nabla
f^{2}(\mathbf{x}_{k})]^{\frac{1}{2}}-\mathbf{I}_{n})[\nabla
f^{2}(\mathbf{x}_{k})]^{-\frac{1}{2}}\nabla f(\mathbf{x}_{k}).$

By construction of $\mathbf{Q}_{h,k}^{-1}$, it holds that

$\displaystyle[\nabla
f^{2}(\mathbf{x}_{k})]^{\frac{1}{2}}\mathbf{Q}_{h,k}^{-1}[\nabla
f^{2}(\mathbf{x}_{k})]^{\frac{1}{2}}-\mathbf{I}_{n}$
$\displaystyle\preceq\sigma_{\max}([\nabla
f^{2}(\mathbf{x}_{k})]^{\frac{1}{2}}\mathbf{Q}_{h,k}^{-1}[\nabla
f^{2}(\mathbf{x}_{k})]^{\frac{1}{2}}-\mathbf{I}_{n})\mathbf{I}_{n}$
$\displaystyle=\max_{i=1,\ldots.n}\left\\{\frac{\sigma_{i}(\nabla
f^{2}(\mathbf{x}_{k}))}{\sigma_{i}(\mathbf{Q}_{h,k})}-1\right\\}=0,$

and thus $\bar{\lambda}(\mathbf{x}_{k})-\lambda(\mathbf{x}_{k})\leq 0$. Using
similar arguments we take

$\displaystyle\lambda(\mathbf{x}_{k})^{2}-\bar{\lambda}(\mathbf{x}_{k})^{2}$
$\displaystyle\leq\max_{i=1,\ldots.n}\left\\{1-\frac{\sigma_{i}(\nabla
f^{2}(\mathbf{x}_{k}))}{\sigma_{i}(\mathbf{Q}_{h,k})}\right\\}\lambda(\mathbf{x})^{2}$
$\displaystyle=\left(1-\frac{\sigma_{k,n}}{\sigma_{k,N+1}}\right)\lambda(\mathbf{x})^{2}$

and thus the left-hand side inequality of (27) follows directly. As for the
result in (28), using the fact that $\mathbf{Q}_{h,k}$ is symmetric positive-
definite, we take

$\displaystyle\|\hat{\mathbf{d}}_{h,k}\|_{\mathbf{x}}^{2}$
$\displaystyle=\hat{\mathbf{d}}_{h,k}^{T}\nabla
f^{2}(\mathbf{x}_{k})\hat{\mathbf{d}}_{h,k}$ $\displaystyle=\nabla
f(\mathbf{x}_{k})^{T}[\mathbf{Q}_{h,k}]^{-\frac{1}{2}}\left([\mathbf{Q}_{h,k}]^{-\frac{1}{2}}\nabla
f^{2}(\mathbf{x}_{k})[\mathbf{Q}_{h,k}]^{-\frac{1}{2}}\right)[\mathbf{Q}_{h,k}]^{-\frac{1}{2}}f(\mathbf{x}_{k})$
$\displaystyle\leq\sigma_{\max}\left([\mathbf{Q}_{h,k}]^{-\frac{1}{2}}\nabla
f^{2}(\mathbf{x}_{k})[\mathbf{Q}_{h,k}]^{-\frac{1}{2}}\right)\bar{\lambda}(\mathbf{x}_{k})^{2},$

and since $\sigma_{\max}\left([\mathbf{Q}_{h,k}]^{-\frac{1}{2}}\nabla
f^{2}(\mathbf{x}_{k})[\mathbf{Q}_{h,k}]^{-\frac{1}{2}}\right)=1$ the result
follows. Last, we prove (29). Using the definitions of
$\mathbf{\hat{d}}_{h,k}$ and $\hat{d}_{h,k}$ we have

$\displaystyle\left\|[\nabla^{2}f({\mathbf{x}_{k}})]^{1/2}\left(\mathbf{\hat{d}}_{h,k}-\mathbf{d}_{h,k}\right)\right\|_{2}$
$\displaystyle=\left\|[\nabla^{2}f(\mathbf{x}_{k})]^{\frac{1}{2}}([\nabla^{2}f(\mathbf{x}_{k})]^{-1})-[\mathbf{Q}_{h,k}]^{-1})\nabla
f(\mathbf{x}_{k})\right\|_{2}$
$\displaystyle\leq\left\|\mathbf{I}_{n}-[\nabla^{2}f(\mathbf{x}_{k})]^{\frac{1}{2}}[\mathbf{Q}_{h,k}]^{-1}[\nabla^{2}f(\mathbf{x}_{k})]^{\frac{1}{2}}\right\|_{2}\bar{\lambda}(\mathbf{x})$
$\displaystyle=\max\left\\{1-\frac{\sigma_{1}}{\sigma_{1}},\ldots,1-\frac{\sigma_{N+1}}{\sigma_{N+1}},\ldots,1-\frac{\sigma_{n}}{\sigma_{N+1}}\right\\}\bar{\lambda}(\mathbf{x})$
$\displaystyle=\left(1-\frac{\sigma_{n}}{\sigma_{N+1}}\right)\bar{\lambda}(\mathbf{x})$

which concludes the proof of the lemma. ∎

######  Lemma A.4.

Let $\lambda(\mathbf{x}_{k})>\eta\ $ for some
$\eta\in(0,\frac{3-\sqrt{5}}{2})$. Then, there exists $\gamma>0$ such that the
coarse direction $\mathbf{\hat{d}}_{h,k}$ will yield reduction in value
function

$f(\mathbf{x}_{k}+t_{k}\mathbf{\hat{d}}_{h,k})-f(\mathbf{x}_{k})\leq-\gamma,$

for any $k\in\mathbb{N}$.

###### Proof.

By Lemma A.1 together with (25) and (26) we have that

$\displaystyle f(\mathbf{x}_{k}+t_{k}\mathbf{\hat{d}}_{h,k})$
$\displaystyle\leq f(\mathbf{x}_{k})+\langle\nabla
f(\mathbf{x}_{k}),\mathbf{x}_{k+1}-\mathbf{x}_{k}\rangle+\omega_{*}(\|\mathbf{x}_{k+1}-\mathbf{x}_{k}\|_{\mathbf{x}_{k}})$
$\displaystyle=f(\mathbf{x}_{k})+t_{k}\nabla
f(\mathbf{x}_{k})^{T}\mathbf{\hat{d}}_{h,k}+\omega_{*}(t_{k}\|\mathbf{\hat{d}}_{h,k}\|_{\mathbf{x}_{k}})$
$\displaystyle\leq
f(\mathbf{x}_{k})-t_{k}\bar{\lambda}(\mathbf{x})^{2}+\omega_{*}(t_{k}\bar{\lambda}(\mathbf{x}))$

where the last inequality follows from (28) and the fact that $\omega_{*}(x)$
is a monotone increasing function. Note that the above inequality is valid for
$t_{h,k}<\frac{1}{\hat{\lambda}(\mathbf{x}_{h,k})}$. Next, the right-hand side
is minimized at $t^{*}_{h}=\frac{1}{1+\hat{\lambda}(\mathbf{x}_{h,k})}$ and
thus

$\displaystyle f(\mathbf{x}_{k}+t^{*}\mathbf{\hat{d}}_{h,k})$
$\displaystyle\leq
f(\mathbf{x}_{k})-\frac{\bar{\lambda}(\mathbf{x})^{2}}{1+\bar{\lambda}(\mathbf{x}_{h,k})}+\omega_{*}\left(\frac{\bar{\lambda}(\mathbf{x})}{1+\bar{\lambda}(\mathbf{x}_{h,k})}\right)$
$\displaystyle=f(\mathbf{x}_{k})-\bar{\lambda}(\mathbf{x}_{k})+\log\left(1+\bar{\lambda}(\mathbf{x}_{k})\right).$

Using the inequality $-x+\log(1+x)\leq-\frac{x^{2}}{2(1+x)}$ for any $x>0$, we
obtain

$\displaystyle f(\mathbf{x}_{k}+t^{*}\mathbf{\hat{d}}_{h,k})$
$\displaystyle\leq
f(\mathbf{x}_{k})-\frac{\bar{\lambda}(\mathbf{x}_{k})^{2}}{2(1+\bar{\lambda}(\mathbf{x}_{k}))}$
$\displaystyle\leq f(\mathbf{x}_{k})-\alpha
t^{*}\bar{\lambda}(\mathbf{x}_{k})^{2}$
$\displaystyle=f(\mathbf{x}_{k})+\alpha t^{*}\nabla
f(\mathbf{x}_{k})^{T}\mathbf{\hat{d}}_{h,k},$

thus $t^{*}$ satisfies the back-tracking line search exit condition and that
it will always return a step size
$t_{k}>\beta/(1+\hat{\lambda}(\mathbf{x}_{k}))$. Therefore,

$f_{h}(\mathbf{x}_{k}+t_{k}\mathbf{\hat{d}}_{h,k})-f(\mathbf{x}_{k})\leq-\alpha\beta\frac{\bar{\lambda}(\mathbf{x}_{k})^{2}}{1+\bar{\lambda}(\mathbf{x}_{k})}.$

Additionally, combining the left-hand side in (27), the fact that the function
$x\rightarrow\frac{x^{2}}{1+x}$ is monotone increasing for any $x>0$ and since
by assumption $\lambda(\mathbf{x}_{k})>\eta$ we obtain

$f(\mathbf{x}_{k}+t_{k}\mathbf{\hat{d}}_{h,k})-f(\mathbf{x}_{k})\leq-\alpha\beta\frac{\eta^{2}}{1+\eta}.$

which concludes the proof by setting $\gamma:=\alpha\beta\eta^{2}/(1+\eta)$. ∎

Next, we estimate the sub-optimality gap of the process.

######  Lemma A.5.

Let $\lambda(\mathbf{x}_{k})<1$. Then,

$\omega(\lambda(\mathbf{x}_{k}))\leq
f(\mathbf{x}_{k})-f(\mathbf{x}^{*})\leq\omega_{*}(\lambda(\mathbf{x}_{k})).$

If in addition $\lambda(\mathbf{x}_{k})<0.68$ then,

$f(\mathbf{x}_{k})-f(\mathbf{x}^{*})\leq\lambda(\mathbf{x}_{k})^{2}.$

###### Proof.

The first result follows directly from Lemma A.2 Further, it holds that
$\omega_{*}(x)\leq x^{2},0\leq x\leq 0.68$ and thus we take the second result.
∎

Therefore, $\lambda(\mathbf{x}_{k})^{2}\leq\epsilon,\epsilon\in(0,0.68^{2})$,
can be used as an exit condition for the Low-Rank Newton method. Next, we
prove that the line search selects the unit step.

######  Lemma A.6.

If $\bar{\lambda}(\mathbf{x}_{k})\leq(1-2\alpha)/2$ then the Low-Rank Newton
method accepts the unit step, $t_{k}=1$.

###### Proof.

By Lemma A.1 we have

$\displaystyle f(\mathbf{x}_{k}+\mathbf{\hat{d}}_{h,k})$ $\displaystyle\leq
f(\mathbf{x}_{k})-\bar{\lambda}(\mathbf{x})^{2}+\omega_{*}(\bar{\lambda}(\mathbf{x}))$
$\displaystyle=f(\mathbf{x}_{k})-\bar{\lambda}(\mathbf{x})^{2}-\bar{\lambda}(\mathbf{x})-\log(1-\bar{\lambda}(\mathbf{x}))$

which holds since by assumption $\bar{\lambda}(\mathbf{x})<1$. Further,
$-x-\log(1-x)$ for all $x\in(0,0.81)$ and hence

$\displaystyle f(\mathbf{x}_{k}+\mathbf{\hat{d}}_{h,k})$ $\displaystyle\leq
f(\mathbf{x}_{k})-\bar{\lambda}(\mathbf{x})^{2}+\frac{1}{2}\bar{\lambda}(\mathbf{x})^{2}+\bar{\lambda}(\mathbf{x})^{3}$
$\displaystyle=f(\mathbf{x}_{k})-\frac{1}{2}(1-2\bar{\lambda}(\mathbf{x}))\bar{\lambda}(\mathbf{x})^{2}.$

Therefore, if $\bar{\lambda}(\mathbf{x}_{k})\leq(1-2\alpha)/2$ we take

$f(\mathbf{x}_{k}+\mathbf{\hat{d}}_{h,k})\leq
f(\mathbf{x}_{k})-\alpha\bar{\lambda}(\mathbf{x})^{2}$

which satisfies the line search condition for $t_{k}=1$. ∎

######  Lemma A.7.

Let $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ be a strictly convex self-
concordant function. If $\bar{\lambda}(\mathbf{x}_{k})<\frac{1}{t_{k}}$, we
have that

  1. (i)

$\nabla^{2}f(\mathbf{x}_{k+1})\preceq\frac{1}{(1-t_{k}\bar{\lambda}(\mathbf{x}_{k}))^{2}}\nabla^{2}f(\mathbf{x}_{k}),$

  2. (ii)

$[\nabla^{2}f(\mathbf{x}_{k+1})]^{-1}\preceq\frac{1}{(1-t_{k}\bar{\lambda}(\mathbf{x}_{k}))^{2}}[\nabla^{2}f(\mathbf{x}_{k})]^{-1}$.

###### Proof.

The proof is analogue to Lemma A.11 below but with
$\bar{\lambda}(\mathbf{x}_{k})$ instead of $\hat{\lambda}(\mathbf{x}_{k})$. ∎

The next lemma shows the local super-linear rate of the Low-Rank Newton
method.

######  Lemma A.8.

Suppose that the sequence is generated by (25) with $t_{k}=1$. Then,

$\lambda(\mathbf{x}_{k+1})\leq\frac{1-\frac{\sigma_{k,n}}{\sigma_{k,N+1}}+\lambda(\mathbf{x}_{k})}{\left(1-\lambda(\mathbf{x}_{k})\right)^{2}}\lambda(\mathbf{x}_{k}).$

###### Proof.

By Lemma A.7, the upper bound in (27) and $t_{k}=1$ we have

$\nabla^{2}f(\mathbf{x}_{k+1})\preceq\frac{1}{(1-\lambda(\mathbf{x}_{k}))^{2}}\nabla^{2}f(\mathbf{x}_{k}).$
(30)

Now, by the definition of the Newton decrement we have that

$\lambda(\mathbf{x}_{k+1})=\left[\nabla
f(\mathbf{x}_{k+1})^{T}[\nabla^{2}f(\mathbf{x}_{k+1})]^{-1}\nabla
f(\mathbf{x}_{k+1})\right]^{1/2},$

and combining this with inequality (30) we take

$\lambda(\mathbf{x}_{k+1})\leq\frac{1}{1-\lambda(\mathbf{x}_{k})}\left\|[\nabla^{2}f({\mathbf{x}_{k}})]^{-1/2}\nabla
f(\mathbf{x}_{k+1})\right\|_{2}.$ (31)

Denote $\mathbf{Z}:=[\nabla^{2}f({\mathbf{x}_{k}})]^{-1/2}\nabla
f(\mathbf{x}_{k+1})$. Using the fact that

$\nabla
f(\mathbf{x}_{k+1})=\int_{0}^{1}\nabla^{2}f({\mathbf{x}_{k}}+y\mathbf{\hat{d}}_{h,k})\mathbf{\hat{d}}_{h,k}\
dy+\nabla f(\mathbf{x}_{k})$

we see that

$\displaystyle\mathbf{Z}$
$\displaystyle=[\nabla^{2}f({\mathbf{x}_{k}})]^{-1/2}\left(\int_{0}^{1}\nabla^{2}f({\mathbf{x}_{k}}+y\mathbf{\hat{d}}_{h,k})\mathbf{\hat{d}}_{h,k}\
dy+\nabla f(\mathbf{x}_{k})\right)$
$\displaystyle=\underbrace{\int_{0}^{1}[\nabla^{2}f({\mathbf{x}_{k}})]^{-1/2}\nabla^{2}f({\mathbf{x}_{k}}+y\mathbf{\hat{d}}_{h,k})[\nabla^{2}f({\mathbf{x}_{k}})]^{-1/2}\
dy}_{\mathbf{T}}\
[\nabla^{2}f({\mathbf{x}_{k}})]^{1/2}\mathbf{\hat{d}}_{h,k}-[\nabla^{2}f({\mathbf{x}_{k}})]^{1/2}\mathbf{d}_{h,k},$

where $\mathbf{d}_{h,k}$ is the Newton direction. Next adding and subtracting
the quantity $\mathbf{T}[\nabla^{2}f({\mathbf{x}_{k}})]^{1/2}\mathbf{d}_{h,k}$
we have that

$\mathbf{Z}=\mathbf{T}[\nabla^{2}f({\mathbf{x}_{k}})]^{1/2}(\mathbf{\hat{d}}_{h,k}-\mathbf{d}_{h,k})+(\mathbf{T}-\mathbf{I}_{N\times
N})[\nabla^{2}f({\mathbf{x}_{k}})]^{1/2}\mathbf{d}_{h,k},$

and thus

$\left\|\mathbf{Z}\right\|\leq\left\|\underbrace{\mathbf{T}[\nabla^{2}f({\mathbf{x}_{k}})]^{1/2}(\mathbf{\hat{d}}_{h,k}-\mathbf{d}_{h,k})}_{\mathbf{Z}_{1}}\right\|_{2}+\left\|\underbrace{(\mathbf{T}-\mathbf{I}_{N\times
N})[\nabla^{2}f({\mathbf{x}_{k}})]^{1/2}\mathbf{d}_{h,k}}_{\mathbf{Z}_{2}}\right\|_{2}\\\
$ (32)

Using again (30) and since, by assumption, $y\hat{\lambda}(\mathbf{x}_{k})<1$
we take

$[\nabla^{2}f({\mathbf{x}_{k}})]^{-1/2}\nabla^{2}f({\mathbf{x}_{k}}+y\mathbf{\hat{d}}_{h,k})[\nabla^{2}f({\mathbf{x}_{k}})]^{-1/2}\preceq\frac{1}{\left(1-y\lambda(\mathbf{x}_{k})\right)^{2}}\mathbf{I}_{N\times
N}.$

We are now in position to estimate both norms in (32). For the first one we
have that

$\displaystyle\left\|\mathbf{Z}_{1}\right\|$
$\displaystyle\leq\left\|\int_{0}^{1}\frac{1}{\left(1-y\lambda(\mathbf{x}_{k})\right)^{2}}dy\right\|_{2}\left\|\
[\nabla^{2}f({\mathbf{x}_{k}})]^{1/2}\left(\mathbf{\hat{d}}_{h,k}-\mathbf{d}_{h,k}\right)\right\|_{2}$
$\displaystyle=\frac{1}{1-\lambda(\mathbf{x}_{k})}\left\|[\nabla^{2}f({\mathbf{x}_{k}})]^{1/2}\left(\mathbf{\hat{d}}_{h,k}-\mathbf{d}_{h,k}\right)\right\|_{2}$
$\displaystyle\leq\left(1-\frac{\sigma_{k,n}}{\sigma_{k,N+1}}\right)\frac{\lambda(\mathbf{x}_{k})}{1-\lambda(\mathbf{x}_{h,k})}.$

where the last inequality follows from (29). Next, the second norm implies

$\displaystyle\left\|\mathbf{Z}_{2}\right\|$
$\displaystyle\leq\left\|\int_{0}^{1}\left(\frac{1}{\left(1-y\lambda(\mathbf{x}_{k})\right)^{2}}-1\right)dy\right\|_{2}\left\|[\nabla^{2}f({\mathbf{x}_{k}})]^{1/2}\mathbf{d}_{h,k}\right\|_{2}$
$\displaystyle=\frac{\lambda(\mathbf{x}_{k})}{1-\lambda(\mathbf{x}_{k})}\lambda(\mathbf{x}_{k}).$

Putting this all together, inequality (31) becomes

$\lambda(\mathbf{x}_{k+1})\leq\frac{1-\frac{\sigma_{k,n}}{\sigma_{k,N+1}}+\lambda(\mathbf{x}_{k})}{\left(1-\lambda(\mathbf{x}_{k})\right)^{2}}\lambda(\mathbf{x}_{k}).$

as required. ∎

Recall that
$\varepsilon:=\min_{k\in\mathbb{N}}\\{\frac{\sigma_{k,n}}{\sigma_{k,N+1}}\\}$
and $\eta:=(3-\sqrt{9-4\varepsilon})/2$. Below we prove Theorem 3.1.

###### proof of Theorem 3.1.

By Lemma A.4, we see that one step of the first phase of the Low-Rank Newton
method decreases the objective by $\gamma:=\alpha\beta\eta^{2}/(1+\eta)>0$. By
Lemma A.4 we also note that the line-search can be applied to determine the
step size parameter. This proves the result of the first phase of the Low-Rank
Newton method. Next, by Lemma A.8, since
$\frac{\sigma_{k,n}}{\sigma_{k,N+1}}\geq\varepsilon$ and
$\lambda(\mathbf{x}_{k})\leq\eta$ by assumption, we take,

$\lambda(\mathbf{x}_{k+1})\leq\frac{1-\varepsilon+\eta}{\left(1-\eta\right)^{2}}\lambda(\mathbf{x}_{k}).$

Moreover, by definition $\eta:=(3-\sqrt{9-4\varepsilon})$ and since
$0<\varepsilon\leq 1$ we obtain
$\frac{\hat{e}+\hat{\eta}}{\left(1-\hat{\eta}\right)^{2}}<1$. This shows that
$\lambda(\mathbf{x}_{k+1})<\lambda(\mathbf{x}_{k})$ and, in particular, this
process converges with super-linear rate for $\lambda(\mathbf{x}_{k})\leq\eta$
and some $\eta\in(0,\frac{3-\sqrt{5}}{2})$. Finally, by Lemma A.6, for
$\lambda(\mathbf{x}_{k})\leq\eta$, the line search accepts the unit step which
concludes the proof of the theorem. ∎

###  A.3 Proof of Theorem 3.4

Theorem 3.4 considers a sequence generated by
$\mathbf{x}_{k+1}=\mathbf{x}_{k}-t_{k}\mathbf{P}\mathbf{Q}_{H,k}^{-1}\mathbf{R}\nabla
f(\mathbf{x}_{k})$, thus for the purposes of this section we denote
$\hat{\mathbf{d}}_{h,k}:=-\mathbf{P}\mathbf{Q}_{H,k}^{-1}\mathbf{R}\nabla
f(\mathbf{x}_{k})$, where $\mathbf{Q}_{H,k}^{-1}$ is defined in (14). Below we
will make use of the approximate decrements $\hat{\lambda}(\mathbf{x})$ and
$\tilde{\lambda}(\mathbf{x})$ defined in (15), respectively.

Before proving Lemma 3.2 we state an upper bound for
$\hat{\lambda}(\mathbf{x})$ which was proved in Tsipinakis and Parpas [2021].

######  Lemma A.9 (Tsipinakis and Parpas [2021]).

Let $\lambda(\mathbf{x}_{h,k})$ be the newton decrement in (21). For any
$k\in\mathbb{N}$ it holds that

$\tilde{\lambda}(\mathbf{x}_{h,k})\leq\lambda(\mathbf{x}_{h,k}).$ (33)

###### Proof.

By the definition of $\tilde{\lambda}(\mathbf{x}_{k})$ in (15) we have that

$\displaystyle\tilde{\lambda}(\mathbf{x}_{k})$ $\displaystyle=\left[\nabla
f(\mathbf{x}_{k})^{T}\mathbf{P}[\mathbf{R}\nabla^{2}f(\mathbf{x}_{k})\mathbf{P}]^{-1}\mathbf{R}\nabla
f(\mathbf{x}_{k})\right]^{1/2}.$ $\displaystyle=\left[\nabla
f(\mathbf{x}_{k})^{T}[\nabla^{2}f(\mathbf{x}_{k})]^{-\frac{1}{2}}\Pi_{[\nabla^{2}f(\mathbf{x}_{k})]^{\frac{1}{2}}\mathbf{P}}[\nabla^{2}f(\mathbf{x}_{k})]^{-\frac{1}{2}}\nabla
f(\mathbf{x}_{k})\right]^{1/2},$

where
$\Pi_{[\nabla^{2}f(\mathbf{x}_{k})]^{\frac{1}{2}}\mathbf{P}}:=([\nabla^{2}f(\mathbf{x}_{k})]^{\frac{1}{2}}\mathbf{P})[([\nabla^{2}f(\mathbf{x}_{k})]^{\frac{1}{2}}\mathbf{P})^{T}([\nabla^{2}f(\mathbf{x}_{k})]^{\frac{1}{2}}\mathbf{P})]^{-1}([\nabla^{2}f(\mathbf{x}_{k})]^{\frac{1}{2}}\mathbf{P})^{T}$
is the orthogonal projection onto the
$\operatorname{span}\\{[\nabla^{2}f(\mathbf{x}_{k})]^{\frac{1}{2}}\mathbf{P}\\}$.
By the idepotency of the orthogonal projection we have that

$\displaystyle\tilde{\lambda}(\mathbf{x}_{k})$
$\displaystyle=\left\langle\Pi_{[\nabla^{2}f(\mathbf{x}_{k})]^{\frac{1}{2}}\mathbf{P}}[\nabla^{2}f(\mathbf{x}_{h,k})]^{-\frac{1}{2}}\nabla
f(\mathbf{x}_{k}),\Pi_{[\nabla^{2}f(\mathbf{x}_{k})]^{\frac{1}{2}}\mathbf{P}}[\nabla^{2}f(\mathbf{x}_{k})]^{-\frac{1}{2}}\nabla
f(\mathbf{x}_{k})\right\rangle^{1/2}$
$\displaystyle\leq\left\|\Pi_{[\nabla^{2}f(\mathbf{x}_{k})]^{\frac{1}{2}}\mathbf{P}}\right\|\left\|[\nabla^{2}f(\mathbf{x}_{k})]^{-\frac{1}{2}}\nabla
f(\mathbf{x}_{k})\right\|.$

Since
$\sigma_{1}(\Pi_{[\nabla^{2}f(\mathbf{x}_{k})]^{\frac{1}{2}}\mathbf{P}})=1$ we
take $\tilde{\lambda}(\mathbf{x}_{k})\leq\lambda(\mathbf{x}_{k})$ as required.
∎

Note that since
$(\nabla^{2}f(\mathbf{x}_{k})^{\frac{1}{2}}\mathbf{P})^{T}(\nabla^{2}f(\mathbf{x}_{k})^{\frac{1}{2}}\mathbf{P}$
is always invertible the above result holds with probability one. Below we
prove Lemma 3.2. We state again Lemma 3.2 for convenience.

######  Lemma A.10.

For all $k\in\mathbb{N}$ we have

$\displaystyle\mathbf{d}_{h,k}^{T}\nabla^{2}f(\mathbf{x}_{k})\hat{\mathbf{d}}_{h,k}=\hat{\lambda}^{2}(\mathbf{x}_{k})$
(34)
$\displaystyle\mathbf{\hat{d}}_{h,k}^{T}\nabla^{2}f(\mathbf{x}_{k})\hat{\mathbf{d}}_{h,k}\leq\hat{\lambda}^{2}(\mathbf{x}_{k})$
(35)
$\displaystyle\|[\nabla^{2}f(\mathbf{x}_{k})]^{\frac{1}{2}}(\mathbf{d}_{h,k}-\hat{\mathbf{d}}_{h,k})\|\leq\hat{e}_{k}$
(36)
$\displaystyle\sqrt{\frac{\sigma_{k,N}}{\sigma_{k,p+1}}}\tilde{\lambda}(\mathbf{x}_{k})\leq\hat{\lambda}(\mathbf{x}_{k})\leq\tilde{\lambda}(\mathbf{x}_{k})\leq\lambda(\mathbf{x}_{k}),$
(37)

where
$\hat{e}_{k}:=\sqrt{\lambda(\mathbf{x}_{k})^{2}-\hat{\lambda}(\mathbf{x}_{k})^{2}}$.

Further, suppose that for all $k\in\mathbb{N}$ such that
$\mathbf{x}_{k}\neq\mathbf{x}^{*}$ it holds $\hat{\lambda}(\mathbf{x}_{k})>0$.
Then, there exists $\mu_{k}\in(0,1]$ such that
$\hat{e}_{k}\leq(1-\frac{\sigma_{k,N}}{\sigma_{k,p+1}}\mu_{k}^{2})^{1/2}\lambda(\mathbf{x}_{k})$.

###### Proof of Lemma 3.2.

Equality (34) follows directly from the definitions of
$\hat{\mathbf{d}}_{h,k}$ and $\mathbf{d}_{h,k}$. By (34) we also get
$\hat{\lambda}^{2}(\mathbf{x}_{k})=-\nabla
f(\mathbf{x}_{k})^{T}\hat{\mathbf{d}}_{h,k}$. For inequality (35), since
$\mathbf{Q}_{H,k}$ is symmetric positive definite, we have

$\displaystyle\mathbf{\hat{d}}_{h,k}^{T}\nabla^{2}f(\mathbf{x}_{k})\hat{\mathbf{d}}_{h,k}$
$\displaystyle=\nabla
f(\mathbf{x}_{k})^{T}\mathbf{P}\mathbf{Q}_{H,k}^{-1}\mathbf{R}\nabla^{2}f(\mathbf{x}_{k})\mathbf{P}\mathbf{Q}_{H,k}^{-1}\mathbf{R}\nabla
f(\mathbf{x}_{k})$
$\displaystyle=\mathbf{z}^{T}\mathbf{Q}_{H,k}^{-\frac{1}{2}}(\mathbf{R}\nabla^{2}f(\mathbf{x}_{k})\mathbf{P})\mathbf{Q}_{H,k}^{-\frac{1}{2}}\mathbf{z}$

where $\mathbf{z}:=\mathbf{Q}_{H,k}^{-\frac{1}{2}}\mathbf{R}\nabla
f(\mathbf{x}_{k})$ and note that
$\mathbf{z}^{T}\mathbf{z}=\hat{\lambda}^{2}(\mathbf{x}_{k})$. By construction
of $\mathbf{Q}_{H,k}^{-1}$, it holds that

$\displaystyle\mathbf{Q}_{H,k}^{-\frac{1}{2}}(\mathbf{R}\nabla^{2}f(\mathbf{x}_{k})\mathbf{P})\mathbf{Q}_{H,k}^{-\frac{1}{2}}$
$\displaystyle\preceq\sigma_{\max}(\mathbf{Q}_{H,k}^{-\frac{1}{2}}(\mathbf{R}\nabla^{2}f(\mathbf{x}_{k})\mathbf{P})\mathbf{Q}_{H,k}^{-\frac{1}{2}})\mathbf{I}_{N}$
$\displaystyle=\max\left\\{\frac{\sigma_{k,1}}{\sigma_{k,1}},\ldots,\frac{\sigma_{k,p+1}}{\sigma_{k,p+1}},\frac{\sigma_{k,p+2}}{\sigma_{k,p+1}},\ldots,\frac{\sigma_{k,N}}{\sigma_{k,p+1}}\right\\}\mathbf{I}_{N}$
$\displaystyle=\mathbf{I}_{N}$

Then,
$\mathbf{z}^{T}\mathbf{Q}_{H,k}^{-\frac{1}{2}}(\mathbf{R}\nabla^{2}f(\mathbf{x}_{k})\mathbf{P})\mathbf{Q}_{H,k}^{-\frac{1}{2}}\mathbf{z}\leq\hat{\lambda}^{2}(\mathbf{x}_{k})$
and (35) is proved. Next,

$\displaystyle\|[\nabla^{2}f(\mathbf{x}_{k})]^{\frac{1}{2}}(\mathbf{d}_{h,k}-\hat{\mathbf{d}}_{h,k})\|^{2}$
$\displaystyle=\mathbf{d}_{h,k}^{T}\nabla^{2}f(\mathbf{x}_{k})\mathbf{d}_{h,k}-2\mathbf{d}_{h,k}^{T}\nabla^{2}f(\mathbf{x}_{k})\hat{\mathbf{d}}_{h,k}+\mathbf{\hat{d}}_{h,k}^{T}\nabla^{2}f(\mathbf{x}_{k})\hat{\mathbf{d}}_{h,k}$
$\displaystyle\leq\lambda(\mathbf{x}_{k})^{2}-\hat{\lambda}(\mathbf{x}_{k})^{2},$

where the last inequality follows from (21), (34) and (35). This proves
inequality (36). For the bounds in (37) we start by proving
$\hat{\lambda}(\mathbf{x}_{k})\leq\tilde{\lambda}(\mathbf{x}_{k})$. We have

$\displaystyle\hat{\lambda}(\mathbf{x}_{k})^{2}-\tilde{\lambda}(\mathbf{x}_{k})^{2}$
$\displaystyle=(\mathbf{R}\nabla
f(\mathbf{x}_{k}))^{T}(\mathbf{Q}_{H,k}^{-1}-[\mathbf{R}\nabla^{2}f(\mathbf{x}_{k})\mathbf{P}]^{-1})(\mathbf{R}\nabla
f(\mathbf{x}_{k}))$
$\displaystyle=(\mathbf{Q}_{H,k}^{-\frac{1}{2}}\mathbf{R}\nabla
f(\mathbf{x}_{k}))^{T}(\mathbf{I}_{N}-\mathbf{Q}_{H,k}^{\frac{1}{2}}[\mathbf{R}\nabla^{2}f(\mathbf{x}_{k})\mathbf{P}]^{-1}\mathbf{Q}_{H,k}^{\frac{1}{2}})(\mathbf{Q}_{H,k}^{-\frac{1}{2}}\mathbf{R}\nabla
f(\mathbf{x}_{k}))$

We also have

$\displaystyle\mathbf{I}_{N}-\mathbf{Q}_{H,k}^{\frac{1}{2}}[\mathbf{R}\nabla^{2}f(\mathbf{x}_{k})\mathbf{P}]^{-1}\mathbf{Q}_{H,k}^{\frac{1}{2}}$
$\displaystyle\preceq\sigma_{\max}(\mathbf{I}_{N}-\mathbf{Q}_{H,k}^{\frac{1}{2}}[\mathbf{R}\nabla^{2}f(\mathbf{x}_{k})\mathbf{P}]^{-1}\mathbf{Q}_{H,k}^{\frac{1}{2}})\mathbf{I}_{N}$
$\displaystyle=\max\left\\{1-\frac{\sigma_{k,1}}{\sigma_{k,1}},\ldots,1-\frac{\sigma_{k,p+1}}{\sigma_{k,p+1}},\ldots,1-\frac{\sigma_{k,p+1}}{\sigma_{k,N}}\ldots\right\\}\mathbf{I}_{N}$
$\displaystyle=0.$

Putting this all together, we take
$\hat{\lambda}(\mathbf{x}_{k})\leq\tilde{\lambda}(\mathbf{x}_{k})$. On the
other hand, with similar arguments we have that

$\displaystyle\tilde{\lambda}(\mathbf{x}_{k})^{2}-\hat{\lambda}(\mathbf{x}_{k})^{2}$
$\displaystyle=\mathbf{z}^{T}(\mathbf{Q}_{H,k}^{\frac{1}{2}}[\mathbf{R}\nabla^{2}f(\mathbf{x}_{k})\mathbf{P}]^{-1}\mathbf{Q}_{H,k}^{\frac{1}{2}}-\mathbf{I}_{N})\mathbf{z}$
$\displaystyle\leq\left(\frac{\sigma_{k,p+1}}{\sigma_{k,N}}-1\right)\mathbf{z}^{T}\mathbf{z}=\left(\frac{\sigma_{k,p+1}}{\sigma_{k,N}}-1\right)\hat{\lambda}(\mathbf{x}_{k})^{2}$

and thus the lower bound of (37) follows. Putting this all together and
combining it with Lemma A.9, (37) has been proved. Moreover, assumption
$\hat{\lambda}(\mathbf{x}_{k})>0$ and (37) imply
$\tilde{\lambda}(\mathbf{x}_{k})>0$. Then there exists
$\mu_{k}\in(0,\frac{\tilde{\lambda}(\mathbf{x}_{k})}{\lambda(\mathbf{x}_{k})}]$.
Thus, we obtain

$\mu_{k}\sqrt{\frac{\sigma_{k,N}}{\sigma_{k,p+1}}}\lambda(\mathbf{x}_{k})\leq\hat{\lambda}(\mathbf{x}_{k})\leq\lambda(\mathbf{x}_{k}).$
(38)

Therefore, by the lower bound of the last result we directly get
$\hat{e}_{k}\leq(1-\frac{\sigma_{k,N}}{\sigma_{k,p+1}}\mu_{k}^{2})^{1/2}\lambda(\mathbf{x}_{k})$,
which concludes the proof. ∎

######  Lemma A.11.

Let $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ be a strictly convex self-
concordant function and
$\mathbf{x}_{k+1}=\mathbf{x}_{k}-t_{k}\mathbf{P}\mathbf{Q}_{H,k}^{-1}\mathbf{R}\nabla
f(\mathbf{x}_{k})$. If $\hat{\lambda}(\mathbf{x}_{h,k})<\frac{1}{t_{h,k}}$, we
have that

  1. (i)

$\nabla^{2}f(\mathbf{x}_{k+1})\preceq\frac{1}{(1-t_{k}\hat{\lambda}(\mathbf{x}_{k}))^{2}}\nabla^{2}f(\mathbf{x}_{k}),$

  2. (ii)

$[\nabla^{2}f(\mathbf{x}_{k+1})]^{-1}\preceq\frac{1}{(1-t_{k}\hat{\lambda}(\mathbf{x}_{k}))^{2}}[\nabla^{2}f(\mathbf{x}_{k})]^{-1}$.

###### Proof.

Consider the case (i). From the upper bound in (24) that arise for self-
concordant functions we have that

$\displaystyle\nabla^{2}f(\mathbf{x}_{k+1})$
$\displaystyle\preceq\frac{1}{(1-t_{k}\|\hat{\mathbf{d}}_{h,k}\|_{\mathbf{x}_{k}})^{2}}\nabla^{2}f(\mathbf{x}_{k})$
$\displaystyle\preceq\frac{1}{(1-t_{k}\hat{\lambda}(\mathbf{x}_{k}))^{2}}\nabla^{2}f(\mathbf{x}_{k}).$

where the last inequality holds from (35). As for the case (ii), we make use
of the lower bound in (24), and thus, by (35), we have

$\nabla^{2}f(\mathbf{x}_{k+1})\succeq(1-t_{k}\hat{\lambda}(\mathbf{x}_{k}))^{2}\nabla^{2}f(\mathbf{x}_{k}).$

Since, further, $f$ is strictly convex we take

$[\nabla^{2}f(\mathbf{x}_{k+1})]^{-1}\preceq\frac{1}{(1-t_{k}\hat{\lambda}(\mathbf{x}_{k}))^{2}}[\nabla^{2}f(\mathbf{x}_{k})]^{-1},$

which concludes the proof. ∎

The next lemma estimates the sub-optimality gad of the process.

######  Lemma A.12.

Let $\lambda(\mathbf{x}_{k})<1$. Then,

$\omega(\lambda(\mathbf{x}_{k}))\leq
f(\mathbf{x}_{k})-f(\mathbf{x}^{*})\leq\omega_{*}(\lambda(\mathbf{x}_{k})).$

If in addition $\lambda(\mathbf{x}_{k})<0.68$ then,

$f(\mathbf{x}_{k})-f(\mathbf{x}^{*})\leq\lambda(\mathbf{x}_{k})^{2}.$

###### Proof.

The first result follows directly from Lemma A.2 Further, it holds that
$\omega_{*}(x)\leq x^{2},0\leq x\leq 0.68$ and thus we take the second result.
∎

Thus $\lambda(\mathbf{x}_{k})$ can be used as an exit condition for Algorithm
1. Combining the result above with (38) one may take

$f(\mathbf{x}_{k})-f(\mathbf{x}^{*})\leq\omega_{*}\left(\frac{\sigma_{k,p+1}^{1/2}}{\mu_{k}\sigma_{k,N}^{1/2}}\hat{\lambda}(\mathbf{x}_{k})\right)$

and thus, for
$\hat{\lambda}(\mathbf{x}_{k})<\mu_{k}\sqrt{\frac{\sigma_{k,N}}{\sigma_{k,p+1}}}$,
$\hat{\lambda}(\mathbf{x}_{k})$ can be used as an exit condition for Algorithm
1.

######  Lemma A.13.

Let $t_{k}:=\frac{1}{1+\hat{\lambda}(\mathbf{x}_{k})}$. Then, for any
$k\in\mathbb{N}$ we have

$f(\mathbf{x}_{k+1})-f(\mathbf{x}_{k})\leq-\omega(\hat{\lambda}(\mathbf{x}))$

where $\omega(x)$ is defined in (22).

###### Proof.

By Lemma A.1 we have that

$\displaystyle f(\mathbf{x}_{k}+t_{k}\mathbf{\hat{d}}_{h,k})$
$\displaystyle\leq f(\mathbf{x}_{k})+\langle\nabla
f(\mathbf{x}_{k}),\mathbf{x}_{k+1}-\mathbf{x}_{k}\rangle+\omega_{*}(\|\mathbf{x}_{k+1}-\mathbf{x}_{k}\|_{\mathbf{x}_{k}})$
$\displaystyle=f(\mathbf{x}_{k})+t_{k}\nabla
f(\mathbf{x}_{k})^{T}\mathbf{\hat{d}}_{h,k}+\omega_{*}(t_{k}\|\mathbf{\hat{d}}_{h,k}\|_{\mathbf{x}_{k}})$
$\displaystyle\leq
f(\mathbf{x}_{k})-t_{k}\hat{\lambda}(\mathbf{x})^{2}+\omega_{*}(t_{k}\hat{\lambda}(\mathbf{x}))$

where for the last inequality we used (35), the monotonicity of
$\omega_{*}(x)$ and $\hat{\lambda}(\mathbf{x})^{2}=-\nabla
f(\mathbf{x}_{k})^{T}\mathbf{\hat{d}}_{h,k}$. Then, using the definition of
$t_{k}$ we take

$\displaystyle f(\mathbf{x}_{k}+t_{k}\mathbf{\hat{d}}_{h,k})$
$\displaystyle\leq
f(\mathbf{x}_{k})-\frac{\hat{\lambda}(\mathbf{x})^{2}}{1+\hat{\lambda}(\mathbf{x}_{h,k})}+\omega_{*}\left(\frac{\hat{\lambda}(\mathbf{x})}{1+\hat{\lambda}(\mathbf{x}_{h,k})}\right)$
$\displaystyle=f(\mathbf{x}_{k})-\hat{\lambda}(\mathbf{x}_{k})+\log\left(1+\hat{\lambda}(\mathbf{x}_{k})\right).$

which concludes the proof. ∎

For the next result we will need the following definitions:
$\hat{M}_{k}:=\mathbb{E}[\mu_{k}\sqrt{\frac{\sigma_{k,N}}{\sigma_{k,p+1}}}|\mathcal{F}_{k}]$,
$\hat{M}:=\min_{k\in\mathbb{N}}\\{\hat{M}_{k}\\}$ and
$\hat{\varepsilon}:=\sqrt{1-\hat{M}^{2}}$.

######  Lemma A.14.

Suppose that for all $k\in\mathbb{N}$ such that
$\mathbf{x}_{k}\neq\mathbf{x}^{*}$ it holds $\hat{\lambda}(\mathbf{x}_{k})>0$
and suppose also that $t_{k}=1$. Then,

$\mathbb{E}\left[\lambda(\mathbf{x}_{k+1})|\mathbf{x}_{k}\right]\leq\frac{\hat{\varepsilon}+\lambda(\mathbf{x}_{k})}{\left(1-\lambda(\mathbf{x}_{k})\right)^{2}}\lambda(\mathbf{x}_{k}),$

which holds almost surely.

###### Proof.

By Lemma A.11, inequality (37) and $t_{k}=1$ we have

$\nabla^{2}f(\mathbf{x}_{k+1})\preceq\frac{1}{(1-\lambda(\mathbf{x}_{k}))^{2}}\nabla^{2}f(\mathbf{x}_{k}),$
(39)

with probability one. Now, by the definition of the Newton decrement we have
that

$\lambda(\mathbf{x}_{k+1})=\left[\nabla
f(\mathbf{x}_{k+1})^{T}[\nabla^{2}f(\mathbf{x}_{k+1})]^{-1}\nabla
f(\mathbf{x}_{k+1})\right]^{1/2},$

and combining this with inequality (39) we take

$\lambda(\mathbf{x}_{k+1})\leq\frac{1}{1-\lambda(\mathbf{x}_{k})}\left\|[\nabla^{2}f({\mathbf{x}_{k}})]^{-1/2}\nabla
f(\mathbf{x}_{k+1})\right\|_{2}.$ (40)

Denote $\mathbf{Z}:=[\nabla^{2}f({\mathbf{x}_{k}})]^{-1/2}\nabla
f(\mathbf{x}_{k+1})$. Using the fact that

$\nabla
f(\mathbf{x}_{k+1})=\int_{0}^{1}\nabla^{2}f({\mathbf{x}_{k}}+y\mathbf{\hat{d}}_{h,k})\mathbf{\hat{d}}_{h,k}\
dy+\nabla f(\mathbf{x}_{k})$

we see that

$\displaystyle\mathbf{Z}$
$\displaystyle=[\nabla^{2}f({\mathbf{x}_{k}})]^{-1/2}\left(\int_{0}^{1}\nabla^{2}f({\mathbf{x}_{k}}+y\mathbf{\hat{d}}_{h,k})\mathbf{\hat{d}}_{h,k}\
dy+\nabla f(\mathbf{x}_{k})\right)$
$\displaystyle=\underbrace{\int_{0}^{1}[\nabla^{2}f({\mathbf{x}_{k}})]^{-1/2}\nabla^{2}f({\mathbf{x}_{k}}+y\mathbf{\hat{d}}_{h,k})[\nabla^{2}f({\mathbf{x}_{k}})]^{-1/2}\
dy}_{\mathbf{T}}\
[\nabla^{2}f({\mathbf{x}_{k}})]^{1/2}\mathbf{\hat{d}}_{h,k}-[\nabla^{2}f({\mathbf{x}_{k}})]^{1/2}\mathbf{d}_{h,k},$

where $\mathbf{d}_{h,k}$ is the Newton direction. Next adding and subtracting
the quantity $\mathbf{T}[\nabla^{2}f({\mathbf{x}_{k}})]^{1/2}\mathbf{d}_{h,k}$
we have that

$\mathbf{Z}=\mathbf{T}[\nabla^{2}f({\mathbf{x}_{k}})]^{1/2}(\mathbf{\hat{d}}_{h,k}-\mathbf{d}_{h,k})+(\mathbf{T}-\mathbf{I}_{N\times
N})[\nabla^{2}f({\mathbf{x}_{k}})]^{1/2}\mathbf{d}_{h,k},$

and thus

$\left\|\mathbf{Z}\right\|\leq\left\|\underbrace{\mathbf{T}[\nabla^{2}f({\mathbf{x}_{k}})]^{1/2}(\mathbf{\hat{d}}_{h,k}-\mathbf{d}_{h,k})}_{\mathbf{Z}_{1}}\right\|_{2}+\left\|\underbrace{(\mathbf{T}-\mathbf{I}_{N\times
N})[\nabla^{2}f({\mathbf{x}_{k}})]^{1/2}\mathbf{d}_{h,k}}_{\mathbf{Z}_{2}}\right\|_{2}\\\
$ (41)

Using again (39) and since, by assumption, $y\hat{\lambda}(\mathbf{x}_{k})<1$
we take

$[\nabla^{2}f({\mathbf{x}_{k}})]^{-1/2}\nabla^{2}f({\mathbf{x}_{k}}+y\mathbf{\hat{d}}_{h,k})[\nabla^{2}f({\mathbf{x}_{k}})]^{-1/2}\preceq\frac{1}{\left(1-y\lambda(\mathbf{x}_{k})\right)^{2}}\mathbf{I}_{N\times
N}.$

We are now in position to estimate both norms in (41). For the first one we
have that

$\displaystyle\left\|\mathbf{Z}_{1}\right\|$
$\displaystyle\leq\left\|\int_{0}^{1}\frac{1}{\left(1-y\lambda(\mathbf{x}_{k})\right)^{2}}dy\right\|_{2}\left\|\
[\nabla^{2}f({\mathbf{x}_{k}})]^{1/2}\left(\mathbf{\hat{d}}_{h,k}-\mathbf{d}_{h,k}\right)\right\|_{2}$
$\displaystyle=\frac{1}{1-\lambda(\mathbf{x}_{k})}\left\|[\nabla^{2}f({\mathbf{x}_{k}})]^{1/2}\left(\mathbf{\hat{d}}_{h,k}-\mathbf{d}_{h,k}\right)\right\|_{2}$
$\displaystyle\leq\frac{\sqrt{\lambda(\mathbf{x}_{k})^{2}-\hat{\lambda}(\mathbf{x}_{k})^{2}}}{1-\lambda(\mathbf{x}_{h,k})}.$

Next, the second norm implies

$\displaystyle\left\|\mathbf{Z}_{2}\right\|$
$\displaystyle\leq\left\|\int_{0}^{1}\left(\frac{1}{\left(1-y\lambda(\mathbf{x}_{k})\right)^{2}}-1\right)dy\right\|_{2}\left\|[\nabla^{2}f({\mathbf{x}_{k}})]^{1/2}\mathbf{d}_{h,k}\right\|_{2}$
$\displaystyle=\frac{\lambda(\mathbf{x}_{k})}{1-\lambda(\mathbf{x}_{k})}\lambda(\mathbf{x}_{k}).$

Putting this all together, inequality (40) becomes

$\lambda(\mathbf{x}_{k+1})\leq\frac{\sqrt{\lambda(\mathbf{x}_{k})^{2}-\hat{\lambda}(\mathbf{x}_{k})^{2}}}{\left(1-\lambda(\mathbf{x}_{k})\right)^{2}}+\frac{\lambda(\mathbf{x}_{k})}{\left(1-\lambda(\mathbf{x}_{k})\right)^{2}}\lambda(\mathbf{x}_{k}).$

Taking now expectation, conditioned on the $\sigma$-algebra $\mathcal{F}_{k}$
we have that

$\mathbb{E}[\lambda(\mathbf{x}_{k+1})|\mathcal{F}_{k}]\leq\frac{\mathbb{E}\left[\sqrt{\lambda(\mathbf{x}_{k})^{2}-\hat{\lambda}(\mathbf{x}_{k})^{2}}|\mathcal{F}_{k}\right]}{\left(1-\lambda(\mathbf{x}_{k})\right)^{2}}+\frac{\lambda(\mathbf{x}_{k})}{\left(1-\lambda(\mathbf{x}_{k})\right)^{2}}\lambda(\mathbf{x}_{k}).$
(42)

By the concavity of the square root and Jensen’s inequality we have that

$\displaystyle\mathbb{E}\left[\sqrt{\lambda(\mathbf{x}_{k})^{2}-\hat{\lambda}(\mathbf{x}_{k})^{2}}|\mathcal{F}_{k}\right]$
$\displaystyle\leq\sqrt{\lambda(\mathbf{x}_{k})^{2}-\mathbb{E}[\hat{\lambda}(\mathbf{x}_{k})|\mathcal{F}_{k}]^{2}}.$

Next, by (38), we have
$\hat{M}\lambda(\mathbf{x}_{k})\leq\mathbb{E}[\hat{\lambda}(\mathbf{x}_{k})|\mathcal{F}_{k}]$
and then

$\displaystyle\lambda(\mathbf{x}_{k})^{2}-\mathbb{E}[\hat{\lambda}(\mathbf{x}_{k})|\mathcal{F}_{k}]^{2}$
$\displaystyle\leq(1-\hat{M}^{2})\lambda(\mathbf{x}_{k})^{2}$

Putting this all together, inequality (42) becomes

$\mathbb{E}\left[\lambda(\mathbf{x}_{k+1})|\mathcal{F}_{k}\right]\leq\frac{\hat{\varepsilon}+\lambda(\mathbf{x}_{k})}{\left(1-\lambda(\mathbf{x}_{k})\right)^{2}}\lambda(\mathbf{x}_{k}),$

which holds almost surely as required. ∎

Now, we are ready to prove Theorem 3.4.

###### proof of Theorem 3.4.

We start by showing reduction in value function of phase (i). Taking
expectation conditioned on $\mathcal{F}_{k}$ in the result of Lemma A.13 we
have that

$\mathbb{E}[f(\mathbf{x}_{k+1})|\mathcal{F}_{k}]-f(\mathbf{x}_{k})\leq-\mathbb{E}[\omega(\hat{\lambda}(\mathbf{x}))|\mathcal{F}_{k}]\leq-\omega(\mathbb{E}[\hat{\lambda}(\mathbf{x})|\mathcal{F}_{k}]),$

where the last inequality holds by the convexity of $\omega$ and Jensen’s
inequality. Next, from (38) it holds that
$\mathbb{E}\left[\hat{\lambda}(\mathbf{x}_{h,k})|\mathcal{F}_{k}\right]\geq\hat{M}\lambda(\mathbf{x}_{h,k})$,
and since $\lambda(\mathbf{x}_{h,k})>\hat{\eta}$ by assumption we get
$\mathbb{E}\left[\hat{\lambda}(\mathbf{x}_{h,k})|\mathcal{F}_{k}\right]>\hat{\eta}\hat{M}$.
Since also $\omega(x)$ is monotone increasing

$\mathbb{E}[f(\mathbf{x}_{k+1})|\mathcal{F}_{k}]-f(\mathbf{x}_{k})\leq-\omega_{*}(\hat{\eta}\hat{M})$

which holds almost surely. Taking total expectation on both sides the claim
follows with $\gamma:=\mathbb{E}[\omega(\hat{\eta}\hat{M})]>0$.

For phase (ii), using Lemma A.14, the fact that
$x\rightarrow\frac{\hat{\varepsilon}+x}{(1-x)^{2}}$ is monotone increasing and
$\lambda(\mathbf{x}_{h,k})\leq\hat{\eta}$ we have that

$\mathbb{E}\left[\lambda(\mathbf{x}_{h,k+1})|\mathcal{F}_{k}\right]\leq\frac{\hat{\varepsilon}+\hat{\eta}}{\left(1-\hat{\eta}\right)^{2}}\lambda(\mathbf{x}_{h,k}).$

Next, by definition $\hat{\eta}:=\frac{3-\sqrt{5+4\hat{e}}}{2}$ and since
$0\leq\hat{e}<1$ we obtain
$\frac{\hat{e}+\hat{\eta}}{\left(1-\hat{\eta}\right)^{2}}<1$. Thus,
$\mathbb{E}\left[\lambda(\mathbf{x}_{h,k+1})|\mathcal{F}_{k}\right]<\lambda(\mathbf{x}_{h,k})$
for any $\hat{\eta}\in(0,\frac{3-\sqrt{5}}{2})$. Taking total expectation on
both sides the claim follows. Finally, by Markov’s inequality and Lemma A.12
we have that

$\operatorname{Prob}(f_{h}(\mathbf{x}_{h,k})-f_{h}(\mathbf{x}_{h}^{*})\geq\epsilon)\leq\frac{\mathbb{E}[f_{h}(\mathbf{x}_{h,k})-f_{h}(\mathbf{x}_{h}^{*})]}{\epsilon}\leq\frac{\mathbb{E}[\lambda(\mathbf{x}_{h,k})^{2}]}{\epsilon}.$

for some tolerance $\epsilon\in(0,0.68^{2})$. Since the right-hand side
converges to zero in the limit then convergence is attained with probability
one. ∎

###  A.4 Proof of Theorem 3.5

If a function is strongly convex, then there exists $m,M>0$ such that

$f(\mathbf{y})\leq f(\mathbf{x})+\nabla
f(\mathbf{x})^{T}(\mathbf{y}-\mathbf{x})+\frac{M}{2}\|\mathbf{y}-\mathbf{x}\|_{2}^{2},$
(43)

and,

$m\mathbf{I}_{n}\leq\nabla^{2}f(\mathbf{x})\leq M\mathbf{I}_{n}.$ (44)

As in Appendix A.3, we denote
$\mathbf{x}_{k+1}=\mathbf{x}_{k}-t_{k}\mathbf{P}\mathbf{Q}_{H,k}^{-1}\mathbf{R}\nabla
f(\mathbf{x}_{k})$,
$\hat{\mathbf{d}}_{h,k}:=-\mathbf{P}\mathbf{Q}_{H,k}^{-1}\mathbf{R}\nabla
f(\mathbf{x}_{k})$ and $\hat{\lambda}(\mathbf{x})$ as in (15). Then, by (43)
we have

$f(\mathbf{x}_{k+1})-f(\mathbf{x}_{k})\leq t_{k}\nabla
f(\mathbf{x})^{T}\hat{\mathbf{d}}_{h,k}+\frac{M}{2}t_{k}^{2}\|\hat{\mathbf{d}}_{h,k}\|_{2}^{2}.$
(45)

Recall that
$f(\mathbf{x})^{T}\hat{\mathbf{d}}_{h,k}=-\hat{\lambda}(\mathbf{x}_{h,k})^{2}$.
Further, combining (35) and (44) we take

$\hat{\lambda}(\mathbf{x}_{h,k})^{2}\geq\hat{\mathbf{d}}_{h,k}^{T}\nabla^{2}f(\mathbf{x})\hat{\mathbf{d}}_{h,k}\geq
m\|\hat{\mathbf{d}}_{h,k}\|_{2}^{2}.$

Therefore, inequality (45) becomes

$f(\mathbf{x}_{k}+t_{k}\hat{\mathbf{d}}_{h,k})-f(\mathbf{x}_{k})\leq-
t_{k}\hat{\lambda}(\mathbf{x}_{h,k})^{2}+\frac{Mt_{k}^{2}}{2m}\hat{\lambda}(\mathbf{x}_{h,k})^{2}$

Minimizing both sides of the above inequality over $t$ we have

$f(\mathbf{x}_{k}+t_{k}^{*}\hat{\mathbf{d}}_{h,k})-f(\mathbf{x}_{k})\leq-\frac{m}{2M}\hat{\lambda}(\mathbf{x}_{h,k})^{2}.$

Taking now expectation conditioned on $\mathcal{F}_{k}$ the above inequality
becomes

$\mathbb{E}[f(\mathbf{x}_{k+1})|\mathcal{F}_{k}]-f(\mathbf{x}_{k})\leq-\frac{m}{2M}\mathbb{E}[\hat{\lambda}(\mathbf{x}_{h,k})|\mathcal{F}_{k}]^{2}.$

Moreover, by (38) it holds
$\mathbb{E}[\hat{\lambda}(\mathbf{x}_{h,k})|\mathcal{F}_{k}]\geq\mathbb{E}[\mu_{k}\sqrt{\frac{\sigma_{k,N}}{\sigma_{k,p+1}}}|\mathcal{F}_{k}]\lambda(\mathbf{x}_{h,k})=\hat{M}_{k}\lambda(\mathbf{x}_{h,k})$
and therefore

$\mathbb{E}[f(\mathbf{x}_{k+1})|\mathcal{F}_{k}]-f(\mathbf{x}_{k})\leq-\frac{m}{2M}\hat{M}_{k}\lambda(\mathbf{x}_{h,k}),$

and taking total expectation the result follows.

###  A.5 Proof of Theorem 3.6

Let $\mathbf{x}_{k+1}=\mathbf{x}_{k}+t_{k}|\hat{\mathbf{d}}_{h,k}|$ and in
addition define $\hat{\lambda}(\mathbf{x})$ be the the approximate decrement
based on the current coarse model, i.e.,

$\hat{\lambda}(\mathbf{x}_{k}):=\left[\nabla
f(\mathbf{x}_{k})^{T}|\hat{\mathbf{Q}}_{h,k}^{-1}|\nabla
f(\mathbf{x}_{k})\right]^{1/2},$

where $|\hat{\mathbf{d}}_{h,k}|$ and $|\hat{\mathbf{Q}}_{h,k}^{-1}|$ are
defined in (17). By construction of the truncated Hessian matrix in (16) we
have that

$\frac{1}{|\sigma_{k,1}|}\mathbf{I}_{n}\preceq|\mathbf{Q}_{h,k}^{-1}|\preceq\frac{1}{|\sigma_{k,p+1}|}\mathbf{I}_{n}.$

Using the above inequality we obtain the following bounds for
$\|\hat{\mathbf{d}}_{h,k}\|_{2}$ and $\hat{\lambda}(\mathbf{x}_{k})$,

$\|\hat{\mathbf{d}}_{h,k}\|^{2}_{2}\leq\frac{1}{\sigma_{k,p+1}^{2}}\|\nabla
f(\mathbf{x}_{k})\|^{2}\ \ \text{and}\ \
\hat{\lambda}(\mathbf{x}_{k})^{2}\geq\frac{1}{|\sigma_{k,1}|}\|\nabla
f(\mathbf{x}_{k})\|^{2}.$

Using the Lipschitz continuity of the gradients (inequality 45) together with
the above inequalities we take

$\displaystyle f(\mathbf{x}_{k+1})-f(\mathbf{x}_{k})$ $\displaystyle\leq-
t_{k}\hat{\lambda}(\mathbf{x}_{k})^{2}+\frac{M}{2}t_{k}^{2}\|\hat{\mathbf{d}}_{h,k}\|_{2}^{2}$
$\displaystyle\leq\left(\frac{Mt_{k}}{2\sigma_{k,p+1}^{2}}-\frac{1}{|\sigma_{k,1}|}\right)t_{k}\|\nabla
f(\mathbf{x}_{k})\|^{2}.$

Letting $t_{k}:=\frac{\sigma_{k,p+1}^{2}}{M|\sigma_{k,1}|}$ the above
inequality becomes

$f(\mathbf{x}_{k+1})-f(\mathbf{x}_{k})\leq-\frac{\sigma_{k,p+1}^{2}}{2M\sigma_{k,1}^{2}}\|\nabla
f(\mathbf{x}_{k})\|^{2}.$

Finally, we obtain the required result by applying the PL inequality (18) on
the last inequality and then taking expectations on both sides.

##  Appendix B Extra Numerical Results and Details

All the datasets used in the experiments can be found in
https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/ and
http://archive.ics.uci.edu/ml/index.php. What follows is a description of the
algorithms used in comparisons.

  1. 1.

Gradient Descent (GD) with Armijo-rule

  2. 2.

Accelerated Gradient Descent (AccGD) with Armijo-rule and momentum $0.5$.

  3. 3.

NewSamp with with Armijo-rule Erdogdu and Montanari [2015]

  4. 4.

Adam with $t_{k}=10^{-3},\beta_{1}=0.9,\beta_{2}=0.99,\epsilon=10^{-8}$ Kingma
and Ba [2014]

  5. 5.

Sigma with Armijo-rule Tsipinakis and Parpas [2021]

  6. 6.

SigmaSVD with Armijo-rule

  7. 7.

Cubic Regularization of the Newton’s method (Cubic Newton) with line search
Nesterov et al. [2018]

To be fair in the comparisons we perform the same step size strategy for all
algorithms but Adam. Note that NewSamp forms the Hessian matrix as in (10) and
although it is not directly comparable to our approach, we include it in our
experiments to show that our approach outperforms sub-sampled Newton methods.
Similar to $N$, we denote $|S_{m}|$ the number of samples that NewSamp uses to
form the Hessian. In addition, the number $p$ of the eigenvalues is selected
the same for both SigmaSVD and NewSamp in all the experiments. The results of
this section were generated on standard laptop machine with Intel i7-10750H
CPU processor. The numerical results of main text were generated on GPU
processor in Google’s colab. The code for generating all the numerical results
will be uploaded on github in the near future.

###  B.1 Non-linear regression

Datasets | Problem | $m$ | $n$ | $N$ | $p$
---|---|---|---|---|---
CTslices | Non-linear least-squares | $53,500$ | $385$ | $0.5n$ | $60$
CovType | Non-linear least-squares | $581,012$ | $54$ | $0.7n$ | $5$
Gisette | Non-linear least-squares | $6,000$ | $5,000$ | $0.5n$ | $350$
W8T | Non-linear least-squares | $14,951$ | $300$ | $0.5n$ | $60$
Table 2: Set-up and dataset details for non-convex, non-linear regression
prblem.

In this section we revisit the non-convex problem of Section 4 to perform
simulations on different datasets. Full details on the datasets used for the
non-linear least-squares problem and how the parameters $N$ and $p$ are
selected are given in Table 2. Figure 2 shows that SigmaSVD, Adam and Cubic
Newton are the best algorithms as they return the lowest training errors, yet
SigmaSVD does better in problems with several saddle points or flat areas. The
three algorithms perform similar in terms of GPU time when the problem
dimensions is small. However, for problems with $n$ large the results are
favourable for SigmaSVD.

(a) CovType

(b) CtSlices

(c) w8t

(d) CovType

(e) CtSlices

(f) w8t

Figure 2: Non-convex minimization. All methods in plots from (a) to (c) are
initialized at the origin while from (e) to (h) the initializer is selected
randomly from A Gaussian $\mathcal{N}(0,1)$.

###  B.2 Log-linear regression

Here we minimize the following self-concordant function,

$\mathbf{x}^{*}=\underset{\mathbf{x}\in{\mathbb{R}}^{n}}{\operatorname{arg\
min}}f(\mathbf{x})=\underset{\mathbf{x}\in{\mathbb{R}}^{n}}{\operatorname{arg\
min}}\left(-\sum_{i=1}^{m}\log(b_{i}-\mathbf{a}_{i}^{T}\mathbf{x})\right).$

We generate two datasets to illustrate the efficiency of the proposed method
for the regimes $m>n$ and $n>m$ as follows: $\\{\mathbf{a_{i}}\\}_{i=1}^{m}$
is generated from the multivariate Gaussian distribution with zero mean and
unit variance and $\\{b_{i}\\}_{i=1}^{m}$ from uniform distribution. Full
details and the set-up for algorithm 1 are given in Table 3.

Datasets | Problem | $m$ | $n$ | $N$ | $|S_{m}|$ | $p$
---|---|---|---|---|---|---
Synthetic | Log linear model | $10,000$ | $1,000$ | $0.5n$ | $0.3m$ | $150$
Synthetic | Log linear | $1,000$ | $10,000$ | $0.1n$ | - | $150$
Table 3: Set-up and dataset details for Log-linear regression.

(a) $m>n$

(b) $m>n$

(c) $n>m$

(d) $n>m$

Figure 3: Log Linear Regression. Plots (a) and (c) show comparisons between
the optimization algorithms for the regime $m>n$ while (b) and (d) for the
regime $m<n$.

Figure 3 shows the performance of the optimization algorithms. Due to the
domain of this problem, the performance of Adam is missing. The performance of
NewSamp is missing from the experiment for the regime $n>m$ due to the large
value of $n$. Clearly, the fastest method is Sigma which achieves a very fast
super-linear convergence rate. As shown in Theorem 3.4, the reason for the
difference in the performance between Sigma and SigmaSVD is that the ratio
$\sigma_{p+1}/\sigma_{n}$ is large which indicates that performing SVD on the
reduce Hessian matrix we discard much of the important second-order
information of the problem. However, even in such a scenario, in contrast to
NewSamp and gradient descent methods, SigmaSVD is able to converge to the
global minimum with a super-linear rate.

###  B.3 Logistic regression

In this section we are concerned with the problem of finding the maximum
likelihood in generalized linear models. In particular, we consider the
regularized logistic regression which as an optimization problem takes the
following form,

$\mathbf{x}^{*}=\underset{\mathbf{x}\in{\mathbb{R}}^{n}}{\operatorname{arg\
min}}f(\mathbf{x})=\underset{\mathbf{x}\in{\mathbb{R}}^{n}}{\operatorname{arg\
min}}\left(\frac{1}{m}\sum_{i=1}^{m}\log(1+e^{-b_{i}\mathbf{a}_{i}^{T}\mathbf{x}_{h}})+\frac{\ell}{2}\|\mathbf{x}\|_{2}^{2}\right).$

Note that the logistic model is a strongly convex function with Lipschitz
continuous Hessian matrix. The details and results for the logistic regression
experiment are given in Table 4 and Figure 4, respectively. Note that for
leukemia dataset, $m$ is too small to perform batch learning and thus the
performance of Adam is omitted for this example. Second-order methods clearly
outperform first-order methods for the CtSlices dataset. Further, NewSamp is
slightly faster than the multilevel methods but obtains slightly worse
convergence rate (Figures 4(a)) and 4(d), respectively). However, as expected,
the efficiency of NewSamp reduces drastically for large values of $n$
(Leukemia dataset). Note that in this example the multilevel methods perform
similarly which indicates that there is no need to use the full spectrum of
the reduced Hessian matrix. Last, Figures 4(c) and 4(f) show the efficiency of
SigmaSVD on a problem with over a million parameters. Such a problem lies at
the heart of large scale machine learning and deep learning. Even when the
coarse dimensions are selected very small ($N=0.001n$ in this experiment),
Figure 4(c) shows that multilevel methods are capable of returning solution
with very good accuracy much faster than first-order methods.

Datasets | Problem | $m$ | $n$ | $N$ | $|S_{m}|$ | $p$ | $\ell$
---|---|---|---|---|---|---|---
CTslices | Logistic model | $53,500$ | $385$ | $0.5n$ | $0.3m$ | $60$ | $10^{-6}$
Leukemia | Logistic model | $38$ | $7,129$ | $0.1n$ | $0.7m$ | $150$ | $10^{-6}$
News20 | Logistic model | $19,996$ | $1,355,129$ | $0.001n$ | $-$ | $100$ | $10^{-12}$
Table 4: Set-up and dataset details for Logistic regression.

(a) CtSlices

(b) Leukemia

(c) News20

(d) CtSlices

(e) Leukemia

(f) News20

Figure 4: Logistic Regression. Plots (a) to (c) show the norm of the gradient
vs cpu time in seconds while (d) to (f) norm of the gradient vs iterations for
three machine learning datasets.

###  B.4 Support Vector Machines

We can train Support Vector Machines (SVMs) using the primal problem over
hinge-q loss or Huber loss function. Since our approach requires twice
differentiable functions we consider the hinge-2 loss function for training
the SVMs,

$\mathbf{x}^{*}=\underset{\mathbf{x}\in{\mathbb{R}}^{n}}{\operatorname{arg\
min}}\frac{1}{2}\|\mathbf{x}\|_{2}^{2}+\frac{\ell}{2}\sum_{i=1}^{m}\max\\{0,1-b_{i}\mathbf{a}_{i}^{T}\mathbf{x}\\}^{2}.$

Note that the objective function of this section is convex, however the
Hessian matrix is not Lipschitz continuous. The details and results for the
logistic regression experiment are given in Table 5 and Figure 5,
respectively. Again, both multilevel methods outperform its counterparts.
Here, as also in the logistic regression example, using only a few eigenvalues
to form the reduced Hessian matrix does not seem to affect the performance of
SigmaSVD. Further, we observed that the performance of NewSamp becomes very
poor, even compared to first-order methods, when highly regularized solutions
are required (Figures 5(b) and 5(d)). On the other hand, both multilevel
methods, no matter how large or small the regularization parameter is chosen,
always outperform the first-order methods.

Datasets | Problem | $m$ | $n$ | $N$ | $|S_{m}|$ | $p$ | $\ell$
---|---|---|---|---|---|---|---
CovType | SVM | $581,012$ | $54$ | $0.5n$ | $0.3m$ | $10$ | $10^{-2}$
W8T | SVM | $14,951$ | $300$ | $0.5n$ | $0.5m$ | $60$ | $10^{3}$
Table 5: Set-up and dataset details for Logistic regression.

(a) W8T

(b) CovType

(c) W8T

(d) CovType

Figure 5: Support Vector Machines. Plots (a) and (b) show the norm of the
gradient vs cpu time in seconds while (c) and (d) norm of the gradient vs
iterations for two machine learning datasets.