8
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ If you mention which tags you would like created either in the post or a comment, someone with higher rep can take care of it for you. $\endgroup$ – David Z Dec 12 '11 at 4:03
  • $\begingroup$ @DavidKetcheson: I'm not quite sure if we really need these specific tags... yet. Maybe when we have a lot of SVD questions... but I think linear-algebra ought to be a sufficient tag for now. $\endgroup$ – J. M. Dec 12 '11 at 8:05
11
$\begingroup$

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. $$

$\endgroup$
  • 2
    $\begingroup$ I think it's better to be explicit and say that the matrix 2-norm is precisely the same as the largest singular value... also, $\mathbf A$ doesn't need to be square to apply the threshold criterion (but one replaces the $n$ with $\max(m,n)$ for an $m\times n$ matrix). $\endgroup$ – J. M. Dec 12 '11 at 7:57
  • 1
    $\begingroup$ instead of A†=Uf(Σ)*VH shouldn't it be A†=Vf(Σ)*U'? $\endgroup$ – drjrm3 Dec 12 '11 at 16:35
  • $\begingroup$ Yes, thank you both. I just wrote a parallel pseudoinverse routine for Hermitian matrices a few days ago, and apparently didn't give enough thought to the differences relative to the general case. $\endgroup$ – Jack Poulson Dec 12 '11 at 16:44
  • 2
    $\begingroup$ Pseudoinverses are not a cure-all for ill-conditioned problems, though they do tend to help. High-order polynomial interpolation is notoriously ill-conditioned and is best avoided. I recommend that you read the Wiki article on Runge's phenomenon and then switch to a different basis for your interpolation (e.g., piecewise polynomials or splines). $\endgroup$ – Jack Poulson Dec 12 '11 at 17:30
  • 1
    $\begingroup$ This discussion is no longer related to your original question; whether or not your implementation of a questionable least squares approximation converges is not relevant to your original question. I suggest you post another question if you are still confused. $\endgroup$ – Jack Poulson Dec 13 '11 at 0:17
11
$\begingroup$

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.

$\endgroup$
  • 3
    $\begingroup$ Can't upvote this enough... this is the single most important point about using SVD techniques correctly. Remember that the singular values of $A^T A$ are the squares of the singular values of $A$! $\endgroup$ – khinsen Dec 12 '11 at 11:14
  • $\begingroup$ Precisely. Squaring a tiny number results in an even tinier number. Considering the fact that the 2-norm condition number is the ratio of the largest to the smallest singular value, this is one way of seeing why the normal equations approach can be a bad idea for large systems. $\endgroup$ – J. M. Dec 12 '11 at 11:39
2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ thanks, i'm aware of most of this. i have searched many resources - the problem is that my solutions become worse and worse as i grow my basis set being used. i have come to the conclusion that this is because the number of singular values being rejected as my basis set increases grows relative to the number of basis functions being used. that is, when i used n = 100 basis functions, i must throw out 95 of the singular values and this yields a worse result than if i use 10 basis sets and only throw out one singular value. $\endgroup$ – drjrm3 Dec 14 '11 at 20:10
  • $\begingroup$ Right - well, as others have suggested, you should look at your basis functions, and perhaps find a basis that better approximates your data. Splines could be an option. So could Kriging, though that is a different approach entirely, which has yielded good results for me on several occasions. $\endgroup$ – OscarB Dec 14 '11 at 20:49

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.