Perturbation Bounds for Low-Rank Inverse Approximations under Noise

Thank you for your valuable comments and suggestions. We are glad that you recognize our improvements over classical results. We address your specific questions and comments below.

Comments:

Though, for instances where gap can be of the order $O(\log n)$ , the bound obtained seems worse compared to the EY-N bound.

Thank you for the comment. We fully acknowledge that our bounds do not apply in regimes where $\delta_{n-p} < \\|E\\|$. In such cases—particularly when the gap is as small as $O(\log n) \ll \Theta(\sqrt{n}) = \\|E\\|$—our bound may indeed be weaker than the classical EY–N bound.

We would also like to note that when $\delta_{n-p}$ is too small relative to the noise level $\\|E\\|$, the low-rank inverse $A_p^{-1}$ itself may cease to be a stable or meaningful approximation to $A^{-1}$. This behavior is illustrated in Appendix J.1, where we provide a concrete example showing that the spectral instability of small eigenvalues fundamentally limits robustness in this regime.

We will revise the manuscript to clarify this limitation more explicitly.

For justifying the assumption of $||E|| \ll \delta_{n-p}$ , the authors provide discussion assuming sub-Gaussian noise, on the variance of noise that ensures this condition w.h.p which is fine as long as $m \gg n$ , otherwise the permissible noise scale can be much smaller - $\tilde{O}(1/\sqrt{n})$ . This also potentially limits the applicability of the bound for high-dimensional data, which is typically the case for modern machine learning applications.

Thank you for this helpful comment. We apologize if our presentation was unclear, and we appreciate the opportunity to clarify.

In the example discussed, we implicitly assume $m \geq n$ so that $A^{-1}$ and $A_p^{-1}$ are well-defined. Otherwise, when $m < n$, the matrix $A = M^\top M$ becomes rank-deficient and both $\lambda_n$ and $\lambda_{n-p}$ are zero for small $p$.

In such cases, one should instead consider the pseudoinverses $A^\dagger$ and $A_p^\dagger$. Our framework and techniques still apply, with $\lambda_m$ and $\delta_{m-p}$ replacing $\lambda_n$ and $\delta_{n-p}$, and the roles of $m$ and $n$ adjusted accordingly. In particular, under sub-Gaussian noise, when $\\|M\\|_F^2 \gg n \log m$, the permissible noise scale becomes $\tilde{O}(1/\sqrt{m})$—which is often reasonable even in high-dimensional settings.

More broadly, we note that this noise tolerance is especially relevant in sparse or structured data scenarios, where the effective energy of $M$ is concentrated. When entries of $M$ are $\Theta(1)$ and $M$ is dense, one should expect a significantly larger permissible noise scale. The assumption remains meaningful and the bound applicable.

We will include a clarification of this point in the final version.

Thank you for your valuable comments. We are glad that you recognize the importance of the fundamental problem we study and appreciate the improvements our work provides over classical results. We address your specific remarks below.

Comments:

In the main theoretical results, besides that the noise matrix bounded by the spectrum, the matrix $A$ is required to be real symmetric and positive semi-definite, and the noise matrix $E$ is required to be symmetric. These requirements will also limit the application of the results.

Thank you for the comment. We apologize if our presentation caused any confusion. We would like to clarify that the main result in Theorem 2.1 assumes $A$ is symmetric positive definite, but the extensions provided in Appendices A and B apply to general symmetric matrices, without requiring $A$ to be positive semidefinite.

Regarding the boundedness of the noise, the condition that $||E||$ is small relative to the spectrum of $A$ is natural and necessary: without such control, the perturbation can overwhelm the small eigenvalues or singular values of $A$, making any spectral-norm guarantee for the inverse ill-posed.

We also note that our approach extends naturally to the case where $A$ is symmetric but rank-deficient (i.e., $\operatorname{rank}(A) = r < n$). In that setting, we work with the pseudoinverse $A^\dagger$ and project onto the subspace corresponding to the nonzero eigenvalues $\tilde\lambda_r, \tilde\lambda_{r-1}, \dots, \tilde\lambda_{r-p+1}$. The contour $\Gamma$ is then reconstructed accordingly using $(\lambda_r, \lambda_{r-p+1}, \delta_{r-p})$.

We will revise the paper to more clearly state these extensions and their assumptions in the final version.

Only several specific matrices are checked in the empirical study, including the 69-by-69 1990 US Census covariance matrix, the 1083-by-1083 BCSSTK09 stiffness matrix, and synthetic matrices from discretised Hamiltonians. Some examples with real-world applications will better demonstrate the value of this work.

