This paper settles an open and challenging question pertaining to the design of simple and optimal high-order methods for solving smooth and monotone variational inequalities (VIs). A VI involves finding x ⋆is an element of X\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x<^>\star \in {\mathcal {X}}$$\end{document} such that ⟨F(x),x-x ⋆⟩>= 0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle F(x), x - x<^>\star \rangle \ge 0$$\end{document} for all x is an element of X\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x \in {\mathcal {X}}$$\end{document}. We consider the setting in which F:Rd -> Rd\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F: {\mathbb {R}}<^>d \rightarrow {\mathbb {R}}<^>d$$\end{document} is smooth with up to (p-1)th\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(p-1)<^>{\text {th}}$$\end{document}-order derivatives. For p=2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p = 2$$\end{document}, the cubic regularization of Newton's method has been extended to VIs with a global rate of O(epsilon-1)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O(\epsilon <^>{-1})$$\end{document} (Nesterov in Cubic regularization of Newton's method for convex problems with constraints, Tech. rep., Universite catholique de Louvain, Center for Operations Research and Econometrics (CORE), 2006). An improved rate of O(epsilon-2/3loglog(1/epsilon))\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O(\epsilon <^>{-2/3}\log \log (1/\epsilon ))$$\end{document} can be obtained via an alternative second-order method, but this method requires a nontrivial line-search procedure as an inner loop. Similarly, the existing high-order methods based on line-search procedures have been shown to achieve a rate of O(epsilon-2/(p+1)loglog(1/epsilon))\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O(\epsilon <^>{-2/(p+1)}\log \log (1/\epsilon ))$$\end{document} (Bullins and Lai in SIAM J Optim 32(3):2208-2229, 2022; Jiang and Mokhtari in Generalized optimistic methods for convex-concave saddle point problems, 2022; Lin and Jordan in Math Oper Res 48(4):2353-2382, 2023). As emphasized by Nesterov (Lectures on convex optimization, vol 137, Springer, Berlin, 2018), however, such procedures do not necessarily imply the practical applicability in large-scale applications, and it is desirable to complement these results with a simple high-order VI method that retains the optimality of the more complex methods. We propose a pth\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p<^>{\text {th}}$$\end{document}-order method that does not require any line search procedure and provably converges to a weak solution at a rate of O(epsilon-2/(p+1))\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O(\epsilon <^>{-2/(p+1)})$$\end{document}. We prove that our pth\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p<^>{\text {th}}$$\end{document}-order method is optimal in the monotone setting by establishing a lower bound of omega(epsilon-2/(p+1))\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega (\epsilon <^>{-2/(p+1)})$$\end{document} under a generalized linear span assumption. A restarted version of our pth\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p<^>{\text {th}}$$\end{document}-order method attains a linear rate for smooth and pth\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p<^>{\text {th}}$$\end{document}-order uniformly monotone VIs and another restarted version of our pth\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p<^>{\text {th}}$$\end{document}-order method attains a local superlinear rate for smooth and strongly monotone VIs. Further, the similar pth\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p<^>{\text {th}}$$\end{document}-order method achieves a global rate of O(epsilon-2/p)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O(\epsilon <^>{-2/p})$$\end{document} for solving smooth and nonmonotone VIs satisfying the Minty condition. Two restarted versions attain a global linear rate under additional pth\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p<^>{\text {th}}$$\end{document}-order uniform Minty condition and a local superlinear rate under additional strong Minty condition.
