Exploring Jacobian Inexactness in Second-Order Methods for Variational Inequalities: Lower Bounds, Optimal Algorithms and Quasi-Newton Approximations

Authors: artem  agafonov
petr  ostroukhov
roman  mozhaev
konstantin  yakovlev
eduard  gorbunov
martin  takáč
alexander  gasnikov
dmitry  kamzolov
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 19:06:04

Abstract: Variational inequalities represent a broad class of problems, including minimization and min-max problems, commonly found in machine learning. Existing second-order and high-order methods for variational inequalities require precise computation of derivatives, often resulting in prohibitively high iteration costs. In this work, we study the impact of Jacobian inaccuracy on second-order methods. For the smooth and monotone case, we establish a lower bound with explicit dependence on the level of Jacobian inaccuracy and propose an optimal algorithm for this key setting. When derivatives are exact, our method converges at the same rate as exact optimal second-order methods. To reduce the cost of solving the auxiliary problem, which arises in all high-order methods with global convergence, we introduce several Quasi-Newton approximations. Our method with Quasi-Newton updates achieves a global sublinear convergence rate. We extend our approach with a tensor generalization for inexact high-order derivatives and support the theory with experiments.

Main paper: 
 introduction: 
 In this paper, we primarily address the problem of solving the Minty Variational Inequality (MVI) [77, 17] . Given a continuous operator F : X → R d , where X ⊆ R d is a closed bounded convex subset with a diameter D = max x,y∈X ∥x -y∥, the objective is to find a point x * ∈ X such that ⟨F (x), x -x * ⟩ ≥ 0, for all x ∈ X .(1) The solution to (1) is referred to as a weak solution of the Variational Inequality (VI) [37] . In contrast, the Stampacchia variational inequality problem [48] consists in finding a point x * ∈ X such that ⟨F (x * ), x -x * ⟩ ≥ 0, for all x ∈ X .(2) This solution is often called a strong solution to the variational inequality. When the operator F is both continuous and monotone, the weak and strong solutions are equivalent [37] .Assumption 1.1 The operator F (x) is called monotone, if ⟨F (x) -F (y), x -y⟩ ≥ 0, for all x, y ∈ X .Another useful assumption is L 1 -smoothness.Assumption 1.2 The operator F (x) is L 1 -smooth, if it has Lipschitz-continuous first-order derivative ∥∇F (x) -∇F (y)∥ op ≤ L 1 ∥x -y∥, for all x, y ∈ X .First-order methods. Variational inequalities encompass a wide range of problems, including minimization [19, 90, 93, 15] , min-max problems, Nash equilibrium, differential equations, and others [37, 16] . The extensive research on VI methods dates back several decades, with a notable breakthrough in the 1970s-the development of the Extragradient method [62, 6] . Subsequently, it was demonstrated that this method achieves global convergence of O ε -1 [80, 43] , matching the convergence rates of other first-order2 methods such as optimistic gradient [92, 78, 63, 44] , forward-backward splitting [95] , and dual extrapolation [81] . These first-order methods collectively exhibit optimal convergence [89] .Second-order and high-order methods. To achieve further notable acceleration of methods for VIs, one can leverage information about higher-order derivatives. For instance, simply incorporating firstorder derivatives (Jacobian) can significantly enhance the convergence speed of the method. Following recent advancements in second-order and high-order methods with global rates for minimization [85, 82, 10, 79, 83, 46, 58, 14, 51, 64, 22] , several high-order methods for VIs have been proposed [21, 68, 52, 88, 70, 84, 72, 5] . However, all these methods involve a line-search procedure, resulting in Õ ε -2/3 convergence for the case of first-order information (Jacobians). Recent works [69, 2] propose methods with improved rates O ε -2/3 and establish the lower bound Ω ε -2/3 , rendering these algorithms optimal.Jacobian's approximation. In the last decade, VIs found new applications in machine learning.There are many problems that could not be reduced to minimization, including reinforcement learning [87, 55] , adversarial training [75] , GANs [42, 31, 41, 76, 28, 67, 91] , classical learning tasks in supervised learning [56, 9] , unsupervised learning [98, 8] , image denoising [35, 27] , robust optimization [13] . Applying second-order methods, without even mentioning high-order ones, described in a previous paragraph to machine learning problems could be a challenging task. Although these methods may theoretically converge faster, computing exact Jacobians and the per-iteration costs can be expensive. Therefore, it seems natural to introduce inexact approximations of the first-order derivatives. In the context of minimization, several works with inexact Hessians were introduced for both convex [40, 3, 7, 4] and nonconvex [24, 25, 23, 61, 99, 94, 74, 11, 12, 34] problems. Regarding VIs, Quasi-Newton (QN) methods [73, 14, 57] can be highlighted, though they, unfortunately, achieve only local convergence in the strongly monotone case [18, 38] . These methods are relatively less advanced for VIs compared to their counterparts in the field of minimization, where they are considered classics in optimization due to their effectiveness and practicality [86] . Modern research on QN approximations for minimization includes methods that exhibit global convergence [57, 54, 53] . A recent work [71] introduces the Newton-MinMax method for convexconcave unconstrained min-max optimization problems, demonstrating an optimal rate under special assumptions on the accuracy of the Jacobian approximation. However, the field of VIs lacks globally convergent inexact second-order methods with an explicit dependence on the accuracy of the Jacobian. This raises several natural questions:What are the lower bounds for methods with inexact Jacobians? Can we construct an optimal method with inexact first-order information? What is the proper way to approximate the Jacobian to ensure global convergence and reduce the iteration complexity?In our work, we attempt to answer these questions in a systematic manner.Optimality measure. Most of our results are stated for the monotone setting (Assumption 1.1). In this context, the optimality of a point x ∈ X is typically measured by a gap function GAP(•) : X → R + [95, 80, 81, 78, 69] , defined byEQUATIONwhere ε ≥ 0 is the accuracy of solution. The boundedness of X and the existence of a strong solution ensure that the gap function is well-defined. If ε = 0, we get by (1) that x is a weak solution of VI.We explore the performance of the proposed algorithm in scenarios involving nonmonotone operators F . In such cases, it is essential to assume that the operator satisfies the Minty condition to ensure that the problem is computationally manageable [32] .Assumption 1. 3 The operator F (x) satisfies Minty condition, if there exists a point x * such that ⟨F (x), x -x * ⟩ ≥ 0, for all x ∈ X .(EQUATIONThe range of applications of nonmonotone VIs satisfying Minty conditions is quite extensive [20, 29, 39, 36, 66, 60] . We note, that this condition is weaker than monotonicity [30, 50, 59 ] and guarantees the existence of at least one strong solution since F is continuous and X is closed and bounded [47] .To measure the optimality of point x we define the residue function RES(•) : X → R + [30, 50, 59 ]EQUATIONThe boundedness of X and the existence of a strong solution ensure that the residual function is well-defined. If ε = 0, by (2), we get that x is a strong solution of VI.Contributions. The main contribution of this paper lies in the development of a new second-order method robust to inexactness in the Jacobian, a common occurrence in machine learning. We demonstrate the algorithm's optimality in the monotone case by establishing a lower bound for this key setting. Expanding further:1. We introduce a novel second-order algorithm, VIJI (Second-order Method for Variational Inequalitues under Jacoibian Inexactness), designed to handle δ-inexact 3 Jacobian information. Specifically, in the context of smooth and monotone VIs, VIJI achieves a convergence rate of O δD 2 T + L1D Notably, when δ ≤ L1D √ T , our method matches the convergence rates of optimal exact second-order methods [2, 69] . 2. We establish the optimal performance of our algorithm on monotone smooth operators by deriving a theoretical complexity lower bound of Ω δD 2 T + L1D 3T 3/2to find weak solution for the case of δ-inexact Jacobians. 3. Our algorithm involves solving a variational inequality subproblem. To tackle this challenge, we introduce an approximation condition, which makes the solution computationally feasible. 4. We introduce a new Quasi-Newton update for approximating the Jacobian, which significantly decreases the per-iteration cost of the algorithm while maintaining a global sublinear convergence rate. Numerical experiments demonstrate the practical benefits of our method. 5. We extend our algorithm for higher-order VIs with inexact high-order derivatives, resulting inO p-1 i=1 δiD i+1 T (i+1)/2 + Lp-1D p+1 T (p+1)/2rate for monotone and smooth VIs with δ i -inexact i-th derivative to find weak solution. Moreover, we extend our proposed high-order method to nonmonotone VIs. 6. We propose a restarted version of VIJI for strongly monotone VIs, which exhibits a linear rate.Comparison with Lin, Mertikopoulos, and Jordan [71] . To the best of our knowledge, the work [71] is the most closely related to our research.The objective of [71] was to develop a method for convex-concave unconstrained min-max optimization under inexact Jacobian information with an optimal convergence rate O(ε -2/3 ) matching the lower bound Ω(ε -2/3 ) [69] . The authors successfully achieve this goal by proposing a second-order algorithm Inexact-Newton-MinMax method based on the Perseus [69] . To attain optimal convergence, they constrained the Jacobian inexactness with a function that decreases as the method converges and bounded the norm of Jacobian from above. While these assumptions might be suitable for randomized sampling in finite-sum and stochastic problems, they may not hold for many approximation strategies, such as Quasi-Newton algorithms. The aim of our work, however, is to study the impact of Jacobian inaccuracy on the convergence of second-order methods for VIs (a special case of which are min-max problems) and to identify the explicit dependence of the convergence rate on the inexactness. Compared to the work [71] , the Jacobian inaccuracy directly affects the step sizes in our algorithm, allowing us to achieve a convergence rate of O(ε -2/3 +δε -1 ) for any given δ. VIJI can be viewed as a generalization of Inexact-Newton-MinMax. With the same assumption on δ as in [71] , our methods for min-max optimization are equivalent. The inexactness of the subproblem and the solution approach proposed in [71] remain valid for our method even with arbitrary large δ. Further details about application of our method to min-max problems can be found in Appendix J.Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?Answer: [Yes]We proposed an algorithm for variational inequalities with inexact Jacobians in Section 3, supported its optimality by establishing the lower bound in Section 4, introduced Quasi-Newton updates in Section 5, studied strongly-convex case in Section 6 and discussed tensor generalization in Section 7. Section 8 provides experimental results to validate our proposed approach.Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper. 
preliminaries: 
 Notation. Let R d be a finite-dimensional vector space with scalar product ⟨•, •⟩. For vector x ∈ R d we denote Euclidean norm as ∥x∥. For X ∈ R d1×...×dp , we defineX[z 1 , • • • , z p ] = 1≤ij ≤dj ,1≤j≤p (X i1,••• ,ip )z 1 i1 • • • z p ip ,and∥X∥ op = max ∥z i ∥=1,1≤j≤p X[z 1 , • • • , z p ].Fixing p ≥ 1 and letting F : R d → R d be a continuous and high-order differentiable operator, we define ∇ (p) F (x) as the p th -order derivative at a point x ∈ R d . To be more precise, letting z 1 , . . . , z k ∈ R d , we have∇ (k) F (x)[z 1 , • • • , z k ] = 1≤i1,...,i k ≤d ∂Fi 1 ∂xi 2 •••∂xi k (x) z 1 i1 • • • z k i k .Taylor approximation and oracle feedback. The starting point for our method is the first-order Taylor polynomial of the operatorF at point v : Φ v (x) = F (v) + ∇F (v)[x -v].Since the computation of Jacobian ∇F (v) could be a quite tiresome task, it seems natural to introduce an inexact approximation J(v). Based on it we introduce inexact Taylor approximation -one of the main building blocks of our algorithmEQUATIONwhere J(x) satisfies the following assumption.Assumption 2.1 For given v ∈ X δ-inexact Jacobian satisfies:EQUATIONAs it was shown, e.g. in [52] , Assumption 1.2 allows to control the quality of approximation of operator F by its Taylor polynomialEQUATIONThe next lemma is counterpart of (9) for the case of inexact Jacobian.Lemma 2.2 Let Assumptions 1.2 and 2.1 hold. Then, for anyx, v ∈ X ∥F (x) -Ψ v (x)∥ ≤ L1 2 ∥x -v∥ 2 + δ∥x -v∥.Question: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes]Justification: the assumptions used in this paper are classical for methods for variational inequalities [69, 21] . Our methods include the solution of auxilarary subproblem (as other globally convergent second-order methods for VIs [69, 21] ), which requires subsolver or addditional subroutine, what we explicitly stated. Moreover, we introduced the way to reduce the complexity and provided the computational cost of the proposed procedure.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
