#  Why Line Search when you can Plane Search?
SO-Friendly Neural Networks allow Per-Iteration Optimization of
Learning and Momentum Rates for Every Layer

Betty Shea<EMAIL_ADDRESS>
University of British Columbia Mark Schmidt<EMAIL_ADDRESS>
University of British Columbia, Canada CIFAR AI Chair (Amii)

###### Abstract

We introduce the class of SO-friendly neural networks, which include several
models used in practice including networks with 2 layers of hidden weights
where the number of inputs is larger than the number of outputs. SO-friendly
networks have the property that performing a precise line search to set the
step size on each iteration has the same asymptotic cost during full-batch
training as using a fixed learning. Further, for the same cost a plane-search
can be used to set both the learning and momentum rate on each step. Even
further, SO-friendly networks also allow us to use subspace optimization to
set a learning rate and momentum rate for each layer on each iteration. We
explore augmenting gradient descent as well as quasi-Newton methods and Adam
with line optimization and subspace optimization, and our experiments indicate
that this gives fast and reliable ways to train these networks that are
insensitive to hyper-parameters.

##  1 Should we use Subspace Optimization in Machine Learning?

To train machine learning (ML) models, we often try to minimize a function $f$
in terms of some parameters $w$. Starting from an initial guess $w_{0}$, a
prototpyical algorithm for performing this minimization is _gradient descent
with momentum_ (GD+M) (Polyak, 1964) which generates a sequence of guesses
according to the formula111With the convention that $w_{-1}=w_{0}$.

$w_{k+1}=w_{k}-\alpha_{k}\nabla f(w_{k})+\beta_{k}(w_{k}-w_{k-1}).$ (1)

This update subtracts the gradient $\nabla f(w_{k})$ multiplied by what we
call the _learning rate_ $\alpha_{k}$, and then adds the momentum term
$(w_{k}-w_{k-1})$ multiplied by what we call the _momentum rate_ $\beta^{k}$.
Many variations exist including stochastic gradient descent (SGD) which uses a
stochastic approximation of the gradient to reduce the iteration cost (Robbins
and Monro, 1951) and quasi-Newton methods which pre-multiply the gradient by a
matrix to improve the convergence rate (see Nocedal and Wright, 2006). The
most popular variant in ML is the Adam optimizer which incorporates a
variation on momentum and uses a scaling of a stochastic approximation of the
gradient (Kingma and Ba, 2014).

In the modern practice of ML, the learning rate $\alpha_{k}$ and momentum rate
$\beta_{k}$ are often set to fixed constants for all $k$. This can result in
poor performance of Adam and related methods (Vaswani et al., 2020), so
practitioners often run the algorithm with different choices of the constant
learning and momentum rates. After running the algorithm with different hyper-
parameter choices, we then choose the run giving the best performance (Zhang
et al., 2022). This “hyper-parameter search” approach to solving an
optimization problem is expensive compared to traditional deterministic
numerical optimization approaches like quasi-Newton methods. Popular quasi-
Newton methods, such as those based on the limited-memory Broyden-Fletcher-
Goldfarb-Shanno (L-BFGS) update (Nocedal, 1980), typically give good
performance “out of the box”. Specifically, deterministic methods often give
good performance without any hyper-parameter tuning using their default hyper-
parameters.

A key ingredient that allows L-BFGS to work well with default parameters is
using a _line search_ (LS) to set the learning rate. In particular, they
search for a value of $\alpha_{k}$ satisfying conditions guaranteeing
sufficient progress in optimizing $f$ such as the strong Wolfe conditions
(Wolfe, 1969; Nocedal and Wright, 2006). This may require testing multiple
learning rates on each iteration, but various heuristics exist that allow us
to quickly find a suitable $\alpha_{k}$ (Moré and Thuente, 1994; Shanno and
Phua, 1978). With these heuristics, the cost of using line search is typically
smaller than running the algorithm for multiple fixed learning rates. Further,
since line searches allow us to find a good learning rate on each iteration,
using line searches has the potential to drastically outperform any fixed
learning rate. Several recent works highlight the advantages of using line
search for modern ML problems within SGD (Truong and Nguyen, 2018; Vaswani et
al., 2019) as well as Adam (Vaswani et al., 2020; Galli et al., 2023).
However, using line search increases the iteration cost of the method and
these methods still use a fixed momentum rate.

In the context of many linear models like least squares and logistic
regression, it is possible to use line searches “for free”. In particular, for
what we call linear-composition problems (LCPs) using a LS or even performing
line optimization (LO) to numerically solve for the best step size does not
change the asymptotic cost of the iterations compared to using a fixed step
size. Further, for this same cost LCPs also allow us to use a _plane search_
to set the learning rate $\alpha_{k}$ and momentum rate $\beta_{k}$ on each
iteration of GD+M (1) and related methods. In Section 2 we review how we can
optimize the learning and momentum rates for LCPs for the same cost as using a
fixed $\alpha_{k}$ and $\beta_{k}$, and how this leads to significantly faster
convergence in practice. Using a plane search is a special case of methods
that use _subspace optimization_ (SO) to search along multiple directions on
each iteration, and in Section 3 we discuss the history of SO methods.

LO and SO are not typically used for neural networks since performing LO or SO
in general neural network model cannot be done “for free” as it can in linear
models with suitable directions. In this work, we introduce the class of SO-
friendly neural networks (Section 4). The key property of these networks is
that performing LO to set the learning rate, or SO to set the learning rate
and momentum rate, does not increase the asymptotic cost of training compared
to using a fixed learning and momentum rate. While SO-friendly networks are a
restricted class, they include some networks used in practice such as networks
with 2 layers of weights where the number of inputs is much larger than the
number of outputs. A further appealing property of SO-friendly networks is
that we can optimize over a separate learning rate and momentum rate for each
layer. For networks with 2 layers of weights this gives us 4 step sizes that
are tuned “for free” on each iteration, and this flexibility can lead to even
faster convergence in practice in some settings. In Section 5 we consider
augmenting L-BFGS and Adam with LO and SO, and show that this leads to faster
training for both LCPs and SO-friendly networks.

###  1.1 List of Contributions by Figure

We have performed a wide variety of experiments that support our main claim
that LO and SO can and should be used to improve the performance of various ML
algorithms for suitable problem structures. Here we give of figures
highlighting various observations from our experiments:

  * •

Figure 2 illustrates how LO and SO can significantly improve the performance
of GD+M on the LCP of logistic regression.

  * •

Figure 5 shows that GD+M with LO and SO can even improve performance over
accelerated gradient methods, but that performing SO over more than 2
directions only results in marginal performance gains.

  * •

Figure 6 shows that LO and SO also lead to improved performance of GD+M for
training 2-layer neural networks.

  * •

Figure 8 shows that using SO to set per-layer step sizes can sometimes
substantially improve performance for training 2-layer networks, but sometimes
makes performance much worse.

  * •

Figure 10 shows that using per-layer step sizes consistently improves
performance if regularization is added.

  * •

Figure 13 and 14 show that LO can improve the performance of a quasi-Newton
method, and that using SO to optimize the step sizes in a quasi-Newton method
with momentum further improves performance.

  * •

Figures 17 and 18 shows that LO improves the performance of Adam for LCPs
though not always for neural networks, but that a multi-direction Adam method
with step sizes set using SO typically improves performance over single-
direction Adam methods.

###  1.2 Comments on Limitations of the Applicability of Subspace Optimization

In this work we restrict attention to deterministic optimization methods, to
the restricted class of SO-friendly networks, and to differentiable objective
functions. Before beginning, we briefly comment on these choices:

  * •

Why focus on deterministic methods? While stochastic methods do indeed become
superior to deterministic as the number of examples increases, in many
practical scenarios we have a finite dataset and we perform multiple passes
through the data to train a model. For example, on the ImageNet dataset it is
standard to perform 60 passes through the data (and this cost may be
multiplied by the number of optimization hyper-parameter settings that are
explored). It is not obvious that stochastic methods necessarily dominate
deterministic methods in this common multi-pass setting, since it is possible
that deterministic methods could exist that find accurate solutions in 10-20
passes through the data. Further, deterministic methods typically do not
require a search over optimization hyper-parameters. The insensitivity of
deterministic methods to hyper-parameter settings also means we can have more
confidence that deterministic optimizers will work robustly in new scenarios.

  * •

Why focus on SO-friendly networks? Most deep neural networks are not SO-
friendly and thus it is not efficient to use SO. Nevertheless, a variety of
SO-friendly networks are used in practice (see Sections 2 and 4) including
many networks with two hidden layers of weights or deep networks where we have
sub-optimized individual layers. There has been recent interest in the ML
literature on improved optimization methods for such networks (Mishkin et al.,
2022), and we believe that practitioners in these settings would substantially
reduce training times (or increase accuracy) by exploiting SO.

  * •

Why focus on differentiable objectives? Some objective functions arising in ML
are non-differentiable, such as support vector machines which are LCPs and
neural networks with rectified linear unit (ReLU) non-linear functions. Note
that SO remains efficient for LCPs and SO-friendly networks even when the
problem is non-differentiable. However, in non-differentiable cases it is less
obvious which directions to use as search directions (since neither the
gradient nor its negation necessarily point in a direction that we can search
in order to improve the objective function).

Nevertheless, in Section 6 we discuss how future work could use SO as a useful
ingredient in stochastic methods and/or methods to train general deep neural
networks.

##  2 Line Search and Plane Search for Linear Composition Problems (LCPs)

Optimizing a function over a low-dimensional subspace is easier than
optimizing the function over the full space. However, for a given problem it
is possible that SO is still too expensive for it to be beneficial.
Nevertheless, for certain problems we can perform LO or SO over a low-
dimensional subspace without increasing the asymptotic iteration cost compared
to using a fixed step size. We call these problems _SO-friendly_. Narkiss and
Zibulevsky (2005a) were the first to highlight an SO-friendly problem, for a
class of LCPs.

LCPs have the form $f(w)=g(Xw)$ for the $n\times d$ data matrix $X$ and a
function $g$. For this problem setting, SO is efficient in the common case
where the matrix multiplication with $X$ is much more expensive than
evaluating $g$. For example, for a target vector $y$ the least squares
objective is $f(w)=\frac{1}{2}\|Xw-y\|^{2}$, which can be written as an LCP
with $g(m)=\frac{1}{2}\|m-y\|^{2}$. In this setting the matrix multiplication
costs $O(nd)$ but evaluating $g$ only costs $O(n)$. Another canonical example
is logistic regression, where $g(m)=\sum_{i=1}^{n}\log(1+\exp(-y_{i}m_{i}))$
for binary labels $y_{i}\in\\{-1,1\\}$. Many other classic models in ML can be
written as LCPs such as robust regression with the Huber loss and support
vector machines. More advanced situations where LCPs arise include include
neural tangent kernels (NTKs) (Jacot et al., 2018) and transfer learning by
re-training the last layer of a pre-trained deep neural network (Donahue et
al., 2014).

###  2.1 Efficient Line Search (LS) and Line Optimization (LO)

The classic gradient descent (GD) iteration takes the form

$w_{k+1}=w_{k}-\alpha_{k}\nabla f(w_{k}).$ (2)

To ensure that the step size $\alpha_{k}$ is sufficient small, standard
deterministic optimizers impose the Armijo sufficient decrease condition

$f(w_{k}-\alpha_{k}\nabla f(w_{k}))\leq f(w_{k})-\alpha_{k}\sigma\|\nabla
f(w_{k})\|^{2},$ (3)

for a positive sufficient decrease factor $\sigma$. LS codes may also ensure
the step size is not too small by imposing an additional curvature condition
as in the strong Wolfe conditions (Nocedal and Wright, 2006). There exist
approaches to finding an $\alpha_{k}$ satisfying these types of conditions
that require only a small number of evaluations of $f$ (Moré and Thuente,
1994). While these methods are effective for general smooth functions, in the
special case of LCPs we can efficiently use LO to set $\alpha_{k}$ in order to
decrease $f$ by a larger amount than LS methods.

For LCPs the gradient has the form $\nabla f(w)=X^{T}\nabla g(Xw)$ and thus
the gradient descent iteration can be written

$w_{k+1}=w_{k}-\alpha_{k}X^{T}\nabla g(Xw_{k}).$

If we use a fixed step size, then for large $n$ and $d$ the _cost of this
update is dominated by the 2 matrix multiplications with $X$ required for each
iteration_. But instead of using a a fixed step size, consider the line
optimization (LO) problem of finding the $\alpha_{k}$ that maximally decreases
$f$,

$\displaystyle\mathop{\rm arg\,min}_{\alpha}f(w_{k}-\alpha X^{T}\nabla
g(Xw_{k})))$ (4) $\displaystyle\equiv$ $\displaystyle\mathop{\rm
arg\,min}_{\alpha}g(X(w_{k}-\alpha X^{T}\nabla g(Xw_{k})))$
$\displaystyle\equiv$ $\displaystyle\mathop{\rm
arg\,min}_{\alpha}g(\underbrace{Xw_{k}}_{m_{k}}-\alpha\underbrace{X(X^{T}\nabla
g(m_{k}))}_{d_{k}})$ $\displaystyle\equiv$ $\displaystyle\mathop{\rm
arg\,min}_{\alpha}g(m_{k}-\alpha d_{k}),$

where we have defined $m_{k}=Xw_{k}$ and $d_{k}=X(X^{T}\nabla g(m_{k}))$. The
key to efficiently optimizing over $\alpha$ is tracking the “memory”
$m_{k}\in\mathbb{R}^{n}$ across iterations. Notice that computing $d_{k}$
requires 2 matrix multiplications, that we have
$m_{k+1}=m_{k}-\alpha_{k}d_{k}$, and that given $m_{k}$ and $d_{k}$ we can
evaluate any value of $\alpha$ in $O(n)$ time plus the inexpensive cost of
evaluating $g$. This allows us to numerically _solve the LO problem with a
total of only 2 matrix multiplications per iteration_. Thus, in terms of the
bottleneck operation the cost of gradient descent with LO for LCPs is the same
as the cost of using a fixed step size. In the case of logistic regression,
the iteration cost with a fixed step size is $O(nd)$ while the iteration cost
if we use bisection to solve the LO to an accuracy of $\epsilon$ over a
bounded domain is only $O(nd+n\log(1/\epsilon))$. Thus, for large $d$ we can
perform LO to a high accuracy “for free”.

