Perturbation Bounds for Low-Rank Inverse Approximations under Noise

AI summary

Overview

Research area: Machine learning / numerical linear algebra — specifically spectral-norm perturbation theory for low-rank approximations of matrix inverses.

Technical level: Advanced. The paper relies on contour integration of matrix functions, resolvent expansions, Riesz projectors, Weyl's inequality, and classical perturbation theory.

Scope: The paper derives and empirically tests non-asymptotic spectral-norm bounds on the error made when a noise-corrupted matrix inverse is replaced by its best rank-p approximation.

What This Paper Is About

Many machine learning and optimization methods avoid inverting a large matrix exactly and instead keep only the p most important directions of the inverse. In practice the matrix is not observed exactly — it arrives with noise from sampling, sketching, quantization, or deliberate perturbation for differential privacy. This paper asks how much the resulting rank-p inverse approximation can change under such noise, and derives an explicit, spectrum-aware bound on that change in spectral norm. The authors study the quantity ‖(Ã⁻¹)_p − A_p⁻¹‖, where A is an n×n symmetric matrix, A_p⁻¹ is the best rank-p approximation of A⁻¹, and à = A + E is the noisy observation.

Key Contributions

1. A new non-asymptotic spectral-norm bound. For real symmetric positive-definite A, if 4‖E‖ ≤ min{λ_n, δ_{n−p}}, then ‖(Ã⁻¹)p − A_p⁻¹‖ ≤ 4‖E‖/λ_n² + 5‖E‖/(λ{n−p} δ_{n−p}) (Theorem 2.1). The abstract and introduction state the same result in the form ≤ 5(‖E‖/λ_n² + ‖E‖/(λ_{n−p} δ_{n−p})).

2. A contour-bootstrapping technique for a non-entire matrix function. The authors apply contour-integral methods to f(z) = 1/z, which is not entire, localizing the resolvent expansion around the smallest eigenvalues and controlling the perturbation of the associated spectral projectors. This yields the bootstrapping lemma F ≤ 2F₁ (Lemma 3.1).

3. Improvement over classical baselines. The bound improves on naive adaptations of classical full-inverse bounds by up to a factor of √n, and provides a quantitative criterion — depending on the eigengap, spectral decay, and alignment of E with low-curvature directions — for when a low-rank inverse remains reliable. The paper states this is the first work to give non-asymptotic spectral-norm bounds for ‖(Ã⁻¹)_p − A_p⁻¹‖ under general additive noise.

4. Extension beyond positive-definite matrices and a concrete downstream application. Theorem 4.1 extends the result to arbitrary real symmetric A (with eigenvalues λ₁ > … > λ_k > 0 > λ_{k+1} > … > λ_n), and Section B applies the framework to the preconditioned conjugate-gradient (PCG) method, where it yields tighter condition-number estimates and a provable n^{1/4} improvement in the guaranteed iteration count relative to the classical Eckart–Young–Mirsky–Neumann-based analysis. Remark 4.7 notes the approach also applies to rank-deficient matrices via their pseudoinverses.

Main Findings

• The bound splits into two interpretable terms. The first term ‖E‖/λ_n² reflects classical perturbation scaling for the full inverse; the second term ‖E‖/(λ_{n−p} δ_{n−p}) captures the extra sensitivity introduced by projecting onto the subspace spanned by the smallest eigenvalues of A.

• The gap condition matters. The spectral condition 4‖E‖ ≤ λ_n guarantees à is invertible; the gap assumption 4‖E‖ ≤ δ_{n−p} guarantees by Weyl's inequality that the eigenvectors of the smallest p eigenvalues of à remain aligned with those of A. When δ_{n−p} ≪ ‖E‖, eigenvalues of à can reorder and the low-rank approximation becomes unstable (Section F.1 gives an example).

• Explicit comparison to the classical baseline. The Eckart–Young–Mirsky–Neumann (EYM–N) bound is ‖(Ã⁻¹)p − A_p⁻¹‖ ≤ 2(‖E′‖ + λ{n−p}⁻¹) ≤ 8‖E‖/(3λ_n²) + 2/λ_{n−p}, valid when 4‖E‖ ≤ λ_n. That bound degrades when λ_{n−p} ≪ λ_n²/‖E‖ and fails to capture the limit ‖(Ã⁻¹)_p − A_p⁻¹‖ → 0 as ‖E‖ → 0.

• Where the two bounds coincide and where the new one wins. The EYM–N bound and the bootstrapped bound coincide in order of magnitude when ‖E‖ ≫ λ_n²/λ_{n−p}. When ‖E‖ ≪ λ_n²/λ_{n−p}, the new bound is smaller by a factor of min{λ_n²/(λ_{n−p}‖E‖), δ_{n−p}/‖E‖}. This gain regime arises when p < sr(A⁻¹) := Σ_{i=1}^n λ_n/λ_i, or when δ_{n−p}λ_{n−p} ≪ λ_n².