Thank you for this helpful comment. We agree that additional real-world application examples could further illustrate the value of our results.

In this paper, our primary goal was to introduce the problem of low-rank inverse approximation under noise, address its key theoretical challenges, and present a clean and general analysis framework. We deliberately selected datasets that are diverse in structure and origin—including real-world matrices from scientific computing and statistics—to provide meaningful validation while keeping the focus on the underlying mathematical phenomena.

That said, we fully agree that exploring richer applications (e.g., in private optimization, kernel methods, or compressed solvers) is a natural next step, and we plan to pursue these directions in future work. We will revise the manuscript to clarify this motivation and scope.

Thank you for your valuable comments and questions. We are glad that you recognize our contributions to perturbation theory for low-rank pseudoinverses and appreciate your interest in the potential connections between our theoretical results and practical applications. We address your specific points below.

Questions:

Can you propose a low-cost procedure (e.g., randomized sketch) to approximate $\delta_{n-p}$ and $\lambda_n$ sufficiently well to certify the gap condition without full eigendecomposition?
[Related comment] The requirement $4\\|E\\| < \delta_{n-p}$ may be hard to check in black-box or data-streaming settings, as estimating tail eigengaps can be as expensive as the inverse itself.

Thank you for this thoughtful question. We agree that efficiently estimating $\delta_{n-p}$ and $\lambda_n$ is important for practical use, particularly in large-scale, streaming, or black-box settings where full eigendecomposition is infeasible.

Fortunately, randomized Krylov subspace or Lanczos methods offer scalable ways to approximate extremal eigenvalues. For symmetric PSD matrices, applying Lanczos with a modest number of iterations can yield reliable estimates of $\lambda_n$ and $\lambda_{n-p}$, and hence of the eigengap $\delta_{n-p} = \lambda_{n-p} - \lambda_{n-p+1}$. These methods require only matrix-vector products and scale as $O(km)$, where $k$ is the number of iterations and $m$ is the cost of a matrix-vector multiplication.

We also note that exact evaluation of the gap condition is not always necessary. In many practical settings—especially when the spectrum exhibits decay or separation—coarse estimates or structural priors (e.g., from known physical models or sketch-based spectrum summaries) are sufficient to verify the condition approximately. For example, it is often possible to use a small number of random projections or polynomial filters to check whether the tail spectrum is well-separated without resolving every eigenvalue.

We will add a brief discussion of these practical considerations in the final version.

How sensitive is the bound to slight misestimation of $\delta_{n-p}$ or $\lambda_1$ when constructing $\Gamma$? [Related comment] While theoretically elegant, the custom rectangular contour $\Gamma$ (Section 3) depends on spectral extremal values $(\lambda_1, \lambda_n, \delta_{n-p})$; practical implementation may require careful numerical tuning and stability analysis.

Thank you for this insightful question. As you noted, our bound and the construction of $\Gamma$ depend on controlling $(\lambda_n, \lambda_{n-p}, \delta_{n-p})$, and we agree that understanding the impact of misestimation is an important direction.

Our method is robust to moderate misestimation: as long as the error in estimating $\lambda_n$ or $\delta_{n-p}$ is no greater than $\\|E\\|$, the bound remains valid up to a constant factor. The key condition is that the contour $\Gamma$ must enclose exactly the eigenvalues $\tilde\lambda_n, \tilde\lambda_{n-1}, \dots, \tilde\lambda_{n-p+1}$ while staying sufficiently away from all perturbed eigenvalues. This can be ensured by adjusting $\Gamma$—e.g., shifting its left edge closer to 0 or expanding the right edge toward $\lambda_{n-p}$. Specifically, it suffices that the distance from $\Gamma$ to $\lambda_n$ is at least $\min \\{c \lambda_n, \\|E\\| \\}$ and to $\lambda_{n-p}$ or $\lambda_{n-p+1}$ is at least $\min \\{c \delta_{n-p}, \\|E\\| \\}$ for some constant $0 < c < 1$.

If the misestimation greatly exceeds $\\|E\\|$, however, key eigenvalues such as $\tilde\lambda_n$ or $\tilde\lambda_{n-p+1}$ may approach $0$ or $\lambda_{n-p}$ respectively, forcing $\Gamma$ too close to these points. This can cause the bound to deteriorate due to a potential blow-up in $\\|(zI - A)^{-1}\\|$.

We will consider including a brief discussion of this stability condition in the final version and agree that further analysis of numerical tuning and adaptive contour selection is a promising direction for future work.

Section A mentions general symmetric $A$; can your approach handle non-square or non-Hermitian matrices (e.g., via pseudospectrum)?

