###### Abstract

Despite the impressive numerical performance of the quasi-Newton and
Anderson/nonlinear acceleration methods, their global convergence rates have
remained elusive for over 50 years. This study addresses this long-standing
issue by introducing a framework that derives novel, adaptive quasi-Newton and
nonlinear/Anderson acceleration schemes. Under mild assumptions, the proposed
iterative methods exhibit explicit, non-asymptotic convergence rates that
blend those of the gradient descent and Cubic Regularized Newton’s methods.
The proposed approach also includes an accelerated version for convex
functions. Notably, these rates are achieved adaptively without prior
knowledge of the function’s parameters. The framework presented in this study
is generic, and its special cases includes algorithms such as Newton’s method
with random subspaces, finite-differences, or lazy Hessian. Numerical
experiments demonstrated the efficiency of the proposed framework, even
compared to the l-BFGS algorithm with Wolfe line-search.

Adaptive Quasi-Newton and Anderson Acceleration Framework with Explicit Global
(Accelerated) Convergence Rates

Damien Scieur Samsung SAIL Montreal

##  1 Introduction

Consider the problem of determining the minimizer $x^{\star}$ of the
unconstrained minimization problem

$f^{\star}\stackrel{{\scriptstyle\text{def}}}{{=}}f(x^{\star})=\min_{x\in\mathbb{R}^{d}}f(x),$

where $d$ is the problem dimension, and the function $f$ has a Lipschitz
continuous Hessian.

######  Assumption 1.

The function $f(x)$ has a Lipschitz continuous Hessian with constant $L$:

$\forall\;\;y,\,z\in\mathbb{R}^{d},\quad\|\nabla^{2}f(z)-\nabla^{2}f(y)\|\leq
L\|z-y\|.$ (1)

In this study, $\|.\|$ stands for the maximal singular value of a matrix and
the $\ell_{2}$ norm for a vector. Numerous twice-differentiable problems, such
as logistic or least-squares regression, satisfy Assumption 1.

The Lipschitz continuity of the Hessian is crucial when analyzing second-order
algorithms because it extends the concept of smoothness to second-order. The
groundbreaking work by [61] revealed the remarkable convergence rate of
Newton’s method for problems satisfying Assumption 1 when augmented with cubic
regularization. For instance, if the problem is convex, the accelerated
gradient descent typically achieves $O(\frac{1}{t^{2}})$, whereas accelerated
second-order methods achieve $O(\frac{1}{t^{3}})$. Recent advancements
achieved even faster convergence rates, up to $\mathcal{O}(\frac{1}{t^{7/2}})$
using hybrid methods [57, 17], or the direct acceleration of second-order
methods [59, 34, 53].

However, second-order methods are not scalable, particularly for the high-
dimensional problems common in machine learning. The limitation is that exact
second-order methods require the solution of a linear system involving the
Hessian $\nabla^{2}f$. This has motivated alternative approaches that balance
the efficiency of second-order methods with the scalability of first-order
methods, such as quasi-Newton methods or Anderson/nonlinear acceleration
methods (which are equivalent to quasi-Newton methods [29]). Due to space
limitation, the results of this study are presented under the prism of quasi-
Newton methods, but their links with Anderson acceleration is explained
extensively in appendix D.

Quasi-Newton (qN) methods minimize differentiable functions by iteratively
updating an approximate Hessian matrix using previous gradients, effectively
balancing scalability and efficiency. This approach makes them highly suitable
for large-scale optimization problems across diverse fields. For instance,
l-BFGS is a widely used and effective optimization method for unconstrained
functions (e.g., fminunc from Matlab) and is often considered state-of-the-art
[1].

Despite the impressive performance of quasi-Newton methods and nonlinear
acceleration schemes, the following long-standing question has remained
unanswered for over 50 years.

What are the nonasymptotic global convergence rates of quasi-Newton and
Anderson acceleration methods?

This question is challenging: Over the years, extensive research has
catapulted the popular l-BFGS algorithm to an exceptionally high level of
efficiency, as attested by various studies (e.g., [56, 89, 58]). However, its
theoretical convergence guarantees are notably lacking and do not accurately
reflect its actual numerical performance. Therefore, additional numerical
improvements or obtaining fast rates without arming the numerical convergence
may be increasingly difficult or infeasible.

###  1.1 Contributions

##### Theoretical guarantees

This study presents generic updates that are novel quasi-Newton methods or
nonlinear acceleration schemes with cubic regularization that meet the
following requirements.

  1. 1.

The assumptions for the theoretical analysis are simple and verifiable
(section 3.1).

  2. 2.

The algorithm is suitable for large-scale problems, as for a fixed memory
budget $N\leq d$, its per-iteration cost is linear in the dimension.

  3. 3.

The algorithm exhibits explicit, global, and nonasymptotic convergence rates
that interpolate between those of the first and second order methods (sections
3 and C):

     * •

Convergence on Nonconvex problems (Theorem 5): $\min_{i\leq t}\|\nabla
f(x_{i})\|\leq O(t^{-\frac{2}{3}}+t^{-\frac{1}{3}})$,

     * •

(Star-)convex problems (Theorems 6 and 7): $f(x_{t})-f^{\star}\leq
O(t^{-2}+t^{-1})$,

     * •

Accelerated rate on convex problems (Theorem 8): $f(x_{t})-f^{\star}\leq
O(t^{-3}+t^{-2})$.

  4. 4.

The algorithm is adaptive to the problem’s constants (algorithms 4 and 7):
both accelerated and classical methods require only an initial estimate of the
Lipchitz constant (section F.1).

##### Novel Analysis

To address these points, this study explores a novel paradigm, rethinking from
scratch the framework underlying qN methods (section 2).

##### Numerical Efficiency

The algorithm outperform l-BFGS in many scenarios (section 4).

##### Practicability

A particular focus have been put on making the method simple, generic and
adaptive, to make it suitable for practical applications. The method is simple
to implement, and requires fewer hyperparameters than classical qN schemes.

##### Generic Framework

The framework can be instantiated as many kind of previously known methods
(section 3.3.2), and recovers the cubic regularization of Newton’s method in
its most extreme case.

###  1.2 Limitations in Current qN Schemes

In most classical qN methods, a (Wolfe) line search algorithm (often in
addition to other techniques such as secant equation filtering or re-scaling)
is required to ensure global convergence. Without such a line search, the
performance of qN methods is poor, even in a simple quadratic case in two
dimensions [64].

Nevertheless, some previous work already attempted to determine rates for qN
methods (or to derive new ones), but often violates at least one of the
previous points: 1) the analysis requires non-verifiable assumptions, 2) the
algorithm is not suitable for large-scale problems as the per-iteration cost
is at least $O(d^{2})$, 3) the rates are local or do not interpolate between
first and second order rates, 4) the algorithm depends on potentially unknown
parameters. See appendix B for a detailed literature review.

Violates 1: For instance, the ARC method [18, 19] or proximal qN methods [92,
68] show accelerated rates for quasi-Newton under similar assumptions.
However, the authors state that the convergence rate is derived under a
nonverifiable assumption and typically their rates do not rely on or exploit
the accuracy of second-order approximations.

Violates 2: By using online algorithms and the Monteiro-Svaiter acceleration
technique, [50] achieves accelerated rates
$O(\min\\{\frac{1}{t^{2}},\frac{1}{t^{2.5}}\\})$ for qN methods, but require
the storage and inversion of a $d\times d$ matrix, which does not scale well
in high dimension, and does not perform well compared to BFGS.

Violates 3: Recent research on quasi-Newton updates has unveiled explicit and
nonasymptotic rates of convergence [65, 67, 66, 54, 55]. Nonetheless, these
analyses suffer from several significant drawbacks, as they are local, full-
memory (hence require storing a $d\times d$ matrix and the per-iteration cost
is $O(d^{2})$) and sometimes requires access to the Hessian matrix.

Violates 4: [52] introduced an adaptive regularization technique combined with
cubic regularization, but the method relies on the knowledge of $L$ in
Assumption 1.

##  2 Rethinking Quasi-Newton Methods

The starting point is the cubic upper bound on the objective function $f$ and
the upper bound on the gradient variation derived from Assumption 1 [61]:

$\displaystyle\textstyle\|\nabla f(y)-\nabla
f(x)-\nabla^{2}f(x)(y-x)\|\leq\frac{L}{2}\|y-x\|^{2},$ (2)
$\displaystyle\textstyle f(y)\leq f(x)+\nabla f(x)(y-x)$
$\displaystyle\qquad+\frac{1}{2}(y-x)^{T}\nabla^{2}f(x)(y-x)+\frac{L}{6}\|y-x\|^{3},$
(3)

which holds for all $x,y\in\mathbb{R}^{d}$. Minimizing eq. 3 over $y$ yields
the cubic regularization of Newton’s method [61].

The main steps in deriving this algorithm are as follows: 1) The minimization
will be constrained to a subspace of dimension $N\leq d$, reducing the per-
iteration computation cost. 2) The Hessian (in the subspace) will be
approximated using differences of gradients. 3) From the previous points, an
upper bound for the objective function and the gradient norm will be
constructed, leading to a Type I and Type II method. 4) To ensure convergence,
the update of the subspace will be chosen such that it spans the gradient
(deterministic) or spans a portion of the gradient in expectation (random
subspace).

###  2.1 First Ingredient: Subspace Minimization

Minimizing the upper bound eq. 3 is expensive in high dimensions because it
requires an eigenvalue decomposition of the Hessian $\nabla^{2}f(x)$ [61].
Rather, let $D_{t}$ be an $d\times N$ matrix of directions (the construction
of $D_{t}$ is described in section 2.4). By constraining the update
$x_{t+1}-x_{t}$ in the columns span of $D_{t}$, that is,

$x_{t+1}=x_{t}+D_{t}\alpha_{t},$ (4)

where $\alpha_{t}$ is a vector of $N$ coefficients, the minimization problem
is simplified to

$\displaystyle\alpha_{t}=\operatorname*{{arg\,min}}_{\alpha\in\mathbb{R}^{N}}$
$\displaystyle\;\;f(x_{t})+\nabla f(x_{t})D_{t}\alpha$
$\displaystyle\;\;\textstyle+\frac{1}{2}(D_{t}\alpha)^{T}\nabla^{2}f(x_{t})D_{t}\alpha+\frac{L}{6}\|D_{t}\alpha\|^{3}.$

Hence, the complexity of the eigenvalue decomposition of
$D_{t}^{T}\nabla^{2}f(x)D_{t}$ is now $O(N^{2}d+N^{3})$.

###  2.2 Second Ingredient: Multisecant Approximation of the Hessian

Typically, quasi-Newton methods approximate the Hessian using the properties
of the secant equation:

$\nabla^{2}f(x_{i})(x_{i}-x_{i-1})\approx\nabla f(x_{i})-\nabla f(x_{i-1}).$

Usually, those updates are performed recursively, that is, by updating an
approximation of the Hessian one secant equation at a time.

This study approximates the Hessian using the secant equations simultaneously
by using a finite-difference approximation between two sets of points $Y_{t}$
and $Z_{t}$.

More formally, let the matrix of directions $D_{t}$ and the matrix of
normalized gradient differences $G_{t}$ be

$\displaystyle D_{t}$
$\displaystyle=\left[\frac{y^{(t)}_{1}-z^{(t)}_{1}}{\|y^{(t)}_{1}-z^{(t)}_{1}\|_{2}},\ldots,\frac{y^{(t)}_{N}-z^{(t)}_{N}}{\|y^{(t)}_{N}-z^{(t)}_{N}\|_{2}}\right],$
(5) $\displaystyle G_{t}$ $\displaystyle=\left[\ldots,\frac{\nabla
f(y^{(t)}_{i})-\nabla
f(z^{(t)}_{i})}{\|y^{(t)}_{i}-z^{(t)}_{i}\|_{2}},\ldots\right].$

and let the matrices $Y_{t}$, $Z_{t}$ be defined as:

$\displaystyle Y_{t}=[y_{1}^{(t)},\ldots,y_{N}^{(t)}],\quad
Z_{t}=[z_{1}^{(t)},\ldots,z_{N}^{(t)}].$ (6)

A naïve approach can be $Y_{t}=[x_{t-N},\ldots,x_{t-1}]$ and
$Z_{t}=[x_{t-N+1},\ldots,x_{t}]$ (this will not be the case in this study; see
section 2.4). Intuitively, the matrix $G_{t}$ is a finite-difference
approximation of the Hessian matrix product $\nabla^{2}f(x_{t})D_{t}$. More
precisely, the next theorem states a bound on the approximation error of this
product as a function of the error vector $\varepsilon_{t}$:

$\displaystyle\varepsilon_{t}\stackrel{{\scriptstyle\text{def}}}{{=}}[e^{(t)}_{1},\ldots,e^{(t)}_{N}],\;e^{(t)}_{i}\stackrel{{\scriptstyle\text{def}}}{{=}}\|y_{i}^{(t)}-z_{i}^{(t)}\|+2\|z^{(t)}_{i}-x_{t}\|$
(7)

######  Theorem 1.

Let the function $f$ satisfy Assumption 1. Let matrices $D_{t},\,G_{t}$ be
defined as in eq. 6 and vector $\varepsilon$ as in eq. 7. Subsequently, for
all $w\in\mathbb{R}^{d}$ and $\alpha\in\mathbb{R}^{N}$

$\displaystyle\left|w^{T}(\nabla^{2}f(x_{t})D_{t}-G_{t})\alpha\right|$
$\displaystyle\leq\textstyle\frac{L\|w\|}{2}|\alpha|^{T}\varepsilon_{t},$ (8)
$\displaystyle\|w^{T}(\nabla^{2}f(x_{t})D_{t}-G_{t})\|$
$\displaystyle\leq\textstyle\frac{L\|w\|}{2}\|\varepsilon_{t}\|.$ (9)

##### Proof sketch

The detailed proof can be found in appendix G. The main idea of the proof is
as follows. The superscript$~{}^{(t)}$ has been removed for clarity. From eq.
2 with $y=y_{i}$, $z=z_{i}$, then Assumption 1,

$\displaystyle\frac{\|\nabla f(y_{i})-\nabla
f(z_{i})-\nabla^{2}f(x_{t})(y_{i}-z_{i})\|}{\|y_{i}-z_{i}\|}$
$\displaystyle\leq\frac{L}{2}\|y_{i}-z_{i}\|+\|\nabla^{2}f(x_{t})-\nabla^{2}f(z_{i})\|\leq\frac{L}{2}e_{i}.$

The first term in $e_{i}$ bounds the error of eq. 2, whereas the second term
originates from the distance between eq. 2 and the current point $x_{t}$ where
the Hessian is estimated. Subsequently, it suffices to combine the
inequalities with the coefficients $\alpha_{t}$ to obtain Theorem 1.

###  2.3 Third Ingredient: Objective Function and Gradient Norm Upper bounds

As the approximation error between $\nabla^{2}f(x_{t})D_{t}$ and $G_{t}$ can
be explicitly bounded, carefully replacing the term
$\nabla^{2}f(x_{t})D_{t}\alpha_{t}$ in eqs. 3 and 2 with $G_{t}\alpha_{t}$,
along with an appropriate regularization, leads to Type I and Type II bounds.

######  Theorem 2.

Let the function $f$ satisfy Assumption 1. Let $x_{t+1}$ be defined as in eq.
4; $D_{t},\,G_{t}$ be defined as in eq. 6 and $\varepsilon_{t}$ be defined as
in eq. 7. Subsequently, $\forall\alpha\in\mathbb{R}^{N}$,

$\displaystyle\textstyle f(x_{t+1})\leq f(x_{t})+\nabla
f(x_{t})^{T}D_{t}\alpha+\frac{\alpha^{T}H_{t}\alpha}{2}+\frac{L\|D_{t}\alpha\|^{3}}{6},$
(Type I bound) $\displaystyle\textstyle\|\nabla f(x_{t+1})\|\leq\left\|\nabla
f(x_{t})+G_{t}\alpha\right\|+\frac{L}{2}\Big{(}|\alpha|^{T}\varepsilon_{t}+\|D_{t}\alpha\|^{2}\Big{)},$
(Type II bound)

where
$H_{t}\stackrel{{\scriptstyle\text{def}}}{{=}}(G_{t}^{T}D_{t}+D_{t}^{T}G_{t}+\mathrm{I}L\|D_{t}\|\|\varepsilon_{t}\|)/2$.

The proof can be found in appendix G. Minimizing eqs. Type I bound and Type II
bound leads to algorithms 4 and 5, respectively, whose constant $L$ is
replaced by a parameter $M$ determined by a backtracking line search. Type I
methods often refer to algorithms that aim to minimize the function value
$f(x)$, whereas type II methods minimize the gradient norm $\|\nabla f(x)\|$
[29, 95, 14].

##### Solving the sub-problems

In algorithms 4 and 5, the coefficients $\alpha$ are computed by solving a
minimization sub-problem in $O(N^{3}+Nd)$, where in practice $N$ is smaller
than $d$ (see all the details in appendix E). In algorithm 4, the subproblem
can be formulated as a convex problem in two variables, by using an eigenvalue
decomposition of the matrix $H\in\mathbb{R}^{N\times N}$ [61], while for
algorithm 5, the subproblem can be cast into a linear-quadratic problem of
$O(N)$ variables and constraints that can be solved efficiently with SDP
solvers.

###  2.4 Fourth Ingredient: Direction Update Rule

Algorithm 1 "Orthogonal forward estimate only"

1:First-order oracle for $f$, step-size $h$, matrices $D_{t-1}$, $G_{t-1}$,
$Y_{t-1}$, $Z_{t-1}$, new point $x_{t}$.

2:If # columns of $D_{t-1}$, $G_{t-1}$, $Y_{t-1}$, $Z_{t-1}$ is larger than
$N$, then remove their first column.

3:Compute $g_{t}=\nabla f(x_{t})$, then compute
$d_{t}=-\frac{\tilde{d}}{\|\tilde{d}\|}$, where
$\tilde{d}=g_{t}-D_{t-1}(D_{t-1}^{T}g_{t})$.

4:Compute $x_{t+\frac{1}{2}}=x_{t}+hd_{t}$, the orthogonal forward estimate.

5:Update $Y_{t}=[Y_{t-1},x_{t+\frac{1}{2}}]$, $Z_{t}=[Z_{t-1},x_{t}]$,
$D_{t}=[D_{t-1},d_{t}]$, $G_{t}=\lx@cref{creftype~refnum}{def:matrices_2}$,
$\varepsilon=\lx@cref{creftype~refnum}{eq:error_vector}$ .

6:return $\nabla f(x_{t})$, $D_{t}$, $G_{t}$, $Y_{t}$, $Z_{t}$,
$\varepsilon_{t}$.

Algorithm 2 "Orthogonal random dir." (example)

1:First-order oracle for $f$, step-size $h$, memory $N$, new point $x_{t}$.

2:Generates $N$ random orthonormal directions, for example,
$[D_{t},]=\texttt{qr}(\texttt{Rand}(d,N))$.

3:Create $Z_{t}=[x_{t},\ldots,x_{t}]$, $Y_{t}=Z_{t}+hD_{t}$, then update
$G_{t}=\lx@cref{creftype~refnum}{def:matrices_2}$, $\varepsilon=[h,\ldots,h]$
.

4:return $\nabla f(x_{t})$, $D_{t}$, $G_{t}$, $Y_{t}$, $Z_{t}$,
$\varepsilon_{t}$.

One critical theoretical property in the analysis is the alignment of the
gradient $\nabla f(x_{t})$ with the directions in $D_{t}$. This section
presents updates that ensure good theoretical properties of $D_{t}$ (see
summary in table 1).

Below are some assumptions for updating $Y_{t},\,Z_{t},\,D_{t}$, which are
called requirements. While not overly restrictive, naïve methods, such as
keeping only the previous iterations, will not satisfy these requirements.

All convergence results rely on _one_ of these conditions on the projector
onto $\textbf{span}(D_{t})$,

$\displaystyle\textstyle
P_{t}\stackrel{{\scriptstyle\text{def}}}{{=}}D_{t}(D_{t}^{T}D_{t})^{-1}D_{t}^{T}.$
(10)

######  Requirement 1a.

For all $t$, projector $P_{t}$ of stochastic matrix $D_{t}$ satisfies
$\mathbb{E}[P_{t}]=\frac{N}{d}\textbf{I}$.

######  Requirement 1b.

For all $t$, $P_{t}\nabla f(x_{t})=\nabla f(x_{t})$.

