🔢 mathematics

Stable High-Order Interpolation on the Grassmann Manifold by Maximum-Volume Coordinates and Arnoldi Orthogonalization

This paper proposes a stable high-order interpolation framework for the Grassmann manifold that combines Maximum-Volume local coordinates to bypass costly SVD-based mappings with Arnoldi-orthogonalized polynomial bases to avoid ill-conditioned linear systems, thereby achieving accurate and efficient geometric approximation.

Original authors: Qiang Niu, Wen Jiang, Jie Fei, Ruoyu Xiong, Yuxuan Li

Published 2026-03-26

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Qiang Niu, Wen Jiang, Jie Fei, Ruoyu Xiong, Yuxuan Li

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are trying to draw a smooth, winding road on a map, but the map isn't a flat piece of paper. Instead, it's a globe (or a complex, curved surface). You have a few specific checkpoints (like cities) where you know exactly where the road is and how fast it's turning. Your goal is to draw the entire road connecting these points so that a driver can follow it perfectly.

This is the problem scientists face when dealing with Grassmann Manifolds. In the real world, this "road" represents things like the changing shape of a fluid flow, the orientation of a robot's arm, or the most important features of a video. The "checkpoints" are data points collected from simulations or sensors.

The paper by Niu, Jiang, and colleagues introduces a new, super-stable way to draw this road, even when the road twists and turns wildly. Here is how they do it, explained simply:

The Problem: The "Flat Map" Trap

Usually, to draw a road on a globe, people try to flatten it out onto a piece of paper (a flat map) to do the math.

The Distortion: When you flatten a globe, things get stretched and distorted. If you try to draw a straight line on the paper, it might look like a curve on the globe.
The Math Explosion: To draw a smooth curve through many points, you need high-level math (polynomials). But on these "flat maps," the math gets incredibly unstable. It's like trying to balance a tower of Jenga blocks where the bottom blocks are wobbly. As you add more points (make the tower taller), the whole thing collapses into nonsense. This is called ill-conditioning.

The Solution: The "MV-(C)V+A" Framework

The authors built a two-part machine to solve this. Think of it as a Stabilizer and a Smart Ruler.

Part 1: The Stabilizer (Maximum-Volume Coordinates)

The Old Way: Imagine trying to flatten the globe by picking a random spot to cut and stretch. Sometimes you pick a spot that stretches the map so much that the cities look huge and distorted. This makes the math impossible.
The New Way (MV): The authors use a clever trick called Maximum-Volume (MV) coordinates. Instead of randomly picking a spot to flatten the map, they use a "Householder" tool (a mathematical mirror) to rotate the globe so that the area they are looking at is perfectly flat and undistorted.
The Analogy: It's like a camera operator who doesn't just zoom in on a wobbly hand; they physically rotate the camera to find the perfect angle where the hand looks steady and square. This ensures the "map" they are drawing on is never stretched or squashed, keeping the math clean.

Part 2: The Smart Ruler (Arnoldi Orthogonalization)

The Old Way: Once the map is flat, they try to draw the curve using a standard ruler (monomial basis). But as the curve gets more complex (higher degree), the ruler starts to shake violently. The numbers get so big and messy that the computer gives up.
The New Way (V+A and CV+A): They switch to a Smart Ruler based on the Arnoldi method.
- Instead of using a shaky ruler, they build a set of perfectly balanced, interlocking blocks (orthogonal polynomials).
- Because these blocks are perfectly balanced, you can stack them as high as you want (high-degree interpolation) without the tower falling over.
- The "Confluent" Twist: Sometimes, you don't just know where the road is; you also know how fast it's turning (velocity). The authors extended their Smart Ruler to handle these "turning speeds" without breaking the balance.

How It Works Together (The Three-Stage Process)

Stage 1: The Setup (Geometric Stabilization)
They take the messy, curved data and use the "Householder" tool to rotate it into a nice, flat, stable coordinate system. No more stretching or distortion.
Stage 2: The Drawing (Algebraic Interpolation)
In this flat, stable world, they use the "Smart Ruler" (Arnoldi method) to draw a perfect, high-degree curve through the points. Because the ruler is balanced, the math stays precise, even with hundreds of points.
Stage 3: The Return (Reconstruction)
They take the perfect curve they drew on the flat paper and "wrap" it back onto the globe. Because the setup was so stable, the road on the globe is smooth, accurate, and doesn't have any weird kinks.

Why This Matters

In the past, if you tried to predict complex systems (like weather patterns or how a bridge vibrates) using many data points, the math would crash, and the predictions would be garbage.

This new method allows scientists to:

Go Higher: They can use much more complex curves (higher degrees) to capture fine details.
Stay Stable: The math doesn't explode, even with noisy data.
Save Time: It avoids the heavy, slow calculations (like SVD) that usually slow down these simulations.

In a nutshell: The authors figured out how to flatten a curved world without stretching it, and then used a super-stable ruler to draw a perfect line through it. This means we can now model complex, changing systems with incredible accuracy, something that was previously too risky to try.

1. Problem Statement

The paper addresses the challenge of performing high-order interpolation (both Lagrange and Hermite) on the Grassmann manifold $\text{Gr}(n, p)$ , which represents the space of $p$ -dimensional subspaces in $\mathbb{R}^n$ . This problem is critical in applications like parametric Model Order Reduction (MOR), computer vision, and machine learning.

The authors identify two primary numerical bottlenecks that hinder high-order interpolation on manifolds:

Geometric Instability: Traditional methods rely on Riemannian normal coordinates, which require Singular Value Decomposition (SVD), matrix logarithms, and exponentials. These operations are computationally expensive and, crucially, require differentiating through matrix decompositions for Hermite interpolation (velocity matching), leading to numerical sensitivity and instability.
Algebraic Instability: Once data is mapped to a local Euclidean coordinate space, standard polynomial interpolation relies on monomial bases (Vandermonde matrices). The condition number of these matrices grows exponentially with the polynomial degree, causing coefficient blow-up and catastrophic loss of precision, especially when derivative data (confluent Vandermonde systems) is included.

2. Methodology: The MV-(C)V+A Framework

The authors propose a decoupled, three-stage framework that integrates Maximum-Volume (MV) coordinates with Arnoldi-orthogonalized polynomial bases. This approach separates geometric nonlinearity from algebraic ill-conditioning.

Stage I: Geometric Stabilization (Householder-MV Coordinates)

Instead of using discrete row permutations (as in previous MV approaches) which can break continuity, the authors introduce a continuous orthogonal transformation:

Reference Selection: A reference subspace $U_m$ is selected (typically the last sample).
Householder QR: A Householder QR decomposition is performed on $U_m$ to generate an orthogonal matrix $Q$ . This acts as a continuous rotation rather than a discrete swap.
Coordinate Extraction: All sample subspaces are transformed by $Q^T$ . The local coordinates $\Xi$ are extracted algebraically as $\Xi = \tilde{U}_2 \tilde{U}_1^{-1}$ , where $\tilde{U}$ is the transformed Stiefel representative.
Derivative Mapping: Tangent vectors (velocities) are mapped to the coordinate domain using a purely algebraic chain rule formula, avoiding the differentiation of SVDs.
Benefit: This ensures the geometric condition number $\kappa_{geo}$ remains near the theoretical lower bound ( $\sqrt{p}$ ), preventing the amplification of coordinate errors during reconstruction.

Stage II: Algebraic Interpolation (V+A and CV+A)

In the Euclidean coordinate domain, the authors replace the unstable monomial basis with the Vandermonde-with-Arnoldi (V+A) method and its Confluent (CV+A) extension:

Krylov Subspace Construction: The method constructs a discrete orthonormal polynomial basis directly from the parameter nodes using the Arnoldi process.
Orthogonalization: Instead of solving an ill-conditioned Vandermonde system $Vc = f$, the method solves a well-conditioned system $Q^T Q = I$ (where $Q$ is the orthonormal basis).
Confluent Extension: For Hermite interpolation, the method uses an augmented operator to simultaneously orthogonalize function values and derivatives, maintaining a condition number of $\kappa_2(Q) = 1$ .
Evaluation: Evaluation at new points is performed via a stable recurrence relation driven by an upper Hessenberg matrix, avoiding explicit power calculations.

Stage III: Geometric Reconstruction

The interpolated Euclidean coordinates are mapped back to the Grassmann manifold:

Cholesky Retraction: The coordinate matrix $\Xi$ is reconstructed into a Stiefel representative $\tilde{U}$ using a Cholesky factorization of $(I + \Xi^T \Xi)$ , avoiding expensive SVDs or matrix square roots.
Inverse Transformation: The global orthogonal map $Q$ is applied to return the subspace to its original orientation.

3. Key Contributions

Stabilized Geometric Mapping: Introduction of a Householder-stabilized MV coordinate system that provides continuous, well-conditioned local representations without discrete permutations or SVD computations.
Algebraic Robustness: Adaptation of the V+A and CV+A methods to manifold interpolation, eliminating the exponential growth of condition numbers associated with high-degree polynomial fitting.
Theoretical Error Bounds: Rigorous derivation of total error bounds showing that the interpolation error decomposes into geometric and algebraic components. The authors prove that the method achieves spectral convergence and limits error growth to machine precision ( $O(\epsilon_{mach})$ ).
Derivative Handling: A purely algebraic formulation for mapping manifold velocities to coordinate derivatives, bypassing the numerical instability of differentiating matrix factorizations.

4. Numerical Results

The paper validates the framework through extensive experiments:

High-Degree Stability: On $\text{Gr}(1000, 10)$ , the proposed method successfully performed degree-35 Hermite interpolation. In contrast, traditional methods (Normal, Local, and original MV coordinates) suffered complete algebraic collapse with errors exploding to $O(100)$ due to condition numbers exceeding $10^{16}$ .
Noise Robustness: In the presence of $10^{-10}$ noise, traditional methods amplified errors significantly due to high condition numbers. The MV-(C)V+A method maintained reconstruction errors near the noise floor ( $\approx 10^{-10}$ ), demonstrating optimal conditioning.
Application to PDEs: The method was applied to a parametric reduced-order model for the 1D Helmholtz equation (degree 23 Hermite interpolation). It maintained relative subspace errors near $10^{-10}$ , whereas baseline methods failed to converge.
Orthogonality Preservation: The reconstructed subspaces strictly preserved orthonormality, with orthogonality loss remaining at machine precision levels.

5. Significance

This work provides a robust, high-order interpolation framework for Grassmann manifolds that overcomes the historical trade-off between geometric fidelity and algebraic stability.

Enables High-Order ROM: It allows for the construction of highly accurate, derivative-enhanced surrogate models for complex physical systems (like the Helmholtz equation) where high-degree polynomials are necessary to capture nonlinear behaviors.
Computational Efficiency: By avoiding SVDs and using efficient recurrence relations, the method is computationally feasible for large-scale problems.
Generalizability: The framework is applicable to any manifold-valued data requiring high-order approximation and can be extended to other matrix manifolds (e.g., Stiefel, Rotation groups) and multivariate parameter domains.

In summary, the MV-(C)V+A framework represents a significant advancement in geometric numerical analysis, enabling reliable, machine-precision interpolation on manifolds where previous methods were limited to low-degree approximations.