Neumann-Series Decomposition Overview

• Neumann-series decomposition is a method for expressing the inverse of an operator or matrix as an infinite geometric series, providing explicit error bounds when truncating the series.
• It is widely applied in high-dimensional regression adjustments, differential equation solvers, and quantum error mitigation to achieve efficient, high-accuracy computations.
• Algorithmic optimizations, such as basis factorization and recursive schemes, enhance the method's scalability in iterative solvers and multigrid approaches.

A Neumann-series decomposition is a systematic expansion founded on the Neumann series for resolvents of linear operators or matrices; it appears across functional analysis, numerical linear algebra, quantum information, regression adjustment for randomized experiments, and the explicit construction of solutions to ordinary differential equations. The central principle is to expand the inverse $(I - M)^{-1}$ (or more generally, the resolvent $(zI - M)^{-1}$) as a geometric series in powers of $M$ (or related operators) that converge under suitable spectral conditions. This decomposition enables explicit control of remainders, facilitates high-accuracy numerical computation, and yields refined theoretical analyses in high-dimensional regimes. The following sections review key instances, structure, algorithms, scaling laws, and implications across major research domains.

1. Algebraic Structure of the Neumann Series

The Neumann series for a bounded linear operator $M$ (typically a matrix) with $\|M\| < 1$ is

$(I - M)^{-1} = \sum_{k=0}^{\infty} M^k,$

which converges absolutely in any induced operator norm if and only if the spectral radius of $M$ is less than one. When truncating after $N$ terms, the remainder is geometrically controlled: $$(I - M)^{-1} = \left(\sum_{k=0}^{N-1} M^k\right) + M^N (I - M)^{-1}.$$ This fundamental result underpins rapid matrix inversion schemes, polynomial preconditioners, and systematic construction of bias corrections.

Explicitly, for invertible matrices,

• If $\|I-M\| < 1$, the inverse $M^{-1}$ is given by

$M^{-1} = \sum_{k=0}^\infty (I - M)^k.$

This form can be actively exploited for computational efficiency and error analysis in truncated schemes.

2. Systematic Corrections in Regression Adjustment for Randomized Experiments

In high-dimensional randomized trials with covariate adjustment, the Neumann-series decomposition yields a hierarchy of bias corrections to OLS-based average treatment effect estimators (Song, 11 Nov 2025). In the finite-population, randomization-based design, for each treatment arm $a \in \{0,1\}$, one writes the matrix inverse appearing in adjustment as

$$\Sigma_a^{-1} = (I_p - \Delta_a)^{-1} = \sum_{k=0}^\infty \Delta_a^k,$$

where $\Sigma_a$ is the sample covariance and $\Delta_a = I_p - \Sigma_a$.

The arm-wise correction is thus decomposed: $$R_a = \bar{x}_a^\top\,\Sigma_a^{-1}\,u_a = \sum_{k=0}^d \bar{x}_a^\top \Delta_a^k u_a + T_a^{(d)},$$ with $T_a^{(d)}$ the Neumann-tail remainder bounded by $\| \Delta_a \|^{d+1} / (1 - \|\Delta_a\| )$.

A degree-$d$ Neumann-corrected ATE estimator is then defined as

$$\hat{\tau}^{(d)} = \hat{\tau}_{\text{OLS}} + \sum_{k=0}^d ( \widehat{R}_1^{[k]} - \widehat{R}_0^{[k]} ),$$

where each $\widehat{R}_a^{[k]}$ is a sample analog constructed from observed residuals and Neumann weights.

This methodology yields a strict enlargement in the admissible dimensionality $p$: the estimator $\hat{\tau}^{(d)}$ is asymptotically normal under the diffuse leverage and Lindeberg-type conditions whenever

$p^{d+3} (\log p)^{d+1} = o(n^{d+2}),$

as opposed to the $p = o(n^{1/2})$ or $o(n^{2/3})$ previously established for classical regression adjustment. The remainder control tightens with higher $d$ due to geometric cancellation of inverse fluctuations up to order $d$, with estimation error from correction terms shown to be $o_P(n^{-1/2})$ for fixed $d$.

3. Neumann-Series Representations in Differential Equations

Neumann-series decompositions give explicit solution representations for linear ODEs, notably the one-dimensional Schrödinger and perturbed Bessel equations (Kravchenko et al., 2015, Kravchenko et al., 2016). Using transmutation operators and a Fourier–Legendre expansion of integral kernels, the solution to equations like

$-y''(x) + q(x) y(x) = \omega^2 y(x)$

is given as a uniformly convergent Neumann series of Bessel functions: $$c(\omega, x) = \cos \omega x + 2\sum_{m=0}^{\infty} (-1)^m \beta_{2m}(x) j_{2m}(\omega x),$$ where $j_n$ are spherical Bessel functions and the coefficients $\beta_n(x)$ are computed via stable recurrence or explicit integrals (SPPS).

Key features:

• Truncated series produce uniform-in-$\omega$ approximations with controlled error, ensuring high-accuracy spectral computations.
• Similar approaches extend to perturbed Bessel problems: $$L\,u = -u''(x) + \left( \frac{l(l+1)}{x^2} + q(x) \right) u(x) = \omega^2 u(x),$$ where the Neumann–Bessel representation yields a regular solution

$$u_l(\omega, x) = d(\omega) \sqrt{\omega x} J_{l+\frac12}(\omega x) + \sum_{n=0}^\infty (-1)^n \beta_n(x) j_{2n}(\omega x).$$

Convergence rates and uniformity follow from fine properties of the expansion kernels and the algebraic decay of coefficients.

The utility of such decompositions is manifest in the computation of vast sets of eigenvalues (“nondeteriorating accuracy”) for Sturm–Liouville spectral problems without loss for high-index modes.