viji algorithm: 
 In this section, extending on recent optimal high-order method for MVIs Perseus [69] , we present our proposed method, dubbed as VIJI and detailed in Algorithm 1.Input: initial point x 0 ∈ X , parameters L 1 , η, sequence {β k }, and opt ∈ {0, 1, 2}.Initialization: set s 0 = 0 ∈ R d . for k = 0, 1, 2, . . . , T do Compute v k+1 = argmax v∈X {⟨s k , v -x 0 ⟩ -1 2 ∥v -x 0 ∥ 2 }. Compute x k+1 ∈ X such that condition (12) holds true. Compute λ k+1 such that 1 32 ≤ λ k+1 L1 2 ∥x k+1 -v k+1 ∥ + β k+1 ≤ 1 22 . Compute s k+1 = s k -λ k+1 F (x k+1 ). Output: x =    xT = 1 T k=1 λ k T k=1 λ k x k , if opt = 0, x T , else if opt = 1, x k T for k T = argmin 1≤k≤T ∥x k -v k ∥, else if opt = 2.The model of objective and subproblem's solution. We begin the description of the algorithm by introducing the inexact model of objectiveEQUATIONwhere η > 0 is given constant. Here, we introduced the additional regularization term ηδ(x -v). As we will demonstrate, this term is crucial for ensuring that the method's subproblem has a solution.For brevity, we use a regularization constant of 5L. Using larger coefficients would yield the same convergence rate. As other dual extrapolation-type methods, VIJI includes the following subproblemEQUATIONFirst of all, the strong solution of this VI exists because Ω η v k+1 (x) is continuous, and X is a closed, bounded, and convex set. Next, we demonstrate that in a monotone setting VI ( 11) is also monotone.Lemma 3.1 Let Assumptions (1.1), (1.2), (2.1) hold. Then for any x, v k+1 ∈ X VI (11) is monotone 1 2 ∇Ω v (x) + ∇Ω v (x) T ⪰ 4L 1 ∥x -v∥I d×d + 5L 1 (x-v)(x-v) T ∥x-v∥ + (η -1)δI d×d .Following the work [69] , one can find a strong solution to such VI using mirror-prox methods from [1] , achieving the following approximate conditionEQUATION)This ensures that the subproblem is computationally solvable in the monotone setting. In specific cases, such as minimax optimization, other efficient subsolvers can be employed [49, 2, 71] .Adaptive dual stepsizes. Adaptive stepsizes in dual space λ k are another core aspect of the algorithm. Due to inaccuracies in the Jacobian, applying the standard adaptive strategy, such as in the Perseus algorithm [69] , can lead to excessively large steps, potentially slowing down the method.To address this issue, an additional term β k is incorporated into the adaptive strategy for selecting λ k . Thus, when ∥x k -v k ∥ is small (indicating proximity to the optimum), β k has a greater influence on the choice of λ k , preventing the method from taking overly aggressive steps. Similar behavior can be observed in accelerated second-order methods for minimization with inexact Hessians from theoretical [3, Lemma 5], [4, Appendix B, Lemma E.3] and practical [4, Section 7] perspectives.Convergence in monotone setting. Now, we are prepared to present the convergence theorem of Algorithm 1 in the monotone case.Theorem 3.2 Let Assumptions 1.1, 1.2, 2.1 hold. Then, after T ≥ 1 iterations of Algorithm 1 with parameters β k = δ, η = 10, opt = 0, we get the following boundEQUATIONThe upper bound (13) consists of two terms. The first one corresponds to exact convergence and matches the lower bound for second-order VI methods [69] . The second term illustrates the impact of the Jacobian's inexactness on the convergence rate, aligning with the lower bound for first-order methods [89] . We also notice that one can make the bound from (13) tighter for the so-called restricted gap-function [81] defined as RGAP(y) = sup x∈C ⟨F (x), y -x⟩, where C ⊆ X and the solution set X * of (2) satisfies X * ⊆ C. In particular, following similar steps as in the proof of Theorem 3.2, one can deriveO L1 D 3 T 3/2 + δ D 2T bound for RGAP(x T ), where D = sup x∈C ∥x -x 0 ∥.In the next theorem, we show that the method can achieve the optimal convergence rate of secondorder methods under additional assumption on δ. (v k ) -J(v k ))[x k -v k ]∥ ≤ δ k ∥x k -v k ∥, δ k ≤ L1 2 ∥x k -v k ∥. (14) Then, after T ≥ 1 iterations of Algorithm 1 with parameters β k = L1 2 ∥x k -v k ∥, η = 10, opt = 0, we get the following boundGAP(x T ) = sup x∈X ⟨F (x), xT -x⟩ ≤ O L1D 3 T 3/2 .Note, that condition ( 14) is verifiable, indicating that the method can adjust to δ k in cases of controllable inexactness. Specifically, at each iteration, we can solve subproblem (12) . If the assumption regarding δ k is not met, we can improve Jacobian approximation and repeat procedure. Moreover, similarly to the previous theorem, one can tighten the above bound for RGAP(x T ) and get the dependence on D = sup x∈C ∥x -x 0 ∥ instead of D.Convergence in nonmonotone setting. To begin with, in the nonmonotone case, the subproblem (11) may not exhibit monotonicity, and solving (12) becomes challenging [32] . Yet, in certain specific scenarios, such as unconstrained minimization tasks, it remains feasible to find a solution by leveraging the cubic structure of the subproblem [26] . The following theorem establishes the convergence of VIJI in the nonmonotone setting. .The convergence rate could be improved by additional assumption on inexact Jacobian.Theorem 3.5 Let Assumptions 1.2, 1.3, 2.1 hold. Let {x k , v k } be iterates generated by Algorithm 1 that satisfy (14) . Then after T ≥ 1 iterations of Algorithm 1 with parametersβ k = L 2 ∥x k -v k ∥, η = 10, opt = 2 we get the following bound RES(x) = sup x∈X ⟨F (x T ), xT -x⟩ = O L1D 3 T .Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof? Answer: [Yes] Justification: For all theoretical results, we provide proofs in the main paper and Appendix. Guidelines:• The answer NA means that the paper does not include theoretical results.• All the theorems, formulas, and proofs in the paper should be numbered and crossreferenced. • All assumptions should be clearly stated or referenced in the statement of any theorems.• The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition. • Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material. • Theorems and Lemmas that the proof relies upon should be properly referenced. 