The first condition guarantees that $D_{t}$ partially spans the gradient
$\nabla f(x_{t})$ in expectation because $\mathbb{E}[P_{t}\nabla
f(x_{t})]=\frac{N}{d}\nabla f(x_{t})$. The second condition requires the
possibility of moving toward the current gradient when taking the step
$x_{t}+D_{t}\alpha_{t}$.

The norm of the relative error must be bounded.

######  Requirement 2.

For all $t$, $(\|\varepsilon_{t}\|/\|D_{t}\|)\leq\delta$.

The Requirement 2 is also nonrestrictive, as it simply prevents the secant
equations at $y_{i}-z_{i}$ and $z_{i}-x_{i}$ from diverging significantly.
Generally, $\delta$ satisfies the (very) crude bound $\delta\leq
O(\|x_{0}-x^{\star}\|)$.

Finally, the condition number of the matrix $D$ must also be bounded. The
following section provides explicit updates for ensuring that this condition
is satisfied.

######  Requirement 3.

For all $t$, the condition number
$\kappa_{D_{t}}\stackrel{{\scriptstyle\text{def}}}{{=}}\sqrt{\|D_{t}^{T}D_{t}\|\|(D_{t}^{T}D_{t})^{-1}\|}$
is bounded as $\kappa_{D_{t}}\leq\kappa$.

| Complexity | # Grads | $\kappa$ | $\delta$
---|---|---|---|---
Forw. est. | $O(Nd)$ | 2 | 1 | $O(R)$
Random | At least $O(d)$ | $N+1$ | 1 | $O(h)$
Pruning | $O(Nd+N^{2})$ | 2 | ? | $O(R)$
Orthogonal. | $O(N^{2}d)$ | $N+1$ | 1 | $O(h)$
Table 1: Comparison between different updates rules: "forward estimates only,"
"orthogonal random directions," pruning or orthogonalization. If the gradient
computation is costly or if $N$ is large, the forward estimate is probably the
best method as it only requires the computation of two new gradients, $\nabla
f(x_{t})$ and $\nabla f(x_{t+\frac{1}{2}})$. Otherwise, the orthogonalization
or random directions might be the methods of choice given their small
constants $\delta$ and $\kappa$.

####  2.4.1 "Orthogonal Forward Estimate Only" Update Rule (Recommended)

The "orthogonal forward estimate only" update maintains $D_{t}$ orthonormal,
that is, $D_{t}^{T}D_{t}=I$ for all $t$, while ensuring that $\nabla f(x_{t})$
belongs to the span of columns of $D_{t}$ (see algorithm 1). These conditions
are ensured thanks to an intermediate iterate $x_{t+\frac{1}{2}}$ which is
used to estimate $\nabla^{2}f(x_{t})\nabla f(x_{t})$; this is called the
orthogonal forward estimate:

$x_{t+\frac{1}{2}}=x_{t}-hd_{t},\;d_{t}=\frac{\nabla
f(x_{t})-\tilde{D}_{t-1}(\tilde{D}_{t-1}^{T}\nabla f(x_{t}))}{\|\nabla
f(x_{t})-\tilde{D}_{t-1}(\tilde{D}_{t-1}^{T}\nabla f(x_{t}))\|},$

where $h>0$ is a small step size and $\tilde{D}_{t-1}$ is simply the matrix
$D_{t-1}$ whose first column has been removed if its number of columns equals
$N$. This estimate is a gradient descent step projected onto the orthogonal
space of $\textbf{span}(\tilde{D}_{t-1})$, and is inexpensive because the
orthogonality of $D_{t}$ is maintained over the iterations.

After computing the forward estimate, the matrices $Y_{t}$ and $Z_{t}$ are
updated as the moving history of the previous forward iterates and previous
iterates:

$Y_{t}=[x_{t-N+\frac{3}{2}},\ldots,x_{t+\frac{1}{2}}],\qquad
Z_{t}=[x_{t-N+1},\ldots,x_{t}].$

Matrices $D_{t}$ and $G_{t}$ then follow eq. 5. See algorithm 1 for a detailed
implementation. This method has several advantages; it ensures good
theoretical performance, particularly because $\kappa=1$ (see Theorem 3), at
the cost of only one extra gradient evaluation.

######  Theorem 3.

The “orthogonal forward estimate only” update described in algorithm 1
satisfies Requirements 3 and 1b with $\kappa=1$.

####  2.4.2 "Random Orthogonal Directions"

The "random orthogonal direction" update (see algorithm 2) generates a batch
of $N$ random orthogonal directions at each iteration, such that

$\mathbb{E}[D_{t}D_{t}^{T}]=\frac{N}{d}I.$

Subsequently, it remains to update the matrices $Y_{t},\,Z_{t},G_{t}$,

$Z_{t}=[\underbrace{x_{t},\ldots,x_{t}}_{N\text{ times}}],\quad
Y_{t}=Z_{t}+hD_{t},\quad G_{t}=\lx@cref{creftype~refnum}{def:matrices_2}.$

For instance, $D_{t}$ could be the $Q$ matrix of a qr decomposition of a
random $N\times d$ matrix (complexity: $O(N^{2}d)$), or even simpler, an
aggregation of random canonical vectors (e.g., see [44], where the complexity
is $O(Nd)$). The major advantages of this approach are that $\kappa=1$ and
$\delta=\sqrt{N}\cdot h$. However, $N$ additional gradient calls are required
to computes $G_{t}$.

####  2.4.3 Other Matrix Updates: Pruning or Orthogonalization

It is possible to create other matrix updates, e.g., the iterates only (stores
the last forward estimate and previous iterates) or greedy (stores all
previous forward estimates and iterates) strategies, as detailed below:

$\displaystyle\begin{cases}Y_{t}&=[x_{t+\frac{1}{2}},x_{t},x_{t-1},\ldots,x_{t-N+2}],\\\
Z_{t}&=[x_{t},x_{t-1},x_{t-2}\ldots,x_{t-N+1}]\end{cases}$ (Iterates only)
$\displaystyle\begin{cases}Y_{t}&=[x_{t+\frac{1}{2}},x_{t},x_{t-\frac{1}{2}},\ldots,x_{t-\frac{N+2}{2}}],\\\
Z_{t}&=[x_{t},x_{t-\frac{1}{2}},x_{t-1},\ldots,x_{t-\frac{N+1}{2}}]\end{cases}$
(Greedy)

where, this time,
$x_{t+\frac{1}{2}}\stackrel{{\scriptstyle\text{def}}}{{=}}x_{t}-h\nabla
f(x_{t})$, with $h$ small. However, the directions in $D_{t}$ are not
orthogonal, hence, $\kappa$ in Requirement 3 may be large [88, 72].
Nevertheless, the condition number can be controlled via pruning or
orthogonalization.

##### Pruning.

It is sufficient to verify the condition number of $D_{t}$ and then prune the
columns of $Y_{t},\,Z_{t},\,D_{t},\,\text{and }G_{t}$ until $\kappa_{D_{t}}$
is sufficiently small, for instance, until $\kappa_{D_{t}}\leq 10^{3}$. The
advantage is that this method requires only one extra gradient computation
$\nabla f(x_{t+\frac{1}{2}})$ to construct $G_{t}$.

##### Orthogonalization

From the matrices $Y_{t},\,Z_{t}$, the matrix $D_{t}$ is computed as
$D_{t}=\texttt{qr}(Z_{t}-Y_{t})$. Subsequently, the rest of the procedure
follows steps $2$ and $3$ from algorithm 2. It ensure the orthogonality of
$D_{t}$, but requires $N$ extra gradients to evaluate $G_{t}$.

The pruning strategy is more cost-effective than orthogonalization but
sacrifices its control over the history size. The orthogonalization technique
resembles the "random orthogonal directions" rule but potentially offers more
relevant directions than random ones.

Algorithm 3 Generic iterative type I method

First-order oracle $f$, initial iterate and smoothness $x_{0},\,M_{0}$, number
of iterations $T$.

for $t=0,\,\ldots,\,T-1$ do

Update $Y_{t}$, $Z_{t}$, $D_{t}$, $G_{t}$, and $\varepsilon_{t}$ (see section
2.4).

$x_{t+1},M_{t+1}\leftarrow\texttt{[Alg.
\ref{alg:type1}]}(f,G_{t},D_{t},\varepsilon_{t},x_{t},\frac{M_{t}}{2})$

end for

return $x_{T}$

Algorithm 4 Type I subroutine with backtracking line-search

1:First-order oracle for $f$, matrices $G,\,D$, vector $\varepsilon$, iterate
$x$, initial smoothness $M_{0}$.

2:Initialize $M\leftarrow\frac{M_{0}}{2}$

3:do

4: $M\leftarrow 2M\;\;$ and
$\;\;H\leftarrow\frac{G^{T}D+D^{T}G}{2}+\mathrm{I}_{N}\frac{M\|D\|\|\varepsilon\|}{2}$

5: $\alpha^{\star}\leftarrow\min\limits_{\alpha}f(x)+\nabla
f(x)^{T}D\alpha+\frac{\alpha^{T}H\alpha}{2}+\frac{M\|D\alpha\|^{3}}{6}$

6: $x_{+}\leftarrow x+D\alpha$

7:while $\textstyle f(x_{+})\geq f(x)+\nabla
f(x)^{T}D\alpha^{\star}+\frac{[\alpha^{\star}]^{T}H\alpha^{\star}}{2}+\frac{M\|D\alpha^{\star}\|^{3}}{6}$

8:return $x_{+}$, $M$

Algorithm 5 Type II subroutine with backtracking line-search

Same as algorithm 4, but minimize and check the upper bound eq. Type II bound
instead of eq. Type I bound on lines 4 and 6.

###  2.5 Miscellaneous

##### Link with qN and Anderson acceleration

The algorithms 4 and 5 are strongly related to known quasi-Newton methods and
Anderson acceleration technique, see section D.3.

##### Backtracking Line-Search

The smoothness parameter is replaced by $M_{t}$, found by backtracking [61].
The parameter $M_{0}$ is estimated by finite-difference, see section F.1.

##  3 Convergence Rates for Iterative Type I Methods

This section analyzes the convergence rates of the methods that use algorithm
4 as a subroutine; see algorithm 3. An analysis of methods that use algorithm
5 is left for future work.

###  3.1 Assumptions

This section lists the important assumptions regarding the function $f$.
Subsequent results require a bound on the radius of the sublevel set
$\\{x:f(x)\leq f(x_{0})\\}$.

######  Assumption 2.

The radius of the sub-level set $\\{x:f(x)\leq f(x_{0})\\}$ is bounded by
$R<\infty$.

Some results require $f$ to be star convex or convex to ensure convergence
toward $f(x^{\star})$.

######  Assumption 3.

The function $f$ is star convex if, for all $x\in\mathbb{R}^{d}$ and
$\forall\tau\in[0,1]$,

$f((1-\tau)x+\tau x^{\star})\leq(1-\tau)f(x)+\tau f(x^{\star}).$

######  Assumption 4.

The function $f$ is convex if, for all $y,\,z\in\mathbb{R}^{d}$, $f(y)\geq
f(z)+\nabla f(z)(y-z)$.

###  3.2 Rates of Convergence

When $f$ satisfies Assumption 1, algorithm 3 ensures a consistent minimal
decrease in function at every step.

######  Theorem 4.

Let $f$ satisfy Assumption 1. Subsequently, at each iteration $t\geq 0$,
algorithm 3 starting at $x_{0}$ with $M_{0}>0$ achieves

$\displaystyle\textstyle f(x_{t+1})\leq
f(x_{t})-\frac{M_{t+1}}{12}\|x_{t+1}-x_{t}\|^{3},$ (11) with $\displaystyle
M_{t+1}<\max\left\\{2L\;;\;\frac{M_{0}}{2^{t}}\right\\}.$

Moreover, the total function evaluation is bounded by
$2t+\log_{2}\left(\frac{M_{0}}{L}\right)$.

Under some mild assumptions, algorithm 3 converges to a critical point for
nonconvex functions and converges to an optimum when the function is star-
convex.

######  Theorem 5.

Let $f$ satisfy Assumption 1 and assume that $f$ is bounded below by
$f^{\star}$. Let Requirements 1b, 2 and 3 hold, and $M_{t}\geq M_{\min}$.
Subsequently, algorithm 3 starting at $x_{0}$ with $M_{0}>0$ achieves

$\displaystyle\min_{i=1,\,\ldots,\,t}\|\nabla f(x_{i})\|\leq$
$\displaystyle\max\Big{\\{}\frac{3L}{t^{2/3}}\left(12\frac{f(x_{0})-f^{\star}}{M_{\min}}\right)^{2/3}\;;$
$\displaystyle\left(\frac{C_{1}}{t^{1/3}}\right)\left(12\frac{f(x_{0})-f^{\star}}{M_{\min}}\right)^{1/3}\Big{\\}},$

$\textstyle C_{1}=\delta
L\left(\frac{\kappa+2\kappa^{2}}{2}\right)+\max_{i\in[0,t]}\|(I-P_{i})\nabla^{2}f(x_{i})P_{i}\|.$

######  Theorem 6.

Let $f$ satisfy Assumptions 1, 2 and 3. Let Requirements 1b, 2 and 3 hold.
Then, for $t\geq 1$, algorithm 3 starting at $x_{0}$ with $M_{0}>0$ achieves

$\displaystyle f(x_{t})-f^{\star}$ $\displaystyle\leq
6\frac{f(x_{0})-f^{\star}}{t(t+1)(t+2)}$
$\displaystyle+\frac{1}{(t+1)(t+2)}\frac{L(3R)^{3}}{2}+\frac{1}{t+2}\frac{C_{2}(3R)^{2}}{4},$

$C_{2}\stackrel{{\scriptstyle\text{def}}}{{=}}\delta
L\frac{\kappa+2\kappa^{2}}{2}+\max_{i\in[0,t]}\|\nabla^{2}f(x_{i})-P_{i}\nabla^{2}f(x_{i})P_{i}\|.$

The next theorem shows that when random directions (that satisfy Requirement
1a) are used, $f(x_{t})$ also converges in expectation to $f(x^{\star})$ when
$f$ is convex.

######  Theorem 7.

Assume $f$ satisfy Assumptions 1, 2 and 4. Let Requirements 1a, 2 and 3 hold.
Then in expectation over the matrices $D_{i}$ and for $t\geq 1$, algorithm 3
starting at $x_{0}$ with $M_{0}>0$ achieves

$\displaystyle\mathbb{E}_{D_{t}}[f(x_{t})-f^{\star}]\leq$
$\displaystyle\frac{1}{1+\frac{1}{4}\left[\frac{N}{d}t\right]^{3}}(f(x_{0})-f^{\star})$
$\displaystyle+\frac{1}{\left[\frac{N}{d}t\right]^{2}}\frac{L(3R)^{3}}{2}+\frac{1}{\left[\frac{N}{d}t\right]}\frac{C_{3}(3R)^{2}}{2},$

$C_{3}\stackrel{{\scriptstyle\text{def}}}{{=}}\textstyle\delta
L\frac{\kappa+2\kappa^{2}}{2}+\frac{(d-N)}{d}\max_{i\in[0,t]}\|\nabla^{2}f(x_{i})\|$.

##### Accelerated method

Due to space limitation, the accelerated algorithm 3 is presented in appendix
A, see algorithms 6 and 7. Indeed, while the algorithm is theoretically
intriguing, it underperforms numerically, likely because it trades-off its
adaptivity for better worst-case convergence rates.

######  Theorem 8.

Assume $f$ satisfies Assumptions 1, 2 and 4. Let Requirements 1b, 2 and 3
hold. Subsequently, for $t\geq 1$, the accelerated algorithm 7 starting at
$x_{0}$ with $M_{0}>0$ achieves

$\displaystyle f(x_{t})-f^{\star}\leq$ $\displaystyle
C_{4}\frac{(3R)^{2}}{(t+3)^{2}}+9\max\left\\{M_{0}\;;\;2L\right\\}\left(\frac{3R}{t+3}\right)^{3}$
$\displaystyle+\frac{\frac{\tilde{\lambda}^{(1)}R^{2}}{2}+\frac{\tilde{\lambda}^{(2)}R^{3}}{6}}{(t+1)^{3}}.$
$\displaystyle\textstyle\tilde{\lambda}^{(1)}=\frac{\delta}{2}\left(L\kappa+M_{1}\kappa^{2}\right)+\|\nabla^{2}f(x_{0})-P_{0}\nabla^{2}f(x_{0})P_{0}\|,$
$\displaystyle\tilde{\lambda}^{(2)}=M_{1}+L,$ $\displaystyle
C_{4}=30\kappa_{D}\left(\delta\max\\{4L,M_{0}\\}+\max_{i\leq
t}\|(I-P_{i})\nabla f(x_{i})P_{i})\|\right)$

###  3.3 Interpretation, Comparison With First-Order Methods and Special Cases

The rates presented in Theorems 5, 6, 7 and 8 combine the rates of the cubic
regularized Newton’s method and gradient descent (or coordinate descent, as in
Theorem 7) for functions with Lipschitz-continuous Hessian. As $C_{1},C_{2}$,
$C_{3}$, and $C_{4}$ decrease, the rates approach those of second-order
methods. For simplicity in this section, those constants are denoted as $C$,
and $\kappa$ is assumed to be equal to $1$.

####  3.3.1 Interpretation and Comparison

The constant $C$ quantifies the estimation error of
$D_{t}^{T}\nabla^{2}f(x_{t})D_{t}$ by $H_{t}$ in (Type I bound) into two
terms:

$\textstyle C\leq O\big{(}\delta L+\textstyle\max_{i\leq
t}\|(I-P_{i})\nabla^{2}f(x_{i})\|\Big{)}.$

The first term is the error caused by approximating $\nabla^{2}f(x)D_{t}$ by
$G_{t}$ and the second is the subspace approximation error of
$\nabla^{2}f(x_{t})$ in the span of the columns of $D_{t}$. This approximation
is more explicit in $C_{3}$, where increasing $N$ reduces the constant to
$N=d$.

The rate associated with the constant $C$ is a perturbed version of the rate
of first order methods (appendix C): if the function has Lipchitz-continuous
gradients with constant $\mathcal{L}$ (appendix C, eq. 12), in the worst case
over $P_{i}$, the constant $C$ is bounded as $C\leq O(\delta L+\mathcal{L})$.
Hence, the rates of Theorems 5, 6 and 8 are perturbed versions of the rates of
gradient descent and accelerated gradient descent, respectively, whereas
Theorem 7 is a perturbed version of the rate of coordinate descent. The
perturbation is not surprising, as this study did not assume that $f$ has a
Lipchitz continuous gradient.

####  3.3.2 Special Cases

The framework presented in this study can be reduced to some known methods,
such as the cubic regularization of Newton’s method. This section explores
some special/extreme cases of this study’s framework when the update rule
reads $Y_{t}=Z_{t}+hD_{t}$ for some $D_{t}$ that satisfy the requirements
listed in section 2.4. Consider the following three extreme cases.

  1. a)

The stepsize $h$ tends to 0, hence $y_{i}\rightarrow z_{i}$ and

$G_{t}\rightarrow[\ldots,\nabla^{2}f(z_{i}^{(t)})d_{i}^{(t)},\ldots]_{i=1\ldots
N}.$

This requires a second-order oracle, but makes $\delta$ smaller, that is,
$\delta=O(\|x_{t}-z_{i}^{(t)}\|)$.

  2. b)

The forward estimates are all centered in $x_{t}$, hence
$Z_{t}=[x_{t},\ldots,x_{t}]$. Then, at each iteration $t$, computing $G_{t}$
requires $N$ additional gradients, but makes $\delta$ smaller, that is,
$\delta=O(h)$.

  3. c)

The algorithm’s memory equals $d$ (full memory version), hence $P_{t}=I_{d}$
for all $t$. This raises the complexity of each iteration to $O(d^{3})$, but
this makes $(I-P_{t})\nabla f(x_{t})=0$.

The combinations of these special cases are summarized in table 2.

| Satisfies a) | Satisfies b) | Satisfies c)
---|---|---|---
Limited Memory QN with guarantees (this paper) |  |  |
Quasi-Newton with Second Order Oracle [67, 66, 65] | ✓ |  |
Subspace Newton with finite-difference (this paper) |  | ✓ |
Full-Memory QN with guarantees (this paper, [50]) |  |  | ✓
Random Subspace [26, 39, 44] | ✓ | ✓ |
Cubic Newton with finite-difference [40] |  | ✓ | ✓
Cubic Newton with Lazy Hessian [25] | ✓ |  | ✓
Cubic Newton [61, 59] | ✓ | ✓ | ✓
Table 2: By considering the extreme cases a), b), and c) from section 3.3.2,
the algorithm presented in this study reduces to known methods, up to the
Cubic regularization of Newton’s method in the most extreme setting.

