[Phys-l] singular value decomposition

[John Denker <jsd@av8n.com>]

On 03/29/2012 12:05 PM, I wrote:
> Once upon a time, a teacher had an ensemble of 2x2 matrices.
> He assigned a team of students to calculate the eigenvalues 
> and produce a histogram thereof.
>
> One student, Bob, worked all day calculating the matrix elements.
> Another student, Carol, worked all night calculating the eigen-
> values and plotting the histogram.

Let me tell you a little more about the calculation Bob did.

He found that the matrix is the sum of two terms, where each
term is a direct product.

The first term (the "big term") is of the form bb · bb^† where 
bb is a column vector and the dagger (†) denotes transpose, so 
bb^† is a row vector.  In terms of components, 
          bb = [b, b]^† 
where      b = 1 +- 0.06%

The other term (the "small term") is of the form ss · ss^† where
in terms of components,
         ss = [s, -s]^†
where     s = 1/sqrt(810) +- 0.06%

Bob dutifully multiplied out the direct products and added the
two terms.  The rest of the story proceeds as described in my
previous note.

==========

Here's yet another part of the story.  There is another student,
Ted, whom I have not previously mentioned.

Ted noticed that it might be helpful to perform a change of
basis.  If we rotate the axes 45 degrees, you can see that bb 
lines up with one axis and ss lines up with the other, exactly.

If we rotate back to the original axes, we can write
       bb = R [b sqrt(2), 0         ]^†
       ss = R [0        , s sqrt(2) ]^†

where R is the 45-degree rotation matrix, and is known exactly,
with no uncertainty whatsoever.

Now multiplication distributes over addition ... and transposition
also distributes over addition, so Ted can write the entire matrix
of interest as

           [ 2 b^2     0     ]
    M = R  [                 ] R^†                          [4]
           [ 0         2 s^2 ]

      = R S R^†

which defines S in the obvious way.

The RHS is called the /singular value decomposition/ (SVD) of the 
matrix M.  Ted could carry out the multiplications called for on 
the RHS ... but he is not obliged to do so.  It would produce the 
same numerical values that Bob got for the matrix elements of M.

  Note that in this case the SVD is identical to the eigenvalue
  decomposition (EVD).

In case you are wondering why Ted might want to carry around three
matrices (R, S, and R^†) instead of just one matrix (M), take
another look at equation [4].
  a) In general, given the SVD you can pick off the eigenvalues
   by inspection.  They are sitting there on the diagonal of S.
   (Remember, in this example the original goal was to find the 
   eigenvalues and histogram them.)
  b) In general, you can also pick off the eigenvectors by 
   inspection.  The columns of R are the eigenvectors.
  c) In this example, the matrix elements of S are (by construction) 
   statistically independent of each other, whereas the matrix elements 
   of M are highly correlated.  Therefore the SVD makes it much, much 
   easier to understand what is going on.
  d) Matrix S is numerically well behaved and easy to use, in contrast 
   to M which gives you terrible problems with "small differences between
   large numbers" as soon as you try to do anything with it.

Et cetera.  There's a lot more I could say.  Fat books have been written
on the subject of SVD and related methods.  This note is just meant to
serve as a reminder that such methods exist.  The next time you are 
playing in this space you can read up on the details.

================

This makes contact to something you already know, as follows:

You presumably learned in connection with quantum mechanics that you
can always diagonalize the Hamiltonian.  Well,
 -- The idea works for all matrices, not just Hamiltonians.
 -- It works in general, not just in quantum mechanics.
 -- It even works for matrices that are singular (not full rank).
 -- It even works for non-square matrices.

The SVD always exists.  Up to trivial re-arrangements, it is unique,
unless there is some degeneracy.

If the matrix is Hermitian, the SVD is the same as the EVD.  This is
good to know, because there are more tricks for finding the EVD than
the general SVD.

In general, there are no closed-form expressions for the SVD or EVD,
but if you have a numerical matrix (as opposed to algebraic abstractions) 
there are efficient algorithms for grinding out the decomposition.  
For example, LAPACK knows how to do it.  GSL knows how to do it.

Sometimes in special cases, as in the example we have been considering, 
you can find a good decomposition analytically (rather than numerically).

====

Also note that there is something I like to call a ballpark decomposition
(BPD).  That is, sometimes you can write M = R S R^† where S is not
necessarily diagonal ... just a lot better-behaved than M.  Especially
if there are more than two dimensions, it might be important to 
diagonalize part of M even if it is not worthwhile to diagonalize the
rest.  Also if M has a problem with "small differences between large
numbers" you might be able to make the problem go away without exactly
diagonalizing the system.  You can treat S as a diagonal matrix plus
a small perturbation.

=============================

Bottom line:  Just because you can calculate the matrix elements of M
doesn't mean you want to.  The most compact representation of M is not
necessarily the most informative.  Furthermore, even if you do calculate 
the matrix elements of M, you might want to keep the SVD representation 
around /also/.