🚀 go-pugleaf

>Newsgroups: sci.math.num-analysis,comp.lang.fortran
>Subject: SVD Troubles
>From: ajayshah@alhena.usc.edu (Ajay Shah)
>Date: 26 Jun 91 13:14:40 GMT
>
>
>I get different results in doing a SVD using Numerical Recipes
>and using Linpack.
>
>The diagonal elements of W using LINPACK:
>
>    8546.9933049439    499.76058704715    237.98494198145    176.87448423245
>    168.17932502747    164.22538937388    155.90921231184    154.22581082297
>    1.1581996038494D-06    7.5680111867615D-07
>
>Using Numerical Recipes:
>
>    8546.9933049439    499.76058704715    237.98494198145    176.87448423245
>    168.17932502747    154.22581082297    155.90921231184    164.22538937388
>    7.5680116056573D-07    1.1581996113596D-06


[LINPACK and Numerical Recipes versions of FORTRAN program deleted.]

[Matrix deleted.]


After carefully reading what the "Numerical Recipes" book has to say (in
section 9 of chapter 2) about Singular Value Decomposition (SVD), I have come
to the conclusion that there is no real discrepancy between the results you
obtained with LINPACK and the results you obtained with the Numerical Recipes
(NR) subroutine.

Suppose that a matrix A has a SVD given by U*W*V', where W is a diagonal
matrix with non-negative entries (the singular values of A), the matrix U is
column-orthogonal, and the square matrix V is orthogonal.  Even when all of
the singular values (SVs) of A are distinct (as they are in your example),
the SVD of A is *not* unique, since permutations of the SVs are permitted (as
long as the corresponding columns of the matricies U and V are permuted in
the same way).  LINPACK was nice enough to list the SVs in decreasing order,
whereas NR did not list them in order.  However, both routines gave you
essentially the same *set* of SVs; the only real question concerns the number
of significant figures to which the corresponding SVs agree.

All of the SVs which are greater than 1 agree to all 14 significant figures.
You can't ask for better agreement than that!  The SVs which are less than 1
(the smallest two SVs) agree to only 7 or 8 significant figures.  Since you
are using double precision, this appears to be cause for concern, at first.
It is not really cause for concern if we understand how particularly small
SVs should be treated in the SVD.  It may be desirable to simply set such SVs
equal to *zero*, in which case, agreement to only 7 or 8 out of 14
significant figures becomes a moot point.  Before we will know whether we
should do this, we must first answer the question, "What constitutes
'particularly small'?"

Compare the largest SV to the smallest SV; the ratio of largest to smallest
is approximately 10**10 (the "condition number" of the matrix A).  Since
double precision is giving you about 14 significant figures, I would estimate
the machine's floating-point precision (machine epsilon) to be about
10**(-13).  If the condition number of A were greater than the reciprocal of
the machine epsilon, 10**13, the matrix A would be considered "ill-
conditioned" (i.e., poorly suited to some numerical computing techniques,
especially to the solving of linear systems by certain algorithms).  You have
about 3 orders of magnitude to go before your matrix A would be considered
ill-conditioned.

It seems to me that the closer the matrix is to being ill-conditioned, the
more one ought to consider the smallest SV to be "particularly small".  Given
a machine epsilon of 10**(-13), a difference of 10 orders of magnitude
between the largest and smallest SVs is significant.  (If the machine epsilon
were 10**(-100), a difference of 10 orders of magnitude would not seem like
much.) In light of the condition number and the machine epsilon, the smallest
two SVs are small enough that you should probably set them equal to zero
before using the results of the SVD in further computations.  With this
adjustment, LINPACK and NR give *exactly* (to *all* reported significant
figures) the same set of SVs.

I do not know how you intend to use the computed SVD of your matrix A, but if
you plan to solve linear systems such as Ax=b, where the vector 'b' is given
and the vector 'x' is unknown, then setting the smallest two SVs to zero may
well give *more* accurate solutions to Ax=b than if you use all of the SVs as
given.  (See the "Numerical Recipes" book, section 9, chapter 2, regarding
this claim.) Here, we are talking about solving Ax=b in the "least squares
sense", which means to find a vector 'x' which minimizes the 2-norm
(Euclidean length) of r=Ax-b (the residual vector).

DISCLAIMER: I am not a full-time numerical analyst, so if a more knowledgable
and experienced reader would like to let me know (via e-mail) whether they
agree/disagree with this treatment of the question, I would appreciate the
feedback.  Corrections of any blatantly incorrect statements should probably
be posted directly to the news group.


~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Frederick W. Chapman (N3IJC)          Campus Phone:  (215) 758-3218
User Services Group                   Network Server UserID:  FC03
Computing Center                      Internet:  FC03@NS.CC.LEHIGH.EDU
Lehigh University                     Bitnet:  FC03@LEHIGH