Note that the the value of $f$ achieved by LO may be substantially smaller
than the value achieved by using the strong Wolfe conditions. Indeed, for non-
quadratic functions the Armijo condition (3) used as part of the strong Wolfe
conditions may exclude the optimal step size (see Figure 1). Indeed, even for
convex problems the decrease in $f$ achieved by LO can be arbitrarily larger
than the largest decrease for any step size satisfying the Armijo condition.
Unfortunately, if $g$ is non-convex then the LO may be difficult to solve.
However, applying a local optimizer to the LO problem is still likely to
achieve a greater decrease in the function than line searches.

Figure 1: Valid step sizes allowed by Armijo condition (in blue) and update
minimizing function (in red) by line optimization. For a fixed $\sigma$, the
improvement of line optimization over the step sizes allowed by the Armijo
condition can be made arbitrarily large even for convex functions.

###  2.2 Efficient Plane Search (PS) for the Memory Gradient Method

The GD+M update (1) for LCPs has the form

$w_{k+1}=w_{k}-\alpha_{k}X^{T}\nabla g(Xw_{k})+\beta_{k}(w_{k}-w_{k-1}).$

With a fixed learning rate $\alpha_{k}$ and momentum rate $\beta_{k}$, notice
that the dominant cost of this update is again the 2 matrix multiplications
with $X$. Consider instead using a SO to set $\alpha_{k}$ and $\beta_{k}$ to
maximally decreases the function $f$,

$\displaystyle\mathop{\rm arg\,min}_{\alpha,\beta}f(w_{k}-\alpha X^{T}\nabla
g(Xw_{k})+\beta(w_{k}-w_{k-1}))$ (5) $\displaystyle\equiv$
$\displaystyle\mathop{\rm arg\,min}_{\alpha,\beta}g(X(w_{k}-\alpha X^{T}\nabla
g(Xw_{k})+\beta(w_{k}-w_{k-1})))$ $\displaystyle\equiv$
$\displaystyle\mathop{\rm
arg\,min}_{\alpha,\beta}g(\underbrace{Xw_{k}}_{m_{k}}-\alpha\underbrace{X(X^{T}\nabla
g(m_{k}))}_{d_{k}}+\beta(m_{k}-m_{k-1}))$ $\displaystyle\equiv$
$\displaystyle\mathop{\rm
arg\,min}_{\alpha,\beta}g(\underbrace{(1+\beta_{k})m_{k}-\alpha d_{k}+\beta
m_{k-1}}_{\text{potential $m_{k+1}$}}).$

Using the GD+M update (1) while optimizing $\alpha_{k}$ and $\beta_{k}$ is
called the _memory gradient_ method (Miele and Cantrell, 1969). For LCPs
adding the momentum term requires us to store $m_{k-1}$, but still only
requires 2 matrix multiplications total per iteration (the same as using a
fixed learning and momentum rate). If we solve the two-dimensional SO to
accuracy $\epsilon$ using a cutting plane method the cost for logistic
regressions iterations with optimized $\alpha_{k}$ and $\beta_{k}$ values is
$O(nd+n\log(1/\epsilon))$. Thus, using PS to optimize the learning and
momentum rate does not increase the asymptotic cost compared to using LO to
optimize the learning rate. However, note that optimizing both rates “for
free” is somewhat specific to the use of momentum as the second direction; if
we use directions that are not in the span of the previous parameters and
gradients this would require additional matrix multiplications.

Note that the memory gradient method does not require $\alpha_{k}$ or
$\beta_{k}$ to be positive. Thus, the method could set $\beta_{k}$ close to 0
to “reset” the momentum or could even use negative momentum. A common
variation on LCPs is functions of the form $f(w)=g(Xw)+g_{0}(w)$ for a
regularization function $g_{0}$. The memory gradient remains efficient in this
setting in the typical case where evaluating $g_{0}$ is much less expensive
than matrix multiplications with $X$. A typical choice is L2-regularization,
$g_{0}(w)=(\lambda/2)\|w\|^{2}$ for a positive regularization parameter, which
only costs $O(d)$ to evaluate.

###  2.3 LO and SO for LCPs in Practice

LCPs are widely used in ML, but is not common to fit these models with LO or
SO. We believe that this is because it is not common knowledge that this leads
to impressive performance gains compared to standard approaches to set the
step size(s). To help give the reader context, in this section we present
empirical results comparing a variety of ways to set the step size(s) for
linear and logistic regression. In particular, our first set of experiments
compares the following methods:

  * •

GD(1/L): gradient descent (2) with a step size of $\alpha_{k}=1/L_{k}$, where
$L_{k}$ is an estimate of the maximum curvature of the function. We initialize
to $L_{0}=1$ and each subsequent iteration is initialized with the previous
estimate $L_{k}=L_{k-1}$. But on each iteration the value $L_{k}$ is doubled
until we satisfy the inequality $f(w_{k}-(1/L_{k})\nabla f(w_{k}))\leq
f(x_{k})-(1/2L_{k})\|\nabla f(w_{k})\|^{2}$. This backtracking approach
corresponds to using the Armijo condition (3) with $\sigma=1/2$. This approach
guarantees a similar theoretical convergence rate to using a fixed step size
of $\alpha_{k}=1/L$, where $L$ is the Lipschitz constant of the gradient (see
Beck and Teboulle, 2009). In particular, this approach is within a constant
factor of the worst-case optimal fixed step size for gradient descent for
minimizing convex functions.

  * •

GD(LS): gradient descent (2) using a line search to find a learning rate
satisfying the strong Wolfe conditions. We use a standard implementation of
the line search (see Nocedal and Wright, 2006, Algorithm 3.5). The Armijo
sufficient decrease parameter $\sigma$ in (3) is set to 0.0001 which allows
larger step sizes than the GD(1/L) method, while the curvature parameter is
set to 0.9 requiring that the step sizes are large enough to slightly decrease
the magnitude of the directional derivative along the line. We initialize the
line search with the previous step size $\alpha_{k-1}$, double the step size
during the initial “bracketing phase”, and then take the midpoint of the
current bracket during the “zoom” phase.

  * •

GD+M(LS): a non-linear conjugate gradient (CG) variation on the gradient
descent with momentum (1) update,

$w_{k+1}=w_{k}+\alpha_{k}d_{k},$

where

$d_{k}=-\nabla f(w_{k})+\eta_{k}(w_{k}-w_{k-1}).$

We set $\eta_{k}$ using the non-negative variant of the formula of Polak and
Ribiere (1969),

