Paper deep dive

Primal Acceleration of Newton's Method

Nikita Doikov

Year: 2026Venue: arXiv preprintArea: math.OCType: PreprintEmbeddings: 60

Intelligence

Status: not_run | Model: - | Prompt: - | Confidence: 0%

Entities (0)

No extracted entities yet.

Relation Signals (0)

No relation signals yet.

Cypher Suggestions (0)

No Cypher suggestions yet.

Abstract

Abstract:We develop a new direct accelerated Newton method for minimizing convex functions with Lipschitz continuous Hessian. The algorithm uses only primal variables and performs just one linear solve per iteration. With a simple predetermined choice of parameters, it achieves the global convergence rate of $O(1/k^3)$ in terms of the functional residual. To the best of our knowledge, this is the first second-order method for this problem class attaining this rate while relying solely on one linear system solve per iteration (without solving auxiliary nonlinear regularized subproblems, such as cubic regularization, performing nonlinear parameter searches, or using dual extragradient corrections). Our method can be implemented in a Hessian-free way, using an inexact linear system solver, while preserving the fast global rate. We further extend our construction to arbitrary geometry through Bregman divergence, and to composite optimization problems.

PDF

Open source PDF →Open local PDF →

Trouble viewing inline? Open PDF directly →

Full Text

57,957 characters extracted from source content.

Expand or collapse full text