Thank you for this thoughtful question. Yes, our approach can be extended beyond the symmetric full-rank setting in two directions:

1. Symmetric but rank-deficient matrices:
   Our method extends naturally to symmetric matrices $A$ of rank $r < n$. In this setting, we use the pseudoinverse $A^\dagger$ and work within the subspace corresponding to the nonzero eigenvalues $\lambda_r, \lambda_{r-1}, \dots, \lambda_{r-p+1}$. The contour $\Gamma$ is then redefined in terms of $(\lambda_r, \lambda_{r-p+1}, \delta_{r-p})$, and the same analysis applies by projecting $\tilde{A}^{-1}$ accordingly.

2. Non-square or non-Hermitian matrices:
   For general rectangular or non-Hermitian matrices $A$, one can embed $A$ into a symmetric matrix

   $$A' = \begin{bmatrix} 0 & A \\\ A^\top & 0 \end{bmatrix},$$

   which is symmetric (and Hermitian if $A$ is complex) but typically rank-deficient. Our contour-based framework then applies to $A'$, and one can recover information about $A^\dagger$ via its singular value decomposition. This connection may also be viewed through the lens of the pseudospectrum, and we believe further exploration in this direction would be interesting.

We will add a brief note about these extensions in the final version.

Is there a practical way to measure or bound $\\|E\\|$ within the low-curvature eigenspace to tighten the bound further in specific applications (e.g., quantization noise)?

This is an excellent observation. In many applications—particularly those involving quantization or structured noise—the perturbation $E$ tends to have low energy in the low-curvature (i.e., small-eigenvalue) directions of $A$. This suggests that a refined bound depending on $\\|\Pi_{>p} E \Pi_{>p}\\|$, where $\Pi_{>p}$ denotes the projection onto the tail eigenspace, could offer a tighter and more informative alternative to the global spectral norm bound. In practice, this quantity can often be estimated using a small number of random projections or empirical measurements when the noise model is known (e.g., isotropic or bounded quantization noise). We agree that incorporating such direction-sensitive refinements is a promising direction and will include a brief discussion in the future work section of the final version.

Other comments:

The gap $\delta_{n-p}$ is sometimes denoted $\delta_{n-p}$ and elsewhere $\delta_{n-p+1}$ — please standardize.

Thank you for pointing out this inconsistency. We will carefully review our notation and ensure that it is used consistently throughout the paper in the final version.

The analysis assumes a fixed matrix $A$ ; extensions to time-varying or adaptive contexts (e.g., iterative solvers) are not addressed.

Thank you for the insightful comment. You are correct that our current analysis focuses on a fixed matrix $A$ and does not directly address time-varying or adaptive settings. Such extensions would indeed require leveraging additional structure specific to each application—for example, spectral stability over iterations or control over update noise in optimization routines.

That said, we believe the tools developed in this work—particularly our contour-based approach and spectrum-sensitive bounds—could serve as a foundation for analyzing robustness in iterative solvers and other adaptive contexts. We view this as a promising direction for future research and will include a note to this effect in the final version.

Thank you for the valuable comments. We address your specific points below.

Comments:

The authors basically prove a sharper perturbation bound using techniques in complex analysis. [...] the new bound on the error takes into account the spectral gap, spectral decay, so it is naturally "sharper", because it is more specific. This is also my view of the experimental section, which shows that the author's bound is tighter than classical ones.

We appreciate this perspective. While we fully acknowledge that the contour integral representation is a well-established tool in functional perturbation theory, our contribution lies in how we extend and apply this tool in a new regime.

To the best of our knowledge, our work is the first to (1) apply the contour method to a non-analytic function like $f(x) = 1/x$ (which lacks a convergent power series around $x = 0$), and (2) develop a contour bootstrapping argument that directly targets the low-curvature part of the spectrum — where classical bounds are typically loose. These aspects are central to our technique and not found in prior literature.

Moreover, the improvement over classical bounds (such as those derived from the Eckart–Young–Neumann theorem or Neumann expansions) can reach up to a $\sqrt{n}$ factor in relevant regimes. We believe this is a meaningful gain, particularly for applications in high-dimensional settings or with structured noise.

That said, we appreciate that the novelty may not have been sufficiently emphasized in the current version, and we will revise the introduction and related work sections to highlight these contributions more clearly.

[...] a purely theoretical result for numerical linear algebra. [...] the baseline result here can be derived in just 3 or 4 lines, based on very classic ideas. [...] But on the other hand, it is also hard to view this as a machine learning paper, because the lack of direct connections to applications.

