Below, we compute with exact real numbers using a realistic / conservative model of computability like TTE.

Assume that there is an algorithm that, given a symmetric real matrix $M$, finds real matrices $P$ and $D$, where $P$ is orthogonal and $D$ is diagonal, such that $M = PDP^T$. I will show that no such algorithm exists:


$$M(\epsilon) = \begin{cases} \left[\begin{matrix}1 & \epsilon\\ \epsilon & 1\end{matrix}\right] ,& \epsilon \geq 0 \\ \left[\begin{matrix}1 - \epsilon & 0\\0 & 1 + \epsilon\end{matrix}\right] ,& \epsilon \leq 0 \end{cases}$$

Define $P(\epsilon)$ and $D(\epsilon)$ analogously.

  • If $\epsilon > 0$ then $P(\epsilon)$ must necessarily be $\frac 1 {\sqrt 2}\left[\begin{matrix}1 & -1\\1 & 1\end{matrix}\right]$.

  • If $\epsilon < 0$, then $P(\epsilon)$ must necessarily be $\left[\begin{matrix}1 & 0\\0 & 1\end{matrix}\right]$.

We know that $P(\epsilon)$ can't be the first matrix if the sum of the absolute values of its entries is less than $2 \sqrt 2 \approx 2.82$; and we know that it can't be the second matrix if the sum of the absolute values of its entries is greater than $2$. It's easily determinable whether the sum of the absolute values of its entries is less than $2 \sqrt 2$ or greater than $2$. Therefore we can decide whether $\epsilon \leq 0$ or $\epsilon \geq 0$, which is impossible!

In a way, this is strange, because many sources say that the Singular Value Decomposition (SVD) and Schur Decomposition (which are generalisations of the Spectral Decomposition) are numerically stable. They're also widely used in numerical applications. But I've just tested the examples above for small $\epsilon$ using SciPy and got incorrect results.

So my question is, how do numerical analysts get around this problem? Or why is this apparently not a problem?

I could venture some guesses: While finding eigenvectors of general matrices may be impossible, it is possible to find their eigenvalues. Also, it's possible to "shift" a problematic matrix by some small $\epsilon$ so that its eigendecomposition is computable.

  • 2
    $\begingroup$ I feel that the first sentence is a bit vague: "Assume that there is an algorthim that..." but you don't say in which computation model you place the question. This seems to me of paramount importance for the rest of the question (BSS will give you a different answer compared to the case where $\epsilon$ must be written down in a real-world computer memory)(also the 0 test is computable in floating-point arithmetic in real-world computers). $\endgroup$ – Loïc Teyssier Aug 24 at 14:33
  • 3
    $\begingroup$ @ogogmad, the word "realistic" doesn't help $\endgroup$ – Matt F. Aug 24 at 14:35
  • 2
    $\begingroup$ @LoïcTeyssier It is. It's a continuous function $\endgroup$ – ogogmad Aug 24 at 14:57
  • 9
    $\begingroup$ @LoïcTeyssier $M(\epsilon) = \left[\begin{matrix}\min\left(0, \epsilon\right) + 1 & \max\left(0, \epsilon\right)\\\max\left(0, \epsilon\right) & 1 - \min\left(0, \epsilon\right)\end{matrix}\right] $, which is clearly computable $\endgroup$ – ogogmad Aug 24 at 15:00
  • 3
    $\begingroup$ I stand corrected ;) $\endgroup$ – Loïc Teyssier Aug 24 at 15:02

The singular value decomposition, when applied to a real symmetric matrix $A = \sum_i \lambda_i(A) u_i(A) u_i(A)^T$, computes a stable mathematical object (spectral measure $\mu_A = \sum_i \delta_{\lambda_i(A)} u_i(A) u_i(A)^T$, which is a projection-valued measure) using a partially unstable coordinate system (the eigenvalues $\lambda_i(A)$ and eigenvectors $u_i(A)$; the eigenvalues are stable, but the eigenvectors are not). The numerical instability of the latter reflects the coordinate singularities of this coordinate system, but does not contradict the stability of the former. But in numerical computations we have to use the latter rather than the former, because standard computer languages have built-in data representations for numbers and vectors, but usually do not have built-in data representations for projection-valued measures.

An analogy is with floating-point arithmetic. The operation of multiplication of two floating point numbers (expressed in binary $x = \sum_i a_i(x) 2^{-i}$ or decimal $x = \sum_i b_i(x) 10^{-i}$) is a stable (i.e., continuous) operation on the abstract real numbers ${\bf R}$, but when viewed in a binary or decimal representation system becomes "uncomputable". For instance, the square of $1.414213\dots$ could be either $1.99999\dots$ or $2.0000\dots$, depending on exactly what is going on in the $\dots$; hence questions such as "what is the first digit of the square of $1.414213\dots$" are uncomputable. But this is an artefact of the numeral representation system used and is not an indicator of any lack of stability or computability for any actual computational problem that involves the abstract real numbers (rather than an artificial problem that is sensitive to the choice of numeral representation used). In contrast, floating point division when the denominator is near zero is a true singularity; regardless of what numeral system one uses, this operation is genuinely discontinuous (in a dramatic fashion) on the abstract reals and generates actual instabilities that cannot be explained away as mere coordinate singularity artefacts.