##  4 Numerical Experiments

This section compares the methods generated by this study’s framework to the
l-BFGS algorithm from minFunc [69], see fig. 1. Additional experiments are
described in appendix F. The tested methods were type-I iterative algorithms
(algorithm 3 using the techniques from section 2.4). The step size of the
forward estimation was set to $h=10^{-9}$, and the condition number
$\kappa_{D_{t}}$ was maintained below $\kappa=10^{9}$ using iterate-only and
greedy techniques. The accelerated algorithm 7 is used only with the  forward
estimate-only technique. These methods are evaluated on a logistic regression
problem on the Madelon UCI dataset [42].

Regarding the number of iterations, the greedy orthogonalized version
outperformed the others owing to the orthogonality of the directions
(resulting in a condition number of one) and the meaningfulness of the
previous gradients/iterates. However, in terms of gradient oracle calls, the
recommended method, orthogonal forward estimate only, achieved the best
performance by balancing the cost per iteration (only two gradients per
iteration) and efficiency (small and orthogonal directions, reducing
theoretical constants). Surprisingly, the accelerated method underperforms,
likely owing to its tightened theoretical analysis reducing its inherent
adaptivity.

##  5 Conclusion and Future work

This study introduces a generic framework for developing novel quasi-Newton
and Anderson/nonlinear acceleration schemes that offer a global convergence
rate in various scenarios, including accelerated convergence on convex
functions, with minimal assumptions.

The proposed approach requires an additional gradient step for the forward
estimate, as discussed in Section 2.4. However, this forward estimate is
crucial for enabling the algorithm’s adaptivity.

Studying the special case $N=d$ (although unsuitable for large-scale problems)
could highlight super-linear convergence rates in future research. Moreover,
using the average-case analysis framework from existing literature [63, 74,
27, 21, 62] can also improve the constants in Theorems 6 and 5 to match those
in Theorem 7. Furthermore, exploring the convergence rates of the type-2
method is also worthwhile.

Ultimately, the results presented in this study open new avenues of research.
It may also provide a potential foundation for investigating additional
properties of existing quasi-Newton methods. This may even lead to the
discovery of convergence rates for adaptive cubic-regularized BFGS variants.

##### Main limitation

Despite its strong theoretical convergence rates and its ability to recover
cubic regularization in the limit, the algorithm’s complexity rises to
$O(d^{3})$ when $N$ approaches $d$ owing to the cubic minimization subproblem.
This contrasts with current full-memory qN methods, which have a complexity of
$O(d^{2})$ but offer local superlinear convergence rates [65, 66, 67]. Solving
the cubic regularized model is an active research topic [16, 15, 51, 33, 48],
and its study is out of the scope of this paper. Nevertheless, it is very
likely that future researchs may reduce the complexity to $O(d^{2})$ by
updating the solution of the subproblem with low-rank perturbations, or by
solving inexactly the cubic subproblem like in [50].

Figure 1: Comparison between the type-1 methods proposed in this study and the
optimized implementation of l-BFGS from minFunc [69] with default parameters,
except for the memory size. All methods use a memory size of $N=25$.

## References

  * [1] Hari Om Aggrawal and Jan Modersitzki  “Hessian Initialization Strategies for $\ell$-BFGS Solving Non-linear Inverse Problems”  In _International Conference on Scale Space and Variational Methods in Computer Vision_ , 2021, pp. 216–228  Springer
  * [2] Donald G Anderson  “Iterative procedures for nonlinear integral equations”  In _Journal of the ACM (JACM)_ 12.4, 1965, pp. 547–560
  * [3] Kimon Antonakopoulos, Ali Kavis and Volkan Cevher  “Extra-Newton: A First Approach to Noise-Adaptive Accelerated Second-Order Methods”  In _arXiv preprint arXiv:2211.01832_ , 2022
  * [4] Claude Brezinski  “Application de l’$\varepsilon$-algorithme à la résolution des systèmes non linéaires”  In _Comptes Rendus de l’Académie des Sciences de Paris_ 271.A, 1970, pp. 1174–1177
  * [5] Claude Brezinski  “Sur un algorithme de résolution des systèmes non linéaires”  In _Comptes Rendus de l’Académie des Sciences de Paris_ 272.A, 1971, pp. 145–148
  * [6] Claude Brezinski and Michela Redivo–Zaglia  “The genesis and early developments of Aitken’s process, Shanks’ transformation, the $\varepsilon$–algorithm, and related fixed point methods”  In _Numerical Algorithms_ 80.1 Springer, 2019, pp. 11–133
  * [7] Claude Brezinski, Michela Redivo-Zaglia and Yousef Saad  “Shanks sequence transformations and Anderson acceleration”  In _SIAM Review_ 60.3 SIAM, 2018, pp. 646–669
  * [8] Claude Brezinski and M Redivo Zaglia  “Extrapolation methods: theory and practice”  Elsevier, 1991
  * [9] Claude Brezinski, Stefano Cipolla, Michela Redivo-Zaglia and Yousef Saad  “Shanks and Anderson-type acceleration techniques for systems of nonlinear equations”  In _arXiv:2007.05716_ , 2020
  * [10] Charles G Broyden  “The convergence of a class of double-rank minimization algorithms: 2. The new algorithm”  In _IMA journal of applied mathematics_ 6.3 Oxford University Press, 1970, pp. 222–231
  * [11] Charles George Broyden  “The convergence of a class of double-rank minimization algorithms 1. general considerations”  In _IMA Journal of Applied Mathematics_ 6.1 Oxford University Press, 1970, pp. 76–90
  * [12] Richard H Byrd and Jorge Nocedal  “A tool for the analysis of quasi-Newton methods with application to unconstrained minimization”  In _SIAM Journal on Numerical Analysis_ 26.3 SIAM, 1989, pp. 727–739
  * [13] Richard H Byrd, Jorge Nocedal and Ya-Xiang Yuan  “Global convergence of a class of quasi-Newton methods on convex problems”  In _SIAM Journal on Numerical Analysis_ 24.5 SIAM, 1987, pp. 1171–1190
  * [14] Marco Canini and Peter Richtárik  “Direct nonlinear acceleration”  In _Operational Research_ 2192, 2022, pp. 4406
  * [15] Yair Carmon and John Duchi  “Gradient descent finds the cubic-regularized nonconvex Newton step”  In _SIAM Journal on Optimization_ 29.3 SIAM, 2019, pp. 2146–2178
  * [16] Yair Carmon and John C Duchi  “Analysis of krylov subspace solutions of regularized non-convex quadratic problems”  In _Advances in Neural Information Processing Systems_ 31, 2018
  * [17] Yair Carmon, Arun Jambulapati, Yujia Jin and Aaron Sidford  “Recapp: Crafting a more efficient catalyst for convex optimization”  In _International Conference on Machine Learning_ , 2022, pp. 2658–2685  PMLR
  * [18] Coralia Cartis, Nicholas IM Gould and Philippe L Toint  “Adaptive cubic regularisation methods for unconstrained optimization. Part I: motivation, convergence and numerical results”  In _Mathematical Programming_ 127.2 Springer, 2011, pp. 245–295
  * [19] Coralia Cartis, Nicholas IM Gould and Philippe L Toint  “Adaptive cubic regularisation methods for unconstrained optimization. Part II: worst-case function-and derivative-evaluation complexity”  In _Mathematical programming_ 130.2 Springer, 2011, pp. 295–319
  * [20] Andrew R Conn, Nicholas IM Gould and Ph L Toint  “Convergence of quasi-Newton matrices generated by the symmetric rank one update”  In _Mathematical programming_ 50.1-3 Springer, 1991, pp. 177–195
  * [21] Leonardo Cunha et al.  “Only tails matter: Average-Case Universality and Robustness in the Convex Regime”, 2022
  * [22] Alexandre d’Aspremont, Damien Scieur and Adrien Taylor  “Acceleration methods”  In _Foundations and Trends® in Optimization_ 5.1-2 Now Publishers, Inc., 2021, pp. 1–245
  * [23] William C Davidon  “Variable metric method for minimization”  In _SIAM Journal on Optimization_ 1.1 SIAM, 1991, pp. 1–17
  * [24] Nikita Doikov, El Mahdi Chayti and Martin Jaggi  “Second-order optimization with lazy Hessians”  In _arXiv preprint arXiv:2212.00781_ , 2022
  * [25] Nikita Doikov and Martin Jaggi  “Second-order optimization with lazy Hessians”  In _International Conference on Machine Learning_ , 2023, pp. 8138–8161  PMLR
  * [26] Nikita Doikov and Peter Richtárik  “Randomized block cubic Newton method”  In _International Conference on Machine Learning_ , 2018, pp. 1290–1298  PMLR
  * [27] Carles Domingo-Enrich, Fabian Pedregosa and Damien Scieur  “Average-case acceleration for bilinear games and normal matrices”  In _arXiv preprint arXiv:2010.02076_ , 2020
  * [28] John R Engels and Hector J Martinez  “Local and superlinear convergence for partially known quasi-Newton methods”  In _SIAM Journal on Optimization_ 1.1 SIAM, 1991, pp. 42–56
  * [29] Haw-Ren Fang and Yousef Saad  “Two classes of multisecant methods for nonlinear acceleration”  In _Numerical Linear Algebra with Applications_ 16.3, 2009, pp. 197–221
  * [30] Roger Fletcher  “A new approach to variable metric algorithms”  In _The computer journal_ 13.3 Oxford University Press, 1970, pp. 317–322
  * [31] Roger Fletcher and Michael JD Powell  “A rapidly convergent descent method for minimization”  In _The computer journal_ 6.2 Oxford University Press, 1963, pp. 163–168
  * [32] William F Ford and Avram Sidi  “Recursive algorithms for vector extrapolation methods”  In _Applied numerical mathematics_ 4.6 Elsevier, 1988, pp. 477–489
  * [33] Yihang Gao, Man-Chung Yue and Michael Ng  “Approximate Secular Equations for the Cubic Regularization Subproblem”  In _Advances in Neural Information Processing Systems_ 35, 2022, pp. 14250–14260
  * [34] Alexander Gasnikov et al.  “Near optimal methods for minimizing convex functions with lipschitz $p$-th derivatives”  In _Conference on Learning Theory_ , 2019, pp. 1392–1393  PMLR
  * [35] Eckart Gekeler  “On the solution of systems of equations by the epsilon algorithm of Wynn”  In _Mathematics of Computation_ 26.118, 1972, pp. 427–436
  * [36] Saeed Ghadimi, Han Liu and Tong Zhang  “Second-order methods with cubic regularization under inexact information”  In _arXiv preprint arXiv:1710.05782_ , 2017
  * [37] Hiva Ghanbari and Katya Scheinberg  “Proximal quasi-Newton methods for regularized convex optimization with linear and accelerated sublinear convergence rates”  In _Computational Optimization and Applications_ 69 Springer, 2018, pp. 597–627
  * [38] Donald Goldfarb  “A family of variable-metric methods derived by variational means”  In _Mathematics of computation_ 24.109, 1970, pp. 23–26
  * [39] Robert Gower, Dmitry Kovalev, Felix Lieder and Peter Richtárik  “Rsn: Randomized subspace newton”  In _Advances in Neural Information Processing Systems_ 32, 2019
  * [40] Geovani Nunes Grapiglia, Max LN Gonçalves and GN Silva  “A cubic regularization of Newton’s method with finite difference Hessian approximations”  In _Numerical Algorithms_ Springer, 2022, pp. 1–24
  * [41] Andreas Griewank and Ph L Toint  “Local convergence analysis for partitioned quasi-Newton updates”  In _Numerische Mathematik_ 39.3 Springer, 1982, pp. 429–448
  * [42] Isabelle Guyon  “Design of experiments of the NIPS 2003 variable selection benchmark”  In _NIPS 2003 workshop on feature extraction and feature selection_ 253, 2003
  * [43] Isabelle Guyon et al.  “Design and analysis of the causation and prediction challenge”  In _Causation and Prediction Challenge_ , 2008, pp. 1–33  PMLR
  * [44] Filip Hanzely, Nikita Doikov, Yurii Nesterov and Peter Richtarik  “Stochastic subspace cubic Newton method”  In _International Conference on Machine Learning_ , 2020, pp. 4027–4038  PMLR
  * [45] K Jbilou and H Sadok  “Vector extrapolation methods. Applications and numerical comparison”  In _Journal of Computational and Applied Mathematics_ 122.1-2 Elsevier, 2000, pp. 149–165
  * [46] Khalide Jbilou and Hassane Sadok  “Analysis of some vector extrapolation methods for solving systems of linear equations”  In _Numerische Mathematik_ 70.1 Springer, 1995, pp. 73–89
  * [47] Khalide Jbilou and Hassane Sadok  “Some results about vector extrapolation methods and related fixed-point iterations”  In _Journal of Computational and Applied Mathematics_ 36.3 Elsevier, 1991, pp. 385–398
  * [48] Ru-Jun Jiang, Zhi-Shuo Zhou and Zi-Rui Zhou  “Cubic regularization methods with second-order complexity guarantee based on a new subproblem reformulation”  In _Journal of the Operations Research Society of China_ 10.3 Springer, 2022, pp. 471–506
  * [49] Ruichen Jiang, Qiujiang Jin and Aryan Mokhtari  “Online Learning Guided Curvature Approximation: A Quasi-Newton Method with Global Non-Asymptotic Superlinear Convergence”  In _arXiv preprint arXiv:2302.08580_ , 2023
  * [50] Ruichen Jiang and Aryan Mokhtari  “Accelerated Quasi-Newton Proximal Extragradient: Faster Rate for Smooth Convex Optimization”  In _arXiv preprint arXiv:2306.02212_ , 2023
  * [51] Rujun Jiang, Man-Chung Yue and Zhishuo Zhou  “An accelerated first-order method with complexity analysis for solving cubic regularization subproblems”  In _Computational Optimization and Applications_ 79 Springer, 2021, pp. 471–506
  * [52] Dmitry Kamzolov, Klea Ziu, Artem Agafonov and Martin Takáč  “Accelerated Adaptive Cubic Regularized Quasi-Newton Methods”  In _arXiv preprint arXiv:2302.04987_ , 2023
  * [53] Dmitry Kovalev and Alexander Gasnikov  “The first optimal acceleration of high-order methods in smooth convex optimization”  In _arXiv preprint arXiv:2205.09647_ , 2022
  * [54] Dachao Lin, Haishan Ye and Zhihua Zhang  “Explicit convergence rates of greedy and random quasi-Newton methods”  In _Journal of Machine Learning Research_ 23.162, 2022, pp. 1–40
  * [55] Dachao Lin, Haishan Ye and Zhihua Zhang  “Greedy and random quasi-newton methods with faster explicit superlinear convergence”  In _Advances in Neural Information Processing Systems_ 34, 2021, pp. 6646–6657
  * [56] Dong C Liu and Jorge Nocedal  “On the limited memory BFGS method for large scale optimization”  In _Mathematical programming_ 45.1-3 Springer, 1989, pp. 503–528
  * [57] Renato DC Monteiro and Benar Fux Svaiter  “An accelerated hybrid proximal extragradient method for convex optimization and its implications to second-order methods”  In _SIAM Journal on Optimization_ 23.2 SIAM, 2013, pp. 1092–1125
  * [58] José Luis Morales  “A numerical study of limited memory BFGS methods”  In _Applied Mathematics Letters_ 15.4 Elsevier, 2002, pp. 481–487
  * [59] Yurii Nesterov  “Accelerating the cubic regularization of Newton’s method on convex problems”  In _Mathematical Programming_ 112.1 Springer, 2008, pp. 159–181
  * [60] Yurii Nesterov  “Introductory lectures on convex optimization”  Springer, 2004
  * [61] Yurii Nesterov and Boris T Polyak  “Cubic regularization of Newton method and its global performance”  In _Mathematical Programming_ 108.1 Springer, 2006, pp. 177–205
  * [62] Courtney Paquette, Bart Merriënboer, Elliot Paquette and Fabian Pedregosa  “Halting Time is predictable for large models: A universality property and average-case analysis”  In _Foundations of Computational Mathematics_ Springer, 2022
  * [63] Fabian Pedregosa and Damien Scieur  “Acceleration through spectral density estimation”  In _Proceedings of the 37th International Conference on Machine Learning (ICML)_ , 2020
  * [64] MJD Powell  “How bad are the BFGS and DFP methods when the objective function is quadratic?”  In _Mathematical Programming_ 34 Springer, 1986, pp. 34–47
  * [65] Anton Rodomanov and Yurii Nesterov  “Greedy quasi-Newton methods with explicit superlinear convergence”  In _SIAM Journal on Optimization_ 31.1 SIAM, 2021, pp. 785–811
  * [66] Anton Rodomanov and Yurii Nesterov  “New results on superlinear convergence of classical quasi-Newton methods”  In _Journal of optimization theory and applications_ 188 Springer, 2021, pp. 744–769
  * [67] Anton Rodomanov and Yurii Nesterov  “Rates of superlinear convergence for classical quasi-Newton methods”  In _Mathematical Programming_ Springer, 2021, pp. 1–32
  * [68] Katya Scheinberg and Xiaocheng Tang  “Practical inexact proximal quasi-Newton method with global complexity analysis”  In _Mathematical Programming_ 160 Springer, 2016, pp. 495–529
  * [69] Mark Schmidt  “minFunc: unconstrained differentiable multivariate optimization in Matlab”  In _Software available at http://www. cs. ubc. ca/~ schmidtm/Software/minFunc. htm_ , 2005
  * [70] Robert B Schnabel  “Quasi-Newton Methods Using Multiple Secant Equations.”, 1983
  * [71] Damien Scieur  “Generalized framework for nonlinear acceleration”  In _arXiv preprint arXiv:1903.08764_ , 2019
  * [72] Damien Scieur, Alexandre d’Aspremont and Francis Bach  “Regularized nonlinear acceleration”  In _Advances in Neural Information Processing Systems (NIPS)_ , 2016
  * [73] Damien Scieur, Alexandre d’Aspremont and Francis Bach  “Regularized nonlinear acceleration”  In _Mathematical Programming_ Springer, 2020
  * [74] Damien Scieur and Fabian Pedregosa  “Universal Asymptotic Optimality of Polyak Momentum”  In _Proceedings of the 37th International Conference on Machine Learning (ICML)_ , 2020
  * [75] Damien Scieur, Lewis Liu, Thomas Pumir and Nicolas Boumal  “Generalization of Quasi-Newton methods: application to robust symmetric multisecant updates”  In _International Conference on Artificial Intelligence and Statistics_ , 2021, pp. 550–558  PMLR
  * [76] Damien Scieur, Edouard Oyallon, Alexandre d’Aspremont and Francis Bach  “Online Regularized Nonlinear Acceleration”  In _arXiv:1805.09639_ , 2018
  * [77] David F Shanno  “Conditioning of quasi-Newton methods for function minimization”  In _Mathematics of computation_ 24.111, 1970, pp. 647–656
  * [78] Avram Sidi  “Convergence and stability properties of minimal polynomial and reduced rank extrapolation algorithms”  In _SIAM Journal on Numerical Analysis_ 23.1 SIAM, 1986, pp. 197–209
  * [79] Avram Sidi  “Efficient implementation of minimal polynomial and reduced rank extrapolation methods”  In _Journal of Computational and Applied Mathematics_ 36.3 Elsevier, 1991, pp. 305–337
  * [80] Avram Sidi  “Extrapolation vs. projection methods for linear systems of equations”  In _Journal of Computational and Applied Mathematics_ 22.1 Elsevier, 1988, pp. 71–88
  * [81] Avram Sidi  “Minimal polynomial and reduced rank extrapolation methods are related”  In _Advances in Computational Mathematics_ 43.1 Springer, 2017, pp. 151–170
  * [82] Avram Sidi  “Vector extrapolation methods with applications”  SIAM, 2017
  * [83] Avram Sidi  “Vector extrapolation methods with applications to solution of large systems of equations and to PageRank computations”  In _Computers & Mathematics with Applications_ 56.1 Elsevier, 2008, pp. 1–24
  * [84] Avram Sidi and Jacob Bridger  “Convergence and stability analyses for some vector extrapolation methods in the presence of defective iteration matrices”  In _Journal of Computational and Applied Mathematics_ 22.1 Elsevier, 1988, pp. 35–61
  * [85] Avram Sidi and Yair Shapira  “Upper bounds for convergence rates of acceleration methods with initial iterations”  In _Numerical Algorithms_ 18.2 Springer, 1998, pp. 113–132
  * [86] Andrzej Stachurski  “Superlinear convergence of Broyden’s bounded $\theta$-class of methods”  In _Mathematical Programming_ 20.1 Springer, 1981, pp. 196–212
  * [87] Alex Toth and CT Kelley  “Convergence analysis for Anderson acceleration”  In _SIAM Journal on Numerical Analysis_ 53.2 SIAM, 2015, pp. 805–819
  * [88] Evgenij E Tyrtyshnikov  “How bad are Hankel matrices?”  In _Numerische Mathematik_ 67.2 Springer, 1994, pp. 261–269
  * [89] Gerhard Venter  “Review of optimization techniques”  John Wiley & Sons Ltd., 2010
  * [90] Homer F Walker and Peng Ni  “Anderson acceleration for fixed-point iterations”  In _SIAM Journal on Numerical Analysis_ 49.4 SIAM, 2011, pp. 1715–1735
  * [91] Zhe Wang, Yi Zhou, Yingbin Liang and Guanghui Lan  “A Note on Inexact Condition for Cubic Regularized Newton’s Method”  In _arXiv preprint arXiv:1808.07384_ , 2018
  * [92] Zengxin Wei, Gaohang Yu, Gonglin Yuan and Zhigang Lian  “The superlinear convergence of a modified BFGS-type method for unconstrained optimization”  In _Computational optimization and applications_ 29 Springer, 2004, pp. 315–332
  * [93] Hiroshi Yabe, Hideho Ogasawara and Masayuki Yoshino  “Local and superlinear convergence of quasi-Newton methods based on modified secant conditions”  In _Journal of Computational and Applied Mathematics_ 205.1 Elsevier, 2007, pp. 617–632
  * [94] Hiroshi Yabe and Naokazu Yamaki  “Local and superlinear convergence of structured quasi-Newton methods for nonlinear optimization”  In _Journal of the Operations Research Society of Japan_ 39.4 The Operations Research Society of Japan, 1996, pp. 541–557
  * [95] Junzi Zhang, Brendan O’Donoghue and Stephen Boyd  “Globally convergent type-I Anderson acceleration for nonsmooth fixed-point iterations”  In _SIAM Journal on Optimization_ 30.4 SIAM, 2020, pp. 3170–3197