the lower bound: 
 In this section, we establish a theoretical lower bound for the complexity of first-order algorithms using inexact Jacobians for monotone MVIs. The proof technique draws inspiration from works [33, 4] and based on lower bounds from [69, 89] .We start by describing the available information and the method's structure. The considered class of algorithms relies on data provided by a first-order δ-inexact oracle, denoted as O :EQUATION)The method is able to generate points {x k } k≥0 that satisfy the following conditions ∈ Lin(F (x 0 ), . . . , F (x k )), x = argmax x∈X {⟨s, x -x 0 ⟩ -1 2 ∥x -x 0 ∥ 2 }, x k+1 ∈ X satisfies that ⟨Ω x(x k+1 -x k ), x -x k+1 ⟩ ≥ 0 for all x ∈ X , where Ω x(h) = a 1 F (x) + a 2 J(x)[h] + b 1 h + b 2 ∥h∥h.Next, we state the primary assumption concerning the method's ability to generate new points. Assumption 4.1 The method generates a recursive sequence of iterates {x k } k≥0 that satisfies the following condition: for all k ≥ 0, we have that x k+1 ∈ X satisfies that ⟨Ω x(x k+1 -x k ), x-x k+1 ⟩ ≥ 0 for all x ∈ X , wherex = argmax x∈X {⟨s, x -x 0 ⟩ -1 2 ∥x -x 0 ∥ 2 } and s ∈ Lin(F (x 0 ), . . . , F (x k )).As highlighted in [69] , Assumption 4.1 is suitably satisfied by various dual extrapolation methods. However, it might not be applicable to alternative methods for variational inequalities, such as extragradient methods and their variants. We leave the generalization of inexact lower bounds for these algorithms to future research, as even lower bounds for exact algorithms [2, 69] do not address this case. Now, let us introduce the generalization of smoothnessEQUATIONFinally, we present the lower bound theorem for first-order methods with inexact Jacobians.Theorem 4.3 Let some first-order method M satisfy Assumption 4.1 and have access only δ-inexact first-order oracle 15. Assume the method M ensures for any L 0 -zero-order smooth and L 1 -first-order smooth monotone operator F the following convergence rateEQUATION)Then for all T ≥ 1 we haveEQUATION)5 Quasi-Newton ApproximationIn this section, inspired by Quasi-Newton (QN) methods for Hessian approximation, we propose QN approximations for Jacobians. Our goal is to create a simple scheme to approximate the first-order derivative and thereby reduce the complexity of the subproblem. We compute J x using a QN update and use it as an inexact Jacobian in the model Ω η v (x) (10) .EQUATION)where r is a rank of approximation,J 0 ∈ R d×d ⪰ 0, u i ∈ R d , v i ∈ R d are known. U ∈ R r×d and V ∈ R r×d are matrices of stacked vectors U = [u 0 , . . . , u r-1 ] and V = [v 0 , . . . , v r-1 ] and C ∈ R r×r is a diagonal matrix C = diag([c 0 , . . . , c r-1 ]). If u i = v i the update becomes symmetric.L-Broyden is a non-symmetric variant of QN approximation [45, 38] of the following formEQUATIONIn a view of (19) , this update is obtained by settingu i = y i -J i s i , v i = s i , c i = 1/(s ⊤ i s i ),and the rank r = m is equal to memory size m.Damped L-Broyden is another option for QN approximation with non-symmetric damped updateEQUATIONBy choosing u i = y i -J i s i , v i = s i , c i = 1/((m + 1)s ⊤ i s i ), r = m, we derive this update from (19) . We define the matrix J r (J 0 , U, V, C) = J r (J 0 , Y, S, C), where Y and S are formed by stacking the vectors [y 0 , . . . , y m-1 ] and [s 0 , . . . , s m-1 ]. The matrix J m (J 0 , Y, S) can be computed for any given pair (Y , S). Next, we describe two strategies for the choice of (s, y) pairs used in (20) , (21) .QN with operator history is the well-known classic variant where operator differences are stored:s i = z i+1 -z i , y i = F (z i+1 ) -F (z i ).This approach is computationally efficient as it does not require additional operator calculations.QN with JVP sampling is based on fast computation of Jacobian-Vector Products (JVP):y i = ∇F (x)s i ,where s i are random vectors uniformly distributed on the unit sphere such that ∥s i ∥ = 1 and s 0 , . . . , s m-1 are linearly independent. Note, for m ≪ d, each s i is linearly independent with high probability. This approach requires only m operator/JVP computations per step, which is considerably fewer than the d JVPs needed for a full Jacobian. Utilizing the current Jacobian information allows to improve the accuracy of the approximation.In the following theorem, we demonstrate that these approximations satisfy Assumption 2.1 and condition (8) for both the QN with operator history and JVP sampling methods.Theorem 5.1 Let F (x) be L 0 -zero-order smooth operator. For m-memory L-Broyden approximation of the Jacobian J x = J m defined iteratively by (20) with 0 ⪯ J 0 ⪯ L 0 I, we have δ ≤ (m + 2)L 0 .For m-memory Damped L-Broyden approximation J x = J m of the Jacobian defined iteratively by (21) with 0 ⪯ J 0 ⪯ L0 m+1 I, the condition δ ≤ 2L 0 holds true.With the primary toolkit for QN approximation in VIs established, we can now discuss efficient way of solving the subproblem (11) , which takes the following formEQUATION)Let us introduce a parameter τ = ∥y -x∥ for a segment search in τ ∈ [0; D]. To solve (22) , we consider another problemEQUATIONwhere 22) and ( 23) are equivalent when τ = ∥y τ -x∥. The subproblem ( 23) can be reformulated as minimization problemA τ = J x + (ηδ + 5L 1 τ )I. Problems (EQUATIONwhere ( 23) is an optimality condition for (24) . The goal is to find y τ such that υ(τ ) = |τ -∥y τ -x∥| ≤ ε. As υ(τ ) is a continuous function of τ , we can find this solution via bisection segment-search with log 2 D ε iterations. This ray-search procedure is similar to the subproblem solution for the Cubic Regularized Newton subproblem.For r-rank QN approximation J r from (19) , we can effectively compute A -1 τ F (x), where A τ = J r + (ηδ + 5L 1 τ )I = U ⊤ CV + J 0 + (ηδ + 5L 1 τ )I = U ⊤ CV + B and B = J 0 + (ηδ + 5L 1 τ )I by using the Woodbury matrix identity [96, 97] .A -1 τ F (x) = B + U ⊤ CV -1 F (x) = B -1 F (x) -B -1 U ⊤ (C -1 + V B -1 U ⊤ ) -1 V B -1 F (x).For computational efficiency, it is better to choose J 0 as a diagonal matrix, then inversionB -1 is computed by O(d) arithmetical operations, C -1 by O(r) operations. V B -1 U ⊤ requires O(r 2 d) for classical multiplication and can be improved by fast matrix multiplication.(C -1 + V B -1 U ⊤ ) -1can be computed by O(r 3 ), as an inverse of r-rank matrix. The rest of the operations are cheaper. Thus, the total number of arithmetic operations is O(r 2 d) instead of O(d 3 ) for Jacobian inversion. We need to perform this inversion logarithmic number of times. Therefore, the total computational cost with the segment-search procedure is Õ r2 d + r 3 log 2 ( D ε ) .6 Strongly monotone setting Assumption 6.1 The operatorEQUATIONTo leverage the strong monotonicity of the objective function and achieve a linear convergence rate, we introduce the restarted version of Algorithm 1 dubbed as VIJI-Restarted. Restart techniques are widely utilized in optimization and typically preserve optimality. In other words, restarting an optimal method for convex functions or monotone problems effectively transforms it into an optimal method for strongly convex or strongly monotone problems. Within iteration of VIJI-Restarted listed as Algorithm 2, we execute VIJI for a predefined number of iterations (26) . Next, the output of this run is used as initial point for next run of VIJI with parameters reset, and this iterative process continues.Input: initial point z 0 ∈ X , D = max x,y∈X , parameters L, δ.Initialization: n = log D ε , R = D. for i=1, ..., n do Set x 0 = z i-1 , R i-1 = R 2 i-1 .Run Algorithm 1 with β k = δ, η = 10, opt = 0 for T i iterations, whereEQUATION)Setz i = x Ti .Theorem 6.2 Let Assumptions 1.2, 2.1, 6.1 hold. Then the total number of iterations of Algorithm 2 to reach desired accuracy ∥z s -x * ∥ ≤ ε, where x * is the solution of (2) isO L1D µ 2 3 + δ µ + 1 log D ε .To make our paper more self-contained, we provide convergence rate of Algorithm 3 in nonmonotone case.Theorem 7.3 Let Assumptions 4.2, 1.3 and 7.1 with i = p -1 hold. Then after T ≥ 1 iterations of Algorithm 3 with parameters η = 5p, opt = 2 we get the following boundRES(x) := sup x∈X ⟨F (x t ), xt -x⟩ = O Lp-1D p+1 T p 2 + p-1 i=1 δiD i+1 T i 2 . Proof. Consider ⟨F (x k ), x k -x⟩. ⟨F (x k ), x k -x⟩ = ⟨F (x k ) -Ω p,v k (x k ), x k -x⟩ + ⟨Ω p,v k (x k ), x k -x⟩ (53) ≤ ∥F (x k ) -Ω p,v k (x k )∥∥x k -x∥ + L p-1 p! ∥x k -v k ∥ p+1 + p-1 i=1 δ i i! ∥x k -v k ∥ i+1 .From definition of Ψ p,v k (x k ) and triangle inequality we get∥F (x k ) -Ψ p,v k (x k )∥ = ∥F (x k ) -Ω p,v k (x k )∥ - 5L p-1 (p -1)! ∥x k -v k ∥ p - p-1 i=1 η i δ i i! ∥x k -v k ∥ i .From this and (55) we getEQUATION≤ L p-1 p! ∥x k -v k ∥ p + p-1 i=1 1 i! δ i ∥x k -v k ∥ i + 5L p-1 (p -1)! ∥x k -v k ∥ p + p-1 i=1 η i δ i i! ∥x k -v k ∥ i = L p-1 p! (5p + 1)∥x k -v k ∥ p + p-1 i=1 δ i i! (η i + 1)∥x k -v k ∥ i . Thus, ⟨F (x k ), x k -x⟩ ≤ L p-1 p! (5p + 1)∥x k -v k ∥ p ∥x k -x∥ + p-1 i=1 δ i i! (η i + 1)∥x k -v k ∥ i ∥x k -x∥ + L p-1 p! ∥x k -v k ∥ p+1 + p-1 i=1 δ i i! ∥x k -v k ∥ i+1 ≤ L p-1 p! ∥x k -v k ∥ p (5p + 2)D + p-1 i=1 δ i i! ∥x k -v k ∥ i (η i + 2)D.From second inequality of (38) and definition of x k T we get∥x k T -v kt ∥ 2 ≡ min 1≤k≤T ∥x k -v k ∥ 2 ≤ 1 T T k=1 ∥x k -v k ∥ 2 ≤ 4 T ∥x * -x 0 ∥ 2 ≤ 4D 2 T ⇔ ∥x k T -v k T ∥ ≤ 4 i 2 D i T i 2.Combining all these results, we getRES(x) = sup x∈X ⟨F (x k ), x k -x⟩ ≤ L p-1 p! ∥x k -v k ∥ p (5p + 2)D + p-1 i=1 δ i i! ∥x k -v k ∥ i (5p + 2)D ≤ 2 p (5p + 2) p! L p-1 D p+1 T p 2 + p-1 i=1 2 i δ i i! (η i + 2) D i+1 T i 2 = O L p-1 D p+1 T p 2 + p-1 i=1 δ i D i+1 T i 2 □Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)? Answer: [Yes] Justification: our algorithm is listed in the paper (Algorithm 1) and the details are fully described (Sections 3 and 5). The expereimantal details for reproducibility are describes in Section 8 and Appendix I. The code is attached. Guidelines:• The answer NA means that the paper does not include experiments.• If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not. • If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable. • Depending on the contribution, reproducibility can be accomplished in various ways.For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed. • While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors).• It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
tensor generalization: 
 In this section, we generalize the results presented in Section 3 to the p-th order case. We consider a higher-order method with inexact high-order derivatives, satisfying the following assumption.EQUATIONBased on inexact (p -1)-th-order Taylor approximationEQUATIONThe tensor generalization of Algorithm 1, referred to as VIHI (High-order Method for Variational Inequalitues under High-order derivatives Inexactness ) and detailed in Appendix G, involves the inexact solution of the subproblem, which satisfies the following condition:sup x∈X ⟨Ω v k+1 (x k+1 ), x k+1 -x⟩ ≤ Lp-1 p! ∥x k+1 -v k+1 ∥ p+1 + p-1 i=1 δi i! ∥x k+1 -v k+1 ∥ i+1 .Another difference in VIHI compared to VIJI (Algorithm 1) is the adaptive strategy for λ k+1 :1 4(5p-2) ≤ λ k Lp-1 p! ∥x k+1 -v k+1 ∥ p+1 + p-1 i=1 δi i! ∥x k+1 -v k+1 ∥ i+1 ≤ 1 2(5p+1) .The other steps of Algorithm 1 remain unchanged for the higher-order method. Now, we are ready to present the convergence properties of VIHI.Theorem 7.2 Let Assumptions 1.1, 4.2 with i = p -1, and 7.1 hold. Then, after T ≥ 1 iterations of VIHI with parameters η i = 5p, opt = 0, we get the following boundGAP(x T ) = sup x∈X ⟨F (x), xT -x⟩ ≤ O Lp-1D p+1 T p+1 2 + p-1 i=1 δiD i+1 T i+1 2.Finally, we extend our tensor generalization to nonmonotone case and obtain the following result.Theorem 7.3 Let Assumptions 1.3, 4.2 with i = p -1, and 7.1 hold. Then after T ≥ 1 iterations of VIHI with parameters η = 5p, opt = 2 we get the following boundRES(x) := sup x∈X ⟨F (x t ), xt -x⟩ = O Lp-1D p+1 T p 2 + p-1 i=1 δiD i+1 T i 2 . 
experiments: 
 In this section, we present numerical experiments to demonstrate the efficiency of our proposed methods. We consider the cubic regularized bilinear min-max problem of the form:min x∈R d max y∈R d f (x, y) = y ⊤ (Ax -b) + ρ 6 ∥x∥ 3 ,where ρ > 0, b = [1, 0, . . . , 0] ∈ R d , and A ∈ R d×d , with all 1 on the main diagonal and all -1 on the upper diagonal, the rest elements are 0. To reformulate it as variational inequality, we defineF (x) = [∇ x f (x, y), -∇ y f (x, y)].This problem is inspired by the first-order lower bound function for variational inequalities and min-max problems and is commonly used to verify the convergence of high-order methods for VI [71, 52] .We implement our second-order method for Variational Inequalities with Quasi-Newton Approximation (VIQA) as a PyTorch optimizer. The code is available in the OPTAMI package4 [58] .VIQA Broyden refers to L-Broyden approximation (20) and VIQA Damped Broyden to (21) , which are used as inexact Jacobians in VIJI (Algorithm 1). We compare them with the Extragradient method (EG) [62] , first-order Perseus (Perseus1), and second-order Perseus with Jacobian (Perseus2). In Figure 1 , one can see that second-order information in Perseus2 significantly accelerates the convergence compared to Perseus1. However, it is expensive to compute Jacobian and solve secondorder subproblems every iteration. VIQA with the proposed Damped Broyden approximation ( 21) is significantly faster than EG, Perseus1, and VIQA Broyden (20) . It shows that this approximation improves the convergence of first-order Perseus1 and confirms the theoretical result that Damped Broyden is a more accurate approximation than classical Broyden from Theorem 5.1. The detailed parameters and setup are presented in Appendix I.Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes] Justification: we describe compute resourses in Appendix I.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
