Generalized Low-Rank Matrix Approximation

• Generalized low-rank matrix approximation is a method that decomposes high-dimensional matrices into low-rank surrogates, reducing storage and speeding up computations.
• It embeds low-rank representations within iterative techniques like Neumann series to replace costly full matrix operations with efficient series evaluations.
• The approach is applied in uncertainty quantification, signal processing, PDE solvers, and wireless communications while maintaining controlled error propagation.

Generalized low-rank matrix approximation (GLRMA) encompasses a class of techniques for representing high-dimensional matrices by low-rank surrogates, often in combination with analytic expansions to solve or accelerate associated large-scale linear algebra problems. Such formulations typically arise in uncertainty quantification, fast signal processing, PDE solvers, wireless communications, and scientific computing when one must invert, factor, or manipulate structured or stochastic matrices at scale. Modern approaches frequently embed low-rank representations within Neumann series expansions or other iterative methods, enabling both memory savings and performance improvements while maintaining mathematically rigorous control on error propagation.

1. Low-Rank Approximation and Matrix Decomposition Principles

GLRMA operates by decomposing an arbitrary matrix $A \in \mathbb{R}^{N \times N}$ (typically arising as the discretized operator in stochastic PDEs, graph Laplacians, or Gram matrices in MIMO) as $A = \bar{A} + \widetilde{A}$, where $\bar{A}$ is a mean or fixed component and $\widetilde{A}$ encodes random or structured perturbations. To reduce storage and computational costs, $\widetilde{A}$ is approximated by a rank-$k$ representation $\widetilde{A}^* = U V^T$ ($U, V \in \mathbb{R}^{N \times k}$, $k \ll N$), as in randomized SVD or generalized low-rank approximation algorithms (Zhu et al., 16 Jan 2026). The associated inversion or solution task is then reduced from an $O(N^3)$ direct solve to a sequence of operations involving $O(Nk)$ matrix-vector products.

In graph signal processing or bandlimited reconstruction, low-rank filters are used to represent ideal or approximate projections onto informative subspaces, bypassing the need for full eigendecomposition by leveraging Chebyshev polynomial approximations or sparse matrix expansions (Wang et al., 2018).

2. Neumann Series Expansion for Matrix Inversion

Whenever the matrix inverse $(I + B)$ is needed and $\|B\| < 1$ (for a consistent operator norm), the classical Neumann series is invoked:

$(I + B)^{-1} = \sum_{r=0}^{\infty} (-B)^r.$

This expansion transforms inversion into repeated application of the low-rank perturbation, which is highly efficient when $B$ is low-rank and small in norm. For practical computation, the series is truncated after $R$ terms to give:

$(I + B)^{-1} \approx \sum_{r=0}^{R} (-B)^r,$

with tail error bounded by

$$\|R_R\| \le \frac{\|B\|^{R+1}}{1 - \|B\|}\|\bar{A}^{-1}\|$$

where $B = \bar{A}^{-1}\widetilde{A}^*$ (Zhu et al., 16 Jan 2026). This principle generalizes to more complex transmission operators (e.g., the Neumann-Poincaré operator in boundary integral equations (Cherkaev et al., 2020, Choi, 2024)) and is also the core of Krylov and algebraic multigrid smoothers (Thomas et al., 2021), where triangular solves and projection steps are replaced by SpMV-based series evaluations.

3. Algorithmic Embedding: Fast Iterative Solvers and Preconditioners

GLRMA is integrated into iterative solvers as follows:

• Precompute the low-rank factors $U, V$ for $\widetilde{A}_m$ (or equivalent perturbation).
• Solve the mean problem once, i.e., compute $\bar{u}_l = \bar{A}^{-1}b_l$ for each right-hand side.
• For each sample or realization, compute $B_m = \bar{A}^{-1}U V^T$ and update the solution as:

$$w^{(0)} = \bar{u}_l,\quad w^{(r+1)} = - Y (V_m^T w^{(r)})\quad (Y = \bar{A}^{-1}U),$$

accumulating the sum $\sum_{r=0}^R (-B_m)^r \bar{u}_l$ (Zhu et al., 16 Jan 2026).

In graph sampling, ranking and selection are based on surrogate objectives derived from truncated Neumann series in terms of appropriate graph filters ($T^{\mathrm{Poly}}$ via Chebyshev expansion), avoiding matrix inversion while providing theoretical control on MSE error (Wang et al., 2018).

In Krylov GMRES and AMG, truncated Neumann series can replace direct solves of $(I + L)$ by $I - L + L^2 - ...$, preserving backward stability and convergence (Thomas et al., 2021).

