Multivariate Monomial Vandermonde Matrices

• Multivariate monomial Vandermonde matrices are defined by evaluating multivariate monomials at a set of nodes, crucial for interpolation, approximation, and spectral estimation.
• They admit exact determinantal factorizations on tensor-product grids, enabling efficient computation and robust analysis in high-dimensional settings.
• Advanced techniques like Arnoldi-based orthogonalization mitigate ill-conditioning, reducing sample complexity while ensuring numerical stability in polynomial approximation.

A multivariate monomial Vandermonde matrix is a canonical linear-algebraic object encoding the evaluation of multivariate monomials at a set of nodes in $\mathbb{R}^d$ or $\mathbb{C}^d$. These matrices arise in polynomial interpolation, least-squares approximation, super-resolution, spectral estimation, and the analysis of polynomial systems. Their conditioning and invertibility are fundamental to both numerical algorithms and theoretical aspects of approximation theory, often dictating the stability and accuracy of associated computational procedures.

1. Definition and Algebraic Structure

Let $d\in\mathbb{N}$, and let $X = \{x_1, \ldots, x_M\} \subset \mathbb{R}^d$ or $\mathbb{C}^d$ be a node set. Given a multi-index set $$\mathcal{A} = \{\alpha \in \mathbb{N}_0^d : |\alpha| \leq n\}$$ (with $|\alpha| := \sum_{j=1}^d \alpha_j$ and $N = n+1$ for single-variable/sup-norm cases), the associated multivariate monomial Vandermonde matrix $V$ is defined by

$$V_{j,\alpha} = x_j^{\alpha} = \prod_{k=1}^d (x_{j,k})^{\alpha_k}, \quad j = 1, \ldots, M,~\alpha \in \mathcal{A}$$

When nodes are mapped onto $(S^1)^d \subset \mathbb{C}^d$ (i.e., $z_j = e^{2\pi i t_j}$ with $t_j \in [0,1)^d$), $V_{j,\alpha} = z_j^\alpha$ encodes trigonometric polynomials. The column-dimension $D = \#\mathcal{A}$, e.g. $D = \binom{n+d}{d}$ for total degree $\leq n$, reflects the algebraic complexity of the system (Kunis et al., 2019, Friedland et al., 20 Jan 2026, Zhu et al., 2023, Zhang et al., 2024).

2. Determinantal Factorization and Tensor-Product Grids

For almost-square tensor product grids $X \times Y$ with $n = m$ or $n = m+1$, and admissible polynomial bases with semiseparable row/column structure, multivariate Vandermonde determinants admit exact factorizations into univariate Vandermonde determinants. Specifically, for polynomial blocks $P(x)$ and $Q(y)$ as in (Marchi et al., 2013),

$$\det [P_{ij}(x_k) Q_{ij}(y_\ell)]_{(k,\ell),(i,j)} = \pm \prod_{j=1}^n \det [P_{ij}(x_k)]_{k,i} \cdot \prod_{i=1}^m \det [Q_{ij}(y_\ell)]_{\ell, j}$$

Such factorization generalizes conjectures from Padua and Padua-like points, and establishes new classes of point sets for which multivariate Vandermonde determinants reduce to products of one-dimensional determinants (Marchi et al., 2013).

3. Conditioning, Singular Values, and Node Geometry

The conditioning of multivariate Vandermonde matrices is strongly governed by the configuration of the nodes. In the analytic context (e.g. nodes on $(S^1)^d$), the minimal singular value $\sigma_{\min}$ admits lower bounds determined by cluster structures and separation statistics (Kunis et al., 2019, Kunis et al., 2021):

• For clustered nodes, introduce a cluster-complexity $C$, the maximal product of inverted pairwise distances within clusters:

$$C := \max_{j} \prod_{k \neq j,\,\|t_j - t_k\| \leq 1/N} \frac{1}{N \|t_j - t_k\|}$$

The smallest singular value satisfies, under suitable separation $\rho$ and maximal cluster size $\lambda$,

$$\sigma_{\min}(V) \gtrsim C^{-1} N^{-d/2} \lambda^{-\lambda + O(d)}$$

• In the univariate case $d=1$, the power-law decay $\sigma_{\min} \sim \tau^{\lambda - 1}$ (with $\tau$ normalized separation) precisely quantifies instability near colliding nodes.

Alternatively, for arbitrary distinct nodes $Z = \{z_1, \ldots, z_s\} \subset B_2^n$ and high-degree regime $N \geq s-1$, stability is quantified by the max-min projection separation