conclusion: 
 In this work, we introduced a second-order method specifically designed to handle Jacobian inexactness. We demonstrated its optimality in the monotone case by introducing a new lower bound and extended its applicability to tensor methods. However, similar to other high-order methods with global convergence properties, our algorithm involves a subproblem that necessitates an additional subroutine for its solution. To address this challenge, we proposed a computationally feasible criterion for solving the subproblem and implemented Quasi-Newton approximations for Jacobians, resulting in a significant reduction in per-iteration cost. Future investigations could explore incorporating inexactness within the operator itself and developing adaptive schemes to dynamically adjust for the level of inexactness encountered during the optimization process. Another open problem is a design of more accurate Quasi-Newton approximations specifically for Jacobians, focusing on non-symmetrical structure of Jacobian and specific inexactness criteria such as Assumption 2.1. Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines? 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?• The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 
new assets: 
 Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?Answer: [NA]Justification: [NA] Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file. 
crowdsourcing and research with human subjects: 
 Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?Answer: [NA]Justification: [NA] Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector. 
institutional review board (irb) approvals or equivalent for research with human subjects: 
 Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?Answer: [NA]Justification: [NA] Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.• Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: EQUATIONProof. For any x, y ∈ X∥F (x) -Ψ v (x)∥ ≤ ∥F (x) -Φ v (x)∥ + ∥Φ v (x) -Ψ v (x)∥ (9) ≤ L 2 ∥x -v∥ 2 + ∥(∇F (v) -J(v))[x -v]∥ ≤ L 2 ∥x -v∥ 2 + δ∥x -v∥. □ Lemma 3.1 Let Assumptions (1.1), (1.2), (2.1) hold. Then for any x, v ∈ X VI (11) is monotone1 2 ∇Ω v (x) + ∇Ω v (x) T ⪰ 4L 1 ∥x -v∥I d×d + 5L 1 (x-v)(x-v) T ∥x-v∥ + (η -1)δI d×d .Proof. For all x, v ∈ X1 2 ∇Ω v (x) + ∇Ω v (x) T (10) = 1 2 J(v) + J(v) T + ηδI d×d + 5L∥x -v∥I d×d + 5L (x -v)(x -v) T ∥x -v∥ (29) ⪰ 1 2 ∇F (x) + ∇F (x) T + 4L∥x -v∥I d×d + 5L (x-v)(x-v) T ∥x-v∥ + (η -1)δI d×d (3) ⪰ 4L∥x -v∥I d×d + 5L (x-v)(x-v) T ∥x-v∥ + (η -1)δI d×d . □ B Proofs of Theorems 3.2, 3.3 B.1 Proof of Theorem 3.2In this section, we let L := L 1 . To show the convergence of Algorithm 1, we define the following Lyapunov functionEQUATION)Lemma B.1 Let Assumption 1.2, 2.1 hold. Then, for every integer T ≥ 1, we haveT k=1 λ k ⟨F (x k ), x k -x⟩ ≤ E 0 -E T + ⟨s T , x -x 0 ⟩ -1 8 T k=1 ∥x k -v k ∥ 2 , for all x ∈ X .Proof. By the definition of Lyapunov function (30) and Step 2 of Algorithm 1, we haveEQUATION)By the update formula for v k+1 , we get⟨x -v k+1 , s k -v k+1 + x 0 ⟩ ≤ 0, for all x ∈ X .Letting x = v k+2 in this inequality and using ⟨a, b⟩ = 1 2 (∥a + b∥ 2 -∥a∥ 2 -∥b∥ 2 ), we have⟨s k , v k+2 -v k+1 ⟩ ≤ ⟨v k+1 -x 0 , v k+2 -v k+1 ⟩ = 1 2 ∥v k+2 -x 0 ∥ 2 -∥v k+1 -x 0 ∥ 2 -∥v k+2 -v k+1 ∥ 2 .(32) Plugging Eq. ( 32) into Eq. ( 31) and using Step 5 of Algorithm 1, we obtain:E k+1 -E k (32) ≤ ⟨s k+1 -s k , v k+1 -x 0 ⟩ + ⟨s k+1 -s k , v k+2 -v k+1 ⟩ -1 2 ∥v k+2 -v k+1 ∥ 2 = ⟨s k+1 -s k , v k+2 -x 0 ⟩ -1 2 ∥v k+2 -v k+1 ∥ 2 ≤ λ k+1 ⟨F (x k+1 ), x 0 -v k+2 ⟩ -1 2 ∥v k+2 -v k+1 ∥ 2 = λ k+1 ⟨F (x k+1 ), x 0 -x⟩ + λ k+1 ⟨F (x k+1 ), x -x k+1 ⟩ + λ k+1 ⟨F (x k+1 ), x k+1 -v k+2 ⟩ -1 2 ∥v k+2 -v k+1 ∥ 2, for any x ∈ X . Summing up this inequality over k = 0, 1, . . . , T -1 and changing the counter k + 1 to k yields thatEQUATIONUsing the update formula for s k+1 and letting s 0 = 0 d ∈ R d , we haveI = T k=1 ⟨λ k F (x k ), x 0 -x⟩ = T k=1 ⟨s k-1 -s k , x 0 -x⟩ = ⟨s 0 -s T , x 0 -x⟩ = ⟨s T , x -x 0 ⟩. (34)Since x k+1 ∈ X satisfies (12), we haveEQUATIONwhere Ω η v (x) : R d → R d is defined in (10) . Letting x = v k+1 in (35) , we haveEQUATION)Then,EQUATION) ≤ L 2 ∥x k -v k ∥ 2 ∥x k -v k+1 ∥ + δ∥x k -v k ∥∥x k -v k+1 ∥ + L 2 ∥x k -v k ∥ 3 + δ∥x k -v k ∥ 2 -5L∥x k -v k ∥⟨x k -v k , x k -v k+1 ⟩ -ηδ⟨x k -v k , x k -v k+1 ⟩ Next, using ⟨x k -v k , x k -v k+1 ⟩ ≥ ∥x k -v k ∥ 2 -∥x k -v k ∥∥v k -v k+1 ∥ and ∥x k -v k+1 ∥ ≤ ∥x k -v k ∥ + ∥v k -v k+1 ∥, we get ⟨F (x k ), x k -v k+1 ⟩ ≤ L 2 ∥x k -v k ∥ 3 + L 2 ∥x k -v k ∥ 2 ∥v k -v k+1 ∥ + δ∥x k -v k ∥ 2 + δ∥x k -v k ∥∥v k -v k+1 ∥ + L 2 ∥x k -v k ∥ 3 + δ∥x k -v k ∥ 2 -5L∥x k -v k ∥ 3 + 5L∥x k -v k ∥ 2 ∥v k -v k+1 ∥ -ηδ∥x k -v k ∥ 2 +ηδ∥x k -v k ∥∥v k -v k+1 ∥ = 11L 2 ∥x k -v k ∥ 2 ∥v k -v k+1 ∥ -4L∥x k -v k ∥ 3 + (η + 1)δ∥x k -v k ∥∥v k -v k+1 ∥ -(η -2)δ∥x k -v k ∥ 2 Next, II ≤ T k=1 11λ k L 2 ∥x k -v k ∥ 2 ∥v k -v k+1 ∥ -4λ k L∥x k -v k ∥ 3 +(η + 1)δλ k ∥x k -v k ∥∥v k -v k+1 ∥ -(η -2)δλ k ∥x k -v k ∥ 2 -1 2 ∥v k -v k+1 ∥ 2 ≤ T k=1 1 2 ∥x k -v k ∥∥v k -v k+1 ∥ -1 4 ∥x k -v k ∥ 2 -1 2 ∥v k -v k+1 ∥ 2 ,where the last inequality is due to the following choice of η = 10 and λ :EQUATION)Plugging ( 34) and ( 37) into (33) yields thatT k=1 λ k ⟨F (x k ), x k -x⟩ ≤ E 0 -E T + ⟨s T , x -x 0 ⟩ -1 8 T k=1 ∥x k -v k ∥ 2 .Lemma B.2 Let Assumptions 1.2, 1.3, 2.1 hold and let x ∈ X . For every integer T ≥ 1, we haveEQUATIONwhere x * ∈ X denotes the weak solution to the VI.Proof. For any x ∈ X , we haveE 0 -E T + ⟨s T , x -x 0 ⟩ = E 0 -max v∈X ⟨s T , v -x 0 ⟩ -1 2 ∥v -x 0 ∥ 2 + ⟨s T , x -x 0 ⟩.Since s 0 = 0 d , we have E 0 = 0 andE 0 -E T + ⟨s T , x -x 0 ⟩ ≤ -⟨s T , x -x 0 ⟩ -1 2 ∥x -x 0 ∥ 2 + ⟨s T , x -x 0 ⟩ = 1 2 ∥x -x 0 ∥ 2 . This together with Lemma B.1 yields that T k=1 λ k ⟨F (x k ), x k -x⟩ + 1 8 T k=1 ∥x k -v k ∥ 2 ≤ 1 2 ∥x -x 0 ∥ 2 , for all x ∈ X ,which implies the first inequality. Since the VI satisfies the Minty condition, there exists x * ∈ X such that ⟨F (x k ), x k -x * ⟩ ≥ 0 for all k ≥ 1. Letting x = x * in the above inequality yields the second inequality. □ Lemma B.3 Let Assumptions 1.2, 1.3, 2.1 hold. For every integer T ≥ 1, we haveEQUATIONwhere x * ∈ X denotes the weak solution to the VI.Proof. Without loss of generality, we assume that x 0 ̸ = x * . We haveT k=1 (λ k ) -2 ( 1 32 ) 2 ≤ T k=1 (λ k ) -2 λ k L 2 ∥x k -v k ∥ + δ 2 = T k=1 L 2 ∥x k -v k ∥ + δ 2 ≤ T k=1 L 2 2 ∥x k -v k ∥ 2 + 2T δ 2 Lemma B.2 ≤ 2L 2 ∥x * -x 0 ∥ 2 + 2T δ 2 .By the Hölder inequality, we haveT k=1 1 = T k=1 (λ k ) -2 1/3 (λ k ) 2/3 ≤ T k=1 (λ k ) -2 1/3 T k=1 λ k 2/3.Putting these pieces together yields thatT ≤ 32 2/3 (2L 2 ∥x * -x 0 ∥ 2 + 2δ 2 T ) 1 3 T k=1 λ k 2/3, Plugging this into the above inequality yields that1 T k=1 λ k 2 ≤ 2048L 2 ∥x * -x 0 ∥ 2 T 3 + 2048δ 2 T 2 □ Theorem 3.2 Let Assumptions 1.2, 1.1, 2.1.Then, after T ≥ 1 iterations of VIJI with parameters β = δ, η = 10, opt = 0, we get the following boundGAP(x T ) = sup x∈X ⟨F (x), xT -x⟩ ≤ 16 √ 2LD 3 T 3/2 + 16 √ 2δD 2 T .Proof. Letting x ∈ X , we derive from the monotonicity of F and the definition of xT (i.e., opt = 0) that⟨F (x), xT -x⟩ = 1 T k=1 λ k T k=1 λ k ⟨F (x), x k -x⟩ .Combining this inequality with the first inequality in Lemma B.2 yields that⟨F (x), xT -x⟩ ≤ ∥x-x0∥ 2 2( T k=1 λ k ), for all x ∈ X . Since x 0 ∈ X , we have ∥x -x 0 ∥ ≤ D and hence⟨F (x), xT -x⟩ ≤ D 2 2( T k=1 λ k ) , for all x ∈ X .Then, we combine Lemma B.3 and the fact that ∥x * -x 0 ∥ ≤ D to obtain that⟨F (x), xT -x⟩ ≤ D 2 2 2048L 2 D 2 T 3 + 2048δ 2 T 2 ≤ 16 √ 2LD 3 T 3/2 + 16 √ 2δD 2 T , for all x ∈ X .By the definition of a gap function. ( 4), we haveEQUATION) □ B.2 Proof of Theorem 3.3We directly follow the steps of the proof of Theorem 3.2. Lemmas B.1, B.2 remain the same.Because of the choice of β k+1 = L1 2 ∥x k+1 -v k+1 ∥ adaptive strategy for λ k+1 looks as follows:EQUATIONwhere x * ∈ X denotes the weak solution to the VI.Proof. Without loss of generality, we assume that x 0 ̸ = x * . We haveT k=1 (λ k ) -2 ( 1 32 ) 2 ≤ T k=1 (λ k ) -2 (λ k (L∥x k -v k ∥)) 2 ≤ T k=1 L 2 ∥x k -v k ∥ 2 Lemma B.2 ≤ 4L 2 ∥x * -x 0 ∥ 2 .By the Hölder inequality, we haveT k=1 1 = T k=1 (λ k ) -2 1/3 (λ k ) 2/3 ≤ T k=1 (λ k ) -2 1/3 T k=1 λ k 2/3.Putting these pieces together yields thatT ≤ 32 2/3 (4L 2 ∥x * -x 0 ∥ 2 ) 1 3 T k=1 λ k 2/3.Plugging this into the above inequality yields that 1T k=1 λ k ≤ 64L∥x * -x 0 ∥ T 3/2 .Then, by following the rest of the proof of Theorem 3.2, we get Theorem 3.3 Let Assumptions 1.2, 1.1 hold. Let {x k , v k } be iterates generated by Algorithm 1 and∥(∇F (v k ) -J(v k ))[z k -v k ]∥ ≤ δ k ∥z k -v k ∥, δ k ≤ L1 2 ∥x k -v k ∥. Then, after T ≥ 1 iterations of Algorithm 1 with parameters β k = L1 2 ∥x k -v k ∥, η = 10, opt = 0, we get the following boundGAP(x T ) = sup x∈X ⟨F (x), xT -x⟩ ≤ 32LD 3 T 3/2 .C Proofs of Theorems 3.4, 3.5In this section, we let L := L 1 .Theorem 3.4 Let Assumptions 1.2, 1.3, 2.1 hold. Then after T ≥ 1 iterations of Algorithm 1 with parameters β k = δ, η = 10, opt = 2 we get the following boundEQUATIONProof. Since we consider opt = 2, then RES(x) = RES(x k T ) = sup x∈X ⟨F (x k T ), x k T -x⟩ .From (36) we getEQUATION≤ ∥FEQUATION)Next, from triangle inequality we have∥F (x k ) -Ψ v k (x k )∥ = F (x k ) -Ω η v k (x k ) + ηδ(x k -v k ) + 5L 2 ∥x k -v k ∥(x k -v k ) ≥ ∥F (x k ) -Ω η v k (x k )∥ -ηδ∥x k -v k ∥ - 5L 2 ∥x k -v k ∥ 2 .From this and (29) we get∥F (x k ) -Ω η v k (x k )∥ ≤ L 2 ∥x k -v k ∥ 2 + δ∥x k -v k ∥ + ηδ∥x k -v k ∥ + 5L 2 ∥x k -v k ∥ 2 = 3L∥x k -v k ∥ 2 + δ(η + 1)∥x k -v k ∥.Now, we can return to (43) :EQUATION) ≤ 3L∥x k -v k ∥ 2 ∥x k -x∥ + δ(η + 1)∥x k -v k ∥∥x k -x∥ + L 2 ∥x k -v k ∥ 3 + δ∥x k -v k ∥ 2 = L∥x k -v k ∥ 2 3∥x k -x∥ + 1 2 ∥x k -v k ∥ + δ∥x k -v k ∥ ((η + 1)∥x k -x∥ + ∥x k -v k ∥) .Since D := max x,y∈X ∥x -y∥,EQUATIONNext, from second inequality from (38) and from definition of x k T in Algorithm 1 we obtainEQUATIONSince (44) holds for any x ∈ X , we get final resultEQUATION) ≤ 14LD∥x * -x 0 ∥ 2 T + 2δ(η + 2)D∥x * -x 0 ∥ √ T ≤ 14LD 3 T + 2δ(η + 2)D 2 √ T = 14LD 3 T + 24D 2 √ T .□ Theorem 3.5 Let Assumptions 1.2, 1.3 hold. Let {x k , v k } be iterates generated by Algorithm 1 that satisfy (14) . Then after T ≥ 1 iterations of Algorithm 1 with parametersEQUATION)The proof of this theorem repeats the proof of Theorem 3.4 with δ ≤ L 2 ∥x k -v k ∥.Theorem 4.3 Let some first-order method M satisfy Assumption 4.1 and have access only δ-inexact first-order oracle 15. Assume the method M ensures for any L 0 -zero-order smooth and L 1 -first-order smooth monotone operator F the following convergence rateGAP(x) ≤ O(1) max δD 2 Ξ1(T ) ; L1D 3 Ξ2(T ) .Then for all T ≥ 1 we haveΞ 1 (T ) ≤ T, Ξ 2 (T ) ≤ T 3/2 .Proof. We prove this Theorem by contradiction. Assume the existence of a method M that satisfies the conditions of Theorem 4.3 and achieves faster rate in one of the terms (18) . First, suppose Ξ 1 (T ) > T . Consider the first-order lower bound from [89] , which is established using a quadratic min-max problem as the worst-case function. In this scenario, the operator has a 0-Lipschitz continuous Jacobian. Applying first-order method M to this lower bound and using an inexact Jacobian J(x) = L 0 I d×d , yields the rate O L0D 2 Ξ1(T ) , Ξ 1 (T ) > T , which is faster than the lower bound Ω L0D 2 T , contradicting our assumption.Secondly, let us assume Ξ 2 (T ) > T 3/2 . The lower bound for exact second-order methods is [69] . By taking exact Jacobian in the method M (δ = 0), we place M in the class of exact second-order methods. Consequently, we obtain a contradiction with the lower bound. □ E Proof of Theorem 5.1Ω L1D 3 T 3/2Proof. The proof is common for both formulas (20) and (21) , where α = 1 for classical L-Broyden approximation and α = m + 1 for Damped L-Broyden approximation. First, from L 0 -zero-order smoothness, get∥∇F (x) -J x ∥ op ≤ ∥∇F (x)∥ op + ∥J x ∥ op ≤ L 0 + ∥J x ∥ opNow, we upper-bound ∥J x ∥ op = ∥J m ∥ op by induction:J i+1 op ≤ J i + (yi-J i si)s ⊤ i αs ⊤ i si op ≤ J i I - sis ⊤ i αs ⊤ i si + yis ⊤ i αs ⊤ i si op ≤ J i I - sis ⊤ i αs ⊤ i si op + yis ⊤ i αs ⊤ i si op ≤ ∥J i ∥ op I - sis ⊤ i αs ⊤ i si op + ∥yi∥∥si∥ αs ⊤ i si ≤ ∥J i ∥ op + L0 α ,where the last inequality is coming from L 0 -zero-order smoothness for operator difference or JVP. By summing up the previous inequality for i in 0, . . . , m -1, we get ∥J m ∥ op ≤ ∥J 0 ∥ op + mL0 α . Finally, we prove the result of Theorem 5.1. □ F Proof of Theorem 6.2In this section, we let L := L 1 .Theorem 6. EQUATION) ≤ 2048L 2 ∥x 0 -x * ∥ 2 T 3 + 2048δ 2 T 2 1 2 • 1 2 ∥x 0 -x * ∥ 2 ≤ 2 max 2048L 2 ∥x 0 -x * ∥ 2 T 3 i , 2048δ 2 T 2 i 1 2 • 1 2 ∥x 0 -x * ∥ 2 .Denote R := D, R i = R 2 i , i ≥ 1. Now we run Algorithm 1 in cycle for i = 1, ..., n and restart it every time its distance to the solution becomes at least twice less than R i-1 . Thus, let T i be number of iterations we run Algorithm 1 inside cycle of Algorithm 2. In other words, letT i be such that xTi-1 -x * ≤ Ri-1 2where xTi+1 is the point, where we restart Algorithm 1. Then the number of iterations before the i-th restart isµ∥x Ti -x * ∥ 2 ≤ 2 max 2048L 2 R 2 i-1 T 3 i , 2048δ 2 T 2 i 1 2 • 1 2 R 2 i-1 ≤ µR 2 i-1 4 ⇔ (47) ⇔ max 2048L 2 R 2 i-1 T 3 i , 2048δ 2 T 2 i 1 2 ≤ µ 2 √ 2 .Deriving T i for each case under max, we get   T i ≥ 2 14 3 L 2 3 R 2 3 i-1 µ 2 3T i ≥ 2 7 δ µ . Thus,EQUATIONNow we calculate the total number of restarts to reach ∥x Tn -x * ∥ ≤ ε. From (47) we can get that∥x Tn -x * ∥ ≤ 1 µ 2 max 2048L 2 R 2 n-1 T 3 n , 2048δ 2 T 2 n 1 2 • 1 2 R 2 n-1 (48) ≤ R n-1 √ 2µ 2 max µ 2 8 , µ 2 8 1 4 = R n-1 2 = R • 2 -(n-1)-1 ≤ ε.Deriving n from last inequality, we getEQUATIONNow we provide auxiliary estimation, that we will need next:EQUATIONFinally, we can get the total number of iterations of Algorithm 1 inside Algorithm 2:EQUATION≤ 2 Input: initial point x 0 ∈ X , parameters L 1 , η, sequence {δ i } p-1 i=1, and opt ∈ {0, 1, 2}.Initialization: set s 0 = 0 ∈ R d . for k = 0, 1, 2, . . . , T do Compute v k+1 = argmax v∈X {⟨s k , v -x 0 ⟩ -1 2 ∥v -x 0 ∥ 2 }. Compute x k+1 ∈ X such that condition (53) holds true. Compute λ k+1 such that 1 4(5p-2) ≤ λ k Lp-1 p! ∥x k -v k ∥ p-1 + p-1 i=1 δi i! ∥x k -v k ∥ i-1 ≤ 1 2(5p+1) . Compute s k+1 = s k -λ k+1 F (x k+1 ). Output: x =    xT =In this section, we provide more details on the generalization of Algorithm 1 with high-order derivatives. We provide the pseudocode of the resulting method in Algorithm 3. To show the convergence of Algorithm 3, we use the Lyapunov function (30) .Define the (p -1)-th order approximations of the FEQUATIONEQUATIONOn each step, our method solves the following subproblem:EQUATIONAdditionally, we will need an auxiliary result from [52] , based on Assumption 4.2. The authors show, that this Assumption allows to control the quality of approximation of operator F by its high-order Taylor polynomial:EQUATIONFirst of all, we provide high-order generalizations of auxiliary lemmas for second-order case from Section A.Lemma G.1 Let Assumptions 7.1 and 4.2 with i = p -1 hold. Then, for any x, v ∈ XEQUATION)Proof.∥F (v) -Ψ p,v (x)∥ ≤ ∥F (v) -Φ p,v (x)∥ + ∥Φ p,v (x) -Ψ p,v (x)∥ (54) ≤ Lp-1 p! ∥x -v∥ p + p-1 i-1 1 i! ∥(∇ i F (v) -G i (v))[x -v] i-1 ∥∥x -v∥ (27) ≤ Lp-1 p! ∥x -v∥ p + p-1 i=1 δi i! ∥x -v∥ i .Lemma G.2 Let Assumptions 1.1, 7.1 and 4.2 with i = p -1 hold. Then for any x, v k+1 ∈ X VI (11) is relatively strongly monotone ifη i ≥ p 1 2 ∇Ω p,v (x) + ∇Ω p,v (x) T 4Lp-1 (p-1)! ∥x -v∥ p-1 I n×n + ∥x -v∥ p-3 (x -v)(x -v) T + p-1 i=1 δi i! ∥x -v∥ i-3 (η i -i)∥x -v∥ 2 I n×n + η i (i -1)(x -v)(x -v) T .Proof.1 2 ∇Ω p,v (x) + ∇Ω p,v (x) T (28) = 1 2 p-1 i=1 1 (i-1)! G i (v)[x -v] i-1 + G i (v)[x -v] i-1 T + p-1 i=1 ηiδi i! ∥x -v∥ i-1 I n×n + (i -1)∥x -v∥ i-3 (x -v)(x -v) T + 5Lp-1 (p-1)! ∥x -v∥ p-1 I n×n (55) ⪰ 1 2 ∇F (x) + ∇F (x) T - Lp-1 (p-1)! ∥x -v∥ p-1 I n×n - p-1 i=1 1 (i-1)! δ i ∥x -v∥ i-1 I n×n + p-1 i-1 ηiδi i! ∥x -v∥ i-1 I n×n + (i -1)∥x -v∥ i-3 (x -v)(x -v) T + 5Lp-1 (p-1)! ∥x -v∥ p-1 I n×n + (p -1)∥x -v∥ p-3 (x -v)(x -v) T = 1 2 ∇F (x) + ∇F (x) T + 4Lp-1 (p-1)! ∥x -v∥ p-1 I n×n + 5Lp-1 (p-2)! ∥x -v∥ p-3 (x -v)(x -v) T + p-1 i=1 δi i! (η i -i)∥x -v∥ i-1 I n×n + p-1 i=1 ηiδi(i-1) i! ∥x -v∥ i-3 (x -v)(x -v) T .From monotonicity of F we know that 1 2 (F (x) -F (x) T ) ⪰ 0. Thus,1 2 ∇Ω p,v (x) + ∇Ω p,v (x) T ⪰ 4Lp-1 (p-1)! ∥x -v∥ p-1 I n×n + 5Lp-1 (p-2)! ∥x -v∥ p-3 (x -v)(x -v) T + p-1 i=1 δi i! (η i -i)∥x -v∥ i-1 I n×n + p-1 i=1 ηiδi(i-1) i! ∥x -v∥ i-3 (x -v)(x -v) T .By rearranging the terms we get1 2 ∇Ω p,v (x) + ∇Ω p,v (x) T 4Lp-1 (p-1)! ∥x -v∥ p-1 I n×n + ∥x -v∥ p-3 (x -v)(x -v) T + p-1 i=1 δi i! ∥x -v∥ i-3 (η i -i)∥x -v∥ 2 I n×n + η i (i -1)(x -v)(x -v) T .Thus, if η i ≥ p, we get that 1 2 ∇Ω p,v k+1 (x) + ∇Ω p,v k+1 (x) T is relatively strongly monotone. □In this subsection we provide theoretical results directly connected to convergence rate of Algorithm 3. Firstly, we need to introduce additional technical lemmas, that generalize corresponding lemmas in Section B.Lemma G.3 Let Assumption 1.1 hold and η i = 5p. Then, for every T ≥ 1, we haveEQUATION)Proof. Since the first steps in the proof of this Lemma are the same as in Lemma B.1, we start our reasoning from (33):T k=1 λ k ⟨F (x k ), x k -x⟩ ≤ E 0 -E T + T k=1 λ k ⟨F (x k ), x 0 -x⟩ I + T k=1 λ k ⟨F (x k ), x k -v k+1 ⟩ -1 2 ∥v k -v k+1 ∥ 2 II .(57) Using the update formula for s k+1 and letting s 0 = 0 d ∈ R d , we haveEQUATION)Now consider II. Using (28) we get⟨F (x k ), x k -v k+1 ⟩ (28) = ⟨F (x k ) -Ψ p,v k (x k ), x k -v k+1 ⟩ + ⟨Ω p,v k (x k ), x k -v k+1 ⟩ - p-1 i=1 ηiδi i! ∥x k -v k ∥ i-1 ⟨x k -v k , x k -v k+1 ⟩ - 5Lp-1 (p-1)! ∥x k -v k ∥ p-1 ⟨x k -v k , x k -v k+1 ⟩ (55) ≤ Lp-1 p! ∥x k -v k ∥ p ∥x k -v k+1 ∥ + p-1 i=1 1 i! δ i ∥x k -v k ∥ i ∥x k -v k+1 ∥ + ⟨Ω p,v k (x k ), x k -v k+1 ⟩ - p-1 i=1 ηiδi i! ∥x k -v k ∥ i-1 ⟨x k -v k , x k -v k+1 ⟩ - 5Lp-1 (p-1)! ∥x k -v k ∥ p-1 ⟨x k -v k , x k -v k+1 ⟩ . Next, using ⟨x k -v k , x k -v k+1 ⟩ ≥ ∥x k -v k ∥ 2 -∥x k -v k ∥∥v k -v k+1 ∥ and ∥x k -v k+1 ∥ ≤ ∥x k -v k ∥ + ∥v k -v k+1 ∥, we get ⟨F (x k ), x k -v k+1 ⟩ ≤ Lp-1 p! ∥x k -v k ∥ p+1 + ∥x k -v k ∥ p ∥v k -v k+1 ∥ + p-1 i=1 1 i! δ i ∥x k -v k ∥ i+1 ∥x k -v k ∥ i ∥v k -v k+1 ∥ - p-1 i=1 ηiδi i! ∥x k -v k ∥ i+1 -∥x k -v k ∥ i ∥v k -v k+1 ∥ - 5Lp-1 (p-1)! ∥x k -v k ∥ p+1 -∥x k -v k ∥ p ∥v k -v k+1 ∥ + ⟨Ω p,v k (x k ), x k -v k+1 ⟩ .From definition of the subproblem (53) we have that⟨Ω p,v k (x k ), x k -v k+1 ⟩ ≤ sup x∈X ⟨Ω p,v k (x k ), x k -x⟩ ≤ Lp-1 p! ∥x k -v k ∥ p+1 + p-1 i=1 δi i! ∥x k -v k ∥ i+1 .Thus,⟨F (x k ), x k -v k+1 ⟩ ≤ Lp-1(1+5p) p! ∥x k -v k ∥ p ∥v k -v k+1 ∥ - Lp-1(5p-2) p! ∥x k -v k ∥ p+1 + p-1 i=1 δi(1+ηi) i! ∥x k -v k ∥ i ∥v k -v k+1 ∥ - p-1 i=1 δi(ηi-2) i! ∥x k -v k ∥ i+1From this we getII = T k=1 λ k ⟨F (x k ), x k -v k+1 ⟩ -1 2 ∥v k -v k+1 ∥ 2 ≤ T k-1 λ k Lp-1(1+5p) p! ∥x k -v k ∥ p ∥v k -v k+1 ∥ - λ k Lp-1(5p-2) p! ∥x k -v k ∥ p+1 +λ k p-1 i=1 δi(1+ηi) i! ∥x k -v k ∥ i ∥v k -v k+1 ∥ -λ k p-1 i=1 δi(ηi-2) i! ∥x k -v k ∥ i+1 -1 2 ∥v k -v k+1 ∥ 2 ηi=5p = T k=1 λ k Lp-1 p! ∥x k -v k ∥ p-1 + p-1 i=1 δi i! ∥x k -v k ∥ i-1 (5p + 1)∥x k -v k ∥∥v k -v k+1 ∥ -λ k Lp-1 p! ∥x k -v k ∥ p-1 + p-1 i=1 δi i! ∥x k -v k ∥ i-1 (5p -2)∥x k -v k ∥ 2 -1 2 ∥v k -v k+1 ∥ 2 .Now, if we choose λ k in a such way thatEQUATIONwe get the followingII (59) ≤ T k=1 1 2 ∥x k -v k ∥∥v k -v k+1 ∥ -1 4 ∥x k -v k ∥ 2 -1 2 ∥v k -v k+1 ∥ 2 ≤ T k=1 max γ 1 2 γ∥x k -v k ∥ -1 2 γ 2 -1 4 ∥x k -v k ∥ 2 ≤ T k=1 1 4 ∥x k -v k ∥ 2 -1 8 ∥x k -v k ∥ 2 -1 4 ∥x k -v k ∥ 2 = -1 8 T k=1 ∥x k -v k ∥ 2 .Finally, combining estimations of I and II with (57), we getT t=1 λ k ⟨F (x k ), x k -x⟩ ≤ E 0 -E T + ⟨s T , x -x 0 ⟩ -1 8 T k=1 ∥x k -v k ∥ 2 .Lemma G.4 Let Assumptions 4.2, 1.3 hold. For every T ≥ 1 we haveEQUATION)Proof. From Hölder inequalityEQUATIONConsider denominator:T k=1 λ -2 p-1 k (59) ≤ T k=1 λ -2 p-1 k (4(5p -2)) 2 p-1 λ 2 p-1 k L p-1 p! ∥x k -v k ∥ p-1 + p-1 i=1 δ i i! ∥x k -v K ∥ i-1 2 p-1 = T k=1 (20p -8) 2 p-1 L p-1 p! ∥x k -v k ∥ p-1 + p-1 i=1 δ i i! ∥x k -v k ∥ i-1 2 p-1 .For p ≥ 2 we have2 p-1 = 2, p = 2, 2 p-1 ≤ 1, p ≥ 3.Consider some nonnegative sequence {a i |a i ≥ 0, ∀i ∈ 1, n}. For p = 2 from Jensen inequality we know that (n i=1 a i ) 2 ≤ n i=1 na 2 i .From Lemma 7 of [65] we know that ∀q ∈ [0, 1], x, y > 0 → (x + y) q ≤ x q + y q . Thus, for p ≥ 3 we can come to the same conclusion:EQUATION)Using this inequality, we getEQUATION≤ p(20p -8)2 p-1 T k=1 L p-1 p! 2 p-1 ∥x k -v k ∥ 2 + p(20p -8) 2 p-1 T k=1 p-1 i=1 δ i i! 2 p-1 ∥x k -v k ∥ 2(i-1) p-1 (38) ≤ 4p(20p -8) 2 p-1 L p-1 p! 2 p-1 ∥x * -x 0 ∥ 2 + p(20p -8) 2 p-1 T k=1 p-1 i=1 δ i i! 2 p-1 ∥x k -v k ∥ 2(i-1) p-1Consider the second factor in this inequality. If we denotea k = ∥x k -v k ∥ 2(i-1) p-1 , b k = 1, c = p-1 i-1 , d = p-1p-i , then we can use Hölder inequality:n k=1 |a k b k | ≤ n k=1 |a k | c 1 c n i=1 |b k | d 1 d , 1 c + 1 d = 1.From this we getT k=1 ∥x k -v k ∥ 2(i-1) p-1 ≤ T k=1 ∥x k -v k ∥ 2 i-1 p-1 T k=1 1 p-i p-1 (38) ≤ 4∥x * -x 0 ∥ 2(i-1) p-1 T p-i p-1 .Thus,T k=1 λ -2 p-1 k ≤ 4p(20p -8) 2 p-1 L p-1 p! 2 p-1 ∥x * -x 0 ∥ 2 + 4p(20p -8) 2 p-1 p-1 i=1 δ i i! 2 p-1 ∥x * -x 0 ∥ 2(i-1) p-1 T p-i p-1 Consider (a + b) p-12 , p ≥ 2. For p = 2 from Lemma 7 in [65] we get (a + b)EQUATION. In other words,EQUATION)Thus,T k=1 λ -2 p-1 k p-1 2 (62) ≤ 2 p 2   (4p) p-1 2 (20p -8) L p-1 p! ∥x * -x 0 ∥ p-1 + (4p) p-1 2 (20p -8) p-1 i=1 δ i i! 2 p-1 ∥x * -x 0 ∥ 2(i-1) p-1 T p-i p-1 p-1 2    ≤ 2 p p p-1 2 (20p -8) L p-1 p! ∥x * -x 0 ∥ p-1 + 2 p 2 p-1 i=1 δ i i! ∥x * -x 0 ∥ i-1 T p-i 2 = 2 p p p-1 2 (20p -8) L p-1 p! ∥x * -x 0 ∥ p-1 + 2 3p 2 p p-1 2 (20p -8) p-1 i=1 δ i i! ∥x * -x 0 ∥ i-1 T p-i 2 .Now, we can get the final result1 T k=1 λ k ≤ T k=1 λ -2 p-1 k p-1 2 T p+1 2 ≤ 2 p p p-1 2 (20p -8) Lp-1 p! ∥x * -x 0 ∥ p-1 T p+1 2 + 2 3p 2 p p-1 2 (20p -8) p-1 i=1 δ i ∥x * -x 0 ∥ i-1 i!T i+1 2 □Finally, we provide the convergence rate of Algorithm 3 in monotone case.Theorem 7.2 Let Assumptions 1.1, 4.2 with i = p -1, and 7.1 hold. Then, after T ≥ 1 iterations of VIHI with parameters η i = 5p, opt = 0, we get the following boundEQUATION.Proof. Consider ∀x ∈ X , opt = 0.EQUATION≤ 1T k=1 λ k 1 2 ∥x -x 0 ∥ 2 ≤ D 2 2 T i=1 λ k (60) ≤ 2 p-1 p p-1 2 (20p -8) L p-1 D p+1 p!T p+1 2 + 2 3p 2 -1 p p-1 2 (20p -8) p-1 i=1 δ i D i+1 i!T i+1 2 . Thus, GAP(x T ) = sup x∈X ⟨F (x), xT -x⟩ = O L p-1 D p+1 T p+1 2 + p-1 i=1 δ i D i+1 T i+1 2 . □ G.For r-rank QN approximation J r from ( 19), we can effectively compute A -1 τ F (x), whereEQUATIONτ )I by using the Woodbury matrix identity [96, 97] .A -1 τ F (x) = B + U ⊤ CV -1 F (x) = B -1 F (x) -B -1 U ⊤ (C -1 + V B -1 U ⊤ ) -1 V B -1 F (x).For J 0 = ιI, the identity could be simplified even more. Hence, B = ιI + (ηδ + 5L 2 τ )I = χI, where χ = ι + ηδ + 5L 2 τ . Then,A -1 τ F (x) = χI + U ⊤ CV -1 F (x) = 1 χ F (x) -1 χ U ⊤ (χC -1 + V U ⊤ ) -1 V F (x).We consider the cubic regularized bilinear min-max problem of the form:EQUATIONwhereρ > 0, b = [1, 0, . . . , 0] ∈ R d , and A ∈ R d×d such that A =          1 -1 0 • • • 0 0 1 -1 . . . . . . . . . . . . . . . . . . 0 0 • • • 0 1 -1 0 • • • 0 0 1 .         To reformulate it as variational inequality, we define F (x) = [∇ x f (x, y), -∇ y f (x, y)]. Following [71] , we plot the restricted primal-dual gap (65) , written in a closed form:gap(z, β) = ρ 6 ∥x∥ 3 + β∥Ax -b∥ + 2 3 2 ρ ∥A ⊤ y∥ 3 2 + b ⊤ y,where z = (x, y).Setup. All methods and experiments were performed using Python 3.11.5, PyTorch 2.1.2, numpy 1.24.3 on a 16-inch MacBook Pro 2023 with an Apple M2 Pro and 32GB memory.Parameters. For Figure 1 , the dimension of the problem d = 50. The regularizer ρ is set to ρ = 1e -3. The starting point (x 0 , y 0 ) = (0, 0) is all zeroes. We present results of last iteration of each method or opt=1 from Algorighm 1. The diameter β for the gap is set as β = 1. For QN methods, we set memory size m = r = 20. We perform a total of 100000 iterations for all methods except Perseus2 with 1000 iterations. We plot each 500 iterations for a better visualisation. In Figure 1 , the left plot refers to the gap decrease per iteration and the right plot refers to the gap decrease per JVP/operator computations. EG, Perseus1, VIQA computes 2 operators per iteration and Perseus2 computes d + 1 JVP/operators.We finetuned learning rate and presented the run with the best results. For EG1: lr = 0.5. For Perseus1: L 0 = 0.2334. For Perseus2: L 1 = 0.0001. For VIQA Broyden: L 1 = 0.001, δ = J 0 = 0.4. For VIQA Damped Broyden: L 1 = 0.001, δ = J 0 = 0.22. Note, that tuned L 1 = 0.001 is actually a theoretical smoothness constant for the problem (63) and can be chosen with such value without finetuning. Also, Perseus2 is working with a smaller constant than theoretical.In this section, we consider the problem of finding a global saddle point of the following min-max optimization problem:EQUATIONi.e., a tuple(x * , y * ) ∈ R m × R n such that f (x * , y) ≤ f (x * , y * ) ≤ f (x, y * ), for all x ∈ R m , y ∈ R n ,where the continuously differentiable objective function f is convex-concave: f (x, y) is convex in x for all y ∈ R n and concave in y for all x ∈ R m .Input: initial point z 0 ∈ X , parameters L 1 , δ, η, τ .Initialization: set s 0 = 0 ∈ R d . for k = 0, 1, 2, . . . , T do Set v k+1 = z 0 + s k+1 . Compute z k+1 ∈ R n+m such that condition (68) holds true. Compute λ k+1 such that 1 16 ≤ λ k+1 L 2 ∥z k+1 -v k+1 ∥ + δ ≤ 1 12 . Compute s k+1 = s k -λ k+1 F (z k+1 ). Output: zT = 1 T k=1 λ k T k=1 λ k z k . Assumption J.1 The function f (x, y) ∈ C 2 has L-Lipschitz-continious second-order derivative if ∥∇ 2 f (z) -∇ 2 (v)∥ ≤ L∥z -v∥, for all z, v ∈ R n+m .Following [71] , we define the restricted gap function to measure the optimality of point ẑ = (x, ŷ) gap(ẑ, β) = maxEQUATIONwhere β is sufficiently large such that ∥ẑ -z * ∥ ≤ β.Problem ( 64) is a special case of VI defined by the following operatorEQUATIONThe Jacobian of F is defined as followsEQUATIONThe following lemma [80, 71] provides properties of operator F .Lemma J.2 Let Assumption J.1 hold. Then 1. The operator F is monotone, i.e. satisfies Assumption 1.1.2. The operator F is L-smooth, i.e. satisfies Assumption 1.2.3. F (z * ) = 0 for any global saddle point z * ∈ R n+m of the function f .We assume that inexact approximation of Jacobian satisfies Assumption 2.1.Since we consider unconstrained minmax optimization, Step 2 of VIJI simplifies to v k+1 = z 0 + s k and the subproblem (11) changes tofind z k+1 ∈ R n+m such that Ω η v k+1 (z k+1 ) = 0.Usually, to find the solution of this subproblem it is necessary to run some addition subroutine. Following the work [71] , we introduce the following approximate condition:EQUATIONwhere τ ∈ (0, 1) is a tolerance parameter.The version of VIJI for unconstrained min-max problems is referred to as VIJI-MinMax and is detailed in Algorithm 4. This subproblem can solved by strategy proposed in [71] , resulting in O (n + m) ω ε -2/3 + (n + m) 2 ε -2/3 log log(1/ε) complexity, where ω ≈ 2.3728 is the matrix multiplication constant.The following theorem provides convergence rate for Algorithm 4.Theorem J.3 Let Assumptions J.1, 2.1 hold. Then, after T ≥ 1 iterations of VIJI-MinMax with parameters η = 5, we get the following boundgap(z T , β) ≤ 1152 √ 2L∥z0-z * ∥ 3 T 3/2 + 576 √ 2δ∥z0-z * ∥ 2 T, where β = 5∥z 0 -z * ∥, z * = (x * , y * ).Now, let us introduce additional assumption on δ, to obtain the convergence with optimal rate O(ε -2/3 ).Theorem J.4 Let Assumptions J.1 hold. Let∥(∇F (v k ) -J(v k ))[z k -v k ]∥ ≤ δ k ∥z k -v k ∥hold for iterates {z k , v k } generated by Algorithm 4, whereEQUATIONThen, after T ≥ 1 iterations of VIJI-MinMax with parameters η = 5, we get the following boundgap(z T , β) ≤ 1944 √ 2L∥z0-z * ∥ 3 T 3/2, where β = 5∥z 0 -z * ∥, z * = (x * , y * ).Note, that it is also possible to change adaptive strategy for λ k in this case to1 16 ≤ 3 2 L∥z k -v k ∥ ≤ 1 12(by introducing parameter β as in Algorithm 1) to eliminate the dependence on δ from the step size while achieving the same convergence rate.Finally, let us show, that we mathch the convergence of [71] under the same assumptions on Jacobian's inexactness and subproblem's solution.Assumption J.5 ([71] ) Let∥(∇F (v k ) -J(v k ))[z k -v k ]∥ ≤ δ k ∥z k -v k ∥, ∥J(v k )∥ ≤ κhold for iterates {z k , v k } generated by Algorithm 4, whereEQUATIONwhere τ 0 < L 4 .Next, we show, that (69) follows from (70) .Let us consider two cases.If ∥z k -v k ∥ ≤ 1, then from (70), we get δ k ≤ τ 0 ∥z k -v k ∥ ≤ L 2 ∥z k -v k ∥. Otherwise, if ∥z k -v k ∥ ≥ 1, we obtain from (68) ∥F (v k )∥ -∥J(v k )∥∥z k -v k ∥ -δ∥z k -v k ∥ -5L∥z k -v k ∥ ≤ ∥∇Ω v k (z k )∥ ≤ τ ∥F (v k )∥.Using our assumptions, we get(1 -τ )∥F (v k )∥ ≤ κ∥z k -v k ∥ + δ∥z k -v k ∥ + 5L∥z k -v k ∥ 2 ≤ (κ + δ + 5L)∥z k -v k ∥ ≤ (κ + 21 4 L)∥z k -v k ∥. Next, δ ≤ L(1-τ ) 4κ+6L ∥F (v k )∥ ≤ L∥z k -v k ∥.Thus, the assumptions of Theorem J.4 hold, and Algorithm 4 achieves O(ε -2/3 ) convergence.J.4 Proof of Theorem J.3Again, as in the proof of Theorem 3.2, we introduce the Lyapunov functionEQUATION)Note that the scalar product ⟨s k , v -z 0 ⟩ can be omitted (see the proof of [71, Theorems 3.1, 4.1]), which would also eliminate the somewhat redundant Step 2 of Algorithm 4. However, we chose to retain it, as it does not affect the method's performance. We retain Lemma J.6 Let Assumptions J.1, 2.1 hold. Then, for every integer T ≥ 1, we haveT k=1 λ k ⟨F (z k ), z k -z⟩ ≤ E 0 -E T + ⟨s T , z -z 0 ⟩ -1 8 T k=1 ∥z k -v k ∥ 2 , for all z ∈ R n+m .Proof. Following the proof of Lemma B.1, we arrive at (35) , where the proofs begin to slightly diverge due to the changes in the subproblem. At this juncture, we have:EQUATION≤ L 2 ∥z k -v k ∥ 2 ∥z k -v k+1 ∥ + δ∥z k -v k ∥∥z k -v k+1 ∥ + τ L 2 ∥z k -v k ∥ 2 ∥z k -v k+1 ∥ +τ δ∥z k -v k ∥∥z k -v k+1 ∥ -5L∥z k -v k ∥⟨z k -v k , z k -v k+1 ⟩ -ηδ⟨z k -v k , z k -v k+1 ⟩ Next, using ⟨z k -v k , z k -v k+1 ⟩ ≥ ∥z k -v k ∥ 2 -∥z k -v k ∥∥v k -v k+1 ∥ and ∥z k -v k+1 ∥ ≤ ∥z k -v k ∥ + ∥v k -v k+1 ∥, we get ⟨F (z k ), z k -v k+1 ⟩ ≤ (τ + 1) L 2 ∥z k -v k ∥ 3 + (τ + 1) L 2 ∥z k -v k ∥ 2 ∥v k -v k+1 ∥ + (τ + 1)δ∥z k -v k ∥ 2 +(τ + 1)δ∥z k -v k ∥∥v k -v k+1 ∥ -5L∥z k -v k ∥⟨z k -v k , z k -v k+1 ⟩ -ηδ⟨z k -v k , z k -v k+1 ⟩ τ <1 ≤ 6L∥z k -v k ∥ 2 ∥v k -v k+1 ∥ -4L∥z k -v k ∥ 3 + (η + 1) δ∥z k -v k ∥∥v k -v k+1 ∥ -(η -1) δ∥z k -v k ∥ 2 Next, II ≤ T k=1 1 2 ∥z k -v k ∥∥v k -v k+1 ∥ -1 4 ∥z k -v k ∥ 2 -1 2 ∥v k -v k+1 ∥ 2 ,where the last inequality is due to the following choice of η = 10 and λ :EQUATIONPlugging ( 73) into (72) yields thatT k=1 λ k ⟨F (z k ), z k -z⟩ ≤ E 0 -E T + ⟨s T , z -z 0 ⟩ -1 8 T k=1 ∥z k -v k ∥ 2 .Next, Lemma B.2 can be applied. Next Lemma is a counterpart of Lemma B.3 for the current choice of λ.Lemma J.7 Let Assumptions J.1, 2.1 hold. For every integer T ≥ 1, we haveEQUATIONwhere x * ∈ X denotes the weak solution to the VI.Proof. Without loss of generality, we assume that z 0 ̸ = z * . We haveT k=1 (λ k ) -2 16 -2 ≤ T k=1 (λ k ) -2 (λ k (L∥z k -v k ∥ + δ)) 2 = T k=1 (L∥z k -v k ∥ + δ) 2 ≤ T k=1 2L 2 ∥z k -v k ∥ 2 + 2T δ 2 Lemma B.2 ≤ 8L 2 ∥z * -z 0 ∥ 2 + 2T δ 2 .By the Hölder inequality, we haveT k=1 1 = T k=1 (λ k ) -2 1/3 (λ k ) 2/3 ≤ T k=1 (λ k ) -2 1/3 T k=1 λ k 2/3.Putting these pieces together yields thatT ≤ 16 2/3 (8L 2 ∥z * -z 0 ∥ 2 + 2δ 2 T ) 1 3 T k=1 λ k 2/3, Plugging this into the above inequality yields that1 T k=1 λ k 2 ≤ 2048L 2 ∥z * -z 0 ∥ 2 T 3 + 512δ 2 T 2Next, by Lemma J.6, we have 0Lem. J.2 ≤ T k=1 λ k ⟨F (z k ), z k -z⟩ + -1 8 T k=1 ∥z k -v k ∥ 2 ≤ E 0 -E T + ⟨s T , z -z 0 ⟩ (71) = 1 2 ∥v 0 -z 0 ∥ 2 - 1 2 ∥v T -z 0 ∥ 2 + ⟨v T -z 0 , v * -v 0 ⟩.By applying Young's ineqaulity and the fact that z 0 = v 0 , we get0 ≤ - 1 2 ∥v k -z 0 ∥ 2 + 1 4 ∥v k -z 0 ∥ 2 + ∥z * -z 0 ∥ 2 = - 1 4 ∥v k -z 0 ∥ 2 + ∥z * -z 0 ∥ 2 .Thus, we have ∥v k -z 0 ∥ ≤ 2∥z * -z 0 ∥. From Lemma B.2 we also have that ∥z k -v k ∥ ≤ 2∥z * -z 0 ∥. Thus,∥v k -z * ∥ ≤ ∥v k -z 0 ∥ + ∥z 0 -z * ∥ ≤ 3∥z 0 -z * ∥ ≤ β, ∥z k -z * ∥ ≤ ∥z k -v k ∥ + ∥z k -z * ∥ ≤ 5∥z 0 -z * ∥ = β. Lemma B.2 also implies T k=1 λ k (z k -z) ⊤ F (z k ) ≤ 1 2 ∥z 0 -z∥ 2 .By Proposition [71, Proposition 2.9], we havef (x T , y) -f (x, ỹT ) ≤ 1 T k=1 λ k T k=1 λ k (z k -z) ⊤ F (z k ) .Putting these pieces together yields f (x T , y) -f (x, ỹT ) ≤1 2( T k=1 λ k ) ∥z 0 -z∥ 2 .This together with Lemma J.7 yieldsf (x T , y) -f (x, ỹT ) ≤ 8 √ 2L∥z 0 -z * ∥∥z 0 -z∥ 2 T 3/2 + 4 √ 2δ∥z 0 -z∥ 2 T .Since ∥z k -z * ∥ ≤ β for all k ≥ 0, we have ∥z T -z * ∥ ≤ β. By the definition of the restricted gap function, we havegap(z T , β) ≤ 32 √ 2L∥z0-z * ∥(∥z0-z * ∥+β) 2 T 3/2 + 16 √ 2δ(∥z0-z * ∥+β) 2 T ≤ 1152 √ 2L∥z0-z * ∥ 3 T 3/2 + 576 √ 2δ∥z0-z * ∥ 2 T. Therefore, we conclude from the above inequality that there exists some T > 0 such that the output ẑ satisfies that gap(ẑ, β) ≤ ϵ. Lemma J.8 Let Assumption J.1. For every integer T ≥ 1, we have 1T k=1 λ k 2 ≤ 54L 2 ∥z * -z 0 ∥ 2 T 3 ,where x * ∈ X denotes the weak solution to the VI.Proof. Without loss of generality, we assume that z 0 ̸ = z * . We haveEQUATION≤ T k=1 9L 2 4 ∥z k -v k ∥ 2 Lemma B.2 ≤ 9L 2 4 ∥z * -z 0 ∥ 2 .By the Hölder inequality, we haveT k=1 1 = T k=1 (λ k ) -2 1/3 (λ k ) 2/3 ≤ T k=1 (λ k ) -2 1/3 T k=1 λ k 2/3.Putting these pieces together yields thatT ≤ 16 2/3 ( 9 4 L 2 ∥z * -z 0 ∥ 2 ) 1 3 T k=1 λ k 2/3, Plugging this into the above inequality yields that 1T k=1 λ k 2 ≤ 54L 2 ∥z * -z 0 ∥ 2 T 3 .Next, directly following the steps of proof of Theorem J.3, we get the following rategap(z T , β) ≤ 1944 √ 2L∥z0-z * ∥ 3 T 3/2 .NeurIPS Paper ChecklistJustification: Yes, the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction).For clarity, let us note that VI methods that use only the information about the operator itself are commonly referred to as first-order methods. Methods that also use information about the Jacobian of the operator (the first-order derivative) are known as second-order methods. This somewhat contradictory notation stems from applications to minimization problems, where the operator is the gradient (i.e., the first-order derivative), and the Jacobian is the Hessian (i.e., the second-order derivative).A formal definition is provided in the following section.https://github.com/OPTAMI/OPTAMIT k=1 λ k T k=1 λ k x k , if opt = 0, x T , else if opt = 1, x k T for k T = argmin 1≤k≤T ∥x k -v k ∥, else if opt = 