Primal Acceleration of Newton’s Method Nikita Doikov Thanks: School of Operations Research and Information Engineering (ORIE), Cornell University, Ithaca, NY, USA. nikita.doikov@cornell.edu August 21, 2026 Abstract We develop a new direct accelerated Newton method for minimizing convex functions with Lipschitz continuous Hessian. The algorithm uses only primal variables and performs just one linear solve per iteration. With a simple predetermined choice of parameters, it achieves the global convergence rate of O⁡(1/k3)O(1/k^3) in terms of the functional residual. To the best of our knowledge, this is the first second-order method for this problem class attaining this rate while relying solely on one linear system solve per iteration (without solving auxiliary nonlinear regularized subproblems, such as cubic regularization, performing nonlinear parameter searches, or using dual extragradient corrections). Our method can be implemented in a Hessian-free way, using an inexact linear system solver, while preserving the fast global rate. We further extend our construction to arbitrary geometry through Bregman divergence, and to composite optimization problems. 1 Introduction Let us consider the standard unconstrained minimization problem: min∈ℝn⁡f⁡(), array[]rcl _ x ^nf( x), array (1) where the objective f:ℝn→ℝf:R^n is convex and several times differentiable. The Fast Gradient Method (FGM) of Nesterov [28] is the optimal first-order algorithm for solving problem (1), which can be implemented in the following momentum form, starting from an arbitrary initialization 0=0∈ℝn x_0= y_0 ^n, we iterate: FGM:k+1=k−αk∇f(k)k+1=k+1+βk(k+1−k),k≥0. .FGM: array[]rcl x_k+1&=& y_k- _k∇ f( y_k)\\[10.0pt] y_k+1&=& x_k+1+ _k( x_k+1- x_k) array \, k≥ 0. (2) Parameter αk>0 _k>0 plays the role of the step-size, and βk∈[0,1] _k∈[0,1] is the momentum coefficient, which specifies how fast we forget the history. The choice βk=0 _k=0 (no momentum) corresponds to the classic gradient descent, while setting, for example βk:=k+3 _k:= kk+3, yields the accelerated rate of convergence of O⁡(1/k2)O(1/k^2) (see Table 1 below for the detailed comparison of rates). In this paper, we study a direct generalization of the classic FGM iterates (2) to second-order optimization methods, incorporating the Hessian (or its approximation) k≈∇2f(k) H_k≈∇^2f( y_k) in the following most natural way, which we call the primal acceleration of Newton method: This paper:k+1=k−(k+1αk)−1∇f(k)k+1=k+1+βk(k+1−k),k≥0. .This paper: array[]rcl x_k+1&=& y_k- ( H_k+ 1 _k I )^-1∇ f( y_k)\\[10.0pt] y_k+1&=& x_k+1+ _k( x_k+1- x_k) array \, k≥ 0. (3) Note that method (3) operates solely in the primal space of variables, without aggregating the dual objects (gradients), in contrast to the estimate sequence framework [30], and without performing any extragradient corrections. Moreover, setting k:= H_k:= 0 immediately recovers the classic FGM. Most importantly, algorithm (3) is very simple to implement as it requires only one linear system solve per iteration. At the same time, using the exact Hessian k:=∇2f(k) H_k:=∇^2f( y_k) and a predefined schedule of second-order step-sizes αk _k equips the algorithm with accelerated O⁡(1/k3)O(1/k^3) rate on convex functions with Lipschitz Hessian. To the best of our knowledge, this is the first second-order method with this rate on our problem class with such simple per-iteration work, as all previously known acceleration techniques required solving auxiliary regularized subproblems (such as cubic regularization) or performing nonlinear parameter searches. We summarize convergence rates for our new approach and for previously known methods in the following table. First-order methods: Hk= H_k= 0 in (3) Method Parameters Convergence Rate, FkF_k Gradient Descent αk=1L1 _k= 1L_1 βk=0 _k=0 O⁡(L1R2k)O ( L_1R^2k ) Fast Gradient Method [28] αk=1L1⋅k+1k+2 _k= 1L_1· k+1k+2 βk=k+3 _k= kk+3 O⁡(L1R2k2)O ( L_1R^2k^2 ) Second-order methods: Hk=∇2f(yk) H_k=∇^2f( y_k) in (3) Method Parameters Convergence Rate, FkF_k Classical Newton’s Method αk=+∞ _k=+∞ βk=0 _k=0 No global convergence Cubic Regularization [33] αk=2L2‖k+1−k‖2 _k= 2L_2\| x_k+1- x_k\|_2 βk=0 _k=0 O⁡(L2D3k2)O ( L_2D^3k^2 ) Gradient Regularization [26] αk=2L2‖∇f(k)‖2 _k= 2 L_2\|∇ f( x_k)\|_2 βk=0 _k=0 O⁡(L2D3k2+e−k)O ( L_2D^3k^2+e^-k ) This paper: αk=13L2R⋅(k+1)(k+2)k+3 _k= 13L_2R· (k+1)(k+2)k+3 βk=k+4 _k= kk+4 O⁡(L2R3k3)O ( L_2R^3k^3 ) Optimal [27, 25, 6] Nonlinear subproblem for each k O⁡(L2R3k7/2)O ( L_2R^3k^7/2 ) Table 1: Summary of optimization methods and their convergence rates in terms of the functional residual Fk:=f⁡(k)−f⋆F_k:=f( x_k)-f , and O⁡(⋅)O(·) hides an absolute numerical constant. We denote by L1L_1 the Lipschitz constant of the gradient, and by L2L_2 the Lipschitz constant of the Hessian, both measured w.r.t. the standard Euclidean norm. We denote by D the diameter of the initial sublevel set, and R≥‖0−⋆‖2R≥\| x_0- x \|_2 is an arbitrary bound for the initial distance to a solution ⋆ x , which we assume to exist. Related Literature. Newton’s method is a powerful and well-developed approach for solving nonlinear optimization problems across various domains. See [35, 34, 30] for the detailed bibliography and for the classic results. One of the fundamental research questions related to the Newton-type methods is the question of its global convergence, starting from an arbitrary initialization. A significant progress in this direction was achieved after the development of the cubic regularization technique with its global convergence guarantees [33]. Adaptive and inexact second-order methods were developed in [8, 9]. Universal algorithms based on the idea of cubic regularization suitable to minimize functions with Hölder continuous Hessian were considered in [22, 17]. Some of the further successful globalization approaches to Newton’s method in the convex case include contracting technique [13], gradient regularization [26, 18, 15]. Notably, the Newton method with gradient regularization possesses the same global O⁡(1/k2)O(1/k^2)-rate as the cubically regularized Newton on convex functions with Lipschitz continuous Hessian (see Table 1). At the same time, this method requires solving only one linear system per iteration, using as the regularization parameter αk _k a power of the gradient norm at the current iteration. In contrast, the cubically regularized Newton method requires solving a univariate nonlinear subproblem of the same complexity as that for trust-region methods [11]. A parallel line of research is devoted to acceleration of the Newton method, with the core initial result of the accelerated cubic Newton [29]. The latter algorithm requires performing one cubic Newton step and one correcting gradient step, which is an auxiliary minimization of the estimating function. Thus, the complexity of each step is of the same cost as for the cubic Newton, while the method achieves an improved O⁡(1/k3)O(1/k^3)-rate on our problem class of convex functions with Lipschitz Hessian. This rate was improved in [27], by replacing the cubic prox function by quadratic one and using inexact proximal-point steps with the total rate of O⁡(1/k7/2)O(1/k^7/2). However, each iteration would require solving an auxiliary line-search subproblem with the total logarithmic cost of second-order oracle calls and non-trivial approximate proximal-point steps. The latter rate was shown to be the optimal one, matching the corresponding lower bounds [4, 2]. Further, using the language of high-order proximal-point steps, it was shown in [31] that the corresponding nonlinear subproblem required to be solved in each iteration of the optimal methods can be made convex, albeit still nonlinear. The optimal second-order scheme from [27] was extended to an arbitrary non-Euclidean geometry in [12]. Notice that the gap between the optimal complexity O⁡(1/ε2/7)O(1/ ^2/7) and the corresponding O⁡(1/ε1/3)O(1/ ^1/3) of the accelerated cubic Newton from [29] is of the factor 1/ε(1/3−2/7)=1/ε1/211/ ^(1/3-2/7)=1/ ^1/21, which seems rather negligible for practical applications and, with a careless implementation, can be offset by the cost of the subproblem. An optimal in terms of the oracle complexity second-order method, omitting an auxiliary logarithm, was developed in [25]. The accelerated scheme in [25] is based on inexact proximal-point view with a predefined optimal rate of convergence, which omits solving an auxiliary nonlinear line-search subproblem as in [27]. At the same time, each iteration might require performing several cubically regularized second-order steps, with amortized optimal complexity. Another approach for obtaining an optimal second-order oracle complexity was proposed in [6] and based on the idea of adaptive Monteiro-Svaiter acceleration from [27]. Note that in both optimal methods from [25] and [6], each step of the accelerated scheme requires solving a certain amortized number of linear systems, as well as performing an auxiliary extragradient step. The idea of the extragradient corrections was further successfully employed in the Extra-Newton method from [3], obtaining an accelerated O⁡(1/k3)O(1/k^3) rate by performing minimization of the quadratic model over the bounded convex set followed by an extragradient correction. In the simplest case of the feasible set provided by the Euclidean ball, the subproblem from [3] is of trust-region type and its cost is the same as that one of the cubically regularized Newton [29]. At the same time, the method [3] admits the use of stochastic gradients and Hessians. The complexity of the accelerated stochastic second-order methods was improved in [1], and for the inexact second-order oracles in [10] based on the inexact proximal-point framework [27]. Therefore, all of the aforementioned acceleration schemes perform an extragradient correction (or an auxiliary estimating function minimization), in order to shift the new proximal center against possible errors caused by approximately solving a certain nonlinear subproblem. Finally, we note that an independent line of work considers more refined problem classes, such as convex functions with Lipschitz or Hölder-continuous third-derivative [21, 15, 23], quasi-self-concordant and generalized-self-concordant function [5, 40, 7, 14, 39], for which it is possible to obtain faster rate beyond those presented in Table 1. Such acceleration comes from more refined smoothness assumptions on the target objective. In this work, we focus on the broad class of convex functions with Lipschitz Hessian, for which our new algorithmic framework seems to be substantially new. Contributions. We propose a very simple second-order acceleration algorithm (3), which, in contrast to previously known schemes, operates solely in the primal space of variables (without extragradient corrections). Our new method can be interpreted as a further development of contracting-proximal acceleration framework from [16]. In contrast to [16], where each iteration requires approximately minimizing an auxiliary contracting subproblem in the logarithmic number of steps of the basic method, we show that it is enough to perform only a single Newton step to ensure that the new point is sufficiently close to the minimizer of the contracted subproblem. In this regard, our viewpoint resembles the classic path-following approach from the theory of the interior-point methods [36, 32, 30, 19, 20]. Another interpretation of our approach is a direct generalization of the fast gradient method in the classic momentum-based form (2). Our new method uses only one linear system solve per iteration, and does not require solving any nonlinear subproblem (as in cubic regularization), performing auxiliary searches, or performing extragradient corrections. At the same time, it achieves the accelerated O⁡(1/k3)O(1/k^3) rate as the original acceleration of the cubic Newton method from [29]. Therefore, to the best of our knowledge, this is the first accelerated second-order method with such a simple implementation, achieving such rate. In Section 2 we show how the momentum-based iterations (3) are related to the more standard prox-center notation used in the literature. For our analysis, the core result about the Newton method which we use is its classic local quadratic convergence, which we briefly review in Section 3. We provide the accelerated analysis of our algorithm and prove the O⁡(1/k3)O(1/k^3) rate in Section 4. In Section 5, we generalize our construction further, showing that we can implement our method for different geometries through the language of Bregman divergence, as well as for composite optimization problems, including simple constraints or non-smooth regularizers, while the corresponding iteration requires solving a nonlinear subproblem. However, we show the same O⁡(1/k3)O(1/k^3) global rate even with approximate solution to the subproblem. In the case of the simplest second-order iterate of the form (3) our result additionally shows that there is no need to solve the linear system exactly. This allows to implement each iteration with efficient large-scale Hessian-free solvers. 2 Algorithm Reparametrization We analyze convergence of method (3) with exact Hessian: k=∇2f(k) H_k=∇^2f( y_k). Therefore, our method would require to solve just one linear system per iteration with the new Hessian computed at the intermediate point. We fix the standard Euclidean norm for vectors ‖2:=⟨,⟩1/2\| x\|_2:= x, x ^1/2, while we show the extension of our analysis to arbitrary norms through Bregman divergence in Section 5. To simplify the reasoning, it is convenient to use the following standard reparameterization. We fix an arbitrary sequence of positive numbers akk≥0\a_k\_k≥ 0 and denote their partial sums, starting from A0:=0A_0:=0, and the normalized ratios by: Ak:⁣=∑i=1kai,γk:=ak+1Ak+1∈(0,1]. array[]rclA_k&:=&Σ _i=1^ka_i, _k\;\;:=\;\; a_k+1A_k+1\;\;∈\;\;(0,1]. array Along with the main sequence of points kk≥0\ x_k\_k≥ 0 we keep track of the auxiliary sequence of prox centers kk≥0\ v_k\_k≥ 0, starting from the same initialization 0=0∈ℝn v_0= x_0 ^n. Then, the prediction points used to collect the oracle information about the function are taken as convex combinations: k:⁣=γkk+(1−γk)k. array[]rcl y_k&:=& _k v_k+(1- _k) x_k. array (4) Our analysis is inspired by the framework of contracting-point methods [16], which suggests to compute the next proximal point k+1 v_k+1 as an inexact minimization of the following subproblem with contracted objective: k+1≈argmin∈ℝn[hk():=Ak+1f(γk+(1−γk)k)+12∥−k∥22]. array[]rcl v_k+1&≈& _ x ^n [\,h_k( x)\;\;:=\;\;A_k+1f ( _k x+(1- _k) x_k )+ 12\| x- v_k\|_2^2\, ]. array (5) Our main observation is that it is enough to perform just a single classical Newton step to compute the next prox center: k+1:⁣=k−∇2hk(k)−1∇hk(k)=k−ak+1(ak+12Ak+1∇2f(k)+)−1∇f(k). array[]rcl v_k+1&:=& v_k-∇^2h_k( v_k)^-1∇ h_k( v_k)\\[10.0pt] &=& v_k-a_k+1 ( a_k+1^2A_k+1∇^2f( y_k)+ I )^-1∇ f( y_k). array (6) It appears that by controlling γk _k, we can ensure that the current point k v_k is in the region of local quadratic convergence of Newton’s method, and one step (6) is enough to ensure a good approximation to the minimum of the contracted subproblem (5). This observation also resembles the classic theory of path-following interior-point methods [30]. After the prox center is updated, we follow the rule of similar triangles to get the next main point: k+1:⁣=γkk+1+(1−γk)k. array[]rcl x_k+1&:=& _k v_k+1+(1- _k) x_k. array (7) It is easy to see that k+1 x_k+1 and k y_k are only different in the displacement of the v-component: k+1−k=(4),(7)γk(k+1−k)=(6)−(∇2f(k)+Ak+1ak+12)−1∇f(k), array[]rcl x_k+1- y_k& ( YkDef),( XkNext)=& _k( v_k+1- v_k)\;\; ( VkNext)=\;\;- (∇^2f( y_k)+ A_k+1a_k+1^2 I )^-1∇ f( y_k), array which gives exactly the first formula in (3) with the regularization parameter equal to αk=ak+12Ak+1. array[]rcl _k&=& a_k+1^2A_k+1. array (8) It remains to obtain the momentum formula for the next prediction point k+1 y_k+1. We notice that k+1:⁣=γk+1k+1+(1−γk+1)k+1=(7)γk+1[k+1γk(k+1−k)]+(1−γk+1)k+1=k+1+γk+1[1γk−1](k+1−k), array[]rcl y_k+1&:=& _k+1 v_k+1+(1- _k+1) x_k+1\;\; ( XkNext)=\;\; _k+1 [ x_k+ 1 _k( x_k+1- x_k) ]+(1- _k+1) x_k+1\\ \\ &=& x_k+1+ _k+1 [ 1 _k-1 ]( x_k+1- x_k), array which gives the prediction update in (3) with the momentum parameter equal to βk=γk+1[1γk−1]=ak+2ak+1⋅AkAk+2. array[]rcl _k&=& _k+1 [ 1 _k-1 ]\;\;=\;\; a_k+2a_k+1· A_kA_k+2. array (9) Therefore, we have a direct correspondence between our new parameterization and the primal form of method (3). 3 Properties of One Newton’s Step Our assumption on the objective f is that it is convex and it has the Lipschitz continuous Hessian with some constant L2>0L_2>0: ‖∇2f()−∇2f()‖≤L2‖−‖2,,∈ℝn, array[]rcl\|∇^2f( y)-∇^2f( x)\|&≤&L_2\| y- x\|_2, x, y ^n, array (10) where in the left hand side of (10) we have the standard spectral norm of a symmetric matrix. Thus, ‖∇f()−∇f()−∇2f()(−)‖2≤L22‖−‖22,,∈ℝn. array[]rcl\|∇ f( y)-∇ f( x)-∇^2f( x)( y- x)\|_2&≤& L_22\| y- x\|_2^2, x, y ^n. array (11) Let us investigate what we can say about one Newton’s step (6) as applied to the contracted objective hk(⋅)h_k(·), defined by (5). We denote by k⋆ v_k the exact minimum of hkh_k: k⋆:⁣=argmin∈ℝn⁡hk(), array[]rcl v_k &:=& _ x ^nh_k( x), array which is unique due to strong convexity of the regularizer; we have ∇2hk()⪰∇^2h_k( v) I for all ∈ℝn v ^n. The optimum satisfies the following stationary condition: ∇hk(k⋆)=ak+1∇f(k+1⋆)+k⋆−k= 0,k+1⋆:⁣=γkk⋆+(1−γk)k. array[]rcl∇ h_k( v_k )&=&a_k+1∇ f( x_k+1 )+ v_k - v_k\;\;=\;\; 0,\\ \\ x_k+1 &:=& _k v_k +(1- _k) x_k. array (12) Therefore, we can relate the distance to the optimum k⋆ v_k after one Newton’s step, as follows: ‖k+1−k⋆‖2≤‖∇2hk(k)(k+1−k⋆)‖2=(12),(6)‖∇hk(k⋆)−∇hk(k)−∇2hk(k)(k⋆−k)‖2≤(11)ℓk2‖k−k⋆‖22, array[]rcl\| v_k+1- v_k \|_2&≤&\|∇^2h_k( v_k)( v_k+1- v_k )\|_2\\ \\ & ( VkStarStat),( VkNext)=&\|∇ h_k( v_k )-∇ h_k( v_k)-∇^2h_k( v_k)( v_k - v_k)\|_2\\ \\ & ( LipHess)≤& _k2\| v_k- v_k \|^2_2, array (13) where ℓk:⁣=ak+13Ak+12L2 array[]rcl _k&:=& a_k+1^3A_k+1^2L_2 array is the Lipschitz constant of the Hessian of the contracted objective. Note that our simple reasoning (13) establishes the classic result about local quadratic convergence of pure Newton’s method. The sufficient condition for the contraction to ensure that the current point k v_k is in the region of local quadratic convergence is ‖k−k⋆‖2≤2ℓk, array[]rcl\| v_k- v_k \|_2&≤& 2 _k, array (14) which we can control by choosing an appropriately small ℓk _k. This choice defines the convergence rate for our accelerated method, as we see in the next section. Together with the distance to the solution, we also need the bound for the new gradient norm: ‖∇hk(k+1)‖2=(6)‖∇hk(k+1)−∇hk(k)−∇2hk(k)(k+1−k)‖2≤(11)ℓk2‖k+1−k‖22. array[]rcl\!\!\!\!\!\|∇ h_k( v_k+1)\|_2&\! ( VkNext)=&\|∇ h_k( v_k+1)-∇ h_k( v_k)-∇^2h_k( v_k)( v_k+1- v_k)\|_2 ( LipHess)≤ _k2\| v_k+1- v_k\|_2^2. array (15) Coupling this bound with an upper bound on the step length: ‖k+1−k‖2≤‖∇hk(k)‖2\| v_k+1- v_k\|_2≤\|∇ h_k( v_k)\|_2, we have the quadratic convergence also in terms of the gradient norm. 4 Accelerated Rate Analysis We are ready to establish the following main convergence result. Theorem 1. Assume that a solution x⋆x to (1) exists and let R≥‖0−⋆‖2R≥\| x_0- x \|_2 be any upper estimate. Let ak+1a_k+1 be chosen such that ak+13Ak+12≤1L2R. array[]rcl a_k+1^3A_k+1^2&≤& 1L_2R. array (16) Then, the following invariant is preserved, for any k≥0k≥ 0: 12‖0−⋆‖22+Akf(⋆)≥12‖k−⋆‖22+Akf(k). array[]rcl 12\| x_0- x \|^2_2+A_kf( x )&≥& 12\| v_k- x \|^2_2+A_kf( x_k). array (17) Proof. We prove (17) by induction. For k=0k=0 it holds as we set A0:=0A_0:=0 and 0=0 x_0= v_0. Assume that (17) holds for a current iterate k≥0k≥ 0. Then, for the next step: 12‖0−⋆‖22+Ak+1f(⋆)=12‖0−⋆‖22+Akf(⋆)+ak+1f(⋆)≥(17)12‖k−⋆‖22+Akf(k)+ak+1f(⋆). array[]rcl 12\| x_0- x \|_2^2+A_k+1f( x )&=& 12\| x_0- x \|_2^2+A_kf( x )+a_k+1f( x )\\ \\ & ( MainInvariant)≥& 12\| v_k- x \|_2^2+A_kf( x_k)+a_k+1f( x ). array (18) By convexity of f, the right hand side can be lower bounded by the contracted objective hk(⋅)h_k(·) at ⋆ x , which is strongly convex. Thus, 12‖0−⋆‖22+Ak+1f(⋆)≥hk(⋆)≥12‖k⋆−⋆‖22+hk(k⋆)=12‖k⋆−⋆‖22+12‖k⋆−k‖22+Ak+1f(k+1⋆)≥12‖k⋆−⋆‖22+12‖k⋆−k‖22+Ak+1f(⋆). array[]rcl 12\| x_0- x \|_2^2+A_k+1f( x )&≥&h_k( x )\;\;≥\;\; 12\| v_k - x \|_2^2+h_k( v_k )\\ \\ &=& 12\| v_k - x \|_2^2+ 12\| v_k - v_k\|_2^2+A_k+1f( x_k+1 )\\ \\ &≥& 12\| v_k - x \|_2^2+ 12\| v_k - v_k\|_2^2+A_k+1f( x ). array We obtain the following useful bound on the distance to the exact minimum k⋆ v_k of the contracted objective: ‖k⋆−⋆‖22+‖k⋆−k‖22≤‖0−⋆‖2≤R2. array[]rcl\| v_k - x \|_2^2+\| v_k - v_k\|_2^2&≤&\| x_0- x \|^2\;\;≤\;\;R^2. array (19) However, in the actual algorithm we do not have access to the exact minimum k⋆ v_k . Instead, we substitute for it one Newton’s step k+1 v_k+1. Applying the convexity lower bound for both function values in (18), f⁡(k)≥f⁡(k+1)+⟨∇f(k+1),k−k+1⟩,f⁡(⋆)≥f⁡(k+1)+⟨∇f(k+1),⋆−k+1⟩, array[]rclf( x_k)&≥&f( x_k+1)+ ∇ f( x_k+1), x_k- x_k+1 ,\\ \\ f( x )&≥&f( x_k+1)+ ∇ f( x_k+1), x - x_k+1 , array we get: 12‖0−⋆‖22+Ak+1f(⋆)≥12‖k−⋆‖22+Ak+1[f⁡(k+1)+⟨∇f(k+1),(γk⋆+(1−γk)k)−k+1⟩]=12‖k−⋆‖22+Ak+1f(k+1)+ak+1⟨∇f(k+1),⋆−k+1⟩. array[]rcl&&\!\!\!\!\!\!\!\!\!\!\!\! 12\| x_0- x \|_2^2+A_k+1f( x )\\ \\ &≥& 12\| v_k- x \|_2^2+A_k+1 [f( x_k+1)+ ∇ f( x_k+1),( _k x +(1- _k) x_k)- x_k+1 ]\\ \\ &=& 12\| v_k- x \|_2^2+A_k+1f( x_k+1)+a_k+1 ∇ f( x_k+1), x - v_k+1 . array Since the norm is Euclidean, we can expand the square as follows: 12‖k−⋆‖22=12‖k+1−⋆‖22+12‖k−k+1‖22+⟨k−k+1,k+1−⋆⟩. array[]rcl 12\| v_k- x \|_2^2&=& 12\| v_k+1- x \|_2^2+ 12\| v_k- v_k+1\|_2^2+ v_k- v_k+1, v_k+1- x . array (20) Substituting this equation we finally arrive at 12‖0−⋆‖22+Ak+1f(⋆)≥12∥k+1−⋆∥22+Ak+1f(k+1)+12∥k−k+1∥22+⟨k+1−k+ak+1∇f(k+1),⋆−k+1⟩=12‖k+1−⋆‖22+Ak+1f(k+1)+12‖k−k+1‖22+⟨∇hk(k+1),⋆−k+1⟩≥(15)12‖k+1−⋆‖22+Ak+1f(k+1)+12‖k−k+1‖22⋅(1−ℓk‖⋆−k+1‖2), array[]rcl&&\!\!\!\!\!\!\!\!\!\!\!\! 12\| x_0- x \|_2^2+A_k+1f( x )\\ \\ &≥& 12\| v_k+1- x \|_2^2+A_k+1f( x_k+1)+ 12\| v_k- v_k+1\|_2^2+ v_k+1- v_k+a_k+1∇ f( x_k+1), x - v_k+1 \\ \\ &=& 12\| v_k+1- x \|_2^2+A_k+1f( x_k+1)+ 12\| v_k- v_k+1\|_2^2+ ∇ h_k( v_k+1), x - v_k+1 \\ \\ & ( GradNormBound)≥& 12\| v_k+1- x \|_2^2+A_k+1f( x_k+1)+ 12\| v_k- v_k+1\|_2^2· (1- _k\| x - v_k+1\|_2 ), array and to establish (17) for the next iterate, it is sufficient to ensure ℓk‖⋆−k+1‖2≤1 _k\| x - v_k+1\|_2≤ 1, which follows from our choice ℓkR=ak+13Ak+12L2R≤ 1, array[]rcl _kR&=& a_k+1^3A_k+1^2L_2R\;\;≤\;\;1, array (21) and the following bound: ℓk‖⋆−k+1‖2≤ℓk(‖⋆−k⋆‖2+‖k+1−k⋆‖)≤(13)ℓk‖⋆−k⋆‖2+ℓk22‖k−k⋆‖22≤(19)max(a,b)∈ℝ2a2+b2≤R2⁡[ℓka+ℓk22b2]≤(21) 1. array[]rcl _k\| x - v_k+1\|_2&≤& _k (\,\| x - v_k \|_2+\| v_k+1- v_k \|\, )\\ \\ & ( NextIterDistance)≤& _k\| x - v_k \|_2+ _k^22\| v_k- v_k \|_2^2\\ \\ & ( BoundContractedMin)≤& _ subarrayc(a,b) ^2\\ a^2+b^2≤ R^2 subarray [ _ka+ _k^22b^2 ]\;\; ( EllKChoice)≤\;\;1. array Therefore, we proved (17) for all k≥0k≥ 0. ∎ From (17), we immediately obtain the global convergence rate: f⁡(k)−f⋆≤R22Akf( x_k)-f ≤ R^22A_k. It remains to choose the controlling coefficients akk≥1\a_k\_k≥ 1 such that (16) is satisfied and that AkA_k is as big as possible. Among various possibilities, we can use the following predetermined choice that leads to simple formulas for our momentum update (3). We set Ak+1=133L2R(k+1)(k+2)(k+3), array[]rclA_k+1&=& 13^3L_2R(k+1)(k+2)(k+3), array and, hence, ak+1=Ak+1−Ak=133L2R(k+1)(k+2)(k+3−k)=132L2R(k+1)(k+2). array[]rcla_k+1\;\;=\;\;A_k+1-A_k&=& 13^3L_2R(k+1)(k+2)(k+3-k)\;\;=\;\; 13^2L_2R(k+1)(k+2). array This choice ensures that condition (16) holds: ak+13Ak+12=1L2R⋅(k+1)3(k+2)3(k+1)2(k+2)2(k+3)2=1L2R⋅(k+1)(k+2)(k+3)2≤1L2R,k≥0. array[]rcl a_k+1^3A_k+1^2&=& 1L_2R· (k+1)^3(k+2)^3(k+1)^2(k+2)^2(k+3)^2\;\;=\;\; 1L_2R· (k+1)(k+2)(k+3)^2\;\;≤\;\; 1L_2R, k≥ 0. array At the same time, we have αk=(8)ak+12Ak+1=13L2R⋅(k+1)(k+2)(k+3)βk=(9)ak+2Akak+1Ak+2=(k+2)(k+3)k(k+1)(k+2)(k+1)(k+2)(k+2)(k+3)(k+4)=k+4. array[]rcl _k& ( AlphaKFormula)=& a_k+1^2A_k+1\;\;=\;\; 13L_2R· (k+1)(k+2)(k+3)\\ \\ _k& ( BetaKFormula)=& a_k+2A_ka_k+1A_k+2\;\;=\;\; (k+2)(k+3)k(k+1)(k+2)(k+1)(k+2)(k+2)(k+3)(k+4)\;\;=\;\; kk+4. array (22) We obtain the following concrete instance of the primal accelerated Newton’s method: k+1=k−(∇2f(k)+3L2Rk+3(k+1)(k+2))−1∇f(k)k+1=k+1+k+4(k+1−k),k≥0. . array[]rcl x_k+1&=& y_k- (∇^2f( y_k)+3L_2R k+3(k+1)(k+2) I )^-1∇ f( y_k)\\ \\ y_k+1&=& x_k+1+ kk+4( x_k+1- x_k) array \, k≥ 0. (23) for which we have the following global rate. Corollary 1. Let R≥‖0−⋆‖2R≥\| x_0- x \|_2 in method (23). Then, we have: f⁡(k)−f⋆≤33L2R32k(k+1)(k+2)=O(L2R3k3),k≥1. array[]rclf( x_k)-f &≤& 3^3L_2R^32k(k+1)(k+2)\;\;=\;\;O ( L_2R^3k^3 ), k≥ 1. array In practice, instead of the rigid choice of parameter L2RL_2R in (23), we might want to choose a tuned second-order stepsize to have the fastest rate for a particular problem. At the same time, from a theoretical perspective, optimization algorithm (23) is the first known to us second-order method with only one linear system solve per iteration and with the global O⁡(1/k3)O(1/k^3) rate of convergence for our problem class of convex functions with Lipschitz continuous Hessian. Instead of using exact solution to the linear system in (23) at each iteration, it is sufficient to have an approximate solution, which we show in the next section. We also further extend method (23) for solving composite optimization problem (allowing to use simple constraints or a non-smooth regularizer in (1)), as well as allow the method to adapt to a different geometry of the problem through Bregman divergence. 5 Bregman Geometry and Composite Problems We show that our approach can be easily generalized to more abstract settings. Let V be a finite-dimensional real vector space with a fixed norm ∥⋅∥\|·\| on it (not necessarily Euclidean). We denote by ∗V^* the dual space, which is the space of linear functions on V. We use the symbol ⟨⋅,⋅⟩ ·,· to denote the value of a linear form ∈∗ s ^* on a vector ∈ x : ⟨,⟩:⁣=⁡()∈ℝ. array[]rcl s, x &:=& s( x)\;\;∈\;\;R. array The dual norm is defined in the standard way: ‖∗:⁣=max∈:‖≤1⟨,⟩,∈∗. array[]rcl\| s\|_*&:=& _ x \;:\;\| x\|≤ 1 s, x , s ^*. array We consider optimization problems in the following composite form: min∈Q[F():⁣=f()+ψ()], array[]rcl _ x∈ Q [\,F( x)&:=&f( x)+ψ( x)\, ], array (24) where ψ is a simple regularizer, dom⁡ψ=Q⊆ ψ=Q , which is a proper closed convex function, possibly non-differentiable. At the same time f is several times differentiable in a neighborhood of every ∈Q x∈ Q. We treat its gradient as the linear form, ∇f()∈∗∇ f( x) ^*, and the Hessian is the linear map ∇2f():→∗∇^2f( x):V ^*, which we assume to be Lipschitz continuous on Q with constant L2>0L_2>0: ‖∇2f()−∇2f()‖:⁣=max∈:‖≤1∥(∇2f()−∇2f())∥∗≤L2‖−‖,,∈Q. array[]rcl\|∇^2f( y)-∇^2f( x)\|&:=& _ h \;:\;\| h\|≤ 1\|(∇^2f( y)-∇^2f( x)) h\|_*\\ \\ &≤&L_2\| y- x\|, x, y∈ Q. array (25) The standard integration argument provides us the bound for the gradient approximation: ‖∇f()−∇f()−∇2f()(−)‖∗≤L22‖−‖2,,∈Q. array[]rcl\|∇ f( y)-∇ f( x)-∇^2f( x)( y- x)\|_*&≤& L_22\| y- x\|^2, x, y∈ Q. array (26) We also fix a 11-strongly convex differentiable distance function d⁡(⋅)d(·). The Bregman divergence is: βd(,):⁣=d⁡()−d⁡()−⟨∇d(),−⟩≥12‖−‖2. array[]rcl _d( x; y)&:=&d( y)-d( x)- ∇ d( x), y- x \;\;≥\;\; 12\| y- x\|^2. array (27) The main algebraic identity which is the key structural property of the Bregman divergence (compare it with (20)) is the next one, for any v, + v^+, and x it holds βd(,)=βd(+,)+βd(,+)+⟨∇d()−∇d(+),+−⟩, array[]rcl _d( v; x)&=& _d( v^+; x)+ _d( v; v^+)+ ∇ d( v)-∇ d( v^+), v^+- x , array (28) which follows directly from the definition. Consider the following general second-order scheme, starting from 0=0∈Q x_0= v_0∈ Q, iterate k≥0k≥ 0: 1. Set k:=γkk+(1−γk)k y_k:= _k v_k+(1- _k) x_k and use k=∇2f(k) H_k=∇^2f( y_k). 2. Compute k+1 v_k+1 such that 1ak+1(∇d(k+1)−∇d(k))+∇f(k)+γkk(k+1−k)+ψ′(k+1)≈∈∗, array[]rcl 1a_k+1 (∇ d( v_k+1)-∇ d( v_k) )+∇ f( y_k)+ _k H_k ( v_k+1- v_k )+ψ ( v_k+1)&≈& 0\;\;∈\;\;V^*, array where ψ′(k+1)∈∂ψ⁡(k+1)ψ ( v_k+1)∈∂ψ( v_k+1) is a subgradient of the regularizer. 3. Set k+1:=γkk+1+(1−γk)k x_k+1:= _k v_k+1+(1- _k) x_k. Note that computation of k+1 v_k+1 corresponds to the approximate minimization of the following subproblem: k+1≈argmin∈Qak+1⟨∇f(k),−k⟩+ak+122Ak+1⟨k(−k),−k⟩+ak+1ψ()+βd(k;). array[]rcl v_k+1&≈& _ x∈ Q \\;a_k+1 ∇ f( y_k), x- v_k + a_k+1^22A_k+1 H_k( x- v_k), x- v_k \\[10.0pt] && \;\;+\;a_k+1ψ( x)+ _d( v_k; x)\; \. array (29) This is minimization of a convex quadratic function, regularized by the Bregman divergence, and augmented by the composite component. Hence, we assume that ψ is simple enough that we can perform this minimization efficiently. Moreover, we need only an inexact solution to (29) with the following simple tolerance condition. We require to compute k+1 v_k+1 which satisfies the inequality: ‖1ak+1(∇d(k+1)−∇d(k))+∇f(k)+γkk(k+1−k)+ψ′(k+1)‖∗≤δk+1,δk+1:=γk2L22min⁡[‖k+1−k‖2,R2], array[]rcl cases& \| 1a_k+1 (∇ d( v_k+1)-∇ d( v_k) )+∇ f( y_k)+ _k H_k( v_k+1- v_k)+ψ ( v_k+1) \|_*\;≤\; _k+1,\\[15.0pt] & _k+1\;:=\; _k^2L_22 [\,\| v_k+1- v_k\|^2,\,R^2\, ], cases array (30) where ψ′(k+1)∈∂ψ⁡(k+1)ψ ( v_k+1)∈∂ψ( v_k+1) is a certain subgradient, and R>0R>0 is a fixed parameter. Note that we can expect from any reasonable numerical method solving subproblem (29) to make the left hand side in the inequality in (30) small, and our condition with δk+1 _k+1 just specifies the required tolerance precisely. Note that it is very important that we admit inexact solution to (29), as in general, this subproblem can be quite complicated due to an arbitrary choice of d and even for simple regularizers ψ. At the same time, the subproblem remains convex, so standard convex optimization techniques can be employed. Clearly, without composite part ψ⁡()≡0ψ( x)≡ 0 and using the Euclidean norm d⁡()=12‖22d( x)= 12\| x\|_2^2, this algorithm gives our previous method with a possibility of using inexact linear solver. The condition (30) becomes: ‖∇f(k)+(γkk+1ak+1)(k+1−k)‖2≤γk2L22min⁡[‖k+1−k‖22,R2], array[]rcl \|∇ f( y_k)+ ( _k H_k+ 1a_k+1 I )( v_k+1- v_k) \|_2&≤& _k^2L_22 [\,\| v_k+1- v_k\|_2^2,\;R^2\, ], array (31) and thus we can apply any Hessian-free method solving the linear system, e.g., the conjugate gradient method, making the left hand side in (31) small. In this case, we do not need to compute the full Hessian. In our solver, we can use only the Hessian-vector products, which is suitable for large-scale models. For our new most general scheme we prove the following convergence result. Theorem 2. Let R≥2βd(0,⋆)R≥ 2 _d( x_0; x ) and let each k+1 v_k+1 satisfy the inexact condition (30). Let ak+1a_k+1 be chosen such that ak+13Ak+12≤25L2R. array[]rcl a_k+1^3A_k+1^2&≤& 25L_2R. array (32) Then, for any k≥0k≥ 0 we have βd(0,⋆)+AkF(⋆)≥βd(k,⋆)+AkF(k). array[]rcl _d( x_0; x )+A_kF( x )&≥& _d( v_k; x )+A_kF( x_k). array (33) Proof. We prove (33) by induction in k≥0k≥ 0. It clearly holds for k=0k=0. Consider any k≥0k≥ 0: βd(0,⋆)+Ak+1F(⋆)=βd(0,⋆)+AkF(⋆)+ak+1F(⋆)≥(33)βd(k,⋆)+AkF(k)+ak+1F(⋆). array[]rcl _d( x_0; x )+A_k+1F( x )&=& _d( x_0; x )+A_kF( x )+a_k+1F( x )\\ \\ & ( NewInvariant)≥& _d( v_k; x )+A_kF( x_k)+a_k+1F( x ). array (34) The right hand side is lower bounded by the following objective with contracted smooth part: hk()=Ak+1f(γk+(1−γk)k)+ak+1ψ()+Akψ(k)+βd(k,). array[]rclh_k( x)&=&A_k+1f ( _k x+(1- _k) x_k )+a_k+1ψ( x)+A_kψ( x_k)+ _d( v_k; x). array Denote the minimum of hk(⋅)h_k(·) by k⋆ v_k , which satisfies the optimality condition (compare with (12)): −ak+1ψ′(k⋆):⁣=ak+1∇f(k+1⋆)+∇d(k⋆)−∇d(k)∈−ak+1∂ψ(k⋆),k+1⋆:⁣=γkk⋆+(1−γk)k. array[]rcl-a_k+1ψ ( v_k )&:=&a_k+1∇ f( x_k+1 )+∇ d( v_k )-∇ d( v_k)\;\;∈\;\;-a_k+1∂ψ( v_k ),\\ \\ x_k+1 &:=& _k v_k +(1- _k) x_k. array (35) Using the key property of the Bregman divergence (28) and convexity of all components, we thus have βd(0,⋆)+Ak+1F(⋆)≥hk(⋆)≥βd(k⋆,⋆)+k(k⋆)≥Ak+1F(k+1⋆)+βd(k,k⋆)+βd(k⋆,⋆), array[]rcl _d( x_0; x )+A_k+1F( x )&≥&h_k( x )\;\;≥\;\; _d( v_k ; x )+ h_k( v_k )\\ \\ &≥&A_k+1F( x_k+1 )+ _d( v_k; v_k )+ _d( v_k ; x ), array where we use the same notation as before, k+1⋆:=γkk⋆+(1−γk)k x_k+1 := _k v_k +(1- _k) x_k. Hence, we obtain the following generalization of (19): βd(k⋆,⋆)+βd(k,k⋆)≤βd(0,⋆). array[]rcl _d( v_k ; x )+ _d( v_k; v_k )&≤& _d( x_0; x ). array (36) It is remarkable that the reasoning works for any choice of the convex distance function d. Taking into account that d is strongly convex (27) with respect to a fixed primal norm, and using the bound for the initial distance, we obtain: ‖k⋆−⋆‖2+‖k⋆−k‖2≤R2. array[]rcl\| v_k - x \|^2+\| v_k - v_k\|^2&≤&R^2. array (37) We can continue the reasoning as follows, using the convexity of f and ψ: βd(0,⋆)+Ak+1F(⋆)≥(34)βd(k,⋆)+ak+1F(⋆)+AkF(k)≥βd(k,⋆)+ak+1ψ(∗)+Akψ(k)+Ak+1f(k+1)+ak+1⟨∇f(k+1),⋆−k+1⟩≥βd(k,⋆)+ak+1ψ(k+1)+Akψ(k)+Ak+1f(k+1)+ak+1⟨∇f(k+1)+ψ′(k+1),⋆−k+1⟩≥βd(k,⋆)+Ak+1F(k+1)+ak+1⟨∇f(k+1)+ψ′(k+1),⋆−k+1⟩. array[]rcl&&\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! _d( x_0; x )+A_k+1F( x )\\ \\ & ( GenProof1)≥& _d( v_k; x )+a_k+1F( x )+A_kF( x_k)\\ \\ &≥& _d( v_k; x )+a_k+1ψ( x^*)+A_kψ( x_k)+A_k+1f( x_k+1)+a_k+1 ∇ f( x_k+1), x - v_k+1 \\ \\ &≥& _d( v_k; x )+a_k+1ψ( v_k+1)+A_kψ( x_k)+A_k+1f( x_k+1)\\ \\ && +\;a_k+1 ∇ f( x_k+1)+ψ ( v_k+1), x - v_k+1 \\ \\ &≥& _d( v_k; x )+A_k+1F( x_k+1)+a_k+1 ∇ f( x_k+1)+ψ ( v_k+1), x - v_k+1 . array Now, we use the key identity (28), which explains the structure of the Bregman divergence: βd(k,⋆)=βd(k+1,⋆)+βd(k,k+1)+⟨∇d(k)−∇d(k+1),k+1−⋆⟩. array[]rcl _d( v_k; x )&=& _d( v_k+1; x )+ _d( v_k; v_k+1)+ ∇ d( v_k)-∇ d( v_k+1), v_k+1- x . array Hence, βd(0,⋆)+Ak+1F(⋆)≥βd(k+1,⋆)+Ak+1F(k+1)+βd(k,k+1)+ak+1⟨∇f(k+1)+ψ′(k+1)+1ak+1(∇d(k+1)−∇d(k)),⋆−k+1⟩≥(30)βd(k+1,⋆)+Ak+1F(k+1)+βd(k,k+1)−ℓk2‖k+1−k‖2⋅‖k+1−⋆‖+ak+1⟨∇f(k+1)−∇f(k)−k(k+1−k),⋆−k+1⟩≥(26)βd(k+1,⋆)+Ak+1F(k+1)+βd(k,k+1)−ℓk‖k+1−k‖2⋅‖k+1−⋆‖≥(27)βd(k+1,⋆)+Ak+1F(k+1)+12‖k−k+1‖2(1−2ℓk‖k+1−⋆‖) array[]rcl&&\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! _d( x_0; x )+A_k+1F( x )\\ \\ &≥& _d( v_k+1; x )+A_k+1F( x_k+1)+ _d( v_k; v_k+1)\\ \\ && +\;a_k+1 ∇ f( x_k+1)+ψ ( v_k+1)+ 1a_k+1(∇ d( v_k+1)-∇ d( v_k)), x - v_k+1 \\ \\ & ( DeltaKDef)≥& _d( v_k+1; x )+A_k+1F( x_k+1)+ _d( v_k; v_k+1)- _k2\| v_k+1- v_k\|^2·\| v_k+1- x \|\\ \\ && +\;a_k+1 ∇ f( x_k+1)-∇ f( y_k)- H_k( x_k+1- y_k), x - v_k+1 \\ \\ & ( GradApprox)≥& _d( v_k+1; x )+A_k+1F( x_k+1)+ _d( v_k; v_k+1)- _k\| v_k+1- v_k\|^2·\| v_k+1- x \|\\ \\ & ( Bregman)≥& _d( v_k+1; x )+A_k+1F( x_k+1)+ 12\| v_k- v_k+1\|^2 (1-2 _k\| v_k+1- x \| ) array where ℓk=ak+13L2Ak+12 _k= a_k+1^3L_2A_k+1^2. Therefore, to complete the proof we need to check the following inequality: ℓk‖k+1−⋆‖≤12. array[]rcl _k\| v_k+1- x \|&≤& 12. array (38) Let us observe the main properties of the inexact composite Newton step equipped with the condition (30). Due to strong convexity of the distance function, and using convexity of the objective, so that k⪰ H_k 0, we have ‖k+1−k⋆‖2≤⟨∇d(k+1)−∇d(k⋆),k+1−k⋆⟩≤⟨∇d(k+1)−∇d(k⋆)+γkak+1k(k+1−k⋆),k+1−k⋆⟩=(35)⟨∇d(k+1)−∇d(k)+ak+1∇f(k)+γkak+1k(k+1−k)+ak+1ψ′(k+1),k+1−k⋆⟩+ak+1⟨∇f(k+1⋆)−∇f(k)−γkk(k⋆−k),k+1−k⋆⟩+ak+1⟨ψ′(k⋆)−ψ′(k+1),k+1−k⋆⟩⏟≤0≤(30),(26)(ak+1δk+1+ℓk2‖k−k⋆‖2)‖k+1−k⋆‖. array[]rcl&&\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\| v_k+1- v_k \|^2\\ \\ &≤& ∇ d( v_k+1)-∇ d( v_k ), v_k+1- v_k \\ \\ &≤& ∇ d( v_k+1)-∇ d( v_k )+ _ka_k+1 H_k( v_k+1- v_k ), v_k+1- v_k \\ \\ & ( BregmanOptimality)=& ∇ d( v_k+1)-∇ d( v_k)+a_k+1∇ f( y_k)+ _ka_k+1 H_k( v_k+1- v_k)\\ \\ && +\;a_k+1ψ ( v_k+1), v_k+1- v_k \\ \\ &&+\;a_k+1 ∇ f( x_k+1 )-∇ f( y_k)- _k H_k( v_k - v_k), v_k+1- v_k \\ \\ &&+\; a_k+1 ψ ( v_k )-ψ ( v_k+1), v_k+1- v_k _≤ 0\\ \\ & ( DeltaKDef),( GradApprox)≤& (a_k+1 _k+1+ _k2\| v_k- v_k \|^2 )\| v_k+1- v_k \|. array Therefore, we justified that ‖k+1−k⋆‖≤ak+1δk+1+ℓk2‖k−k⋆‖2≤(30)ℓk2‖k−k⋆‖2+ℓk2R2. array[]rcl\| v_k+1- v_k \|&≤&a_k+1 _k+1+ _k2\| v_k- v_k \|^2\;\; ( DeltaKDef)≤\;\; _k2\| v_k- v_k \|^2+ _k2R^2. array (39) It remains to notice that ℓk‖k+1−⋆‖≤ℓk‖k⋆−⋆‖+ℓk‖k+1−k⋆‖≤(39)ℓk‖k⋆−⋆‖+ℓk22‖k−k⋆‖2+ℓk22R2≤(37)max(a,b)∈ℝ2a2+b2≤R2⁡[ℓka+ℓk22b2]+ℓk22R2≤ℓkR+ℓk22R2≤12, array[]rcl _k\| v_k+1- x \|&≤& _k\| v_k - x \|+ _k\| v_k+1- v_k \|\\ \\ & ( VkNewBregmanBound)≤& _k\| v_k - x \|+ _k^22\| v_k- v_k \|^2+ _k^22R^2\\ \\ & ( BoundGenMin2)≤& _ subarrayc(a,b) ^2\\ a^2+b^2≤ R^2 subarray [ _ka+ _k^22b^2 ]+ _k^22R^2\;\;≤\;\; _kR+ _k^22R^2\;\;≤\;\; 12, array where we used our selection of the parameters (32). ∎ Corollary 2. We see that up to a numerical constant in the right hand side, condition (32) coincides with that one from Theorem 1. For example, we can set Ak+1=2335L2R(k+1)(k+2)(k+3),ak+1=2325L2R(k+1)(k+2) array[]rclA_k+1&=& 23^35L_2R(k+1)(k+2)(k+3), a_k+1\;\;=\;\; 23^25L_2R(k+1)(k+2) array and γk=ak+1Ak+1=3k+3 _k= a_k+1A_k+1= 3k+3 in our general algorithm, which gives the following accelerated rate, for k≥1k≥ 1: F⁡(k)−F⋆≤(33)βd(0,⋆)Ak≤335L2R3k(k+1)(k+2)=O⁡(L2R3k3). array[]rclF( x_k)-F & ( NewInvariant)≤& _d( x_0; x )A_k\;\;≤\;\; 3^35L_2R^3k(k+1)(k+2)\;\;=\;\;O ( L_2R^3k^3 ). array 6 Discussion In this work, we show that the direct inclusion of second-order information into the classic fast gradient method in the momentum form (3) achieves an improved convergence rate of O⁡(1/k3)O(1/k^3) for convex functions with Lipschitz continuous Hessian. This is the same rate as for the accelerated cubic Newton from [29]. However, in our approach we use only one (possibly inexact) linear solve per iteration, and do not require performing correcting extragradient steps. To the best of our knowledge, this is the first second-order method with such simple iteration complexity for this problem class. One of the main open questions is how to perform an adaptive estimation of the parameters from (23) such as L2L_2 and R. While currently we can treat their product γ≡L2Rγ≡ L_2R as a single scalar second-order step-size that we fine-tune, it would be interesting to develop an adaptive second-order scheme with only one linear solve per iteration. It would be also interesting to explore optimal second-order schemes with an improved O⁡(1/k7/2)O(1/k^7/2) rate and simple complexity per iteration. This seems to require replacing our global parameter R by a local estimate of the proximal step, which immediately leads to complications of the subproblem, which is no longer linear. Finally, it would be interesting to establish connections to classic quasi-Newton methods, by combining our analysis with recent developments of their non-asymptotic analysis [38, 37, 24]. We keep these directions for future research. Acknowledgment We are very grateful to Yurii Nesterov and Anton Rodomanov for useful comments on a preliminary version of this manuscript, which improved the presentation. References [1] Artem Agafonov, Dmitry Kamzolov, Alexander Gasnikov, Ali Kavis, Kimon Antonakopoulos, Volkan Cevher, and Martin Takáč. Advancing the lower bounds: An accelerated, stochastic, second-order method with optimal adaptation to inexactness. arXiv preprint arXiv:2309.01570, 2023. [2] Naman Agarwal and Elad Hazan. Lower bounds for higher-order convex optimization. In Conference On Learning Theory, pages 774–792. PMLR, 2018. [3] Kimon Antonakopoulos, Ali Kavis, and Volkan Cevher. Extra-newton: A first approach to noise-adaptive accelerated second-order methods. Advances in Neural Information Processing Systems, 35:29859–29872, 2022. [4] Yossi Arjevani, Ohad Shamir, and Ron Shiff. Oracle complexity of second-order methods for smooth convex optimization. Mathematical Programming, 178:327–360, 2019. [5] Francis Bach. Self-concordant analysis for logistic regression. 2010. [6] Yair Carmon, Danielle Hausler, Arun Jambulapati, Yujia Jin, and Aaron Sidford. Optimal and adaptive Monteiro-Svaiter acceleration. Advances in Neural Information Processing Systems, 35:20338–20350, 2022. [7] Yair Carmon, Arun Jambulapati, Qijia Jiang, Yujia Jin, Yin Tat Lee, Aaron Sidford, and Kevin Tian. Acceleration with a ball optimization oracle. Advances in Neural Information Processing Systems, 33:19052–19063, 2020. [8] Coralia Cartis, Nicholas IM Gould, and Philippe L Toint. Adaptive cubic regularisation methods for unconstrained optimization. part i: motivation, convergence and numerical results. Mathematical Programming, 127(2):245–295, 2011. [9] Coralia Cartis, Nicholas IM Gould, and Philippe L Toint. Adaptive cubic regularisation methods for unconstrained optimization. part i: worst-case function-and derivative-evaluation complexity. Mathematical programming, 130(2):295–319, 2011. [10] Lesi Chen, Chengchang Liu, Luo Luo, John Lui, and Jingzhao Zhang. Optimal convex optimization with inexact second-order oracles. arXiv preprint arXiv:2607.24520, 2026. [11] Andrew R Conn, Nicholas IM Gould, and Philippe L Toint. Trust region methods. SIAM, 2000. [12] Juan Pablo Contreras, Cristóbal Guzmán, and David Martínez-Rubio. Non-euclidean high-order smooth convex optimization. arXiv preprint arXiv:2411.08987, 2024. [13] Nikita Doikov. New second-order and tensor methods in Convex Optimization. PhD thesis, Ph. D. thesis, Université catholique de Louvain, 2021. [14] Nikita Doikov. Minimizing quasi-self-concordant functions by gradient regularization of newton method: N. doikov. Mathematical Programming, pages 1–39, 2025. [15] Nikita Doikov, Konstantin Mishchenko, and Yurii Nesterov. Super-universal regularized Newton method. SIAM Journal on Optimization, 34(1):27–56, 2024. [16] Nikita Doikov and Yurii Nesterov. Contracting proximal methods for smooth convex optimization. SIAM Journal on Optimization, 30(4):3146–3169, 2020. [17] Nikita Doikov and Yurii Nesterov. Minimizing uniformly convex functions by cubic regularization of Newton method. Journal of Optimization Theory and Applications, pages 1–23, 2021. [18] Nikita Doikov and Yurii Nesterov. Gradient regularization of Newton method with Bregman distances. Mathematical Programming, pages 1–25, 2023. [19] Pavel Dvurechensky and Yurii Nesterov. Global performance guarantees of second-order methods for unconstrained convex minimization. 2018. [20] Pavel Dvurechensky and Yurii Nesterov. Improved global performance guarantees of second-order methods in convex minimization. Foundations of Computational Mathematics, pages 1–41, 2025. [21] Geovani Nunes Grapiglia and Yu Nesterov. Tensor methods for minimizing convex functions with hölder continuous higher-order derivatives. SIAM Journal on Optimization, 30(4):2750–2779, 2020. [22] Geovani Nunes Grapiglia and Yurii Nesterov. Regularized Newton methods for minimizing functions with Hölder continuous Hessians. SIAM Journal on Optimization, 27(1):478–506, 2017. [23] Slavomír Hanzely, Farshed Abdukhakimov, and Martin Takáč. Newton method revisited: Global convergence rates up to o (1/kˆ 3) for stepsize schedules and linesearch procedures. 2026. [24] Qiujiang Jin and Aryan Mokhtari. Non-asymptotic superlinear convergence of standard quasi-Newton methods. Mathematical Programming, 200(1):425–473, 2023. [25] Dmitry Kovalev and Alexander Gasnikov. The first optimal acceleration of high-order methods in smooth convex optimization. Advances in Neural Information Processing Systems, 35:35339–35351, 2022. [26] Konstantin Mishchenko. Regularized Newton method with global ⁡(1/k2)O(1/k^2) convergence. SIAM Journal on Optimization, 33(3):1440–1462, 2023. [27] Renato DC Monteiro and Benar F Svaiter. An accelerated hybrid proximal extragradient method for convex optimization and its implications to second-order methods. SIAM Journal on Optimization, 23(2):1092–1125, 2013. [28] Yurii Nesterov. A method for solving the convex programming problem with convergence rate O(1/kˆ2). In Dokl. akad. nauk Sssr, volume 269, pages 543–547, 1983. [29] Yurii Nesterov. Accelerating the cubic regularization of Newton’s method on convex problems. Mathematical Programming, 112(1):159–181, 2008. [30] Yurii Nesterov. Lectures on convex optimization, volume 137. Springer, 2018. [31] Yurii Nesterov. Inexact accelerated high-order proximal-point methods. Mathematical Programming, 197(1):1–26, 2023. [32] Yurii Nesterov and Arkadi Nemirovski. Interior-point polynomial algorithms in convex programming. SIAM, Philadelphia, 1994. [33] Yurii Nesterov and Boris T Polyak. Cubic regularization of Newton method and its global performance. Mathematical programming, 108(1):177–205, 2006. [34] Jorge Nocedal and Stephen J Wright. Numerical optimization. Springer, 2006. [35] Boris T Polyak. Newton’s method and its use in optimization. European Journal of Operational Research, 181(3):1086–1096, 2007. [36] James Renegar. A mathematical view of interior-point methods in convex optimization. SIAM, 2001. [37] Anton Rodomanov. Quasi-Newton methods with provable efficiency guarantees. PhD thesis, UCLouvain, 2022. [38] Anton Rodomanov and Yurii Nesterov. New results on superlinear convergence of classical quasi-Newton methods. Journal of optimization theory and applications, 188(3):744–769, 2021. [39] Andrei Semenov, Martin Jaggi, and Nikita Doikov. Gradient-normalized smoothness for optimization with approximate hessians. In International Conference on Learning Representations, volume 2026, pages 13974–14029, 2026. [40] Tianxiao Sun and Quoc Tran-Dinh. Generalized self-concordant functions: a recipe for newton-type methods. Mathematical Programming, 178(1):145–213, 2019.

Instant research discovery