A Multilevel Low-Rank Newton Method with Super-linear Convergence Rate and its Application to Non-convex Problems
Authors: Nick Tsipinakis
Panagiotis Tigas
Panos Parpas
Venue: Accepted by TMLR ()

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

Published in Transactions on Machine Learning Research (01/2026)
A Multilevel Low-Rank Newton Method with Super-linear Convergence Rate and its Application to Non-convex Prob- lems
Nick Tsipinakis
nikolaos. tsspinakkssuunddittaee...
Department of Mathematics and Computer Science
UniDistance Suisse Brig, Switzerland
Panagiotis Tigas
ptigas@robots.or.co
Department of Compuler Science Oxford University Oxford, UK
Panos Parpas
panospparpo@@@mpprial..ccc uk
Department of Computing Imperial College London London, UK
Reviewed on OpenReview: https: // openreview. net/ f id= PKakPeVVja
Abstract
Second-order methods can address the shortcomings of first-order methods for the optimiza- tion of large-scale machine learning models. However, second-order methods have signifi- cantly higher computational costs associated with the computation of second-order infor- mation. Subspace methods that are based on randomization have addressed some of these computational costs as they compute search directions in lower dimensions. Even though super-linear convergence rates have been empirically observed, has not been possible to rigorously show that these variants of second-order methods can indeed achieve such fast rates. Also, it is not clear whether subspace methods efficient for non-convex settings. To address these shortcomings, we develop a link between multigrid optimization methods and low-rank Newton methods that enables us to prove the super-linear rates of stochastic low-rank Newton methods rigorously. Our method does not require any computations in the original model dimension. We further propose truncated version of the method that capable of solving high-dimensional non-convex problems. Preliminary numerical experi- ments show that method has better escape rate from saddle points compared to the state-of-thee ahr first-order methods and thus returns lower training errors.
1
Introduction
When it comes to applying second-order optimization methods in machine learning, there two open ques- tions: 1. Can second-order methods be implemented efficiently? 2. Can second-order methods outperform standard first-order methods for traditional ML metrics such as generalization error? Regarding the latter question, several recent. articles have argued the elficacy of second-order methods in deep learning (Singh & Alistarh 2020: Pascanu & Bengio, 2013 Dauphin et al.) 2014), reinforcement learning (Wu et al., 2017] and variational inference (Regier results obtained via second-order methods, the first question regarding the efficiency of these methods has become more urgent. In this paper we provide some answers to the first question.
2017) (to name but a few). Given the recent promising experimental

---