###### Contents

  1. 1 Introduction
    1. 1.1 Contributions
    2. 1.2 Limitations in Current qN Schemes
  2. 2 Rethinking Quasi-Newton Methods
    1. 2.1 First Ingredient: Subspace Minimization
    2. 2.2 Second Ingredient: Multisecant Approximation of the Hessian
    3. 2.3 Third Ingredient: Objective Function and Gradient Norm Upper bounds
    4. 2.4 Fourth Ingredient: Direction Update Rule
      1. 2.4.1 "Orthogonal Forward Estimate Only" Update Rule (Recommended)
      2. 2.4.2 "Random Orthogonal Directions"
      3. 2.4.3 Other Matrix Updates: Pruning or Orthogonalization
    5. 2.5 Miscellaneous
  3. 3 Convergence Rates for Iterative Type I Methods
    1. 3.1 Assumptions
    2. 3.2 Rates of Convergence
    3. 3.3 Interpretation, Comparison With First-Order Methods and Special Cases
      1. 3.3.1 Interpretation and Comparison
      2. 3.3.2 Special Cases
  4. 4 Numerical Experiments
  5. 5 Conclusion and Future work
  6. A Accelerated Algorithm
  7. B Related work
    1. B.1 Inexact, Subspace, and Stochastic Newton Methods
    2. B.2 Nonlinear and Anderson Acceleration
    3. B.3 Quasi-Newton Methods
    4. B.4 Close Related Work
      1. B.4.1 (Accelerated) Quasi-Newton with Secant Inexactness
      2. B.4.2 ARC: Adaptive Regularization algorithm using Cubics
      3. B.4.3 Proximal Quasi-Newton Methods
      4. B.4.4 Proximal Extragradient Quasi-Newton Methods with Online Estimation
  8. C Known rates of convergence and Comparison
    1. C.1 (Accelerated) Gradient Descent
    2. C.2 (Accelerated) Cubic Regularized Newton’s Method
    3. C.3 Relation Between Parameters
    4. C.4 Comparing rates of convergence
  9. D Link with quasi-Newton and Anderson/Nonlinear Acceleration
    1. D.1 Anderson Acceleration and Nonlinear Acceleration
    2. D.2 Single-secant and Multisecant Quasi-Newton Methods
      1. D.2.1 The Ideas Behind Single-Secant and Multisecant Hessian Approximation
      2. D.2.2 Davidon-Fletcher-Powell (DFP) Formula
      3. D.2.3 Multisecant Broyden Methods
      4. D.2.4 Link with Anderson Acceleration
    3. D.3 Links with Algorithms 4 and 5
  10. E Solving the sub-problems
    1. E.1 Solving the Type 1 Subproblem
    2. E.2 Solving the Type 2 Subproblem
      1. E.2.1 Fundamentals of SOCP
      2. E.2.2 SOCP Formulation of the Type 2 Subproblem
  11. F Additional Numerical Experiments
    1. F.1 Initial Parameter for the Backtracking Line search
    2. F.2 Scalability w.r.t. Dimension and Memory
    3. F.3 Influence of h
    4. F.4 Impact of the memory parameter N
    5. F.5 Nonconvex optimization
    6. F.6 Comparison of Type 1 Methods on Convex Problems
      1. F.6.1 Square loss and cubic regularization
      2. F.6.2 Logistic regression
    7. F.7 Comparison of Type 2 Methods on Convex Problems
      1. F.7.1 Square loss and cubic regularization
      2. F.7.2 Logistic regression
  12. G Missing proofs
    1. G.1 Technical Result: Hessian Approximation
    2. G.2 Technical Results: Cubic Subproblem
    3. G.3 Technical Results: Decrease Guarantees
    4. G.4 Technical Results: Accelerated Algorithm
    5. G.5 Missing proofs from Sections 2 and 3

## Supplementary Materials

##  Appendix A Accelerated Algorithm

This section introduces algorithm 7, an accelerated variant of algorithm 3 for
convex functions, designed using the estimate sequence technique from [59]. It
consists in iteratively building a function $\Phi_{t}(x)$, that reads

$\textstyle\Phi_{t}(x)=\frac{1}{\sum_{i=0}^{t}b_{i}}\left(\sum_{i=0}^{t}b_{i}\left(f(x_{i})+\nabla
f(x_{i})(x-x_{i})\right)+\lambda^{(1)}_{t}\frac{\|x-x_{0}\|^{2}}{2}+\lambda^{(2)}_{t}\frac{\|x-x_{0}\|^{3}}{6}\right).$

