Let $\F$ be a finite field. Given $N$ distinct points $\bm{\alpha}_1,\ldots,\bm{\alpha}_N \in \F^m$ and associated values $\beta_1,\ldots,\beta_N \in \F$, we study the problem of computing an $m$-variate polynomial $f(\vy) \in \F[\vy]$ of individual degree less than $d$, with $d^m \ge N$, that satisfies $f(\bm{\alpha}_i)=\beta_i$ for all $i\in[N]$, whenever such an $f$ exists. This is known as the\textit{ multivariate polynomial interpolation} problem. In this work, given $N=d^m$, we present a deterministic algorithm for multivariate polynomial interpolation over an arbitrary finite field $\F$ that runs in near-quadratic time $N^{2+o(1)}\poly(m,d,\log|\F|)$, assuming $d$ is $\omega(1)$. This bound holds whenever the associated \emph{multivariate Vandermonde} matrix is nonsingular, has constant co-rank, or admits some special structure. In the general case, we give a randomized Las Vegas algorithm attaining the same near-quadratic running time, provided the finite field is sufficiently large and $d$ grows similar to the previous case. Unlike other multivariate interpolation strategies, we do not assume black-box access to the polynomial. In particular, we are not permitted to query arbitrary input points to obtain their evaluations.

Our interpolation result follows by combining Wiedemann's classical approach~\cite{w86} for solving sparse linear systems with known bounds for fast \emph{multivariate multipoint evaluation} over finite fields~\cite{umans08, KU11, BGKM23, BGGKU24}. Specifically, our deterministic interpolation strategy relies on a deterministic variant of Wiedemann's approach, which is more commonly used in the randomized setting.