$\eta_{k}=\max\left\\{0,\frac{\nabla f(w_{k})^{T}(\nabla f(w_{k})-\nabla
f(w_{k-1})}{\|\nabla f(w_{k})\|^{2}}\right\\}.$

We set $\alpha_{k}$ using a line search to satisfy the strong Wolfe
conditions. We reset the momentum rate to $\eta_{k}=0$ on iterations where the
directional derivative $\nabla f(w_{k})^{T}d_{k}$ of the search direction
$d_{k}$ is not negative.

  * •

GD(LO): gradient descent (2) using LO (4) to set the learning rate on each
iteration. We discuss how we numerically solve the LO in Appendix A.

  * •

GD+M(LO): the non-linear CG method GD+M(LS) but using LO to set the step size
$\alpha_{k}$. Note that with exact computation this is equivalent to the
linear conjugate gradient method for the special case of strongly-convex
quadratic functions.

  * •

GD+M(SO): gradient descent with momentum (1) where the learning rate and
momentum rate are optimized on each iteration using SO (5). We again refer to
Appendix A for how we numerically solve the SO problems. Note that with exact
computation this method is also equivalent to the linear conjugate gradient
method for strongly-convex quadratic functions.

Figure 2: Performance of different gradient-based methods for fitting logistic
regression models. Each plot is a different dataset. The black line only
backtracks, the blue lines use a line search that can decrease or increase the
step size to satisfy the strong Wolfe conditions, the orange lines use LO, and
the magenta line uses SO. The GD methods use the gradient direction, the
GD+M(L*) methods use the gradient direction and momentum with the non-linear
conjugate gradient relationship between the parameters, and the GD+M(SO)
method optimizes the learning rate and momentum rate. We see that LS methods
tend to dominate the 1/L method, GD+M methods tend to dominate GD methods, LO
methods tend to dominate LS methods, and the best performance on every dataset
was achieved with SO.

We explored the performance of these methods for logistic regression on 16
datasets drawn from the libsvm (Chang and Lin, 2011) collection and the UCI ML
Repository (Dua and Graff, 2017). These datasets were selected in Shea and
Schmidt (2023) to represent a variety of optimization challenges. We did not
add a bias term or standardize the datasets, which makes these datasets
particularly-challenging for optimizers to fit (an exception is the “madelon”
dataset, where we standardized the features as they were poorly scaled and the
optimizers were not able to make any reasonable progress). Figure 2 shows the
sub-optimality against the number of iterations of the different methods on
these datasets.222We approximated the optimal solution by running an variation
on the non-monotonic Barzilai-Borwein method of Raydan (1997) for 5,000
iterations. While in our experiments we plot the performance in terms of the
number of iterations, we remind the reader that for LCPs the iteration costs
of all of these methods are dominated by the cost of the 2 matrix
multiplications required per iteration. Thus, as the size of the data grows
the runtimes of the iterations of all methods we compare will only differ by a
small constant factor. This makes the comparison fair since we are most
interested in performance on larger datasets.

We observed several trends across the datasets in this experiment:

  1. 1.

GD(1/L) was always the worst method, even though theoretically it is close to
optimal. The superiority of practical step sizes highlights the limitations of
selecting algorithms purely based on theoretical convergence rates.

  2. 2.

LO methods tended to outperform LS methods, indicating that a more precise
line search improves performance. The GD(LO) method outperformed GD(LS) on
most datasets while the GD+M(LO) method outperformed GD+M(LS) on all datasets.

  3. 3.

Adding momentum tended to help, despite momentum not improving the theoretical
rate for this class of problems, and momentum helped more when using LO.
GD+M(LS) outerperformed GD(LS) on most problems, while after the first
iteration GD+M(LO) dominated GD(LO) across all iterations across all problems.

  4. 4.

The memory gradient method GD+M(SO), that optimizes both the learning rate and
momentum rate, performed as well or better than all other methods across every
iteration on every dataset.

To investigate why LO and SO improve performance, we plotted the step sizes
used by the different methods. In Figure 3 we show the step sizes used by the
different methods that only use the gradient direction. In these plots we see
that the theoretically-motivated 1/L method tended to use much smaller step
sizes than the other methods, either always accepting the initial step size or
decreasing the step size within the first few iterations and then using a
small constant step size. In contrast, the LO step sizes tended to oscillate
(with a period of 2) while often slowly increasing. Only two exceptions to
this pattern were observed, and on these problems the LO method quickly solved
the problem using large step sizes. The LS method tended to start with small
step sizes before growing into a similar range to the LO method’s
oscillations, and tended to also have oscillations but with a longer period.

Figure 3: Step sizes of different gradient descent methods for fitting
logistic regression models. Note that the GD(LO) step sizes tended to lead to
the fastest convergence while the GD(1/L) step sizes always converged slowest.

In Figure 4, we plot the learning rate and momentum rate for the methods that
incorporate momentum. Note that the optimal learning and momentum rates found
by the GD+M(SO) method tended not be constant, often differing by orders of
magnitude between iterations. The optimal learning rate and momentum rate do
seem to show oscillations, but the oscillations seem less regular than the
GD(LO) method where momentum is not used. We see that the non-linear conjugate
gradient method with a precise line optimization GD+M(LO) often closely tracks
the learning rate and momentum rates found by the GD+M(SO) method, while the
imprecise line search GD+M(LS) finds very different values (and indeed often
resets the momentum to zero on early iterations). We finally note that on one
iteration, the optimal momentum rate used by the GD+M(SO) method was negative
(absolute value of the rate is shown).

Figure 4: Step sizes of different gradient plus momentum methods for fitting
logistic regression models. The solid lines are the learning rates and the
dashed lines are the momentum rates. Note that the GD+M(SO) step sizes tended
to lead to the fastest convergence while the GD+M(LS) step sizes usually
converged slowest. A star is used to indicate the iteration where the momentum
rate was negative.

##  3 The Scattered 50+ Year History of Subspace Optimization Methods

The use of SO in both theoretical and practical works has a long history.
Indeed, the idea of SO arises in a variety of contexts in the numerical
optimization literature. In this section we review a variety of the most
closely related works.

Memory Gradient Method: the memory gradient method GD+M(SO) where SO is used
to set the learning rate and momentum rate in the GD+M update (1) was first
explored by Miele and Cantrell (1969). They proposed a particular method to
approximately solve the SO problem and presented numerical results showing
that SO could improve performance on a simple test function. However, Miele
and Cantrell’s work had relatively little impact. This is likely because
numerical optimization work in the 1960s did not have the modern focus on
exploiting problem structure (like LCPs) to quickly solve the SO.

Conjugate Gradient Methods: the conjugate gradient (CG) method was proposed in
1952 for solving positive-definite linear systems (Hestenes and Stiefel,
1952), but it can equivalently be written as using the memory-gradient method
to minimize strongly-convex quadratic functions (Cantrell, 1969). Due to its
outstanding empirical and theoretical properties, various non-linear CG
methods have been proposed for general smooth optimization. The first non-
linear CG method is due to Fletcher and Reeves (1964), but numerous variants
have appeared over the years (see Hager and Zhang, 2006). Non-linear CG
methods typically use a LS along a direction that is a weighted combination of
the gradient and momentum directions, and have been explored in some neural
network settings (LeCun et al., 2002). However, by fixing the weighting
between the learning rate and momentum rate, CG methods may make less progress
per iteration than the memory gradient method. Current non-asymptotic analyses
of non-linear CG methods show that their worst case performance is not better
than gradient descent with a line search (Gupta et al., 2023). Further, it was
already identified by Fletcher and Reeves that on some problems we need to
“reset” the momentum rate to 0 for non-linear CG to obtain good practical
performance. It is known that using LO and resetting non-linear CG every $n$
iterations leads to asymptotic quadratic convergence (as in Newton’s method)
in terms of $n$-iteration cycles (Cohen, 1972). The fact that non-linear CG’s
performance is improved by this resetting mechanism suggests that we might
improve performance by using the memory gradient method’s approach of
optimizing the momentum rate on each iteration (where some iterations could
use small or zero or even negative momentum rates).

Accelerated Gradient Methods: SO was important in the development of first-
order methods achieving the accelerated $O(1/k^{2})$ error rate for minimizing
convex functions (Nemirovski, 1982; Nemirovski and Yudin, 1983). Nemirovski
showed the $O(1/k^{2})$ rate could be achieved using SO over a subspace
containing the current gradient and specific weighted combinations of previous
gradients. Nemirovski later showed that we could achieve this rate with a LS
instead of SO, while Nesterov (1983) eventually showed the rate could be
achieved with a fixed learning rate and a modified momentum. Nesterov’s
accelerated gradient (NAG) algorithm can also be used to achieve accelerated
convergence rates for strongly-convex (Nesterov, 2013) and non-convex (Carmon
et al., 2017) objectives.

However, convergence rates only describe the worst case performance and it has
been highlighted that NAG often underperforms non-linear CG methods in
practice (Narkiss and Zibulevsky, 2005a). In Figure 5, we add a NAG variant to
the previous experiment. In particular, the NAG(1/L) method implements the
backtracking version of NAG described by Beck and Teboulle (2009). In this
experiment we see that the accelerated NAG(1/L) dominates the non-accelerated
GD(1/L), but that NAG(1/L) is consistently outperformed by the non-accelerated
non-linear CG method with LO GD+M(LO) and the non-accelerated memory gradient
method GD+M(SO).

Figure 5: Performance of different gradient-based methods for fitting logistic
regression models. The black lines only backtrack, the orange line uses LO,
the magenta line uses a two-dimensional SO, the light green line uses a three-
dimensional SO, and the dark green line uses a four-dimensional SO. The
NAG(1/L) is an implementation of Nesterov’s accelerated gradient method. The
NAG(SO) adds “gradient momentum” to the memory gradient method GD+M(SO), while
the SNAG(SO) method further adds scaling of the parameter vector. We see that
acceleration on its own tends to be less effective than using LO with an
appropriate direction. We also see that adding additional directions improves
performance but that only small gains tend to be observed by optimizing over
more than 2 directions.

Sequential Subspace Optimization: Narkiss and Zibulevsky (2005a) proposed
sequential subspace optimization (SESOP), which uses SO over multiple
directions in order to obtain the $O(1/k^{2})$ rate for convex optimization
but also obtain CG-like practical performance. They highlight LCPs as a
problem class where SO over the gradient and momentum directions can be
performed without significantly increasing the cost of the method. SESOP has
been applied to various LCP problems (Narkiss and Zibulevsky, 2005a, b; Elad
et al., 2007), and has been combined with diagonal scaling and a Hessian-free
Newton approach (Zibulevsky, 2008).

Supermemory Method: the supermemory gradient method (Cragg and Levy, 1969)
uses the most-recent momentum direction $(w_{k}-w_{k-1})$ and momentum
directions $(w_{t}-w_{t-1})$ from previous iterations $t$. A variation on this
idea is to use several previous gradient directions, and examples of this type
of method include the methods of Elad et al. (2007) and Wang and Yuan (2006).
It is known that using SO over all previous gradient directions has appealing
convergence properties (Drori and Taylor, 2020), but this can impose
significant additional computational costs. We note that for LCPs these
supermemory methods also only requires 2 matrix multiplications per iteration.
However, in our experiments we did not find that additional momentum terms
improved performance enough to justify the increased memory requirements and
increased SO cost.

Multi-Linear Problems: The linear structure in LCPs is what allows us to use
LO and SO within the gradient and memory gradient methods without changing
their iteration cost. A related problem is multi-linear problems where the
function is linear in certain subsets of the variables. Sorber et al. (2016)
discuss performing LO and SO for a supermemory method in the special case of
tensor factorization models. In Appendix B we consider variations on the
memory gradient method for matrix factorization, while Appendix C discusses
using LO or SO for log-determinant problems if we use low-rank directions.
Unfortunately, it does not seem possible to implement LO and SO for generic
multi-linear models without changing the iteration cost. We also note that
typical neural networks are not multi-linear. In the case of Sorber et al.
(2016), they only perform an SO every 4th iteration to keep costs low in their
application.

Acceleration and SO: NAG can be written as augmenting the GD+M method with a
third direction that is a form of “gradient momentum”, $(\nabla
f(w_{k})-\nabla f(w_{k-1}))$ (Shi et al., 2021). A CG method based on these
three directions was proposed by Andrei (2014) that gave excellent numerical
performance. But we could also consider a version of NAG that uses these 3
directions,

$x_{k+1}=x_{k}-\alpha_{k}\nabla
f(w_{k})+\beta_{k}(w_{k}-w_{k-1})+\gamma_{k}(\nabla f(w_{k})-\nabla
f(w_{k-1})),$ (6)

and performs 3-dimensional SO (this still only requires 2 matrix
multiplications per iteration for LCPs). This is appealing since different
choices of $\\{\alpha_{k},\beta_{k},\gamma_{k}\\}$ allow NAG to achieve
accelerated convergence rates in convex, strongly-convex, and non-convex
settings (Nesterov, 2013; Li and Lin, 2022). We thus might expect in some
cases that SO would adapt to the problem difficulty. Indeed, several authors
have shown that SO allows a single algorithm to achieve optimal rates in
several settings (Drori and Taylor, 2020; Nesterov et al., 2021; Guminov et
al., 2023).

In Figure 5 we consider using (6) with a three-dimensional subspace
optimization, under the name NAG(SO). In this experiment we found that NAG(SO)
typically only gave a small improvement over the memory gradient method. This
lack of improvement is consistent with the experiments of Karimi and Vavasis
(2017). They propose a restarted variant of the memory gradient method called
CGSO that can add a correction step to achieve the accelerated rate in the
strongly-convex setting. This correction step adds directions to the memory
gradient method, but they found that in practice the method applied the memory
gradient step in almost every case.

Scaled Gradient and Memory Methods: Sorber et al. (2016) consider a scaled
gradient method,

$x_{k+1}=\delta_{k}x_{k}-\alpha_{k}\nabla f(w_{k}),$

where on each iteration we estimate a scaling factor $\delta_{k}$ for the
original variable. For LCPs we can search for an optimal $\alpha_{k}$ and
$\delta_{k}$ while still only requiring 2 matrix multiplications per
iteration. However, Sorber et al. (2016) found that adding the scaling factor
$\delta_{k}$ did not significantly improve performance in their tensor
factorization experiments. In preliminary experiments we found that this
scaled gradient method improved performance over the gradient method with LO,
but did not tend to outperform the memory gradient method which also uses 2
directions.

We could also use a scaled variant of the memory gradient method

$x_{k+1}=\delta_{k}x_{k}-\alpha_{k}\nabla f(w_{k})+\beta_{k}(w_{k}-w_{k-1}),$
(7)

and solving for the three scaling factors still only requires 2 matrix
multiplications. In preliminary experiments we found that this 3-direction
approach can improve the performance of the memory gradient method, but tended
to be outperformed by the 3-direction NAG(SO) method that has momentum and
gradient momentum (6). On the other hand, a further small performance gain
could be obtained by combining all 4 directions into the update

$x_{k+1}=\delta_{k}x_{k}-\alpha_{k}\nabla
f(w_{k})+\beta_{k}(w_{k}-w_{k-1})+\gamma_{k}(\nabla f(w_{k})-\nabla
f(w_{k-1})).$

We refer to this method as SNAG(SO) in Figure 5.

Second-Order SO: Using LO to set the step size in Newton’s method does not
increase the iteration cost for many problem structures and preserves the
superlinear convergence rate of Newton’s method (Shea and Schmidt, 2023).
Recent work gives explicit non-asymptotic convergence rates of BFGS with LO
(Jin et al., 2024). Conn et al. (1994) discuss using SO with the gradient
direction and the Newton direction. This method improves on the global
convergence properties of Newton’s method while maintaining its superlinear
convergence rate (Zibulevsky, 2013; Shea and Schmidt, 2023). We note that it
is common to use SO to minimize quadratic approximations of the function, as
in trust region methods (see Nocedal and Wright, 2006), but in this work we
focus on methods that use SO to directly minimize the function. An approximate
second-order method has also been proposed in a majorization-minorization
framework for certain problem structures (Chouzenoux et al., 2010). However,
second-order information typically significantly increases the iteration cost.
Thus, in Section 5.1 we consider using SO within quasi-Newton methods that
cheaply approximate second-order information.

Global Optimization with SO: Wang et al. (2016) consider global optimization,
and search for an approximate solution by optimizing over a random subspace.
Recent works have explored iterative variations where a different subspace is
used on each iteration (Cartis and Otemissov, 2022; Cartis et al., 2023a) or
growing subspaces are used (Cartis et al., 2023b). But in this work we
restrict attention to the typical ML setting where the dimension is high and
only a local minimum may be found.

SGD with SO: Richardson et al. (2016) consider augmenting an SGD method with
SO. In particular, they alternate between taking a large number of SGD steps
with a small mini-batch and then taking an SO step with a large mini-batch.
The SO step uses the current gradient and differences between previous
iterates (as in the supermemory method). This method improved the performance
of three different SGD methods (SGD with momentum, NAG, and AdaGrad) in terms
of runtime on several deep learning tasks. In this work we restrict focus to
deterministic methods since these are the key ingredient behind the most
effective SGD methods. In particular, SGD with momentum is a stochastic
version of Polyak’s deterministic GD+M method while it has recently been shown
that the performance advantage of Adam on language models is due to its
deterministic properties (Kunstner et al., 2022). We expect that by improving
the deterministic ingredients of stochastic methods that this work will lead
to faster stochastic methods in the future.

##  4 Line Search and [Hyper-]Plane Search for Neural Networks

While SO has a long history and numerous applications, the work of Richardson
et al. (2016) is the only previous work we are aware of that considers using
SO to train neural networks. The lack of works using SO to train neural
networks is sensible because for general neural networks we cannot evaluate
the objective at many points in a subspace for a lower cost than evaluating
the objective at arbitrary points. Nevertheless, there exist neural networks
that are used in practice where a linear operator with the parameters is the
bottleneck operation. These networks allow us to efficiently implement SO to
set the learning rate and momentum rate, and indeed allow us to optimize over
separate rates for each layer. In these settings we may expect to see the same
performance gains that SO offers for LCPs.

###  4.1 Definition: SO-Friendly Neural Networks

We can write typical neural network objective functions in the form
$f(W,v)=g(h(XW)v)$ for non-linear functions $g$ and $h$. Here, $X$ is the data
matrix (in $\mathbb{R}^{n\times d}$), $W$ is the set of weights in the first
hidden layer (in $\mathbb{R}^{d\times r}$ with $r$ hidden units), and $v$ is
the concatenation of the parameters in the second and higher layers of the
network. We say that a neural network is SO-friendly if the cost of evaluating
$f$ is dominated by the cost of matrix multiplications with $X$. The dominant
computational cost in SO-friendly networks is similar to LCPs, and this allows
us to use LO and SO efficiently even though the parameters $v$ apply to a non-
linear function of the input.

There are a variety of possible SO-friendly neural network structures, but we
will first discuss the simple case of a network with two layers of hidden
weights. This type of network has been popularized in practice under the name
extreme learning machine (ELM) (Huang et al., 2006), and ELMs are among the
best performing out-of-the-box classifiers (Fernández-Delgado et al., 2014).
However, in ELMs the first layer of weights is set randomly while we consider
using SO to fit both layers efficiently. Two-layer networks are also popular
in reinforcement learning within the trust-region policy optimization and
proximal policy optimization frameworks (Schulman et al., 2015, 2017), where
simple networks have been sufficient to solve a number of continuous control
tasks (Rajeswaran et al., 2017). Finally, we note that having an efficient
method for fitting 2-layer networks would allow us to replace 1-layer networks
in settings where these degenerate linear models are used. For example, for
transfer learning we could train the last two layers of a pre-trained deep
network instead of only training the last layer.

###  4.2 Example: 2-Layer Networks with a Single Output - Tied Step Size(s)

A classic example of an SO-friendly neural network is a fully-connected neural
network with two layers of weights and a single output. In this setting $v$
has $r$ elements and $f$ has the simplified form $f(W,v)=g(h(XW)v)$, where the
activation function $h$ applies a non-linear operation such as a sigmoid
function element-wise and the loss function $g$ maps the $n$ predictions to a
scalar loss. The gradient descent update for a 2-layer neural network with a
single output takes the form

$\displaystyle W_{k+1}$
$\displaystyle=W_{k}-\alpha_{k}X^{T}\underbrace{\text{diag}(\nabla
g(h(XW_{k})v_{k}))h^{\prime}(XW_{k})\text{diag}(v_{k})}_{R_{k}},$
$\displaystyle v_{k+1}$
$\displaystyle=v_{k}-\alpha_{k}\underbrace{h(XW_{k})^{T}\nabla
g(h(XW_{k})v_{k})}_{\nabla_{v}f(W_{k},v_{k})},$

where $h^{\prime}$ is the matrix containing the derivative with respect to
each input and “diag” makes a diagonal matrix from a vector. With a fixed step
size the computational cost of this update is dominated by the two
multiplications with $X$ per iteration: one to compute $XW_{k}$ and another to
compute $X^{T}R_{k}$. Consider now performing LO to set the step size,

$\displaystyle\mathop{\rm arg\,min}_{\alpha}f(W_{k}-\alpha
X^{T}R_{k},v_{k}-\alpha\nabla_{v}f(W_{k},v_{k}))$ $\displaystyle\equiv$
$\displaystyle\mathop{\rm arg\,min}_{\alpha}g(h(X(W_{k}-\alpha
X^{T}R_{k})(v_{k}-\alpha\nabla_{v}f(W_{k},v_{k}))))$ $\displaystyle\equiv$
$\displaystyle\mathop{\rm
arg\,min}_{\alpha}g(h(\underbrace{XW_{k}}_{M_{k}}-\alpha\underbrace{X(X^{T}R_{k})}_{D_{k}})(v_{k}-\alpha_{k}\nabla_{v}f(W_{k},v_{k})))$
$\displaystyle\equiv$ $\displaystyle\mathop{\rm
arg\,min}_{\alpha}g(h(\underbrace{M_{k}-\alpha D_{k}}_{\text{potential
$M_{k+1}$}})(v_{k}-\alpha\nabla_{v}f(W_{k},v_{k}))).$

Analogous to LCPs, if we track the $n\times r$ product $M_{k}=XW_{k}$ then we
can perform LO using only two of the bottleneck matrix multiplications per
iteration: one to compute the $d$ by $r$ product $X^{T}R_{k}$, and another to
pre-mutiply this matrix by $X$. These two matrix multiplications cost
$O(ndr)$, and given $M_{k}$ and $D_{k}$ evaluating the LO objective only costs
$O(nr)$. Similar to the LCP case, we can also optimize the momentum rate
without additional matrix multiplications with $X$.

For neural networks the LO and SO problems are non-convex, and we did find
that it was possible to find poor local optima of the SO problem in
particular. Thus, care should be taken in how we try to solve the SO problem
and we should at least use a method that guarantees we decrease $f$ compared
to using a step size (or step sizes) of zero (otherwise LO/SO can result in
worse performance). Nevertheless, we found that the simple generic method used
in our LCP experiments (see Appendix A) typically performed well and we used
this method in all our neural network experiments.

###  4.3 LO and SO for 2-Layer Networks in Practice (Tied Step Sizes)

Figure 6: Performance of different gradient-based methods for 2-layer neural
networks. The black line only backtracks, the blue lines use a line search
that can decrease or increase the step size to satisfy the strong Wolfe
conditions, the orange lines use LO, and the magenta line uses SO. The GD
methods use the gradient direction, the GD+M(L*) methods use the gradient
direction and momentum with the non-linear conjugate gradient relationship
between the parameters, and the GD+M(SO) method optimizes the learning rate
and momentum rate. We see that the 1/L method has better performance for
training neural networks than for training linear models, but that the
GD+M(LO) and GD+M(SO) methods that exploit momentum and optimized step sizes
tend to outperform the other methods..

In this section we consider training 2-layer networks using the methods and
datasets from Section 2.3. We emphasize that the iteration costs of all
methods in this setting are dominated by the cost of the two products with $X$
on each iteration, and thus on large datasets all methods will have similar
runtimes. For these networks we use $r=100$ hidden units, use $h=\tanh$ as the
activation function, use the squared loss objective
$g(h(XW)v)=\|h(XW)v-y\|_{2}^{2}$, and initialize all elements of $W$ and $v$
with a sample from a standard normal distribution divided by the total number
of weights in the network $r(d+1)$.

We plot the results of the 2-layer training experiment in Figure 6. We
highlight two observations from this experiment:

  1. 1.

The essentially-fixed step size method GD(1/L) performed much better for
neural networks than for it did for the linear model. In particular, GD(1/L)
was often not substantially worse than the line search GD(LS) method and in
some cases performed better than the line search method. This counter-
intuitive behaviour is perhaps explained by the edge of stability phenomenon
(Cohen et al., 2020); unlike linear models where a constant step size must
adapt to the maximum curvature of the function, in neural networks it appears
empirically that the curvature of the function can adapt to the step size. In
Figure 7 we plot the step sizes of the different gradient descent methods, and
we see that the GD(1/L) method often found steps on a similar scale to the
GD(LS) method (unlike for linear models where the GD(1/L) step sizes were
typically too small). That the curvature of the network may be adapting to the
step size may explain why classic line searches have had little impact on the
training of neural networks.

  2. 2.

The methods using LO and SO with momentum, GD+M(LO) and GD+M(SO), consistently
outperformed the other methods. Thus, it appears that neural network training
can be sped up by per-iteration step size tuning but that it must be done
precisely and must use an appropriate per-iteration momentum. We did not find
obvious patterns in the step sizes used by GD+M(LO) and GD+M(SO) that were not
seen in the LCP experiments, although we note that the GD+M(SO) method used
negative momentum rates slightly more frequently than for linear models.

Figure 7: Step sizes of different gradient descent methods for fitting 2-layer
neural networks. The basic GD(1/L) tended to find step sizes closer to the
line search and LO methods on these problems.

###  4.4 Example: 2-Layer Networks with a Single Output - Per-Layer Step Sizes

While LO and SO appear to be helpful for training neural networks, we can
further take advantage of the layer structure to efficiently use different
step sizes in different layers. First consider the case of a 2-layer network
and a variant of gradient descent where we have a step size for each layer:

$\displaystyle W_{k+1}$ $\displaystyle=W_{k}-\alpha_{k}^{1}X^{T}R^{k},$
$\displaystyle v_{k+1}$
$\displaystyle=v_{k}-\alpha_{k}^{2}\nabla_{v}f(W_{k},v_{k}),$

where $\alpha_{k}^{1}$ is the step size for the first layer and
$\alpha_{k}^{2}$ is the step size for the second layer. We might expect better
performance by updating the different layers at different rates. Observe that
using SO to set both step sizes requires only the same two matrix
multiplications with $X$ that we need to set a single step size,

$\mathop{\rm
arg\,min}_{\alpha^{1},\alpha^{2}}g(h(\underbrace{M_{k}-\alpha^{1}D_{k}}_{\text{potential
$M_{k+1}$}})(v_{k}-\alpha^{2}\nabla_{v}f(W_{k},v_{k}))).$

The GD+M update (1) with a separate learning and momentum rate for each layer
has the form

$\displaystyle W_{k+1}$
$\displaystyle=W_{k}-\alpha_{k}^{1}X^{T}R_{k}+\beta_{k}^{1}(W_{k}-W_{k-1})$
$\displaystyle v_{k+1}$
$\displaystyle=v_{k}-\alpha_{k}^{2}\nabla_{v}f(W_{k},v_{k})+\beta_{k}^{2}(v_{k}-v_{k-1}),$

and optimizing over all 4 step sizes again only requires the same two matrix
multiplications,

$\displaystyle\mathop{\rm
arg\,min}_{\alpha^{1},\alpha^{2},\beta^{1},\beta^{2}}f(W_{k}-\alpha^{1}X^{T}R_{k}+\beta^{1}(W_{k}-W_{k-1}),v_{k}-\alpha^{2}\nabla_{v}f(W_{k},v_{k})+\beta^{2}(v_{k}-v_{k-1}))$
$\displaystyle\equiv$ $\displaystyle\mathop{\rm
arg\,min}_{\alpha^{1},\alpha^{2},\beta^{1},\beta^{2}}g(h(X(W_{k}-\alpha^{1}X^{T}R_{k}+\beta^{1}(W_{k}-W_{k-1}))(v_{k}-\alpha^{2}\nabla_{v}f(W_{k},v_{k})+\beta^{2}(v_{k}-v_{k-1}))))$
$\displaystyle\equiv$ $\displaystyle\mathop{\rm
arg\,min}_{\alpha^{1},\alpha^{2},\beta^{1},\beta^{2}}g(h(\underbrace{(1+\beta^{1})M_{k}-\alpha^{1}D_{k}-\beta^{1}M_{k-1}}_{\text{potential
$M_{k+1}$}})((1+\beta^{2})v_{k}-\alpha^{2}\nabla_{v}f(W_{k},v_{k})-\beta^{2}v_{k-1})).$

Thus, at each iteration we can numerically search for a learning rate for each
layer and a momentum rate for each layer for the same asymptotic cost as using
fixed step sizes. We could also consider adding previous iterations or
gradients as additional directions as in supermemory methods, or use a scaling
of each layer as in the scaled memory gradient method (7).

###  4.5 LO and SO for 2-Layer Networks in Practice (Per-Layer Step Sizes)

In this section, we repeat our previous experiment but include methods with
per-layer step sizes:

  * •

GD(SB): Gradient descent where SO is used to set a learning rate for each
layer (2 step sizes for 2-layer networks).

  * •

GD+M(SB): Gradient descent with momentum using the non-linear CG direction,
but where SO is used to set a step size for each layer (2 step sizes for
2-layer networks). We use the PR+ formula separately for each layer, but only
reset if the combined direction across layers is not a descent direction.

  * •

GD+M(SO+SB): Gradient descent with momentum where SO is used to set a separate
learning rate and momentum rate for each layer (4 step sizes for 2-layer
networks).

We once again repeat that the iteration costs of all methods are dominated by
the cost of the two products with $X$ on each iteration, although we note that
a 4-dimensional SO problem may take more time to solve than a 2-dimensional or
1-dimensional problem.

Figure 8: Performance of different gradient-based methods for 2-layer neural
networks. The red GD(SB) and GD+M(SB) use SO to a step size for each layer
while the purple GD+M(SO+SB) method set a learning and momentum rate for each
layer. We see that in some cases per-layer step sizes substantially improve
performance while in some cases per-layer step sizes harm performance.

Figure 8 shows the performance of these methods with per-layer step sizes and
the corresponding methods with step sizes tied across layers. In these plots
we observed two conflicting behaviours:

  * •

On some datasets using per-layer step sizes significantly improves
performance.

  * •

On some datasets using per-layer step sizes significantly harms performance.

Unfortunately, based on this performance it does not appear that per-layer
step sizes are a reliable option if we want a method that is free of hyper-
parameters. A potential source of the poor performance of per-layer step sizes
is that the optimized per-layer step sizes are often very large. We compare
the learning rates for gradient descent with optimized tied and per-layer
learning rates in Figure 9, where we see that the per-layer rates are often
much larger for at least one of the layers and have larger oscaillations
(sometimes oscillating between positive and negative step sizes). These large
step sizes lead to a large initial reduction in the objective, but they make
the norm of the parameters much larger than with tied step sizes. This large
norm seems to be the cause of the slow progress on later iterations for some
datasets.

Figure 9: Step sizes of gradient descent methods for fitting 2-layer neural
networks, optimizing a single step size for both layers (GD(LO)) or optimizing
a step size for each layer (GD(SB)). For the per-layer method, the solid line
is for the output layer and the dashed line is for the input layer. We plot
the absolute value, and include markers to indicate iterations where the step
was negative. We see that both methods tend to have cyclic behaviour, we see
that the per-layer method tends to use larger step sizes (for for at least one
of the layers) and tends to have larger oscillations. We also see that the
output layer often uses negative step sizes, in some cases using a negative
step size on every iteration but in many cases using a negative step size on
every second iteration.

To explore whether controlling the norm of the parameters would make per-layer
step sizes viable, we repeated this experiment using L2-regularization of the
parameters. In particular, we add to the objective a term of the form
$(\lambda/2)(\|W\|_{F}^{2}+\|v\|^{2})$ where we set $\lambda=1/n$. Figure 10
shows the performance of the same methods on the regularized objective. In
this plot, we see that per-layer step sizes typically lead to a performance
improvement (and in many cases a large improvement). We plot the four learning
and momentum rates used by best-performing method (GD+M(SO+SB)) in Figure 11,
showing the wide range of per-layer learning and momentum rates that lead to
this strong performance. It is possible that other strategies for controlling
the scale of the parameters, such as trust region methods, would also allow
per-layer step sizes to be effective if we do not want to use a regularized
objective.

Figure 10: Performance of different gradient-based methods for regularized
2-layer neural networks. The red GD(SB) and GD+M(SB) use SO to a step size for
each layer while the purple GD+M(SO+SB) method set a learning and momentum
rate for each layer. We see that per-layer step sizes consistently improve
performance.

Figure 11: Step sizes of gradient descent methods for fitting 2-layer neural
networks, optimizing a learning rate and momentum rate for each layer
(GD+M(SO+SB)). The solid lines are the learning rates and the dashed lines are
the momentum rates. We plot the absolute value, and include markers to
indicate iterations where the step was negative.

###  4.6 Other Examples of SO-Friendly Networks

Above we consider fully-connected 2-layer neural networks with a single
output. But there exists a variety of other settings where neural networks are
SO-friendly. Here we give three examples.

Fully-connected 2-layer networks with a small number of outputs: if we have a
2-layer network with $c$ outputs, the $r\times 1$ vector $v$ is replaced by a
$r\times c$ matrix $V$ and the the cost of forward and backward propagation
through the last layer of weights changes from $O(r)$ to $O(rc)$ for each of
the $n$ examples. In order to remain SO efficient, the $O(nrc)$ cost
associated with applying the last layer must be dominated by the $O(ndr)$ cost
of applying the first layer. Thus, 2-layer networks with multiple outputs are
SO-friendly provided that $c\ll d$, and we can use SO in cases where the
number of outputs is significantly smaller than the number of inputs. Note
that such situations arise in common ML datasets such as MNIST (LeCun et al.,
1998) which has 784 inputs and 10 outputs, CIFAR-10 (Krizhevsky et al., 2009)
which has1024 inputs and 10 outputs, and ImageNet (Deng et al., 2009) which
has 181,503 inputs for an average image and 1,000 outputs. Note that if we fix
the weights in the first layer of a 2-layer network, then the optimization
problem becomes an LCP. Thus, even if the number of labels is large it could
still make sense to do SO on only the second layer on a subset of the
iterations.

Wide-then-narrow deep networks: consider a fully-connected network with $\ell$
layers where the first layer has $r$ outputs and the remaining layers have at
most $\rho$ outputs. In this setting the cost of applying the network to all
examples is $O(ndr+nr\rho+n\rho^{2}\ell)$. Such networks are SO-friendly if we
have $\rho\ll d$ and $\rho\ell\ll r$. Such networks must have a large number
of inputs and can have a large number of hidden units in the first layer, but
must have a relatively smaller number of units in subsequent layers. Thus,
they can perform a large-dimensional linear transformation followed by a
number of lower-dimensional non-linear transforms. We give an example of a
5-layer network that is SO-friendly in Figure 12. While being SO-friendly
limits the possible structures available, note that this structure allows
efficient per-layer optimization of the learning and momentum rates for all
the layers.

Figure 12: Example structure of a 5-layer SO-friendly neural network
($d=10000,r=500,\rho=20,\ell=5$). For the same asymptotic cost as using a
fixed learning rate and momentum rate, we could numerically optimize a
learning rate and momentum rate for each of the 5 layers on each step (10 step
sizes).

Convolutional neural networks (CNNs) with large strides and/or pooling
operators: convolutional neural networks are a popular variant of neural
networks that include convolution and pooling layers (see LeCun et al., 1998;
Krizhevsky et al., 2012). Convolution layers may include a “stride” variable
that controls the size of the output of a layer compared to the input, and
pooling layers often also result in a dimensionality reduction. If a CNN is
setup with large-enough strides or pooling areas, it is possible to make SO-
friendly CNNs where applying the first layer is the dominant cost.

Non-linear neural tangent kernels: the NTK approach to training deep neural
network optimizes a linearized approximation to the network (Jacot et al.,
2018). This linearization allows us to efficiently use SO via the linear CG
method if we use the squared loss function. However, the linearization may
lose important information about the non-linearity present in the network. We
could alternately consider training with a linearization of only the first $l$
layers, and maintaining the non-linearity in the deep layers. If $l$ is chosen
large enough, the network will be SO-friendly (with the standard NTK
corresponding to setting $l$ to the number of layers). This would give a
better approximation than the NTK, and would allow optimizing shared step
sizes for the first $l$ layers and then separate step sizes for the layers
beyond layer $l$.

##  5 Augmenting Quasi-Newton and Adam with Subspace Optimization

The previous sections considered using SO to set step sizes within the GD+M
update. However, SO is also suitable for many of the common variants on this
update. In this section, we discuss using SO to improve on the performance of
two of the most empirically effective methods.

###  5.1 Augmenting Quasi-Newton with Subspace Optimization

Quasi-Newton updates can be written in the form

$w_{k+1}=w_{k}-\alpha_{k}B_{k}\nabla f(w_{k}),$ (8)

where the learning rate is chosen using the strong Wolfe conditions, and
$B_{k}$ is chosen to satisfy a variant of the secant equations, $B_{k}(\nabla
f(w_{k})-\nabla f(w_{k-1}))=(w_{k}-w_{k-1})$. One of the most successful
quasi-Newton methods is the L-BFGS method, which does not explicitly store
$B_{k}$ but implicitly forms the matrix based on the $p$ most-recent values of
$(w_{k}-w_{k-1})$ and $(\nabla f(w_{k})-\nabla f(w_{k-1}))$ (Nocedal, 1980).
The L-BFGS method only requires $O(dp)$ memory and allows multiplications with
$B_{k}$ to be performed in $O(dp)$, so with small $p$ the L-BFGS update has a
similar cost to the GD+M update. But in practice the approximate second-order
information in $B_{k}$ often allows the L-BFGS update to make substantially
more progress per iteration. For precise details on the method, see Nocedal
and Wright (2006, Chapters 6-7). The computational savings for LO in LCPs and
SO-friendly networks can also be used with updates of the form (8),

$\displaystyle\mathop{\rm arg\,min}_{\alpha}f(w_{k}-\alpha B_{k}X^{T}\nabla
g(Xw_{k}))$ $\displaystyle\equiv$ $\displaystyle\mathop{\rm
arg\,min}_{\alpha}g(X(w_{k}-\alpha B_{k}X^{T}\nabla g(Xw_{k})))$
$\displaystyle\equiv$ $\displaystyle\mathop{\rm
arg\,min}_{\alpha}g(m_{k}-\alpha\underbrace{X(B_{k}X^{T}\nabla
g(m_{k}))}_{d_{k}})$ $\displaystyle\equiv$ $\displaystyle\mathop{\rm
arg\,min}_{\alpha}g(\underbrace{(1+\beta)m_{k}-\alpha d_{k})}_{\text{potential
$m_{k+1}$}}).$