The parameters $b_{i}\geq 0$, $\lambda^{(1)}_{t}$, $\lambda_{t}^{(2}$ and the
iterates $X_{t}$ are designed by theory to ensure the following properties,

$B_{t}f(x_{t})\leq\min_{x}\phi_{t}(x),\qquad\phi(x)\leq
B_{t}f(x)+\frac{\tilde{\lambda}^{(1)}+\lambda^{(1)}_{t}}{2}\|x-x_{0}\|^{2}+\frac{\tilde{\lambda}^{(2)}+\lambda^{(2)}_{t}}{6}\|x-x_{0}\|^{3},$

where $B_{t}=\sum_{i=0}^{t}b_{i}$ and $\tilde{\lambda}^{(1)}$,
$\tilde{\lambda}^{(2)}$ are constants determined by the theory.

Once the parameters are set, the accelerated algorithm operates as follow:

  1. 1.

The accelerated algorithm combines linearly $v_{t}$, the optimum of
$\Phi_{t}$, and the previous iterate $x_{t}$.

  2. 2.

It uses a slight modified version of algorithm 4, see algorithm 6.

  3. 3.

There is a distinction between small and large step sizes, identifying which
$\lambda$ needs to be updated. The step size is considered "large" if it
resembles a cubic-Newton step.

Algorithm 6 Type-I subroutine with backtracking for the accelerated method

First-order oracle $f$, matrices $G,\,D$, vector $\varepsilon$, iterate $x$,
initial smoothness parameter $M_{0}$

Initialize $M\leftarrow\frac{M_{0}}{2}$,
$\texttt{ExitFlag}\leftarrow\texttt{None}$

Define
$\gamma_{M}\stackrel{{\scriptstyle\text{def}}}{{=}}\frac{\kappa_{D}}{\|D\|}\left(\frac{3}{2}\|\varepsilon\|+2\frac{(\|I-P)G\|}{M}\right)$

do

$M\leftarrow 2\cdot M\;\;$ and
$\;\;H_{\gamma}\leftarrow\frac{G^{T}D+D^{T}G}{2}+D^{T}D\frac{M\gamma_{M}}{2}$

$\alpha^{*}\leftarrow\operatorname*{{arg\,min}}_{\alpha}f(x)+\nabla
f(x)^{T}D\alpha+\frac{1}{2}\alpha^{T}H_{\gamma}\alpha+\frac{M\|D\alpha\|^{3}}{6}$

$x_{+}\leftarrow x+D\alpha$

if $\frac{2}{3^{3/4}}\frac{\|\nabla f(x_{+})\|^{3/2}}{\sqrt{M}}\leq-\nabla
f(x_{+})^{T}D\alpha$ then

$\texttt{ExitFlag}\leftarrow\texttt{LargeStep}$

end if

if $\frac{\|\nabla f(x_{+})\|^{2}}{M(\gamma_{M}+\|D\alpha\|)}\leq-\nabla
f(x_{+})^{T}D\alpha$ And $\|D\alpha\|\leq(\sqrt{3}-1)\gamma_{M}$ then

$\texttt{ExitFlag}\leftarrow\texttt{SmallStep}$

end if

while ExitFlag is None

return $x_{+}$, $\alpha$, $M$, $\gamma_{M}$, ExitFlag

Algorithm 7 Adaptive Accelerated Type-I Iterative Algorithm

First-order oracle $f$, initial iterate and smoothness $x_{0},\,M_{0}$, number
of iterations $T$.

$\lambda^{(1)}_{0}\leftarrow 0$, $\lambda^{(2)}_{0}\leftarrow 0$

Initialize $G_{0}$, $D_{0}$, $\varepsilon_{0}$ (see section 2.4)

$\\{x_{1},\,M_{1}\\}\leftarrow\texttt{[\lx@cref{creftype~refnum}{alg:type1}]}(f,G_{0},D_{0},\varepsilon_{0},x_{0},M_{0})$

Initialize $\ell_{0}^{(0)}=f(x_{1})$, $\quad\ell_{0}^{(1)}=0$

for $t=1,\,\ldots,\,T-1$ do

Update $G_{t}$, $D_{t}$, $\varepsilon_{t}$ (see section 2.4)

Set $b_{t}\leftarrow\frac{(t+1)(t+2)}{2}$,
$B_{t}\leftarrow\frac{t(t+1)(t+2)}{6}$, $\beta_{t}\leftarrow\frac{3}{t+3}$.

Update $\ell_{t}^{(0)}\leftarrow\ell_{t-1}^{(0)}+b_{t-1}[f(x_{t})-\nabla
f(x_{t})^{T}x_{t}]$,
$\quad\ell_{t}^{(1)}\leftarrow\ell_{t-1}^{(1)}+b_{t-1}\nabla f(x_{t})$

do

$\texttt{ValidBound}\leftarrow\texttt{True}$

Set $v_{t}\leftarrow\operatorname*{{arg\,min}}_{v}\phi_{t}(v)$ (See
proposition 1).

Let $y_{t}\leftarrow\frac{3}{t+3}v_{t}+\frac{t}{t+3}x_{t}$
$\\{x_{t+1},\alpha_{t},\,M_{t+1},\,\gamma_{t},\texttt{ExitFlag}\\}\leftarrow\texttt{[Alg.\ref{alg:subroutine_acc}]}(f,G_{t},D_{t},\varepsilon_{t},y_{t},\frac{M_{t}}{2})$

%% Check if the next $\phi$ is still a lower bound for $B_{t}f(x_{t+1})$

Define $\phi_{+}(x)=\phi_{t}(x)+b_{t}[f(x_{t+1}+\nabla
f(x_{t+1})(x-x_{t+1})]$.

Set $v_{+}\leftarrow\operatorname*{{arg\,min}}_{v}\phi_{+}(v)$ (See
proposition 1).

if $\Phi_{+}(v_{+})\leq B_{t}f(x_{t+1})$ then %% Parameters adjustment if
needed

$\texttt{ValidBound}\leftarrow\texttt{False}\quad$ %% Unsuccessful iteration:
$\phi_{t+1}(v_{t+1})\geq f(x_{t+1})$.

if ExitFlag is LargeStep then

If  $\lambda_{t}^{(2)}=0$ then
$\displaystyle\lambda_{t}^{(2)}\leftarrow\frac{4}{\sqrt{3}}\frac{b_{t+1}^{3}}{B^{2}_{t}}M_{t+1}$.
Else, $\lambda_{t}^{(2)}\leftarrow 2\lambda_{t}^{(2)}$.

else %% Exitflag is SmallStep

If $\lambda_{t}^{(1)}=0$ then
$\displaystyle\frac{b_{t+1}^{2}}{B_{t}}M_{t+1}\left(\gamma_{t}+\|D_{t}\alpha_{t}\|\right)$.
Else,  $\lambda_{t}^{(1)}\leftarrow 2\lambda_{t}^{(1)}$.

end if

else

$\\{\lambda_{t+1}^{(1)},\lambda_{t+1}^{(2)}\\}\leftarrow\\{\lambda_{t}^{(1)},\lambda_{t}^{(2)}\\}$
%% Successful iteration

end if

while ValidBound is False

end for

return $x_{T}$

######  Proposition 1.

Let $v_{t}$ be the minimizer of

$\phi_{t}(v)=\ell_{t}^{(0)}+\left[\ell_{t}^{(1)}\right]^{T}v+\frac{\lambda_{t}^{(1)}}{2}\|v-x_{0}\|^{2}+\frac{\lambda_{t}^{(2)}}{6}\|v-x_{0}\|^{3}.$

where $\lambda_{t}^{(1)}\geq 0$, $\lambda_{t}^{(2)}\geq 0$. Let
$r_{t}=\|v_{t}-x_{0}\|$. Then,

$\displaystyle r_{t}$
$\displaystyle=\textstyle\|v_{t}-x_{0}\|=\begin{cases}0&\texttt{if
}\lambda_{t}^{(1)}=\lambda_{t}^{(2)}=0\\\
\frac{\|\ell_{t}^{(1)}\|}{\lambda_{t}^{(1)}}&\texttt{if
}\lambda_{t}^{(1)}>0\texttt{ and }\lambda_{t}^{(2)}=0\\\
\frac{-\lambda^{(1)}_{t}+\sqrt{[\lambda^{(1)}_{t}]^{2}+2\lambda^{(2)}_{t}\|\ell_{k}\|}}{\lambda^{(2)}_{2}}&\texttt{if
}\lambda_{t}^{(2)}>0\\\ \end{cases}$ $\displaystyle v_{t}$
$\displaystyle\textstyle=\operatorname*{{arg\,min}}\Phi_{t}(x)=x_{0}-r_{t}\frac{\ell_{t}^{(1)}}{\|\ell_{t}^{(1)}\|}$

##  Appendix B Related work

###  B.1 Inexact, Subspace, and Stochastic Newton Methods

Instead of explicitly computing the Hessian matrix and the Newton step,
inexact methods compute an approximation using sampling [3], inexact Hessian
computation [36, 24], or random subspaces [26, 39, 44]. These approaches
substantially reduce per-iteration costs without significantly compromising
the convergence rate. The convergence speed in such cases often represents an
interpolation between the rates observed in gradient descent methods and
(cubic) Newton’s method.

###  B.2 Nonlinear and Anderson Acceleration

Nonlinear acceleration techniques, including Anderson acceleration [2], have a
long standing history [4, 5, 35]. Driven by their promising empirical
performance, they recently gained interest in their convergence analysis [80,
32, 79, 47, 85, 83, 90, 87, 72, 81, 82, 7, 76, 9, 73]. In essence, Anderson
acceleration is an optimization technique that enhances convergence by
extrapolating a sequence of iterates using a combination of previous gradients
and corresponding iterates. Comprehensive reviews and analyses of these
techniques can be found in notable sources such as [47, 8, 46, 45, 6, 22].
However, these methods do not generalize well outside quadratic minimization
and their convergence rate can only be guaranteed asymptotically when using a
line-search or regularization techniques [78, 84, 72].

###  B.3 Quasi-Newton Methods

Quasi-Newton schemes are renowned for their exceptional efficiency in
continuous optimization. These methods replace the exact Hessian matrix (or
its inverse) in Newton’s step with an approximation updated iteratively during
the method’s execution. The most widely used algorithms in this category
include DFP [23, 31] and BFGS [77, 38, 30, 11, 10]. Most of the existing
convergence results predominantly focus on the asymptotic super-linear rate of
convergence [86, 41, 13, 12, 20, 28, 94, 92, 93]. However, recent research on
quasi-Newton updates has unveiled explicit and non-asymptotic rates of
convergence [65, 67, 66, 54, 55]. Nonetheless, these analyses suffer from
several significant drawbacks, such as assuming an infinite memory size and/or
requiring access to the Hessian matrix. These limitations fundamentally
undermine the essence of quasi-Newton methods, typically designed to be
Hessian-free and maintain low per-iteration cost through their low memory
requirement and low-rank structure.

###  B.4 Close Related Work

####  B.4.1 (Accelerated) Quasi-Newton with Secant Inexactness

Recently, [52] introduced an adaptive regularization technique combined with
cubic regularization, with global, explicit (accelerated) convergence rates
for any quasi-Newton method. Based on the secant inexactness inequality, the
technique introduces a quadratic regularization whose parameter is found by a
backtracking line search. However, this algorithm relies on prior knowledge of
the Lipschitz constant specified in Assumption 1. Unfortunately, the paper
does not provide an adaptive method to find jointly the Lipschitz constant as
well, as it is a priory too costly to know which parameter to update. This
aspect makes the method impractical in real-world scenarios.

####  B.4.2 ARC: Adaptive Regularization algorithm using Cubics

In [18, 19] is proposed a generic framework for inexact cubic regularized
Newton’s steps,

$x_{t+1}=\min_{x}f(x_{t})+\nabla
f(x_{t})(x-x_{t})+\frac{1}{2}(x-x_{t})H_{t}(x-x_{t})+\frac{M_{t}}{6}\|x-x_{t}\|^{3},$

where $H_{t}$ is assumed to be an approximation of the Hessian
$\nabla^{2}f(x_{t})$. However, the theoretical analysis presents numerous
problems, in particular, the assumption that the norm of the current step
bounds the approximation

$\|\nabla^{2}f(x_{t})-H_{t}\|\leq C\|x_{t+1}-x_{t}\|,$

for some constant $C$. Follow up works, such as [91], relaxed this assumption
into

$\|\nabla^{2}f(x_{t})-H_{t}\|\leq C\|x_{t}-x_{t-1}\|,$

which is much weaker since it can be verified while computing the step
$x_{t+1}$. Nevertheless, those are assumptions on the matrix $H_{t}$, but
those works do not explicitly construct such a matrix. Even worse - the
assumption might not be met in practice, especially if $H_{t}$ is a subspace
estimation of the matrix $\nabla^{2}f(x_{t})$.

####  B.4.3 Proximal Quasi-Newton Methods

The work of [68, 37] combined qN methods with proximal schemes and provided
sublinear and accelerated convergence rates. However, the rates in [68] are
based on a technical assumption [68, Assumption 2], for which the authors
commented that "Exploring different conditions on the Hessian approximations
that ensure Assumption 2 is a subject of a separate study", and acknowledge in
their conclusion that "Our framework does not rely on or exploit the accuracy
of second-order information, and hence we do not obtain fast local convergence
rates."

In a follow-up work, [37] proposed accelerated convergence rates under similar
assumptions. However, the authors acknowledge the following: "In our numerical
results, we construct $H_{k}$ via L-BFGS and ignore condition
$\sigma_{k+1}H_{k+1}\preceq\sigma_{k}H_{k}$, since enforcing it in this case
causes a very rapid decrease in $\sigma$. It is unclear, however, if a
practical version of Algorithm 5, based on L-BFGS Hessian approximation, can
be derived, which may explain why the accelerated version of our algorithm
does not represent any significant advantage." In addition, their theoretical
convergence results are based on an upper bound on the sequence $\sigma_{k}$,
which current qN schemes cannot ensure.

####  B.4.4 Proximal Extragradient Quasi-Newton Methods with Online Estimation

Based on the technique in [49], [50] developed a novel quasi-Newton method
with the global accelerated rate of convergence of
$O(\min\\{\frac{1}{t^{2}};\frac{\sqrt{d\log t}}{t^{2.5}}\\})$. The main ideas
are as follows: the authors used the framework of inexact proximal method from
[57], used an online algorithm to estimate the Hessian, and then solved a
linear system involving this approximation using conjugate gradients.

The paper focuses on a different regime than this study: [50] explicitly show
that it is possible to break the $O(\frac{1}{t^{2}})$ barrier for first order
methods using full memory qN methods but this implies storing a full $d\times
d$ matrix, and using it in a linear system, leading to per-iteration
complexities of at least $O(d^{2})$.

From a practical point of view, the algorithm requires numerous
hyperparameters such as $\alpha_{1}$, $\alpha_{2}$, $\beta$,$\ldots$, whose
impact on the efficiency is rather unclear. Moreover, numerically, the
algorithm improves over Nesterov’s acceleration but is slower than l-BFGS on
toy experiments.

##  Appendix C Known rates of convergence and Comparison

###  C.1 (Accelerated) Gradient Descent

This section study the rate of gradient decent when function is smooth (i.e.,
has Lipschitz continuous gradients):

$f(y)\leq f(x)+\nabla f(x)(y-x)+\frac{\mathcal{L}}{2}\|y-x\|^{2},$ (12)

Note that the class of functions considered in this paper is not the class of
smooth functions. However, if the function satisfies Assumption 1, the
Lipchits constant can be bounded as

$\mathcal{L}\leq\|\nabla^{2}f(x)\|+LR\qquad\text{for all $x\in\\{x:f(x)\leq
f(x_{0})\\}$.}$ (13)

The rates of plain gradient descent and its accelerated version read [60]
(after replacing $\mathcal{L}$)

$\displaystyle\min_{0\leq i\leq t}\|\nabla f(x_{i})\|$
$\displaystyle\leq\sqrt{\frac{\left[\|\nabla^{2}f(x)\|+LR\right](f(x_{0})-f^{\star})}{t+1}},$
(plain, non-convex) (14) $\displaystyle f(x_{t})-f(x^{\star})$
$\displaystyle\leq\left[\|\nabla^{2}f(x)\|+LR\right]\frac{2}{t+4}R^{2},$
(plain, convex) (15) $\displaystyle f(x_{t})-f(x^{\star})$
$\displaystyle\leq\left[\|\nabla^{2}f(x)\|+LR\right]\frac{4}{(t+2)^{2}}R^{2}.$
(accelerated) (16)

###  C.2 (Accelerated) Cubic Regularized Newton’s Method

When the function has a Lipschitz-continuous Hessian, the cubic regularized
Newton method and its accelerated version converge with the following rates
[61, 59, 44]:

$\displaystyle\min_{0\leq i\leq t}\|\nabla f(x_{i})\|$
$\displaystyle\leq\frac{16L}{9}\left(\frac{3(f(x_{0})-f^{\star})}{2tM_{\min}}\right)^{2/3},$
(plain, non-convex) (17) $\displaystyle f(x_{t})-f(x^{\star})$
$\displaystyle\leq\quad 9L\frac{R^{3}}{(t+4)^{2}},$ (plain, convex) (18)
$\displaystyle\mathbb{E}[f(x_{t})]-f(x^{\star})$
$\displaystyle\leq\left(\frac{d-N}{N}\right)\frac{\mathcal{L}(3R)^{2}}{2t}+\left(\frac{d}{N}\right)^{2}\frac{L(3R)^{3}}{3t^{2}}+O\left(\frac{1}{t^{3}}\right),$
(Random Subspace, convex) (19) $\displaystyle f(x_{t})-f(x^{\star})$
$\displaystyle\leq L\frac{14R^{3}}{t(t+1)(t+2)}.$ (accelerated) (20)

###  C.3 Relation Between Parameters

Given that this paper does not make the assumption of Lipchitz-continuous
gradients, it becomes necessary to establish connections between various
quantities to facilitate the comparison of rates. To streamline the notation,
all numeric constants are substituted with the big $O$ notation, and the
subsequent equations are derived for the "orthogonal forward estimate only"
update rule, hence $\|D\|=1$ and $\kappa=1$.

##### Relation between $\delta$ and $R$.

The constant $\delta$ represents the upper bound on the relative error (see
Requirement 2):

$\forall t,\;\;\frac{\|\varepsilon_{t}\|}{\|D_{t}\|}\leq\delta.$

For a fixed memory, and assuming $h$ small, since $\varepsilon$ is the norm
between iterates, $\delta$ is upper-bounded as

$\delta\leq O(R).$ (21)

##### Relation between the different $C_{i}$ and $\mathcal{L}$

The $C_{1},C_{2}$, and $C_{4}$ in Theorems 5, 6 and 8 quantifies the
estimation error of $D_{t}^{T}\nabla^{2}f(x_{t})D_{t}$ by $H_{t}$ in (Type I
bound) into two terms:

$C_{i}\leq O\big{(}\delta L+\textstyle\max_{i\leq
t}\|(I-P_{i})\nabla^{2}f(x_{i})\|\Big{)}.$

The first term is the error caused by approximating $\nabla^{2}f(x)D_{t}$ by
$G_{t}$, and the second is the subspace approximation error of
$\nabla^{2}f(x_{t})$ in the span of the columns of $D_{t}$.

Intuitively, the constants $C_{i}$ can be seen as an approximation of an upper
bound on $\mathcal{L}$ in a neighborhood of size $\delta$. This is similar to
eq. 13 but the norm of the Hessian is taken in a subspace, hence the $C_{i}$’s
are smaller. Indeed, using eq. 21, in the worst case, if all iterates
satisfies $\|x_{i}-x^{\star}\|<R$,

$C_{i}=O(RL+\textstyle\max_{i\leq t}\|(I-P_{i})\nabla^{2}f(x_{i})\|).$ (22)

##### Other updates

Note that eqs. 21 and 22 are valid only for the "orthogonal forward estimate
only" update rule. If the random orthogonal forward estimate, or the
orthogonalization of the "greedy" or "iterates only" update rules were used,
the results would have been

$\delta=O(h),\qquad C_{i}=O(hL+\textstyle\max_{i\leq
t}\|(I-P_{i})\nabla^{2}f(x_{i})\|),$

where $h$ is small. However, the comparison with gradient descent or Newton’s
method wouldn’t have been fair as the orthogonalization update rules requires
$N$ additional gradient calls.

###  C.4 Comparing rates of convergence

##### Non convex

The rate from Theorem 5 reads

$\displaystyle\min_{i=1,\,\ldots,\,t}\|\nabla f(x_{i})\|\leq$
$\displaystyle\max\Big{\\{}\frac{3L}{t^{2/3}}\left(12\frac{f(x_{0})-f^{\star}}{M_{\min}}\right)^{2/3}\;;\;\left(\frac{C_{1}}{t^{1/3}}\right)\left(12\frac{f(x_{0})-f^{\star}}{M_{\min}}\right)^{1/3}\Big{\\}},$

where $\textstyle C_{1}=\frac{3\delta
L}{2}+\max_{i\in[0,t]}\|(I-P_{i})\nabla^{2}f(x_{i})P_{i}\|.$ In the case where
$C_{1}$ is small, the rate matches exactly eq. 17. In the other case, using
the approximation from eq. 22,

$\displaystyle\min_{i=1,\,\ldots,\,t}\|\nabla
f(x_{i})\|\leq\left(\frac{O(RL+\textstyle\max_{i\leq
t}\|(I-P_{i})\nabla^{2}f(x_{i})\|)}{t^{1/3}}\right)\left(12\frac{f(x_{0})-f^{\star}}{M_{\min}}\right)^{1/3}$

which differs significantly from eq. 14, as the rate is
$O(\frac{1}{\sqrt{t}})$. However, this might be an artifact of the theoretical
analysis, since the function was not assumed to be smooth.

##### Star convex

After using the approximation from eq. 22, the rate from Theorem 6 reads

$\displaystyle f(x_{t})-f^{\star}\leq
O\left(\frac{f(x_{0})-f^{\star}}{t^{3}}\right)+O\left(\frac{LR^{3}}{t^{2}}\right)+O\left(\frac{[RL+\textstyle\max_{i\leq
t}\|(I-P_{i})\nabla^{2}f(x_{i})\|]R^{2}}{t}\right)$ (23)

The term in $t^{-2}$ is exactly the one from eq. 18, while the term is
$t^{-1}$ has the same dependency in $R^{3}$ compared to eq. 15. However,
$\|(I-P)\nabla^{2}f(x_{i})\|$ could be much smaller than $\|\nabla^{2}f(x)\|$.

##### Convex with random coordinates or random subspace

The rate from Theorem 7 reads

$\displaystyle\mathbb{E}_{D_{t}}[f(x_{t})-f^{\star}]\leq\frac{1}{1+\frac{1}{4}\left[\frac{N}{d}t\right]^{3}}(f(x_{0})-f^{\star})+\frac{1}{\left[\frac{N}{d}t\right]^{2}}\frac{L(3R)^{3}}{2}+\frac{1}{\left[\frac{N}{d}t\right]}\frac{[O(\delta
L)+\frac{(d-N)}{d}\max\limits_{i\in[0,t]}\|\nabla^{2}f(x_{i})\|](3R)^{2}}{2}.$

The rate is similar to eq. 19, up to an additional $O(\delta L/t)$ term. This
extra term comes from the estimation of the Hessian with finite difference,
while the method presented in [44] uses exact Hessian-vector products.

##### Convex, accelerated rates

After using the approximation from eq. 22, and ignoring the terms
$\tilde{\lambda}^{(1)}$, $\tilde{\lambda}^{(2)}$ for clarity, the rate from
Theorem 8 reads

$\displaystyle f(x_{t})-f^{\star}\leq[RL+\max_{i=0\ldots t}\|(I-P_{i})\nabla
f(x_{i})\|]\frac{(3R)^{2}}{(t+3)^{2}}+9\max\left\\{M_{0}\;;\;2L\right\\}\left(\frac{3R}{t+3}\right)^{3}$

The rate is exactly a combination of eq. 20 and eq. 16, but the constant
ascociated to the $1/t^{2}$ rate is smaller in practice: eq. 21 is a
conservative bound and $\|(I-P_{i})\nabla^{2}f(x)\|\leq\|\nabla^{2}f(x)\|$.

##  Appendix D Link with quasi-Newton and Anderson/Nonlinear Acceleration

This section presents the fundamentals of Anderson/nonlinear acceleration
(section D.1), quasi-Newton schemes (section D.2), and their relationship with
the method proposed in this paper (section D.3).

###  D.1 Anderson Acceleration and Nonlinear Acceleration

Anderson acceleration, also known as nonlinear acceleration, is a powerful
technique that enhances the convergence speed of fixed point iterations and
optimization algorithms. Initially developed for solving linear systems,
Anderson acceleration has gained popularity due to its effectiveness in
accelerating iterative methods, including the ones in optimization. The method
leverages previous iterations to construct an improved estimate of the
objective function’s minimizer.

The Anderson acceleration algorithm employs the following approximation to
compute weights:

$\nabla
f\left(\sum_{i=0}^{N}\beta_{i}x_{i}\right)\approx\sum_{i=0}^{N}\beta_{i}\nabla
f(x_{i}),\;\;\sum_{i=0}^{N}\beta_{i}=1.$

When the function $f$ is quadratic, this approximation becomes an equality.
The underlying idea is as follows: since the optimum satisfies $\nabla
f(x^{\star})=0$,

$\sum_{i=0}^{N}\beta_{i}\nabla f(x_{i})\approx 0\;\;\Rightarrow\nabla
f\left(\sum_{i=0}^{N}\beta_{i}x_{i}\right)\approx
0\quad\Rightarrow\sum_{i=0}^{N}\beta_{i}x_{i}\approx x^{\star}.$

The Anderson acceleration steps are thus given by

$x_{t+1}=\sum_{i=0}^{N}\beta_{i}^{\star}x_{t-i+1},\quad\beta^{\star}=\operatorname*{{arg\,min}}_{\beta}\|\sum_{i=0}^{N}\beta_{i}\nabla
f(x_{t-i+1})\|^{2}$

Over the past decades, the ideas behind Anderson acceleration have been
refined. For example, the constraint can be eliminated by considering the step
$x_{t+1}-x_{t}$ instead:

$\displaystyle x_{t+1}-x_{t}$
$\displaystyle=\left(\sum_{i=0}^{N}\beta_{i}x_{t-i+1}\right)-x_{t}$
$\displaystyle=\sum_{i=0}^{N}\tilde{\beta}_{i}x_{t-i+1}.$

The vector $\tilde{\beta}_{i}$ has the property that its sum equals zero.
Hence, it can be rewritten as

$\displaystyle x_{t+1}-x_{t}$
$\displaystyle=\sum_{i=1}^{N}\alpha_{i}(x_{t-i+1}-x_{t-i})$
$\displaystyle\alpha$
$\displaystyle=\operatorname*{{arg\,min}}_{\alpha}\left\|\nabla
f(x_{t})+\sum_{i=1}^{N}\alpha_{i}(\nabla f(x_{t-i+1})-\nabla
f(x_{t-i}))\right\|$

where $\alpha\in\mathbb{R}^{N}$ has no constraint. By writing
$d_{i}=x_{t-i+1}-x_{t-i}$, $g_{i}=\nabla f(x_{t-i+1})-\nabla f(x_{t-i})$, and
$D=[d_{t},\ldots,d_{t-N+1}],\,G=[g_{t},\,\ldots,\,g_{t-N+1}]$, the step
becomes

$\displaystyle x_{t+1}-x_{t}$
$\displaystyle=D_{t}\alpha,\quad\alpha=\operatorname*{{arg\,min}}_{\alpha}\|\nabla
f(x_{t})+G_{t}\alpha\|.$

However, this version of Anderson acceleration is non-convergent because there
is no contribution from $\nabla f(x_{t})$ in the step $x_{t+1}-x_{t}$. The
most popular solution to this problem is introducing a mixing parameter that
combines gradient steps, resulting in the following expression:

$\displaystyle x_{t+1}$ $\displaystyle=x_{t}-h\nabla
f(x_{t})+(D-hG)\alpha,\quad\alpha=\operatorname*{{arg\,min}}_{\alpha}\|\nabla
f(x_{t})+G\alpha\|.$ (AA Type II)

Following a similar idea, recent works have introduced a type I variant of the
algorithm [29, 90, 95, 14] that minimizes the function value instead of the
gradient norm:

$x_{t+1}=x_{t}-h\nabla
f(x_{t})+(D-hG)\alpha,\quad\alpha=\operatorname*{{arg\,min}}f(x_{t})+\nabla
f(x_{t})D_{t}\alpha+\frac{1}{2}\alpha^{T}D_{t}^{T}G_{t}\alpha,$ (AA Type I)

By incorporating regularization [72, 14], globalization techniques [95], or
performing a line search on the parameter $h$, the algorithm converges towards
$x^{\star}$.

###  D.2 Single-secant and Multisecant Quasi-Newton Methods

Quasi-Newton methods, such as the Broyden-Fletcher-Goldfarb-Shanno (BFGS)
method, approximate the Hessian matrix to solve unconstrained optimization
problems efficiently. These methods avoid the expensive computation of the
exact Hessian by using iterative updates based on previous iterates and
gradients of the objective function.

This section focuses on other commonly used quasi-Newton methods: the Davidon-
Fletcher-Powell (DFP) and Broyden type-1 and type-2 updates.

####  D.2.1 The Ideas Behind Single-Secant and Multisecant Hessian
Approximation

In quasi-Newton methods, the Hessian approximation is updated using the secant
equation, which relates the gradients and Hessian at two different points. For
a twice continuously differentiable function, the secant equation is given by:

$\nabla f(y)-\nabla f(x)=\nabla^{2}f(\xi)(y-x),$

where $\xi$ is a point on the line segment connecting $x$ and $y$. This
equation serves as the basis for updating the Hessian approximation.

Based on this remarkable identity, quasi-Newton methods update an
approximation of the Hessian $B_{t}$ or its inverse $H_{t}$ such that the
approximation satisfies

$\nabla f(x_{t})-\nabla f(x_{t-1})=B_{t}(x_{t}-x_{t-1}),\quad
H_{t}\left(\nabla f(x_{t})-\nabla f(x_{t-1})\right)=x_{t}-x_{t-1}.$

What distinguishes the different updates is how to fix the remaining degrees
of freedom. For instance, the simple SR-1 method updates $H_{t}$ such that

$\min_{H}\|H-H_{t-1}\|_{F}\quad:H=H^{T},\;\;H\left(\nabla f(x_{t})-\nabla
f(x_{t-1})\right)=x_{t}-x_{t-1}.$ (24)

Those methods are called single-secant as they update $H_{t}$ only one secant
equation at a time. Hence, in general, $H_{t}$ only satisfies the latest
secant equation.

Multisecant updates, on the other hand, approximate the Hessian using a batch
of secant equations. By introducing matrices
$D_{t}=[x_{t-N+1}-x_{t-N},\ldots,x_{t}-x_{t-1}]$ and $G_{t}=[\nabla
f(x_{t-N+1})-\nabla f(x_{t-N}),\ldots,\nabla f(x_{t})-\nabla f(x_{t-1})]$, the
multisecant updates satisfy

$G_{t}=B_{t}D_{t},\quad\text{or}\quad H_{t}G_{t}=D_{t}.$

Unfortunately, when imposing symmetry, it is impossible to satisfy multiple
secants at a time [70]. However,it is possible to enforce symmetry while
approximating the secant equation in a least square sense [71, 75].

When symmetry is not imposed, the solution for $B_{t}$ and $H_{t}$ can be
obtained as:

$B_{t}=G_{t}[D_{t}]^{\dagger}+B_{0}(I-D_{t}D_{t}^{\dagger}),\quad
H_{t}=D_{t}[G_{t}]^{\dagger}+H_{0}(I-G_{t}G_{t}^{\dagger}),$ (25)

where $B_{0}$ and $H_{0}$ are the initial approximations, and $[A]^{\dagger}$
denotes the pseudo-inverse of matrix $A$. Different choices of pseudo-inverse
lead to different methods.

The inversion of $B_{t}$ can be computed using the Woodbury matrix identity,
which provides an efficient way to compute the inverse. The update for
$B_{t}^{-1}$ is given by:

$B_{t}^{-1}=B_{0}^{-1}\left(I-G_{t}\left(D_{t}^{\dagger}B_{0}^{-1}G_{t}\right)^{-1}D_{t}^{\dagger}B_{0}^{-1}\right)+D_{t}\left(D_{t}^{\dagger}B_{0}^{-1}G_{t}\right)^{-1}D_{t}^{\dagger}B_{0}^{-1}.$

This update is equivalent to the update for $H_{t}$, given that

$B_{0}^{-1}=H_{0},\quad\text{and}\quad
G_{t}^{\dagger}=\left(D_{t}^{\dagger}B_{0}^{-1}G_{t}\right)^{-1}D_{t}^{\dagger}B_{0}^{-1}.$
(26)

In summary, quasi-Newton methods update the Hessian approximation using the
secant equation. Single-secant methods update the approximation using the
secant equation one by one, while multisecant methods use a batch of secant
equations. The choice of updating strategy and pseudo-inverse affects the
behavior of the method.

####  D.2.2 Davidon-Fletcher-Powell (DFP) Formula

The DFP formula is a Quasi-Newton update rule used to iteratively refine an
approximation of the inverse Hessian matrix. It is defined as follows:

$H_{t}=H_{t-1}+\frac{d_{t}d_{t}^{T}}{d_{t}^{T}g_{t}}-\frac{H_{t-1}g_{t}g_{t}^{T}H_{t-1}}{g_{t}^{T}H_{t-1}g_{t}},$
(27)

In the above equation, $g_{t}=\nabla f(x_{t})-\nabla f(x_{t-1})$ represents
the difference in gradients, and $d_{t}=x_{t}-x_{t-1}$ denotes the difference
in parameter values. The DFP formula updates the matrix $H_{t}$ using a rank-
two matrix such that it remains symmetric and positive definite.

####  D.2.3 Multisecant Broyden Methods

The multisecant Broyden methods utilize the update equation from eq. 25, where
$A^{\dagger}$ is chosen as the Moore-Penrose pseudo-inverse of $A$, given by
$A^{\dagger}=(A^{T}A)^{-1}A$. In this equation, $B_{0}$ and $H_{0}$ are scaled
identity matrices. After simplification, the two types of updates can be
expressed as follows:

$\displaystyle B_{t}^{-1}$
$\displaystyle=D_{t}\left(D_{t}^{\dagger}G_{t}\right)^{-1}D_{t}^{\dagger}+B_{0}^{-1}\left(I-G_{t}\left(D_{t}^{\dagger}G_{t}\right)^{-1}D_{t}^{\dagger}\right),$
(28) $\displaystyle H_{t}$
$\displaystyle=D_{t}(G_{t}^{T}G_{t})^{-1}G_{t}^{T}+H_{0}\left(I-G_{t}\left(G_{t}^{T}G_{t}\right)^{-1}G_{t}^{T}\right).$
(29)

Both updates are quite similar, differing mainly in the choice of the pseudo-
inverse of the matrix $G$.

####  D.2.4 Link with Anderson Acceleration

The connection between quasi-Newton methods and Anderson Acceleration is
strong, as, for instance, Broyden methods and Anderson acceleration are
equivalent. To illustrate this, let’s closely examine the update of $\alpha$
in eq. AA Type I:

$\displaystyle x_{t+1}=x_{t}-h\nabla
f(x_{t})+(D_{t}-hG_{t})\alpha,\quad\alpha=\operatorname*{{arg\,min}}f(x_{t})+\nabla
f(x_{t})D_{t}\alpha+\frac{1}{2}\alpha^{T}D_{t}^{T}G_{t}\alpha$
$\displaystyle\Leftrightarrow$ $\displaystyle x_{t+1}=x_{t}-h\nabla
f(x_{t})+(D_{t}-hG_{t})\alpha,\quad\alpha:D_{t}^{T}\nabla
f(x_{t})+D_{t}^{T}G_{t}\alpha=0$ $\displaystyle\Leftrightarrow$ $\displaystyle
x_{t+1}=x_{t}-h\nabla
f(x_{t})+(D_{t}-hG_{t})\alpha,\quad\alpha:\alpha=-(D_{t}^{T}G_{t})^{-1}D_{t}^{T}\nabla
f(x_{t})$ $\displaystyle\Leftrightarrow$ $\displaystyle x_{t+1}=x_{t}-h\nabla
f(x_{t})-(D_{t}-hG_{t})(D_{t}^{T}G_{t})^{-1}D_{t}^{T}\nabla f(x_{t}).$
$\displaystyle\Leftrightarrow$ $\displaystyle
x_{t+1}=x_{t}-\left(D_{t}(D_{t}^{T}G_{t})^{-1}D_{t}^{T}+h\left(I-G_{t}(D_{t}^{T}G_{t})^{-1}D_{t}^{T}\right)\right)\nabla
f(x_{t})$

The above step is precisely the quasi-Newton step
$x_{t+1}=x_{t}-B_{t}^{-1}\nabla f(x_{t})$, where $B_{t}^{-1}$ corresponds to
the Broyden update given by Equation 28, with $B_{0}^{-1}=hI$. A similar
reasoning can be applied to Equation 29.

When considering the single-secant updates, following the same reasoning as in
Section 3 leads to the same conclusion for the SR-1 and DFP updates.

This result is expected since the approximations $H_{t}$ or $B_{t}^{-1}$
satisfy the single or multisecant equation:

$H_{t}G_{t}=D_{t}.$

This indicates that the matrix $H_{t}$ maps vectors from the span of previous
gradients to the span of previous directions. This observation justifies the
construction in eq. 4.

###  D.3 Links with Algorithms 4 and 5

Both Algorithms 4 and 5 can be viewed as quasi-Newton and Anderson/nonlinear
acceleration schemes. The update formulas are

$\displaystyle\min_{\alpha}f(x_{t})+\nabla
f(x_{t})^{T}D_{t}\alpha+\frac{\alpha^{T}H_{t}\alpha}{2}+\frac{M\|D_{t}\alpha\|^{3}}{6},\quad
H_{t}\stackrel{{\scriptstyle\text{def}}}{{=}}\frac{G_{t}^{T}D_{t}+D_{t}^{T}G_{t}+\mathrm{I}M\|D_{t}\|\|\varepsilon_{t}\|}{2}.$
(Type I) $\displaystyle\min_{\alpha}\left\|\nabla
f(x_{t})+G_{t}\alpha\right\|+\frac{M}{2}\Big{(}\sum_{i=1}^{N}|\alpha_{i}|[\varepsilon_{t}]_{i}+\|D_{t}\alpha\|^{2}\Big{)},$
(Type II)

The resemblance with Anderson/nonlinear acceleration is strong, as the
objective function is similar. If the function is quadratic, $L=0$ and
therefore $M$ can also be set to $0$; hence, the coefficients $\alpha$ are
exactly the type I and type II Anderson steps eqs. AA Type I and AA Type II.

The same idea holds when compared to quasi-Newton methods. In both cases, the
optimal solution $\alpha^{\star}$ can be written implicitly:

$\displaystyle\alpha^{\star}=-\left(H_{t}+\frac{MD_{t}^{T}D_{t}\|D_{t}\alpha^{\star}\|}{6}\right)^{-1}D_{t}^{T}\nabla
f(x_{t}),$ (Type I - solution)
$\displaystyle\alpha^{\star}=-\left(G_{t}^{T}G_{t}+\tilde{M}D_{t}^{T}D_{t}\right)^{-1}\left(G_{t}^{T}\nabla
f(x)+\frac{\tilde{M}\|\varepsilon_{t}\|}{2}\partial(|\alpha^{\star}|)\right),$
(Type II - solution)

where $\tilde{M}\stackrel{{\scriptstyle\text{def}}}{{=}}\|\nabla
f(x_{t})+G_{t}\alpha\|M$ and $\partial(|\alpha^{\star}|)$ is a subgradient of
$|\alpha^{*}|$. The step then reads

$\displaystyle x_{t+1}$ $\displaystyle=x_{t}+D\alpha^{\star}$ (Generic step)
$\displaystyle x_{t+1}$
$\displaystyle=x_{t}-D_{t}\left(H_{t}+\frac{MD_{t}^{T}D_{t}\|D_{t}\alpha^{\star}\|}{6}\right)^{-1}D_{t}^{T}\nabla
f(x_{t}),$ (Type I - step) $\displaystyle x_{t+1}$
$\displaystyle=x_{t}-D_{t}\left(G_{t}^{T}G_{t}+\tilde{M}D_{t}^{T}D_{t}\right)^{-1}\left(G_{t}^{T}\nabla
f(x)+\frac{\tilde{M}\|\varepsilon_{t}\|}{2}\partial(|\alpha^{\star}|)\right),$
(Type II - step)

Type I is a quasi-Newton step with a symmetrization of $G^{T}D$ and a
regularization. In contrast, the type II step can be seen as a quasi-Newton
method with a regularization on $G^{\dagger}$, with a correction term on the
gradient. Therefore the Hessian approximation reads

$B_{t}^{-1}=D_{t}\left(H_{t}+\frac{MD_{t}^{T}D_{t}\|D_{t}\alpha^{\star}\|}{6}\right)^{-1}D^{T},\quad
H_{t}=D_{t}\left(G_{t}^{T}G_{t}+\tilde{M}D_{t}^{T}D_{t}\right)^{-1}G_{t}^{T}.$

Again, when the objective function is quadratic, $L=0$ and therefore $M=0$.
Moreover, when $f$ is quadratic, the matrix multiplication $D^{T}G$ satisfies
$D^{T}G+G^{T}D=2D^{T}G$ as $D^{T}G$ becomes symmetric. Hence,

$\displaystyle x_{t+1}$
$\displaystyle=x_{t}-D_{t}\left(D_{t}^{T}G_{t}\right)^{-1}D_{t}^{T}\nabla
f(x_{t}),$ (Type I - quadratic) $\displaystyle x_{t+1}$
$\displaystyle=x_{t}-D_{t}\left(G_{t}^{T}G_{t}\right)^{-1}G_{t}^{T}\nabla
f(x_{t}),$ (Type II quadratic)

The steps are exactly the type I and type II multisecant Broyden methods from
eqs. 28 and 29, with the only difference that there is no initialization
$H_{0}$ or $B_{0}$.

##  Appendix E Solving the sub-problems

##### Solving the Type 1 Subproblem

The Type 1 subproblem is a well-studied problem that involves minimizing a
specific objective function. A method proposed by [61] has proven to be
efficient for solving this problem. The method utilizes eigenvalue
decomposition on a matrix to find the optimal solution. In this paper, the
matrix involved in this problem is relatively small, therefore eigenvalue
decomposition is not a concern even for large-scale problems. The subproblem
aims to determine the norm of the solution, and this can be achieved through
solving one nonlinear equation using bisection or secant method.

##### Solving the Type 2 Subproblem

The Type 2 subproblem can be formulated as a Second-Order Cone Program (SOCP).
The objective function of this subproblem consists of three terms: a norm
term, a sum of absolute values term, and a quadratic term. The norm term can
be transformed using singular value decomposition, and the sum of absolute
values term can be expressed as with linear constraints. The quadratic term
can be simplified using a rotated quadratic cone. By utilizing these
techniques, the Type 2 subproblem can be effectively solved using existing
SOCP solvers.

###  E.1 Solving the Type 1 Subproblem

The Type 1 subproblem can be expressed as follows:

$\min_{\alpha}\nabla
f(x)D\alpha+\frac{1}{2}\alpha^{T}H\alpha+\frac{M}{6}\|D\alpha\|^{3},$

where $H$ is symmetric but not necessarily positive definite. This problem has
been well-studied, and [61] proposed an efficient method to solve it using
eigenvalue decomposition on the matrix $H$. Although eigenvalue decomposition
may be challenging for large-scale problems, it is not a concern here since
$H\in\mathbb{R}^{N\times N}$, with a relatively small $N$ (e.g., $N=25$ in the
experiments).

In essence, the subproblem involves determining the norm of the solution
$r=\|\alpha\|$. This can be accomplished through a simple bisection on the
following system of nonlinear equations:

$\left(H+\frac{MD^{T}Dr}{2}I\right)\alpha=-D^{t}\nabla
f(x),\quad\|\alpha\|=r,\quad r\geq-\lambda_{\min}(H).$ (30)

Interestingly, this problem is equivalent to the following formulation, as
shown in Proposition 2:

$\left(\Lambda+\frac{Mr}{2}I\right)\tilde{\alpha}=-V^{T}(D^{T}D)^{-1/2}D^{t}\nabla
f(x),\;\|\alpha\|=r,\;r\geq-\lambda_{\min}(H),\;\tilde{\alpha}=V^{T}(D^{T}D)^{1/2}\alpha,$
(31)

which involves the eigenvalue decomposition
$(D^{T}D)^{-1/2}H(D^{T}D)^{-1/2}=V\Lambda V^{T}$.

######  Proposition 2.

Problems eq. 30 and eq. 31 are equivalent.

###### Proof.

The first step is to split $D^{T}D=(D^{T}D)^{1/2}(D^{T}D)^{1/2}$ and then
employ an eigenvalue decomposition on
$(D^{T}D)^{-1/2}H(D^{T}D)^{-1/2}=V\Lambda V^{T}$ (where $V$ is orthonormal due
to the symmetry of the matrix):

$\displaystyle\left(H+\frac{MD^{T}Dr}{2}I\right)\alpha=-D^{t}\nabla f(x)$
$\displaystyle\Leftrightarrow$
$\displaystyle(D^{T}D)^{1/2}\left((D^{T}D)^{-1/2}H(D^{T}D)^{-1/2}+\frac{Mr}{2}I\right)(D^{T}D)^{1/2}\alpha=-D^{t}\nabla
f(x)$ $\displaystyle\Leftrightarrow$
$\displaystyle(D^{T}D)^{1/2}V\left(\Lambda+\frac{Mr}{2}I\right)V^{T}(D^{T}D)^{1/2}\alpha=-D^{t}\nabla
f(x)$ $\displaystyle\Leftrightarrow$
$\displaystyle\left(\Lambda+\frac{Mr}{2}I\right)V^{T}(D^{T}D)^{1/2}\alpha=-V^{T}(D^{T}D)^{-1/2}D^{t}\nabla
f(x)$ $\displaystyle\Leftrightarrow$
$\displaystyle\left(\Lambda+\frac{Mr}{2}I\right)\tilde{\alpha}=-V^{T}(D^{T}D)^{-1/2}D^{t}\nabla
f(x).$

∎

Once the eigenvalue decomposition is performed, the subproblem eq. 31 becomes
relatively simple since it involves solving a diagonal system of equations for
a fixed value of $r$. The main objective is to find an interval
$[r_{\min},r_{\max}]$ that encompasses the optimal value $r=\|\alpha\|$. Once
this interval is identified, a straightforward bisection or secant method can
be employed to obtain the optimal solution.

##### Finding initial bounds

Starting with $r_{\min}=\max\\{0,-\lambda_{\min(H)}\\}$ and
$r_{\max}=\max\\{2r_{\min},1\\}$,

$\text{do }r_{\max}\leftarrow 2r_{\max}\quad\text{while
}\left\|\tilde{\alpha}\right\|\geq r_{\max}.$

where
$\tilde{\alpha}=-\left(\Lambda+\frac{Mr_{\max}}{2}I\right)^{-1}V^{T}(D^{T}D)^{-1/2}D^{t}\nabla
f(x)$. Increasing $r_{\max}$ increases the regularization, hence reduces the
norm of $\tilde{\alpha}$.

##### Finding $\alpha$

After $r^{\star}$ has been found such that $|r^{\star}-\|\tilde{\alpha}\||$ is
sufficiently small, the best $\alpha$ is simply

$\alpha=(D^{T}D)^{-1/2}V\tilde{\alpha}=-(D^{T}D)^{-1/2}V\left(\Lambda+\frac{Mr^{\star}}{2}I\right)^{-1}V^{T}(D^{T}D)^{-1/2}D^{t}\nabla
f(x).$

In the case where the diagonal matrix is not invertible, which happens when
$r^{\star}=r_{\min}$, it suffices to use the pseudo-inverse instead.

Note that $D^{T}D$ is an $N\times N$ matrix, where $N$ is small, therefore,
computing its inverse is inexpensive. Moreover, when $D$ is orthogonal,
$D^{T}D=I$, therefore there is no need to invert it. In addition,
$(\Lambda+\frac{Mr^{\star}}{2}I)^{-1}$ can be computed in $O(N)$ complexity
since the matrix is diagonal.

###  E.2 Solving the Type 2 Subproblem

The Type 2 subproblem is given by:

$\min_{\alpha}\underbrace{\left\|\nabla
f(x)+G\alpha\right\|}_{\textbf{(a)}}+\frac{L}{2}\Big{(}\underbrace{\sum_{i=1}^{N}|\alpha_{i}|\varepsilon_{i}}_{\textbf{(b)}}+\underbrace{\|D\alpha\|^{2}}_{\textbf{(c)}}\Big{)}.$
(32)

Although it may not be immediately apparent, this subproblem can be formulated
as a Second-Order Cone Program (SOCP) with $O(N)$ variables and constraints.

####  E.2.1 Fundamentals of SOCP

SOCP solvers handle the following conic problems:

$\displaystyle\min_{x,t_{i},\omega_{i}}c_{0}x+\sum_{i}c_{i}[t_{i};\omega_{i}]\quad\text{subject
to}$ $\displaystyle A_{0}x+\sum_{i=1}^{k}A_{i}[t_{i};\omega_{i}]=b$ (SOCP
Standard Matrix Form) $\displaystyle x\geq 0$
$\displaystyle(t_{i},\omega_{i})\in\mathcal{K}_{i}\;\;\Leftrightarrow
t_{i}\geq\|\omega_{i}\|,\;\;t\geq 0.$

Here, $k$ represents the number of cones, and the cone $\mathcal{K}$ refers to
the second-order cone, also known as the Lorenz cone.

A useful transformation is the rotated quadratic cone, defined as follows:

$\displaystyle[a,b,c]\in\mathcal{K}_{q}\quad\Leftrightarrow\quad
2ab\geq\|c\|^{2}.$

The rotated quadratic cone can be reformulated as a second-order cone using a
linear transformation:

$\text{if}\quad\begin{bmatrix}a\\\ b\\\
c\end{bmatrix}=\begin{bmatrix}\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}}&0\\\
\frac{1}{\sqrt{2}}&-\frac{1}{\sqrt{2}}&0\\\
0&0&I_{K}\end{bmatrix}\begin{bmatrix}t\\\ \omega^{(0)}\\\
\omega\end{bmatrix}\quad\text{then}\quad(t,[\omega^{(0)};\,\omega])\in\mathcal{K}\quad\Leftrightarrow\quad[a,b,c]\in\mathcal{K}_{q}.$