Decision: The authors present second-order methods for variational inequalities that can handle inexact Jacobian information, along with lower bounds that depend on the amount of inexactness. The reviewers all agree that the work is interesting and impactful, and thereby worthy of acceptance.

Reviews: review: paper_summary: paper_summary: This paper presents a new second-order algorithm for variational inequalities that is applicable to inexact Jacobian. The complete convergence analysis and a theoretical lower bound are provided in the paper.
strength_weakness: strength_weakness: The advantages of the proposed method over some related works are not stated clearly (see question 1).
questions: questions: 1. What is the advantage of the proposed method over [60] since the convergence rate is not improved? To my understanding, the assumption of Jacobian inexactness in [60] is not suitable for quasi-newton. But why "suitable to quasi-newton" can be claimed as an advantage? 
2. Will inexact Jacobian causes more iterations to converge than exact Jacobian, which eventually leads to a higher complexity than exact Jacobian?
limitations: limitations: I did not see any negative societal impact.
paper_summary: paper_summary: This paper studies the impact of inexactness Jacobian in second-order methods for solving variational inequalities.

This article is very informative， including a novel inexact second-order framework which achieves a convergence rate that matches the proposed lower bound, feasible solution to the sub problem, combination of quasi-Newton methods, and the generalization to the higher-order inexactness.