We can consider adding directions to quasi-Newton methods and using SO to set
the corresponding step sizes. A natural update to consider is a hybrid with
the gradient method and quasi-Newton method

$w_{k+1}=w_{k}-\alpha_{k}^{1}\nabla f(w_{k})-\alpha_{k}^{2}B_{k}\nabla
f(w_{k}).$

This method would guarantee at least as much progress as gradient descent on
each iteration, which could lead to more progress than a quasi-Newton update
if the matrix $B_{k}$ is badly scaled. However, this update would require 3
matrix multiplications with $X$ for LCPs and SO-friendly neural networks.
Thus, on large problems adding the gradient direction would increase the
iteration cost by around 50% compared to the standard quasi-Newton method.

Alternately, we could consider adding a momentum term to the quasi-Newton
update,

$w_{k+1}=w_{k}-\alpha_{k}B_{k}\nabla f(w_{k})+\beta_{k}(w_{k}-w_{k-1}).$ (9)

The computational savings for SO in LCPs and SO-frinedly networks are
preserved for this “quasi-Newton with momentum” method. In the case of LCPs
solving for $\alpha_{k}$ and $\beta_{k}$ corresponds to solving

$\displaystyle\mathop{\rm arg\,min}_{\alpha,\beta}f(w_{k}-\alpha
B_{k}X^{T}\nabla g(Xw_{k})+\beta(w_{k}-w_{k-1}))$ $\displaystyle\equiv$
$\displaystyle\mathop{\rm arg\,min}_{\alpha,\beta}g(X(w_{k}-\alpha
B_{k}X^{T}\nabla g(Xw_{k})+\beta(w_{k}-w_{k-1})))$ $\displaystyle\equiv$
$\displaystyle\mathop{\rm
arg\,min}_{\alpha,\beta}g(m_{k}-\alpha\underbrace{X(B_{k}X^{T}\nabla
g(m_{k}))}_{d_{k}}+\beta(m_{k}-m_{k-1}))$ $\displaystyle\equiv$
$\displaystyle\mathop{\rm
arg\,min}_{\alpha,\beta}g(\underbrace{(1+\beta)m_{k}-\alpha d_{k}+\beta
m_{k-1}}_{\text{potential $m_{k+1}$}}),$