Thanks to this transformation, the rotated quadratic cone can be included in
SOCP solvers.

####  E.2.2 SOCP Formulation of the Type 2 Subproblem

The SOCP of eq. 32 is composed of the three terms a, b, and c.

##### Term (a)

Let $U_{G}\Sigma_{G}V^{T}_{G}$ be the singular value decomposition of $G$.
Write $P_{G}=U_{G}U_{G}^{T}$ as the projector onto the columns of $G$. Then,

$\displaystyle\left\|\nabla f(x)+R\alpha\right\|$
$\displaystyle=\left\|P_{G}\nabla f(x)+P_{G}G\alpha+(I-P_{G})\nabla
f(x)\right\|$ $\displaystyle=\sqrt{\left\|P_{G}\nabla
f(x)+R\alpha\right\|^{2}+\|(I-P_{G})\nabla f(x)\|^{2}}$
$\displaystyle=\sqrt{\left\|U_{G}\left(U_{G}^{T}\nabla
f(x)+\Sigma_{G}V^{T}_{G}\alpha\right)\right\|^{2}+\|(I-P_{G})\nabla
f(x)\|^{2}}$ $\displaystyle=\sqrt{\left\|U_{G}^{T}\nabla
f(x)+\Sigma_{G}V^{T}_{G}\alpha\right\|^{2}+\|(I-P_{G})\nabla f(x)\|^{2}}$

Let the vector $\omega_{1}=\left[U_{G}^{T}\nabla
f(x)+\Sigma_{G}V\alpha;\,\|(I-P_{G})\nabla f(x)\|\right]$. Hence,

$\left\|\nabla
f(x)+G\alpha\right\|=\min_{t_{1},\,\alpha,\,\omega_{1}}t_{1}:(t_{1},\omega_{1})\in\mathcal{K}_{L},\quad\omega_{1}=\left[U_{G}^{T}\nabla
f(x)+\Sigma_{G}V\alpha;\,\|(I-P_{G})\nabla f(x)\|\right].$

##### Term (b)

This term is standard in linear programming. Let
$\alpha=\alpha_{+}-\alpha_{-}$, with $\alpha_{+},\alpha_{-}\geq 0$,

$\sum_{i=1}^{N}|\alpha_{i}|\varepsilon_{i}=\sum_{i=1}^{N}(\alpha_{+}+\alpha_{-})\varepsilon_{i}.$

##### Term (c)

Let $U_{D}\Sigma_{D}V^{T}_{D}$ be the singular value decomposition of $D$.
Using the rotated cone, the constraint can be written as

$2t_{3}b\geq\|U_{D}\Sigma_{D}V_{D}\alpha\|^{2}=\|\Sigma_{D}V_{D}\alpha\|^{2},\quad
b=\frac{1}{2}.$

Using the transformation into a Lorenz cone, this is equivalent to

$\begin{bmatrix}1&0&0\\\ 0&1&0\\\
0&0&\Sigma_{D}V_{D}^{T}\end{bmatrix}\begin{bmatrix}t_{3}\\\ b\\\
\alpha\end{bmatrix}=\begin{bmatrix}\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}}&0\\\
\frac{1}{\sqrt{2}}&-\frac{1}{\sqrt{2}}&0\\\
0&0&I_{k}\end{bmatrix}\begin{bmatrix}t_{2}\\\ \omega^{(0)}_{2}\\\
\omega_{2}\end{bmatrix},\quad
b=\frac{1}{2},\quad(t_{2},\,[\omega^{(0)}_{2},\,\omega_{2}])\in\mathcal{K}.$

Simplification. Note that, since $b=\frac{1}{2}$, the value can be immediately
replaced. Same idea with $t_{3}$: the constraint is written as

$t_{3}=\frac{t_{2}+\omega^{(0)}_{2}}{\sqrt{2}},\quad t_{3}\geq 0.$