We appreciate the comment and the concerns about positioning. However, we believe the contribution is both conceptually and technically novel, and relevant for ML-adjacent problems, as we outline below.

While classical techniques can yield rough upper bounds in a few lines, they do not offer sharp, spectrum-aware guarantees in the spectral norm—particularly in the regime where the inverse is most sensitive (i.e., near small eigenvalues). To our knowledge, no prior work provides non-asymptotic, low-rank inverse perturbation bounds of the kind we establish.

In particular, the problem of approximating $A_p^{-1}$ via $\tilde{A}_p^{-1}$ presents several distinct and nontrivial challenges:

1. Controlling the smallest singular values is delicate, as they are highly sensitive to noise. Even small Gaussian perturbations can significantly shift or destroy the least singular value. For instance, when $E$ is a standard GOE, it is well-known that $||\tilde{A}^{-1}|| = O(n)$ with high probability for any fixed real $n \times n$ matrix $A$.

2. The error $||\tilde{A}_p^{-1} - A_p^{-1}||$ is not scale-invariant, unlike many standard perturbation measures. Determining tolerable noise levels while preserving utility is a meaningful and technically subtle question, especially in systems involving sketching, quantization, or privacy.

3. The transition from spectral gaps in $A$ to bounds on the inverse error is analytically subtle, particularly near the tail of the spectrum, where eigengaps are narrow and easily erased by noise.

Our contour bootstrapping technique was developed to address these issues directly. To the best of our knowledge, this is the first work to provide sharp, non-asymptotic spectral-norm bounds for low-rank inverse approximations—with up to a $\sqrt{n}$ improvement over classical estimates in structured regimes.

As for application relevance: although the paper focuses on foundational results, these are directly motivated by practical challenges in ML pipelines—such as preconditioned solvers, private optimization, and low-rank kernel methods. We will add a brief discussion of such connections in the final version to make this link more explicit.

We focused this paper on resolving the core theoretical question clearly and precisely, as this problem had not previously been addressed at this level of generality. We believe this kind of rigorous, broadly applicable contribution is well aligned with NeurIPS.

Thank you for your valuable comments and suggestions. We are glad that you appreciate our novel technique related to the smallest eigenvalues and that you find interested in some potential extensions of our work. Below, we address your specific question and concerns.

Questions:

Equation (2) is derived in Section F. Upon checking, I noticed that the authors are using Weyl's inequality to bound the eigengap between $\\| \tilde A^{-1}-A^{-1} \\|$ as $\\|\tilde A^{-1} -A^{-1}\\| \ge\tilde\lambda_{n-p}^{-1} - \lambda_{n-p}^{-1}$ . However, as far as I know Weyl's is $|\lambda_i^{A} -\lambda_{i}^B| \leq \\|A -B\\|$. Are the authors using some specific Weyl's inequality? If yes, please cite it accordingly, or if there is any assumption to this statement, it should be stated in the theorem statement and not assumed in the proofs.

Thank you for this question and the careful reading. We apologize if the presentation was unclear and appreciate the opportunity to clarify. In Section F, we use the standard version of Weyl’s inequality as you stated. The inequality $\\|\tilde A^{-1} -A^{-1}\\| \geq |\tilde\lambda_{n-p}^{-1} - \lambda_{n-p}^{-1}|$ is a direct consequence of applying Weyl's inequality to the matrices $\tilde{A}^{-1}$ and $A^{-1}$.

Since $4\\|E\\| < \lambda_n$ and $A$ is positive definite, $\tilde{A}$ remains positive definite as well. Therefore, $\tilde\lambda_{n-p}^{-1}$ and $\lambda_{n-p}^{-1}$ are the $(p+1)^{th}$ eigenvalues of $\tilde{A}^{-1}$ and $A^{-1}$, and Weyl's inequality applies. We will add this clarification in the proof of Section F in the final version.

While the authors state the bounds are for PSDs in the main paper, the first line in Setup (line 95) assumes for $A \in \mathbb{R}^{n \times n}$ that $\lambda_1 \geq \lambda_2 \geq \dots \geq \lambda_n > 0$, which defines a PD matrix.

Thank you for pointing this out. You are correct—the setup assumes that $A$ is positive definite. Indeed, when $A$ is PSD, its inverse $A^{-1}$ is well-defined only when $\lambda_n > 0$. We will clarify this in the final version.