$$\rho(Z, j) = \max_{\|v\|_2 = 1} \min_{i \neq j} |\langle v, z_j - z_i \rangle|, \quad \kappa(Z) = \min_j \rho(Z, j)$$

One then has

$$\sigma_{\min}(V_N(Z)) \gtrsim \frac{\kappa(Z)^{s-1}}{(4n)^{s-1} s \sqrt{s \nu(n,N)}}, \quad \nu(n,N) = \binom{N + n}{N}$$

This projection-based approach avoids a priori separation assumptions and applies to all distinct node sets (Friedland et al., 20 Jan 2026).

4. Separation Criteria and Stability Regimes

A Shannon–Nyquist-type separation, $q n \geq C(d)$ with $q$ minimal node separation and $C(d)$ scaling linearly or logarithmically with $d$, yields explicit spectral bounds (Kunis et al., 2021):

• Linear regime ($C(d) \sim d$): $\sigma_{\min}(V) \gtrsim n^{d/2}$ and $\mathrm{cond}(V) = O(1)$, i.e., uniformly bounded in dimension.
• Logarithmic regime ($C(d) \sim \log d$): $\sigma_{\min}(V)/n^{d/2}$ decays rapidly, $\mathrm{cond}(V)$ grows slightly more than exponentially.

These separation conditions underpin robust stability for interpolation and super-resolution methods. Notably, as node separation decreases (especially in higher dimensions), extreme ill-conditioning becomes inevitable (Kunis et al., 2019, Kunis et al., 2021).

5. Arnoldi-Based Remedies and Polynomial Bases

The multivariate Vandermonde-Arnoldi (V+A) method overcomes exponential ill-conditioning by replacing the monomial basis with a discrete orthonormal basis constructed via a Stieltjes/Arnoldi orthogonalization process (Zhu et al., 2023, Zhang et al., 2024):

• Construct an $M \times N$ orthonormal matrix $Q$ such that $(1/M) Q^T Q = I$ on sample points $X$.
• Transform the original least-squares/interpolation problem $V c \approx f$ into $Q d \approx f$, where coefficient recovery is stable to machine precision.
• Sample complexity is $M = O(N^2)$ for deterministic admissible meshes and $M = O(N^2 \log N)$ for randomized sampling, with near-optimal error bounds.
• Leverage-score weighting further reduces the sample complexity to $M = O(N \log N)$.

For derivative-augmented data, the multivariate confluent Vandermonde with $G$–Arnoldi orthogonalization generalizes to build $G$-orthonormal bases facilitating Hermite interpolation, boundary-value enforcement, and PDE solution. The coefficient matrix in the transformed least-squares system is orthonormal, and only the recurrence matrix $H$ governs conditioning (Zhang et al., 2024).

6. Geometric and Cluster Effects Beyond Pairwise Distances

In higher dimensions ($d > 1$) and cluster sizes $\lambda \geq 3$, the minimal singular value depends nontrivially on the geometric configuration of nodes, not just pairwise distances. For example, the angle between points within clusters can alter $\sigma_{\min}$'s decay exponent—from linear ($\tau$) to quadratic ($\tau^2$) as demonstrated in specific three-point cluster constructions (Kunis et al., 2019). A plausible implication is that full geometric invariants (including angles, higher-order determinants) are necessary to precisely model conditioning in multivariate clusters.

7. Numerical and Theoretical Implications

The analysis of multivariate monomial Vandermonde matrices illuminates several foundational insights:

• Exponential ill-conditioning is generic without careful node selection or basis transformation.
• Clustered nodes, low separation, or inadequate geometric diversification of sample sets directly cause instability in polynomial interpolation, super-resolution, ESPRIT, MUSIC, and related algorithms.
• Factorization results for special grids (e.g. Padua-like, almost-square product grids) facilitate fast determinant computation and explicit construction of optimal interpolation sets (Marchi et al., 2013).
• Arnoldi-based approaches, including $G$–Arnoldi, provide robust computational frameworks for stable polynomial approximation, even in high dimension or for derivative-enriched datasets (Zhu et al., 2023, Zhang et al., 2024).
• Practical sample complexity results enable stable spectral approximation on irregular domains with minimal oversampling.

In summary, multivariate monomial Vandermonde matrices exhibit rich algebraic, geometric, and numerical properties, combining subtle combinatorics of multi-index systems, geometry-driven stability estimates, and modern computational strategies for polynomial approximation and inverse problems.