SANIA: Polyak-type Optimization Framework Leads to Scale Invariant Stochastic Algorithms

## Abstract

Abstract Adaptive optimization methods are widely recognized as among the most popular approaches for training Deep Neural Networks (DNNs). Techniques such as Adam, AdaGrad, and AdaHessian utilize a preconditioner that modifies the search direction by incorporating information about the curvature of the objective function. However, despite their adaptive characteristics, these methods still require manual fine-tuning of the step-size. This, in turn, impacts the time required to solve a particular problem. This paper presents an optimization framework named SANIA to tackle these challenges. Beyond eliminating the need for manual step-size hyperparameter settings, SANIA incorporates techniques to address poorly scaled or ill-conditioned problems. We also explore several preconditioning methods, including Hutchinson’s method , which approximates the Hessian diagonal of the loss function. We conclude with an extensive empirical examination of the proposed techniques across classification tasks, covering both convex and non-convex contexts.

## 1 INTRODUCTION

Machine Learning (ML), especially Deep Neural Networks (DNNs), has emerged as a transformative tool, setting the stage for unprecedented advances across many disciplines, including computer vision
(Krizhevsky et al.,, ; Simonyan and Zisserman,, ; He et al.,, ) and natural language processing (Wolf et al.,, ; Mikolov et al.,, ; Devlin et al.,, ; Radford et al.,, ), as well as science
(Xie and Grossman,, ; Gómez-Bombarelli et al.,, )
and engineering (Bello et al.,, ; LeCun et al.,, )
to name a few.
The enormous potential of these models is enabled through the efficacy of the optimization methods that train them.
In the domain of ML the training task can be expressed as solving the following problem

$$ $\textstyle{\min}_{w\in\mathbb{R}^{d}}f(w):=\tfrac{1}{n}\sum_{i=1}^{n}f_{i}(w),$ (1) $$