4. Convergence, Error Analysis, and Complexity Bounds

Convergence of the Neumann series requires $\rho(B) < 1$ (spectral radius); in low-rank settings, this is assured when random/stochastic perturbations are bounded, or when diagonal dominance holds (as in massive MIMO scenarios (Zhu et al., 2015, Dimitrov et al., 2017)). Truncation error decays exponentially with the number of terms, and for sufficiently small $\|B\|$,

$\text{Error} \sim \|B\|^{R+1}$

allowing controlled trade-off between runtime and accuracy (Zhu et al., 16 Jan 2026, Wang et al., 2018). In massive MIMO, closed-form MSE error formulas in terms of antenna-user ratio $\beta = M/K$ and beta-function terms $B_{a,M}$ enable precise choice of $N$ (Zhu et al., 2015).

Complexity can be further reduced using optimized factorization strategies for Neumann series evaluation, such as prime-base splitting and mixed-basis recursion. These algorithms lower the multiplication count from $O(N)$ to $O(\log N)$, with practical speedups of $1.8\times$ to $2.5\times$ over classical Horner evaluation for matrix sizes of several hundred or more (Dimitrov et al., 2017).

5. Applications Across Scientific and Engineering Domains

GLRMA with Neumann-series embedding is utilized in:

• Uncertainty quantification and PDE solvers, for efficient inverse computation in stochastic systems (Zhu et al., 16 Jan 2026).
• Graph signal processing: optimal node sampling and robust bandlimited signal reconstruction, leveraging low-rank filter surrogates (Wang et al., 2018).
• Large-scale power-flow studies in distribution networks, via accelerated probabilistic solvers for Newton steps (Chevalier et al., 2020).
• Wireless communications: matrix inversion approximation (MIA) for precoding/detection in massive MIMO systems, with rigorous performance-complexity analysis (Zhu et al., 2015, Dimitrov et al., 2017).
• Preconditioning and smoothers in iterative solvers (GMRES/AMG): replacing triangular solves by SpMV-based Neumann expansions, preserving stability and improving parallel scaling (Thomas et al., 2021).

6. Implementation Considerations and Numerical Performance

Efficient GLRMA implementations rely on:

• Exploiting low-rank structure for storage/memory savings: $O(Nk)$ instead of $O(N^2)$ per sample.
• Precomputing and reusing core operations (e.g., $\bar{A}^{-1}$, filter polynomials, and $Y = \bar{A}^{-1}U$).
• Careful control of truncation ($R$) and rank ($k$) for accuracy-runtime tradeoff.
• Hardware-optimized factorization trees and SpMV pipelines, amenable to GPUs and many-core architectures (Thomas et al., 2021).
• Reordering to exploit sparsity and symmetry, thus reducing non-normality of factors for swift convergence.

Reported empirical results include $20$–$40\%$ runtime savings and $10^{-4}$–$10^{-5}$ relative error for stochastic PDEs (Zhu et al., 16 Jan 2026), up to $10\times$–$167\times$ speedup in probabilistic power-flow studies (Chevalier et al., 2020), and $2\times$ multiplication cost reduction in Neumann-series matrix inversion (Dimitrov et al., 2017).

7. Theoretical Connections to Spectral Analysis and Boundary Operators

GLRMA is conceptually unified with geometric series expansions in spectral theory; for example, boundary-integral Neumann–Poincaré operators admit explicit infinite-matrix forms in specialized bases, with exponential off-diagonal decay and block-diagonalization under symmetry (Choi, 2024, Cherkaev et al., 2020). Such representations inform both theoretical bounds (e.g., extremal conductivity via conformal coefficients, spectral monotonicity under domain deformation) and practical approximations in spectral computations.

---

References cited:

(Zhu et al., 16 Jan 2026): "An efficient solver based on low-rank approximation and Neumann matrix series for unsteady diffusion-type partial differential equations with random coefficients" (Wang et al., 2018): "A-Optimal Sampling and Robust Reconstruction for Graph Signals via Truncated Neumann Series" (Zhu et al., 2015): "On the Matrix Inversion Approximation Based on Neumann Series in Massive MIMO Systems" (Dimitrov et al., 2017): "On the Computation of Neumann Series" (Chevalier et al., 2020): "Accelerated Probabilistic Power Flow in Electrical Distribution Networks via Model Order Reduction and Neumann Series Expansion" (Thomas et al., 2021): "Neumann Series in GMRES and Algebraic Multigrid Smoothers" (Choi, 2024): "Matrix representation of the Neumann-Poincaré operator for a torus" (Cherkaev et al., 2020): "Geometric series expansion of the Neumann-Poincaré operator: application to composite materials"