• Worked gain example. For a matrix with spectrum {n, 2n, …, 10n, 20n, 20n, …, 20n} and p = 10, with standard Gaussian noise, the EYM–N bound evaluates to O(√n/n² + 1/n) = O(1/n) while the new bound gives O(√n/n²) = O(n^{−3/2}) — the expected √n-level gain.

• Random-noise specialization. If E is a Wigner matrix (symmetric with i.i.d. sub-Gaussian entries), then ‖E‖ = (2 + o(1))√n with high probability, giving ‖(Ã⁻¹)p − A_p⁻¹‖ = O(√n/λ_n² + √n/(λ{n−p}δ_{n−p})), versus O(√n/λ_n² + 1/λ_{n−p}) for the EYM–N bound, which is larger when √n ≪ δ_{n−p}.

• Recovery of the full-inverse case. Setting p = n gives δ_{n−p} = λ₀ = +∞, so the second term vanishes and the bound reduces to Θ(‖E‖/λ_n²), matching the Neumann bound.

• Empirical agreement. On real and synthetic matrices — sample covariance matrices, discretized elliptic operators, and sparse structural stiffness matrices such as BCSSTK09 — the bound tracks the true error within a small constant factor (Section 7.2 reports typically below 10), while classical estimates often overpredict by one to two orders of magnitude.

• Assumption admissibility. Section 7.1 evaluates the condition 4‖E‖ < min{λ_n, δ_{n−p}} on the 1990 US Census covariance and BCSSTK09 stiffness matrices, finding safety margins that comfortably exceed noise levels common in differential-privacy and structural-engineering applications. For a sample covariance A = MᵀM with M ∈ ℝ^{m×n}, m ≥ n, and E symmetric with i.i.d. zero-mean sub-Gaussian entries of variance Δ², both conditions hold with high probability if ‖M‖_F² ≥ m log n and m > Cn^{3/2}Δ/log n for some constant C > 0.

• Robustness to misestimation. The factor 1/2 in the contour construction may be replaced by any fixed constant c ∈ (0,1) by adjusting the contour, changing the estimate only up to a constant. As long as estimation errors in λ_n and δ_{n−p} are within ‖E‖, the bound remains valid up to a constant factor.

• A stronger but unevaluated bound. In the intermediate regime λ_n²/λ_{n−p} ≪ ‖E‖ ≪ min{δ_{n−p}, λ_n}, the bound in Theorem 5.2 offers an asymptotic improvement over both Theorem 2.1 and EYM–N, but depends on additional structural quantities such as the alignment of E with the low-curvature eigenspace and is not empirically evaluated.

• Extension result. Theorem 4.1 states that if 4‖E‖ < min{δ_{k−k₁}, δ_{k+p−k₁}, σ_n} and σ_{n−p} − σ_{n−p+1} > 2‖E‖, then the error is bounded by 4‖E‖/λ_k² + 5‖E‖/(λ_{k−k₁}δ_{k−k₁}) + 4‖E‖/|λ_{k+1}|² + 5‖E‖/(|λ_{k+p−k₁+1}|δ_{k+p−k₁}).

• Why the non-entire case is hard. For deformed Wigner perturbations, ‖Ã⁻¹‖ = O(n) with high probability for any fixed real A, so the smallest singular values of A are effectively destroyed by noise. The paper also reports (Section F.2) a counterexample showing that a bound on ‖Ã_p − A_p‖ does not control the inverse approximations.

Methodology in Plain English

The authors express the rank-p inverse approximation as a contour integral. They draw a closed curve in the complex plane that surrounds the p smallest eigenvalues of A, deliberately avoiding the singularity of 1/z at zero and excluding the larger eigenvalues. Integrating z⁻¹(zI − A)⁻¹ around this curve reproduces exactly the best rank-p approximation of A⁻¹; doing the same with à produces the noisy counterpart.

The difference between the two integrals is bounded by an integral over the curve, which makes the interaction between A and the noise E explicit rather than hiding it. The problem is that the perturbed resolvent (zI − Ã)⁻¹ is awkward to handle directly. The authors expand it as an infinite series in E and show, through a bootstrapping argument, that the whole series is at most twice its first term — a relationship previously known for entire functions such as exp(z) and the constant function 1, which the authors extend here to f(z) = 1/z under the condition 4‖E‖ ≤ min{λ_n, δ_{n−p}}.