Returning back to matrices, whereas the individual eigenvectors $u_i(A)$ of a real symmetric matrix $A$ are not uniquely defined (there is a choice of sign for $u_i(A)$, even when there are no repeated eigenvalues) or continuously dependent on $A$, the spectral measure $\mu_A := \sum_i \delta_{\lambda_i(A)} u_i(A) u_i(A)^T$ is unambiguous; it is the unique projection-valued measure for which one has the functional calculus $$ f(A) = \int_{\bf R} f(E)\ d\mu_A(E)$$ for any polynomial $f$ (or indeed for any continuous function $f \colon {\bf R} \to {\bf R}$). The spectral measure $\mu_A$ depends continuously on $A$ in the vague topology; indeed one has the inequality $$ \| f(A) - f(B) \|_F \leq \|f\|_\text{Lip} \|A-B\|_F$$ for any real symmetric $A,B$ and any Lipschitz $f$, where $\|\|_F$ denotes the Frobenius norm (also known as the Hilbert-Schmidt norm or 2-Schatten norm). This allows for the possibility for stable computation of this measure, and indeed standard algorithms such as tridiagonalisation methods using (for instance) the QR factorisation and Householder reflections do allow one to compute this measure in a numerically stable fashion (e.g., small roundoff errors only lead to small variations in any test $\int_{\bf R} f(E)\ d\mu_A(E)$ of the spectral measure $\mu_A$ against a given test function $f$), although actually demonstrating this stability rigorously for a given numerical SVD algorithm does require a non-trivial amount of effort.

The practical upshot of this is that if one uses a numerically stable SVD algorithm to compute a quantity that can be expressed as a numerically stable function of the spectral measure (e.g., the inverse $A^{-1}$, assuming that the spectrum is bounded away from zero), then the computation will be stable, despite the fact that the representation of this spectral measure in eigenvalue/eigenvector form may contain coordinate instabilities. In examples involving eigenvalue collision such as the one you provided in your post, the eigenvectors can change dramatically (while the eigenvalues remains stable), but when the time comes to apply the SVD to compute a stable quantity such as the inverse $A^{-1}$, these dramatic changes "miraculously" cancel each other out and the algorithm becomes numerically stable again. (This is analogous to how a stable floating point arithmetic computation (avoiding division by very small denominators) applied to an input $x = 1.99999\dots$ and an input $x' = 2.00000\dots$ will lead to outcomes that are very close to each other (as abstract real numbers), even though all the digits in the representations of $x$ and $x'$ are completely different; the changes in digits "cancel each other out" at the end of the day.)