Published in Transactions on Machine Learning Research (01/2026)
The development of randomized methods such as sketching and sub-sampling has enabled the application of second-order methods to large scale ML. models. Super-linear or linear-quadratic (composite) convergence for sub-sampled methods has been shown in various recent works (Berahas et al., 2017 Bollapragada et al.
2019 Erdogdu & Montanari 2015 Pilanci & Wainwright 2017 Khoraa-Ki & Mahoney 2019 Na
jet al. 2023). But all the existing methods require O(n") operations to compute the Newlon direction. Stochastic multilevel or sub-space methods, such as the one proposed in this paper, have ((nN) cost for computing and storing the Hessian, and O(n²N) cost for computing the Newton direction (where N < n is the dimension of the sub-space). However, it is unclear if stochastic multilevel methods can still converge super-linearly when the hessian is approximated via a random method. Although it is clear from numerical experiments that randomized Newton methods can obtain super-linear convergence, rigorous proof under general conditions has not appeared before Hanzely et al. 2020); Gower el al. (2019)
In terms of the thcoretical contributions of this paper, the algorithms of Pilanci & Wainwright (017 and Tsipinakis & Parpas| (2024) are particularly relevant since they are analyzed using self-concordant functions. In Pilancii & Wainwright/ (2017), restrictiiee assumptions are made regarding the Hessian. In particular, the authors assume that the square root of the Hessian is available or easily computable, and the method requires access to the full Hessian and both the proof and hence the algorithm cannot be extended to the non-convex case. In this paper we relax all these assumptions without compromising the convergence rate. In Tsipinakis & Parpas (2024), super-linear convergence was established, but the algorithm cannot be applied to the non-convex case without non-trivial modifications. This paper links multilevel optimization methods and randomized Newton methods. We exploit this link and establish the global convergence of the algorithm and local super-linear convergence rate for self-concordant functions. Because of the way the low-rank Hessian approximation is performed, our method can casily be extended to the non-convex case and to very high dimensions (even when the Hessian is dense).
Regarding the stochastic subspace methods, the work of Gower et al. 2019 is quite relevant, however the method only achieves a lincar rate, whereas we achiove a super-linear rate for the same iteration complexity. In addition, the method in Gower et al. (2019) cannot be exlended to the non-convex case without non-trivial modifications. In Hanzely et al. (2020). the authors develop a related method that is based on the cubic Newton method. But the method has lincar rate, and it is unclear if the method can be efficiently applicd to non-convex problems. Shang et al. (024 develop an accelerated double-sketching subspace Newton method that also attains a local linear rate. Furthermore, in Fuji et al. (2025) and Cartis et al. (2025), the authors establish global convergence for subspace-reguularzeed Newton methods. However, fast local convergence rates are not available, and these methods have only been tested on relatively small non-convex problems.
Here, order to apply the Newton method to non-convex problems, we take a truncated SVD of the Hessian to keep the N + most informative eigenvalues and replace the rest with the (N + 1)th eigenvalue while negative cigenvalues are replaced by their absolute value. Therefore, as in the convex case, enforcing the Hessian matrix to be positive definite and premultiplying the negative gradient by its inverse, we perform local change of coordinates and thus we should expect an accelerated convergence behavior compared to the first-order methods. In addition, when dealing with non-convex functions, the fact that we do not allow for sufficiently small eigenvalues means that slow and flat manifolds around saddle points will turn into saddles whose eigenvalues are large and we therefore achieve a faster escape rate from the unstable region of the saddle points. Similar approaches when computing the Newton direction have alrcady been explored (O'Leary-Ruseberr
al.. 2019 Paternain et. al., 2019 Erdogdu & Montanari 2015). However, they all require the computation of the full Hessian matrix and hence they are inefficient for high dimensional problems. A work that scales to high dimensional problems which also approximates the Hessian matrix was explored by Marteau-Fert Ferey et al. (2019). The difference to our approach is that it regularizes the Hessian matrix instead of performing a Truncated SVD. However, selecting the regularization parameter such that the Hessian matrix becomes positive definite is a difficult task and thus the method is suitable to convex problems. In addition, their method has limited applicability since it relies on the existence of an lz regularizer to work. Such an assumption is not present. in this paper.
To this end, all the aforementioned advantages of our approach are demonstrated in our numerical experi- ments. We show that our method is efficient in terms of computational requirements and it can be applied to models with millions of parameters even when the Hessian is dense. In particular, the experiments suggest
2

---

Published in Transactions on Machine Learning Research (01/2026)
that the proposed methods are well-suited for problems with low-rank Hessian matrix an assumption that is satisfied in most practical problems, especially in non-convex settings. Our numerical experiments suggest that the proposed method behaves similarly to the cubic Newton method (in terms of efficiency and ability to escape saddle points) but with significantly less computational requirements. Morcover, we apply the method to the minimization of a deep auloencoder model, which is known to exhibit numerous flat regions and saddle points. The results demonstrate that the proposed method escapes saddle points and flat regions significantly more efficiently than Adam, despite constructing its iterates in a much smaller subspace.
2 Background and Method
In this section, we briefly discuss the main components of multilevel optimization methods, and introduce the coarse-grained model, which we later modify to compute search directions. Tn the optimization literature, multilevel methods are also referred to as multigrid methods. Since there is notion of "grid" in the proposed algorithm we use the term multilevel, but in the context of this paper the two terms are equivalent. Multilevel optimization methods have two components: (1) A hierarchy of coarse models, and (2) Linear prolongation and restriction operators used to "transfer" information up and down the hicrarchy of models. In this section, we describe these two components; however, for simplicity, we consider a two-level hierarchy. Our results can be extended to more than two levels without substantial changes.
2.1 The Coarse- poorse-grained Model and Main Assumptions
We first develop the proposed methodology for a convex optimization model. The non-convex case is dis- cussed in Scction 3.1.1 Consider the following model,
x
arg min f(x). xER"
(1)
where we assume that the function f is strictly convex and colforancor..nt. Moreover, it has a closed sub- level set and is bounded below that x* exists. The definition and properties of self-concordant functions can be found in appendix A. and we refer the interested reader to Boyd & Vandenberghe (2004); : Nesterov (2004); Nesterov et al. (2018); ; Nesterov & Nemirovskil (1944) for additional background information. We will adopt the terminology used in the multilevel literature and call the model in equation the fine model. We will also assume that we can construct a lower dimensional model called the coarse model. In traditional multilevel methods the coarse model is typically derived by varying a discretization parameter. In machine learning applications coarse model can be derived by varying the number of pixels in image processing applications (Galun et al., 2015). or by varying the dictionary size and fidelity in statistical pattern recogni- tion (see c.g. Hovhannisyan et al. (2016) for examples in face recognition, and background extraction from video). We denote the coarse model as,
arg min F(y). YERN
(2)
We assume that N < n and typically N < and that F also a strictly convex, self-concordant function. The property that N < justifies the use of the terms coarse-grained reduced order model. The role of the coarse model is to help generate high quality scarch directions from the current incumbent point but at reduced compulational cost. Let R& E N be the restriction operulor with which one may transfer information from the fine to coarse model and thus we define the initial point in the coarse model as yo == R&X; at itcration k. For simplification purposes from now onwards we will omit the subscript k from the restriction operator. In order for the search direction to be useful (e.g. a descent direction), we assume that first-order coherency condition below holds,
RV = VF
See Wen & Goldfarb (2009) for a discussion of the first-order coherency condition multilevel optimization. In addition to the restriction operator we also assume that a prolongation operator, P Rx is available. We make the following assumption about these two operators.
3

---

Published in Transactions on Machine Learning Research (01/2026)
Assumption 1. The restriction and prolongation operators R and P satisfy, P= =RR,, where > 0, and P has full column rank, Le.,, rank N.
The assumptions regarding the restriction and prolongation operators are standard (Briggs et al. 2000 and, without loss of generality, assume that 1. A simple way to construct the prolongation and restriction operators, which satisfies Assumption 1 arises from the naive Nyström method (see Drineas & Mahoney (2005) for an introduction). In particular, we construct R and P as described in the definition below.
Definition 1. Let S,, {1, :.....n]] n. and denote SN C Sa, with the property that the elements N are uniformly selected by the set Sn without replacement. Purthermore, assume that 8; is the ith element of Sx- Then the prolongation operator P is generated as follows: the jth column of P is the 8t column of Taxn and, furthermore, we set R as P
Beyond definition several ways to select the restriction operators are available in the literature, e.g., Gaussian, 1- or s-hashing, and sampling matrices; see Cartis et al. (2022). However, for our experiments, we prefer uniform sampling for constructing these operators because it yields efficient iterations by avoiding expensive matrix multiplications.
In addition to the first-order coherency condition, we will also assume that the coarse model also second- order coherent,
(4)
A simple and surprisingly effective method to construct the coarse model in equation 2 while satisfying both the first and second order coherency conditions in equation 3 and equation 4 is the so-called Galerkin model,
- F(y): = (RV - yo) + RV (V' -o)) -oo).
In the context of multilevel optimization methods, the Galerkin model was experimentally tested in Gratton et al. (2010) where the authors found that it compared favorably with other methods to construct coarse- grained models. Motivated by the excellent numerical results in Gratton et al. (2010), and its simplicity, we will adopt the model in equation 5 as the coarso-grained model in the proposed algorithm. The model in equation also has close links with the recently proposed randomized Newlon methods which we briefly discuss below (for more details see Tsipinakis & Parpas [0024]; Ho et al. (2019)).
We compute the coorse direction by p y= Yo. Furthermore, in the case of the Galerkin model defined in equation we have the closed form solution,
H = arg min
[VV (V) -
(6)
dRERN
Note that the coarse direction is a vector in R. To correct the incumbent solution X we must prolongate it to Rn using the prolongation operator,
d = PdHA = -\V)
(7)
where the h and H subscripts denote directions related to the fine and coarse model respectively. If we naively set P as the identity matrix, the search direction in equation 7 becomes the Newton direction. When P is selected as in Definition 0 then we obtain the randomized Newton method which is based on uniform sampling. In this case, the cost to compute the reduced Hessian is O(nN) which is much cheaper than Onn2) of the Newton method (for details on how to compute the reduced Hessian matrix see Remark 2.1 in Tsipinakis & Parpas (2024)). The general multilevel method of this section achieves local super- linear convergence for both self-coneordant and strongly convex functions "Tsipinakis & Parpas (2024) Ho et al. (2019) However, when assuming self-concordant function, one can additionally show that the general multilevel method enjoys a global and scale invariant convergence analysis (Tsipinakis & Parpas 2024). In the next section, we develop a variant of the general multilevel method which achieves a global analysis with local super-linear convergence rate, which is also applicable to non-convex functions.
4

---

Published in Transactions on Machine Learning Research (01/2026)
3 Low-rank Multilevel Newton Methods
We saw above that the general multilevel method can be viewed randomized Newton method. In this section we discuss connections of the multilevel scheme with the low-rank Newlon method through the naive Nyström method and propose constructing the coarse model using a Truncated Singular Value Decomposition (T-VVD). The version of our method which is suitable for non-convex optimization problems is given in Algorithm !
Let A R be positive definite matrix and Y E ROXN with rank N < n. The Nyström method builds a rank-N approximation of A, namely AN, as follows,
AN = AY(Y"AY) A
(8)
The Nyström method has been extensively studied and has been shown to be efficient low-rank approx- imation method for different sampling techniques which gives practitioners the freedom to choose from a wide range of random matrices (for more details see Drineas & Mahoney (2005); Gittens (2011); ; [moola & Schölkopf (0000); Williams & Seeger (2000)) ). Then, one may obtain the naive Nyström method when Y is selected as in Definition I Similarly, we may obtain a rank-N approximation A by T-SVD,
Ax = UNENUT
For A being a symmetric positive definite matrix, EN E NNN a diagonal matrix containing the N-largest eigenvalues of A, with σi 02 2 2 ON, and Un E ROxN contains the corresponding eigenvectors. Although the naive Nyström method is more efficient compared to the T-SVD method, for the algorithms propose below we employ the T-SVD as it gives direct access to the eigenvalues of the Hessian matrix.
3.1 Convergence Analysis of Low-rank Newton Method for Self-concordant functions
Based on the Nyström method in equation g substituting A with v²](x;), Y with P and multiplying right and left with [V2 [xx)(x we obtain,
V [(x)) P(RV)[(x)PP))
(9)
which implies that the coarse direction equation can be viewed as an approximation of the Newton direction that is based on a low-rank approximation of the Hessian matrix. The fact that equation g offfers a low rank approximation of the Hessian matrix gives us the freedom to design search directions using different low-rank approximation approaches. Here, similar to the approach proposed in Erdogdu & Montanari (015
we perform a T-SVD approximation method (note that the work of Erdogdu & Montanari (2015) is based on sub-sampling and thus it is different to ours and not directly comparable). In particular, we construct an approximation of the inverse Hessian, namely Qi,kk by computing the (N + 1)sh T-SVD f2f(x; and
then replace all the cigenvalues after the (N + 1)th cigenvalue with ON+1- Formally, Q is constructed follows:
(10) where UR and Σa,N are the matrices of eigenvectors and eigenvalues obtained by the T- SVD process at Xs, respectively. Thus, by definition, the eigenvalues of Q have the form diag "i.o,,,,..,,,/,,.,.
Q = OF + U (E - INxN)
Then, du, = - V((x( is a descent direction. As result, for convex (stronlgy self-concordant) problems, this approach is efficient when the valuable second-order information is concentrated on the first, say N, eigenvalues. In particular, eigenvalues of large magnitude correspond to more informative directions since in these directions the objective function has a large curvaluur,, whereas for eigenvalues that are close to zero the curvature becomes almost flat. Therefore, by employing equation [0 in the computation of the search direction equation 7 we aim to determine the subspace spanned by the dominant cigenvectors associated with the largest eigenvalues, which are those that yield the fast convergence rate of the Newton method. Such problem structures are typical in machine learn- ing problem (Berahas et al. 2017); e.g., the Hessian matrix is (nearly) low-rank and/or there big gap
5

---

Published in Transactions on Machine Learning Research (01/2026)
between the nN and NN+1 eigenvalues. In the convergence analysis, we make use of the Newton decrement, which is given by,
((x+) = (V 1/2
We also define the following quantities
k,n
E/ :
6: lim inf ik R=+00
Smin := inf & EENN.
0,N411
We are now position to present the convergence analysis of the method. The convergence split inlo two phases according to the magnitude of {(x+). If X(xx) is sufficiently small, then the method enters the fast region of convergence. The fast region is governed by 17 = (3 v - 42min//. As the theorem shows, the rate of convergence depends on the magnitude of E.
Theorem 3.1. Let be strictly conver self-concortont function and suppose that the sequence (xx)NEN is generated by X = X3= 1. Q V fxxx), where Q in 10). Suppose also that e & 0. Then, there exist Y 0 and 17 (0, 3-vB) such that
(1) if {(x+) 7,, then f(xx+1) - f fxx)) --7
(ti) if {(x+)= S n, then the line search selects the unit step and
(a) if (0, 1), then the method achieves a composite convergence rate, i.e.,
A(x++1) S 1 1 + X(X+) M(xs) < A(xs)
Emiu
(1 (1x)(xx))²
(b) of 1, then the method converges with at least super-linens rate, i.e.,
A(xx+1) 1 1- 1 + A(x+)
X(xx) (1 (1- x(xx))
The proof of the theorem appears in Appendix The convergence rate depends on the value of E. If (0, 1); then the method attains a composite rate. If = 1, then the sequence (Ex)KeN converges to one since lim sumsup--oe & 1. In this case the method converges at least super-lincarly. Quadratic convergence rate will be achieved æ and there exists some Ko E N such that F for all k: ko- The theorem does not discuss explicitly the case
0. However, it is readily seen that if this scenario occurs, then 7 0 and thus the method converges in phase [i) Furthermore, unlike the classical theory for strongly convex functions, the advantage of using colf-conccrd as our main assumption results in obtaining a global analysis that has an intuitive local fast convergence rate, which only depends on the ratio between on and dN+1- The cost per itcration of the proposed method O(Nn?). This is much better than ( of the full Newton method.
3.1.1 Extension Non-convex Problems
For self-concordaant and other convex functions finding x* can in many instances be considered as an easy task (c.g., when descent methods are applicable) since all local minima are global. In this case, the unique global minimum can be attained at an x for which [|V(())] = 0. On the other hand, for non-convex problems, finding x" in general an NP hard problem. For the latter class of problems, a point x for which [[V = 0 can be one of the many local minima or a critical point such matrix is indefinite. Here, non-convex function f. To achieve this, we propose a variant of the low-rank Newton method which we conjecture will have a better escape rate from saddles compared to (stochastic) first-order methods (our conjecture is validated by numerical experiments).
saddle, for which the Hessian concerned with the task of finding local minimum of a general possible
Constructing the coarse dircction using the definition in equation 10 is particularly suitable for convex problems since all eigenvalues are positive. Then, for any i E f1. 1) we have o, > which ensures the descent property of dz However, when dealing with non-convex functions such guarantee may not

---

Published in Transactions on Machine Learning Research (01/2026)
Algorithm SigmaSVD
Input: p < N, E (0, E (0.0.)8). e (0.0.5), E (0, Ps E RaXN X4.0 R"
Compute [Q by equation 16 using the randomized T-SVD (Halko al. 2011
Form | by equation 17
Quit if (V (Xn,))())))) d) S
5: Armijo search: while fi (x + tu x > (xn,)) + atn, (nnk tn + Sta,
Update: Xh,4+1 == Xb, + th tn.dddnsl go
Return X4,k
be true, e.g., at a neighborhood of saddle point. Since for the Newton method the negative gradient is pre-multiplied by the inverse Hessian matrix, negative eigenvalue will result in changing the sign of the corresponding gradient entry which may yield the Newton method converging to saddle point or a local maximum. On the other hand, the Newton method breaks down when there are zero eigenvalues. To address these shortcomings, we compute the (N + 1)th T-SVD of the Hessian as above, but here we replace the negative eigenvalues by their absolute value. Further, sufficiently small eigenvalues are replaced by positive scalar to ensure the non-singularity of the approximated Hessian matrix. Formally, we define 9k.i N + RR, 1,2,..... N, such that,
Ek, 4> 0,
(11)
and [R == diag(ge, (( ..... NENN [k,,))). Then we obtain the truncated low-rank approxima- tion of the Hessian matrix as follows,
Q k == _o(ok.N+1) "Ixn + UR.N (gh -
(12)
Defining the approximated inverse Hessian matrix as above necessary order to obtain a descent direction, which together with an appropriately selected step-size parameter we can guarantee the descent nature of the Low-Rank Newton method for non-convex problems. Despite the lower iteration cost compared to the Newton method, there still may be cases that forming the Hessian matrix and computing the T-SVD is too expensive for some applications (when is extremely large). In the next section, we address this issue.
3.2
Sharse-grained Low-rank Newton Method with Analysis for Seff-concordant Functions
The computational bottleneck of the procedure described in the previous section arises from the fact that the computations are performed over the full Hessian. To address this, we perform the T-SVD on the Hessian matrix of the reduced order model. We begin by describing a low-rank multilevel method which is suitable for couvex optimization. Later, we will present the extension for non-convex functions.
Let QH, be T-SVD-based rank-p approximation of the reduced Hessian matrix in equation where p < N < n. Then, as in equation 9 we take
1 PQ H. PIRV P Qi ( f(x))
which in fact is a rank-p approximation of the inverse Hessian matrix with the difference that the computa- tional cost of forming Q is significantly reduced since we compute [U,,p41, Ex,ttll (by T-SVD) over the reduced Hessian matrix equation equation 10 thus, we define
We wish our Hessian approximation to have a similar structure as in
:=::= o kkpti INNN + UR CEE -
kkpt I
(13)
Therefore, dr, - V is an approximation of the Newton direction. We call this method SigmaSVD as it is approximation of Sigma in Tsipinakis & Parpas (2024). Let us now define the approximate decrements for both SigmaSVD and the general multilevel method of scction 2
where URP Nxp.
A(xx) = ( V 1/2
(14)
Á(x+) = (V

---

Published in Transactions on Machine Learning Research (01/2026)
respectively. Both quantities analogous to the Newton decrement ((x+) and serve the same purpose. In this section, we will construct the operators randomly at each iterations according to Definition 1 thus we slightly change the notation of random operators to Pa and RR. This way we enhance the applicability of the method. The following assumption emerges naturally.
Assumption 2. For all EEN N il holds {(xx) > O with some probability 6 > 0.
The above assumplion is Lypical when analyzing probabilistic subspace methods (Cartis et al. 2022 Tsip-
inakis & Parpas, 2024 Ho et al., (019)). Since {(x;) is norm of R2 (xxe), assumption 2 effectively prevents the latter from becoming zero with some probability. This assumptions is needed when analyzing these methods because (x =
can occur when R E mxll(Vf(x))).
When assuming self-concordant functions, the convergence analysis of the proposed method is similar that in Tsipinakis & Parpas (2024), with the difference that the lerm that controls its convergence rate is given by é == (x²(x)) (xx)) instead of ég == (X²(x;) (xxx)) It has been shown that EN (1 x2)2/ 0 ph S Tsipinakis & Parpas 2024). Intuitively, it should be expected éN to further incorporate the information carried out by the T-SVD. Self-concordant functions allow to prove such informative upper bound.
Lemma 3.1. For all & N we have that
1. d. k (2f((x) dn. {x(xx)
a fxx )dn, S (2(x
f[ (dm. - an,r)l S êk
4. 0epp+1
V
Oc. {(xx) s ((xx) s }(xx) S A(xx)
Further, lel Assumplion hold. Then, there erists pa E (0 such that éx < (1 probability 6.
2/22\()))))
with
The proof of the theorem below is based on the proof technique of Theorem 3.1 quantities
We define the following
Eg Ok.N Emin inffEx le E NN, =:= liminffe Ek.
Tkpt+1
k-+00
(15)
1 == lim inf pk: Pmin == inf(ue k E
k
We will show that the region the fast convergence is governed by f == 3- 3-v5+4 where := /I - Emiin
Theorem 3.2. Suppose that f is a strictly self- coff-concordant function, Px selected as in Definition ! and
Emin Mmin * 0. Suppose also that Assumption 2 holds and the sequence (λą)жeN is generated by X
ls Q V (xxx). Then, there exist constants > O and i E (0, 3-155 such that
(i) if A(xn,)) fi then
fn(x),4+1) - fn(xn,,) S --9,
with probability at each E N,
(i) if ((xx,) I f, and
(a) if either E e 10,1) or μ E (0, 1), then
-I- Emiu + A(X*) ((x+) ((xx),
X(xe+1)
(1 xxxx))
and thus the sequence (A(x)))kEN achieves composite rate with probability 8 at each k N,

---

Published in Transactions on Machine Learning Research (01/2026)
(b) of 1, then
1(x++1) X(xx)
ENPL + A(X1)
(1 - )(xx))
and thus the sequence (A(x)))KEN converyes super-linearly with probability at each k E N.
The proof of the theorem appears in Appendix C The difference between theorem 3.1 and 3.2 is that the latter involves Mk; which effectively is the quantity that controls the level of approximation bet ween the coarse and fine directions. Therefore, a slower convergence rate is expected compared to the low-rank Newton method section 3.1 μa approaches one, then SigmaSVD will achieve the fast rate of the low-rank Newton method. The method is also an approximation the multilevel method in [Tsiinakis & Parpas (2024). Theorem 3.2 shows that SigmaSVD will converge similar to the standard multilevel method in Tsipinakis & Parpas [0024 if ε is large enough. As in theorem 3.1 if E 0 μ 0, then and the
method will converge in phase 11)
3.3 Analysis for Non-convex Functions and Polyak-Lojasiewicc Inequality
In this section we analyze SigmaSVD) abandoning self-concordance and convexity from our assumptions. Specifically, assuming Lipschitz continuous gradients, we show reduction in the value of the objective function for sequences (xx)keN generated using the truncated corse direction equation I7 below. If we further assume that our objective function satisfies the Dolyak-.ojasiewice (PI.) inequality, then the proposed method will converge with a lincar rate.
Based on the discussion in section 3.1.1 and given the eigenvalue decomposition which, as in 3.2 E and Ukp the truncated reduced Hessian matrix is defined follows
returns
1 == 0k,++ (Ekg+1) INNN + (9np()((Ep)) - -h.p -E
(16)
then the truncated coarse direction given by
d == - - =
- -P/
(17)
The full algorithm together with step-size strategy is specified in Algorithm When the objective function is convex or sell-concordant Algorithm | coincides with SigmaSVD in section 3.2
Assumption 3. The gradient of is Lipschitz continuous on R", that is, there ewists M such that for all x,y E T
M 2 lly x|
f(y) < f(x) Vf(x)) x)
A diroct consequence the above assumption is the boundness of the Hcssian matrix from above: v2f(x) = MI,, for all R
Assumption 4. There exists ğ E (0, such that R4 V fx(x)) 2 Vf with probability 0.
Assumption 4 is typical when analyzing subspacc methods, see Cartis ct al. 2022 Cartis & Roberts (2024)
Cartis et al. 2025). In this section we use it in place of assumption to prevent inellective coarse steps with
probability 8.
Theorem 3.3. Let assumptions
3 and 4 hold. Suppose also that the sequence (xi)eeN is generated by l 1 Mar
X4+1 = + 4/ dd. and 1o <
Then
A2v
f(xx+1) - f(xx) S 2w02 M
[Vf(xx))
with probability 8.

---

Published in Transactions on Machine Learning Research (01/2026)
The above result shows that can obtain descent method when the function has Lipschitz gradients. The result is global and remains true with probability 8. The method will not diverge since, with probability 1 ò- f may be zero, and thus f(x++1) = f(x)). However, due to the randomness of the method, at future iteration, it is expected that the coarse model will not be ineffective (and thus j * 0), which effectively yields a convergent algorithm. In order derive a rale of convergence, we make use of the following assumption.
Assumption 5. The PL inequality holds; there ezists 5 > 0 such that
[ 2
((((=)) f(x")
(18)
The PL condition is weaker than convexity (see Karimi et al. 2016) for more details) and has been (empiri- cally) shown that it is often satisfied for over- wee--arameterized neural networks (see Belkin (2021) and references thercin). We think that this result goes some way to explaining the behavior of the algorithm (at least close enough to a minimum).
Theorem 3.4. Let assumptions and |oupe [dn,el- Then
hold.
Assume that (xx)REN is generated by Xk-1 = Xx +
f(xx+1)- - f(x))= s (1 (1- 2w² (f(x=) (oxs))-((x')) f(xx))
VEA²
with probability 8.
Refer to Appendix D for the proof. The theorem presented is global, demonstrating that convergence to f(x*) can be achieved with at least a linear rate. Theorems 3.3 and 3.4 become particularly elfective when both f and are sufficiently large. Objective functions with high f and often include those with low-rank Hessians or concentrated second-order information in the dominant eigenvalues. Functions satisfying the PL inequality typically exhibit low-rank Hessians. Given the method's low per-iteration computational cost (O(N + pN²) in parallel), algorithm is expected to escape saddle points efficiently and converge rapidly to local minimum. This will be empirically validated in the following section.
Remark 3.1. Since the Euclidean norm is used to oblain the converyence results in this section, we are in a position to quantify the parameters P and å depending on the choice of the operator R. 1f R is Goussian matric whose entries are drawn from N(, N-)), then for any P E (0, 1) we have é e (2- (Cartis et al., (022). Another possible choice that gields sporse matrices, and is therefore less computationally expensive than using dense Gaussian matrices, is given by s-hashing matrices (Kane & Nelson) 2014)_ When using s-hashing malrices to construct the coarse models, then for any fa E (0, we have 8 E N(1 NNL-=2) /C, problem-dependent
where = C2(1 - N2)N,, and C1, C2,C2 are problem- dependent constants. To the best of our knowledge, such
explicit quantitative values for f and under defivition 1 are not available in the literature.
Remark 3.2. In praclice, we select N and as follows: When R is constructed as in definition 1 these parameters should be selected according to the effective rank of the Hessian at the optimum, R= tr(V= f(x")) The quantity R indicates how many directions really matter near the solution, where most methods slow down significantly. If R is known, in proctice N and p should be selected larger and slightly smaller than R, respectively. Alternatively, they can be selected in similar fashion according to the quantity maxfnn Rlun), suggested in Erdogdu & Montanari] 2015). On the other hand, if R. is Gaussian s-hashing N and p should be selected according to Remark 3.1 i.e., by specifying ja and 8 one obtains N directly-
Amar f(x)))
4 Numerical results
We are now ready to validate the efficiency of the proposed method on different machine learning models. We illusratte that the method is efficient and compares favorably with other stato-of-the-art algorithms broad class of problems. In particular, apply algorithm problems minimizing self-concordant.
10

---

Published in Transactions on Machine Learning Research (01/2026)
EigaatyO NERSE
SigtadeD
10-4
CLBl REWion
Norton
1N0)
Humaeii p015
11
GPU time (ooe)
CPU Sizzo (Dec)
GPU time (sce)
(a) Gisette
(b) Gisette
(c) Gisette with Xo
Figure 1: Non-convex minimization. All the methods in plot (a) are initialized at the origin, while in plot (b) the initializer selected randomly by (). Plot (c) shows the convergence behavior of SigmaSVD for different values of p.
fanctions, strongly conver functions with or without Lipschitz continuity and a non-conver functions. this section we only report the results for the non-convex cases. The additional numerical results together with detailed description of the algorithms, training datasets, and the setup used to obtain the results appear in the Appendix
4.1 Non-linear least-squares
Suppose that are given training set (a,,b, with a, Rh and be [0, and consider solving the non-linear least-squares problem,
(b: o(a{x))²
1 exp(w)--I
where is a real function, known as sigmoid. Since we apply the mean squared error over the sigmoid function, the non-linear least-squares is non-convex problem. We compare the performance of algorithm [ against gradient descent (GD), accelerated gradient descent (AGD), cubic Newton and Adam. Although Adam, as batch algorithm, is not directly comparable to SigmaSVD, we include it in our comparisons to demonstrate the advantages of second-order methods compared to the stochastic first-order method on problems with several saddle points and flat areas. For algorithm [ GD and AGD the Armijo rule is applied to determinc the step size parameter with constants a = 0.001, 3 = 0.7. A line search strategy applicd select the regularization parameter the Cubic Newton method (for details see Nesterov et al. (2018)). In all experiments we select the fixed eigenvalue threshold, ie., V 10-10 10 in equation 11 The momentum
parameters for Adam are selected suggested in Kingma & Ba 2014 while for AGD the momentum is selected as 0.5.
In Figure ! we demonstrate the performance of the different optimization methods on the non-linear least- squares problem for the Gisette dataset (see also Figure 3 in Appendix E for simulations on different data regimes). In Figures Ia and Ic we illustrate the reduction of the value function over GPU time. Observe that SigmaSVD able to return a better solution compared to the first-order methods. Clearly, first-order methods are trapped in a flat area (Hessian and gradient are zero) and thus they are stuck far from the local or global minima. This is not surprise since in such areas the gradient is almost zero and thus first- order methods unable to progress and although in theory they have been shown to always escape even illconditioned saddles, here we see that they require an extremoly large amount of iterations which makes them inefficient for practical applications. On the other hand, Cubic Newton is able to escape saddle points in one iteration, however we observed that many experiments it converges rapidly to early local minima probably due to the lack of randomnoss. Morcover, as suggested Theorem 3.2 Figure I illustrates that SigmaSVD enjoys a fast local convergence rate and is faster than its counterparts. It is also faster and has the ability to achieve lower training errors for different initialization points.
11

---

Published in Transactions on Machine Learning Research (01/2026)
Table 1: Probability of successfully escaping from saddle points and convergence to the global minimum for various values of N and fixed p = 450, X = 0. Each row shows the probability over 50 trials.
Escape Rate Probability - Gisette
N
Probability
0.1n
18%
0.13n
46%
0.26n
52%
0.36n
66%
0,42n
80%
0.46n
92%
In Table aan Figure Tc we demonst how the choice of N and p affects the behavior of SigmaSVD around saddle or flat area. We revisit the experiment of Figure [b (i.e., Xo =0) 0) in which SigmaSVI escapes from such the flat area in one iteration (similar to Cubic Nowton). Figure shoosw the behavior of SigmaSVD for different values of p with fixed 0.5m. Observe that for very small values of p SigmaSVD is trapped the same saddle as the first-order, however for p 2 350 it converges to the global minimum. Next, Table we fix p 450 and vary N to show the escape rate probability from this region, otherwise the probability of convergence the global minimum over 50 trials. Similarly to p, see that the escape rate probability is proportional to N. Both experiments verify Theorem 3.2 which indicates that for very small values in p and N, one should expect SigmaSVD to behave as a first-order method. Furthermore, Table | and Figure show that SigmaSVD is able to reach the behavior of the Cubie Newlon method (that is, escape from the saddle point one iteration) using only half the dimensions of the problem and about 20% of the cigenvalues of the reduced Hessian, which constitutes a significant improvement in total complexity as also depicted in Figure
4.2 MNIST deep autoencoder
In this section we consider training the MNIST deep autoencoder which is based on the work presented in Hinton Salakhutdinoy (2006). The MNIST auto-encoder is considered a benchmark optimization problem because the presence of saddle points along its optimization trajectory poses significant challenges to optimization algorithms, see Reddi et al. (2018) and reference therein. The MNIST auto-encoder consists of an encoder with layers-size 28 28, 1000, 500, 250, 30 that it is followed by a symmetric decoder and in total has 2.8M paramcters. Further, the sigmoid activation function is applied to all layers.
Here, we compare SigmaSVD against Adam. Since the MNIST autoencoder problem is large, to avoid memory issues, we adapt SigmaSVD for batch learning. In addition we use first-order momentum, namely 81 E (0, 1), similar to that of Adam. In particular, our new itcrations have the following form
x = X - ttP Rmx+1,
where
m++1 1 3
m
m = B1 m. + (1 31) V
As a result, SigmaSVD with momentum effectively applies Adam's updates, substituting the Hessian for the second moment of the gradients (which a diagonal matrix). Our goal demonstrate that the Hessian can be essential in deep architectures and that its diagonal approximations may severely slow down optimization algorithms or lead to suboptimal performance in the presence of saddles points and flat areas.
To obtain the results the batch size is set to 128. Furthermore, for Adam, we set the learning rate to 0.0025 (yields the best convergence rate) while keeping the momentum parameters at their default values (Kingma & Ba 2014). For SigmaSVD, we use two coarse models of different sizes: one with N = 2,800 and tx = 0.01, and another with 1,400 and & = 0.02, where for both we used 31 = 0.9. In both cases, we found that setting v= 10-8 yielded the best performance. Additionally, we use a full eigenvalue decomposition instead of truncated SVD. This is feasible because the coarse models are small, allowing for efficient computation.
12

---

Published in Transactions on Machine Learning Research (01/2026)
Adam
Adam
SigmaSVD, n=2800
SigmaSVD, n-2800
SigmaSVD,
SigmaSVD, n-1400
10t1
10
Test error
Thaining
0.0
2.5
5.0
10.0
12.5
15.0
17.5
20.0
0.0
2.5
5.0
7.5
10.0
12.5
15.0
17.5
20.0
Epochs
Epochs
MNIST
MNIST
Figure 2: Minimizing the MNIST deep autoencoder. Performance comparison between SigmaSVD and Adam. The left plot shows the reduction in the value of the objective function and the right plot the generalization performance versus epochs.
A full SVD is also advantagcous as it climinates the need to tune extra hyperparameter (p in this case). More importantly, selecting a larger generally leads to better escape rates near saddle points, as observed in Figure 1c. We ran the experiments on a single A100 GPU with 40Gb RAM.
The results of this experiment are shown in Figure 2 The left plot illustrates how quickly the algorithms reduce (x) over epochs. Clearly, both SigmaSVD variants outperform Adam in terms of convergence speed in the first 20 epochs. This significant difference in convergence rate stems from the presence of multiple saddle points and flat regions along the optimization trajectory. In particular, during the early stages of training, we observed instances where the gradient norm was nearly zero (indicating saddle points) and cases where all eigenvalues of the reduced Hessian were zero (suggesting a possible flat region). SigmaSVD handles these situations efficiently, escaping quickly by leveraging second-order information, whereas Adam, relying only on a diagonal approximation, struggles in such ill-conditioned settings. Additionally, the right plot compares the error, showing that the rapid convergence in training does not translate to overfitting, Additionally, note that the proposed methods achieved these results while updating only 1, 10400 and 2, 800 parameters per iteration, respectively, whereas Adam updates all 2.8M parameters. This highlights the significance of employing rich and informative preconditioners, such as the reduced Hessian, in optimizing ill-conditioned deep neural networks. Further, regarding the wall-clock time comparison, Adam completes an epoch significantly faster (18 seconds) compared to the SigmaSVD variants (600 and 1, 100 seconds, respectively). However, comparing convergence based on GPU time would be unfair this stage, as Adam is a highly optimized industrial algorithm, whereas SigmaSVD has been developed primarily for academic research. Nonetheless, this experiment suggests that, order to design efficient algorithms for minimizing large neural networks with flat regions and saddle points, one should consider hybrid approaches, combining computationally inexpensive method such as Adam when gradients are large, with more sophisticated method such as SigmaSVD near saddle points and flat regions.
5 Conclusion
We develop stochastic multilevel low-rank Newton lype method that enjoys a super-linear convergence rate for self-concordant functions with low per-iteration cost. We further propose a variant that is applicable to non-convex problems and establish its lincar when the PI. incquality is satisfied. Preliminary numerical experiments show that our method is efficient for large-scale optimization problems with millions of dimen- sions. The experiments also show that the method is well-suited to problems with low-rank Hessian matrices.
13

---

Published in Transactions on Machine Learning Research (01/2026)
It is also faster and has improved escape rate from saddles and flat-areas compared to first-order methods in non-convex cases. As future direction, we aim to analyze the batch variant of our method and develop a hybrid approach that is efficient when training deep neural networks. We also plan to provide convergence analysis for general non-convex functions.
References
Mikhail Belkin. Fit without fear: remarkable mathematical phenomena of deep learning through the prism of interpolation. Acta Numerica, 30:203-248 2021.
Albert S Berahas, Raghu Bollapragada, and Jorge Nocedal. An investigation of Newton-sketch and subsam- pled Newlon methods. Xiv preprinl urXiv: 1eee10050022111 2017.
Raghu Bollapragada, Richard H Byrd, and Jorge Nocedal. Exact and inexact subsampled Newton met for optimization. IMA Journal of Numerical Analysis, 39(2):545 578, 2019.
Stephen Boyd and Lieven Vandenberghe. Conver optimization. Cambridge University Press, Cam- bridge, 2004. ISBN 0-521-83378-7. doi: 10.101//B0078051180444 URL httpt:://oo..oooo org/10..017/ CB09780511804441
htttp:://bd..org/10.11177
William L. Briggs, Van Emden Henson, and Steve F. McCormick. multigrid tutorial. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, second edition, 2000. ISBN 0-89871-462-1. doi: 0111137/..9708987711
URL httpss///aaii.orr//1113//////// 97888997790099
Coralia Cartis and Lindon Roberts. Randomized subspace derivative-free optp:aiiieeeeeeeooiiiiiiiai optimization with quadratic models and sccond-order convergence. arXiv preprint arXiu:241.144991, 2024.
Coralia Cartis, Jaroslav Fowkes, and Zhen Shao. Randomised subspace methods for non-convex optimization, with applications to nonlinear least-squares. arXiv preprinl arXiv:221100987. 2022.
Coralia Cartis, Zhen Shao, and Edward Tansley. Random subspace cubic-regularization methods, with applications to low-rank functions. arXiv preprint arXiv:2501.099 .09734, 2025.
Yann N Dauphin, Razvan Pascanu, Caglar Gulcehre, Kyunghyun Cho, Surya Ganguli, and Yoshua Bengio. Identifying and attacking the saddle point problem in high- high-dimensional non-convex optimization. Advances in neural information processing systems, 27, 2014.
Petros Drineas and Michael W. Mahoney, On the Nyström method for approximating a Gram matrix for improved ternel-based learning. Mach. Learn. Res., 6:2153-2177 2005. ISSN 1532-4435.
Murat A Erdogdu and Andrea Montanari. Convergence rates of sub-sampled Newton methods. In Proceedings of the 28th International Conference Neural Information. Processing Systems- Volume 2, pp. 4052-3060. MIT Press, 2015.
Terunari Fuji, Pierre-Louis Poirion, and Akiko Takeda. Theoretical analysis of the randomized subspace regularized Newton method for non-convex optimization. Open Journal of Mathematical Oplimizattoo, 6: 1 1-35, 2025.
Meirav Galun, Ronen Basri, and Irad Yavneh. Review of methods inspired by image analysis applications. Numer. Math. Theory Methods Appl. 8(2):283- 312, 2015. ISSN 1004-8979. doi: 10.4428/nmtma... w1si.. URL stttps://ooi.orgee 42000/442200000 11200/nmtman nmtma. 2015. w14si
for data and
Alex Gittens. The spectral norm error of the naive Nystrom extension. arXiv preprint arXio: 110.5555, 2011.
Robert Gower, Dmitry Kovalev, Felix Lieder, and Peter Richtárík. RSN: Randomized Subspace Newton. Advances in Neural Information Processing Systems, 32, 2019.
14

---

Published in Transactions on Machine Learning Research (01/2026)
Serge Gratton, Mélodie Mouffe, Annick Sartenaer, Philippe L. Toint, and Dimitri Tomanos. Numerical experience with a recursive trust-region method for multilevel nonlinear bound-consttrained optimization. Optim. Methods Softw., 25(3):359 386, 2010. ISSN 1055-6788, doi: 10.1060/106566780009333339 URL
Nathan Halko, Per-Gunnat Martinsson, and Joel A Tropp. Finding structure with randomness: probabilistic
algorithms for constructing approximate matrix decompositions. SIAM Rev., 53(2):217-22 2011. ISSN
0036-1445. doi: 2011370090771 806. URL ttttps://doii.org/10... 11t..011117/9907788000
Filip Hanzely, Nikita Doikov, Yurii Nesterov, and Peter Richtarik. Stochastic subspace cubic Newton method. In International Conference on Machine Learning, pp. 4027-4038. PMLR, 2020.
Geoffrey E Hinton and Ruslan R Salakhutdinov. Reducing the dimensionality of data with neural networks. science, 313(5786):504-507, 2006.
Chin Pang Ho, Michal Kočvara, and Panos Parpas. Newton-type multilevel optimization method. Oplimiza- tion Methods and Software, pp. 1 1-34, 2019.
Roger A. Horn and Charles R. Johnson. Malrix analysis. Cambridge University Press, Cambridge, second edition, 2013. ISBN 978-0-521-54823-6.
Vahan Hovhannisyan, Panos Parpas, and Stefanos Zafeiriou. MAGMA: multilevel accelerated gradient mirror descent algorithm for large-scale convex composite minimization. SIAM J. Imaging Scí., (4):18829 1857, 2016. ISSN 1936-4954. doi: 10.1137/15M100
URL http://....rr/11111111////10..00
Daniel M Kane and Jelani Nelson. Sparser
transforms. Journal of the ACM (JACM),
(1(1)11-23, 23, 2014.
Hamed Karimi, Julie Nutini, and Mark Schmidt, Linear convergence of gradient and proximal-gradien meth- ods under the Polyak-Loojssiewiz condition. In Machine Learning and Knowledge Discovery in Databases: European Conference, ECML PKDD 2016, Riva del Garda, Italy, Seplember 19-23, 2016, Proceedings, Part I 16, pp. 795-811. Springer, 2016.
Diederik Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXio preprinl arXiv:1412. 6980, 2014.
Ulysse Marteau-Ferey, Francis Bach, and Alessandro Rudi. Globally convergent Newton methods for ill- conditioned generalized self-concordant losses. Advances in Neural Information Processing Systems, 32, 2019.
Sen Na, Michal Derezińikki, and Michael W Mahoney. Hessian averaging in stochastic Newton methods achieves superlinear convergence. Mathematical Programming, 201(1):473-55) 520, 2023.
Yurii Nesterov. Introductory lectures on conver oplimization, volume 87 of Applied Oplimization. Kluwer Academic Publishers, Boston, MA, 2004. ISBN 1-4020-7553-7. doi: 10.1007/978-1-40 8853-9. URL basic course.
kttpps///oi: 1ttp://0000000//1111111-----9-00
Yurii Nesterov and Arkadii Nemirovskii. Interior-point polynomial algorithmss in conver programming, volume 13 of STAM Studies in Applied Mathemalics. Society for Industrial and Applied Mathemai- ics (SIAM), Philadelphia, PA, 1994. ISBN 0-89871-319-6. doi: 10.118771979111117777111 URL
https://ádi... org/ 10. 1117/11 9781611970791
Yurii Nesterov al. Lectures on conver optimization, volume 137. Springer, 2018.
Thomas O'Leary-Roseberry Nick Alger, and Omar Ghattas. Incxact Newton methods for stochastic non- convex optimization with applications to neural network training, arXiv preprinl arXiv:1905.00 067388, 2019.
Razvan Pascanu and Yoshua Bengio. Revisiting natural gradient for deep networks. 13X 2013.
arviv preprint
15

---

Published in Transactions on Machine Learning Research (01/2026)
Santiago Paternain, Aryan Mokhtari, and Alejandro Ribeiro. A Newton- based method for nonconvex opti-
mization with fast evasion of saddle points. STAM Journal on Oplimizalãon, 29(1):343-35 -388, 2019.
Mert Pilanci and Martin Wainwright. Newton sketch: A near linear-time optimization algorithm with lincar-quadratic convergence. STAM Journal on Optimization, 27(1):205-24) 245, 2017.
Sashank Reddi, Manzil Zaheer, Suvrit Sra, Barnabas Poczos, Francis Bach, Ruslan Salakhutdinov, and Alex Smola. A generic approach for escaping saddle points. In International conference on artificial intelligence and statistics, pp. 2233-1242. PMLR, 2018.
Jeflrey Regier, Michael Jordan, and Jon McAuliffe. Fast black-box variational inference through stochastic trust-region optimization. Advonces in Neural Information Processing Systems, 30, 2017.
Farbod Poosta-Khorasan and Michael W Mahoney. Sub-sampled Newton methods. Mathematical Program- ming, 174(1):293 326, 2019.
Jun Shang, Haishan Ye, and Xiangyu Chang. Accelerated doublc-skotching subspace Newton. European
Journal of Operational Research, 319(2):484-498, 2024.
Sidak Pal Singh and Dan Alistarh. Woodfisher: Efficient second-order approximation for neural network compression. Advances in Neural Information Processing Systems, 11118098-11000 2020.
AJ Smola. and B Schölkopf. Sparse greedy matrix approximation for machine learning. Seventeenth International Conference on Machine Learning (ICML 2000), pp. 911 918. Morgan Kaufmann, 2000.
Nick Tsipinakis and Panos Parpas. A multilevel method for self-concordant minimization. Journal of Opti- mization Theory and Applications, pp. 1-51, 2024.
Zaiwen Wen and Donald Goldfarb. A line search multigrid method for large-scale nonlinear optimization.
SIAM J. Optim., 20(3):1478 -508-10,3, 2009. ISSN 1052-6234. doi: /0.1137/08071524X. URL hutpps://oii, org/ 10. 11371/0717/8874
Christopher Williams and Matthias Seeger. Using the Nyström method to speed up kernel machines. Ad- vances in neural information processing systems, 13, 2000.
Yuhuai Wu, Elman Mansimov, Roger B Grosso, Shun Liao, and Jimmy Ba. Sccond-order optimization for deep reinforcement learning using kronecker-factored approximation. In NIPS, 2017.
A Background
Throughout this paper, all vectors are denoted with bold lowercase letters, i.e., x E IR. and all matrices with bold uppercase letters, i.e., A e _m=x The function ||x||2 (x", x22 is the fa- or Euclidean norm of x. The spectral norm of A the norm induced by the f--norm on R" and it is defined as ||A||2 ||x 12
max
|| 2|
It can be shown that ||A||2 = 01 (A), where 0ı (A) (or simply 01) is the largest singular value of A, see Horn & Johnson 2013 Section 5.6). For two symmetric matrices A and B we write A > B when xT(AA B 2 0 for all x E R", or otherwise when the matrix A- - B is positive semi-definite. Below we present the main properties and inequalities for self-concordant. functions. For a complete analysis see Nesterov (2004); Boyd
& Vandenberghe (2004). A univariate convex function : R + R called self-concordant with constant
M, 2 0 if,
(0"( Mz0"(x)3/2,
(19)
Examples of self-concordant functions include but not. limited to linear, convex quadratic, negative loga- rithmic and negative log-determinant function. Moreover, the class of self-concordant functions includes p-strongly convex functions with L-Lipschitz continuous Hessians, as the latter are sclf-concordant with constant Ma L/22)3/2). Therefore, if f is throc-times differentiable, the self-coneordance assumption is less restrictive than the classical assumptions, see Nesterov et al. (0018.
16

---

Published in Transactions on Machine Learning Research (01/2026)
Based on eq. 19 multivariate convex function f R" - R is called self-concordant if ((t) := f(x tu) satisfies equation is preserved under composition with any affine function. In addition, for any convex self-concordant function M²
19 for all x dom f; E IR" and / E IR such that x x E dom f. Further, self-concordance
è with constant Ma 2 0 it can be shown that ((z) := 7 (x)) is self-concordant with constant Ma = 2. Next, given x E dom f and assuming that f(x) positive-defintt we can define the following norms
||u||x := (xx)( w///
1/2 and l|vll= (v.)
(20)
for which it holds that (u, v) S | Therefore, the Newton decrement can be written as
Xy(x) = 1 [V
(21)
In addition, take into consideration two auxiliary functions, both introduced in Nesterov (2004) Define the functions w and wa such that
w(x) = x + and W- (x) = - 1 -oogoo xx) x),
(22)
where dom + far R: 2 0} and dom w. fa E 0 < 1}, respectively, and log(æ) denotes the natural logarithm of æ. Moreover, note that both functions are convex and their range is the set of positive real numbers. Further, from the definition equation T9 for Ma 2, have that
d dt
(1))
s 1,
from which, after integration, we obtain the following bounds
)0 + (0)2222 S d((1) (1) S (1 (0)
(23)
where the lower and the uppers bounds hold for t 2 0 and t E [0, (00) with t E dom ø, respectively (see also Boyd & Vandenberghe (2004)). Consider now self-coneordant functions on IR" and let S(x) =
{y E RR"
ly x < 1J. For any x dom f and y E S(x), we have that Nesterov (2004))
(1 - lly x x )xx
V [1 ly x f(x).
(24)
The following results will be used later for the convergence analysis of the algorithms. For their proofs see Nesterov (2004).
Lemma A.1 ((Nesterov 2004)). Let x, E dom f. If lly - x|1= < 1, then,
f(y) S fx) (V/(x),) - -x) + w. (lyy - x1|x)
Lemma A.2 Nesterov (2004)). Let x.y E dom f. Then,
f(y) 2 f(x) (Vf(x),y x) +wIVf((y) - v((x))))).
If in addition [[V [x())] < 1, then,
f(y) S f(x) + (V/(x), - x) + +w. (V/)
Throughout this paper, refer to notions such linear and quadratic convergence rates. Denote X the iterate generated by an iterative process at the tt iteration. The sub-optimality gap of the Newton method for self-concordant function satisfies the bound f(xx) - f(x') s (x which holds for (y(x)) S 0.68 [Boyd & Vandenberghe (2004h) and thus one can estimate the convergence rate in terms of the local norm of the gradient. It is known that the Newton method achieves a local quadratic convergence rate. In this setting, we say that a process converges quadratically if X (xx++)) I Rxy(xx)?, for R 0. In addition to quadratic converge rate, we say that a process converges with super- -ineer rate if \y(xx+1)\\((xx) s R(k), where R(k) O.
17

---

Published in Transactions on Machine Learning Research (01/2026)
B Proof of Theorem 31
The update rule of the Low-Rank Newton method at iteration & is
Xk+1 = Xx: + tx tg
(25)
where for the purposes of this section we use a == - Qã,kk Vi.kvv(se (x+e), where Q is defined in equation 10
Define also the corresponding approximate Newton decrement X(x) = (V)(x)) ) (x)))
Lemma B.1. For any k E N it holds that
= d
(26)
okin A(xx) S )(xx) S A(x+)
(27)
0N,N+1
Id S {(x+)
(28)
[(>x(x)]] (an. - dr.it) dnix)1 s 5(1-
Tknn
A(XA).
(29)
JR,N+1
where || llx defined in equation 20
Proof. Equality equation follows directly by the definition of dr.x. For the right-hand side equation we get
x(xx)2 )(xx)² vf(x)) 7
-f((x) // -|
By construction of Q holds that
1, In 3 maax ( (xx)] ( I,,
o (VF
1 = 0, (2(
max I
and thus İ(x+)= - A(xx)= 1 0. Using similar arguments we take
ov(V ([V](((x))) X(x)2
)(xx)² )(xx)² < max
i= l=.....m
(Q(Qnn))
Th,n
X(x)²
JK.N+1
and thus the left-hand side inequality of equation 27 follows directly, As for the result in equation 28 using the fact that Qi, is symmetric positive-definite we take
id = a V )
T((x))](QnnN]
[ (
I Fmax (IQu.el (xx)))) )(xx)
18

---

Published in Transactions on Machine Learning Research (01/2026)
and since max definitions of dn,& and d have
the result follows. Last, we prove equation 29 Using the
v [(xx)) 1/2 (dn, - dn
= 22 (xxx (1V [[xx](xx)) - Quai ) /
(xx)
I
[ fxxx)
QQ 1 [xxx) 1(x)
max
1 -
x++1
X(x)
01
NN+1
N
1
On
X(x)
ON+1
which concludes the proof of the lemma.
Lemma B.2. Let {(xx) > for some yj E (0, 3-v5). Then, there erists >> 0 such that the coarse direction
an will gield reduction value function
f(x: + (gdn,)) - f(xx) S -%,
for any & E N.
Proof. By Lemma A.1 logether with equation 25 and equation 26 we have that
f(x* + txdn,x) & f(xx)= + (V fxxe) ) xh+1 - xx) + we (1x0+1 - x xxllxx)
= fxxx) + Vf (x+)) ) a + wa (dn
S fxx)) (xx + w. (tx
where the last inequality follows from equation 28 and the fact that w. (x) is monotone increasing function. Note that the above inequality is valid for ln ((x,,)
Next, the right-hand side is minimized at t
and thus x+
X(x)2
A(x)
f(x* + + d), S f(xx)
+ wk I+ A(x(.k)
1 {(x1,,))
= f(xx - A(x+) + log (1 + {(x+))
Using the incquality =x + log(1 + x) S 2(1+x)
x² for any x > O, we obtain
{(xx)²
f(xx + t' d'd,))) S f((xx) 2(1 + ((xx))
$ f(xx) - at" {(xx) )2
=xxx) at V fxxx) Tank.
thus t* satisfics the back-tracking line search exit condition and that it will always return stop size t > 3/(1 + {(x))). Therefore,
fn(x, + txd),.) f(xx) S -a, )(xx)²
1 {(xx)
Additionally, combining the lent-hand side in equation
27 the fact that the function a It
is monotone
increasing for any 0 and since by assumption A(x+) 17 obtain
ox 1/
f(x: + (,.n.k)) - [(xx)
1
which concludes the proof by setting = a:31/2/1 (11 + 7))
19

---

Published in Transactions on Machine Learning Research (01/2026)
Next, estimate the sub-optimality gap of the process. Lemma B.3. Let X(xx) 1. Then,
w(A(xa) S [(xx) m(x') S as (M(x)).
If in addition (xxx) < 0.68 then,
f(x() - f(x²)((==)))))))>))) - A(xx)²
Proof. The first result follows directly from Lemma A.2 Further, it holds that w. S 2²,0 2 x < 0.68 and thus we take the second result.
Therefore, )(xx)² & e,e E (0, 0.68"), can be used as an exit condition for the Low-Rank Newton method. Next, we prove that the line search selects the unit step.
Lemma B.4. If ((xx) s (1 22 then the Low-Rank Newton method accepts the unit step, ta 1.
Proof. By Lemma A.1 we have
f(x: + dn,x) S f(xx) X(x)² + w. (A(x))
f((x)) - X(x)2 1(x) - log(1 - X(x))
which holds since by assumption X(x) <. Further, - log(1 -x) for all E (0, 0.81) and hence
I(xx dm,)) f(xx) T(x)² _ [X(x)² X(x)*
f((x)) 2
(1 - (12xxx))((xx²²
Therefore, if ((xx) S (1 122)/2 /2 we take
f(x: + an,s) S f(xx) ax(x)²
which satisfies the line search condition for 1.
Lemma B.5. Let f _ be a strictly conver self-concordant function. If ((xx) 2, we have that
v2f(xx+1) 1 (1xx))(x)))= 2 f(xx),
(ti) [V2 (1=+
V (2xxx)]
Proof. The proof is analogue to Lemma
C.3 below but with X(xx) instead of )(xx).
The next lemma shows the local super-linear rate of the Low-Rank Newton method.
Lemma B.6. Suppose that the sequence is generated by equation 25 with to 1. Then,
0ha A(x++1) & o
+ A(xx)
(1-X(xx))² A(X6).
Proof. By Lemma B.5 the upper bound in equation 27 and to = 1 we have
v2f(xx+1) = {(xx))²
22 v²f(xx).
(30)
Now, by the definition of the Newton decrement we have that
X(x++1) /
20

---

Published in Transactions on Machine Learning Research (01/2026)
and combining this with inequality equation 30 take
|| [xx)
((x++1)
(31)
A(x1.)
Denote Z = [xx((x) (x) Using the fact that
[(((++1)= S/¹ v2I(xx + (dnx)dnn dy + (x(x))
we that
(2'v2f(xx 1/2 (L/ fix + ydnx)) dy + f
Z = [V* 1-13/
[xx( 1/2 f(x: + VV 1/2 dy [xx)))) /.,
- dn.k.
where da,k is the Newton direction. Next adding and subtracting the quantity TIV we have that
z =
- da.x) + (T [(xx) dd,
and thus
|z|| <
T fxx)) -nax) +
T INXN) V2 f(xx) d2
(32)
21
2,
Using again equation 30 and since, by assumption, yd(xa) < 1 we take
V ((x) + xl f(xx)
(1 y)(x+)) INNN.
We now in position to estimate both norms in equation 32 For the first one have that
2 fxx)) 1/2 (dn. - dnx))l
7 <
dy
(1 - y)(xx))
1
[xx((
(an.k -dxx)) 11
1 {(xx)
Tkn
A(x+)
X(Xn,))
where the last inequality follows from equation 29 Next, the second norm implies
Zz |1 <
dy
[xx))(x
(xxx)
A(xx)
)(xx).
- A(xx)
Putting this all logether, inequality equation becomes
Jh.n
+ *(xx))
X(x++1) 5
(1 - )(xx))² 2(xx).
av.N+1
required.
21

---

Published in Transactions on Machine Learning Research (01/2026)
Next prove theorem 3.1 Recall the
Proof of theorem 3.1 By lemmma B.2 see that one step of the first phase of the Low-Rank Newton method
decreases the objective by := aBn (1 + 77) > O. By Lemma B.2 we also note that the line-search can be applied to determine the step size parameter. This proves the result of phase [i) of the Low-Rank Newton
method.
[ii we have {(x+) S 7. Further, by lemma B.4 and since ((xx) S 7, the line search accepts the unit step. Moreover, by strict convexity, EQ 0 for all l E I (0, 1), then Emin E (0, 1] Hence, by
Next, in phase
Lemma B.6. since
dkn 2 Emin and A(x+) 7), take
02,N+1
1 1 Emin + x(xx). ((x++1)
(1 n1)
Moreover, by definition 17 = (3 -9-9-4 4e) and since Emi S 1 we obtain -- 1. This shows (1-1)) that ((x++1) < X(xx) and, in particular, this process converges with composite rate for A(x+) S 1) and some
E (0, 3-y5)
Let ε = Then, 1 Ek 1 since (Ek)KeN boundod and & S lim supr EN 1. Using again the result of lemma B.6 together with A(xx) S 1) we have that 1-5etý 1, for all & N. Therefore, A(x++1) < A(xx). [1 (1-1)*
This mcans
- + {(x+) (1-A(xx)) = 0,
lim k0000
which proves the super-linear rate for the low-rank Newton method required.
C Proof of Theorem 3.2
Let be strictly convex self-concordant function. Theorem 3.2
considers the sequence generated by
X = X5: - PR Q R. (xx)).
For the purposes of this section denote d == - P& P V((x)) where Q is defined in equation 13
and P& in definition Below we will make use of the approximate decrements {(x) and X(x) defined in equation 14. respectively. Before proving Lemma 3.1 state upper bound for X(x) which was proved in Tsipinakis & Parpas 2024).
Lemma C.1 (Tsipinaksss & Parpas 2024 ) Let A(xh,b) be the newton decrement in equation 21 For any EEN it holds that
i(x),)) I A(xn.)).
(33)
Note that the above result holds with probability one. Below we prove Lemma 3.1 3.1 for convenience.
We state again Lemma
Lemma C.2. For all ls E N we have
T V )dn. {2(xx)
(34)
ahe V2 f(xx )dn. S {2(x))
(35)
2 f(xx)) (dh. - -dn.k)l s Ex
(36)
OK.N T(xx) s )(xx) s x(x+) & X(xx),
(37)
Tkp++1
where éx v ((xx))2 }(xx))
Further, let Assumption 2 hold. Then, there crists Pa E (0, 1] such that é S (1 probability 8.
ONN
1/2x)))(xx) with
ahpll
22

---

Published in Transactions on Machine Learning Research (01/2026)
Proof of Lemma 3.1
Equality equation 34 follows directly from the definitions of d and d, By equa-
tion 34 we also get (xxx)) (xx(xx)) d For inequality equation 35 since QHæ is symmetric positive definite, we have
a V2 ((x))dn, V PP
(((xx))PPPAA) Qn (xx, R. R(()())
= =2" (Ra V (xx) )Q
= (x(x). By construction of Q it holds that where z : QH R& V fxx)) and note that z=z=='=''>
QHe R V2 PA) PA) I max QH (R (Ax
max
0k11
0k.p+1
0kpt2
OK.N
IN
0k11
0kpt1
0kptl
0kptt1
IN
Then, z² Q²QQ
H1 (R (R+ V2xx QH S (xxx and equation 35 is proved. Next,
|[V2 (dh, dx,x))12 dh 2 f(xx) da. 2dh v2 fxxx)dn + ah. V2 f(x))da
- A(xx)² - {(xx}?
where the last inequality follows from equation
21 equation and equation 35 This proves inequality
equation 36 For the bounds in equation 37 we start by proving ((xx) I ((x+).
We have
x(xx)² -xxx)= (RV RR BR
RR Ra Q ((x)) Px. Q ) ( RV
We also have
IN RV - H max (IN
Q R 2
f(x)
Q IN
0k.1
Kkptl
Økpt1
-max
... IN
0k.1
0kpt1
OK,N
0.
Putting this all together, we take )(xx)= S A(xx). On the other hand, with similar arguments we have that
A(xx)² -xxx)=(x=====-) =z= (R. Q IN)z
[k.pt1 -1 xT
(Okpt1-))
-1 )(xx)²
OK,N
OK.N
and thus the lower bound of equation 37 follows. Putting this all together and combining it with Lemma
C.1 equation 37 has been proved. Moreover, assumption 2 and equation 37 imply X(xx) 0. Then there
exists E (0, x(xa) with probability 5. Thus, we obtain
OK,N ((xa) & ((xx) S A(xx),
pk Okptt
(38)
with probability 8. Therefore, by the lower bound of the last result. we directly get ex (1 Ju.. (/2x)) with probability too, which concludes the proof.
Tkall
Lemma C.3. Let R" R be strictly conver celf-concordant function and X P R. fx(x). If ((xn.)) < tn, have that
Xx.
23

---

Published in Transactions on Machine Learning Research (01/2026)
42 f²((x))) (i) f(((x+11) 1 (xx1))(x))²
(ii [((((++1)) (x=()((x)))² V2 f(x+))
Proof. Consider the case we have that
From the upper bound in equation 24 that arises for self-coneordanl functions
(xx++))) 3
v2 sxx))
(1 ts |d llxa) )2
1
V2 v(x(x
(1 1xx(xx))2
where the last inequality holds from equation 35
As for the case
ii
we make use of the lower bound in
equation 24 and thus, by equation 35 we have
v f(x(x+))) = (1 (1x f(xx)
Since, further, f is strictly convex we take
((++1)))) (1 {xxxx)))²
[xx(x=]
which concludes the proof.
The next lemma estimates the sub-optimality gap of the process.
Lemma C.4. Let A(xx) < 1. Then,
w(((x)) s f(xx) - f(xx) s w. (A(xx))
If in addition A(xx) < 0.68 then,
-(x))- - {(xo))
Proof. The first result follows directly from Lemma A.2 Further, it holds that (. S S x S 0.68, and thus we take the second result.
Thus X(xx) can be used as an exit condition for Algorithm one may take
Combining the result above with equation 38
1/2
Ckgtt {(x()) 1/2
f(xx) f(x") S we
WKOKN
and thus, for )(xx)= < Mry 0n.t41
{(xx) be used exit condition for Algorithm 1
Lemma C.5. Let to == I+xx=))
Then, for any l E N we have
f(xx+1) - f(xx) S -w((xx))
where w(x) is defined in equation 22
Proof. By Lemma A.1 we have that
f(xx + (odn.)
S fxxx) + (V/(xe), Xh+1 xx) + wa ([xx+11 - xxolxx))
= f(xx) + tu (xxx) Tan (ttll//dn)))))))
< f(xx) txx(x)) ww (txA()))
24

---

Published in Transactions on Machine Learning Research (01/2026)
where for the last inequality we used equation 35 the monotonicity of (x) and 1(x)² (xx(x)) "ai,k Then, using the definition of l we take
)(xx)2 + w. I + {(xen.))
{(XCE)
f(x: + txdn,)) S f(xx))
1 1 )(xx,1))
f((x)) {(x+) log (1+1(x+))
which concludes the proof.
Lemma C.6. Let A(xx) < 1 and set th 1. Further, suppose that assumplion holds. Then
A(Xa-1) VI VI- ENFÉ 1(x1), A(xx);
- a(xx))²
with probability 6.
Proof. By Lemma
C33
inequality equation 37 and la have
V ('I(x++)) (1 - )xxx)) V (xxx;
(39)
with probability one. Now, by the definition of the Newton decrement we have that
X(x+11)= / 1/2
and combining this with inequality equation 39 we take
A(x++1) 1 A(x+)
| (/x))
(40)
Denote Z = f(xx+ Using the fact that
v((x++1) 1 V2 + (ddn))dn. dy + (x(x))
we that
z = [Vx((( 1/2 (f'v2 (22((x )dn dy V [((xx)
[ [(x)) 1/22 f(x: + v dy f(xx))
1/2 2d -xx2f(xx)) 1/2 dn..
where da,k is the Newton direction. Next adding and subtracting the quantity TIV? f(xx)] dr. we have that
Z T dr.x) + (T INxN [VxN)(V2 d
and thus
|z|| <
((x)) dhi)) +
T INXN) ) 42 fxx: 1/2 d
(41)
2,
Using again equation 39 and since, by assumption, yd(xx)) < 1 we take
[xx)) -xxxx)))11/222 1/2025(xx + [V2 (1- v((xx))) INNN-
25

---

Published in Transactions on Machine Learning Research (01/2026)
We now in position to estimate both norms in equation 41
For the first one have that
I. v)(xx)) dy
| -
[Vx)(x)) (dn. - dnw))11
1
|V2 [xxx)]
(an. du 1
1 {(x+)
)(xx)² )2- - )(xx)²
{(X1,k)
Next, the second norm implies
Zz v
xx(x)xxx))2 - dy
|(x](xx)]
X(xx)
1 - ((xx) 1(xx).
Putting this all logether, inequality equalion
40 becomes
((x++1) s V A(xx)² {(xa)²
A(xx)
+
(1 - )(xx)) A(xx).
(1- -=(xxx))
Therefore, by lemma C.2 we take
v EkPl
1(x+)
A(x++1)
(1- -xxx))))) )(xx):
(1 (1- )xx+))
with probability
Proof of theorem 3.2
From eq.
38
we have that
JK.N X(xx) S )(xx).
Mmin v V (xxx) S /
0kpt+1
with probability d for each & E N. Using this inequality, the result lemma C.5 and the monotonicity of w we take
f(x+11) fxxx) S -w (Hmiu v (Ex)) S -- (Mmin v Emin
where the second inequality uses the assumption A(x) 9) The above bound holds with probability 8, and, further, = w (Mmin /Emian) > 0 as required.
Additionally, if X(x)= < f/, then from lemma
C.6 we get that
A(x++1) < vI II Emin/ min ((x+) ( 5
V Emink min + f ((xx);
(1>((x)))
(1 - 1))
with probability 8. Moreover, let E E 10, and (0,1). 1). Then Emin, E (0,1). By n ;= 3-V5+4E that n E (0, 3- Then, from the definition of 9, it holds
we have
2
Emin/ + t
1,
and thus X(xx+1) < A(xB). Therefore, this process converges with composite rate and a probability 8 each &s E N.
26

---

Published in Transactions on Machine Learning Research (01/2026)
If, the other hand, E μ then from lemma C.6 we have that
((x++1) / Expl + f {(xx),
(1 1))
with probability 8. Since / - EXP E (0. and using the definition of i as above, we obtain X(x++1) < ((xx). Hence, since ē = μr 1, then both (Ex) and (Jk)keN converge to and thus
M(x++1) 1(x+) < V1 VI-EMH ((xa)
11 + A(xx) O. (1 {(xx))²
and lim lo-ao
(1 {(x+))
Therefore, we conclude that the process converges with super- Seper-linear rate and probability 8 at each & e
N.
D Proofs of Results of Section 3.3
Proof of Theorem 3.5 Let ((x) be the the approximate decrement based on the current coarse model
[xxx)) In 1/2
{(xx)=
Q H are defined in equation 16
The following results will be used later in the proof. By construction of Q H.kl and assumption 3 for all & E N we have
vI, 3 lokpti L, 3 1QH = MI.,
where the left-hand side inequality holds by This implies
d PA Q H. R& vf(x)) {(xx)
(42)
where, here, for simplicity, we use d instead of | in equation
We also obtain
1 M
{(x))²
R. V((x)))
[VS(x)]1?
(43)
M
with probability >> O. Furthermore, by assumption 3 we take
M
f(x++1) f(xx) < (xxx)
1 12lldn.ll3
43
Mw?
4 1
2v 1/
)(xx))
$ 2M {(xx)
v
2
where the last inequality follows from our assumption to lg. Finally, using eq. 43 into the last inequality we get the desired result with probability 8 8>0 0.
Proof of Theorem 3.4 Combining the result of theorem 3.3 and that for all N we have that lok.pt1| 2
obtain
a3v [
f(x++1) f(xx) <
2 M
with probability 8. Using assumption 5 we get
N. (f(x)) f(x')).
f(xx+1) f(xx) <
WPM
The result. follows by adding and subtracting f(x") in the above inequality.
27

---

Published in Transactions on Machine Learning Research (01/2026)
Datascts
Problem
171
72
NV
CTTslices
Non-linear least-squares
53, 500
385
0.5m
60
CovType
Non-linear least-squares
581,012
54
0.7n
Gisctte
Non-lincar least-squares
6, 000
5, 5,000
0.5m
350
W8T
Non-linear least-squares
14, 951
300
0.5m
60
Table Set-up and dataset details for non-convex, non-linear regression problem.
E Extra Numerical Results and Details
All the datasets used in the experiments can be found in httpss//www cwwe. ntu. edu. tw/-cjlin/ abbswttools/datasettae and htttp//arccive. ics uci, edu/ml/index, php What follows is a description of the algorithms used in comparisons.
1. Gradient Descent (GD) with Armijo-rule
2. Accelerated Gradient Descent (AccGI) with Armijo-rule and momentum 0.5.
3. NewSamp with Armijo-rule (Erdogdu & Montanari 2015
Adam with 4 = 10-3,,40 81 = 0.9,,82 = 0.99, = 10-8 Kingma & Ba 2014)
5. Sigma with Armijo-rule (Tsipinakis & Parpas 2024
SigmaSVD with Armijo-rule
Cubic Regularization of the Newton's method (Cubic Newton) with line search (Neserovv et al. 20188)
To be fair in the comparisons we perform the same step size stratogy for all algorithms but Adam. Noto that NewSamp forms the Hessian matrix as in equation 10 and although it not directly comparable to our approach, we include it in our experiments show that our approach outperforms sub-sampled Newton methods. Similar to N, we denote |Smll the number of samples that NewSamp uses form the Hessian. In addition, the number of the cigenvalues is selected the same for both SigmaSVD and NewSamp in all the experiments. The results of this section were generated on standard laptop machine with Intel 17-10750H CPU processor. The numerical results of main text were generated on GPU processor in Google's colab, The code for generating all the numerical results will be uploaded on github in the near future.
E.1 Non-linear regression
In this section we revisit the non-convex problem of Section g to perform simulations on different datasets. Full details on the datasets used for the non-linear least-squares problem and how the parameters and p are selected are given in Table ? Figure 3 shows that SigmaSVD, Adam and Cubic Newton are the best algorithms as they return the lowest training errors, yet SigmaSVD does better in problems with several saddle points or flat areas. The three algorithms perform similar in terms of GPU time when the problem dimensions is small. However, for problems with n large the results are favourable for SigmaSVD.
E.2 Linelinear regression
Here we minimize the following self-concordant function,
= arg mins(x = arg min xER"
log( - a. x)
x
We generate two datasets illustrate the efficiency of the proposed method for the regimes m > and follows: {a mm is generated from the multivariate Gaussian distribution with zero mean and unit
28

---

Published in Transactions on Machine Learning Research (01/2026)
SigmaSVD
Sigma5VD
ACCED
Rec013
ADAM
ADMM
Cublic Newton
Cubiic Nowiton
22-00-1-
f(8)
19
16
30
10
GPU time (snc)
GPU time (sec)
(a) CovType
(b) CtSlices
SigmasVI
RccGD
ADAM
10-1
Cubli Newton
7:8]
3x010
Sigma5VD
Acecio
2xx
ADAM
Cubic Newton
CPU time (sec)
GPU time (sec)
(c) w8t
(d) Cov'T
SigmasVD
Siqma5VD
AccGD
AccGO
ADAM
ADAM
Newton
Cubic Newton
10-1
0-1
(00
1(0
10
GPU time (sec)
GPU time (s0x)
(e) CtSlices
w8t
Figure 3: Non-convex minimization. All methods in plots from (a) to (c) are initialized at the origin while from (d) to (f) the initializer is selected randomly from a Gaussian A(0,1 1).
29

---

Published in Transactions on Machine Learning Research (01/2026)
Datascts
Problem
7
71
N
[5m
Synthetic
Log linear model
10, 000
1,000
0.5m
0.3m
150
Synthetic
Log linear
1, 000
10, 000
0.1m
150
Table 3: Set-up and dataset details for Log-liness regression.
SignaSVD
10
SigmaSvD
Signa
Sigma
AccGD
10"
ACcGD
18
10
1.0
:0
-0-3
20
CPU time (smm)
CPU time (sre)
(b)
Sigmasvi
10
Sigmasvo
Sigma
Sigmia
ApcCi
10
ACCOD
104
NewSamp
11F
-0-1
10-1
-0-3
10-1
20-1
20-9
200
400
3000
108
200
3101
iterations
iterations
(c) n
(d) m m
Figure 4: Log Linear Regression. Plots (a) and (c) show comparisons between the optimization algorithms for the regime m while (b) and (d) for the regime
variance and (m from uniform distribution. Full details and the set-up for algorithm are given in Table
Figure 4 shows the performance of the optimization algorithms. Due to the domain of this problem, the performance of Adam is missing. The performance of NewSamp is missing from the experiment for the regime n > m due to the large value of n. Clearly, the fastest method is Sigma which achieves very fast super-linear convergence rate. As shown in Theorem 3.2 the reason for the difference in the performance between Sigma and SigmaSVD is that the ratio opt1/o,, /o is large which indicates that performing SVD on the reduced Hcssian matrix we discard much of the important sccond-order information of the problem. However, even in such scenario, in contrast to NewSamp and gradient descent methods, SigmaSVD is able to converge to the global minimum with super-linear rate.
30

---

Published in Transactions on Machine Learning Research (01/2026)
Datascts
Problem
772
71
N
[5ml
p
CTTslices
Logistic model
53, 500
385
0.5m
0.3m
60
10
Leukemia
Logistic model
38
129
0.1n
0.7m
150
10
News20
Logistic model
19, 996
1,355, 129
0.001n
100
10 12
Table 4: Set-up and dataset details for Logistic regression.
Datascts
Problem
772
n
N
Sm
p
CovType
SVM
581, 012
54
0.5n
0.3m
10-2
W8T
SVM
14, 951
300
0.5n
0.5m
60
10a
Table 5: Set-up and dataset details for Support Vector Machines.
E.3 Logistic regression
In this section we concerned with the problem of finding the maximum likelihood in generalized linear models. In particular, consider the regularized logistic regression which as optimization problem takes the following form,
x = arg minf(x) = arg min
Σ 2=1
log(1 + dia Xx
2 |
XER
xER=
Note that the logistic model a strongly convex function with Lipschitz continuous Hcssian matrix. The details and results for the logistic regression experiment are given in Table and Figure 5 respectively. Note that for leukemia dataset, m is too small to perform batch learning and thus the performance of Adam is omitted for this example. Second-order methods clearly outperform first-order methods for the CtSlices dataset. Further, NewSamp is slightly faster than the multilevel methods but obtains slightly worse convergence rate (Figures 5ab and 5d respectively). However, expected, the efficiency of NewSamp reduces drastically for large values o (Leukemia dataset). Note that in this example the multilevel methods perform similarly which indicates that there no need to use the full spectrum of the reduced Hcssian matrix. Last, Figures 5c and 5i show the efficiency of SigmaSVD on a problem with over a million parameters. Such a problem lies at the heart of large scale machine learning and deep learning. Even when the coarse dimensions are selected very small (N 0.001 in this experiment), Figure shows that multilevel methods capable of returning solution with very good accuracy much faster than first-order methods.
E.4
Support Vector Machines
We can train Support Vector Machines (SVMs) using the primal problem over hinge-q loss Huber loss function. Since our approach requires twice differentiable functions consider the hinge-2 loss function for training the SVMs,
m
x arg min xER
x11 + max b,a, x}2
Note that the objective function of this section is convex, however the Hessian matrix is not Lipschitz continuous. The details and results for the logistic regression experiment are given in Table and Figure g respectively. Again, both multilevel methods outperform its counterparts. Here, as also in the logistic regression example, using only a few eigenvalues to form the reduced Hessian matrix does not seem to affect the performance of SigmaSVD. Further, we observed that the performance of NewSamp becomes very poor, even compared to first-order methods, when highly regularized solutions are required (Figures Gb and 6d] On the other hand, both multilevel methods, matter how large small the regularization parameter is chosen, always outperform the first-order methods.
31

---

Published in Transactions on Machine Learning Research (01/2026)
SignaSvD
Sigmasvo
Sigma
Sigma
AcctD
Accii
10
12
NewSamp
NewSamp
(172 /1/(33))
10-1
10-4
10"1
10-9
10-10
2.5
5.6
10.0
12:5
5.0
17.5
70.5
16
15
CPU time (snc)
CPU time (scc)
(a) CtSlices
(b) Leukemia
SigmasVo
10
Sigma
ActGD
10
NewSamp
00:4
10
10-6
50
100
100
200
200
fime (ven))
ilreaiions
(c) News20
(d) CtSlices
SigmasvD
Sigma
ACcGD
NowSamp
((%.))ll Alla
10
10
SigmasvDD Sgma
Arcii
10
ADAM
200
(00)
600
3000
100m
12000
100
350
200
200
Berations
(e)
Leukemia
(f) News20
Figure 5: Logistic Regression. Plots (a) to (c) show the norm of the gradient vs cpu time in seconds while
(d) to (f) show the norm of the gradient vs iterations for three machine learning datasets.
32

---

Published in Transactions on Machine Learning Research (01/2026)
SigmaSub
SigmasvD
Sigma
Sigma
AceGO
Zecra
31-
NewSamp
NewSemm
12
CPU time (s0c)
CPU siman [soc)
(a) W8T
(b) CovTTye
Signation Simn
Sgmatvv
Sigma
AccO0
10-1
.00
200
260
300
156
269
20
Atrratons
(c) W8T
(d) CovType
Figure 6: Support Vector Machines. Plots (a) and (b) show the norm of the gradient epu time in seconds while (c) and (d) norm of the gradient vs iterations for two machine learning datascts.
33


Reviews:
[{"review_id": "ffR1Jvi8b9", "reviewer": "Reviewer_VAxm", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors consider the numerical solution of large-scale optimization problems related to machine learning models. They introduce a novel variant of two-level Newton methods that relies on randomized methods for the computation of approximate second order information related to the reduced Hessian matrix. Exploiting such subspace information allows to reduce the global computational cost of the method. The resulting algorithm named SigmaSVD (Algorithm 1) is presented in Section 2 and analyzed at the theoretical level in Section 3 and Appendices B, C and D, respectively. A super-linear convergence rate is shown for such a low rank Newton method for self-concordant functions. The authors propose a variant applicable to non-convex problems and establish a linear rate of convergence when the Polyak-Lojasiewicz inequality is satisfied. Numerical experiments on various optimization problems related to machine learning are provided in Section 4 and Appendix E showing the efficiency and applicability of SigmaSVD.\n\nStrengths\n========\n- The manuscript is well written and covers both theoretical and numerical aspects.\n\n- Convincing numerical experiments given in Section 4 and Appendix E are provided. These numerical experiments support the claims related to the convergence rate and demonstrate the efficiency of the proposed algorithm. The ability of SigmaSVD to solve large-scale non-convex problems efficiently is also highlighted.\n\nWeaknesses\n============\n- SigmaSVD relies on randomized methods for the computation of approximate second order information related to the reduced Hessian matrix. The authors employ the Nyström method with uniform sampling as a randomized method. Would it be possible to consider more advanced sampling methods that could offer a lower computational complexity for a similar accuracy? The authors are expected to provide more information and details on this aspect.\n\n- The efficiency of the algorithm proposed by the authors is partly based on the fact that the reduced Hessian matrix exhibits a low rank structure for applications related to the optimization of machine learning models. This statement should be highlighted more clearly in the manuscript (e.g. in the Introduction and Conclusion). Reference articles related to this claim should be given and discussed.\n\n- The efficiency of SigmaSVD is partly governed by two hyperparameters: the size of the coarse problem $N$ and the order of the low-rank approximation of the reduced Hessian matrix $p$. The influence of both parameters is shown at the experimental level only in Section 4 and Appendix E, respectively. Comments concerning the determination of the size of the coarse problem are expected at the theoretical level.", "claims_and_evidence_yn": "Yes", "claims_explanation": "- One of the key strengths of the manuscript is that the authors provide a theoretical proof of the convergence rate of their novel algorithm  and also extensive  numerical experiments related to the optimization of machine learning models. Comments related to numerical experiments are convincing and numerical  results support the theoretical claims related to the main theorems of the manuscript.  \n\n- To improve the quality of the manuscript, I have proposed three points to address. The main goal is to explain if the proposed algorithm is of large applicability or not. These are important points to address for increasing the impact of the manuscript.", "audience_yn": "Yes", "audience_explanation": "The authors introduce a novel algorithm for the optimization of large-scale machine learning models. This algorithm is a second order method exhibiting a super linear convergencE rate and is also applicable when solving non-convex problems.\n\nThe manuscript offers two major contributions:\n- a theoretical proof related to the convergence rate of the method (given in  Appendices B and C),\n- detailed numerical experiments on different machine learning models that provide a convincing demonstration of the efficiency of their method (Section 4 and Appendix E).\n\nI am convinced that both contributions would be of interest for TMLR's audience.", "requested_changes": "Critical\n========\n- The authors should provide additional comments on the choice of the restriction and prolongation operators based on the literature related to multilevel methods for optimization problems. \n\n- The authors should explain more carefully how they determine in practice the size of the coarse problem $N$ and the order of the low rank approximation of the reduced Hessian matrix $p$.  \nBoth parameters have a strong impact on the convergence rate of SigmaSVD and detailed comments are expected at the theoretical level in the core of the manuscript. \n\n- Additional comments are required concerning the use of randomized methods as discussed above. See the PhD thesis of Nathan Epperly (Caltech, 2025) \"Make the Most of What You Have: Resource-Efficient Randomized Algorithms for Matrix Computations\" for an accurate overview of the latest contributions in this field. \n\nLess critical\n=========\n- The authors should incorporate in Section 4 or Appendix E comments related to the computational cost of SigmaSVD in practice. It would be of interest to know the different contributions in terms of CPU or GPU time.\n\n- The bibliography should be extended by considering latest contributions in the field. In the current version of the manuscript, no references on multilevel methods for optimization problems dated later than 2023 are provided.\n\n- To improve readability, Figures 1, 3, 4 and 5 should be larger.\n\n- Section 3.2, page 8, please use \"SigmaSVD will converge...\"\n- Section 3.3, page 8, please use \"Polyak-Lojasiewicz\" instead of \"PL\" in the title of this section\n- Section 4.1, page 10, please use \"selected as 0.5.\"\n- Section 4.1, page 10, please use \"In Figures 1a, 1c ...\"\n- Section 4.2, page 12, please use \"one with  N = 2,800...\"\n- Section 4.2, page 12, please use \"another with  N = 1,400...\"\n- Section 4.2, page 12, please use \"as observed in Figure 1c.\"\n- References, please check if arXiv references have been published.\n- Appendix B, page 16, please use \" where $Q_{h,k}^{-1}$ is defined in equation 10\".\n- Proof of Lemma B.1, page 17, please replace at several places $\\nabla f^2(x_k)$ with $\\nabla^2 f(x_k)$.\n- Proof of Lemma B.1, page 17, please use $\\nabla f(x_k)$ instead of $\\nabla f(x_k)^T$ in the last term of the right-hand side in the first equation of the proof.\n- Proof of Lemma B.9, page 19, equation just before relation (32), please use $I_{n \\times n}$ instead of $I_{N \\times N}$.\n- Page 20, just before the proof of Theorem 3.1, please remove \"Recall the\".\n- Page 21, Lemma C.1, please use \"Newton decrement\".\n- Page 23, Lemma C.3, please use \"...that arises...\".\n- Page 24, Lemma C.6, please use $I_{n \\times n}$ instead of $I_{N \\times N}$.\n- Page 25, Appendix D, last relation at the bottom of the page, please use $I_{N \\times N}$ instead of $I_{n}$.\n- Page 26, Appendix E, third item, please modify \"Newton with with\".\n- Page 28, Appendix E, please use \"reduced Hessian matrix\".", "broader_impact_concerns": "No comments required in my opinion.", "additional_comments": "The authors should provide a link to their software repository. It is indeed important to help TMLR's audience to reproduce the numerical experiments provided in Section 4 and Appendix E, respectively. This would be a nice asset that may increase the impact of the manuscript.", "extra_scores": {}, "extra_text": {}}, {"review_id": "plDvUwJPvV", "reviewer": "Reviewer_VqxF", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors introduce the SigmaSVD method applied to smooth functions for unconstrained optimization. They show SigmaSVD to be superlinearly convergent in the convex setting under the assumption of self concordance and in the nonconvex setting under the assumption of self concordance and a PL assumption. The analysis is then extended to the nonconvex and stochastic setting where only a partial gradient and hessian vector product information is available. The proofs are done well, build off of previous work, and provide novel results. \n\nThe computational cost is shown to be less than a full Newton-like method, but the use of a truncated singular value decomposition (T-SVD) still requires considerably more computation than well implemented first order methods. Also, tuning how `coarse' (selecting $N<n$) the T-SVD should be is likely problem dependent. The authors do not discuss these complications in much detail. \n\nWhen applied to a nonlinear least squares, SigmaSVD shows a considerable computational efficiency advantage. In this example, full gradient and hessian vector products are available. The second example applies SigmaSVD to training an autoencoder. Notably, in this example SigmaSVD is modified to incorporate some aspect of momentum inspired by ADAM. No substantive proofs are given for the algorithm under these modifications. Furthermore, a full SVD is used in this example. While the authors say this is reasonable given the size of the coarse model, the numerical results demonstrate a questionable computational advantage of SigmaSVD over ADAM.", "claims_and_evidence_yn": "Yes", "claims_explanation": "The background information, proofs and numerical experiments give clear evidence that SigmaSVD is superlinearly convergent as claimed. They demonstrate that the cost per iteration is considerably cheaper than other Newton-like variants.", "audience_yn": "Yes", "audience_explanation": "While many of the building blocks of this algorithm are not particularly new, in conjunction the SigmaSVD algorithm shows considerable promise. The idea of extracting only enough curvature information to have a computational edge over gradient based methods appears to be a fruitful area of research and should be explored more.", "requested_changes": "My main critique is the lack of discussion surrounding how to choose the size of the T-SVD so that there still is an advantage over gradient based methods. Obviously, tuning was done to select $N$ for the examples given. It is not clear how sensitive the method's advantage is to the size of the T-SVD. Furthermore, one of the examples choose $N=n$ making the T-SVD a full SVD. This may indicate that truncating is not viable in practice.", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "1IkWAohLCk", "reviewer": "Reviewer_ATKW", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper presents a reduced-rank multi-level Newton's method for convex and non-convex objectives. Some simple theoretical properties of the presented methods are established, and numerical results are obtained on MNIST data sets.", "claims_and_evidence_yn": "Yes", "claims_explanation": "descriptions of methods, and imposed assumptions are accurate, proper, and numerical results seem reasonable.", "audience_yn": "No", "audience_explanation": "At this moment, I feel both theoretical and numerical results are too limited to be of broad interest - however, with changes made (see below), the paper has the potential to reach many ML/optimization researchers for its relevance.", "requested_changes": "I have two major comments about results in this paper.\n\n1. On the convex side, the paper assumes the objective function is self-concordant. While this is (somewhat) a standard assumption for Newton's method, I find it a bit concerning: usually, self-concordant functions arise as barrier functions of convex constraints, which are then optimized (using Newton's method) as an intermediate step in interior point methods. In this scenario, it is easy to verify self concordance because it is by construction of barrier functions. In this paper, however, it is directly assumed that the objective function is self concordant: how do we know this is the case? Are common loss functions in ML, e.g. least squares/cross entropy, self-concordant functions? If so, what are their self-concordant parameters (functions of dimension n/N) and how are they related to the properties of the training data?\n\nThe authors must clarify this, preferably with real-world examples in machine learning or statistics, to justify the relevance of proposed method which requires self-concordance of the objective function.\n\n2. On the non-convex side, I have a few comments.\n\n2-(a). I do not find a similar assumption to self concordance when the function is non-convex. The authors state the PL condition, which is well-known to be similar to strong convexity and sufficient for the purpose of first-order methods, but how about self concordance? Without self concordance, is Newton's method very useful, compared simply with fixed step size first-order methods? This also leads me to my second question...\n\n2-(b). The convergence rates established in Theorem 3.4 are linear convergence, not quadratic convergence. Then why should we care about Newton's method? With PL condition, gradient descent achieves linear convergence. Even with stochastic gradients, variants of SVRG with finite sum structures also achieve linear convergence.\n\nOverall, for my second criticism, I think quadratic convergence must be proved for non-convex functions to establish the relevance of Newton's method. If only linear convergence is obtained, then SVRG/Gradient Descent can already do that, with much smaller iteration complexity, on non-convex functions that satisfy the PL condition. If the argument is that Newton's method has a smaller constant in linear convergence, then it should be compared, theoretically, with constants from GD, SVRG, Accelerated GD, to see which one has smaller constants in which scenarios.\n\n3. On the numerical side, MNIST/UCI might be too much a toy example to demonstrate the effectiveness of the proposed method. I recommend the following changes, to make the numerical experiments more interesting.\n\n3-(a). Alongside with MNIST/UCI, the authors should study their methods on more synthetic data using simulations. I have a feelings that synthetic examples might be more interesting than MNIST, because by controlling various parameters of a testing function, different aspects (e.g. dependency on dimension, dependency on self-concordance parameter) could be revealed more clearly.\n\n3-(b). In addition to MNIST/UCI, the authors should report performance of their methods on a more serious ML data set. I am not familiar with those, but I know MNIST is kind of too \"toy\".", "broader_impact_concerns": "N/A", "additional_comments": "N/A", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)