Overall, I found this paper well-written and the results are neat and  interesting.
strength_weakness: strength_weakness: There may lack some discussion on the related work:

In terms of the inexact Jacobian: It is very necessary to discuss the relation with one that considers the inexact Jacobian for convex optimization, especially the technical difficulties over them [3, 4].

In terms of the quasi-Newton methods for solving VI, the author may need to discuss [A] which does not require the symmetric QN update for saddle point problems.

I will also be happy to see discussion on the **technical difficulties** in comparison with optimal second-order methods for solving VI, especially [59], where Algorithm 1 use very similar framework to.


[A]. Liu C, Luo L. Quasi-Newton methods for saddle point problems. NeurIPS 2022.
questions: questions: See weakness.

In related work for convex optimization [3, 4], they also show the influence of inexact gradient on the convergence, it seems also important to discuss the behavior when $F(\cdot)$ is inexact.
limitations: limitations: N/A
paper_summary: paper_summary: The paper develops 2nd-order methods for variational inequalities under inexact Jacobian information. Existing methods requires exact derivative information which can be prohibitively expensive for large-problems limiting their practicality. To address this issue, the paper analyzes the impact of inexact derivatives on the convergence rate. When the objective is smooth and monotone, the paper establishes a lower-bound that depends upon the level inaccuracy and develop an optimal algorithm that attains this bound. The paper develops several inexact schemes based on a Quasi-Newton approximation that achieve a gloabl sublinear convergence rates.
strength_weakness: strength_weakness: The main weakness that I can see is that the paper focuses on deterministic inaccuracy in the Jacobian. Most practical ML algorithms use stochastic derivative information. Moreover, the batch size doesn't typically increase as the algorithm progresses, so generally only a noisy neighborhood of the optimum is reached. Thus, it is not immediately clear how relevant the algorithms and analysis in the paper are to this setting.
questions: questions: 1) If the algorithms in this paper only have access to stochastic Jacobian information  and assuming they don't minibatch, do the authors believe the conclusions of the paper still hold modulo the fact that convergence will only occur up to the noise level in the Jacobian? If the answer is yes, then I don't think the weakness above is that significant. 

2)  For the experiments results in Figure 1, how well do the proposed methods perform in terms of wall-clock time? Although they require more iterations then Perseus, I would expect them to be faster due to their use of inexactness. I think inclusion of this information in a table would make a nice addition to the paper.
limitations: limitations: I think the authors have addressed the main limitations of the work at the end of the paper.

Comments: title: Thanks for the rebuttal
comment: Thank you for your response. We appreciate your feedback and are glad that we could address your concerns.
Thank you for your prompt response and for raising the score. We are pleased to hear that our rebuttal addressed your concerns.
I thank the authors for their response and have no further questions.
I appreciate the authors' response, in particular to their thoughtful answer to my question about stochastic Jacobians. I agree this sounds like an interesting direction for future work. Overall, the rebuttal has addressed my questions/concerns, I'm raising my score from 7 to 8.