where $m_{k}=Xw_{k}$ and $d_{k}=X(B_{k}X^{T}\nabla g(m_{k}))$. Similar to the
GD+M update this only requires 2 matrix multiplications with $X$ per
iteration. Thus, for LCPs we can optimally set the learning rate of the L-BFGS
method and use a non-zero momentum term with an optimally set momentum rate.
Further, for SO-friendly networks we could additionally set these step sizes
on a per-layer basis.

In Figures 13 and 14 we plot the performance of methods that use L-BFGS
directions where the step size is set using the Wolfe conditions (QN(LS)), the
step size is set using LO (QN(LO)), or we add a momentum term and set both
step sizes using SO (QN+M(SO)). For all methods, we use the method of Shanno
and Phua (1978) to improve the scaling of $B_{k}$, and for the QN(LS) method
we initialize the line search with a step size of $\alpha_{k}=1$. While the
differences are smaller than what we see for gradient methods, in these
figures we see that using LO often improves performance and that using SO with
the additional momentum term often further improves performance.

Figure 13: Performance of different quasi-Newton methods for fitting logistic
regression models. The blue line uses a line search initialized with 1, the
orange line uses LO, and the magenta line uses SO to optimize the learning and
momentum rates.

Figure 14: Performance of different quasi-Newton methods for 2-layer neural
networks. The blue line uses a line search initialized with 1, the orange line
uses LO, and the magenta line uses SO to optimize the learning and momentum
rates

It may appear surprising that we can improve the performance of quasi-Newton
methods with LO and SO, since asymptotically if quasi-Newton methods converge
to a strict minimizer then the choice of $\alpha_{k}=1$ is optimal and
momentum is not needed. However, we note that for any finite iteration quasi-
Newton can benefit from using step sizes other than 1 and can benefit from
momentum. In Figures 15 and 16 we plot the step sizes. In these plots we that
the line search QN(LS) method typically accepts the initial step size of 1,
while the LO and SO methods typically use step sizes that are slightly larger
than 1. We also see that the momentum rate used by the QN+M(SO) method is
typically much smaller than the step size, and that negative momentum steps
are often used (we were surprised that $\alpha_{k}$ was also negative on one
iteration on the phishing dataset).

Figure 15: Step sizes of different quasi-Newton methods for fitting logistic
regression models. We plot the absolute value, and include markers to indicate
iterations where the step was negative. The dashed line is the momentum rate
for the SO method. We see that the line search typically uses a learning rate
of 1, while the LO and SO methods typically use slightly larger values. Note
that the momentum rate was frequently negative in the SO method.

Figure 16: Step sizes of different quasi-Newton methods for fitting = 2-layer
neural networks. We plot the absolute value, and include markers to indicate
iterations where the step was negative. The dashed line is the momentum rate
for the SO method. We see that the line search typically uses a learning rate
of 1, while the LO and SO methods typically use slightly larger values. Note
that the momentum rate was frequently negative in the SO method.

###  5.2 Augmenting Adam with Subspace Optimization

The Adam optimizer uses the update

$\displaystyle\mu_{k+1}$
$\displaystyle=\beta_{k}^{1}\mu_{k}+(1-\beta_{k}^{1})\nabla f(w_{k})$ (10)
$\displaystyle v_{k+1}$
$\displaystyle=\beta_{k}^{2}v_{k-1}+(1-\beta_{k}^{2})\nabla
f(w_{k})\circ\nabla f(w^{k})$ $\displaystyle w_{k+1}$
$\displaystyle=w_{k}-\alpha_{k}\underbrace{V_{k+1}^{-1}\mu_{k+1}}_{d_{k}},$

for a learning rate $\alpha_{k}$, gradient momentum rate $\beta_{k}^{1}$,
squared-gradient momentum rate $\beta_{k}^{2}$, and $V_{k}$ being a diagonal
matrix where diagonal $i$ given by $(\sqrt{(v_{k})_{i}}+\epsilon)$ for a small
$\epsilon$ (we use $\circ$ for element-wise multiplication).333We omit the
bias correction factors in the update of $w_{k}$ for simplicity and because
they made the performance of the method worse in our experiments. But these
factors do not affect the efficiency of SO. In many applications the three
learning rates $\\{\alpha_{k},\beta_{k}^{1},\beta_{k}^{2}\\}$ are set to fixed
values across iterations, with a default choice being $\alpha_{k}=0.001$,
$\beta_{k}^{1}=0.99$, $\beta_{k}^{2}=0.999$, and $\epsilon=10^{-8}$. However,
for over-parameterized problems recent works show that the out-of-the-box
performance of Adam is improved by using a line search to set $\alpha_{k}$ on
each iteration (Vaswani et al., 2020; Galli et al., 2023). For LCPs and SO-
friendly networks we can further improve the performance of Adam using LO. For
example, with 2 matrix multiplications per iteration we could use LO to set
$\alpha_{k}$: one to compute $\nabla f(w^{k})$ and one to compute
$XV_{k+1}^{-1}\mu_{k}$.

Unfortunately, even for LCPs it does not appear efficient to use SO to set the
learning and momentum rates in Adam. It would require 3 matrix multiplications
per iteration to optimize $\alpha_{k}$ and $\beta_{k}^{1}$: one to compute
$\nabla f(w^{k})$, one to compute $XV_{k+1}{-1}\mu_{k}$, and one to compute
$XV_{k+1}^{-1}\nabla f(w_{k})$. Unfortunately, this increase the iteration by
50% in terms of the bottleneck matrix multiplication operations. Further, it
does not appear that the structure of the algorithm allows to efficiently use
SOs to optimize all three parameters even for LCPs.

Nevertheless, for amenable problem structures it is possible to add additional
directions to Adam that can be efficiently optimized with SO. For example, we
explored adding an additional momentum term $\beta_{k}(w_{k}-w_{k-1})$ to Adam
and found this improved performance. But we found that a larger performance
increase could be obtained by using multiple Adam directions,

$w_{k+1}=w_{k}-\alpha_{k}^{1}d_{k}-\alpha_{k}^{2}d_{k-1},$ (11)

where $d_{k}$ is defined by the Adam algorithm (10) Despite Adam usually being
considered as a method that works with stochastic gradients, we found that
this update was often an effective deterministic optimization method when
using SO to set the step sizes.

We plot the performance of a variety of Adam variants for logistic regression
and 2-layer networks in Figures 17 and 18. The methods we compare are using a
step size of $\alpha=10^{-3}$ as in PyTorch (Adam(Default)), using the strong
Wolfe line search to set $\alpha_{k}$ initialized with the previous step size
(Adam(LS)), using LO to set $\alpha_{k}$ (Adam(LO)), or using SO for update
(11) to optimize over two Adam step sizes (Adam2(SO)). The other Adam hyper-
parameters are fixed at their default values. The Adam direction $d_{k}$ may
not be a descent direction, and on iterations where this happened we made the
Adam(LS) search along the direction $-d_{k}$ instead of $d_{k}$. For logistic
regression we see that the default step size tends to perform poorly, while
the more clever step sizes perform better and the multi-direction method with
SO tends to perform the best. For neural networks the performance differences
were not as clear; on some datasets the default sep size performed poorly
while on other datasets it was competitive. However, on neural networks the SO
method still tended to have the best performance.444An exception to the good
performance of Adam(SO) was on the splice dataset, where on the second
iteration the method converged numerically to a stationary point and halted.
We found that this degenerate behaviour was not observed if we added
regularization.

Figure 17: Performance of different variations on Adam for fitting logistic
regression models. We see that the default Adam settings (black line) tend to
perform poorly, while using a line search (blue) or line optimization (orange)
improves the performance. The multi-direction method using SO(magenta)
dominated the single-direction methods.

Figure 18: Performance of different variations on Adam for fitting 2-layer
neural networks. We see that the default Adam settings (black line) are more
competitive for neural networks than for linear models. The multi-direction
method using SO(magenta) often outperformed the single-direction methods.

The relative performance differences of the different Adam variants can
partially be explained by looking at the step size plots Figures 19 and 20.
For logistic regression we see that the default step is often too small while
for neural networks the default step size appears more reasonable. We also see
that the line search method often uses the same step size for many iterations
(and occasionally used a negative step size), while the LO method was
constantly changing the step size. Periodic behaviour is observed in the LO
method, but the period appears larger than it did for gradient decent with LO.
Finally, the SO method tended to use step sizes with similar magnitudes on
$d_{k}$ and $d_{k-1}$, but on many iterations it used a positive step size on
$d_{k}$ and a negative step size on $d_{k-1}$. Thus, the Adam(SO) method often
seemed to approximate the difference between successive Adam directions in
order to make more progress in decreasing the objective.

Figure 19: Step sizes of different variations on Adam for fitting logistic
regression models. We plot the absolute value, and include markers to indicate
iterations where the step was negative. The dashed line is for the step size
on $d_{k-1}$ in the SO method. We see that the default step size tends to be
too small, that the line search tends to use larger constant values over
several iterations, that the LO seems to show periodic behaviour with a period
longer than 2, and that the SO method typically shows periodic behaviour with
a period of 2 where positive and negative step sizes of similar magnitudes are
used.

Figure 20: Performance of different first-order methods for 2-layer networks.
We plot the absolute value, and include markers to indicate iterations where
the step was negative. The dashed line is for the step size on $d_{k-1}$ in
the SO method. We see that the default step size tends to be in a similar
range to the other methods, that the line search tends to use constant values
over several iterations, the LO is consistently changing the step size, and
that the SO method often shows periodic behaviour with a period of 2 where
positive and negative step sizes of similar magnitudes are used.

###  5.3 Gradient, Quasi-Newton, and Adam with[out] Subspace Optimization