[The situation is a bit more interesting when applying the SVD to a non-symmetric matrix $A = \sum_i \sigma_i(A) u_i(A) v_i(A)^T$. Now one gets two spectral measures, $\mu_{(A^* A)^{1/2}} = \sum_i \delta_{\sigma_i(A)} v_i(A) v_i(A)^T$ and $\mu_{(AA^*)^{1/2}} = \sum_i \delta_{\sigma_i(A)} u_i(A) u_i(A)^T$ which are numerically stable, but these don't capture the full strength of the SVD (for instance, they are not sufficient for computing $A^{-1}$). The non-projection-valued spectral measure $\mu_A = \sum_i \delta_{\sigma_i(A)} u_i(A) v_i(A)^T$ does capture the full SVD in this case, but is only stable using the vague topology on the open half-line $(0,+\infty)$, that is to say $\int_0^\infty f(E)\ d\mu_A(E)$ varies continuously with $A$ as long as $f$ is a test function compactly supported in $(0,+\infty)$, but is unstable if tested by functions that do not vanish at the origin. This is ultimately due to a genuine singularity in the polar decomposition of a non-selfadjoint matrix when the matrix becomes singular, which in one dimension is simply the familiar singularity in the polar decomposition of a complex number near the origin.]

| cite | improve this answer | |
  • 3
    $\begingroup$ Are you sure about the Lipschitz estimate in the operator norm? On $B(\ell^2)$ already the absolute value map is not Lipschitz continuous on the self-adjoint part, so I would expect the constant in this inequality to blow up as the size of the matrices increases. It is fine for the Hilbert-Schmidt norm, though. $\endgroup$ – Mateusz Wasilewski Aug 24 at 19:41
  • 1
    $\begingroup$ Oops, you are right, thanks! $\endgroup$ – Terry Tao Aug 24 at 20:03
  • 1
    $\begingroup$ So the purpose of the spectral decomposition is to extend a continuous function $f:\mathbb R \to \mathbb R$ to matrices, by using the identity $f(U\Sigma U^*) = U f(\Sigma) U^*$. Have I understood that correctly? $\endgroup$ – ogogmad Aug 25 at 7:46
  • 12
    $\begingroup$ I would just like to remark to this very nice answer that there is no a priori reason why we could not directly implement the spectral measure on a computer, nor is there any reason, other than historic accident, that the floating-point format uses numerically unstable representation of the reals. Signed binary (with digits $-1, 0, 1$) is numerically stable. $\endgroup$ – Andrej Bauer Aug 25 at 8:12
  • 1
    $\begingroup$ Perhaps it is helpful to think about some of the standard things you would want to do with an SVD decomposition from this perspective. For example, you often replace a matrix $\sum \lambda_i u_i v_i^T$ by a low rank approximation $\sum_{\lambda_i > \kappa} \lambda_i u_i v_i^T$ for some cut off $\kappa$. This is integrating the spectral measure against $x \chi_{\mu}(x)$ where $\chi_{\mu}$ is the characteristic function of $(\kappa, \infty)$. (continued) $\endgroup$ – David E Speyer Aug 25 at 14:26

The SVD decomposition falls under the family of phenomena where discontinuity implies non-computability. (Intuitively, this is because at the point of discontinuity infinite precissions is required.)

In this particular case we speak of the (dis)continuity of a multivalued function which takes a matrix to any of its decompositions, or better, the non-existence of a realizer for the statement "For every matrix $M$ there exist suitable $U$, $\Sigma$, $V$ yielding an SVD decomposition of $M$." I belive this statement has no conitnuous realizer in function realizability, and hence no computable one either.

Some other examples of this phenomenon are:

  • The sign function $\mathrm{sgn} : \mathbb{R} \to \{-1,0,1\}$ is discontinuous, therefore non-computable. In fact, every computable map $\mathbb{R} \to \{0,1\}$ is constant.
  • The rank of a matrix is non-computable.
  • Gaussian elimination (as taught in school) cannot be performed because testing for zero is non-computable.
  • The number of distinct zeroes of a polynomial is non-computable.

So why are these, along with your observation, a problem? There are several answers, depending on the context.

In floating-point numerics, all computations are done with a fixed finite precision and numerical errors are simply unavoidable. In this setting the non-computability manifests itself as numerical instability. In your case, we might simply compute the wrong decomposition.

In some situations we can restrict to computation in a subring of $\mathbb{R}$ in which the problem disappears. For example, many of the above problems are non-existent when we restrict to $\mathbb{Q}$ or the algebraic numbers.

In exact real-arithmetic there are no numerical errors, as precision always adapts automatically to achieve the desired result. In this setting non-computability really is non-computability. The algorithm will diverge at points of discontinuity. In your case, it will just run forever trying to determine in which of the two cases it is.

There are models of real-number computation that pretend we can perform exact zero-testing, notably the Blum-Shub-Smale model. They are often used in computational geometry to side-step questions about non-computability. There are various theorems guaranteeing that a small perturbation of the input can get us out of trouble, at the price of possibly computing the wrong result.

| cite | improve this answer | |
  • 4
    $\begingroup$ I thought your question was, and I quote: "So my question is, how do numerical analysts get around this problem? Or why is this apparently not a problem?" What about SVD? Are you asking whether it is computable? $\endgroup$ – Andrej Bauer Aug 24 at 12:40
  • 1
    $\begingroup$ Ok, at this point I do not know what I am supposed to answer, sorry. I mean, I can answer "yeah, it seems so", but that's not what you're expecting. $\endgroup$ – Andrej Bauer Aug 24 at 12:41
  • 4
    $\begingroup$ I really liked this answer. $\endgroup$ – Nik Weaver Aug 24 at 14:21
  • 5
    $\begingroup$ @alephzero: stability and comptuabilithy don't make sense for a fixed matrix. There needs to be some sort of a parameter somewhere, because for a fixed matrix there is nothing to compute, as the SVD decomposition is fixed as well. Perhaps I misunderstand your point. $\endgroup$ – Andrej Bauer Aug 24 at 18:54
  • 2
    $\begingroup$ @AndrejBauer The confusion here is backwards vs. forwards stability. When people say that the SVD is stable, they mean that it is backwards stable, whereas here you are talking about forwards stability. $\endgroup$ – Nick Alger Aug 25 at 3:01

This is primarily an issue of backwards vs. forwards stability. Good SVD algorithms are backwards stable in the sense that the computed singular values and singular vectors are the true singular values and singular vectors of a slightly perturbed problem. You may see this by noting that while $P$ may change drastically as you change $\epsilon$, the product $PDP^T$ changes negligibly.

The SVD is not forwards stable when the singular values have small spectral gap, as your example demonstrates and other answers here discuss in more detail.

For more on backwards and forwards stability, see, e.g., this post and the links therein: https://math.stackexchange.com/a/78907/3060

SCIPY uses LAPACK; some details on the stability of the algorithm are provided here: https://www.netlib.org/lapack/lug/node97.html

| cite | improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.