Since, by construction, $t_{2}\geq\omega^{(0)}_{2}$ and $t_{2}\geq 0$, $t_{3}$
always satisfies the condition, which means both $t_{3}$ and its constraint
can be removed. The constraints thus simplify into

$\begin{bmatrix}\frac{1}{2}\\\ 0\end{bmatrix}+\begin{bmatrix}0\\\
\Sigma_{D}V_{D}^{T}\end{bmatrix}\begin{bmatrix}\alpha\end{bmatrix}=\begin{bmatrix}\frac{1}{\sqrt{2}}&-\frac{1}{\sqrt{2}}&0\\\
0&0&I_{k}\end{bmatrix}\begin{bmatrix}t_{2}\\\ \omega^{(0)}_{2}\\\
\omega_{2}\end{bmatrix},\quad(t_{2},\,[\omega^{(0)}_{2},\,\omega_{2}])\in\mathcal{K}.$

##### Final formulation

Gathering all terms, the final SOCP formulation reads

minimize $\displaystyle
t_{1}+\frac{L}{2}\left((\alpha_{+}+\alpha_{-})^{T}\varepsilon+t_{2}\right)$
subject to $\displaystyle\omega_{1}=\left[U_{G}^{T}\nabla
f(x)+\Sigma_{G}V_{G}^{T}\alpha\;;\;\|(I-P_{G})\nabla f(x)\|\right],$
$\displaystyle\alpha_{+},\alpha_{-}\geq 0$
$\displaystyle\alpha=\alpha_{+}-\alpha_{-}$
$\displaystyle\begin{bmatrix}\textbf{0}_{1\times
N}&-\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}}&0\\\
\Sigma_{D}V_{D}^{T}&\textbf{0}_{N\times 1}&\textbf{0}_{N\times
1}&-I_{N}\end{bmatrix}\begin{bmatrix}\alpha\\\ t_{2}\\\ \omega^{(0)}_{2}\\\
\omega_{2}\end{bmatrix}=\begin{bmatrix}0\\\ -\frac{1}{2}\\\
\textbf{0}_{N\times 1}\end{bmatrix}$
$\displaystyle(t_{1},\,\omega_{1})\in\mathcal{K},\quad(t_{2},[\omega^{(0)}_{2};\omega_{2}])\in\mathcal{K}_{L},\quad
t_{2}\geq 0.$

##### Standard matrix formulation

The SOCP can be written under the standard matrix form eq. SOCP Standard
Matrix Form. Let the variables

$\alpha_{+},\,\alpha_{-}\geq
0,\quad(t_{1},\omega_{1})\in\mathcal{K}_{1},\quad(t_{2},[\omega_{2}^{(0)}\omega_{2}])\in\mathcal{K}_{2},$

where $t_{1}$, $t_{2}$, and $\omega_{2}^{(0)}$ are scalars, $\omega_{2}$,
$\alpha_{+}$, and $\alpha_{-}$ are vectors of size $N$, and $\omega_{1}$ is a
vector of size $N+1$. The SOCP matrices read

$\displaystyle
c_{0}=\begin{bmatrix}\frac{L\varepsilon^{T}}{2}&\frac{L\varepsilon^{T}}{2}\end{bmatrix}\quad
c_{1}=\begin{bmatrix}1&\textbf{0}_{1\times N+1}\end{bmatrix}\quad
c_{2}=\begin{bmatrix}\frac{L}{2\sqrt{2}}&\frac{L}{2\sqrt{2}}&\textbf{0}_{1\times
N}\end{bmatrix}$ $\displaystyle
A_{0}=\begin{bmatrix}-\Sigma_{G}V_{G}^{T}&\Sigma_{G}V_{G}^{T}\\\
\textbf{0}_{2\times N}&\textbf{0}_{2\times N}\\\
\Sigma_{D}V_{D}^{T}&-\Sigma_{D}V_{D}^{T}\end{bmatrix}$ $\displaystyle
A_{1}=\begin{bmatrix}\textbf{0}_{N+1\times 1}&I_{N+1\times N+1}\\\
\textbf{0}_{N+1\times 1}&\textbf{0}_{N+1\times N+1}\\\ \end{bmatrix}$
$\displaystyle A_{2}=\begin{bmatrix}\textbf{0}_{N+1\times
1}&\textbf{0}_{N+1\times 1}&\textbf{0}_{N+1\times N}\\\
-\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}}&\textbf{0}_{1\times N}\\\
\textbf{0}_{N\times 1}&\textbf{0}_{N\times 1}&-I_{N\times N}\\\ \end{bmatrix}$
$\displaystyle b=\begin{bmatrix}\nabla f(x)^{T}U_{G}&\|(I-P_{R})\nabla
f(x)\|&-\frac{1}{2}&\textbf{0}_{N\times 1}\end{bmatrix}^{T}.$

This completes the SOCP formulation of the type 2 subproblem.

##  Appendix F Additional Numerical Experiments

This section presents additional numerical experiments.

##### Methods

The methods compared are the type 1 and type 2 steps with the following
strategies: Iterate only, Forward estimate only, Greedy (refer to section
2.4), and the accelerated type 1 method with the strategy forward estimate
only. The batch methods are not included as they perform poorly regarding the
number of Oracle calls. The baseline is the l-BFGS method from minFunc [69].

##### Method parameters

In all experiments, the memory of the methods is set to $N=25$ and the $h$ for
the forward estimates is set to $h=10^{-9}$. The parameters of the l-BFGS are
left untouched except for the memory. The initial point is $x_{0}=\nabla
f(0_{d})$.

##### Functions

The minimized problems are square loss with cubic regularization, logistic
loss with small quadratic regularization, and the generalized Rosenbrock
function. The regularization parameter of the square loss is set to $1e-3$
times the norm of the Hessian, and the regularization of the logistic loss is
set to $1e-10$ times the square norm of the feature matrix.

##### Dataset

The datasets for the square and the logistic loss are Madelon [42], Sido0
[43], and Marti2 [43] datasets.

##### Post-processing

The dataset matrix is normalized by its norm, then a vector of ones is
concatenated to the data matrix.

###  F.1 Initial Parameter for the Backtracking Line search

The backtracking line search was used in all experiments. The estimation of
the initial value $M_{0}$ (see eq. 33) is based on the following observation.
Since the function satisfies Assumption 1,

$\|\nabla f(y)-\nabla f(x)-\nabla^{2}f(x)(y-x))\|\leq\frac{L}{2}\|y-x\|^{2},$

for some $x,\,y$. Hence, the parameter $L$ can be estimated as

$L\approx 2\frac{\|\nabla f(y)-\nabla
f(x)-\nabla^{2}f(x)(y-x))\|}{\|y-x\|^{2}}.$

Now, define

$\displaystyle s_{h}\stackrel{{\scriptstyle\text{def}}}{{=}}h\nabla
f(x_{0}),\qquad$

for some small $h$ and $\tau>1$, and let $x=x_{0}$ and $y=x_{0}+s_{\tau h}$.
Indeed, if $h$ is small, then

$\tau\left[\nabla f(x_{0}+s_{h})-\nabla
f(x_{0})\right]\approx\tau\nabla^{2}f(x)s_{h}=\nabla^{2}f(x)s_{\tau h}.$

Therefore,

$\|\nabla f(x_{0}+s_{\tau h})-\nabla f(x_{0})-\tau\left[\nabla
f(x_{0}+s_{h})-\nabla f(x_{0})\right]\|\approx\|\nabla f(x_{0}+s_{\tau
h})-\nabla f(x_{0})-\nabla^{2}f(x)s_{\tau h}\|,$

and hence, the Lipchitz constant can be estimated as

$\displaystyle M_{0}=\frac{2}{\|s_{\tau h}\|^{2}}\|\nabla f(x_{0}+s_{\tau
h})-\nabla f(x_{0})-\tau\left[\nabla f(x_{0}+s_{h})-\nabla f(x_{0})\right]\|.$
(33)

In the experiments, $h$ is the same as the algorithm, and $\tau=10$. Various
choices of $\tau,\,h$ have been tested without significantly impacting the
numerical convergence.

###  F.2 Scalability w.r.t. Dimension and Memory

Figure 2: Scaling of the Type 1 method with the "orthogonal forward estimates
only" updates rules w.r.t. $N$ and $d$ to minimize a random logistic
regression function. As predicted by the theory, the scaling is linear in the
dimension and quadratic w.r.t. $N$. The proposed method is suitable for large-
scale problems, as it can quickly solve problems with $d\approx 10^{6}$.

Figure 3: Distribution of the per-iteration time for three methods. The memory
parameter of l-BFGS and the type I method is set to (left to right)
$N=5,25,100$. The time required by the l-BFGS algorithm increases slightly
when $N$ grows, and the per-iteration computation time is approximately two
times faster than the type I method. Surprisingly, the total computation time
of the type-1 method remains constant for different $N$ because the condition
in the backtracking line search is more often satisfied. Note that the $\times
2$ factor between l-BFGS and the type 1 method is expected since the type 1
method requires at least 2 gradient calls.

###  F.3 Influence of h

Figure 4: Influence of the step size $h$ to compute the forward estimate
$x_{+\frac{1}{2}}$ in the "orthogonal forward estimates only" updates rules on
the Madelon dataset to minimize a (left) quadratic and (right) a logistic
loss. The range of acceptable $h$ is rather large. For instance, this range is
$[10^{-9},10^{-1}]$ when minimizing the logistic loss.

###  F.4 Impact of the memory parameter N

Figure 5: Impact of the memory size $N$ on the convergence rate of the type 1
method with the "Orthogonal forward estimate" update rule to minimize a
logistic loss on the Madelon dataset. Left: number of iterations versus
suboptimality, right: time versus suboptimality. Overall, it is always better
to increase the memory parameter in terms of the number of iterations, but
there is an effect of diminishing returns.

###  F.5 Nonconvex optimization

Figure 6: Comparison of type 1 methods on the Generalized Rosenbrock function
in $\mathbb{R}^{100}$.

Figure 7: Comparison of type 2 methods on the Generalized Rosenbrock function
in $\mathbb{R}^{100}$

###  F.6 Comparison of Type 1 Methods on Convex Problems

####  F.6.1 Square loss and cubic regularization

Figure 8: Comparison of type 1 methods: Square loss and cubic regularization
on Madelon dataset

Figure 9: Comparison of type 1 methods: Square loss and cubic regularization
on sido0 dataset

Figure 10: Comparison of type 1 methods: Square loss and cubic regularization
on marti2 dataset

####  F.6.2 Logistic regression

Figure 11: Comparison of type 1 methods: Logistic loss and cubic
regularization on Madelon dataset

Figure 12: Comparison of type 1 methods: Logistic loss and cubic
regularization on sido0 dataset

Figure 13: Comparison of type 1 methods: Logistic loss and cubic
regularization on marti2 dataset

###  F.7 Comparison of Type 2 Methods on Convex Problems

The type-2 method was not the focus of this study. Its prototypical
implementation is rather slow, hence, the time VS suboptimality graph are not
showed.

####  F.7.1 Square loss and cubic regularization

Figure 14: Comparison of type 2 methods: Square loss and cubic regularization
on Madelon dataset

Figure 15: Comparison of type 2 methods: Square loss and cubic regularization
on sido0 dataset

Figure 16: Comparison of type 2 methods: Square loss and cubic regularization
on marti2 dataset

####  F.7.2 Logistic regression

Figure 17: Comparison of type 2 methods: Logistic loss and cubic
regularization on Madelon dataset

Figure 18: Comparison of type 2 methods: Logistic loss and cubic
regularization on sido0 dataset

Figure 19: Comparison of type 2 methods: Logistic loss and cubic
regularization on marti2 dataset

##  Appendix G Missing proofs

In this section, when not needed, the subscript $t$ has been removed for
clarity. The following definitions simplify the notations:

$\displaystyle D_{\dagger}$ $\displaystyle=(D^{T}D)^{-1}D^{T},$ (34)
$\displaystyle D^{T}_{\dagger}$ $\displaystyle=D(D^{T}D)^{-1},$ (35)
$\displaystyle\kappa_{D}$ $\displaystyle=\|D_{\dagger}\|\|D\|,$ (36)

Note that the pseudo inverse $D_{\dagger}$ exists under Requirement 3. Note
that

$D_{\dagger}D=I,\qquad DD_{\dagger}=P_{D}=P.$

###  G.1 Technical Result: Hessian Approximation

This section presents technical results related to the approximation of the
Hessian $\nabla^{2}f(x)$. To simplify notations, let the matrices $H_{0}$ and
$\tilde{H}_{0}$ be

$\displaystyle
H_{0}=\frac{D^{T}R+R^{T}D}{2},\qquad\tilde{H}_{0}=D_{\dagger}^{T}H_{0}D_{\dagger}=\frac{PGD_{\dagger}+D_{\dagger}^{T}G^{T}P}{2}.$
(37)

Intuitively, $\tilde{H}_{0}$ is the Hessian approximation, while $H_{0}$ is
the approximation of the quadratic form $D^{T}\nabla^{2}f(x)D$.

######  Proposition 3 (Subspace Hessian Approximation Error).

Assume $D$ satisfies Requirement 1b. Then, the following holds:

$\left\|\left(\tilde{H}_{0}-P\nabla^{2}f(x)P\right)D\alpha\right\|\leq\frac{L}{2}\|D_{\dagger}\|\|\varepsilon\|\|D\alpha\|$

###### Proof.

Since $D^{T}D_{\dagger}=D_{\dagger}^{T}D^{T}=P$, $D_{\dagger}D=\mathrm{I}$,
$PD=D$, $\|P\|=1$, and using eq. 37,

$\displaystyle\left\|\left[\frac{PGD_{\dagger}+D^{T}_{\dagger}G^{T}P}{2}-P\nabla^{2}f(x)P\right]D\alpha\right\|$
$\displaystyle\leq\quad$
$\displaystyle\frac{1}{2}\big{(}\left\|(PGD_{\dagger}-P\nabla^{2}f(x)P)D\alpha\right\|+\left\|(D^{T}_{\dagger}G^{T}P-P\nabla^{2}f(x)P)D\alpha\right\|\big{)}$
$\displaystyle\leq\quad$
$\displaystyle\frac{1}{2}\Big{(}\left\|G\alpha-\nabla^{2}f(x)D\alpha\right\|+\|D_{\dagger}\|\left\|(G^{T}-D^{T}\nabla^{2}f(x))D\alpha\right\|\Big{)}$

Using inequality eq. 8 for the first term and eq. 9 for second gives

$\displaystyle\left\|\left[\frac{PGD_{\dagger}+D^{T}_{\dagger}G^{T}P}{2}-P\nabla^{2}f(x)P\right]D\alpha\right\|\leq\frac{1}{2}\left(\frac{L}{2}|\alpha|^{T}\varepsilon+\|D_{\dagger}\|\frac{L\|D\alpha\|}{2}\|\varepsilon\|\right)$

Because
$|\alpha|^{T}\varepsilon\leq\|\alpha\|\|\varepsilon\|\leq\|D_{\dagger}\|\|D\alpha\|\|\varepsilon\|$,

$\displaystyle\left\|\left[\frac{PGD_{\dagger}+D^{T}_{\dagger}G^{T}P}{2}-P\nabla^{2}f(x)P\right]D\alpha\right\|\leq\frac{L}{2}\|D_{\dagger}\|\|\varepsilon\|\|D\alpha\|.$

∎

######  Proposition 4.

[Out-of-subspace Error Estimation] Let the function $f$ satisfy Assumption 1.
Let the matrices $D,\,G$ be defined as in eq. 6 and vector $\varepsilon$ as in
eq. 7. Then, for all $\alpha\in\mathbb{R}^{N}$,

$\frac{\|(I-P)\nabla^{2}f(x)D\alpha\|}{\|D\alpha\|}\leq\left(\|(I-P)G\|+L\|\varepsilon\|\right)\frac{\kappa_{D}}{\|D\|}.$

###### Proof.

Indeed, using eq. 9,

$\displaystyle\|(I-P)\nabla^{2}f(x)D\alpha\|$
$\displaystyle=\|(I-P)(G-G+\nabla^{2}f(x)D)\alpha\|$
$\displaystyle\leq\|(I-P)(\nabla^{2}f(x)D-G)\alpha\|+\|(I-P)G\alpha\|$
$\displaystyle\leq\|(\nabla^{2}f(x)D-G)\alpha\|+\|(I-P)G\alpha\|$
$\displaystyle\leq\|\nabla^{2}f(x)D-G\|\|\alpha\|+\|(I-P)G\alpha\|$
$\displaystyle\leq\left(\frac{L\|\varepsilon\|}{2}\|\alpha\|+\|(I-P)G\alpha\|\right)$

Hence,

$\frac{\|(I-P)\nabla^{2}f(x)D\alpha\|}{\|D\alpha\|}\leq\frac{\left(\frac{L\|\varepsilon\|}{2}\|\alpha\|+\|(I-P)G\alpha\|\right)}{\|D\alpha\|}.$

Moreover,

$\displaystyle\frac{\|(I-P)\nabla^{2}f(x)D\alpha\|}{\|D\alpha\|}$
$\displaystyle\leq\left(\frac{L\|\varepsilon\|}{2}+\|(I-P)G\|\right)\|\alpha\|.$
$\displaystyle\leq\left(\frac{L\|\varepsilon\|}{2}+\|(I-P)G\|\right)\frac{\|\alpha\|}{\|D\alpha\|}$
$\displaystyle\leq\max_{\alpha}\left(\frac{L\|\varepsilon\|}{2}+\|(I-P)G\|\right)\frac{\|\alpha\|}{\|D\alpha\|}$
$\displaystyle=\left(\frac{L\|\varepsilon\|}{2}+\|(I-P)G\|\right)\sigma^{-1}_{\min}(D).$

The desired result follows from the fact that
$\kappa_{D}=\frac{\|D\|}{\sigma_{\min}(D)}$. ∎

###  G.2 Technical Results: Cubic Subproblem

This section presents results on the properties of the solution of the cubic
subproblem

$\alpha^{\star}\stackrel{{\scriptstyle\text{def}}}{{=}}\operatorname*{{arg\,min}}_{\alpha}\nabla
f(x)^{T}(D\alpha)+\frac{1}{2}(D\alpha)^{T}\tilde{H}_{\Gamma}(D\alpha)+\frac{M}{6}\|D\alpha\|^{3},\qquad
x_{+}=x+D\alpha^{\star}$ (38)

where $\tilde{H}_{\Gamma}\in\mathbb{R}^{d\times d}$ is a rank $N$ matrix such
that

$\displaystyle\tilde{H}=D_{\dagger}^{T}H_{\gamma}D_{\dagger},\qquad\Leftrightarrow\;\;H=D^{T}\tilde{H}_{\Gamma}D,\qquad
H_{\gamma}=\frac{R^{T}D+D^{T}R+\Gamma}{2},$ (39)

and $\Gamma$ is a $N\times N$ matrix. For instance, setting
$\Gamma=M\|\varepsilon\|\|D\|I$ gives the $H$ in algorithm 4.

######  Proposition 5.

The first-order and second-order conditions of the subproblem eq. 38 read

$\displaystyle D^{T}\nabla
f(x)+H_{\Gamma}\alpha+\frac{M}{2}D^{T}D\alpha\|D\alpha\|$ $\displaystyle=0,$
(40) $\displaystyle H_{\Gamma}+\frac{M}{2}D^{T}D\|D\alpha\|$
$\displaystyle\succeq 0.$ (41)

###### Proof.

See [59], equation (3.3), and [61], equation (2.7). ∎

######  Proposition 6.

Let $f$ satisfies Assumption 1 and $B\in\mathbb{R}^{d\times d}$ be any matrix.
Assume the matrix $D$ satisfies Requirement 1b, and $\alpha$ satisfies the
first-order condition eq. 40. Let $\tilde{H}_{\Gamma}$ be defined in eq. 39.
Then,

$\displaystyle\|\nabla f(x)+BD\alpha-\nabla f(x_{+})\|$
$\displaystyle=\|(\tilde{H}_{\Gamma}-B+\frac{M\|D\alpha\|}{2})D\alpha+\nabla
f(x_{+})\|$ (42)
$\displaystyle\leq\frac{L}{2}\|D\alpha\|^{2}+\|[B-\nabla^{2}f(x)]D\alpha\|.$
(43)

Then, the following equation follows from the optimality condition multiplied
by $D(D^{T}D)^{-1}$, writing $P=DD_{\dagger}=D_{\dagger}^{T}D^{T}$, assuming
$P\nabla f(x)=\nabla f(x)$,

$\nabla f(x)+(\tilde{H}_{\Gamma}+\frac{M\|D\alpha\|}{2})D\alpha=0.$