Our experiments indicate that gradient descent with momentum benefit from LO
and SO, that quasi-Newton methods benefit from LO while SO allows us to add an
optimized momentum term, and that Adam often benefits from LO while using SO
to combine multiple Adam directions yields further improvements. However, up
to this point we have not compared the three styles of methods (gradient,
quasi-Newton, and Adam) to each other. In Figures 21 and 22 we compare the
performance of what we view as the best-performing method of each of the 3
types if we exploit LO/SO and the best performing methods if we do not exploit
the problem structure.

These “best methods” plots indicate that the choice regarding whether or not
to use SO seems more important than choosing between gradient/quasi-
Newton/Adam methods. Among the methods not exploiting SO, only the QN(LS)
method seems to consistently come close to the performance of the SO methods
(this could be due to its natural step size of $\alpha_{k}=1$). On the
logistic regression problems the QN+M(SO) method seems to be the most
effective, always performing at least as well or better than all other
methods. On the neural network problems the QN+M(SO) was also consistently
among the best methods, but on a small number of datasets the multi-direction
Adam2(SO) method performed the best.

Figure 21: Performance of different gradient-based, quasi-Newton, and Adam
methods for fitting logistic regression models. The blue lines use a classic
line search that does not exploit the problem structure while the purple lines
exploit the problem structure to perform SO. We see that the SO methods tend
to dominate, with only the QN(LS) method coming close to the performance of
the SO methods. Among the SO methods, the QN+M(SO) method has among the best
performance across datasets.

Figure 22: Performance of different gradient-based, quasi-Newton, and Adam
methods for fitting 2-layer neural networks. The blue lines use a classic line
search that does not exploit the problem structure while the purple lines
exploit the problem structure to perform SO. We see that SO methods tend to
perform better across most datasets. On most datasets the QN+M(SO) method has
among the best performance, except on some datasets while the Adam2(SO) method
performs better.

##  6 Beyond Deterministic, SO-Friendly, and Differentiable Problems

In this work we show that SO can be efficiently applied to a restricted class
of neural networks, and it seems to generically improve the performance of
deterministic optimization algorithms. However, many deep learning
architectures do not follow these structures and are trained with SGD
variants. In this section we discuss directions that future work could take to
exploit SO within stochastic algorithms and to use SO to speed up training
neural networks that are not SO-friendly.

###  6.1 Using Subspace Optimization within SGD

SGD methods use variants of the GD+M algorithm where $\nabla f(w_{k})$ is
computed based on a subset of the training examples, which can be advantageous
if the number of training examples is large. If we use LO or SO to set the SGD
step size(s), SGD does not converge in general. However, it has recently been
shown that some LS variants do converge for over-parameterized models and that
this gives impressive performance in various settings (Vaswani et al., 2019;
Galli et al., 2023). We believe that SO could be used to further improve the
performance of SGD for over-parameterized models. For under-parameterized
models SO could be used with growing batch sizes (Friedlander and Schmidt,
2012; Bollapragada et al., 2018; De et al., 2016), with alternating SGD steps
(Richardson et al., 2016), or to implement better models of the objective (Asi
and Duchi, 2019). We could alternately simply use a very-large fixed batch
size, such as the 3.2 million used by GPT-3 (Brown et al., 2020), since large
batch sizes makes SGD behave like a deterministic algorithm up to a fixed
accuracy.

###  6.2 Using Subspace Optimization for Deep Learning

Many deep learning models are not SO-friendly. Nevertheless, we may be able to
use SO to implement methods that make more progress per iteration than SGD.
For example, prox-linear methods partially linearize functions that are
written as compositions (such as neural networks) and can guarantee more
progress per iteration than gradient descent (Nesterov, 2007; Lewis and
Wright, 2016; Duchi and Ruan, 2018; Drusvyatskiy and Paquette, 2019). The
iterations of the prox-linear method are LCPs, and thus this is a natural
direction to explore to obtain better optimization algorithms for deep
learning.

An alternative way to use SO for deep learning is with layer-wise training. It
has been shown that layer-wise training allows building deep models that are
competitive on ImageNet (Belilovsky et al., 2019). Layer-wise training can be
formulated as an LCP and may be a way to remove the frustrations associated
with training deep networks (due to sensitivity to hyper-parameters). While
the default layer-wise method would focus on optimizing a single layer, it may
be possible to quickly optimize SO-friendly sub-networks of a deep network to
yield a fast and robust training method.

Acknowledgments and Disclosure of Funding

We would like to thank Michael Zibulevsky for pointing us towards additional
references on the topic. Betty Shea is funded by an NSERC Canada Graduate
Scholarship. The work was partially supported by the Canada CIFAR AI Chair
Program and NSERC Discovery Grant RGPIN-2022-036669.

##  Appendix A Solving SO Problems

Conn et al. (1994) discuss the details of a way to numerically solve the SO
problem at every iteration. They suggest using a full-memory BFGS method,
which finds a local optimum with a local superlinear convergence rate under
suitable conditions. Using the full-memory variant is feasible since the SO
problem is low-dimensional, and allows us to easily modify the BFGS matrix to
prevent it becoming unrelated to the gradient direction (which can cause BFGS
to not converge). Conn et al. (1994) suggest terminating the BFGS iterations
when a lack of potential progress is detected. Despite solving the SO problem
inexactly, they prove that the overall algorithm is convergent if the learning
rate $\alpha$ is initialized with a line search and other step sizes are
initialized to 0 (assuming that the BFGS iterations monotonically decrease
$f$).

In our experiments, we found that a simple approach to solving the SO problem
led to good performance in practice. In particular, we solved the SO using the
simple Barzilai-Borwein (BB) (Barzilai and Borwein, 1988) approximate quasi-
Newton method, augmented with a non-monotonic Armijo line-search (Grippo et
al., 1986) as in Raydan (1997). The BB method has a superlinear convergence
rate for two-dimensional quadratic functions, so we expect it to work
efficiently for low-dimensional problems. To initialize the method, we simply
set all step sizes to zero.

As opposed to quasi-Newton methods, SESOP uses a Newton-like method to
numerically solve the SO (Narkiss and Zibulevsky, 2005a). But this approach
requires extra implementation effort for new problems to compute the second
derivatives on the subspace. Alternately, cutting plane methods generalize the
bisection-style methods used in line-search codes, but this requires
additional memory and these methods require careful implementation in cases
where the SO is non-convex (Hinder, 2018). We note that the best performance
in terms of runtime might be obtained by increasing the accuracy that we solve
the SO problem as the algorithm runs; when far from a solution it is unlikely
that a high accuracy is needed but we can benefit from an accurate solution as
we approach an optimum. However, it is not obvious how to optimally choose the
accuracy at each step555One practical approach could be to vary the maximum
number of iterations allowed for the subproblem solver depending on the stage
of the main optimization problem. For example, if we are at the $k$th
iteration of the optimization process, we can set the maximum number of
iterations allowed in the inner loop maxIterInnerLoop to be some multiple of
$k$..

##  Appendix B Matrix Factorization Problems

Matrix factorization (MF) problems can be written in the form
$f(U,W)=g(UW^{T})$ for a function $g$ and matrices $U\in\mathbb{R}^{n\times
r}$ and $W\in\mathbb{R}^{d\times r}$ for some rank $r\leq\min\\{n,d\\}$. For
this problem setting, SO becomes appealing if the matrix multiplication
$UW^{T}$ is the dominant cost, which happens as the rank $r$ grows (in the
typical case where evaluating $g$ is less expensive than the matrix
multiplication). A classic example is the principal component analysis (PCA)
problem, which can be written $f(U,W)=\frac{1}{2}\|UW^{T}-X\|_{F}^{2}$ for the
data matrix $X$ where $\|\cdot\|_{F}$ is the matrix Frobenius norm. For the
PCA setting we have $g(M)=\frac{1}{2}\|M-X\|^{2}$, and thus evaluating $g$
becomes cheaper than computing the matrix multiplication $UW^{T}$ as the rank
$r$ grows.

The MF gradient matrices have the form $\nabla_{U}f(U,W)=\nabla g(UW^{T})W$
and $\nabla_{W}f(U,W)=\nabla g(UW^{T})^{T}U$. Thus, the GD+M update has the
form

$\displaystyle U_{k+1}$ $\displaystyle=U_{k}-\alpha_{k}^{1}\nabla
g(U_{k}W_{k}^{T})W_{k}+\beta_{k}^{1}(U_{k}-U_{k-1}),$ $\displaystyle W_{k+1}$
$\displaystyle=W_{k}-\alpha_{k}^{2}\nabla
g(U_{k}W_{k}^{T})^{T}U_{k}+\beta_{k}^{2}(W_{k}-W_{k-1}).$

We have written this update in terms of four step sizes: a learning rate
$\alpha_{k}^{1}$ on $U$, a learning rate $\alpha_{k}^{2}$ on $W$, and separate
a momentum rate $\beta_{k}^{1}$ and $\beta_{k}^{2}$ for each matrix. The
bottleneck in this update is matrix multiplications costing $O(ndr)$. With all
rates fixed to constant values this update requires 3 matrix multiplications
with this cost: $U_{k}W_{k}^{T}$ and the matrix $\nabla g(U_{k}W_{k}^{T})$
multiplied by both $W_{k}$ and $U_{k}$. A common variation is to update only
one of the two matrices on each iteration (“alternating minimization”) (Bell
and Koren, 2007), and in this case only 2 matrix multiplications are
required.666If we only update one matrix on each iteration, it can make sense
to replace $U_{k-1}$ and $W_{k-1}$ in the momentum term with the matrices from
earlier iterations.

There are several ways we could incorporate SO into a MF optimization method:

  1. 1.

Alternating Minimization: if we only update one matrix on each iteration, then
the problem considered on each iteration is an LCP. Thus, we can use SO to set
the learning rate and momentum rate with only 2 matrix multiplications.

  2. 2.

Simultaneous Gradient Descent: if we do not use momentum terms, we can use SO
to optimally set the two step sizes $\alpha_{k}^{1}$ and $\alpha_{k}^{2}$,

$\displaystyle\mathop{\rm
arg\,min}_{\alpha^{1},\alpha^{2}}f(U_{k}-\alpha^{1}\nabla_{U}f(U_{k},W_{k}),W_{k}-\alpha^{2}\nabla_{W}f(U_{k},W_{k}))$
$\displaystyle\equiv$ $\displaystyle\mathop{\rm
arg\,min}_{\alpha^{1},\alpha^{2}}g((U_{k}-\alpha^{1}\nabla
g(U_{k}W_{k}^{T})W_{k})(W_{k}-\alpha^{2}\nabla g(U_{k}W_{k})^{T}U_{k})^{T})$
$\displaystyle\equiv$ $\displaystyle\mathop{\rm
arg\,min}_{\alpha^{1},\alpha^{2}}g(U_{k}W_{k}^{T}-\alpha^{1}\underbrace{\nabla
g(M_{k})W_{k}W_{k}^{T}}_{D_{k}^{1}}-\alpha^{2}\underbrace{U_{k}U_{k}^{T}\nabla
g(M_{k})}_{D_{k}^{2}}+\alpha^{1}\alpha^{2}\underbrace{\nabla
g(M_{k})W_{k}U_{k}^{T}\nabla g(M_{k})}_{D_{k}^{3}})$ $\displaystyle\equiv$
$\displaystyle\mathop{\rm
arg\,min}_{\alpha^{1},\alpha^{2}}g(M_{k}-\alpha^{1}D_{k}^{1}-\alpha^{2}D_{k}^{2}+\alpha^{1}\alpha^{2}D_{k}^{3}).$

If we have stored $M_{k}$ from the previous iteration, this requires a total
of 5 matrix multiplications (two to form $D_{k}^{1}$, two to form $D_{k}^{2}$,
and one more to form $D_{k}^{3}$). This is 2 more multiplications than using
fixed a step size, but we would expect that using two step sizes would lead to
more progress per iteration. Further, since the asymptotic cost is increased
by a factor of less than 2, using SO to set both step sizes on each iteration
is faster than re-running with two different sets of fixed step sizes.
Further, compared to methods that only update one matrix per iteration,
updating both matrices with an optimal learning rate only requires one
additional matrix multiplication for both matrices to be updated. Finally, we
note that if only one learning rate was used ($\alpha^{1}=\alpha^{2}$), the LS
would still require up to 5 matrix multiplications. Thus, when updating both
matrices, SO is not “free” in the sense that it does require additional matrix
multiplications.

  3. 3.

Momentum on One Matrix: If we use momentum on only one of the matrices, it
requires a total of 7 matrix multiplications to set all three step sizes. If
we choose the $U$ matrix, this leads to an update of the form

$\displaystyle U_{k+1}$ $\displaystyle=U_{k}-\alpha^{1}\nabla
g(U_{k}W_{k}^{T})W_{k}+\beta(U_{k}-U_{k-1}),$ $\displaystyle W_{k+1}$
$\displaystyle=W_{k}-\alpha^{2}\nabla g(U_{k}W_{k}^{T})^{T}U_{k}).$

The function value for a given choice of $\\{\alpha^{1},\alpha^{2},\beta\\}$
is given by

$\displaystyle
f(U_{k}-\alpha^{1}\nabla_{U}f(U_{k},W_{k})+\beta(U_{k}-U_{k-1}),W_{k}-\alpha^{2}\nabla_{W}f(U_{k},W_{k}))$
$\displaystyle=$ $\displaystyle g(((1+\beta)U_{k}-\beta
U_{k-1}-\alpha^{1}\nabla g(U_{k}W_{k}^{T})W_{k})(W_{k}-\alpha^{2}\nabla
g(U_{k}W_{k})^{T}U_{k})^{T})$ $\displaystyle=$ $\displaystyle
g((1+\beta)U_{k}W_{k}^{T}-\alpha^{1}\nabla
g(U_{k}W_{k}^{T})W_{k}W_{k}^{T}-\alpha^{2}(1+\beta)U_{k}U_{k}^{T}\nabla
g(U_{k}W_{k}^{T})$ $\displaystyle+\alpha^{1}\alpha^{2}\nabla
g(U_{k}W_{k}^{T})W_{k}U_{k}^{T}\nabla
g(U_{k}W_{k}^{T}){\color[rgb]{1,0,0}-\beta U_{k-1}W_{k}^{T}+\alpha^{2}\beta
U_{k-1}U_{k}^{T}\nabla g(U_{k}W_{k}^{T}))}).$

These 7 matrix multiplications are two more than the 5 we require if we do not
include the momentum term (the extra matrix multiplications come from the
terms highlighted in red). If we used momentum on $W$ instead of $U$, it would
similarly require 7 matrix multiplications.

  4. 4.