What remains is to choose the contour so the resulting integral is finite and computable. They use a rectangle with vertical sides at x₀ = λ_n/2 and x₁ = λ_{n−p+1} + δ_{n−p}/2 and height T = 2λ₁. Each point on the contour sits at least half of δ_{n−p} or λ_n away from the spectrum of A. Evaluating the four sides gives the two main terms, one scaling as ‖E‖/λ_n² from the left vertical side and one scaling as ‖E‖/(λ_{n−p}δ_{n−p}) from the right vertical side. Unlike prior contour arguments where only one side dominates, here both vertical sides contribute essentially equally, while the horizontal sides contribute a smaller term.

Why This Matters

Impact on research. Classical matrix perturbation theory addresses the full inverse and ignores truncation to rank p; results on low-rank perturbations under Schatten or Frobenius norms do not provide spectral-norm guarantees for pseudoinverses under general noise; and the Sherman–Morrison–Woodbury formula covers only structured low-rank updates. This paper fills that gap with a spectrum-aware guarantee that becomes increasingly sharp as noise decreases, and it introduces contour bootstrapping for the non-entire function 1/z as a reusable technique.

Real-world applications (as identified in the paper):

• Kernel methods, Gaussian processes, covariance-based inference, and solvers for structured systems and graph Laplacians, all of which use inverse approximations.
• Fast solvers, adaptive preconditioners, and randomized linear algebra, including the PCG application in Section B.
• Differentially private computation, where noise is added deliberately, and where the paper shows the assumptions hold with comfortable margin.
• Structural engineering, where stiffness matrices such as BCSSTK09 are used and where noise levels stay below the derived thresholds.

Industry relevance. Practitioners who use sketching, quantization, or privacy mechanisms to build approximate second-order information get a quantitative criterion for deciding when a rank-p inverse stays trustworthy. The n^{1/4} improvement in guaranteed PCG iteration count translates into a concrete statement about solver effort, and the paper notes that exact verification of the assumptions is often unnecessary because the bound degrades only by a constant factor under moderate misestimation of λ_n and δ_{n−p}.

Future Directions

• Empirically evaluating the refined bound. Remark 2.2 notes that the Theorem 5.2 bound may be tighter in the moderate-noise, slow-spectral-decay regime but states explicitly that it is not empirically evaluated; testing it against the true error would clarify when the extra complexity pays off.

• Completing the general symmetric extension. The provided content for Theorem 4.1 is truncated mid-sentence — the issue of the perturbed eigenvalues not corresponding to the p least singular values of à is described as being resolved by "enforcing the sing…". The full treatment of rank-deficient and indefinite cases via pseudoinverses (Remark 4.7) is a natural area for further development.

• Extending contour bootstrapping to other non-entire functions. The paper argues that the technique was previously applied mainly to entire functions such as exp(z) and to f(z) = 1 for eigenspace projections. Other rational or non-entire matrix functions used in machine learning could be analyzed with the same machinery.

• Relaxing the noise model. The strongest stated result covers symmetric additive noise E, with a Wigner specialization where ‖E‖ = (2 + o(1))√n. Whether the analysis extends to non-symmetric, structured, sparse, or adversarial perturbations is not reported.

Target Audience

Researchers and graduate students in numerical linear algebra, randomized numerical linear algebra, and theoretical machine learning who work on matrix perturbation theory, preconditioning, kernel methods, Gaussian processes, or private and approximate optimization. Also relevant to practitioners who build large-scale solvers and need a certificate for when a low-rank inverse approximation remains reliable under noise. Readers without a background in complex analysis and matrix perturbation theory will find the technical core difficult; the empirical comparisons and the statement of Theorem 2.1 are the most accessible entry points.

Authors’ abstract

Low-rank pseudoinverses are widely used to approximate matrix inverses in scalable machine learning, optimization, and scientific computing. However, real-world matrices are often observed with noise, arising from sampling, sketching, and quantization. The spectral-norm robustness of low-rank inverse approximations remains poorly understood. We systematically study the spectral-norm error $\| (\tilde{A}^{-1})_p - A_p^{-1} \|$ for an $n\times n$ symmetric matrix $A$, where $A_p^{-1}$ denotes the best rank-\(p\) approximation of $A^{-1}$, and $\tilde{A} = A + E$ is a noisy observation. Under mild assumptions on the noise, we derive sharp non-asymptotic perturbation bounds that reveal how the error scales with the eigengap, spectral decay, and noise alignment with low-curvature directions of $A$. Our analysis introduces a novel application of contour integral techniques to the \emph{non-entire} function $f(z) = 1/z$, yielding bounds that improve over naive adaptations of classical full-inverse bounds by up to a factor of $\sqrt{n}$. Empirically, our bounds closely track the true perturbation error across a variety of real-world and synthetic matrices, while estimates based on classical results tend to significantly overpredict. These findings offer practical, spectrum-aware guarantees for low-rank inverse approximations in noisy computational environments.