  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^\star \in \mathcal{X}$
such that $\langle F(x), x - x^\star\rangle \geq 0$ for all $x \in
\mathcal{X}$. We consider the setting in which $F$ is smooth with up to
$(p-1)^{th}$-order derivatives. For $p = 2$, the cubic regularized Newton
method was extended to VIs with a global rate of $O(\epsilon^{-1})$. An
improved rate of $O(\epsilon^{-2/3}\log\log(1/\epsilon))$ can be obtained via
an alternative second-order method, but this method requires a nontrivial
line-search procedure as an inner loop. Similarly, high-order methods based on
line-search procedures have been shown to achieve a rate of
$O(\epsilon^{-2/(p+1)}\log\log(1/\epsilon))$. As emphasized by Nesterov,
however, such procedures do not necessarily imply practical applicability in
large-scale applications, and it would be desirable to complement these results
with a simple high-order VI method that retains the optimality of the more
complex methods. We propose a $p^{th}$-order method that does \textit{not}
require any line search procedure and provably converges to a weak solution at
a rate of $O(\epsilon^{-2/(p+1)})$. We prove that our $p^{th}$-order method is
optimal in the monotone setting by establishing a matching lower bound under a
generalized linear span assumption. Our method with restarting attains a linear
rate for smooth and strictly monotone VIs and a local superlinear rate for
smooth and strongly monotone VIs. Our method also achieves a global rate of
$O(\epsilon^{-2/p})$ for solving smooth and nonmonotone VIs satisfying the
Minty condition and when augmented with restarting it attains a global linear
and local superlinear rate for smooth and nonmonotone VIs satisfying the
strictly/strong Minty condition.
