Paper deep dive
Halpern Iteration Achieves $\tilde{\mathcal{O}}(ε^{-1/p})$ $p$th-Order Oracle Complexity for Monotone Variational Inequalities
Lesi Chen, Xinliang Zhang, Hengyu Wang, Chengchang Liu, Yongchao Chen, Jingzhao Zhang
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 92%
Last extracted: 8/12/2026, 1:41:25 AM
Summary
This paper proposes novel high-order methods for solving smooth monotone variational inequalities (MVI) by combining the Halpern iteration with Newton proximal extragradient (NPE) and a new Anchored Tensor Method (ATM). The proposed Halpern-NPE method achieves a convergence rate of O~(T^-2) for second-order methods, and the generalized Halpern-ATM method achieves O~(T^-p) for p-th order methods, significantly improving upon prior bounds for p >= 2.
Entities (8)
Relation Signals (6)
Halpern-NPE → achievesrate → O~(T^-2)
confidence 95% · Halpern-NPE method that achieves an even faster rate of O~(T^-2) for solving MVIs.
Halpern-ATM → achievesrate → O~(T^-p)
confidence 95% · Halpern-ATM achieves the fast convergence rate of O~(T^-p).
NPE → proposedby → Monteiro and Svaiter
confidence 92% · Monteiro and Svaiter (SIAM J. Optim., 2012) showed that a second-order method, NPE...
Halpern-NPE → combines → Halpern Iteration
confidence 90% · by using a large-step inexact Halpern iteration, we propose a novel Halpern-NPE method
Halpern-ATM → combines → Anchored Tensor Method
confidence 90% · combine it with the Halpern iteration to achieve a faster convergence rate
Halpern-NPE → improvesupon → NPE
confidence 85% · This improves all prior results for p >= 2
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:We study second- and higher-order methods for solving smooth monotone variational inequalities (MVI). Monteiro and Svaiter (SIAM J. Optim., 2012) showed that a second-order method, NPE, converges at the rate of $\mathcal{O}(T^{-1.5})$. For convex-concave minimax optimization, a subset of MVI problems, Chen, Liu, Luo, and Zhang (COLT 2025) recently improved the complexity to $\tilde{\mathcal{O}}( T^{-1.75})$ . However, it is open whether the conjectured complexity for MVI can be improved. In this paper, by using a large-step inexact Halpern iteration, we propose a novel Halpern-NPE method that achieves an even faster rate of $\tilde{\mathcal{O}}(T^{-2})$ for solving MVIs. We also provide the $p$th-order generalization of our method. We first introduce an Anchored Tensor Method (ATM) that achieves the rate of $\mathcal{O}(T^{-(p-1)})$, and then combine it with the Halpern iteration to achieve a faster convergence rate of $\tilde{\mathcal{O}}(T^{-p})$. This improves all prior results for $p \ge 2$ and matches the classical extragradient method for $p=1$.
Tags
Links
- Source: https://arxiv.org/abs/2608.08463v1
- Canonical: https://arxiv.org/abs/2608.08463v1
Trouble viewing inline? Open PDF directly →
Full Text
67,192 characters extracted from source content.
Expand or collapse full text
Halpern Iteration Achieves ~(ϵ−1/p) O(ε^-1/p) ppth-Order Oracle Complexity for Monotone Variational Inequalities Lesi Chen * 1 Xinliang Zhang * 1 Hengyu Wang * 3 4 Chengchang Liu 5 Yongchao Chen 2 3 Jingzhao Zhang 1 2 1IIIS, Tsinghua University 2 College of AI, Tsinghua University 3 Apex Intelligence 4 School of Mathematical Sciences, Tongji University 5 Department of Artificial Intelligence, Westlake University chenlc23, xinliang23@@mails.tsinghua.edu.cn, wanghengyu@apexin.ai liuchengchang@westlake.edu.cn, cyc@apexin.ai, jingzhaoz@mail.tsinghua.edu.cn Abstract We study second- and higher-order methods for solving smooth monotone variational inequalities (MVI). Monteiro and Svaiter (SIAM J. Optim., 2012) showed that a second-order method, NPE, converges at the rate of (T−1.5)O(T^-1.5). For convex-concave minimax optimization, a subset of MVI problems, Chen, Liu, Luo, and Zhang (COLT 2025) recently improved the complexity to ~(T−1.75) O(T^-1.75) . However, it is open whether the conjectured complexity for MVI can be improved. In this paper, by using a large-step inexact Halpern iteration, we propose a novel Halpern-NPE method that achieves an even faster rate of ~(T−2) O(T^-2) for solving MVIs. We also provide the ppth-order generalization of our method. We first introduce an Anchored Tensor Method (ATM) that achieves the rate of (T−(p−1))O(T^-(p-1)), and then combine it with the Halpern iteration to achieve a faster convergence rate of ~(T−p) O(T^-p). This improves all prior results for p≥2p≥ 2 and matches the classical extragradient method for p=1p=1. AI Usage. Lesi Chen, Chengchang Liu, and Jingzhao Zhang shortlisted multiple open problems in optimization and submitted them to the auto-research platform of Apex Intelligence founded by Yongchao Chen. The team conducted extensive large-scale searches; Hengyu Wang first discovered a proof of the rate (T−(p−1))O(T^-(p-1)) with Claude Opus 4.6. Upon verifying the result, the authors conjectured a better result of ~(T−p) O(T^-p), which Xinliang Zhang later found a proof with GPT 5.6 Sol. These results are then verified by Lesi Chen, Xinliang Zhang, Chengchang Liu, and Jingzhao Zhang. The authors also express their gratitude to independent researcher Junwei Zhou (zjw330501@gmail.com), as well as Xiaoyu Cao (caoxiaoyu@apexin.ai) and Huan Wang (wanghuan@apexin.ai) from Apex for their invaluable support during the proof search process. NoHyper**footnotetext: Equal contributions. 1 Introduction Let ⊆ℝdX ^d be a nonempty compact convex set and let F:→ℝdF:X→R^d be a continuously differentiable monotone operator. We study the monotone variational inequality (MVI) problem [24], which targets at finding an solution ∗∈ x^*∈X such that ⟨F(∗),−∗⟩≥0∀∈. F( x^*), x- x^* ≥ 0 ∀ x∈X. (1) Let ()N_X( x) be the normal cone of X at point x, equivalently, the subdifferential of the indicator function of X. The MVI problem is equivalent to monotone inclusion problem ∈(∗)+(∗) 0∈ F( x^*)+N_X( x^*), which captures a lot of optimization problems, especially solving game-theoretical equilibria [37, 27, 56, 7, 36]. A typical example is the following convex-concave minimax optimization problem [28, 16]: min∈max∈ϕ(,), _ u∈U _ v∈Vφ( u, v), (2) which is an MVI problem for the gradient operator (,)=[∇ϕ(,),−∇ϕ(,)] F( u, v)=[ _ uφ( u, v),- _ vφ( u, v)]. The minimax optimization has wide applications in many machine learning problems such as adversarial training [65], AUC maximization [63], and distributionally robust optimization [15]. Following the oracle model established in classical textbooks [49, 54], we consider the ppth-order algorithm class that can query the derivative information of F up to order p−1p-1: ((),D(),…,Dp−1()).( F( x),D F( x),…,D^p-1 F( x)). (3) We also assume that the operator is ppth-order LpL_p-smooth, i.e., Dp−1D^p-1 F is LpL_p-Lipschitz continuous. The optimal oracle complexity for convex optimization has been settled, where the monotone operator ()=∇f() F( x)=∇ f( x) is the gradient of a convex function f. Let T be the total number of oracle calls. For first-order methods (p=1p=1), Nemirovskij and Yudin [49] showed a lower bound of Ω(L1/T2) (L_1/T^2), and Nesterov [51] proposed the optimal method that achieves the matching convergence rate of (L1/T2)O(L_1/T^2). For second- and higher-order methods (p≥2p≥ 2), Monteiro and Svaiter [47] proposed near-optimal methods that converge to ~(L2/T3.5) O(L_2/T^3.5), then Kovalev and Gasnikov [39] and Carmon et al. [14] independently improved it to (L2/T3.5)O(L_2/T^3.5) in p=2p=2 and (Lp/T(3p+1)/2)O(L_p/T^(3p+1)/2) for general p≥2p≥ 2. Importantly, Arjevani et al. [6] provided matching lower bounds of Ω(Lp/T(3p+1)/2) (L_p/T^(3p+1)/2) for all ppth-order methods. However, the complexity of MVI beyond convex optimization is not fully understood besides first-order methods. For p=1p=1, Korpelevich [38] proposed the extragradient method with a convergence rate of (L1/T)O(L_1/T), and Nemirovski [48] showed a matching lower bound of Ω(L1/T) (L_1/T) in 2004. In contrast, the optimal oracle complexity for p≥2p≥ 2 has remained open for many years. In 2012, Monteiro and Svaiter [46] proposed the Newton Proximal Extragradient (NPE) method with a fast second-order convergence rate of (L2/T1.5)O(L_2/T^1.5) for p=2p=2, where the cost at each iteration is nearly the same as matrix inversion/multiplication [23, 21]. Moreover, the ppth-order generalization of NPE achieves the rate of (Lp/T(p+1)/2)O(L_p/T^(p+1)/2) in general [9, 1, 31, 43, 35, 42, 55], which is speculated to be optimal in Adil et al. [1], Lin and Jordan [42]. Very recently, Chen et al. [17] proposed a faster method tailored to minimax Problem (2) that can achieve a new convergence rate of ~(L2/T1.75) O(L_2/T^1.75) using second-order oracles and ~(Lp/T(3p+1)/4) O(L_p/T^(3p+1)/4) using ppth-order oracles [19]. However, it is unknown whether this upper bound is optimal, as it persists a gap compared to the lower bound of Ω(L2/T2.5) (L_2/T^2.5) and Ω(Lp/T(3p−1)/2) (L_p/T^(3p-1)/2) for second-order and ppth-order oracles, respectively [19]. More importantly, it is also open whether the classical (Lp/T(p+1)/2)O(L_p/T^(p+1)/2) rate [46, 42] can be improved for general MVI, without leveraging the minimax problem structure. Contributions. In this paper, we propose a superfast high-order method that solves the MVI problem at the convergence rate of ~(L2/T2) O(L_2/T^2) for p=2p=2 and ~(Lp/Tp) O(L_p/T^p) for general p≥2p≥ 2. Our new upper bound achieves a significant improvement in the exponent compared to both the classical rate of (Lp/T(p+1)/2)O(L_p/T^(p+1)/2) for MVI and the recently established fast rate of ~(Lp/T(3p+1)/4) O(L_p/T^(3p+1)/4) for minimax problems, because p>3p+14>p+12,∀p≥2.p> 3p+14> p+12, ∀ p≥ 2. Notations. We use ∥⋅∥\|\,·\,\| to denote the Euclidean norm for vectors and the spectral norm for matrices and tensors in a unified way. We hide logarithmic factors in the notation ~(⋅) O(\,·\,). Also, we use the notations p(⋅)O_p(\,·\,) and ~p(⋅) O_p(\,·\,) to hide the costants that depend on p. We denote DqD^q F as the qqth-order derivative of an operator :→ℝd F:X→R^d. We also let Proj Proj_X be the projection operator of x onto the set X. 1.1 Technical Overview The algorithm and analysis is surprisingly simple. We achieve the new result based on the Halpern iteration [29] for finding the fixed point of a non-expansive operator :ℝd→ℝ P:R^d→R, which takes the form of t+1=βt0+(1−βt)(t), x_t+1= _t x_0+(1- _t) P( x_t), (4) where βtt=0T−1∈(0,1)\ _t\_t=0^T-1∈(0,1) is a decreasing sequence and the initial point 0 x_0 is also called the anchor [64, 40]. A simple and common choice of βt=(t−1) _t=O(t^-1) can lead to the optimal convergence rate to the fixed point ‖T−(T)‖≤(T−1)\| x_T- P( x_T)\|≤O(T^-1) [41, 29, 3, 10]. When solving MVI problems, a natural candidate operator P that satisfies non-expansiveness is the resolvent/proximal operator η(+)=(Id+η(+))−1) P_η( F+N_X)=( Id+η( F+N_X))^-1) [60, 62]. Let res()=‖−η(+)()‖/η res( x)=\| x- P_η( F+N_X)( x)\|/η be the normalized proximal residual then the Halpern iteration guarantees the convergence rate of res(T)=((ηT)−1) res( x_T)=O((η T)^-1). For first-order methods, we can take η=1/L1η=1/L_1, and the resolvent can be solved in (1)O(1) gradient steps. As a result, it can convert a method with the optimal (L1/T)O(L_1/T) convergence rate in the gap function [52, 38] to a counterpart with the same rate in terms of the stronger notions, known as the gradient norm [64, 40] or tangent residual [11, 13, 12]. Maybe surprisingly, in the high-order setting, we found that the power of Halpern iteration [22, 10] or anchoring technique [61, 64, 40] is not limited to converting solution concepts as in p=1p=1, but also can lead to a faster convergence rate for p≥2p≥ 2. We briefly introduce how to prove the new results as follows. For the case p=2. We show that a simple algorithm called Halpern-NPE can achieve the fast convergence rate of ~(T−2) O(T^-2) by using the following two observations: 1. Since the proximal subproblem is L2L_2-second-order smooth and η−1η^-1-strongly monotone, the NPE method [46, 31, 42, 1] achieve the complexity of Nt=~((ηL2Rt)2/3)N_t= O((η L_2R_t)^2/3) for solving the ttth subproblem, where Rt=‖t−η(+)(t)‖=ηres(t)R_t=\| x_t- P_η( F+N_X)( x_t)\|=η res( x_t) is the distance between the initialization t x_t and the proximal point η(+)(t) P_η( F+N_X)( x_t). 2. Recall that the guarantee of the Halpern iteration immediately give Rt=(t−1)R_t=O(t^-1). Therefore, we can select a very large stepsize of η=(T/L2)η=O(T/L_2), while ensuring the total costs remains nearly unchanged since ∑t=0T−1Nt=~(T) _t=0^T-1N_t= O(T). Finally, substituting the setting of η into res(T)=((ηT)−1) res( x_T)=O((η T)^-1) leads to the claimed convergence rate of ~(L2/T2) O(L_2/T^2). The anchor point 0 x_0 in the Halpern iteration is indispensable for our new results. Without the anchor point, i.e., setting βt=0 _t=0 in equation 4, the Halpern iteration reduces to the proximal point iteration [60]. Although the latter also converges at the rate of ((ηT)−1)O((η T)^-1) in the gap function [45], we can only prove ∑t=0T−1Rt2=(1) _t=0^T-1R_t^2=O(1), thus the total cost of subproblem solving is ηL2∑t=0T−1Rt=(ηL2T)η L_2 _t=0^T-1R_t=O(η L_2 T). Hence, one can only choose η=Θ(T/L2)η= ( T/L_2) and achieve the classical rate of (L2/T1.5)O(L_2/T^1.5) as prior works [46, 34, 5]. For the case p≥2p≥ 2. To generalize the above result to higher-order optimization, we can easily follow the same analysis and know that, as long as there exists a basic tensor method that has the complexity of (ηLpdt)1/(p−1)(η L_pd_t)^1/(p-1) for solving the proximal subproblem, then the large-step Halpern iteration with η=Θ(Tp−1/Lp)η= (T^p-1/L_p) yields a method converges at the rate of ~(T−p) O(T^-p). However, the high-order generalization of NPE [9, 1, 31, 42] can only achieve a complexity of ~(ηLpdt)2/(p+1)) O(η L_pd_t)^2/(p+1)) that is not fast enough for p≥4p≥ 4. To address this issue, we propose an Anchored Tensor Method (ATM), which iteratively adds regularization on the original operator to ensure that every tensor step lies in the local superlinear convergence region. We show that ATM achieves a rate of (T−(p−1))O(T^-(p-1)) for p≥2p≥ 2, which implies the required complexity of (ηLpdt)1/(p−1)(η L_pd_t)^1/(p-1) for the proximal subproblem in the Halpern iteration. Finally, the algorithm called Halpern-ATM achieves the fast convergence rate of ~(T−p) O(T^-p). 1.2 Related Works Convex optimization. When the operator ()=∇f() F( x)=∇ f( x) is the gradient of a convex function f, Nesterov and Polyak [50] proposed the cubic regularized Newton (CRN) method, which is the first globally convergent second-order method and achieves a rate of (L2/T2)O(L_2/T^2) for p=2p=2. Nesterov [53] proposed the accelerated CRN to achieve a fast rate of (L2/T3)O(L_2/T^3). Monteiro and Svaiter [47] proposed the accelerated Newton proximal extragradient (A-NPE) method that converges at an even faster rate of ~(L2/T3.5) O(L_2/T^3.5). For p≥2p≥ 2, Gasnikov et al. [26], Bubeck et al. [8], Jiang et al. [32] proposed ppth-order generalization of A-NPE that converges at the rate of ~(Lp/T(3p+1)/2) O(L_p/T^(3p+1)/2). Very recently, Kovalev and Gasnikov [39], Carmon et al. [14] removed the bisection sub-procedure in A-NPE and achieved the optimal rate of (Lp/T(3p+1)/2)O(L_p/T^(3p+1)/2). On the lower bound part, Agarwal and Hazan [2] showed a lower bound of Ω(Lp/T(5p+1)/2) (L_p/T^(5p+1)/2) for randomized algorithms. Concurrently, Arjevani et al. [6] showed the optimal lower bound of Ω(Lp/T(3p+1)/2) (L_p/T^(3p+1)/2) for deterministic algorithms. Recently, Garg et al. [25] improved the lower bound of randomized and quantum algorithms to Ω~(Lp/T(3p+1)/2) (L_p/T^(3p+1)/2), which has no gap between the upper bounds up to logarithmic factors. Monotone variational inequalities. For a general operator F, Monteiro and Svaiter [46] proposed the Newton proximal extragradient (NPE) method that globally converges at the rate of (L2/T1.5)O(L_2/T^1.5) for p=2p=2, using T second-order oracle calls and (TlogT)O(T T) matrix inversion operations. Bullins and Lai [9] generalized NPE to ppth order and showed a convergence rate of (Lp/T(p+1)/2)O(L_p/T^(p+1)/2). Subsequently, many simpler analyses or alternative algorithms have been found [31, 1, 43, 35, 55, 42], but all the shown convergence rates are (Lp/T(p+1)/2)O(L_p/T^(p+1)/2). Moreover, Jiang et al. [34], Alves and Svaiter [5] proposed bisection-free methods for p=2p=2 that also converge at the rate of (L2/T1.5)O(L_2/T^1.5) and only require a single matrix inversion at each iteration. When F is the gradient operator of a convex-concave function ϕ:→ℝφ:X→R, Chen et al. [17] applied a primal-dual Monteiro-Svaiter acceleration [47] to the proximal function and achieved a fast rate of (L2/T1.75)O(L_2/T^1.75) for second-order minimax optimization (p=2p=2). In the full version, Chen et al. [19] showed the ppth-order generalization achieves the rate of (Lp/T(3p+1)/4)O(L_p/T^(3p+1)/4) and also established a lower bound of Ω(Lp/T(3p−1)/2) (L_p/T^(3p-1)/2). This paper reduces the gap by proposing better upper bounds of (Lp/Tp)O(L_p/T^p). 2 Preliminaries 2.1 Main Assumptions We study the MVI Problem (1) under the following standard assumptions [31, 42, 9, 35]. Assumption 2.1. ⊆ℝdX ^d is a nonempty compact convex set. Assumption 2.2. We assume there exists ⋆, x , such that 0∈(+)(⋆)0∈( F+N_X)( x ). Assumption 2.3. :→ℝd F:X ^d is continuous and monotone: ⟨()−(),−⟩≥0,∀,∈. F( x)- F( y), x- y ≥ 0, ∀ x, y∈X. Assumption 2.4. Assume F:→ℝdF:X→R^d is ppth-order LpL_p smooth: ‖Dp−1()−Dp−1()‖≤Lp‖−‖,∀,∈. \|D^p-1 F( x)-D^p-1 F( y) \|≤ L_p \| x- y \|, ∀ x, y∈X. (5) 2.2 Tensor Steps A basic operation to leverage the ppth-order oracle in equation 3 is the following tensor step [31, 42, 9, 35], which solves the MVI problem/monotone inclusion induced by a local Taylor approximation. Definition 2.1. Under Assumption 2.3 and 2.4, for an input point ∈ x∈X the ppth-order tensor step with regularization parameter M≥LpM≥ L_p outputs =p(;M) y=T_ F^p( x;M) such that ∈¯p()+Mp!‖−‖p−1(−)+(), 0∈ F_ x^p( y)+ Mp! \| y- x \|^p-1( y- x)+N_X( y), (6) where ¯p():=∑k=0p−1Dk()[−]k/k! F_ x^p( y):= _k=0^p-1D^k F( x)[ y- x]^k/k! is (p−1)(p-1)th-order Taylor expansion for F at the center point x. When p=1p=1, the above tensor step is exactly the (projected) gradient step that can be conducted in vector addition operators; When p=2p=2, it can be solved in the same spirit of cubic regularized Newton subproblem [50] using a similar binary search, whose costs is nearly the same as matrix multiplication/inversion time [23, 21]; In general (p≥2p≥ 2), the MVI problem in equation 6 can be solved in polynomial time using the interior point method [59, 58] or the cutting plane method [33]. 2.3 Resolvent Operators and Residuals To apply the Halpern iteration to MVI problems, we make use of the following resolvent/proximal operator [60, 20, 24]. Definition 2.2. Under Assumption 2.3, for η>0η>0, we can define the unique-valued resolvent operator as η()=(Id+η(+))−1(). P_η( x)=( Id+η( F+N_X))^-1( x). It is well-known that η P_η is non-expansive if F is monotone [62, 60]. Proposition 2.1 (Ryu and Boyd [62, Section 6]). Under Assumption 2.3, the operator η P_η is non-expansive for any η>0η>0, ‖η()−η(′)‖≤‖−′‖,∀,′∈.\| P_η( x)- P_η( x )\|≤\| x- x \|, ∀ x, x ∈X. Moreover, it is easy to see that ∗ x^* solved the MVI problem, i.e., ∈(∗)+(∗) 0∈ F( x^*)+N_X( x^*), if and only if ∗ x^* is the fixed point of η P_η, i.e., ∗=η(∗) x^*= P_η( x^*). It motives the definition of the following proximal residual [22, 10, 3]. Definition 2.3. Under Assumption 2.3, we define the proximal residual for a point ∈ x∈X as res():=‖−η()‖/η. res( x):=\| x- P_η( x)\|/η. In addition, finding a point with a small proximal residual is sufficient for approximating the MVI problem: Proposition 2.2 (Cai et al. [10, Proposition 4.1]). Under Assumption 2.3, then for =η() y= P_η( x), dist(,()+())≤res(). dist( 0, F( y)+N_X( y))≤ res( x). (7) The left-hand side in equation 7 is also called the tangent residual [11, 13, 12] at the point y, which implies a strong solution/Stampacchia variational inequality solution [30]. A strong solution also implies a weak solution/Minty variational inequality solution [44] measured by the restricted gap function [52] gap():=sup′∈⟨(),−′⟩ gap( y):= _ y ∈X F( y), y- y on a compact set, because we have gap()≤dist(,()+())⋅diam(D) gap( y)≤ dist( 0, F( y)+N_X( y))· diam(D) by the monotonicity of F and Cauchy–Schwarz. Finally, although the right-hand side in equation 7 depends on the proximal point =η() y= P_η( x), it is straightforward to produce ^∈ y∈X such that dist(,(^)+(^))≤2res() dist( 0, F( y)+N_X( y))≤ 2 res( x) by approximately solving the proximal subproblem. See, e.g. Cai et al. [10, Lemma C.4] for the case p=1p=1, and for general p≥2p≥ 2 is similar. 3 Main Result The main theorem below shows a new complexity upper bound that of T=~(ϵ−1/p)T= O(ε^-1/p) for finding an ϵε-solution such that res()≤ϵ res( x)≤ε, which is equivalent to the convergence rate res()=~(T−p) res( x)= O(T^-p) as claimed. Theorem 3.1. Under Assumptions 2.1-2.4, for any integer p≥2p≥ 2, there exists an algorithm (Algorithm 2) that can return an ϵε-solution ∈ x∈X satisfying res()≤ϵ res( x)≤ε in the ppth-order oracle complexity of T=~(D(Lpϵ)1/p),T= O (D ( L_pε )^1/p ), where D=‖0−∗‖D= \| x_0- x^* \| is the distance of the initial point 0∈ x_0∈X to the optimal solution ∗ x^*. For p=2p=2, we obtain a Newton method for MVI problems that has the complexity of ~(ϵ−1/2) O(ε^-1/2), improving the classical result of (ϵ−2/3)O(ε^-2/3) for MVI [46] and the recent result of (ϵ−4/7)O(ε^-4/7) for minimax problems [17]. Compared to the Ω(ϵ−2/5) (ε^-2/5) lower bound proved in [19], a gap of ϵ−1/10ε^-1/10 persists. Therefore, the oracle complexity remains open even for p=2p=2. For the more general setting p≥2p≥ 2, our new upper bounds also improves the (ϵ−2/(p+1))O(ε^-2/(p+1)) for MVI [46] and the (ϵ−4/(3p+1))O(ε^-4/(3p+1)) for minimax problems [19]. Compared to the lower bound of Ω(ϵ−2/(3p−1)) (ε^-2/(3p-1)), the gap becomes larger as p grows. 4 Halpern-NPE Achieves the ~(T−2) O(T^-2) Rate for p=2p=2 In this section, we first prove Theorem 3.1 for p=2p=2 by introducing a simple second-order method that achieves the fast convergence rate of ~(T−2) O(T^-2). Following the ideas in Section 1.1, we apply a large-step inexact Halpern iteration on the resolvent operator, which is solved by a restarted version [31, 42] of the NPE method [46]. We formally introduce our double-loop Halpern-NPE method algorithm in the following. 4.1 Outer Loop: Inexact Halpern Iteration In equation 4, a simple choice of the anchoring coefficient βt _t for optimal convergence rate is βt=1/(t+2) _t=1/(t+2) [41, 22]. Then, we obtain the following scheme by applying the Halpern iteration βt=1/(t+2) _t=1/(t+2) on the approximate resolvent t≈η(t) y_t≈ P_η( x_t): t+1=1t+20+t+1t+2t,wheret≈η(t). x_t+1= 1t+2 x_0+ t+1t+2 y_t, where y_t≈ P_η( x_t). (8) This inexact Halpern iteration has been analyzed in [22, 3, 10]. This iteration converges at the rate of res(T)=((ηT)−1) res( x_T)=O((η T)^-1) if the approximation error of the resolvent is sufficiently small at every step Lemma 4.1 (Alacaoglu et al. [3, Theorem 2.1]). Under Assumptions 2.1-2.3, if the resolvent approximation error satisifes ‖t−η(t)‖≤δt:=‖t−η(t)‖98t+2log(t+2),t=0,⋯,T−1,\| y_t- P_η( x_t)\|≤ _t:= \| x_t- P_η( x_t)\|98 t+2 (t+2), t=0,·s,T-1, (9) then the inexact Halpern iteration in equation 8 guarantees a convergence rate of res(T)≤4‖0−∗‖η(T+1). res( x_T)≤ 4\| x_0- x^*\|η(T+1). (10) The above result has been extensively used in first-order methods, including parameter-free methods [22], finite-sum problems [10], and structured non-monotone problems [3]. In the high-order settings considered in this paper, we show that it allows a larger stepsize of η=(T/L2)η=O(T/L_2) than the choice of η=(T/L2)η=O( T/L_2) in the proximal point method [46, 9, 42, 31]. 4.2 Inner Loop: Approximating the Resolvent Recall that the resolvent in Definition 2.2 is η()=(Id+η(+))−1() P_η( x)=( Id+η( F+N_X))^-1( x). Equivalently, the proximal point =η() y= P_η( x) is the unique solution of ∈()+1η(−)+(). 0∈ F( y)+ 1η( y- x)+N_X( y). (11) For a fixed x, the above problem is equivalent to solving the MVI problem induced by the regularized operator ():=()+(−)/η G( y):= F( y)+( y- x)/η, which is μ-strongly monotone for μ=1/ημ=1/η: ⟨()−(′),−′⟩≥μ‖−′‖2,∀,′∈. G( y)- G( y ), y- y ≥μ\| y- y \|^2, ∀ y, y ∈X. (12) The complexity of obtaining an approximate solution follows from existing results for smooth strongly MVIs [46, 35, 31, 57]. For instance, we can use the restarted NPE/ARE method [46, 31], which is essentially a second-order extragradient method [38] with large adaptive stepsize τk=(1/(L2‖k+1/2−k‖)) _k=O(1/(L_2\| x_k+1/2- x_k\|)) [46, 31]. See Algorithm 1 for the complete procedure, whose guarantee is stated as follows. Lemma 4.2 (Huang and Zhang [31, Theorem 3.2]). Let the operator :→ℝd G:X→R^d be L2L_2-second-order smooth satisfying equation 5 and η−1η^-1-strongly monotone. Apply Algorithm 1 with the parameters M=Θ(L2),K=((ηL2R)2/3),S=(loglog(R/δ))M= (L_2), K=O ((η L_2R)^2/3 ), S=O ( (R/δ) ) can return a point S∈ y_S∈X such that ‖S−∗‖≤δ\| y_S- y^*\|≤δ in the total ppth-order oracle complexity of K×SK× S, where ∗=(+)−1() y^*=( G+N_X)^-1( 0) is the unique solution to the MVI induced by G and R=‖0−∗‖R=\| y_0- y^*\|. Algorithm 1 Restarted-NPE(,0,M,K,S)( G, y_0,M,K,S) 1:for s=0,⋯,S−1s=0,·s,S-1 2: s,0=s y_s,0= y_s 3: for k=0,⋯,K−1k=0,·s,K-1 4: Perform a cubic regularized Newton step s,k+1/2=2(s,k;M) y_s,k+1/2=T_ G^2( y_s,k;M) 5: Set the adaptive stepsize τs,k=1/(M‖s,k+1/2−s,k‖) _s,k=1/(M\| y_s,k+1/2- y_s,k\|)) 6: Perform an extragradient step s,k+1=1(s,k;τs,k−1) y_s,k+1=T_ G^1( y_s,k; _s,k^-1) 7: end for 8: Restart at s+1=∑s=0S−1τs,ks,k+1/2/∑s=0S−1τs,k y_s+1= _s=0^S-1 _s,k~ y_s,k+1/2~ / _s=0^S-1 _s,k 9:end for 10:return S y_S We remark that the above complexity of K×SK× S can be reduced to K+SK+S with a refined analysis, as in Jiang and Mokhtari [35, Section 7.2] or Huang and Zhang [31, Theorem 4.1]. However, the slightly looser bound of K×SK× S is enough for our subsequent analyses. 4.3 Final Algorithm and Total Complexity Algorithm 2 Halpern-NPE(,0,η,T,M,Ktt=0T−1,Stt=0T−1)( F, x_0,η,T,M,\K_t\_t=0^T-1,\S_t\_t=0^T-1) 1:for t=0,⋯,T−1t=0,·s,T-1 2: Let t()=()+(−t)/η G_t( y)= F( y)+( y- x_t)/η and approximately solved ∈t()+() 0∈ G_t( y)+N_X( y) by t=Restarted-NPE(t,t,M,Kt,St) y_t= Restarted-NPE( G_t, x_t,M,K_t,S_t) 3: Perform the Halpern iteration t+1=1t+20+t+1t+2t x_t+1= 1t+2 x_0+ t+1t+2 y_t 4:end for 5:return T x_T By applying Algorithm 1 to solve the resolvent operator in equation 8, we obtain our final high-order method to solve MVI problems in Algorithm 2. Then, by appropriately selecting the stepsize η=Θ(T)η= (T) as well as sub-solver parameter sequences KtK_t and StS_t, we can obtain the total complexity of Algorithm 2 as follows. Theorem 4.1 (Theorem 3.1 for p=2p=2). Under Assumptions 2.1-2.4 for p=2p=2, apply Algorithm 2 with parameters T= T= (D(L2ϵ)1/2),η=Θ(TL2D),M=Θ(L2), O (D ( L_2ε )^1/2 ), η= ( TL_2D ), M= (L_2), Kt= K_t= ((Tt+1)2/3),St=Θ(loglog(t+2)), O ( ( Tt+1 )^2/3 ), S_t= ( (t+2)), then the algorithm can ensure res(T)≤ϵ res( x_T)≤ε in the total second-order oracle complexity of (TloglogT)=~(D(L2/ϵ)1/2),O(T T)= O(D (L_2/ε )^1/2), (13) where D=‖0−∗‖D= \| x_0- x^* \| is the distance of the initial point 0∈ x_0∈X to the optimal solution ∗ x^*. Proof. Let us prove by induction that our parameter setting ensures the convergence rate of res(T):=‖T−η(T)‖η≤4Dη(T+1). res( x_T):= \| x_T- P_η( x_T)\|η≤ 4Dη(T+1). (14) For t=0t=0, we know from the non-expansiveness of η P_η and the fixed-point property ∗=η(∗) x^*= P_η( x^*) that ‖0−η(0)‖≤‖0−η(∗)‖+‖η(∗)−η(0)‖≤2‖0−∗‖=2D,\| x_0- P_η( x_0)\|≤\| x_0- P_η( x^*)\|+\| P_η( x^*)- P_η( x_0)\|≤ 2\| x_0- x^*\|=2D, (15) which established the induction base of equation 14 by dividing η on both sides. Now, we assume equation 14 holds for all t≤T−1t≤ T-1. Then, Rt:=‖t−η(t)‖R_t:=\| x_t- P_η( x_t)\|, the initial distance of Algorithm 1 satisfies that Rt≤4Dt+1,∀t=0,⋯,T−1.R_t≤ 4Dt+1, ∀ t=0,·s,T-1. (16) Therefore, by Lemma 4.2, at each step the resolvent can be solved up to the approximation error δt _t in equation 9 under the parameter setting M=Θ(L2)M= (L_2), Kt=Θ((ηL2Rt)2/3)=Θ((TRtD)2/3)=((Tt+1)2/3),K_t= ((η L_2R_t)^2/3 )= ( ( TR_tD )^2/3 )=O ( ( Tt+1 )^2/3 ), and, by the definition of δt=Θ(Rt/(t+2log(t+2)) _t= (R_t/( t+2 (t+2)) in equation 9, St=Θ(loglog(Rt/δt))=Θ(loglog(t+2)).S_t= ( (R_t/ _t) )= ( (t+2) ). Therefore, the inexact condition in Lemma 4.1 is satisfied for t≤T−1t≤ T-1, which implies the convergence rate of equation 14 for t=Tt=T by the lemma. This completes the induction. Now, substituting our choices of η in equation 14 yields res(T)=(L2D2/T2), res( x_T)=O (L_2D^2/T^2 ), which is smaller than ϵε by our setting of T. Finally, the total number of oracle calls can be bounded by ∑t=0T−1KtSt=(∑t=0T−1(Tt+1)2/3loglogT)=(TloglogT), _t=0^T-1K_tS_t=O ( _t=0^T-1 ( Tt+1 )^2/3 T )=O(T T), (17) which is equivalent to the bound of ~(D(L2/ϵ)1/2) O(D (L_2/ε )^1/2) in equation 13 by substituting the choice of T. ∎ This theorem shows that the proposed Halpern-NPE converges at a fast rate of ~(T−2) O(T^-2), which improves both the classical rate of (T−1.5)O(T^-1.5) by NPE [46] and the rate of ~(T−1.75) O(T^-1.75) via primal-dual A-NPE [17] for minimax problems. In the next section, we generalize this theorem to all p≥2p≥ 2 and achieve the upper bound of ~(T−p) O(T^-p) claimed by the main Theorem 3.1. 5 Generalizing the Result to All p≥2p≥ 2 To achieve the rate of ~(T−p) O(T^-p) for all p≥2p≥ 2, a natural idea is to also use the Halpern iteration with a large stepsize η=Θ(Tp−1)η= (T^p-1) and reuse the analysis in the proof of Theorem 4.1. To ensure that each subproblem is solvable in ~(1) O(1) costs like in equation 24, the sub-solver for the resolvent operator should converges at least ~(exp(−μ1/(p−1)T)) O( (-μ^1/(p-1)~T)) for μ-strongly monotone operators, or equivalently (under black-box reductions [4, 57, 31, 42]), the ~(T−(p−1)) O(T^-(p-1)) for monotone operators. However, the classical convergence rate of (T−(p+1)/2)O(T^-(p+1)/2) achieved by high-order NPE [9, 31] does not satisfy the above requirement when p≥4p≥ 4. In the following, we first introduce an Anchored Tensor Method (ATM) that achieves the required convergence rate of (T−p)O(T^-p) in Section 5.1, then we further accelerate it with the Halpern iteration as in p=2p=2 to achieve the ~(T−p) O(T^-p) convergence rate in Section 5.2. 5.1 ATM Achieves the (T−(p−1))O(T^-(p-1)) Rate Recall Section 4.2 that approximating the resolvent is η()=(Id+η(+))−1() P_η( x)=( Id+η( F+N_X))^-1( x) is equivalent to solving the MVI problem induced by the η−1η^-1-strongly monotone operator ():=()+(−)/η G( y):= F( y)+( y- x)/η. We introduce a novel Anchored Tensor Method (ATM) that can solve the problem in the ppth-order oracle complexity of (η1/(p−1))O(η^1/(p-1)). It also implies a complexity of (ϵ−1/(p−1))O(ε^-1/(p-1)) for finding an ϵε-solution in the monotone case by applying the same algorithm on an ϵε-regularized operator [19, Lemma 4.1]. Algorithm 3 Anchored-Tensor-Method(,0,μ,M,R,K,S)( G, y_0,μ,M,R,K,S) 1:for k=0,⋯,K−1k=0,·s,K-1 2: Select regularization μk _k according to equation 19. 3: Define the regularized operator k()=()+(μk−μ)(−0) G_k( y)= G( y)+( _k-μ)( y- y_0) 4: Perform a tensor step k+1=kp(k;M) y_k+1=T_ G_k^p( y_k;M) 5:end for 6:for s=K+1,⋯,K+S−1s=K+1,·s,K+S-1 7: Perform a tensor step k+1=p(k;M) y_k+1=T_ G^p( y_k;M) 8:end for 9:return K+S y_K+S The starting point of the anchored tensor method is the following local convergence guarantee for the tensor step [31, 42]. Lemma 5.1 (Lin and Jordan [42, Theorem 3.5]). Let the operator :→ℝd G:X→R^d be LpL_p-pthpth-order smooth satisfying equation 5 and μ-strongly monotone satisfying equation 12. There exists a numerical constant C such that The tensor step +=p(;CLp) y^+=T_ G^p( y;CL_p) satisfies ‖+−∗‖≤CpηLp‖−∗‖(p+1)/2,Cp:=2p(5p−2)p!,\| y^+- y^*\|≤ C_pη L_p\| y- y^*\|^(p+1)/2, C_p:= 2^p(5p-2)p!, where ∗=(+)−1() y^*=( G+N_X)^-1( 0) is the unique solution to the MVI induced by G. Define the local region radius ρ(μ)ρ(μ) and contraction factor θp _p indicated by this lemma as ρ(μ):=12(μCpLp)1/(p−1),θp:=2−(p−1)/2<1.ρ(μ):= 12 ( μC_pL_p )^1/(p-1), _p:=2^-(p-1)/2<1. Consequently, Lemma 5.1 implies hat ‖−∗‖≤ρ(μ)⟹‖+−∗‖≤θpρ(μ).\| y- y^*\|≤ρ(μ) \| y^+- y^*\|≤ _pρ(μ). (18) It means that whenever y enters the local region, the tensor update converges superlinearly to the minimizer. Based on this observation, we formally propose a two-phase Anchored Tensor Method in Algorithm 3: 1. In the first phase with K iterations, the algorithm iteratively performs the tensor step on the anchored/regularized operator k():=()+(μk−μ)(−0) G_k( y):= G( y)+( _k-μ)( y- y_0). The anchoring coefficient μk _k is chosen such that every tensor step lies in the local superlinear convergence region and μk _k decreases over iterations as the algorithm approaches the minimizer. 2. At the beginning of the second phase, the anchoring coefficient μK _K has decreased to μK=μ _K=μ and the algorithm has entered the local region for the unregularized operator G. Therefore, the second phase reaches an δ-solution in additional S=(loglogδ)S=O( δ) iterations. To analyze the ATM method, we first give two lemmas for the anchored operator. Lemma 5.2 (Chen and Luo [18, Lemma A.1]). Let the operator G with MVI solution ∗ y^* satisfies Assumptions 2.1-2.3. For ν>0ν>0, we define ν()=()+ν(−0) G_ν( y)= G( y)+ν( y- y_0) be the regularized operator, and let ν∗ y_ν^* be the unique solution of 0∈ν()+()0∈ F_ν( y)+N_X( y). Then we have ‖ν∗−0‖2+‖ν∗−∗‖2≤‖0−∗‖2 \| y_ν^*- y_0 \|^2+ \| y_ν^*- y^* \|^2≤ \| y_0- y^* \|^2 Lemma 5.3 (Stability of anchoring). Under Assumptions 2.1-2.3, for every ν1>ν2>0 _1> _2>0, ‖ν1∗−ν2∗‖≤(1−ν2ν1)‖0−∗‖.\| y_ _1^*- y_ _2^*\|≤ (1- _2 _1 )\| y_0- y^*\|. Lemma 5.3 follows from the first-order optimality condition of ν∗ y_ν^* and the monotonicity of G. The formal proof is contained in Appendix A. With the help of these two lemmas, we carefully design the schedule of μk _k such that the local region rkr_k contracts fast. Formally, let CpC_p and θp _p be the constants implied by Lemma 5.1, we define rk:=maxρ(μ),(1R+k(1−θp)2R(p−1))−1,μk:=ρ−1(rk)=CpLp(2rk)p−1,r_k:= \ρ(μ), ( 1R+ k(1- _p)2R(p-1) )^-1 \, _k:=ρ^-1(r_k)=C_pL_p(2r_k)^p-1, (19) where R>0R>0 satisfies ‖0−∗‖≤R\| y_0- y^*\|≤ R. By this definition, rkr_k is decreasing from r0r_0 that satisfies r0≤Rr_0≤ R to ρ(μ)ρ(μ). Consequently, μj _j is decreasing from μ0 _0 such that μ0≥ρ−1(R) _0≥ρ^-1(R) to μ. We summarize the properties of these sequences in the following. Lemma 5.4. The sequences rkk=0K−1\r_k\_k=0^K-1 and μkk=0K−1\ _k\_k=0^K-1 defined in equation 19 satisfy the following: 1. rk+1/rk≥(1+θp)/2r_k+1/r_k≥(1+ _p)/2. 2. 1−μk+1/μk≤(1−θp)rk/(2R)1- _k+1/ _k≤(1- _p)r_k/(2R). Lemma 5.4 follows from basic algebraic calculations, and the full proof is contained in Appendix B. Now, we are ready to give the complexity bound of the ATM (Algorithm 3). Theorem 5.1. Let the operator :→ℝd G:X→R^d be LpL_p-pthpth-order smooth satisfying equation 5 and μ-strongly monotone satisfying equation 12. Apply Algorithm 3 with the parameters M=Θ(Lp),K=(R(Lp/μ)1/(p−1)),S=(loglog(D/δ))M= (L_p), K=O (R(L_p/μ)^1/(p-1) ), S=O( (D/δ)) can return a point K+S∈ y_K+S∈X such that ‖K+S−∗‖≤δ\| y_K+S- y^*\|≤δ in the total ppth-order oracle complexity of K+SK+S, where ∗=(+)−1() y^*=( G+N_X)^-1( 0) is the unique solution to the MVI induced by G and R=‖0−∗‖R=\| y_0- y^*\|. Proof. Let us prove by induction that our parameter setting ensures that ‖k−k∗‖≤rk=ρ(μk),\| y_k- y_k^*\|≤ r_k=ρ( _k), (20) where k∗ y_k^* is the unique solution of the MVI induced by the anchored operator k G_k that satisfies 0∈k()+()0∈ G_k( y)+N_X( y). For k=0k=0, we know from Lemma 5.2 that ‖0−0∗‖≤R\| y_0- y_0^*\|≤ R, which established the induction base since r0=Rr_0=R by definition. Now, we assume equation 20 holds for the kkth iteration. Then, by equation 18, ‖k+1−k∗‖≤θprk.\| y_k+1- y_k^*\|≤ _pr_k. Moreover, by Lemma 5.3 and Lemma 5.4 (part 2), we have ‖k+1∗−k∗‖≤(1−μk+1μk)R≤(1−θp)rk2.\| y_k+1^*- y_k^*\|≤ (1- _k+1 _k )R≤ (1- _p)r_k2. Combining the above two inequalities and using the triangle inequality, we obtain ‖k+1−k+1∗‖≤‖k+1−k∗‖+‖k+1∗−k∗‖≤(1+θp)rk2≤rk+1,\| y_k+1- y_k+1^*\|≤\| y_k+1- y_k^*\|+\| y_k+1^*- y_k^*\|≤ (1+ _p)r_k2≤ r_k+1, where the last step uses Lemma 5.4 (part 1). This completes the induction. Consequently, all the tensor step in phase 1 lies in the local superlinear convergence region in Lemma 5.1. We let the total number of iterations of phase 1 be the smallest number such that μk _k decays to μ. By equation 19, such a K satisfies K=(R(Lp/μ)1/(p−1)).K=O (R(L_p/μ)^1/(p-1) ). After phase 1 ends, we have μK=μ _K=μ and thus K= G_K= G. Since the initialization of phase 2 already lies in the local region for G, by equation 18, phase 2 can find a δ-solution in S=(loglog(D/δ)S=O( (D/δ) iterations. ∎ Treating the (loglog(D/δ))O( (D/δ)) factor as a constant, the above theorem shows that ATM achieves a new complexity upper bound of ((Lp/μ)1/(p−1))O((L_p/μ)^1/(p-1)) for LpL_p-ppth-order smooth and μ-strongly monotone operators. This improves the classical ((Lp/μ)2/(p+1))O((L_p/μ)^2/(p+1)) complexity achieved by high-order NPE [46, 9, 42, 31, 1] for any p≥4p≥ 4. As we have discussed, it also implies a complexity of (ϵ−1/(p−1))O(ε^-1/(p-1)) for finding an ϵε-solution in the monotone case by applying the same algorithm on an ϵε-regularized operator [19, Lemma 4.1]. 5.2 Halpern-ATM Achieves the ~(T−p) O(T^-p) Rate Algorithm 4 Halpern-ATM(,0,η,T,M,D,Ktt=0T−1,Stt=0T−1)( F, x_0,η,T,M,D,\K_t\_t=0^T-1,\S_t\_t=0^T-1) 1:for t=0,⋯,T−1t=0,·s,T-1 2: Let t()=()+(−t)/η G_t( y)= F( y)+( y- x_t)/η and approximately solved ∈t()+() 0∈ G_t( y)+N_X( y) by t=Anchored-Tensor-Method(t,t,1/η,M,Rt,Kt,St), y_t= Anchored-Tensor-Method( G_t, x_t,1/η,M,R_t,K_t,S_t), where Rt=4D/(t+1)R_t=4D/(t+1) according to equation 16. 3: Perform the Halpern iteration t+1=1t+20+t+1t+2t x_t+1= 1t+2 x_0+ t+1t+2 y_t 4:end for 5:return T x_T The ATM introduced in the previous section achieves the complexity upper bound of ~(ϵ−1/(p−1)) O(ε^-1/(p-1)) for monotone problems. In this section, we further show an accelerated complexity of ~(ϵ−1/p) O(ε^-1/p) by applying ATM to solve the resolvent operator in the Halpern iteration in Section 4.1. We present the resulting method in Algorithm 4. In light of Theorem 4.1, we choose a large stepsize η=Θ(Tp−1)η= (T^p-1), which indicates that the outer Halpern iteration converges at the rate of ((ηT)−1)=(T−p)O((η T)^-1)=O(T^-p). We can then complete the proof by providing a global upper bound of the total oracle complexity via a similar analysis to the proof of Theorem 4.1 in p=2p=2. This leads to the following main theorem. Theorem 5.2 (Full version of Theorem 3.1 ). Under Assumptions 2.1-2.4, for any integer p≥2p≥ 2, apply Algorithm 4 with parameters T= T= (D(Lpϵ)1/p),η=Θ(Tp−1LpDp−1),M=Θ(Lp), O (D ( L_pε )^1/p ), η= ( T^p-1L_pD^p-1 ), M= (L_p), Kt= K_t= (Tt+1),St=Θ(loglog(t+2)), O ( Tt+1 ), S_t= ( (t+2)), then the algorithm can ensure res(T)≤ϵ res( x_T)≤ε in the total second-order oracle complexity of (TlogT)=~(D(Lp/ϵ)1/p),O(T T)= O(D (L_p/ε )^1/p), (21) where D=‖0−∗‖D= \| x_0- x^* \| is the distance of the initial point 0∈ x_0∈X to the optimal solution ∗ x^*. The proof can be done by replacing the sub-solver in p=2p=2 with ATM for general p≥2p≥ 2 and following the analysis in the proof of Theorem 4.1. The formal proof is contained in Appendix C. This theorem shows that the proposed Halpern-ATM converges at a fast rate of ~(T−p) O(T^-p), which improves both the classical rate of (T−2/(p+1))O(T^-2/(p+1)) by high-order NPE [9, 1, 42, 31] and the rate of ~(T−4/(3p+1)) O(T^-4/(3p+1)) via high-order primal-dual A-NPE [19] for minimax problems. 6 Conclusion and Future Work In this paper, we propose a novel Halpern-NPE method that achieves a fast convergence rate of ~(T−2) O(T^-2) for solving MVI problems, which improves the classical rate of (T−1.5)O(T^-1.5) by NPE [46] and the rate of ~(T−1.75) O(T^-1.75) via primal-dual A-NPE [17] for minimax problems. We also provide the ppth-order generalization of our method, which can achieve the fast rate of ~(T−p) O(T^-p). However, the optimality of our methods remains open for p≥2p≥ 2, where the lower bound in [19] is Ω(T−2.5) (T^-2.5) for p=2p=2 and Ω(T−(3p−1)/2) (T^-(3p-1)/2) in general. We hope the gap can be closed in future work. Acknowledgment As noted in the abstract, the rates are first proved with the assistance of AI. The initial resolvent sub-solver given by AI has triple loops, and the author simplified it with human-readable analyses to Algorithm 1 for p=2p=2 and the single-loop method Algorithm 3 for p≥2p≥ 2. See Appendix D for the full algorithm and analyses provided by AI. References [1] D. Adil, B. Bullins, A. Jambulapati, and S. Sachdeva (2022) Optimal methods for higher-order smooth monotone variational inequalities. arXiv preprint arXiv:2205.06167. Cited by: item 1, §1.1, §1.2, §1, §5.1, §5.2. [2] N. Agarwal and E. Hazan (2018) Lower bounds for higher-order convex optimization. In COLT, Cited by: §1.2. [3] A. Alacaoglu, D. Kim, and S. J. Wright (2024) Revisiting inexact fixed-point iterations for min-max problems: stochasticity and structured nonconvexity. In ICML, Cited by: §1.1, §2.3, §4.1, §4.1, Lemma 4.1. [4] Z. Allen-Zhu and E. Hazan (2016) Optimal black-box reductions between optimization objectives. In NeurIPS, Cited by: §5. [5] M. M. Alves and B. F. Svaiter (2024) A search-free (1/k3/2)O(1/k^3/2) homotopy inexact proximal-newton extragradient algorithm for monotone variational inequalities. SIAM Journal on Optimization 34 (4), p. 3235–3258. Cited by: §1.1, §1.2. [6] Y. Arjevani, O. Shamir, and R. Shiff (2019) Oracle complexity of second-order methods for smooth convex optimization. Mathematical Programming 178 (1), p. 327–360. Cited by: §1.2, §1. [7] W. Ba, T. Lin, J. Zhang, and Z. Zhou (2025) Doubly optimal no-regret online learning in strongly monotone games with bandit feedback. Operations Research 73 (6), p. 3219–3244. Cited by: §1. [8] S. Bubeck, Q. Jiang, Y. T. Lee, Y. Li, and A. Sidford (2019) Near-optimal method for highly smooth convex optimization. In COLT, Cited by: §1.2. [9] B. Bullins and K. A. Lai (2022) Higher-order methods for convex-concave min-max optimization and monotone variational inequalities. SIAM Journal on Optimization 32 (3), p. 2208–2229. Cited by: §1.1, §1.2, §1, §2.1, §2.2, §4.1, §5.1, §5.2, §5. [10] X. Cai, A. Alacaoglu, and J. Diakonikolas (2024) Variance reduced halpern iteration for finite-sum monotone inclusions. In ICLR, Cited by: §1.1, §1.1, §2.3, §2.3, Proposition 2.2, §4.1, §4.1. [11] Y. Cai, A. Oikonomou, and W. Zheng (2022) Finite-time last-iterate convergence for learning in multi-player games. In NeurIPS, Cited by: §1.1, §2.3. [12] Y. Cai, A. Oikonomou, and W. Zheng (2024) Accelerated algorithms for constrained nonconvex-noncancave min-max optimization and comonotone inclusion. In ICML, Cited by: §1.1, §2.3. [13] Y. Cai and W. Zheng (2023) Accelerated single-call methods for constrained min-max optimization. In ICLR, Cited by: §1.1, §2.3. [14] Y. Carmon, D. Hausler, A. Jambulapati, Y. Jin, and A. Sidford (2022) Optimal and adaptive monteiro-svaiter acceleration. In NeurIPS, Cited by: §1.2, §1. [15] Y. Carmon and D. Hausler (2022) Distributionally robust optimization via ball oracle acceleration. In NeurIPS, Cited by: §1. [16] N. Cesa-Bianchi and G. Lugosi (2006) Prediction, learning, and games. Cambridge university press. Cited by: §1. [17] L. Chen, C. Liu, L. Luo, and J. Zhang (2025) Solving convex-concave problems with (ϵ−4/7)O(ε^-4/7) second-order oracle complexity. In COLT, Cited by: §1.2, §1, §3, §4.3, §6. [18] L. Chen and L. Luo (2024) Near-optimal algorithms for making the gradient small in stochastic minimax optimization. JMLR 25 (387), p. 1–44. Cited by: Lemma 5.2. [19] L. Chen, X. Zhang, C. Liu, J. Li, L. Luo, and J. Zhang (2026) Solving convex-concave problems with (ϵ−4/(3p+1))O(ε^-4/(3p+1)) ppth-order oracle complexity. arXiv preprint arXiv:2604.19462. Cited by: §1.2, §1, §3, §5.1, §5.1, §5.2, §6. [20] P. L. Combettes (2018) Monotone operator theory in convex optimization: pl combettes. Mathematical Programming 170 (1), p. 177–206. Cited by: §2.3. [21] J. Demmel, I. Dumitriu, and O. Holtz (2007) Fast linear algebra is stable. Numerische Mathematik 108 (1), p. 59–91. Cited by: §1, §2.2. [22] J. Diakonikolas (2020) Halpern iteration for near-optimal and parameter-free monotone inclusion and strong solutions to variational inequalities. In COLT, Cited by: §1.1, §2.3, §4.1, §4.1, §4.1. [23] R. Duan, H. Wu, and R. Zhou (2023) Faster matrix multiplication via asymmetric hashing. In FOCS, Cited by: §1, §2.2. [24] F. Facchinei and J. Pang (2003) Finite-dimensional variational inequalities and complementarity problems. Springer. Cited by: §1, §2.3. [25] A. Garg, R. Kothari, P. Netrapalli, and S. Sherif (2021) Near-optimal lower bounds for convex optimization for all orders of smoothness. In NeurIPS, Cited by: §1.2. [26] A. Gasnikov, P. Dvurechensky, E. Gorbunov, E. Vorontsova, D. Selikhanovych, and C. A. Uribe (2019) Optimal tensor methods in smooth convex and uniformly convexoptimization. In COLT, Cited by: §1.2. [27] F. Giannessi, A. Maugeri, et al. (1995) Variational inequalities and network equilibrium problems. Springer. Cited by: §1. [28] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio (2020) Generative adversarial networks. Communications of the ACM 63 (11), p. 139–144. Cited by: §1. [29] B. Halpern (1967) Fixed points of nonexpanding maps. Cited by: §1.1, §1.1. [30] P. Hartman and G. Stampacchia (1966) On some non-linear elliptic differential-functional equations. Cited by: §2.3. [31] K. Huang and S. Zhang (2025) An approximation-based regularized extra-gradient method for monotone variational inequalities. SIAM Journal on Optimization 35 (3), p. 1469–1497. Cited by: Appendix D, item 1, §1.1, §1.2, §1, §2.1, §2.2, §4.1, §4.2, §4.2, Lemma 4.2, §4, §5.1, §5.1, §5.2, §5, §5. [32] B. Jiang, H. Wang, and S. Zhang (2021) An optimal high-order tensor method for convex optimization. Mathematics of Operations Research 46 (4), p. 1390–1412. Cited by: §1.2. [33] H. Jiang, Y. T. Lee, Z. Song, and S. C. Wong (2020) An improved cutting plane method for convex optimization, convex-concave games, and its applications. In SIGACT, Cited by: §2.2. [34] R. Jiang, A. Kavis, Q. Jin, S. Sanghavi, and A. Mokhtari (2024) Adaptive and optimal second-order optimistic methods for minimax optimization. In NeurIPS, Cited by: §1.1, §1.2. [35] R. Jiang and A. Mokhtari (2025) Generalized optimistic methods for convex-concave saddle point problems. SIAM Journal on Optimization 35 (3), p. 2066–2097. Cited by: §1.2, §1, §2.1, §2.2, §4.2, §4.2. [36] M. Jordan, T. Lin, and Z. Zhou (2025) Adaptive, doubly optimal no-regret learning in strongly monotone and exp-concave games with gradient feedback. Operations Research 73 (3), p. 1675–1702. Cited by: §1. [37] D. Kinderlehrer and G. Stampacchia (2000) An introduction to variational inequalities and their applications. Cited by: §1. [38] G. M. Korpelevich (1976) The extragradient method for finding saddle points and other problems. Matecon 12, p. 747–756. Cited by: §1.1, §1, §4.2. [39] D. Kovalev and A. Gasnikov (2022) The first optimal acceleration of high-order methods in smooth convex optimization. In NeurIPS, Cited by: §1.2, §1. [40] S. Lee and D. Kim (2021) Fast extra gradient methods for smooth structured nonconvex-nonconcave minimax problems. In NeurIPS, Cited by: §1.1, §1.1. [41] F. Lieder (2020) On the convergence rate of the halpern-iteration. Optimization letters 15 (2), p. 405–418. Cited by: §1.1, §4.1. [42] T. Lin and M. I. Jordan (2024) Perseus: a simple high-order regularization method for variational inequalities. Mathematical Programming, p. 1–42. Cited by: item 1, §1.1, §1.2, §1, §1, §2.1, §2.2, §4.1, §4, §5.1, §5.1, §5.2, Lemma 5.1, §5. [43] T. Lin and M. I. Jordan (2023) Monotone inclusions, acceleration, and closed-loop control. Mathematics of Operations Research 48 (4), p. 2353–2382. Cited by: §1.2, §1. [44] G. J. Minty (1962) Monotone (nonlinear) operators in hilbert space. Cited by: §2.3. [45] A. Mokhtari, A. Ozdaglar, and S. Pattathil (2020) A unified analysis of extra-gradient and optimistic gradient methods for saddle point problems: proximal point approach. In AISTATS, Cited by: §1.1. [46] R. D. Monteiro and B. F. Svaiter (2012) Iteration-complexity of a newton proximal extragradient method for monotone variational inequalities and inclusion problems. SIAM Journal on Optimization 22 (3), p. 914–935. Cited by: item 1, §1.1, §1.2, §1, §1, §3, §4.1, §4.2, §4.3, §4, §5.1, §6. [47] R. D. Monteiro and B. F. Svaiter (2013) An accelerated hybrid proximal extragradient method for convex optimization and its implications to second-order methods. SIAM Journal on Optimization 23 (2), p. 1092–1125. Cited by: §1.2, §1.2, §1. [48] A. Nemirovski (2004) Prox-method with rate of convergence (1/t)O(1/t) for variational inequalities with lipschitz continuous monotone operators and smooth convex-concave saddle point problems. SIAM Journal on Optimization 15 (1), p. 229–251. Cited by: §1. [49] A. S. Nemirovskij and D. B. Yudin (1983) Problem complexity and method efficiency in optimization. Cited by: §1, §1. [50] Y. Nesterov and B. T. Polyak (2006) Cubic regularization of newton method and its global performance. Mathematical Programming 108 (1), p. 177–205. Cited by: §1.2, §2.2. [51] Y. Nesterov (1983) A method for solving the convex programming problem with convergence rate (1/k2)O(1/k^2). In Dokl akad nauk Sssr, Vol. 269, p. 543. Cited by: §1. [52] Y. Nesterov (2007) Dual extrapolation and its applications to solving variational inequalities and related problems. Mathematical Programming 109 (2-3), p. 319–344. Cited by: §1.1, §2.3. [53] Y. Nesterov (2008) Accelerating the cubic regularization of newton’s method on convex problems. Mathematical Programming 112 (1), p. 159–181. Cited by: §1.2. [54] Y. Nesterov (2018) Lectures on convex optimization. 137. Cited by: §1. [55] Y. Nesterov (2023) High-order reduced-gradient methods for composite variational inequalities. arXiv preprint arXiv:2311.15154. Cited by: §1.2, §1. [56] J. Nocedal and S. J. Wright (1999) Numerical optimization. Cited by: §1. [57] P. Ostroukhov, R. Kamalov, P. Dvurechensky, and A. Gasnikov (2020) Tensor methods for strongly convex strongly concave saddle point problems and strongly monotone variational inequalities. arXiv preprint arXiv:2012.15595. Cited by: §4.2, §5. [58] L. Qi and D. Sun (2002) Smoothing functions and smoothing newton method for complementarity and variational inequality problems. Journal of Optimization Theory and Applications 113, p. 121–147. Cited by: §2.2. [59] D. Ralph and S. J. Wright (2000) Superlinear convergence of an interior-point method despite dependent constraints. Mathematics of Operations Research 25 (2), p. 179–194. Cited by: §2.2. [60] R. T. Rockafellar (1976) Monotone operators and the proximal point algorithm. SIAM journal on control and optimization 14 (5), p. 877–898. Cited by: §1.1, §1.1, §2.3, §2.3. [61] E. K. Ryu, K. Yuan, and W. Yin (2019) Ode analysis of stochastic gradient methods with optimism and anchoring for minimax problems. arXiv preprint arXiv:1905.10899. Cited by: §1.1. [62] E. K. Ryu and S. Boyd (2016) Primer on monotone operator methods. Appl. comput. math 15 (1), p. 3–43. Cited by: §1.1, §2.3, Proposition 2.1. [63] Y. Ying, L. Wen, and S. Lyu (2016) Stochastic online auc maximization. In NeurIPS, Cited by: §1. [64] T. Yoon and E. K. Ryu (2021) Accelerated algorithms for smooth convex-concave minimax problems with (1/k2)O(1/k^2) rate on squared gradient norm. In ICML, Cited by: §1.1, §1.1. [65] B. H. Zhang, B. Lemoine, and M. Mitchell (2018) Mitigating unwanted biases with adversarial learning. In Proceedings of the 2018 AAAI/ACM Conference on AI, Ethics, and Society, p. 335–340. Cited by: §1. Appendix A Proof of Lemma 5.3 Proof. The first-order optimality conditions of ν1∗ y_ _1^* and ν2∗ y_ _2^* indicate that −ν1(ν1∗−0)∈(ν1∗)+(ν1∗)- _1( y_ _1^*- y_0)∈ G( y_ _1^*)+N_X( y_ _1^*) and −ν2(ν2∗−0)∈(ν2∗)+(ν2∗)- _2( y_ _2^*- y_0)∈ G( y_ _2^*)+N_X( y_ _2^*). Then, the monotonicity of + G+N_X gives ⟨−ν1(ν1−0)+ν2(ν2−0),ν1−ν2⟩≥0. - _1( y_ _1- y_0)+ _2( y_ _2- x_0), y_ _1- y_ _2 ≥ 0. Using ν1−0=(ν1−ν2)+(ν2−0) y_ _1- y_0=( y_ _1- y_ _2)+( y_ _2- y_0) and Lemma 5.2, we obtain (ν1−ν2)‖0−∗‖‖ν1−ν2‖≥ν1‖ν1−ν2‖2.( _1- _2)\| y_0- y^*\| \| y_ _1- y_ _2 \|≥ _1 \| y_ _1- y_ _2 \|^2. ∎ Appendix B Proof of Lemma 5.4 Proof. Note that we only need to prove the inequalities hold before rkr_k reaches ρ(μ)ρ(μ), or, equivalently, before μk _k reaches μ. Then, the truncated rkr_k and μk _k after that trivially satisfy the same inequalities. Proof of part 1. From the definition of rkr_k in equation 19, we have 1rk+1−1rk≤1−θp2R(p−1). 1r_k+1- 1r_k≤ 1- _p2R(p-1). Rearranging, we have rk+1rk≥(1+(1−θp)rk2R(p−1))−1 r_k+1r_k≥ (1+ (1- _p)r_k2R(p-1) )^-1 (22) Finally, using the inequality (1+x)−1≥1−x(1+x)^-1≥ 1-x for x≥0x≥ 0 and the facts rk≤Rr_k≤ R and p≥2p≥ 2 gives rk+1rk≥1−(1−θp)rk2R(p−1)≥1−1−θp2=1+θp2. r_k+1r_k≥ 1- (1- _p)r_k2R(p-1)≥ 1- 1- _p2= 1+ _p2. Proof of part 2. From the definition of μk _k and rkr_k in equation 19, we have 1−μk+1μk=1−(rk+1rk)p−11- _k+1 _k=1- ( r_k+1r_k )^p-1 Using equation 22 and then applying the inequality 1−(1+x)−(p−1)≤(p−1)x1-(1+x)^-(p-1)≤(p-1)x for x≥0x≥ 0, we obtain that 1−μk+1μj=1−(rk+1rk)p−1≤1−(1+(1−θp)rk2R(p−1))p−1≤(1−θp)rk2R1- _k+1 _j=1- ( r_k+1r_k )^p-1≤ 1- (1+ (1- _p)r_k2R(p-1) )^p-1≤ (1- _p)r_k2R ∎ Appendix C Proof of Theorem 5.2 Proof. Following the proof of Theorem 4.1, we show by induction that our parameter setting ensures the convergence rate of res(T):=‖T−η(T)‖η≤4Dη(T+1). res( x_T):= \| x_T- P_η( x_T)\|η≤ 4Dη(T+1). (23) The induction base of t=0t=0 follows from equation 15 by dividing η on both sides. Now, we assume equation 23 holds for all t≤T−1t≤ T-1. Then, Rt:=‖t−η(t)‖R_t:=\| x_t- P_η( x_t)\|, the initial distance of Algorithm 3 satisfies Rt≤4Dt+1,∀t=0,⋯,T−1.R_t≤ 4Dt+1, ∀ t=0,·s,T-1. Therefore, by Theorem 5.1, at each step the resolvent can be solved up to the approximation error δt _t in equation 9 under the parameter setting M=Θ(Lp)M= (L_p), Kt=Θ(Rt(ηLp)1/(p−1)))=Θ(TRtD)=(Tt+1),K_t= (R_t(η L_p)^1/(p-1)) )= ( TR_tD )=O ( Tt+1 ), and, by the definition of δt=Θ(Rt/(t+2log(t+2)) _t= (R_t/( t+2 (t+2)) in equation 9, St=Θ(loglog(Rt/δt))=Θ(loglog(t+2)).S_t= ( (R_t/ _t) )= ( (t+2) ). Therefore, the inexact condition in Lemma 4.1 is satisfied for t≤T−1t≤ T-1, which implies the convergence rate of equation 23 for t=Tt=T by the lemma. This completes the induction. Now, substituting our choices of η in equation 23 yields res(T)=(LpDp/Tp), res( x_T)=O (L_pD^p/T^p ), which is smaller than ϵε by our setting of T=(D(Lp/ϵ)1/p)T=O(D (L_p/ε )^1/p). Finally, the total number of oracle calls can be bounded by ∑t=0T−1Kt+St=(∑t=0T−1(Tt+1+loglog(t+2)))=(TlogT), _t=0^T-1K_t+S_t=O ( _t=0^T-1 ( Tt+1+ (t+2) ) )=O(T T), (24) which is equivalent to the bound of ~(D(Lp/ϵ)1/p) O(D (L_p/ε )^1/p) in equation 21 by substituting the choice of T. ∎ Appendix D Initial Algorithm and Analysis by AI The initial resolvent sub-solver given by AI has triple loops whose full procedure is presented in Algorithm 5. The initial proof by AI is also accessible through the following link: https://chatgpt.com/share/6a6c372c-6efc-83e8-a378-c6a1ea3b527 The authors simplified the algorithm and its proof by using the Restarted-NPE algorithm [31] for p=2p=2 and the single-loop Anchored-Tensor-Method introduced in Algorithm 3 for general p≥2p≥ 2. Algorithm 5 Triple-Looped-Anchored-Tensor-Method(,0,R,S,K)( G, y_0,R,S,K) 1:R0=R_0=R 2:for s=0,⋯,S−1s=0,·s,S-1 3: s,0=s y_s,0= y_s 4: for k=0,⋯,K−1k=0,·s,K-1 5: Define the anchoring coefficient μk=LpRsp−1(hp+kΔp)p−1,wherehp=12cp−1/(p−1),Δp=hp8(p−1),cp=2p(5p−2)p!. _k= L_pR_s^p-1(h_p+k _p)^p-1, where h_p= 12c_p^-1/(p-1),~~ _p= h_p8(p-1),~~c_p= 2^p(5p-2)p!. 6: Define the regularized operator k()=()+μk(−s) G_k( y)= G( y)+ _k( y- y_s) 7: s,k,0=s,k y_s,k,0= y_s,k 8: for j=0,⋯,5j=0,·s,5 9: Perform one tensor step s,k,j+1=kp(s,k,j;5Lp) y_s,k,j+1=T^p_ G_k( y_s,k,j;5L_p) 10: end for 11: s,k+1=s,k,6 y_s,k+1= y_s,k,6 12: end for 13: Change the anchor point to s+1=s,K y_s+1= y_s,K and contract the region by Rs+1=Rs/4R_s+1=R_s/4 14:end for 15:return S y_S