What should be the criteria for accepting/rejecting singular values?

Question

I am solving a system using singular value decomposition. The singular values (before scaling) are:

1.82277e+29
1.95011e+27
1.15033e+23
1.45291e+21
4.79336e+17
7.48116e+15
8.31087e+12
1.71838e+11
5.63232e+08
2.17863e+08
9.02783e+07
1.72345e+07
1.73889e+05
8.09382e+02
2.16644e+00

I have found that accepting all the singular values and their associated contribution to my solution vector yields poor results. I scale them all by the largest number, yielding singular values of:

1.0
1.06986e-02
6.31091e-07
7.97089e-09
2.62971e-12
4.10428e-14
4.55948e-17
9.42732e-19
3.08998e-21
1.19523e-21
4.95281e-22
9.45510e-23
9.53980e-25
4.44040e-27
1.18854e-29

The best solution only starts to become bad if I include the last two, and only become good around the $10^{-19}$ term.

There is a sharp drop in accuracy when I include the last 2 terms. Why is that? What are the criteria for including/not including singular values?

My matrix equation comes from a linear least squares fitting where I am using a polynomial basis set to fit some noisy data I created. I am solving the standard overdetermined system ($m \times n$ matrix where $m \gg n$) by multiplying each side ($A\cdot X = B$) by the transpose of $A$ ($A^\top A X = A^\top X$) and performing SVD on that.

I am judging the answers to my solutions by how well it approximates my noisy data.

I have also noticed that, even on the 'good' fits, I am not fitting very well near zero (my data ranges from $-10$ to $10$). Why is that?

Answer 1

The standard tolerance for forming a pseudoinverse is to only invert singular values that are at least $\max(m,n) \epsilon \|A\|_2$, where $A$ is $m \times n$, $\epsilon$ is the machine precision, and $\|A\|_2$ coincides with the largest singular value of $A$.

With that said, as J.M. mentioned, it is much more stable to avoid forming $A^H A$:

First, we need to compute the singular value decomposition $$ [U,\Sigma,V] = \mathrm{svd}(A) $$ then, we can define the pseudoinverse through $$ A^\dagger = V\; f(\Sigma)\; U^H, $$ where $$ f(\sigma) = \begin{cases}1/\sigma,\quad\sigma \ge \max(m,n) \epsilon \|A\|_2,\\ 0,\quad\text{otherwise}\end{cases} $$

The solution can then be computed as $$ X = A^\dagger B. $$

Answer 2

Augh!! No, no, a thousand times, no!

The reason people use SVD is precisely to avoid having to form the cross-product matrix $\mathbf A^\top\mathbf A$, since the formation of this matrix is a nice recipe for forming ill-conditioned linear systems! The decomposition is meant to be applied directly to $\mathbf A$. (See also some of my previous answers.)

I have mentioned to you before that the usual criterion for zeroing out singular values is to compare them with the product of the largest singular value and machine epsilon. However, this is rendered moot by your formation of the cross-product matrix. Please do try running the decomposition again, but this time, on the design matrix itself instead of the cross-product matrix. Any other way is flagrant abuse of the decomposition.

Answer 3

I think that several people here have provided valuable tips for your problem.

For future reference, however, your question of how to solve an ill-posed linear least squares problem could be answered by looking at the immense body of litterature on this problem.

Specifically, you could use the TSVD (Truncated Singular Value Decomposition) as a simple method of obtaining the solution: $$ \mathbf{x}_k = \sum_{i=1}^{k} \frac{\mathbf{u}_i^H \mathbf{b}}{\sigma_i}\mathbf{v}_i $$ where $\sigma_i$ is the i'th singular value, $\mathbf{u}_i$ and $\mathbf{v}_i$ are the i'th column in the matrices $U$ and $V$ from the factorization $USV^H=A$, $\mathbf{b}$ is the right-hand side of your problem, and the notation $\mathbf{u}^H$ means the complex conjugate of the entries and then turned to a row vector, such that $\mathbf{u}^H \mathbf{b}$ yields a scalar (dotproduct). Your solution is thus the vector $\mathbf{x}_k$.

The main problem in this setting, aside from being forced to compute the SVD which is quite expensive, is how to choose the number of singular values to use, i.e. $k$. Again, there are a number of ways to do that, but the most popular ones would be the Discrepancy principle, the Generalized Cross-Validation method, and the L-curve.

All of this (and much more) is implemented, in Matlab, in the excellent toolbox Regularization Tools, written by Prof. Per Christian Hansen, who has also published several papers and a few books on inverse problems. The toolbox is easy to use and should be quite easy to translate to other programming languages.

In conclusion, while others have provided important insights on your application that suggest other approaches are more appropriate, the above is a quick summary on how you could solve the problem if you still need to.