where $w\in\mathbb{R}^{d}$ represents the weight parameter, and each $f_{i}:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is a sufficiently smooth function.
To provide a practical context, consider a dataset denoted as $\{(x_{i},y_{i})\}_{i=1}^{n}$, where $x_{i}\in\mathbb{R}^{d}$ is the data sample and $y_{i}\in\mathbb{R}$ represents the label corresponding to that sample. If $f_{i}(w)=\frac{1}{2}({x_{i}}^{T}w-y_{i})^{2}$, this optimization problem gives rise to the well-known least squares problem. Similarly, if $f_{i}(w)=\log(1+e^{-y_{i}{x_{i}}^{T}w})$,
we get logistic regression problem.
Stochastic Gradient Descent.
To address problem ([1](#S1.E1)), one of the fundamental techniques employed is Stochastic Gradient Descent (SGD)
(Robbins and Monro,, ; Polyak,, ; Polyak and Juditsky,, ; Nemirovski et al.,, ; Bottou et al.,, ).
This method iteratively updates the weight parameter $w$ according to the following scheme:

$$ $w_{t+1}=w_{t}-\gamma_{t}\nabla f_{i}(w_{t}),$ (2) $$

where $\gamma_{t}$ is the step-size schedule and $i\subset[n]:=\{1,2,\dots,n\}$ is chosen uniformly as random.
Unfortunately, the optimal step-size schedule often relies on problem-specific parameters, such as the Lipschitz-smoothness constant and the level of stochastic gradient noise, which are frequently not accessible. Consequently, achieving an optimal step-size typically demands a substantial amount of tuning, which can be quite costly in practical applications. The influence of the step-size can be observed in Figure  [1](#S1.F1) where our proposed method achieves 100% accuracy rate after 10 epochs of training while other algorithms require fine-tuning their step-size in order to achieve the same result.

Figure: Figure 1: Accuracy of the models optimized by SANIA and other popular algorithms with different learning rates $lr\in[2^{-20},2^{5}]$ after $10$ epochs of training on colon-cancer with batch-size$=16$ and logistic regression criterion.
Refer to caption: x1.png

Numerous methodologies have been developed to tackle this issue. The first approach that reduces the number of parameters to tune is the AdaGrad method by Duchi et al., (); Li and Orabona, (); Ward et al., ().
An additional challenge arises from the fact that using the same learning rate for each feature
$j\in[d]$ might not yield the best performance. To address this, diagonal preconditioning techniques have been employed in the SGD setting by methods such as AdaGrad by Duchi et al., (), RMSProp by Tieleman et al., (), Adam by Kingma and Ba, (), AMSGrad by Reddi et al., (), AdamW by Loshchilov and Hutter, (),
AdaHessian by Yao et al., (), AdaDelta by Zeiler, (), and OASIS by Jahani et al., (). However, all of these methods still require a considerable degree of parameter tuning to achieve optimal performance. Another approach is associated with parameter-free regret minimization for online learning problems, as discussed in various papers (Mcmahan and Streeter,, ; McMahan and Orabona,, ; Orabona and Pál,, ; Orabona and Tommasi,, ; Orabona,, ; Carmon and Hinder,, ; Ivgi et al.,, ; Defazio and Mishchenko,, ; Cutkosky et al.,, ; Mishchenko and Defazio,, ). Finally, in our paper, we explore the Stochastic Polyak step-size approach as an adaptive parameter-free method.
Stochastic Polyak step-size (SPS) Methods. Polyak step-size method was first proposed for the subgradient by Boris T. Polyak, (, ). Recently, stochastic Polyak step-size was proposed by Oberman and Prazeres, (); Berrada et al., (); Loizou et al., (); Gower et al., (); Orvieto et al., (). Subsequently, lots of variants of SPS have emerged, such as mSPS by D’Orazio et al., () and AdaSLS by Jiang and Stich, (). To further relax the requirements for interpolation condition in SPS, many attempts have been made by Gower et al., (); Orvieto et al., (); Garrigos et al., (); Schaipp et al., (). A variant of second-order expansion for SPS was presented by Li et al., (). Next we describe the main idea of Polyak step-size in more detail.
To derive the deterministic Polyak step-size, let us consider a convex function $f(w)$ and the step ([2](#S1.E2)). We obtain the step-size from the following upper-bound on the distance from the current point $w_{t+1}$ to the minimum $w^{\ast}$:

$$ $\displaystyle\|w_{t+1}-w^{\ast}\|^{2}=\|w_{t}-w^{\ast}\|^{2}+\|\gamma_{t} \nabla f(w_{t})\|^{2}$ $\displaystyle\hskip 68.28644pt-2\gamma_{t}\left\langle\nabla f(w_{t}),w_{t}-w^ {\ast}\right\rangle$ $\displaystyle\leq\|w_{t}-w^{\ast}\|^{2}+\gamma_{t}^{2}\|\nabla f(w_{t})\|^{2}- 2\gamma_{t}(f(w_{t})-f(w^{\ast})).$ $$

Minimizing the right hand side by $\gamma_{t}$, we get

$$ $\gamma_{t}=\tfrac{f(w_{t})-f(w^{\ast})}{\|\nabla f(w_{t})\|^{2}}$ (3) $$

Similarly, in the stochastic case, the Stochastic Polyak step-size (SPS) is defined as

$$ $w_{t+1}=w_{t}-\frac{f_{i}(w_{t})-f_{i}^{\ast}}{\|\nabla f_{i}(w_{t})\|^{2}} \nabla f_{i}(w_{t}),$ (4) $$

where $f_{i}^{\ast}$ is a minimal value of function $f_{i}(w)$.
Another way to derive this formulation is by solving the following optimization problem:

$$ $\displaystyle w_{t+1}=\textstyle{\arg\min}_{w\in\mathbb{R}^{d}}\|w-w_{t}\|_{2} ^{2}$ s.t. $\displaystyle f_{i}(w_{t})+\langle\nabla f_{i}(w_{t}),w-w_{t}\rangle=f_{i}^{ \ast},$ (5) $$

where ([4](#S1.E4)) is an explicit formulation of ([5](#S1.E5)).
In case where the value of $f_{i}^{\ast}$ is known and set to 0 0 for all $i$ (a common scenario in over-parameterized deep neural networks) we obtain this simplified expression for ([4](#S1.E4)):

$$ $w_{t+1}=w_{t}-\tfrac{f_{i}(w_{t})}{\|\nabla f_{i}(w_{t})\|^{2}}\nabla f_{i}(w_ {t}).$ (6) $$

This condition, referred to as the ”interpolation condition”, is expressed as $f_{i}^{\ast}=0$.
Preconditioning / Feature scaling.
Preconditioning is a technique used to improve the convergence rate of algorithms applied to data that may exhibit poor scaling or ill-conditioning. Algorithms leveraging preconditioning typically follow a generic update rule, which can be expressed as

$$ $w_{t+1}=w_{t}-\gamma_{t}B_{t}^{-1}m_{t},$ (7) $$

where $B_{t}\in\mathbb{R}^{d\times d}$ is an invertible positive definite matrix, and $m_{t}$ is a gradient or its approximation. The origin of such a step is Newton method by Newton, (); Raphson, (); ; ; Kantorovich, () which uses the exact Hessian to precondition the gradient of the objective function, i.e. $B_{t}=\nabla^{2}f(w_{t})$ and $m_{t}=\nabla f_{i}(w_{t})$. Newton method can be very effective for minimizing convex objectives. However, the prohibitive cost of computing and inverting the Hessian matrix, together with issues around negative eigenvalues, makes this approach impractical for machine learning tasks. To address this issue, one can use methods that never define the Hessian of the objective function explicitly but rather use its approximation or solve the Newton system using iterative algorithms (Martens et al.,, ).
Quasi-Newton methods (QN).
Methods that construct an approximation of the (inverse) Hessian date back to the 70s such as BFGS by (Broyden,, ; Fletcher,, ; Goldfarb,, ; Shanno,, ), L-BFGS by (Nocedal,, ; Liu and Nocedal,, ), and SR-1 by (Conn et al.,, ; Khalfan et al.,, ).
These optimization methods take advantage of a cheap way to build (inverse) Hessian matrix estimation algorithms based on past gradient information.
One of the most prominent QN method is Symmetric Rank 1(SR-1) which recursively approximates the Hessian as follows

$$ $B_{t+1}=B_{t}+\tfrac{(y_{t}-B_{t}s_{t})(y_{t}-B_{t}s_{t})^{\top}}{(y_{t}-B_{t} s_{t})^{\top}s_{t}},$ (8) $$

where $s_{t}=w_{t+1}-w_{t}$ and $y_{t}=\nabla f_{i}(w_{t+1})-\nabla f_{i}(w_{t})$. Although, SR-1 update only makes a rank-1 change to the previous Hessian approximation and evidently has a simple form, in practice it displays better convergence to the true Hessian than other similar methods like BFGS. Another useful property of this approximation is self-complementarity, which means that we can find the inverse Hessian approximation $B_{t}^{-1}$ using the same vector pair $s_{t}$ and $y_{t}$: $B_{t+1}^{-1}=B_{t}^{-1}+\frac{(s_{t}-B_{t}^{-1}y_{t})(s_{t}-B_{t}^{-1}y_{t})^{
\top}}{(s_{t}-B_{t}^{-1}y_{t})^{\top}y_{t}}$. Note, that this approximation method does not necessarily generate a positive definite matrix.
Contributions.
Before delving into the details, we outline the primary contributions of this work:

- •
We present the General Framework for Preconditioned and Second-order Polyak methods. This framework covers classical optimization methods, provides valuable insights into Polyak step-size methods, and enables the development of novel Polyak step-size methods.
- •
We propose the first Stochastic Cubic Newton method with Polyak step-size.
- •
We introduce the new scale invariant versions of AdaGrad and Adam, which make them invariant to some basis transformations.
- •
We conduct comprehensive experiments encompassing a diverse range of scenarios, including both convex and non-convex settings.

Organisation.
In this paper, we have consolidated our findings and integrated them into a comprehensive framework presented in Section [2](#S2). Additionally, Section [3](#S3) offers a detailed presentation of the results from our experiments.
Notation and Assumptions.
We introduce the notation used throughout the paper and state the underlying assumptions that guide our analysis.

We equip the primal space $w\in\mathbf{E}$ and the dual space $g\in\mathbf{E}^{\ast}$ with the conjugate norms $\|w\|$ and $\|g\|_{\ast}$, respectively. As a special case, for a positive definite matrix $B\in\mathbb{R}^{d\times d}$, we introduce the conjugate Euclidean norms as follows: $\|w\|_{B}=\left\langle Bw,w\right\rangle^{1/2}$ and $\|g\|_{B^{-1}}=\left\langle g,B^{-1}g\right\rangle^{1/2}$. As an example, $\nabla f(w)\in\mathbf{E}^{\ast}$ and $\nabla^{2}f(w)h\in\mathbf{E}^{\ast}$ for $h\in\mathbf{E}$. We define the operator $\odot$ as a component-wise product between two vectors, also known as the Hadamard product. For the vector $w$, $w^{2}$ and $\sqrt{w}$ means component-wise square and square root, respectively. We represent $\operatorname{diag}(w)$ as a diagonal matrix of a given vector $v$ and a vector $\operatorname{diagonal}(H)\in\mathbb{R}^{d}$ as the diagonal of a matrix $H\in\mathbb{R}^{d\times d}$. For simplicity, we denote $g_{t}=\nabla f_{i}(w_{t})$ and $H_{t}=\nabla^{2}f_{i}(w_{t})$ if it is not defined differently. Also, we denote an action of the linear operator as $B[h]^{2}=\left\langle Bh,h\right\rangle$. Unless otherwise stated, our default assumption is that Assumption [1](#S1) holds true.

## 2 SANIA – GENERAL FRAMEWORK

### 2.1 General Framework

In this section, we propose a general framework ([10](#S1.E10)) for preconditioning stochastic Polyak step-size methods. This framework generalizes some well-known first-order, second-order, and Quasi-Newton methods from Polyak step-size perspective. The main feature of the framework is that it highlights some insights about SPS and provides an instrument to generalize existing methods as Polyak step-size methods. It makes them adaptive and parameter-free in the SPS setting. The generality of this framework makes it difficult to propose an explicit step. Therefore, we will focus on the most promising cases and provide their explicit formulations to introduce new methods.
In the following section, we will demonstrate the problem settings required to derive existing and proposed methods using SANIA ([10](#S1.E10)). We note that if any particular variable from the General Framework is not mentioned explicitly it is assumed to be fixed at zero.

### 2.2 Existing methods

SGD. Let us first derive an update rule for the most frequently used variant of Stochastic Gradient Descent (SGD) method using SANIA ([10](#S1.E10)).
We set the parameters as follows:
$\displaystyle\tau_{t}=\gamma_{t},\,m_{t}=\nabla f_{i}(w_{t}),\,D_{t}=0,\,B_{t}
=\mathit{I}.$
The explicit method ([2](#S1.E2)) is the solution of the following implicit problem:
$\begin{gathered}w_{t+1},\alpha_{t+1}=\mathop{\arg\min}_{w\in\mathbb{R}^{d},
\alpha\in\mathbb{R}}\tfrac{1}{2}\|w-w_{t}\|^{2}_{2}+\gamma_{t}\alpha\\
\quad\text{s.t.}\;f_{i}(w_{t})+\langle\nabla f_{i}(w_{t}),w-w_{t}\rangle\leq
\alpha.\end{gathered}$
(11)
The proof is presented in Appendix [B.1](#A2.SS1). Note, that normally $\alpha$ is an upper bound for $f_{i}^{\ast}$. Hence, if $f_{i}^{\ast}$ is known, we can fix $\alpha=f_{i}^{\ast}$. This leads us to Stochastic Polyak step-size method.
Stochastic Polyak step-size (SPS). The update rule for Stochastic Gradient Descent with Polyak step-size can be derived as follows:
We set the parameters as follows:
$\displaystyle\alpha=f_{i}^{\ast},\,\tau_{t}=\gamma_{t},\,m_{t}=\nabla f_{i}(w_
{t}),\,D_{t}=0,\,B_{t}=\mathit{I},$
and solve the following problem:
$\begin{gathered}w_{t+1}=\textstyle{\arg\min}_{w\in\mathbb{R}^{d}}\tfrac{1}{2}
\|w-w_{t}\|^{2}_{2}\\
\hskip 42.67912pt\text{s.t.}\;f_{i}(w_{t})+\langle\nabla f_{i}(w_{t}),w-w_{t}
\rangle\leq f_{i}^{\ast}.\end{gathered}$
(12)
We show that ([4](#S1.E4)) is an explicit version of ([12](#S2.E12)) in Appendix [B.2](#A2.SS2). If $f_{i}^{\ast}$ is known for example in case of interpolation with Assumption [1](#S1), then the method is adaptive and parameter-free. Otherwise, we replaced step-size parameter $\gamma_{t}$ with parameter $f_{i}^{\ast}$. Next, we show that a similar transition could be done for different methods.
Preconditioned SGD.
Preconditioning is used to introduce curvature information into SGD ([7](#S1.E7)). We precondition the stochastic gradient approximation, denoted as $m_{t}$, with a positive definite matrix $B_{t}\succ 0$. There are many methods that fit this description, ranging from the classical Damped Newton method and Quasi-Newton methods (like BFGS) to modern diagonal preconditioning techniques such as Adam, AdaGrad, and Hutchinson method. We can derive Preconditioned SGD from ([10](#S1.E10)).
With $0<\gamma_{t}\leq 1$ as a step-size, we choose the parameters as follows:
$\displaystyle\tau_{t}=\gamma_{t},\,m_{t}=m_{t},\,D_{t}=0,\,B_{t}=B_{t},$
and solve the following problem:
$\begin{gathered}w_{t+1},\alpha_{t+1}=\textstyle{\arg\min}_{w\in\mathbb{R}^{d},
\alpha\in\mathbb{R}}\tfrac{1}{2}\|w-w_{t}\|^{2}_{B_{t}}+\gamma_{t}\alpha\\
\quad\text{s.t.}\;f_{i}(w_{t})+\langle m_{t},w-w_{t}\rangle\leq\alpha.\end{gathered}$
(13)
We get the next explicit step:
$w_{t+1}=w_{t}-\gamma_{t}B_{t}^{-1}m_{t}.$
(14)
Next, we describe some preconditioning methods.
AdaGrad is an adaptive optimization method that approximates the Hessian of the objective function using the cumulative squared gradient information to scale the learning rates. Accumulation of all previous gradients in the preconditioner $B_{t}$ leads to decay in the learning rate $\gamma_{t}$ which increases performance for sparse settings (non-frequent features) at the cost of degrading in case of dense settings.

Adam is incorporating both adaptive learning rates and momentum. The update rule involves the computation of the moving average of both the first and second moments of the gradients. The first moment ($\beta_{1}$) is the mean of the gradients, and the second moment ($\beta_{2}$) is the uncentered variance of the gradients.
The Adam preconditioning is derived by:
$\displaystyle m_{t}=\tfrac{(1-\beta_{1})\sum_{j=1}^{t}\beta_{1}^{t-i}g_{j}}{1-
\beta_{1}^{t}},$
$\displaystyle B_{t}=\operatorname{diag}\left(\textstyle{\sqrt{\tfrac{(1-\beta_
{2})\sum_{j=1}^{t}\beta_{2}^{t-j}g_{j}^{2}}{1-\beta_{2}^{t}}}}\right),$
where $0<\beta_{1},\beta_{2}<1$ are two hyperparameters referred to as first and second moment coefficients. The biased estimates are corrected by dividing them by the bias correction terms, which are powers of the decay rates $\beta_{1}$ and $\beta_{2}$, respectively.
Hutchinson’s method is employed to estimate the diagonal of the Hessian matrix (Hutchinson,, ). To achieve this, the method utilizes only a handful of Hessian-vector products, which can be efficiently computed using backpropagation (Christianson,, ).
Specifically, the product of a Hessian matrix
$\nabla^{2}f(w)$ and a vector $h$ can be computed
through a directional derivative of the gradient, given by
$\tfrac{d}{dt}\left.\nabla f(w+th)\right|_{t=0}=\nabla^{2}f(w)h.$
Hutchinson’s method leverages Hessian-vector products to estimate the diagonal through
$\operatorname{diag}(\nabla^{2}f(w))=\mathbb{E}[h\odot(\nabla^{2}f(w)h)],$
where $h$ is a random vector with Rademacher distribution(^1^11$h_{j}\in\{-1,+1\}$ with equal probability.) or a normal distribution as discussed in (Bekas et al.,, ) and Lemma [B.4](#A2.SS4) in Appendix.
Utilizing this identity, we can estimate the Hessian diagonal by a weighted average of each iteration’s result:

$$ $B_{t}=\beta B_{t-1}+(1-\beta)\operatorname{diag}(h\odot\nabla^{2}f_{i_{t}}(w_{ t})h),$ $$

where $\beta\in(0,1)$ is a momentum parameter, $i_{t}$ is a number of a random function on the step $t$, and
$B_{0}=\tfrac{1}{k}\sum_{j=1}^{k}\operatorname{diag}(h_{j}\odot\nabla^{2}f_{j}(
w_{0})h_{j}),$
where $k$ is a number of functions for initialization of the approximation. To ensure
$B_{t}$ remains positive definite, especially in the face of potential non-convexities in the loss functions, we apply truncation by positive number $\mu$ and retain only the absolute values of elements given by $(B_{t})_{j,j}=\max\{\mu,|B_{t}|_{j,j}\}$.
Preconditioned SPS.
Similarly to SGD and SPS, Polyak step-size could be introduced for Preconditioned SGD methods. Preconditioned SPS (PSPS) was presented by Abdukhakimov et al., (). It can be also derived from SANIA for $B_{t}\succ 0$.
We choose the parameters as follows:
$\displaystyle\alpha=f_{i}^{\ast},\,\tau_{t}=\gamma_{t},\,m_{t}=m_{t},\,D_{t}=0
,\,B_{t}=B_{t},$
and solve the following problem:
$\begin{gathered}w_{t+1}=\textstyle{\arg\min}_{w\in\mathbb{R}^{d}}\tfrac{1}{2}
\|w-w_{t}\|^{2}_{B_{t}}\\
\text{s.t.}\;f_{i}(w_{t})+\langle m_{t},w-w_{t}\rangle\leq f_{i}^{\ast}.\end{gathered}$
(15)
We get the next explicit step:
$w_{t+1}=w_{t}-\frac{f_{i}(w_{t})-f_{i}^{\ast}}{\|m_{t}\|_{B_{t}^{-1}}^{2}}B_{t
}^{-1}m_{t}.$
(16)
In PSPS, the norm in the projection is changed to a weighted norm based on the preconditioning matrix $B_{t}\succ 0$, it helps to improve the convergence rate in case of badly scaled/ill-conditioned datasets.
Gradient regularized Newton method.
One of the main issues of Newton method is a lack of global convergence. To solve it with provably fast convergence, Cubic Regularized Newton method was proposed by Nesterov and Polyak, (). Later, to simplify subproblem solution, the gradient regularization was proposed by Mishchenko, (); Doikov and Nesterov, (). Next, we present a formulation of a Stochastic Cubic Newton Method with gradient regularization from ([10](#S1.E10)).
With $L_{2}$ as a Lipschitz-continuous constant for Hessian, we choose the parameters as follows:
$\displaystyle\tau_{t}=\sqrt{\tfrac{3}{L_{2}\|g_{t}\|}},\,m_{t}=g_{t}=\nabla f_
{i}(w_{t}),$
$\displaystyle D_{t}=H_{t}=\nabla^{2}f_{i}(w_{t}),\qquad B_{t}=\mathit{I},$
and solve the following problem:
$\displaystyle w_{t+1},\alpha_{t+1}=\mathop{\arg\min}_{w\in\mathbb{R}^{d},
\alpha\in\mathbb{R}}\tfrac{1}{2}\|w-w_{t}\|^{2}_{2}+\alpha\sqrt{\tfrac{3}{L_{2
}\|g_{t}\|}}$
(17)
$\displaystyle\text{s.t.}\;f_{i}(w_{t})+\langle g_{t},w-w_{t}\rangle+\tfrac{1}{
2}H_{t}[w-w_{t}]^{2}\leq\alpha.$
We get the next step:
$w_{t+1}=w_{t}-\left[H_{t}+I\sqrt{\tfrac{L_{2}}{3}\|g_{t}\|}\right]^{-1}g_{t}.$
(18)
SP2.
In (Li et al.,, ), the constraint of SPS ([4](#S1.E4)) was extended for the second-order information, aimed at incorporating additional curvature information to accelerate the convergence rate. Next, we present the implicit formulation of SP2 under Assumption [1](#S1):

$$ $\begin{gathered}w_{t+1}=\textstyle{{\arg\min}}_{w\in\mathbb{R}^{d}}\tfrac{1}{2 }\|w-w_{t}\|^{2},\\ \mbox{s.t.}\,f_{i}(w_{t})+\langle g_{t},w-w_{t}\rangle+\tfrac{1}{2}H_{t}[w-w_{ t}]^{2}=0.\end{gathered}$ (19) $$

The explicit formulation was presented only for generalized linear models. In the next sections, we will propose a variant of the explicit solution for SP2 with a connection to Cubic Newton.

### 2.3 Proposed methods

Gradient regularized Newton method with Polyak step-size.
Similarly to SGD and SPS, we propose a new version of Cubic Newton method with Polyak step-size and its stochastic version. If $f_{i}^{\ast}$ is known, for example in the case of interpolation with Assumption [1](#S1), then the method is parameter-free. This result is new both in deterministic and stochastic cases. Similarly to SGD, we fix $\alpha=f_{i}^{\ast}$ in ([17](#S2.E17)) and get the next method.
We choose the parameters as follows:
$\displaystyle\alpha=f_{i}^{\ast},\,\tau_{t}=\sqrt{\tfrac{3}{L_{2}\|g_{t}\|}},
\,m_{t}=g_{t}=\nabla f_{i}(w_{t}),$
$\displaystyle D_{t}=H_{t}=\nabla^{2}f_{i}(w_{t}),\,B_{t}=\mathit{I},$
and solve the following problem:
$\displaystyle w_{t+1}=\textstyle{\arg\min}_{w\in\mathbb{R}^{d}}\tfrac{1}{2}\|w
-w_{t}\|^{2}_{2}$
(20)
$\displaystyle\text{s.t.}\;f_{i}(w_{t})+\langle g_{t},w-w_{t}\rangle+\tfrac{1}{
2}H_{t}[w-w_{t}]^{2}\leq f_{i}^{\ast}.$
The explicit step is formulated as follows:
$w_{t+1}=w_{t}-(1-\kappa_{t})\left[\kappa_{t}\mathit{I}+(1-\kappa_{t})H_{t}
\right]^{-1}g_{t},$
(21)
where $\kappa_{t}=0$ if $f_{i}(w_{t})-f_{i}^{\ast}>\tfrac{1}{2}\|g_{t}\|_{H_{t}^{-1}}^{2}$, otherwise $\kappa_{t}$ is computed by Cubic Newton-type line-search.
Similarly to classical Cubic Newton, the explicit step requires an additional line-search. We present this line-search and provide more details in Appendix [C.1](#A3.SS1). Additionally, we would like to highlight that ([20](#S2.E20)) follows the same formulation as SP2 (Li et al.,, ). Surprisingly, we have discovered that SP2 is a Stochastic Polyak-step size version of the Cubic Newton method.
SANIA Quasi-Newton for $B_{t}\succ 0$.
Similarly to PSPS ([15](#S2.E15)), this approach covers AdaGrad, Adam, Hutchinson’s method, Quasi-Newton methods with $B_{t}\succ 0$, and Newton method for convex functions with $H_{t}\succ 0$. The method is inspired by Affine-Invariant Cubic Newton from (Hanzely et al.,, ). Note, the Hessian approximation $B_{t}$ is used both in scaling of objective norm and in constraint model. We derive it from ([10](#S1.E10)).
The parameters are chosen as follows:
$\displaystyle\alpha=f_{i}^{\ast},\,\tau_{t}=\gamma_{t},\,m_{t}=m_{t},\,D_{t}=B
_{t},\,B_{t}=B_{t},$
and solve the following problem:
$\displaystyle w_{t+1}=\textstyle{\arg\min}_{w\in\mathbb{R}^{d}}\tfrac{1}{2}\|w
-w_{t}\|_{B_{t}}^{2}\qquad\text{s.t.}$
(22)
$\displaystyle f_{i}(w_{t})+\langle m_{t},w-w_{t}\rangle+\tfrac{1}{2}B_{t}[w-w_
{t}]^{2}\leq f_{i}^{\ast}.$
The explicit step is :
$w_{t+1}=w_{t}-\lambda_{t}B_{t}^{-1}m_{t},$
(23)
where for $\upsilon_{t}=\tfrac{2(f_{i}(w_{t})-f_{i}^{\ast})}{\|m_{t}\|_{B_{t}^{-1}}^{2}}$ we define
$\lambda_{t}=\begin{cases}1-\sqrt{1-\upsilon_{t}},&\text{if}\quad\upsilon_{t}
\leq 1,\\
1,&\text{otherwise}.\end{cases}$
(24)
The main difference between PSPS ([16](#S2.E16)) and SANIA-Quasi-Newton ([23](#S2.E23)) is the parameter $\lambda_{t}$. For ([23](#S2.E23)), we can show that $\lambda_{t}\leq 1$ in ([24](#S2.E24)), while in contrast for ([16](#S2.E16)) $\lambda_{t}$ could be much bigger than $1$. For Quasi-Newton and Newton methods, the step-size $\lambda_{t}$ is naturally bounded by $1$, which makes SANIA-Quasi-Newton step-size safer than the step-size of PSPS.
SANIA AdaGrad-SQR. We propose a new preconditioning method, called AdaGrad-SQR, by removing the square root from AdaGrad update. In Section [2.4](#S2.SS4), we will prove that the improved algorithm have ”scale invariance” property. Figure  [2](#S2.F2) shows that the proposed algorithm behaves the same both on original and scaled versions of datasets.
We define $m_{t},B_{t},D_{t}$ for ([23](#S2.E23)) as follows:
$m_{t}=g_{t},\quad\text{and}\quad B_{t}=D_{t}=\operatorname{diag}\left(
\textstyle{\sum_{j=1}^{t}g_{j}^{2}}\right).$
(25)
SANIA Adam-SQR. Along with SANIA AdaGrad-SQR, we propose another ”scale-invariant” method. Following the same idea, it removes the square root from the preconditioning matrix of Adam.
We define $m_{t},B_{t},D_{t}$ for ([23](#S2.E23)) as follows:
$\displaystyle m_{t}$
$\displaystyle=\tfrac{(1-\beta_{1})\sum_{j=1}^{t}\beta_{1}^{t-j}g_{j}}{1-\beta_
{1}^{t}},$
(26)
$\displaystyle B_{t}=D_{t}$
$\displaystyle=\operatorname{diag}\left(\tfrac{(1-\beta_{2})\sum_{j=1}^{t}\beta
_{2}^{t-j}g_{j}^{2}}{1-\beta_{2}^{t}}\right).$
SANIA PCG for Newton method for non-convex functions.
In cases where the functions $f_{i}(w)$ are non-convex, the Hessian matrix $H_{t}$ may not be positive definite but invertible. This characteristic renders the approach not applicable, as $\|g_{t}\|_{H_{t}^{-1}}$ is no longer a norm. To address this issue, we propose a solution based on the rank-$1$ SR-1 approximation.

We present more details in Appendix [4](#algorithm4). In practice, we solve $H_{t}^{-1}g_{t}$ by using Conjugate Gradient method, which allows to compute only Hessian-vector products without computing and storing the full Hessian $H_{t}$.

### 2.4 Affine and Scale Invariance

The family of Stochastic Gradient Methods with Polyak step-size offers an update rule that alleviates the need of fine-tuning the learning rate of an optimizer. However, existing first-order algorithms, whether stochastic or deterministic, perform poorly on ill-conditioned datasets. One possible reason for this is their strong dependence on the chosen basis. This is why, in machine learning, it is common practice to normalize data, as it makes the optimization space and basis more amenable. In the case of generalized linear models (GLM), the choice of basis is directly linked to the handling of ill-conditioned datasets. Changing the basis leads to improvement of conditioning.
Affine invariance is one of the key features of the Newton method, which makes it basis-independent (Nesterov and Nemirovski,, ; Nesterov,, ). Let $A\in\mathbb{R}^{d\times d}$ be a non-degenerate matrix. We consider function $\phi(y)=f(Ay)$. By affine transformation, we denote $f(w)\rightarrow\phi(y)=f(Ay),w\rightarrow A^{-1}y$.
Now, we discuss what is affine invariant friendly and what is not. First of all, the local Hessian norm $\|h\|_{\nabla^{2}f(w)}$ is affine-invariant:

$$ $\displaystyle\|z\|_{\nabla^{2}\phi(y)}^{2}=\left\langle\nabla^{2}\phi(y)z,z \right\rangle=\left\langle A^{\top}\nabla^{2}f(Ay)Az,z\right\rangle=$ $\displaystyle\left\langle\nabla^{2}f(w)h,h\right\rangle=\|h\|^{2}_{\nabla^{2}f (w)}.$ $$

However, the norm $\|z\|_{I}^{2}$ is not affine invariant.
Second of all, Damped Newton method is affine invariant (Lemma 5.1.1 (Nesterov,, )). It means that for the function $f(w)$ Damped Newton method with affine invariant step-size $\gamma_{t}$ generates $w_{t+1}=w_{t}-\gamma_{t}[\nabla^{2}f(w_{t})]^{-1}\nabla f(w_{t})$. For a function $\phi(y)$, Damped Newton method generates $y_{t+1}=y_{t}-\gamma_{t}[\nabla^{2}\phi(y_{t})]^{-1}\nabla\phi(y_{k})$, where $y_{0}=A^{-1}w_{0}$, then $y_{t}=A^{-1}w_{t}$. Essentially, we get a bijection between $y_{t}$ and $w_{t}$. Also, the function values during the optimization are the same $\phi(y_{t})=f(w_{t})$. It means that for GLM, we will automatically get the best basis. Finally, we can show that SANIA Newton and SANIA CG are affine invariant, because the step-size $\lambda_{t}$ in ([24](#S2.E24)) is affine-invariant friendly. All proofs are presented in Appendix [D.2](#A4.SS2.SSS0.Px1).
Scale invariance is a special case of affine invariance, where the matrix $A$ is a diagonal matrix. This implies the removal of rotations from the transformations, allowing only diagonal transformations. To distinguish scale invariance from affine invariance, we denote the transformation $V\in\mathbb{R}^{d\times d}$ as a non-degenerate diagonal matrix. It’s evident that the diagonal preconditioning from AdaGrad, Adam, and Hutchinson is not affine invariant because it does not adapt to rotations. However, they could be scale invariant. It turns out that classical AdaGrad and Adam are not scale invariant, but if we remove the square root, they become scale invariant. We propose the new scale invariant SANIA AdaGrad-SQR in ([25](#S2.E25)) and new scale invariant SANIA Adam-SQR in ([26](#S2.E26)). All proofs are presented in Appendix [D.2](#A4.SS2.SSS0.Px1). Scale invariance property of SANIA Adam-SQR and SANIA AdaGrad-SQR is shown in Figure  [2](#S2.F2), where SANIA Adam-SQR and SANIA AdaGrad-SQR are converging identically for both original and badly scaled versions of the datasets, while using classical Adam and AdaGrad preconditioners result in different convergence steps.

Figure: Figure 2: Comparison of classical and modified Adam and AdaGrad preconditioners used in SANIA optimizing logistic regression for synthetically generated dataset and Mushrooms (scaling factor $k=2$).
Refer to caption: x2.png

Figure  [3](#S2.F3) illustrates that SANIA is able to become scale invariant with various preconditioners. Note that SANIA $\mathit{I}_{d}$, SANIA $(V^{-1})^{2}$, and SANIA $\operatorname{diag}(H^{-1})$ are preconditioned by Identity matrix (i.e. no preconditioning), squared inverse of the scaling vector used to obtain the scaled version of the dataset, and inverse of the Hessian diagonal of the objective function, respectively. One of the most noteworthy observations from this figure is that using the vector employed to transform the dataset for scaling, as a preconditioner, results in a scale invariant method. This essentially leads to convergence in a similar manner as non-preconditioned SANIA applied to the original dataset. In practice, obtaining such information is typically unattainable and often not even approximable. However, by utilizing the curvature of the objective function, we can achieve the same scale invariance property. This is also demonstrated in Figure  [3](#S2.F3) by comparing SANIA preconditioned with the diagonal of the Hessian (SANIA $\operatorname{diag}(H_{t}^{-1})$) on both the original and scaled data. This method results in improved convergence while maintaining scale invariance, albeit with minor numerical instabilities. Nevertheless, SANIA $\operatorname{diag}(H_{t}^{-1})$ is still impractical for large problems involving demanding calculations of Hessian. For reference, in the same figure we display performance of Adam with a constant step size, which deteriorates when scaled data is introduced.

Figure: Figure 3: Observation of scale invariance of SANIA and performance of Adam on 2 datasets with the scaling factor $k=2$ and logistic regression objective function.
Refer to caption: x3.png

Figure: Figure 4: Performance of SANIA and other adaptive methods on 4 datasets (original and badly scaled with scaling factor $k=6$) with logistic regression loss.
Refer to caption: x4.png

## 3 EXPERIMENTS

We test our methods on multiple binary classification datasets from LibSVM (^2^22[https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/](https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/))
, namely colon-cancer, duke and leukemia. Also, we use synthetically generated dataset. To simulate badly scaled data we introduce scaled version of each datasets where its columns are multiplied by a vector $e=\{\exp(a_{i})\}_{i=1}^{d}$ where $a_{i}$ is generated from a uniform distribution on the interval $[-k,k]$. For comparison, we also train Adam, AdaGrad, and Adadelta from PyTorch framework with constant step-sizes. Additional findings and details for the experiments (synthetic dataset generation, learning rates and etc.) can be found in Appendix [E](#A5).

In Figures  [4](#S2.F4) and  [5](#S3.F5) we can see that all presented variations of SANIA outperform other adaptive optimization methods on binary classification task. Once again, note that while other methods require step-size fine-tuning and require multiple runs of experiments, SANIA only requires one run for one set of configurations (i.e. scaling factor, batch-size, and etc.). SANIA finds global minimum of the objective function in both original and scaled versions of datasets.

### 3.1 Convex

Let ${\{(x_{i},y_{i})\}}_{i=1}^{n}$ be our dataset, where $x_{i}\in\mathbb{R}^{d}$ and $y_{i}\in\{-1,+1\}$, then
Logistic regression is defined as
$f_{LogReg}(w)=\tfrac{1}{n}\sum_{i=1}^{n}\log(1+\exp(-y_{i}x_{i}^{T}w)).$

### 3.2 Non-Convex

Let ${\{(x_{i},y_{i})\}}_{i=1}^{n}$ be our dataset, where $x_{i}\in\mathbb{R}^{d}$ and $y_{i}\in\{0,1\}$, then
Non-linear least squares is given by
$f_{NLLSQ}(w)=\tfrac{1}{n}\sum_{i=1}^{n}(y_{i}-\tfrac{1}{(1+\exp(-x_{i}^{T}w))}
)^{2}.$

Figure: Figure 5: Performance of SANIA and other adaptive methods on 2 LibSVM datasets (original and badly scaled with scaling factor $k=6$) with non-linear least squares loss.
Refer to caption: x5.png

## 4 CONCLUSION

In this paper, we introduced a versatile and inclusive framework that not only encompasses classical optimization techniques but also sheds valuable light on Polyak step-size methods. Our research introduce the first Cubic Newton method with Polyak step-size which combines the efficiency of stochastic methods and the robustness of Newton methods. We have presented innovative variants of AdaGrad and Adam optimization algorithms that are scale invariant. Our proposed methods are affine or scale invariant, and this important development ensures the invariance of these methods to basis transformation, expanding their applicability and reliability in various scenarios. Our work is supported by comprehensive experiments including both convex and non-convex settings.

## Appendix A RELATED WORK

Second-order methods have played a crucial role in contemporary optimization since their inception in classical works focused on root-finding algorithms by Newton, (), Raphson, (), Simpson, (), and Bennett, (). Subsequent significant advancements in the Newton method and its local quadratic convergence rates were made by ; ; Kantorovich, (); ; ; Kantorovich, (, ). These methods have been extensively researched, refined, and enhanced in various works, with notable contributions from Moré, (),Griewank, (), Nesterov and Polyak, (). Today, they are widely employed in both industrial and scientific computing. For a comprehensive historical overview of the Newton method, Boris T. Polyak’s paper Polyak, () provides more in-depth insights. Compared to first-order algorithms, second-order methods typically yield faster convergence. However, it’s important to note that the per-iteration computational cost of second-order methods is considerably higher. An example of the classical Newton method can be expressed as follows:

$$ $x_{t+1}=x_{t}-\left[\nabla^{2}f(x_{t})\right]^{-1}\nabla f(x_{t})$ $$

It exhibits quadratic local convergence, but it becomes impractical for large-scale optimization problems due to the necessity of computing the complete Hessian and matrix inversion at each iteration. It also lack of global convergence properties and could diverge if far from the solution.

The Cubic Regularized Newton method by Yurii Nesterov and Boris T. Polyak (Nesterov and Polyak,, ) is one of the main approaches to globalize the Newton method.
This algorithm achieves global convergence with the convergence rate $O(\varepsilon^{-1/2})$ for convex functions.
Nonetheless, a notable limitation of the Cubic Regularized Newton method lies in the auxiliary problem, which typically requires running a separate optimization algorithm to solve it.
Several research papers have proposed regularization techniques based on the gradient norm, aiming to derive an explicit regularized Newton step Polyak, (, ). In Mishchenko, (); Doikov and Nesterov, (), the convergence rate was improved up to $O(\varepsilon^{-1/2})$ for convex functions, under higher assumptions on smoothness it accelerates up to $O(\varepsilon^{-1/3})$ Doikov et al., (). Affine-Invariant Cubic Regularized Newton method with local Hessian norms has the convergence rate $O(\varepsilon^{-1/2})$ and the same subproblem as a classical Newton step Hanzely et al., ().

## Appendix B Proofs

### B.1 Stochastic Gradient Descent with SANIA

Lemma 1

#### Proof.

Denote the Lagrangian as $\mathcal{L}(w,\alpha,\lambda)=f(w)+\tau\alpha+\lambda(g(w)-\alpha)$, where $\lambda\geq 0$ is the Lagrange multiplier. We know that $\frac{\partial\mathcal{L}}{\partial\alpha}=\tau-\lambda$ should be 0 0, which means $\lambda=\tau>0$. According to the complementary slackness, the condition $\lambda(g(w)-\alpha)=0$ should hold. Thus, $\alpha=g(w)$, which means solving problem [27](#A2.E27) is the same as solving problem [28](#A2.E28).

#### Proof.

From Lemma [B.1](#A2.SS1), we know that solving problem [29](#A2.E29) is the same as solving the following problem:

$$ $\displaystyle w_{t+1}=\mathop{\arg\min}\limits_{w\in\mathbb{R}^{d}}\tfrac{1}{2 }\|w-w_{t}\|^{2}_{2}+\gamma_{t}(f_{i}(w_{t})+\langle\nabla f_{i}(w_{t}),w-w_{t }\rangle).$ (31) $$

By taking the derivative of the objective function, we get the solution right away.

### B.2 Stochastic Polyak step-size with SANIA

#### Proof.

Denote the Lagrangian as $\mathcal{L}(w,\lambda)=\tfrac{1}{2}\|w-w_{t}\|^{2}_{2}+\lambda(f_{i}(w_{t})+
\langle\nabla f_{i}(w_{t}),w-w_{t}\rangle-f_{i}^{\ast}),$ and we can get Karush–Kuhn–Tucker(KKT) conditions as below:

$$ $\begin{cases}\frac{\partial\mathcal{L}}{\partial w}=w-w_{t}+\lambda\nabla f_{i }(w_{t})=0\\ \lambda(f_{i}(w_{t})+\langle\nabla f_{i}(w_{t}),w-w_{t}\rangle-f_{i}^{\ast})=0 \\ f_{i}(w_{t})+\langle\nabla f_{i}(w_{t}),w-w_{t}\rangle-f_{i}^{\ast}\leq 0\\ \lambda\geq 0.\end{cases}$ (34) $$

$\lambda\in\mathbb{R}_{+}$ is called Lagrange multiplier, and if $\lambda=0$, then the constrain is not active. We consider these two cases as follows.

$$ $\displaystyle\textit{(i) $\lambda=0$:\,\,}\begin{cases}w_{t+1}=w_{t}\\ f_{i}(w_{t})-f_{i}^{\ast}\leq 0,\textit{It's only true when they are equal.} \end{cases}\textit{(ii) $\lambda>0$:\,\,}\begin{cases}w_{t+1}=w_{t}-\lambda \nabla f_{i}(w_{t})\\ \lambda=\frac{f_{i}(w_{t})-f_{i}^{\ast}}{\|\nabla f_{i}(w_{t})\|^{2}}.\end{cases}$ $$

### B.3 Preconditioned SGD with SANIA

#### Proof.

From Lemma [B.1](#A2.SS1), we know problem [35](#A2.E35) is equivalent to:

$$ $\displaystyle w_{t+1}=\mathop{\arg\min}\limits_{w\in\mathbb{R}^{d}}\tfrac{1}{2 }\|w-w_{t}\|^{2}_{B_{t}}+\gamma_{t}(f_{i}(w_{t})+\langle m_{t},w-w_{t}\rangle).$ (37) $$

Take derivative of $w$ and get solution:

$$ $w_{t+1}=w_{t}-\gamma_{t}B_{t}^{-1}m_{t}.$ (38) $$

By plugging in $m_{t}$ and $B_{t}$, we get formulas for AdaGrad:
$w_{t+1}=w_{t}-\gamma_{t}\frac{g_{t}}{\sqrt{\sum_{j=1}^{t}g_{j}\odot g_{j}}}$,

and for Adam:
$w_{t+1}=w_{t}-\gamma_{t}\frac{\frac{(1-\beta_{1})\sum_{j=1}^{t}\beta_{1}^{t-j}
g_{j}}{1-\beta_{1}^{t}}}{\sqrt{\frac{(1-\beta_{2})\sum_{j=1}^{t}\beta_{2}^{t-j
}g_{j}\odot g_{j}}{1-\beta_{2}^{t}}}}$

### B.4 Hutchinson’s Lemma

#### Proof.

Taking expectation the Hadamard product we have that

$$ $\mathbb{E}[z\odot Hz]=\mathbb{E}[\sum_{i}z_{i}(\sum_{j}H_{ij}z_{j})e_{i}]=\sum _{i}\sum_{j}H_{ij}\mathbb{E}[z_{j}z_{i}]e_{i}.$ (41) $$

Since $\mathbb{E}[z_{j}z_{i}]=I$ we have that $\mathbb{E}[z_{j}z_{i}]=\delta_{ij}=\begin{cases}1\quad\textit{if}\;i=j\\
0\quad\textit{if}\;i\neq j.\end{cases}$

Using the above in [41](#A2.E41) we have that

$$ $\mathbb{E}[z\odot Hz]=\sum_{i}H_{ii}e_{i}$ (42) $$

which is the diagonal of the Hessian matrix.

Let $z$ be a Rademacher random varaible. That is $z_{i}=\begin{cases}1\quad\text{with probability $\frac{1}{2}$}\\
-1\quad\text{with probability $\frac{1}{2}$.}\end{cases}$
Thus for $i,j\in{1,\dots,d}$ and $i\neq j$, we have that $\mathbb{E}[z_{i}]=0$, $\mathbb{E}[z_{i}^{2}]=1$ and $\mathbb{E}[z_{i}z_{j}]=0$.
The same results follow for $z\in\aleph(0,1)$.

## Appendix C Proposed methods

### C.1 Gradient regularized Newton method with Polyak step-size

#### Proof.

For problem ([43](#A3.E43)),
the Lagrangian could be written as follows:

$$ $\displaystyle\mathcal{L}(w,\lambda)=\tfrac{1}{2}\|w-w_{t}\|^{2}_{2}+\lambda \left(f_{i}(w_{t})+\langle g_{t},w-w_{t}\rangle+\tfrac{1}{2}H_{t}[w-w_{t}]^{2} -f_{i}^{\ast}\right).$ $$

Then, we get the next KKT conditions:

$$ $\begin{cases}\frac{\partial\mathcal{L}}{\partial w}=I(w-w_{t})+\lambda\left(g_ {t}+H_{t}(w-w_{t})\right)=0\\ \lambda(f_{i}(w_{t})+\langle g_{t},w-w_{t}\rangle+\tfrac{1}{2}H_{t}[w-w_{t}]^{ 2}-f_{i}^{\ast})=0\\ f_{i}(w_{t})+\langle g_{t},w-w_{t}\rangle+\tfrac{1}{2}H_{t}[w-w_{t}]^{2}-f_{i} ^{\ast}\leq 0\\ \lambda\geq 0.\end{cases}$ $$

Similarly, to previous proofs, the case of inactive constraint with $\lambda=0$ us trivial and we focus on active constraint case.

$$ $\begin{cases}I(w-w_{t})+\lambda\left(g_{t}+H_{t}(w-w_{t})\right)=0\\ f_{i}(w_{t})+\langle g_{t},w-w_{t}\rangle+\tfrac{1}{2}H_{t}[w-w_{t}]^{2}-f_{i} ^{\ast}=0,\\ \lambda>0.\end{cases}$ $$

First, we find $w_{t+1}$ as

$$ $\displaystyle I(w-w_{t})+\lambda\left(g_{t}+H_{t}[w-w_{t}]\right)=0$ $\displaystyle w_{t+1}=w_{t}-\lambda\left[\lambda H_{t}+I\right]^{-1}g_{t}.$ $$

Now, we substitute its new form in the active constraint and get

$$ $\displaystyle f_{i}(w_{t})-f_{i}^{\ast}-\lambda g^{\top}\left[\lambda H_{t}+I \right]^{-1}g_{t}+\tfrac{\lambda}{2}g^{\top}\left[\lambda H_{t}+I\right]^{-1} \lambda H_{t}\left[\lambda H_{t}+I\right]^{-1}g_{t}=0$ $\displaystyle f_{i}(w_{t})-f_{i}^{\ast}-\lambda g^{\top}\left[\lambda H_{t}+I \right]^{-1}g_{t}+\tfrac{\lambda}{2}g^{\top}\left[\lambda H_{t}+I\right]^{-1}( \lambda H_{t}+I)\left[\lambda H_{t}+I\right]^{-1}g_{t}-\tfrac{\lambda}{2}\| \left[\lambda H_{t}+I\right]^{-1}g_{t}\|^{2}_{2}=0$ $\displaystyle f_{i}(w_{t})-f_{i}^{\ast}-\tfrac{\lambda}{2}g^{\top}\left[ \lambda H_{t}+I\right]^{-1}g_{t}-\tfrac{\lambda}{2}\|\left[\lambda H_{t}+I \right]^{-1}g_{t}\|^{2}_{2}=0.$ $$

To simplify the line-search by $\lambda\in[0,+\infty]$, we transform it to $\kappa=\tfrac{1}{1+\lambda}$, which is now $\kappa\in[0,1]$.

$$ $\displaystyle f_{i}(w_{t})-f_{i}^{\ast}-\tfrac{1-\kappa}{2}g^{\top}\left[(1- \kappa)H_{t}+\kappa I\right]^{-1}g_{t}-\tfrac{\kappa(1-\kappa)}{2}\left\|\left [(1-\kappa)H_{t}+\kappa I\right]^{-1}g_{t}\right\|^{2}_{2}=0.$ $$

To simplify the multiple computations of the inverse matrix, one can apply SVD for $H_{t}$ and get the next simplified formulation:

$$ $\displaystyle H_{t}=U_{t}S_{t}U_{t}^{\top}$ $\displaystyle\left[(1-\kappa)H_{t}+\kappa I\right]^{-1}=\left[(1-\kappa)U_{t}S _{t}U_{t}^{\top}+\kappa U_{t}IU_{t}^{\top}\right]^{-1}=U_{t}\left[(1-\kappa)S_ {t}+\kappa I\right]^{-1}U_{t}^{\top}$ $\displaystyle f_{i}(w_{t})-f_{i}^{\ast}-\tfrac{1-\kappa}{2}g^{\top}U_{t}\left[ (1-\kappa)S_{t}+\kappa I\right]^{-1}U_{t}^{\top}g_{t}-\tfrac{\kappa(1-\kappa)} {2}\left\|\left[(1-\kappa)S_{t}+\kappa I\right]^{-1}U_{t}^{\top}g_{t}\right\|^ {2}_{2}=0$ $\displaystyle\tilde{g_{t}}=U_{t}^{\top}g_{t}$ $\displaystyle f_{i}(w_{t})-f_{i}^{\ast}-\tfrac{1-\kappa}{2}\tilde{g_{t}}^{\top }\left[(1-\kappa)S_{t}+\kappa I\right]^{-1}\tilde{g_{t}}-\tfrac{\kappa(1- \kappa)}{2}\left\|\left[(1-\kappa)S_{t}+\kappa I\right]^{-1}\tilde{g_{t}} \right\|^{2}_{2}=0,$ $$

where $S_{t}$ is a diagonal matrix. Note, that this type of line-search is pretty common for Cubic Newton Methods. It adds only additional logarithmic inversions $O(\log\varepsilon^{-1})$ compared to classical Newton.

### C.2 SANIA AdaGrad-SQR for Quasi-Newton.

The following is the explicit implementation of the Quasi-Newton algorithm when choosing AdaGrad-SQR as preconditioning matrix. We add some insurance $\epsilon$ to avoid numerical collapse.

Figure: Algorithm 1 SANIA AdaGrad-SQR

### C.3 SANIA Adam-SQR for Quasi-Newton.

The following is the explicit implementation of the Quasi-Newton algorithm when choosing Adam-SQR as preconditioning matrix. We add some insurance $\epsilon$ to avoid numerical collapse.

Figure: Algorithm 2 SANIA Adam-SQR

### C.4 SANIA PCG for Newton method on convex functions.

For convex setting where Hessian is positive definite, we can choose $B_{t}$ in ([23](#S2.E23)) as Hessian or the approximation matrix of diagonal Hessian. We present detailed algorithm when $B_{t}=\nabla^{2}f_{i}(w_{t})$(we denote as $H_{k}$) as below.

Figure: Algorithm 3 SANIA PCG for convex setting

In practice, we solve this matrix-vector product $(\nabla^{2}f_{i}(w_{t}))^{-1}\nabla f_{i}(w_{t})$
using Conjugate Gradient method. Furthermore, we can incorporate curvature information from Hessian approximation using Hutchinson’s method, Adam or AdaGrad, which allows us to benefit from preconditioned system. In Conjugate Gradient method preconditioning is required to ensure faster convergence and the system can be preconditioned by a matrix $M^{-1}$ that is symmetric and positive-definite. Preconditioned Conjugate Gradient is equivalent to solving the following system:

$$ $\displaystyle E^{-1}\nabla^{2}f_{i}(w_{t})(E^{-1})^{T}E^{T}x=E^{-1}\nabla f_{i }(w_{t}),$ $$

where

$$ $\displaystyle EE^{T}=M.$ $$

If matrix $M_{k}=H_{k}$, then SANIA PCG is affine invariant; if $M_{k}=\operatorname{diag}(H_{k})$, then this method is scale invariant. In experiments you can choose $M_{k}$ as AdaGrad-SQR[25](#S2.E25) or Adam-SQR[26](#S2.E26).

### C.5 SANIA PCG for Newton method on non-convex functions.

For non-convex settings, we cannot use conjugate gradient method to solve this $Hx=g$ (Hessian is not positive definite) linear system of equations anymore. We try to combine Polyak step-size and line searxch Newton-CG method together to get good performance. The following is our specific implementation of the algorithm.

Figure: Algorithm 4 SANIA PCG for Non-convex setting

Since product $B_{t}^{+}\nabla f_{i}(w_{t})$ results in the same direction as $(\nabla^{2}f_{i}(w_{t}))^{-1}\nabla f_{i}(w_{t})$ , and now the algorithm stops once it detects negative curvature, otherwise it still takes CG steps until it hits stopping criteria. You can choose matrix $M_{k}$ to be AdaGrad-SQR[25](#S2.E25) or Adam-SQR[26](#S2.E26) to attain the scale-invariance property and we name them as SANIA PCG AdaGrad-SQR and SANIA PCG Adam-SQR. Notice that the names for the convex and non-convex setting are the same, but the implementation of these methods are slightly different due to the effectiveness of conjugate gradient methods.

## Appendix D Affine and Scale Invariance

### D.1 Affine Invariance

#### Proof.

Let $y_{k}=B^{-1}x_{k}$ for some $k\geq 0$. Then

$$ $\displaystyle y_{k+1}=$ $\displaystyle y_{k}-[\nabla^{2}\phi(y_{k})]^{-1}\nabla\phi(y_{k})=y_{k}-[B^{T} \nabla^{2}f(By_{k})B]^{-1}B^{T}\nabla f(By_{k})$ $\displaystyle=$ $\displaystyle B^{-1}x_{k}-B^{-1}[\nabla^{2}f(x_{k})]^{-1}\nabla f(x_{k})=B^{-1 }x_{k+1}.\quad\square$ $$

Thus, the Newton’s method is affine invariant with respect to affine transformations of variables.

#### Proof.

We define $\lambda_{k}=\begin{cases}1-\sqrt{1-\upsilon_{k}},&\text{if}\quad\upsilon_{k}
\leq 1,\\
1,&\text{otherwise},\end{cases}$ where $\upsilon_{k}=\tfrac{2(f_{i}(x_{k})-f_{i}^{\ast})}{\|\nabla f_{i}(x_{k})\|_{{
\nabla^{2}f(x_{k})}^{-1}}^{2}}$, and for $\hat{\lambda_{k}}$, $\hat{\upsilon_{k}}=\tfrac{2(\phi_{i}(y_{k})-\phi_{i}^{\ast})}{\|\nabla\phi_{i}
(y_{k})\|_{{\nabla^{2}\phi(y_{k})}^{-1}}^{2}}$. Let $y_{k}=B^{-1}x_{k}$ for some $k\geq 0$. We have this condition $\hat{\upsilon_{k}}=\tfrac{2(\phi_{i}(y_{k})-\phi_{i}^{\ast})}{\|\nabla\phi_{i}
(y_{k})\|_{{\nabla^{2}\phi(y_{k})}^{-1}}^{2}}=\tfrac{2(f_{i}(By_{k})-f_{i}^{
\ast})}{\|B^{T}\nabla f_{i}(By_{k})\|_{[B^{T}\nabla^{2}f(By_{k})B]^{-1}}^{2}}=
\tfrac{2(f_{i}(x_{k})-f_{i}^{\ast})}{\|\nabla f_{i}(x_{k})\|_{{\nabla^{2}f(x_{
k})}^{-1}}^{2}}=\upsilon_{k}$ holds, which means $\hat{\lambda_{k}}=\lambda_{k}$.

Then

$$ $\displaystyle y_{k+1}=$ $\displaystyle y_{k}-\lambda_{k}[\nabla^{2}\phi(y_{k})]^{-1}\nabla\phi(y_{k})=y _{k}-\lambda_{k}[B^{T}\nabla^{2}f(By_{k})B]^{-1}B^{T}\nabla f(By_{k})$ $\displaystyle=$ $\displaystyle B^{-1}x_{k}-\lambda_{k}B^{-1}[\nabla^{2}f(x_{k})]^{-1}\nabla f(x _{k})=B^{-1}x_{k+1}.\quad$ $$

Thus, the SANIA Newton method is affine invariant with respect to affine transformations of variables. $\square$

### D.2 Scale Invariance

Kempka et al., (); Zhuang et al., () illustrate this important but overlooked property of an optimization algorithm. It is widely recognized that the convergence rate of minimizing a twice continuously differentiable function $f$ through a first-order optimization algorithm depends heavily on the condition number. To mitigate the impact of the condition number, one effective approach is the use of preconditioners relying on Hessian of the function which yields affine invaraince. Consider the Hessian cannot be easily estimated, Zhuang et al., () shows that scale invariance gives similar advantages to the use of an optimal diagonal preconditioner.

They also showed why algorithms like SGD and Adam have such excellent performances in DNNs even though they are not scale invariant. Because they are intensively linked to the use of batch normalization which normalizes the gradients. Without BN, using SGD with momentum and weight decay, even a tiny learning rate will lead to divergence while training a deep neural network. But for the upgraded version of Adam–AdamW which enjoys scale invariance outperforms Adam when both are finely tuned.

Now, we will show the classical AdaGrad and Adam are not scale invariant but AdaGrad-SQR and Adam-SQR.

#### Proof.

We define $\lambda_{k}=\begin{cases}1-\sqrt{1-\upsilon_{k}},&\text{if}\quad\upsilon_{k}
\leq 1,\\
1,&\text{otherwise},\end{cases}$ where $\upsilon_{k}=\tfrac{2(f_{i}(x_{k})-f_{i}^{\ast})}{\|m_{k}\|_{B_{k}^{-1}}^{2}}$, and for $\hat{\lambda_{k}}$, $\hat{\upsilon_{k}}=\tfrac{2(\phi_{i}(y_{k})-\phi_{i}^{\ast})}{\|\hat{m_{k}}\|_
{\hat{B_{k}}^{-1}}^{2}}$.

Let $y_{k}=V^{-1}x_{k}$ for some $k\geq 0$. We have

$$ $\begin{cases}\hat{B_{k}}=\sum_{j=1}^{k}\nabla\phi_{i_{j}}(y_{k})^{2}=\sum_{j=1 }^{k}[V^{T}\nabla f_{i_{j}}(Vy_{k})]^{2}=V^{T}[\sum_{j=1}^{k}\nabla f_{i_{j}}( x_{k})^{2}]V=V^{T}B_{k}V,\\ \hat{m_{k}}=\nabla\phi_{i_{j}}(y_{k})=V^{T}\nabla f_{i_{j}}(Vy_{k})=V^{T} \nabla f_{i_{j}}(x_{k})=V^{T}m_{k},\\ \hat{\upsilon_{k}}=\tfrac{2(\phi_{i}(y_{k})-\phi_{i}^{\ast})}{\|\hat{m_{k}}\|_ {\hat{B_{k}}^{-1}}^{2}}=\tfrac{2(f_{i}(Vy_{k})-f_{i}^{\ast})}{\|V^{T}m_{k}\|_{ (V^{T}B_{k}V)^{-1}}^{2}}=\tfrac{2(f_{i}(x_{k})-f_{i}^{\ast})}{\|m_{k}\|_{B_{k} ^{-1}}^{2}}=\upsilon_{k}.\end{cases}$ $$

Then

$$ $\displaystyle y_{k+1}=$ $\displaystyle y_{k}-\hat{\lambda_{k}}\hat{B_{k}}^{-1}\hat{m_{k}}=y_{k}-\lambda _{k}[V^{T}B_{k}V]^{-1}V^{T}m_{k}$ $\displaystyle=$ $\displaystyle V^{-1}x_{k}-\lambda_{k}V^{-1}B_{k}^{-1}m_{k}=V^{-1}x_{k+1}.\quad$ $$

Thus, the AdaGrad-SQR method is scale invariant. $\square$
And for Adam-SQR setting where $m_{k}=\frac{(1-\beta_{1})\sum_{j=1}^{k}\beta_{1}^{k-j}\nabla f_{i_{j}}(x_{k})}
{1-\beta_{1}^{k}},B_{k}=\frac{(1-\beta_{2})\sum_{j=1}^{k}\beta_{2}^{k-j}\nabla
f
_{i_{j}}(x_{k})^{2}}{1-\beta_{2}^{k}}$, and $\hat{m_{k}}=\frac{(1-\beta_{1})\sum_{j=1}^{k}\beta_{1}^{k-j}\nabla\phi_{i_{j}}
(y_{k})}{1-\beta_{1}^{k}},\hat{B_{k}}=\frac{(1-\beta_{2})\sum_{j=1}^{k}\beta_{
2}^{k-j}\nabla\phi_{i_{j}}(y_{k})^{2}}{1-\beta_{2}^{k}}$. Similarly, we can get $\hat{B_{k}}=V^{T}B_{k}V,\;\hat{m_{k}}=V^{T}m_{k}.$ Rest proofs are the same. From proofs above we can know for simple AdaGrad and Adam they are not scale invariant , because $\hat{B_{k}}=V^{T}B_{k}\neq V^{T}B_{k}V$.

### D.3 GLM

Suppose $f_{i}$ is the loss over a linear model with

$$ $f_{i}(w)=\psi_{i}(x_{i}^{T}w-y_{i}),$ (50) $$

where $\psi_{i}:\mathbb{R}\rightarrow\mathbb{R}$ is the loss function, and $x_{i}$ is the $i^{th}$ data and $y_{i}$ is the corresponding label. Let the sequence $\{w_{k}\}$ be generated by method as applied to the function $f$. Consider the sequence $\{\hat{w_{k}}\}$, generated by the same method but for function $\phi$
where $\phi(\hat{w_{k}})=f(B\hat{w_{k}})=\psi_{i}(x_{i}^{T}B\hat{w_{k}}-y_{i})$.

Take $x_{i}^{T}B$ as a whole, it can be seen as we are doing linear transformation to the data. When matrix $B$ is badly scaled, it will lead to a ill-conditioning dataset. And it inhibits the performance of the general algorithms, which is specifically reflected in the need for more iterations to converge, or even diverge on the worst case. But if the algorithm enjoys affine invariant property, that is, $\hat{w_{k}}=B^{-1}x_{k}$. Then we have $\psi_{i}(x_{i}^{T}B\hat{w_{k}}-y_{i})=\psi_{i}(x_{i}^{T}BB^{-1}x_{k}-y_{i})=f_
{i}(w)$, which means we automatically have the same function value as the original one as every iteration goes.

## Appendix E Additional Experiments and Details

All experiments were run with $5$ different seeds $(0,1,2,3,4)$ using PyTorch 2.0.1+cu118 on a computing machine with AMD EPYC 7402 24-Core Processor with 2.8GHz of base clock. Default datatype in PyTorch is set to
torch.float64. Code used for the experiments is publicly available (^3^33[https://github.com/fxrshed/SANIA](https://github.com/fxrshed/SANIA)).

### E.1 Datasets

Table [1](#A5.T1) refers to dataset and training specifications used throughout all experiments. In order to simulate badly scaled datasets we use scaling precedure shown in ([51](#A5.E51)).

**Table 1: Description and setting of datasets used for experiments.**
| Dataset | # of datapoints | # of features | Batch-size |
| --- | --- | --- | --- |
| colon-cancer | 62 | 2000 | 16 |
| duke breast-cancer | 44 | 7129 | 16 |
| leukemia | 38 | 7129 | 16 |
| synthetic | 1000 | 1000 | 200 |
| mushrooms | 8124 | 112 | 256 |

$$ $A^{n\times d}=\begin{pmatrix}a_{1,1}&a_{1,2}&\dots&a_{1,d}\\ a_{2,1}&a_{2,2}&\dots&a_{2,d}\\ \vdots&\ddots&\ddots&\vdots\\ a_{n,1}&a_{n,2}&\dots&a_{n,d}\\ \end{pmatrix}\quad\xrightarrow[]{scale}\quad\hat{A}^{n\times d}=\begin{pmatrix }a_{1,1}\times v_{1}&a_{1,2}\times v_{2}&\dots&a_{1,d}\times v_{d}\\ a_{2,1}\times v_{1}&a_{2,2}\times v_{2}&\dots&a_{2,d}\times v_{d}\\ \vdots&\ddots&\ddots&\vdots\\ a_{n,1}\times v_{1}&a_{n,2}\times v_{2}&\dots&a_{n,d}\times v_{d}\\ \end{pmatrix},$ (51) $$

where $v_{i}=e^{b_{j}},\quad b_{j}\in\mathrm{Uniform}[-k,k]$.

Learning rates of algorithms used for experiments are not chosen randomly. To avoid overoptimized learning rates obtained using special algorithms and at the same time to adhere to some fairness of the results we conducted experiments with a series of learning rates $\gamma=2^{n}$ where $n\in range(-2,-16,2)$. Next, we used the best performing step size as the main result for certain optimizer. All of them are presented in Table  [2](#A5.T2).

**Table 2: Learning rates of Adam, AdaGrad and Adadelta.**
| Loss | Dataset | $\gamma$ Adam | $\gamma$ AdaGrad | $\gamma$ Adadelta |
| --- | --- | --- | --- | --- |
| Logistic Regression | colon-cancer | $2^{-10}$ | $2^{-6}$ | $2^{-2}$ |
|  | duke breast-cancer | $2^{-12}$ | $2^{-8}$ | $2^{-2}$ |
|  | leukemia | $2^{-12}$ | $2^{-8}$ | $2^{-2}$ |
|  | synthetic | $2^{-6}$ | $2^{-6}$ | $2^{-4}$ |
| Non-linear least squares | colon-cancer | $2^{-10}$ | $2^{-6}$ | $2^{-2}$ |
|  | duke breast-cancer | $2^{-14}$ | $2^{-12}$ | $2^{-4}$ |
|  | leukemia | $2^{-14}$ | $2^{-2}$ | $2^{-2}$ |
|  | synthetic | $2^{-6}$ | $2^{-6}$ | $2^{-4}$ |

### E.2 Synthetic Dataset

import numpy as np
def generate_synthetic(n, d):
data = np.random.randn(n, d)
w_star = np.random.randn(d) # optimal point

target = data @ w_star
target[target ¡ 0.0] = -1.0
target[target ¿ 0.0] = 1.0

return data, target

### E.3 More findings

In Figure  [6](#A5.F6) we can see that proposed SANIA CG and SP2 for Generalized Linear Models presented in Li et al., () generate identical steps towards the minimum given the exact same set of observations $x_{i}$. However, disadvantage of SP2 in this case is that it has a closed form solution only for GLMs.

Figure: Figure 6: Performance of SANIA and SP2 GLM on 2 datasets (batch size $=1$) with logistic regression.
Refer to caption: x6.png

Figure  [7](#A5.F7) shows that unlike other classical adaptive methods, SANIA with Newton step is scaling invariant. The same behaviour can be observed in Figure  [8](#A5.F8) where SANIA AdaGrad-SQR is not only scaling invariant but also displays significantly better performance compared to Adam, AdaGrad and Adadelta with a constant learning rate.

Figure: Figure 7: SANIA Newton compared to other adaptive methods on original and badly scaled ($k=5$) synthetic binary classification dataset (batch size $=100$) with logistic regression objective function.
Refer to caption: x7.png

Figure: Figure 8: SANIA AdaGrad-SQR compared to other adaptive methods on original and badly scaled ($k=10$) mushrooms dataset (batch size $=256$) with non-linear least squares objective function.
Refer to caption: x8.png

### E.4 Experiments with Cubic Newton with Polyak step-size

In this subsection, we present results for Cubic Newton with Polyak step-size from ([21](#S2.E21)). In Figure [9](#A5.F9), we compare classical Cubic Newton from (Nesterov and Polyak,, ), Gradient Regularized Newton from (Mishchenko,, ; Doikov and Nesterov,, ) and our Cubic Newton with Polyak step-size on full-batch logistic regression with $\tfrac{\mu}{2}\|w\|_{2}^{2}$-regularization, where $\mu=1e-4$. To show globalization properties, we choose the starting point far from the solution $x_{0}=3e$, where $e$ is a vector of all ones. We present Cubic Newton with theoretical parameter $L_{2}=0.1$, with fine-tuned parameter $L_{2}=0.0004$; Gradient Regularized Newton with fine-tuned parameter $L_{2}=0.0004$. There is a huge difference between fine-tuned and theoretical choice. It means that the method is pretty sensitive to the choice of the parameter $L_{2}$. For Cubic Newton with Polyak step-size, we denote approximate $f^{\ast}$ as $\hat{f}$. Then, we present the precise approximation $\hat{f}=f^{\ast}=0.3361$, close lower approximation $\hat{f}=0.3$, and the very simple and naive lower bound $\hat{f}=0$. For all three cases, the convergence is almost the same. It also shows that Cubic Newton with Polyak step-size is very robust to the parameter $\hat{f}$, where even the most naive choice works perfectly fine. Finally, we highlight that Cubic Newton with Polyak step-size significantly overperform other Cubic methods even with fine-tuned parameters.

Figure: Figure 9: Gradient regularized(Cubic) Newton with Polyak step-size vs Cubic Newton methods for $\tfrac{\mu}{2}\|w\|^{2}$-regularized logistic regression
Refer to caption: x9.png

## Peer Reviews

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: Could the authors provide a formal link between scale invariance and using SANIA? It seems that there is some overlap between the two, but SANIA does not ensure scale invariance by itself. It would be interesting for potential users of SANIA to provide conditions for their algorithm to be scale-invariant.

Addition to Fig. 5: how does SANIA compare to concurrent methods in terms of training loss?
soundness: 3
strengths: # Originality

As far as I know, the work presented in this paper is original.


# Clarity

Overall, the paper is easy to follow.


# Quality

I am very grateful to the authors for providing additional experiments (Appendix E), which show not only test accuracy, but also test loss. I would have been even better to provide the training loss, since the proposed algorithms have been designed to improve the optimization process (not the generalization).

Overall, the paper seems to be correct.

# Significance

This paper provides practical uses of *SP2: A Second Order Stochastic Polyak*, Li et al., 2023. 
So the present paper is relevant for the community.
summary: This paper proposes a training algorithm generalizing the "Polyak step-size" to second-order algorithms, such as cubic Newton or quasi-Newton. 

Let $f$ be a function to minimize and $w\_t$ the current estimate of the argmin.
For an order-1 method, the Polyak step-size comes down to choosing the step-size such that we would attain the minimum $w^*$ in one step if $f$ was affine between $w\_t$ and $w^*$. For order-2 methods, the authors distinguish the metric used for the parameter space, denoted by $B_t$, and the metric related to the curvature of the local approximation of $f$, represented by a matrix $D_t$.

Then, the authors propose variations of existing algorithms based on their generalization of Polyak step-size.
weaknesses: # Originality

This paper can be seen as a practical application of *SP2: A Second Order Stochastic Polyak*, Li et al., 2023.


# Clarity

At the end of Section 2, the authors use the notations: "SANIA $I_d$", "SANIA $(V^{-1})^2$" and "SANIA $\mathrm{diag}(H^{-1})$", while writing that there is no preconditioning. Thus, I understand that "SANIA $A$" refers to Eqn. (6) with $D_t = A$. But there is no formal definition of "SANIA $A$". It should appear somewhere.

Just above Eqn. (9), one can read: $m_t = m_t$, which should be "$m_t = g_t$" (?).

Several mistakes, Figure 1, first plot: 
 * the legend is difficult to understand, since several curves have the same label;
 * some curves do not seem to correspond to their label: apparently, light green should be "SANIA $(V^{-1})^2$";
 * the purple curve is dotted, while it is not dotted in the legend.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: What are your thoughts on how learning rate schedules such as cosine decay etc compose with $\lambda_t$ schedule defined in (19)?

How does $\lambda_t$ change during training?
soundness: 3
strengths: 1. Thorough theoretical analysis.
2. Clear writing
summary: This paper introduces cubic Newton method with adaptive Polyak step-size, enhancing robustness and efficiency in non-convex tasks. It also proposes scale-invariant variants of AdaGrad and Adam, which improve optimizer performance on poorly scaled data. Extensive experiments validate the effectiveness of these methods across diverse optimization scenarios.
weaknesses: 1. It would be great to have theoretical/practical comparisons with Sophia [1], which also uses Hutchinson's based method to compute their preconditioner.
2. KATE [2] removes square root to ensure scale invariance of Adagrad. It would be great to have theoretical/empirical comparison with KATE. 

I am willing to improve the score if the above are addressed.

[1] Liu, Hong, et al. "Sophia: A scalable stochastic second-order optimizer for language model pre-training." arXiv preprint arXiv:2305.14342 (2023).
[2] Choudhury, Sayantan, et al. "Remove that Square Root: A New Efficient Scale-Invariant Version of AdaGrad." arXiv preprint arXiv:2403.02648 (2024).

---

confidence: 5
contribution: 3
initial_rating: 6
presentation: 4
questions: 1. Instability of the Euler method and noises are believed to be helpful for generalization (see Ref[1]). Relatively large step size is helpful for the training process to jump out of bad local minima. In the new framework, when the steps are spontaneously obtained, the instability are also almost inhibited spontaneously. This may be a severe problem for the new frame work. 

2. From the numerical results, it seems that the test accuracy of Sania is slightly worse than existing methods. The author should also compare the results with SGD, since it usually has competitive generalization ability.

3. In Line 209, what is the meaning of "Otherwise, we replaced step-size parameter γt to parameter fi∗?". The authors should explain more clearly.

4. Due to the anisotropic property of the loss landscape, the estimate f* may not be a good guess for the local optimization. How will this affect the performance of the new method?

[1] The Implicit Regularization of Dynamical Stability in Stochastic Gradient Descent, Lei Wu, Weijie J. Su, ICML 2023;
soundness: 3
strengths: The new method can be naturally incorporated with existing methods and the realization is not difficult.
summary: The work tries to incorporate a Polyak-type optimization framework with existing optimizers for machine learning. In the new method, the optimization direction and step size are spontaneously determined by a simple local optimization for each step. The new method is not difficult in realization and appears to work. From the numerical results, the test accuracy of network trained with the new method is comparable to existing methods (but not better than them).
weaknesses: From the numerical methods, the new method does not show improvements compared to previous methods.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: In Page 2, SPS has been derived, why not just cite it? The difference between g_t and m_t is unclear to me, these seem to often mean the same thing? Why not keep the notation consistent? In (6) you have written \| w \|_{B_t} is a Euclidean norm, would it be better to say " \| \cdot \|_{B_t} is a Euclidean norm"?
soundness: 3
strengths: Provides a novel formulation of stochastic optimization algorithms and a novel method that does not require fine-tuning of the learning rate parameter. Investigation into scale invariance is novel. Theorems, equations and ideas are presented clearly. Numerical experiments show that the SANIA methodology achieves similar performance to algorithms with fine-tuned learning rates, and improved robustness as the training curves fluctuate less.
summary: This paper presents SANIA, which is a method that doesn't require manual fine-tuning of the learning rate in commonly used stochastic optimization algorithms, leading to faster optimization. The authors present a framework which generalises common stochastic optimization algorithms. The authors also consider affine and scale invariance which seeks to address poorly scaled or ill-conditioned problems. The authors compare performance of SANIA with commonly used algorithms (which have had fine-tuned parameters) for training classifiers for MNIST, FashionMNIST, CIFAR10 and SVHN.
weaknesses: Numerical experiments lack a comparison to other options for tuning-free methods.

Minor comments:

Occasional incorrect grammar, including after equation (7): "This leads us to Stochastic Polyak step-size method." should read "This leads us to the Stochastic Polyak step-size method." Also bad grammar in the statement of Theorem 1. "Another way to derive this formulation is by solving" (4), I think it may be better to reference appendix B.2 to make it clear why this holds. "interpolation condition" isn't very clearly defined in my opinion, it may be better to be more clear. "in practice it displays better convergence to the true Hessian than other
similar methods like BFGS" - I think this needs a reference.

Meta-review: This paper studies a general framework for Preconditioned and Second-order Polyak methods. In particular, the paper proposes a first Stochastic Cubic Newton method with polyak step-size and also introduces the new scale invariant versions of AdaGrad and Adam, which make them invariant to some basis transformations. The paper also provides a few experiments to demonstrate the performance of the algorithm. The improvements (if any) are not very significant but the authors highlight that the key benefit of this work is need for  limited tuning of hyperparameters.

The reviews for the paper are very much on the borderline. The reviewers while acknowledging that the work is interesting also highlighted the weakness of lack of comparison with schedule free optimization and somewhat limited empirical evidence. After reading through the paper, I was extremely disappointed with the empirical evidence and I believe these weakness are very important to address before publication. The paper's core claim is about the eliminating the need for manual stepsize hyperparameter setting, however the experiments are done on extremely small settings (where hyperparameter tuning is of not much value) and lack proper comparisons. I think this paper would greatly benefit from thorough empirical analysis on larger settings (e.g. ImageNet, ResNet, VIT). I recommend rejection in the current form.

Decision: Reject