We would also like to note that our approach naturally extends to the case where $A$ is symmetric PSD of rank $r < n$. In this setting, we replace $A^{-1}, A_p^{-1}$, and $\tilde A_p^{-1}$ with the pseudoinverse of $A$ and the projection of $\tilde A^{\dagger}$ onto the subspace corresponding to the nonzero eigenvalues $\tilde\lambda_r, \tilde\lambda_{r-1}, \dots, \tilde\lambda_{r-p+1}$. The contour $\Gamma$ can then be constructed with respect to $(\lambda_r, \lambda_{r-p+1}, \delta_{r-p})$, and the analysis proceeds similarly.

While $\tilde{A}$ is assumed to be symmetric, in Step 1 of the Proof Overview, its eigenvalues are also assumed to be positive. This was confusing—can $E$ not perturb the eigenvalues of $A$ to make $\tilde{A}$ indefinite?

Thank you for this helpful question—we appreciate the opportunity to clarify. The assumption in Theorem 2.1 implicitly ensures that $\tilde{A}$ remains positive definite, though we should have stated this more explicitly.

Indeed, under the condition $4\\|E\\| < \lambda_n$ from Theorem 2.1, Weyl’s inequality guarantees that all eigenvalues of $\tilde{A}$ are at least $\lambda_n - \\|E\\| > 3\\|E\\| > 0$, which ensures that $\tilde{A}$ remains positive definite. We will include this clarification in the final version of the paper.

In the same line in Step 1 of Proof overview, the authors use this nice trick to swap the norm and the integrals. However $dz$ should be $|dz|$ , which complicates the whole proof (I am unsure if the proof breaks though). The same issue exists in appendix . Is $dz$ and $|dz|$ interchangable in this context? Surely not, right, as one has to hand two cases $\geq 0$ and $< 0$ separately?

Thank you for catching this—we appreciate the opportunity to clarify and correct the notation.

You are absolutely right: the correct form of the inequality should use $|dz|$, in line with the standard bound

In Step 1, we mistakenly wrote $dz$ instead of $|dz|$. Starting from line 178, the derivation should indeed switch to $|dz|$, and this is what we intended throughout. Later, in Step 3, where the rectangular contour $\Gamma$ is explicitly constructed, we estimate each integral segment using $dt$ in place of $|dz|$, as appropriate for real-line integrals.

We will correct the notation in both the main text and the appendix in the final version. Thank you again for pointing this out.

In line 803, the authors state that there exists some permutation of $\pi$ such that $\sigma_{n-p}=\lambda_{\pi(p)}$ . The $\lambda$ here should be an absolute value, as otherwise, yes there is some perturbation when this is true, but such perturbation might not lead to consecutive values be $> 2||E||$ guaranteed to be $> 2||E||$, which is a required assumption in Theorem A.1

Thank you for the question. As you correctly noted, the proper statement is $\sigma_{n-p} = |\lambda_{\pi(p)}|$. We apologize for the typo and any confusion it caused.

We will fix it in the final version.

Other comments:

Given those are resolved, I would like to understand what the authors think this result would mean for problems like the Schatten-1 norm or spectral density estimation under noisy environments.

Thank you for the insightful question and suggestion. We are happy to discuss possible implications of our work in these directions.

1. Schatten-1 Norm Error: Suppose the goal is to estimate $\\| \tilde A_p^{-1} - A_p^{-1}\\|^{S_1}$. Our contour-based approach can be extended to this setting. Note that \\|\tilde A_p^{-1} - A_p^{-1}\\|^{S_1} = \left| \sum_{i=1}^n e_i^\top (\tilde A_p^{-1} - A_p^{-1}) e_i \right| = \left| \int_{\Gamma} \frac{1}{2\pi i z} \sum_{i=1}^n e_i^\top \left[ (zI - \tilde{A})^{-1} - (zI - A)^{-1} \right] e_i \hspace{1mm} dz \right|.$$ We can apply our contour bootstrapping technique to this integral, reducing the problem to bounding a term of the form F_1 = \int_{\Gamma} \left| \frac{1}{z} \sum_{i=1}^n e_i^\top (zI - A)^{-1} E (zI - A)^{-1} e_i \right| \hspace{1mm} |dz|. $$While obtaining a sharp bound for this trace term may require additional assumptions or techniques, our framework provides a natural starting point for extending spectral-norm results to Schatten-1 norms. We would be happy to explore this further if there is a particular application or setting you have in mind.$$
2. Spectral Density Estimation: We are not yet certain which specific formulation or model for spectral density estimation under noise you are referring to. However, our contour-based analysis and sensitivity bounds could potentially inform error guarantees in estimators that rely on smoothed spectral traces (e.g., via Chebyshev expansions or resolvent averaging). We would welcome the opportunity to discuss this direction further.