Replacing $\nabla f(x)$ gives

$\|\nabla f(x)+BD\alpha-\nabla
f(x_{+})\|=\|-(\tilde{H}_{\Gamma}+\frac{M\|D\alpha\|}{2})D\alpha+BD\alpha-\nabla
f(x_{+})\|,$

which is the desired result.

###### Proof.

The inequality follows directly from eq. 2,

$\displaystyle\|\nabla f(x)+BD\alpha-\nabla f(x_{+})\|$
$\displaystyle\leq\|\nabla f(x)+\nabla^{2}f(x)D\alpha-\nabla
f(x_{+})\|+\|BD\alpha-\nabla^{2}f(x)D\alpha\|$
$\displaystyle\leq\frac{L}{2}\|D\alpha\|^{2}+\|[B-\nabla^{2}f(x)]D\alpha\|.$

∎

######  Proposition 7.

Assume $D$ satisfies Requirement 1b. Let $\tilde{H}$ be defined in eq. 39.
Then, for all $\tilde{\Gamma}$, if

$B=\tilde{H}_{\Gamma}-\frac{1}{2}D_{\dagger}\tilde{\Gamma}D_{\dagger}^{T}$

in proposition 6, the following holds:

$\displaystyle\left\|\left(\frac{1}{2}D_{\dagger}\tilde{\Gamma}D_{\dagger}^{T}+\frac{M\|D\alpha\|}{2}\right)D\alpha+\nabla
f(x_{+})\right\|\leq\frac{L}{2}\|D\alpha\|^{2}+\|[B-\nabla^{2}f(x)]D\alpha\|,$
(44)

where

$\displaystyle\|[B-\nabla^{2}f(x)]D\alpha\|\leq\|D\alpha\|\left(\frac{L}{2}\|D_{\dagger}\|\|\varepsilon\|+\frac{\|(I-P)\nabla^{2}f(x)D\alpha\|}{\|D\alpha\|}+\frac{1}{2}\left\|D_{\dagger}(\Gamma-\tilde{\Gamma})D_{\dagger}\right\|\right)$

###### Proof.

From proposition 6,

$\|(\tilde{H}_{\Gamma}-B+\frac{M\|D\alpha\|}{2})D\alpha+\nabla
f(x_{+})\|\leq\frac{L}{2}\|D\alpha\|^{2}+\|[B-\nabla^{2}f(x)]D\alpha\|.$

Replacing $B$ in the left-hand-side gives

$\|(\tilde{H}_{\Gamma}-B+\frac{M\|D\alpha\|}{2})D\alpha+\nabla
f(x_{+})\|=\|(\frac{D_{\dagger}\Gamma
D_{\dagger}^{T}}{2}+\frac{M\|D\alpha\|}{2})D\alpha+\nabla f(x_{+})\|$

Since

$\nabla^{2}f(x)D\alpha=P\nabla^{2}f(x)PD\alpha+(I-P)\nabla^{2}f(x)PD\alpha,$

where $P=D(D^{T}D)^{-1}D^{T}$, and because $PD=D$, the inequality becomes

$\displaystyle\|[B-\nabla^{2}f(x)]D\alpha\|=$
$\displaystyle\|\left[\tilde{H}_{\Gamma}-\frac{1}{2}D_{\dagger}\tilde{\Gamma}D_{\dagger}^{T}-\nabla^{2}f(x)\right]D\alpha\|$
(45) $\displaystyle=$
$\displaystyle\|[P+(I-P)]\left[\tilde{H}_{\Gamma}-\frac{1}{2}D_{\dagger}\tilde{\Gamma}D_{\dagger}^{T}-\nabla^{2}f(x)\right]PD\alpha\|$
(46) $\displaystyle\leq$
$\displaystyle\left\|\left(\tilde{H}_{0}-P\nabla^{2}f(x)P\right)D\alpha\right\|$
(47)
$\displaystyle+\left(\frac{1}{2}\left\|D_{\dagger}^{T}(\Gamma-\tilde{\Gamma})D_{\dagger}\right\|+\frac{\left\|(I-P)\nabla^{2}f(x)D\alpha)\right\|}{\|D\alpha\|}\right)\|D\alpha\|$
(48)

∎

######  Corollary 1 (Bound depending on $\tilde{\Gamma}$).

In proposition 7,

  * •

if $\tilde{\Gamma}=0$ and $\Gamma=M\|D\|\|\varepsilon\|I$,

$\displaystyle\left\|\frac{M\|D\alpha\|}{2}D\alpha+\nabla
f(x_{+})\right\|\leq\frac{L}{2}\|D\alpha\|^{2}+\|D\alpha\|\left(\frac{\|\varepsilon\|}{\|D\|}\left(\frac{L+M\kappa_{D}}{2}\right)\kappa_{D}+\|(I-P)\nabla^{2}f(x)P\|\right)$
(49)

  * •

if $\tilde{\Gamma}=\Gamma$,

$\displaystyle\left\|\left(\frac{1}{2}D_{\dagger}\Gamma
D_{\dagger}^{T}+\frac{M\|D\alpha\|}{2}\right)D\alpha+\nabla
f(x_{+})\right\|\leq\frac{L}{2}\|D\alpha\|^{2}+\|D\alpha\|\left(\frac{L}{2}\frac{\|\varepsilon\|}{\|D\|}\kappa_{D}+\frac{\|(I-P)\nabla^{2}f(x)D\alpha\|}{\|D\alpha\|}\right)$
(50)

  * •

If $\tilde{\Gamma}=D(M\|D\alpha\|)D^{T}$ and $\Gamma=M\|D\|\|\varepsilon\|I$,

$\displaystyle\|\nabla f(x_{+})\|$
$\displaystyle\leq\frac{L+M}{2}\|D\alpha\|^{2}+\|D\alpha\|\left(\frac{\|\varepsilon\|}{\|D\|}\left(\frac{L+M\kappa_{D}}{2}\right)\kappa_{D}+\|(I-P)\nabla^{2}f(x)P\|\right)$
(51)

###  G.3 Technical Results: Decrease Guarantees

This section presents two technical results on the minimal decrease of the
function $f$.

######  Proposition 8.

Let Assumptions 1, 1b, 2 and 3 hold. Then, $\forall y\in\mathbb{R}^{d}$,
algorithm 4 ensures

$f(x_{+})\leq
f(y)+\frac{M+L}{6}\|y-x\|^{3}+\frac{\|y-x\|^{2}}{2}\left(\|\nabla^{2}f(x)-P\nabla^{2}f(x)P\|+\delta\frac{L\kappa+M\kappa^{2}}{2}\right)$

###### Proof.

The output of algorithm 4 ensures that

$\displaystyle f(x_{+})\leq\min_{\alpha}f(x)+\nabla
f(x)^{T}D\alpha+\frac{1}{2}(D\alpha)^{T}\nabla^{2}f(x)D\alpha+\frac{1}{2}\alpha^{T}\left(H-D^{T}\nabla^{2}f(x)D\right)\alpha+\frac{M}{6}\|D\alpha\|^{3}$

However, by the definition of $H$ eq. Type I bound,

$\displaystyle\frac{1}{2}\alpha^{T}\left(H-D^{T}\nabla^{2}f(x)D\right)\alpha$
$\displaystyle\leq$
$\displaystyle\frac{1}{2}\left(\alpha^{T}\left(\frac{G^{T}D+D^{T}G}{2}-D^{T}\nabla^{2}f(x)D\right)\alpha+\|\alpha\|^{2}\frac{M\|D\|\|\varepsilon\|}{2}\right)$
$\displaystyle\leq$
$\displaystyle\frac{1}{2}\left(\alpha^{T}\left(\frac{G^{T}D+D^{T}G}{2}-D^{T}\nabla^{2}f(x)D\right)\alpha+\|D^{\dagger}\|^{2}\|D\alpha\|\frac{M\|D\|\|\varepsilon\|}{2}\right)$
$\displaystyle=$
$\displaystyle\frac{1}{2}\left((D\alpha)^{T}\left(G-\nabla^{2}f(x)D\right)\alpha+\|D^{\dagger}\|^{2}\|D\alpha\|\frac{M\|D\|\|\varepsilon\|}{2}\right).$

The last equality comes from the fact that

$\alpha^{T}\left(D^{T}G\right)\alpha=\alpha^{T}\left(\frac{D^{T}G+G^{T}D}{2}+\frac{D^{T}G-G^{T}D}{2}\right)\alpha=\alpha^{T}\left(\frac{D^{T}G+G^{T}D}{2}\right)\alpha.$

Now, using eq. 8 with $w=D\alpha$ gives

$\displaystyle\frac{1}{2}\alpha^{T}\left(H-D^{T}\nabla^{2}f(x)D\right)\alpha\leq\frac{L\|D\alpha\|}{4}\sum_{i=1}^{N}|\alpha_{i}|\varepsilon_{i}+\|D^{\dagger}\|^{2}\|D\alpha\|\frac{M\|D\|\|\varepsilon\|}{4}.$

Finally, since

$\sum_{i=1}^{N}|\alpha_{i}|\varepsilon_{i}\leq\|\alpha\|\|\varepsilon\|\leq\|D^{\dagger}\|\|D\alpha\|\|\varepsilon\|,$

the inequality becomes

$\displaystyle\frac{1}{2}\alpha^{T}\left(H-D^{T}\nabla^{2}f(x)D\right)\alpha$
$\displaystyle\leq\frac{\|D\alpha\|^{2}}{4}\left(L\|D^{\dagger}\|\|\varepsilon\|+M\|D^{\dagger}\|^{2}\|D\|\|\varepsilon\|\right)$
$\displaystyle=\frac{\|D\alpha\|^{2}}{4}\frac{\|\varepsilon\|}{\|D\|}\left(L\kappa_{D}+M\kappa_{D}^{2}\right).$

All together,

$\displaystyle f(x_{+})$ $\displaystyle\leq$
$\displaystyle\min_{\alpha}f(x)+\nabla
f(x)^{T}D\alpha+\frac{1}{2}(D\alpha)^{T}\nabla^{2}f(x)D\alpha+\frac{1}{2}\alpha^{T}\left(H-D^{T}\nabla^{2}f(x)D\right)\alpha+\frac{M}{6}\|D\alpha\|^{3}$
$\displaystyle\leq$ $\displaystyle\min_{\alpha}f(x)+\nabla
f(x)^{T}D\alpha+\frac{1}{2}(D\alpha)^{T}\nabla^{2}f(x)D\alpha+\frac{\|D\alpha\|^{2}}{4}\frac{\|\varepsilon\|}{\|D\|}\left(L\kappa_{D}+M\kappa_{D}^{2}\right)+\frac{M}{6}\|D\alpha\|^{3}$

Now, by Requirement 3, for all $y$, one can find $\alpha$ such that

$D\alpha=P(y-x)=DD^{\dagger}(y-x).$

Indeed, multiplying both sides by $D^{\dagger}$ gives

$\alpha=D^{\dagger}(y-x).$

Therefore, the minimum can be written as a function of $y$ instead of
$\alpha$,

$\displaystyle f(x_{+})\leq\min_{y\in\mathbb{R}^{d}}\;\;$ $\displaystyle
f(x)+\nabla f(x)^{T}P(y-x)+\frac{1}{2}(P(y-x))^{T}\nabla^{2}f(x)P(y-x)$
$\displaystyle+\frac{\|P(y-x)\|^{2}}{4}\frac{\|\varepsilon\|}{\|D\|}\left(L\kappa_{D}+M\kappa_{D}^{2}\right)+\frac{M}{6}\|P(y-x)\|^{3}.$
(52)

Since $P\nabla f(x)=\nabla f(x)$ by Requirement 1b, and using the crude bound
$\|P(y-x)\|\leq\|y-x\|$,

$\displaystyle f(x_{+})\leq\min_{y\in\mathbb{R}^{d}}\;\;$ $\displaystyle
f(x)+\nabla f(x)^{T}(y-x)+\frac{1}{2}(y-x)^{T}\nabla^{2}f(x)(y-x)$
$\displaystyle+\frac{1}{2}(y-x)\left[\nabla^{2}f(x)-P\nabla^{2}f(x)P\right](y-x)$
$\displaystyle+\frac{\|y-x\|^{2}}{4}\frac{\|\varepsilon\|}{\|D\|}\left(L\kappa_{D}+M\kappa_{D}^{2}\right)+\frac{M}{6}\|y-x\|^{3}.$

Using the lower bound eq. 3,

$f(x)+\nabla
f(x)^{T}(y-x)+\frac{1}{2}(y-x)^{T}\nabla^{2}f(x)(y-x)-\frac{L}{6}\|y-x\|^{3}\leq
f(y),$

the crude bound
$(y-x)\left[\nabla^{2}f(x)-P\nabla^{2}f(x)P\right](y-x)\leq\|\nabla^{2}f(x)-P\nabla^{2}f(x)P\|\|y-x\|^{2}$,
and Requirements 3 and 2 lead to the desired result,

$\displaystyle f(x_{+})\leq
f(y)+\frac{M+L}{6}\|y-x\|^{3}+\frac{\|y-x\|^{2}}{2}\left(\|\nabla^{2}f(x)-P\nabla^{2}f(x)P\|+\delta\frac{L\kappa+M\kappa^{2}}{2}\right)$

∎

######  Proposition 9.

Let Assumptions 1, 1a, 2 and 3 hold. Then, $\forall y\in\mathbb{R}^{d}$,
algorithm 4 ensures

$\displaystyle\mathbb{E}f(x_{+})\leq$
$\displaystyle\left(1-\frac{N}{d}\right)f(x)+\frac{N}{d}f(y)+\frac{N}{d}\frac{(M+L)}{6}\|y-x\|^{3}$
$\displaystyle+\frac{N}{d}\frac{\|y-x\|^{2}}{2}\left(\delta\frac{L\kappa+M\kappa^{2}}{2}+\frac{(d-N)}{d}\|\nabla^{2}f(x)\|\right)$

###### Proof.

The proof is the same as for proposition 8, until equation eq. 52,

$\displaystyle f(x_{+})\leq\min_{y\in\mathbb{R}^{d}}\;\;$ $\displaystyle
f(x)+\nabla f(x)^{T}P(y-x)+\frac{1}{2}(P(y-x))^{T}\nabla^{2}f(x)P(y-x)$
$\displaystyle+\frac{\|P(y-x)\|^{2}}{4}\frac{\|\varepsilon\|}{\|D\|}\left(L\kappa_{D}+M\kappa_{D}^{2}\right)+\frac{M}{6}\|P(y-x)\|^{3}.$

With Requirement 1a, the following relations hold (see [44, lemma 5.7])

$\displaystyle\mathbb{E}[\|P(y-x)\|^{2}]$
$\displaystyle=(y-x)^{T}\mathbb{E}[P](y-x)=\frac{N}{d}\|y-x\|^{2},$ (53)
$\displaystyle\mathbb{E}[\|P(y-x)\|^{3}]$
$\displaystyle\leq\mathbb{E}[\|P(y-x)\|^{2}]\|y-x\|=\frac{N}{d}\|y-x\|^{2},$
(54) $\displaystyle\mathbb{E}[(y-x)^{T}P\nabla^{2}f(x)P(y-x)]$
$\displaystyle\leq\frac{N^{2}}{d^{2}}(y-x)\nabla^{2}f(x)(y-x)+\frac{N(d-N)}{d^{2}}\|\nabla^{2}f(x)\|\|y-x\|^{2}$
(55)

Hence, removing the minimum and taking the expectation of eq. 52 gives

$\displaystyle\mathbb{E}f(x_{+})\leq$ $\displaystyle f(x)+\frac{N}{d}\nabla
f(x)^{T}(y-x)$
$\displaystyle+\frac{1}{2}\left(\frac{N^{2}}{d^{2}}(y-x)\nabla^{2}f(x)(y-x)+\frac{N(d-N)}{d^{2}}\|\nabla^{2}f(x)\|\|y-x\|^{2}\right)$
$\displaystyle+\frac{N}{d}\frac{\|y-x\|^{2}}{4}\frac{\|\varepsilon\|}{\|D\|}\left(L\kappa_{D}+M\kappa_{D}^{2}\right)+\frac{N}{d}\frac{M}{6}\|y-x\|^{3}.$

Using the lower bound from eq. 3

$\frac{1}{2}(y-x)\nabla^{2}f(x)(y-x)\leq
f(y)+\frac{L}{6}\|y-x\|^{3}-f(x)-\nabla f(x)(y-x)$

in the inequality over the expectation gives

$\displaystyle\mathbb{E}f(x_{+})\leq$ $\displaystyle f(x)+\frac{N}{d}\nabla
f(x)^{T}(y-x)$
$\displaystyle+\frac{N^{2}}{d^{2}}\left(f(y)+\frac{L}{6}\|y-x\|^{3}-f(x)-\nabla
f(x)(y-x)\right)$
$\displaystyle+\frac{1}{2}\frac{N(d-N)}{d^{2}}\|\nabla^{2}f(x)\|\|y-x\|^{2}$
$\displaystyle+\frac{N}{d}\frac{\|y-x\|^{2}}{4}\frac{\|\varepsilon\|}{\|D\|}\left(L\kappa_{D}+M\kappa_{D}^{2}\right)+\frac{N}{d}\frac{M}{6}\|y-x\|^{3}.$

After simplification,

$\displaystyle\mathbb{E}f(x_{+})\leq$
$\displaystyle\left(1-\frac{N^{2}}{d^{2}}\right)f(x)+\frac{N^{2}}{d^{2}}f(y)+\frac{N}{d}\left(1-\frac{N}{d}\right)\nabla
f(x)^{T}(y-x)$
$\displaystyle+\frac{1}{2}\frac{N(d-N)}{d^{2}}\|\nabla^{2}f(x)\|\|y-x\|^{2}$
$\displaystyle+\frac{N}{d}\frac{\|y-x\|^{2}}{4}\frac{\|\varepsilon\|}{\|D\|}\left(L\kappa_{D}+M\kappa_{D}^{2}\right)+\left(\frac{N^{2}L}{6d^{2}}+\frac{NM}{6d}\right)\|y-x\|^{3}.$

To simplify the expression, since $N\leq d$,

$\left(\frac{N^{2}L}{6d^{2}}+\frac{NM}{6d}\right)\|y-x\|^{3}\leq\frac{N(M+L)}{6d}\|y-x\|^{3}.$

Finally, since the function is convex,

$\frac{N}{d}\left(1-\frac{N}{d}\right)\nabla
f(x)^{T}(y-x)\leq\frac{N}{d}\left(1-\frac{N}{d}\right)(f(y)-f(x)).$

From this last relation, Requirement 2 and Requirement 3 comes the desired
result,

$\displaystyle\mathbb{E}f(x_{+})\leq$
$\displaystyle\left(1-\frac{N}{d}\right)f(x)+\frac{N}{d}f(y)+\frac{N(M+L)}{6d}\|y-x\|^{3}$
$\displaystyle+\frac{\|y-x\|^{2}}{2}\left(\frac{N}{d}\delta\frac{L\kappa+M\kappa^{2}}{2}+\frac{N(d-N)}{d^{2}}\|\nabla^{2}f(x)\|\right)$

∎

###  G.4 Technical Results: Accelerated Algorithm

##### Notations

The following functions define the estimate sequence,

$\displaystyle\ell_{t}(x)$
$\displaystyle=\sum_{i=2}^{t}b_{i-1}\left(f(x_{i})+\nabla
f(x_{i})(x-x_{i})\right),$ (56) $\displaystyle\phi_{t}(x)$
$\displaystyle=f(x_{1})+\ell_{t}(x)+\frac{\lambda_{t}^{(1)}}{2}\|x-x_{0}\|^{2}+\frac{\lambda_{t}^{(2)}}{6}\|x-x_{0}\|^{3}$
(57) $\displaystyle\Phi_{t}(x)$ $\displaystyle=\frac{\phi_{t}(x)}{B_{t}},$
(58)

where $\lambda_{t}^{(1,2)}$ are non-negative and increasing, and the sequences
$b_{t}$, $B_{t}$ are

$\displaystyle B_{t}$
$\displaystyle=\frac{t(t+1)(t+2)}{6}=\sum_{i=1}^{t}b_{i},$ (59) $\displaystyle
b_{t}$ $\displaystyle=\frac{(t+1)(t+2)}{2}=B_{t+1}-B_{t}.$ (60)

Moreover, the following quantities will be important later,

$\displaystyle v_{t}$
$\displaystyle=\operatorname*{{arg\,min}}_{x}\phi_{t}(x)=\operatorname*{{arg\,min}}_{x}\Phi_{t}(x),$