Momentum on Both Matrices - Exact SO: the update if we use a momentum term on
both matrices is given by

$\displaystyle U_{k+1}$ $\displaystyle=U_{k}-\alpha^{1}\nabla
g(U_{k}W_{k}^{T})W_{k}+\beta^{1}(U_{k}-U_{k-1}),$ $\displaystyle W_{k+1}$
$\displaystyle=W_{k}-\alpha^{2}\nabla
g(U_{k}W_{k}^{T})^{T}U_{k}+\beta^{2}(W_{k}-W_{k-1})).$

The function value for a given choice
$\\{\alpha^{1},\alpha^{2},\beta^{1},\beta^{2}\\}$ is given by

$\displaystyle
f(U_{k}-\alpha^{1}\nabla_{U}f(U_{k},W_{k})+\beta^{1}(U_{k}-U_{k-1}),W_{k}-\alpha^{2}\nabla_{W}f(U_{k},W_{k})+\beta^{2}(W_{k}-W_{k-1})$
$\displaystyle=$ $\displaystyle
g(((1+\beta^{1})U_{k}-\beta^{1}U_{k-1}-\alpha^{1}\nabla
g(U_{k}W_{k}^{T})W_{k})((1+\beta^{2})W_{k}-\beta^{2}W_{k-1}-\alpha^{2}\nabla
g(U_{k}W_{k})^{T}U_{k})^{T})$ $\displaystyle=$ $\displaystyle
g((1+\beta^{1})(1+\beta^{2})U_{k}W_{k}^{T}-\alpha^{1}(1+\beta^{2})\nabla
g(U_{k}W_{k}^{T})W_{k}W_{k}^{T}-\alpha^{2}(1+\beta^{1})U_{k}U_{k}^{T}\nabla
g(U_{k}W_{k}^{T})$ $\displaystyle+\alpha^{1}\alpha^{2}\nabla
g(U_{k}W_{k}^{T})W_{k}U_{k}^{T}\nabla
g(U_{k}W_{k}^{T})+\beta^{1}\beta^{2}U_{k-1}W_{k-1}^{T}$
$\displaystyle{\color[rgb]{1,0,0}-\beta^{1}(1+\beta^{2})U_{k-1}W_{k}^{T}-(1+\beta^{1})\beta^{2}U_{k}W_{k-1}^{T}+\alpha^{2}\beta^{1}U_{k-1}U_{k}\nabla
g(U_{k}W_{k}^{T})+\alpha_{1}\beta^{2}\nabla
g(U_{k}W_{k}^{T})W_{k}W_{k-1}^{T}}$

This requires 9 matrix multiplications, which is 4 more than the 5 we require
if not include any momentum terms.

  5. 5.

Momentum on Both Matrices - Inexact SO: nine matrix multiplications to perform
SO compared to 3 matrix multiplications to use fixed learning and momentum
rates may be a significant increase in the iteration cost. Rather than
computing all 9 required matrix multiplications and precisely optimizing the
step sizes, we could instead compute 2 matrix multiplications plus one matrix
multiplication for each set of step sizes that we try. Specifically, as before
we track $M_{k}=U_{k}W_{k}^{T}$ and use 2 matrix multiplications to compute
$\nabla_{U}f(U_{k},W_{k})=\nabla g(M_{k})W_{k}$ and
$\nabla_{W}f(U_{k},W_{k})=\nabla g(M_{k})^{T}U_{k}$. We could then compute

$\displaystyle
f(U_{k}-\alpha^{1}\nabla_{U}f(U_{k},W_{k})+\beta^{1}(U_{k}-U_{k-1}),W_{k}-\alpha^{2}\nabla_{W}f(U_{k},W_{k})+\beta^{2}(W_{k}-W_{k-1})$
$\displaystyle=$ $\displaystyle
g((U_{k}-\alpha^{1}\nabla_{U}f(U_{k},W_{k})+\beta^{1}(U_{k}-U_{k-1}))(W_{k}-\alpha^{2}\nabla_{W}f(U_{k},W_{k})+\beta^{2}(W_{k}-W_{k-1}))^{T}),$

which requires one matrix multiplication for each prospective
$\\{\alpha^{1},\alpha^{2},\beta^{1},\beta^{2}\\}$ we test. Thus, this may be
advantageous over the precise approach if we can find suitable step sizes
using less than 7 guesses.

Modern variations on MF problems often include regularizers on the matrices
$U$ and $W$, or consider matrix completion settings where $g$ only depends on
a subset of the entries in $UW^{T}$. For example, both of these variations
appear in the probabilistic matrix factorization framework that is widely used
for recommender systems (Mnih and Salakhutdinov, 2007). Regularization does
not complicate the use of SO in the typical case where (given $UW^{T}$)
evaluating the regularizer is much less expensive than evaluating $UW^{T}$
such as penalizing $\|UW^{T}\|_{F}^{2}$. Matrix completion similarly supports
SO since the sparse dependency on $UW^{T}$ can be exploited in all the
required matrix multiplications.

##  Appendix C Log-Determinant Problems

We say that a problem is a log determinant (LD) problem if we have a symmetric
matrix variable $V\in\mathbb{R}^{d\times d}$ that is constrained to be
positive definite and the dominant cost in evaluating the objective is
computing the logarithm of the determinant of $V$, $\log|V|$. A classic
example is fitting the precision matrix of a multivariate Gaussian, where
$f(V)=\text{Tr}(SV)-\log|V|$ for the empirical covariance matrix $S$. In this
example the trace term in this objective is a linear composition but the log
determinant term is not. Note that we also obtain trace and log-determinant
terms in the more complex case of fitting a multivariate student $t$ (Cornish,
1954; Lindsey and Lindsey, 2006), and indeed log-determinants problems tend to
arise when fitting distributions under a change of variables. In this section
we discuss how SO is efficient for this type of problem under certain types of
updates.

A common approach to solving log-determinant problems is coordinate
optimization and block coordinate descent (Banerjee et al., 2006; Friedman et
al., 2008; Scheinberg and Rish, 2009; Hsieh et al., 2014). On each iteration,
these methods update either a single entry or a single column (and
corresponding row) of the matrix. In this setting of updating a single column
in each iteration, SO is efficient due to the multilinearity property of
determinants: if all columns are fixed except one, then the determinant is
linear with respect to that column. Thus, in block coordinate updates where
all columns are fixed except one, the update acts as if the problem was an
LCP.

However, rather than restricting to updating one element/column at a time, we
can alternately consider SO along general rank-1 search directions. For
example, if on iteration $k$ we consider the rank-1 direction
$u_{k}u_{k}^{T}$, by applying the matrix determinant lemma we have

$\displaystyle|V_{k+1}|$ $\displaystyle=|V_{k}+\alpha_{k}u_{k}u_{k}^{T}|$
$\displaystyle=(1+\alpha_{k}u_{k}^{T}V_{k}^{-1}u_{k})|V_{k}|$
$\displaystyle=(1+\alpha_{k}u_{k}^{T}\tilde{u}_{k})|V_{k}|,$

where we compute $\tilde{u}_{k}$ by solving $V_{k}\tilde{u}_{k}=u_{k}$. By
tracking the determinant $|V_{k}|$, we can thus evaluate $|V_{k+1}|$ for many
potential values of $\alpha_{k}$ at the cost of solving a single linear
system. We can additionally use the matrix inversion lemma to perform SO on
two rank-1 directions,

$\displaystyle|V_{k+1}|$
$\displaystyle=|V_{k}+\alpha_{k}^{1}u_{k}u_{k}^{T}+\alpha^{2}_{k}v_{k}v_{k}^{T}|$
$\displaystyle=(1+\alpha^{2}_{k}v_{k}^{T}(V_{k}+\alpha^{1}_{k}u_{k}u_{k}^{T})^{-1}v_{k})|V_{k}+\alpha^{1}_{k}u_{k}u_{k}^{T}|$
$\displaystyle=(1+\alpha^{2}_{k}v_{k}^{T}(V_{k}^{-1}-\frac{\alpha^{1}_{k}}{1+\alpha^{1}_{k}u_{k}^{T}V_{k}^{-1}u_{k}}V_{k}^{-1}u_{k}u_{k}^{T}V_{k}^{-1})v_{k})(1+\alpha^{1}_{k}u_{k}^{T}V_{k}^{-1}u_{k})|V_{k}|$
$\displaystyle=(1+\alpha^{2}_{k}v_{k}^{T}\tilde{v}_{k}-\frac{\alpha^{1}_{k}\alpha^{2}_{k}}{1+\alpha^{1}_{k}u_{k}^{T}\tilde{u}_{k}}v_{k}^{T}\tilde{u}_{k}u_{k}^{T}\tilde{v}_{k})(1+\alpha^{1}_{k}u_{k}^{T}\tilde{u}_{k})|V_{k}|$
$\displaystyle=((1+\alpha^{1}_{k}u_{k}^{T}\tilde{u}_{k})(1+\alpha_{k}^{2}v_{k}^{T}\tilde{v}_{k})-\alpha_{k}^{1}\alpha_{k}^{2}v_{k}^{T}\tilde{u}_{k}u_{k}^{T}\tilde{v}_{k})|V_{k}|,$

where $\tilde{v}_{k}$ solves $V_{k}\tilde{v}_{k}=v_{k}$. Thus, performing an
SO to set the step sizes $\alpha^{1}$ and $\alpha^{2}$ for two rank-1 search
directions can be done for the cost of solving two linear systems. Solving
these linear systems has a similar or smaller cost than computing/updating a
Cholesky factorization to compute the determinant and check positive-
definiteness.

## References

  * Andrei (2014) N. Andrei.  An accelerated subspace minimization three-term conjugate gradient algorithm for unconstrained optimization.  _Numerical Algorithms_ , 65:859–874, 2014.
  * Asi and Duchi (2019) H. Asi and J. C. Duchi.  The importance of better models in stochastic optimization.  _Proceedings of the National Academy of Sciences_ , 116(46):22924–22930, 2019.
  * Banerjee et al. (2006) O. Banerjee, L. E. Ghaoui, A. d’Aspremont, and G. Natsoulis.  Convex optimization techniques for fitting sparse gaussian graphical models.  In _ICML_ , pages 89–96, 2006.
  * Barzilai and Borwein (1988) J. Barzilai and J. M. Borwein.  Two-point step size gradient methods.  _IMA Journal of Numerical Analysis_ , 8:141–148, 1988.
  * Beck and Teboulle (2009) A. Beck and M. Teboulle.  A fast iterative shrinkage-thresholding algorithm for linear inverse problems.  _SIAM journal on imaging sciences_ , 2(1):183–202, 2009.
  * Belilovsky et al. (2019) E. Belilovsky, M. Eickenberg, and E. Oyallon.  Greedy layerwise learning can scale to imagenet.  In _International conference on machine learning_ , pages 583–593. PMLR, 2019.
  * Bell and Koren (2007) R. M. Bell and Y. Koren.  Scalable collaborative filtering with jointly derived neighborhood interpolation weights.  In _Seventh IEEE international conference on data mining (ICDM 2007)_ , pages 43–52. IEEE, 2007.
  * Bollapragada et al. (2018) R. Bollapragada, J. Nocedal, D. Mudigere, H. J. Shi, and P. Tang.  A progressive batching L-BFGS method for machine learning.  In _ICML_ , pages 620–629. PMLR, 2018.
  * Brown et al. (2020) T. Brown, B. Mann, N. Ryder, M. Subbiah, J. D. Kaplan, P. Dhariwal, A. Neelakantan, P. Shyam, G. Sastry, A. Askell, et al.  Language models are few-shot learners.  _NeurIPS_ , 33:1877–1901, 2020.
  * Cantrell (1969) J. W. Cantrell.  Relation between the memory gradient method and the Fletcher-Reeves method.  _Journal of Optimization Theory and Applications_ , 4(1):67–71, 1969.
  * Carmon et al. (2017) Y. Carmon, J. C. Duchi, O. Hinder, and A. Sidford.  “convex until proven guilty”: dimension-free acceleration of gradient descent on non-convex functions.  In _ICML_ , pages 654–663. PMLR, 2017.
  * Cartis and Otemissov (2022) C. Cartis and A. Otemissov.  A dimensionality reduction technique for unconstrained global optimization of functions with low effective dimensionality.  _Information and Inference: A Journal of the IMA_ , 11(1):167–201, 2022.
  * Cartis et al. (2023a) C. Cartis, E. Massart, and A. Otemissov.  Bound-constrained global optimization of functions with low effective dimensionality using multiple random embeddings.  _Mathematical Programming_ , 198(1):997–1058, 2023a.
  * Cartis et al. (2023b) C. Cartis, E. Massart, and A. Otemissov.  Global optimization using random embeddings.  _Mathematical Programming_ , 200(2):781–829, 2023b.
  * Chang and Lin (2011) C. C. Chang and C. J. Lin.  LIBSVM: A library for support vector machines.  _ACM Transactions on Intelligent Systems and Technology_ , 2:27:1–27:27, 2011.
  * Chouzenoux et al. (2010) E. Chouzenoux, J. Idier, and S. Moussaoui.  A majorize–minimize strategy for subspace optimization applied to image restoration.  _IEEE Transactions on Image Processing_ , 20(6):1517–1528, 2010.
  * Cohen (1972) A. I. Cohen.  Rate of convergence of several conjugate gradient algorithms.  _SIAM Journal on Numerical Analysis_ , 9(2):248–259, 1972.
  * Cohen et al. (2020) J. Cohen, S. Kaur, Y. Li, J. Z. Kolter, and A. Talwalkar.  Gradient descent on neural networks typically occurs at the edge of stability.  In _International Conference on Learning Representations_ , 2020.
  * Conn et al. (1994) A. R. Conn, P. L. Toint, A. Sartenaer, and N. Gould.  On iterated-subspace minimization methods for nonlinear optimization.  Technical report, P00024646, 1994.
  * Cornish (1954) E. A Cornish.  The multivariate t-distribution associated with a set of normal sample deviates.  _Australian Journal of Physics_ , 7:531, 1954.
  * Cragg and Levy (1969) E. E. Cragg and A. V. Levy.  Study on a supermemory gradient method for the minimization of functions.  _Journal of Optimization Theory and Applications_ , 4(3), 1969.
  * De et al. (2016) S. De, A. Yadav, D. Jacobs, and T. Goldstein.  Big batch SGD: automated inference using adaptive batch sizes.  _arXiv:1610.05792_ , 2016.
  * Deng et al. (2009) J. Deng, W. Dong, R. Socher, L.-J. Li, K. Li, and F. F. Li.  Imagenet: A large-scale hierarchical image database.  In _2009 IEEE conference on computer vision and pattern recognition_ , pages 248–255. Ieee, 2009.
  * Donahue et al. (2014) J. Donahue, Y. Jia, O. Vinyals, J. Hoffman, N. Zhang, E. Tzeng, and T. Darrell.  Decaf: a deep convolutional activation feature for generic visual recognition.  In _ICML_ , pages 647–655. PMLR, 2014.
  * Drori and Taylor (2020) Y. Drori and A. B. Taylor.  Efficient first-order methods for convex minimization: a constructive approach.  _Mathematical Programming_ , 184(1):183–220, 2020\.
  * Drusvyatskiy and Paquette (2019) D. Drusvyatskiy and C. Paquette.  Efficiency of minimizing compositions of convex functions and smooth maps.  _Mathematical Programming_ , 178:503–558, 2019.
  * Dua and Graff (2017) D. Dua and C. Graff.  UCI machine learning repository, 2017.
  * Duchi and Ruan (2018) J. C. Duchi and F. Ruan.  Stochastic methods for composite and weakly convex optimization problems.  _SIAM Journal on Optimization_ , 28(4):3229–3259, 2018.
  * Elad et al. (2007) M. Elad, B. Matalon, and M. Zibulevsky.  Coordinate and subspace optimization methods for linear least squares with non-quadratic regularization.  _Applied and Computational Harmonic Analysis_ , 23(3):346–367, 2007.
  * Fernández-Delgado et al. (2014) M. Fernández-Delgado, E. Cernadas, S. Barro, and D. Amorim.  Do we need hundreds of classifiers to solve real world classification problems?  _JMLR_ , 15(1):3133–3181, 2014.
  * Fletcher and Reeves (1964) R. Fletcher and C. M. Reeves.  Function minimization by conjugate gradients.  _The computer journal_ , 7(2):149–154, 1964.
  * Friedlander and Schmidt (2012) M. P. Friedlander and M. Schmidt.  Hybrid deterministic-stochastic methods for data fitting.  _SIAM Journal on Scientific Computing_ , 34(3):A1380–A1405, 2012.
  * Friedman et al. (2008) J. Friedman, T. Hastie, and R. Tibshirani.  Sparse inverse covariance estimation with the graphical lasso.  _Biostatistics_ , 9(3):432–441, 2008.
  * Galli et al. (2023) L. Galli, H. Rauhut, and M. Schmidt.  Don’t be so monotone: Relaxing stochastic line search in over-parameterized models.  _arXiv preprint arXiv:2306.12747_ , 2023.
  * Grippo et al. (1986) L. Grippo, F. Lampariello, and S. Lucidi.  A nonmonotone line search technique for Newton’s method.  _SIAM journal on Numerical Analysis_ , 23(4):707–716, 1986.
  * Guminov et al. (2023) S. Guminov, A. Gasnikov, and I. Kuruzov.  Accelerated methods for weakly-quasi-convex optimization problems.  _Computational Management Science_ , 20(1):1–19, 2023.
  * Gupta et al. (2023) S. D. Gupta, R. M. Freund, A. Sun, and A. Taylor.  Nonlinear conjugate gradient methods: worst-case convergence rates via computer-assisted analyses.  _arXiv preprint arXiv:2301.01530_ , 2023.
  * Hager and Zhang (2006) W. W. Hager and H. Zhang.  A survey of nonlinear conjugate gradient methods.  _Pacific journal of Optimization_ , 2(1):35–58, 2006.
  * Hestenes and Stiefel (1952) M. R. Hestenes and E. Stiefel.  Methods of conjugate gradients for solving linear systems.  _Journal of Research of the National Bureau of Standards_ , 49(6):409–436, 1952.
  * Hinder (2018) O. Hinder.  Cutting plane methods can be extended into nonconvex optimization.  In _COLT_ , pages 1451–1454. PMLR, 2018.
  * Hsieh et al. (2014) C. J. Hsieh, M. A. Sustik, I. S. Dhillon, P. Ravikumar, et al.  Quic: quadratic approximation for sparse inverse covariance estimation.  _JMLR_ , 15(1):2911–2947, 2014.
  * Huang et al. (2006) G. B. Huang, Q.Y. Zhu, and C. K. Siew.  Extreme learning machine: theory and applications.  _Neurocomputing_ , 70(1-3):489–501, 2006.
  * Jacot et al. (2018) A. Jacot, F. Gabriel, and C. Hongler.  Neural tangent kernel: Convergence and generalization in neural networks.  _NeurIPS_ , 31, 2018.
  * Jin et al. (2024) Q. Jin, R. Jiang, and A. Mokhtari.  Non-asymptotic global convergence rates of bfgs with exact line search.  _arXiv preprint arXiv:2404.01267_ , 2024.
  * Karimi and Vavasis (2017) S. Karimi and S. Vavasis.  A single potential governing convergence of conjugate gradient, accelerated gradient and geometric descent.  _arXiv:1712.09498_ , 2017.
  * Kingma and Ba (2014) D. P. Kingma and J. Ba.  Adam: A method for stochastic optimization.  _arXiv:1412.6980_ , 2014.
  * Krizhevsky et al. (2009) A. Krizhevsky, G. Hinton, et al.  Learning multiple layers of features from tiny images.  2009\.
  * Krizhevsky et al. (2012) A. Krizhevsky, I. Sutskever, and G. E. Hinton.  Imagenet classification with deep convolutional neural networks.  _Advances in neural information processing systems_ , 25, 2012.
  * Kunstner et al. (2022) F. Kunstner, J. Chen, J. W. Lavington, and M. Schmidt.  Noise is not the main factor behind the gap between sgd and adam on transformers, but sign descent might be.  In _The Eleventh International Conference on Learning Representations_ , 2022.
  * LeCun et al. (1998) Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner.  Gradient-based learning applied to document recognition.  _Proceedings of the IEEE_ , 86(11):2278–2324, 1998.
  * LeCun et al. (2002) Y. LeCun, L. Bottou, G. B. Orr, and K.-R. Müller.  Efficient backprop.  In _Neural networks: Tricks of the trade_ , pages 9–50. Springer, 2002.
  * Lewis and Wright (2016) A. S. Lewis and S. J. Wright.  A proximal method for composite minimization.  _Mathematical Programming_ , 158:501–546, 2016.
  * Li and Lin (2022) H. Li and Z. Lin.  Restarted nonconvex accelerated gradient descent: no more polylogarithmic factor in the $O(\varepsilon^{-7/4})$ complexity.  In _ICML_ , pages 12901–12916. PMLR, 2022.
  * Lindsey and Lindsey (2006) J. K. Lindsey and P. J. Lindsey.  Multivariate distributions with correlation matrices for nonlinear repeated measurements.  _Computational statistics & data analysis_, 50(3):720–732, 2006.
  * Miele and Cantrell (1969) A. Miele and J. W. Cantrell.  Study on a memory gradient method for the minimization of functions.  _Journal of Optimization Theory and Applications_ , 3(6), 1969.
  * Mishkin et al. (2022) A. Mishkin, A. Sahiner, and M. Pilanci.  Fast convex optimization for two-layer relu networks: Equivalent model classes and cone decompositions.  In _International Conference on Machine Learning_ , pages 15770–15816. PMLR, 2022.
  * Mnih and Salakhutdinov (2007) A. Mnih and R. R. Salakhutdinov.  Probabilistic matrix factorization.  _NeurIPS_ , 20, 2007.
  * Moré and Thuente (1994) J. J. Moré and D. J. Thuente.  Line search algorithms with guaranteed sufficient decrease.  _ACM Transactions on Mathematical Software (TOMS)_ , 20(3):286–307, 1994.
  * Narkiss and Zibulevsky (2005a) G. Narkiss and M. Zibulevsky.  Sequential subspace optimization method for large-scale unconstrained problems.  Technical report, Technion - Israel Institute of Technology, 2005a.
  * Narkiss and Zibulevsky (2005b) G. Narkiss and M. Zibulevsky.  Support vector machine via sequential subspace optimization.  Technical report, Technion - Israel Institute of Technology, 2005b.
  * Nemirovski (1982) A. Nemirovski.  Orth-method for smooth convex optimization.  _Engineering Cybernetics_ , 20(2):937–947, 1982\.
  * Nemirovski and Yudin (1983) A Nemirovski and D Yudin.  Information-based complexity of mathematical programming.  _Izvestia AN SSSR, Ser. Tekhnicheskaya Kibernetika (the journal is translated to English as Engineering Cybernetics. Soviet J. Computer & Systems Sci.)_, 1, 1983.
  * Nesterov (1983) Y. Nesterov.  A method for unconstrained convex minimization problem with the rate of convergence $o(1/k^{2})$.  In _Doklady AN USSR_ , pages 543–547, 1983.
  * Nesterov (2007) Y. Nesterov.  Modified Gauss–Newton scheme with worst case guarantees for global performance.  _Optimisation Methods and Software_ , 22(3):469–483, 2007.
  * Nesterov (2013) Y. Nesterov.  Gradient methods for minimizing composite functions.  _Mathematical programming_ , 140(1):125–161, 2013\.
  * Nesterov et al. (2021) Y. Nesterov, A. Gasnikov, S. Guminov, and P. Dvurechensky.  Primal–dual accelerated gradient methods with small-dimensional relaxation oracle.  _Optimization Methods and Software_ , 36(4):773–810, 2021.
  * Nocedal (1980) J. Nocedal.  Updating quasi-Newton matrices with limited storage.  _Mathematics of computation_ , 35(151):773–782, 1980.
  * Nocedal and Wright (2006) J. Nocedal and S. J. Wright.  _Numerical Optimization, 2nd Ed._ Springer, 2006.
  * Polak and Ribiere (1969) E. Polak and G. Ribiere.  Note sur la convergence de méthodes de directions conjuguées.  _Revue française d’informatique et de recherche opérationnelle. Série rouge_ , 3(16):35–43, 1969.
  * Polyak (1964) B. Polyak.  Some methods of speeding up the convergence of iteration methods.  _USSR Computational Mathematics and Mathematical Physics_ , 4:1–17, 1964.
  * Rajeswaran et al. (2017) A. Rajeswaran, K. Lowrey, E. V. Todorov, and S. M. Kakade.  Towards generalization and simplicity in continuous control.  _Advances in Neural Information Processing Systems_ , 30, 2017.
  * Raydan (1997) M. Raydan.  The Barzilai and Borwein gradient method for the large scale unconstrained minimization problem.  _SIAMOpt_ , 7(1):26–33, Feb 1997.
  * Richardson et al. (2016) E. Richardson, R. Herskovitz, B. Ginsburg, and M. Zibulevsky.  Seboost-boosting stochastic learning using subspace optimization techniques.  _Advances in Neural Information Processing Systems_ , 29, 2016.
  * Robbins and Monro (1951) H. Robbins and S. Monro.  A stochastic approximation method.  _The annals of mathematical statistics_ , pages 400–407, 1951.
  * Scheinberg and Rish (2009) K. Scheinberg and I. Rish.  Sinco-a greedy coordinate ascent method for sparse inverse covariance selection problem.  _preprint_ , 2009.
  * Schulman et al. (2015) J. Schulman, S. Levine, P. Abbeel, M. Jordan, and P. Moritz.  Trust region policy optimization.  In _International conference on machine learning_ , pages 1889–1897. PMLR, 2015.
  * Schulman et al. (2017) J. Schulman, F. Wolski, P. Dhariwal, A. Radford, and O. Klimov.  Proximal policy optimization algorithms.  _arXiv preprint arXiv:1707.06347_ , 2017.
  * Shanno and Phua (1978) D. F. Shanno and K. H. Phua.  Matrix conditioning and nonlinear optimization.  _Mathematical Programming_ , 14:149–160, 1978.
  * Shea and Schmidt (2023) B. Shea and M. Schmidt.  Greedy Newton: Newton’s method with exact line search.  _NeurIPS Workshop on Optimization for Machine Learning_ , 2023.
  * Shi et al. (2021) B. Shi, S. S. Du, M. I. Jordan, and W. J. Su.  Understanding the acceleration phenomenon via high-resolution differential equations.  _Mathematical Programming_ , pages 1–70, 2021.
  * Sorber et al. (2016) L. Sorber, I. Domanov, M. Van Barel, and L. De Lathauwer.  Exact line and plane search for tensor optimization.  _Computational Optimization and Applications_ , 63:121–142, 2016.
  * Truong and Nguyen (2018) T. T. Truong and T. H. Nguyen.  Backtracking gradient descent method for general $C^{1}$ functions, with applications to deep learning.  _arXiv:1808.05160_ , 2018.
  * Vaswani et al. (2019) S. Vaswani, A. Mishkin, I. Laradji, M. Schmidt, G. Gidel, and S. Lacoste-Julien.  Painless stochastic gradient: interpolation, line-search, and convergence rates.  In _NeurIPS_ , volume 32, 2019.
  * Vaswani et al. (2020) S. Vaswani, I. Laradji, F. Kunstner, S. Y. Meng, M. Schmidt, and S. Lacoste-Julien.  Adaptive gradient methods converge faster with over-parameterization (but you should do a line-search).  _arXiv:2006.06835_ , 2020.
  * Wang et al. (2016) Z. Wang, F. Hutter, M. Zoghi, D. Matheson, and N. De Feitas.  Bayesian optimization in a billion dimensions via random embeddings.  _Journal of Artificial Intelligence Research_ , 55:361–387, 2016.
  * Wang and Yuan (2006) Z. H. Wang and Y. X. Yuan.  A subspace implementation of quasi-Newton trust region methods for unconstrained optimization.  _Numerische Mathematik_ , 104:241–269, 2006.
  * Wolfe (1969) P. Wolfe.  Convergence conditions for ascent methods.  _SIAM review_ , 11(2):226–235, 1969.
  * Zhang et al. (2022) Y. Zhang, C. Chen, N. Shi, R. Sun, and Z. Q. Luo.  Adam can converge without any modification on update rules.  _arXiv:2208.09632_ , 2022.
  * Zibulevsky (2008) M. Zibulevsky.  Sesop-tn: combining sequential subspace optimization with truncated Newton method, 2008.
  * Zibulevsky (2013) M. Zibulevsky.  Speeding-up convergence via sequential subspace optimization: Current state and future directions.  _arXiv preprint arXiv:1401.0159_ , 2013.