4. Efficient Computation and Algorithmic Optimization

The evaluation of truncated Neumann series, particularly for large matrices, is computationally sensitive. Structured decompositions based on series factorization can substantially reduce operation counts (Dimitrov et al., 2017). Classical approaches use binary basis and Horner's method with $N-1$ multiplications for an $N$-term series. However, factorization using bases of size five (quinary) or recursively defined blocks with $y_{n+1} = y_n^2 + 1$ can reduce the cost to approximately $1.72\log_2(N) - 2$ or even $1.70\log_2(N) - 2$ multiplications in the asymptotic limit.

A representative complexity comparison:

Basis	Multiplication Exponent $\alpha$
Binary (2)	2.000
Ternary (3)	1.893
Size-5 (5)	1.722
Asymptotic	1.7016

When $N$ is a power of $5$, the size-5 basis is typically optimal. For arbitrary $N$, mixing basis strategies (e.g., size-5 and size-2) yields practical savings. These optimizations are relevant for MIMO/Massive MIMO systems and image rendering, where truncated Neumann expansions of size $N \approx 5$–$10$ are common.

5. Measurement Error Mitigation in Quantum Computation

Neumann-series decompositions facilitate measurement error mitigation on quantum devices without requiring explicit noise structure or calibration (Wang et al., 2021). Given a stochastic measurement noise matrix $M$, with $\|I - M\| < 1$, the inverted effect is approximated by truncating

$M^{-1} = \sum_{k=0}^K (I-M)^k + R_{K+1},$

where $R_{K+1}$ is the geometric remainder.

This expansion is operationalized by mapping $M^{-1}$ to a combination of sequential noisy measurements. Specifically, for measurement of an observable $O$ on quantum state $\rho$:

• Sequential reads are performed $k = 1, ..., K+1$ times, forming empirical estimates $\eta^{(k)}$.
• Coefficients $c_K(k-1) = (-1)^{k-1} \binom{K+1}{k}$ give the linear combination approximating the noise-free expectation.
• The overhead in state preparations and measurement rounds is independent of system size (qubit count), as long as noise resistance $\xi = \|I - M\|_1 < 1$.

Bias and mean-square-error bounds follow directly from the truncated remainder and concentration inequalities: $|\Tr[O \rho] - \sum_{k=1}^{K+1} c_K(k-1) \eta^{(k)}| \leq \xi^{K+1},$ with $K$ tunable to achieve target error. The method is both system-size independent and strictly model-agnostic.

6. Neumann-Series Acceleration in Iterative Solvers and Multigrid

Krylov subspace methods (e.g., GMRES) and algebraic multigrid (AMG) smoothers exploit the Neumann-series for efficient inner solves (Thomas et al., 2021). In low-synchronization modified Gram-Schmidt GMRES, the correction matrix

$T_m = (I + L_m )^{-1}$

with strictly lower-triangular $L_m$ admits a Neumann expansion,

$T_m = I - L_m + L_m^2 - L_m^3 + \ldots,$

truncated at the appropriate finite order (since $L_m$ is nilpotent for strictly lower-triangular matrices).

Replacing backward/forward triangular solves by a small number of sparse matrix-vector multiplications with $I - L_m$ or similar, one reduces synchronization and improves parallel scalability with negligible loss in convergence or stability (as shown via matrix perturbation and backward stability analyses).

In multigrid, the Gauss–Seidel smoother or ILU preconditioner leverages the Neumann-series for the inverse of the $(D + L)$ or $(I + U_s)$ factors, achieving substantial speedup (25–50$\times$ on GPU for AMG smoothing, and up to $10\times$ in strong-scaling experiments) compared to conventional triangular solves.

Structural graph reordering (e.g., AMD, symAMD) and scaling decrease the "departure from normality" of $L$ and $U$ factors, accelerating the convergence of the Neumann-series truncated inner iterations.

7. Convergence Properties, Error Control, and Scaling Laws

All practical Neumann-series decompositions rely on spectral or norm constraints for convergence and effective remainder control:

• For matrices or operators $M$, convergence is assured if $\|M\| < 1$ in the induced norm relevant for the application.
• The geometric decay of the remainder, $O(\|M\|^{K+1})$, yields explicit error bounds and guides the selection of truncation order $K$.

In regression adjustment, the analytic control achieved by Neumann truncation directly translates into improved scaling laws. For degree-$d$ correction, admissibility of $p$ covariates expands to $p^{d+3}(\log p)^{d+1} = o(n^{d+2})$ (Song, 11 Nov 2025).

In quantum measurement mitigation, the error $\leq \xi^{K+1}$ for $\xi = \|I - M\|_1$ can be tightly controlled. In ODE spectral methods, the uniform-in-$\omega$ convergence allows for computation of large spectral sets with nondeteriorating (non-polluting) accuracy.

A general theme is that geometric cancellation in the Neumann expansion systematically peels off higher-order bias or approximation error, and, provided the norm constraint is met, each additional degree of truncation propagates a relaxation in the necessary problem regularity or dimension.

---

In summary, Neumann-series decomposition provides a unifying analytic and algorithmic tool for resolving inverse operators in diverse areas—yielding explicit bias corrections, analytic structure for differential equations, efficient error mitigation strategies, and high-performance numerical algorithms with precise error and complexity guarantees. Its efficacy is demonstrated in contemporary research on regression adjustments for high-dimensional experiments (Song, 11 Nov 2025), spectral solution of ODEs (Kravchenko et al., 2015, Kravchenko et al., 2016), quantum device error mitigation (Wang et al., 2021), and numerical linear algebra (Dimitrov et al., 2017, Thomas et al., 2021).