Source-linked AI summary
Near-Quadratic Time Multivariate Polynomial Interpolation
Chandra Kanta Mohapatra
TL;DR
Dense multivariate interpolation at fixed points must handle Vandermonde systems that may be singular, while generic dense solving lacks a truly quadratic-time algorithm. The paper combines fast multivariate matrix-vector multiplication with Wiedemann-style system solving to compute interpolants in near-quadratic time. The algorithm is deterministic for specified rank or structure cases and Las Vegas randomized in the general case.
Problem
Dense multivariate interpolation at arbitrary fixed points requires solving a Vandermonde system that may be singular, unlike the unique-solution univariate case.
Source: §1 · pp. 1–2Method
The algorithm combines Wiedemann-style iterative linear-system solving with fast multivariate Vandermonde matrix-vector multiplication.
Source: §1 · p. 2, §6.2 · p. 9Results
N^{2+o(1)} poly(m, d, log |F|) time is achieved for N = d^m, deterministically in full-rank, constant-co-rank, or specially structured cases and otherwise with a Las Vegas algorithm.
Source: §1 · p. 3Takeaways & Limitations
The algebraic structure of generalized Vandermonde matrices yields an improvement over generic dense linear algebra even for unstructured evaluation sets.
Source: §1 · p. 3In the randomized singular-system procedure, a consistent system may be reported inconsistent if rank computation fails, but an output solution is never incorrect.
Source: §6.5.2 · pp. 18, 20
Abstract
from arXiv · showhide
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.
1 Introduction
The paper studies dense multivariate interpolation at fixed evaluation points, where distinct points need not yield a unique solution and generic dense linear-system solving lacks a truly quadratic-time algorithm. It combines fast multivariate matrix-vector multiplication with Wiedemann-style linear algebra to obtain near-quadratic interpolation algorithms, with determinism depending on the system’s rank and structure.
1 Introduction
- Multivariate interpolation seeks a degree-bounded polynomial matching N distinct point-value constraints, with a solution required only when one exists.The individual degree is less than d, and the interpolation space has dimension d^m.Source: §1 · pp. 1–2
- Unlike univariate interpolation, distinct multivariate evaluation points can produce a singular Vandermonde matrix, so a solution may be nonunique or nonexistent.For simplicity, the paper takes N = d^m, padding the system with vacuous equations when N is smaller.Source: §1 · p. 2
- Fixed-point dense interpolation differs from black-box methods, which can query chosen points, and from fast grid methods, which require structured evaluation points.The paper also notes that no near-linear algorithm is known for the general bivariate case over arbitrary evaluation points.Source: §1 · p. 2
- The approach exploits the multivariate Vandermonde matrix’s algebraic structure to improve on generic dense solving, whose best-known time is O(N^ω).A fast matrix-vector product enables iterative linear algebra to solve the structured system more quickly.Source: §1 · pp. 1–2
- Fast multivariate multipoint evaluation supplies matrix-vector products, while Wiedemann’s method reduces solving the resulting system to computing a Krylov-sequence minimal polynomial.The paper combines these existing ingredients rather than claiming that the component ideas are new.Source: §1 · pp. 3–4
- For a finite field, given N evaluation constraints for an unknown m-variate polynomial of individual degree less than d, a solution, if one exists, is computable in N2+o(1) poly(m, d, log |F|) time.The algorithm is deterministic when the associated multivariate Vandermonde matrix has full rank, constant co-rank, or additional structure; otherwise it is Las Vegas randomized for the general case.Source: §1 · p. 3
2 Proof Overview
The proof overview turns interpolation into structured linear-system solving: fast matrix-vector products make Wiedemann-style methods efficient, while separate procedures handle nonsingular and singular systems. A deterministic coordinate-sequence variant covers the nonsingular case; recursive reductions and randomized preconditioning address singular cases.
2 Proof Overview
- Fast matrix-vector multiplication lets Wiedemann’s iterative framework solve a system without explicitly inverting its matrix.The method is useful whenever matrix-vector multiplication is faster than the general dense product.Source: §2 · pp. 4–5
- An annihilating polynomial p with p(A) = 0 yields a linear relation among b, Ab, ..., A^d b, from which a solution can be recovered using d matrix-vector products.The characteristic polynomial guarantees an annihilating polynomial of degree at most n, but a lower-degree one reduces the required products.Source: §2 · p. 4
- The minimal polynomial of the Krylov sequence {A^i b} has degree at most n and supports solution recovery for nonsingular A.Using this sequence avoids computing the matrix’s minimal polynomial directly.Source: §2 · p. 5
Nonsingular system.
- For nonsingular systems, the deterministic variant computes each coordinate sequence’s minimal polynomial and takes their least common multiple to recover the vector-sequence minimal polynomial.Berlekamp–Massey is applied independently to the coordinate sequences; combined with fast multivariate Vandermonde products, this gives near-quadratic interpolation.Source: Nonsingular system · p. 5
Singular system.
- For singular A, a Krylov-sequence minimal polynomial with nonzero constant term still yields a valid solution, because b lies in A’s column space.The resulting candidate can be verified directly, so no separate invertibility certificate is required.Source: Singular system · pp. 5–6
- When that constant term is zero, the method derives a column dependence and recursively examines the transpose to remove redundant rows and columns.Constant co-rank requires only constantly many reductions before solving a full-rank subsystem; larger co-rank uses randomized triangular-Toeplitz preconditioning and a consistency check, still in near-quadratic time.Source: Singular system · p. 6
3 Related Work
Prior multivariate interpolation methods commonly rely on structured points, black-box access, or other specialized conditions. This work instead treats fixed arbitrary points and the resulting Vandermonde system directly, including rank-deficient cases.
3 Related Work
- Grid interpolation uses recursive univariate interpolation for near-linear time, while lower-set methods exploit structured Newton-type bases; both assume point geometry fixed in advance.This paper makes no assumption about the input points’ structure.Source: §3 · p. 6
- Sparse black-box interpolation methods can choose evaluation points adaptively, unlike the fixed-point setting considered here.The literature also includes Hermite interpolation and ideal-theoretic methods, which may involve costly Gröbner-basis or large linear-algebra computations.Source: §3 · p. 7
- The paper trades the black-box setting’s appeal for determinism by using quadratic space, since the input matrix can be explicitly determined from the points and monomial basis.It describes its contribution as combining existing ideas and necessary tools for this interpolation setting.Source: §3 · p. 7
- For fixed evaluation points, the paper reports no prior near-quadratic algorithm for dense multivariate interpolation without sparsity, black-box access, or special point structure.The cited closest fixed-point approach selects a nonsingular monomial subset and uses O(t^{ω+1}) arithmetic operations.Source: §3 · pp. 7–8
4 Organization
The section outlines the paper’s progression from multivariate Vandermonde matrix-vector products to interpolation-system solving.
4 Organization
- Sections 6.1 and 6.2 address computing matrix-vector products and transposed matrix-vector products for a multivariate Vandermonde matrix.Source: §4 · p. 8
- Section 6.3 presents the deterministic Wiedemann computation.Source: §4 · p. 8
- Sections 6.4 and 6.5 address interpolation-system solving for nonsingular and singular systems, respectively.Source: §4 · p. 8
5 Preliminaries and Notation
The preliminaries establish notation and introduce linear recurrences, matrix polynomials, Krylov sequences, and Berlekamp–Massey.
5 Preliminaries and Notation
- The paper defines soft-O notation as O(m log^c m), where c is a constant.Source: §5 · p. 8
- A sequence has a linear recurrence of order d when an annihilating polynomial of degree d gives the stated zero-sum relation for all i ≥ 0.Source: §5 · p. 8
- For a matrix A, the characteristic polynomial is det(xI − A), while the minimal polynomial is the least-degree monic polynomial satisfying μA(A) = 0.Source: §5 · p. 8
- The Krylov sequence K(A, v) is the set of vectors A^i v for i ∈ N, and its minimal polynomial is denoted μA.Source: §5 · p. 8
- Berlekamp–Massey computes the minimal polynomial from the first 2d terms of a scalar recurrence of degree at most d in eO(d) field operations [GG03, Chapter 12].Source: §5 · p. 8
- The paper points to [GG03, Chapter 12] for an exposition of the introduced concepts.Source: §5 · p. 8
6 Main Results
The section develops near-quadratic multivariate interpolation using fast multivariate Vandermonde products and Krylov-sequence methods. It gives deterministic methods for full-rank and certain singular systems, including constant co-rank, and randomized methods for larger co-rank that never return an incorrect solution, with a deterministic method when rank is sufficiently close to full.
6.1 Matrix-vector multiplication complexity
- A multivariate Vandermonde matrix or any submatrix supports deterministic matrix-vector multiplication in N^{1+o(1)} poly(m, d, log |F|) time when N = d^m and d = ω(1).The multiplication is a restricted multipoint-evaluation problem, using the algorithm of Bhargava et al..Source: §6.1 · p. 9
6.2 Transposed multivariate multipoint evaluation
- The transpose of multivariate multipoint evaluation, including its submatrices, has a deterministic running time of (d^m + N)^{1+o(1)} poly(m, d, log |F|) for sufficiently large d.The method transposes stages of the non-algebraic algorithm of Bhargava et al., avoiding parameter constraints of algebraic alternatives.Source: §6.2 · pp. 10–11
6.3 Deterministic Krylov minimal polynomial
- The minimal polynomial annihilating a Krylov sequence can be computed deterministically in eO(n · M(n) + n^2) time.The method applies Berlekamp–Massey to coordinate sequences and takes the least common multiple of their minimal polynomials.Source: §6.3 · pp. 11–13
6.4 Nonsingular system
- For full-rank multivariate Vandermonde systems, the unique interpolant is computable deterministically in N^{2+o(1)} poly(m, d, log |F|) time when d = ω(1).The result applies to N = d^m distinct data points and polynomials of individual degree less than d.Source: §6.4 · p. 14
6.5 Singular system
- For a singular system, a nonzero constant term in the Krylov-sequence minimal polynomial guarantees a solution computable deterministically in N^{2+o(1)} poly(m, d, log |F|) time.The solution need not be unique; the condition implies the right-hand side lies in the matrix's column space.Source: §6.5 · pp. 14–15
6.5.1 Case 1 : co-rank(A) = O(1)
- For rank n−O(1) systems with fast transposed matrix-vector multiplication, the section introduces a co-rank O(1) case that prunes redundant rows and columns toward an invertible system.The cited theorem statement begins with these conditions; its conclusion is not present in this part.Source: §6.5.1 · p. 15
- For co-rank O(1), interpolation is deterministic and takes N^2+o(1) poly(m, d, log |F|) time, assuming d = ω(1) and a solution exists.The theorem assumes the Krylov-sequence minimal polynomial has zero constant term and the Vandermonde matrix has rank N−O(1).Source: §6.5.1 · pp. 16–17
6.5.2 Case 2 : co-rank = ω(1)
- If a solution exists and d = ω(1), the algorithm computes it in N^2+o(1) poly(m, d, log |F|) time using 4N−2 random elements from F*, and never outputs an incorrect solution.The stated condition is that the Krylov-sequence minimal polynomial has a zero constant term.Source: §6.5.2 · pp. 20–21
- Taking |F| = Ω(N^2) makes the failure probability small; for smaller fields, the method uses a field extension, and repetition reduces failure exponentially on consistent systems.The only possible error is declaring a consistent system inconsistent.Source: §6.5.2 · p. 21
- When rank is guaranteed sufficiently close to full, the system can instead be solved deterministically without knowing its exact rank in advance.This regime generalizes the constant co-rank case.Source: §6.5.2 · p. 21