The SVD subroutine computes the singular value decomposition for a numerical matrix. The input to the SVD subroutine is as follows:

a

is the input matrix that is factored as described in the following discussion.

The SVD subroutine returns the following output arguments:

u

is an orthonormal matrix

q

is an vector that contains the singular values

v

is an orthonormal matrix

If , the SVD subroutine factors a real matrix into the form

where

and contains the singular values of . The columns of contains of the orthonormal eigenvectors of , and contains the orthonormal eigenvectors of . contains the square roots of the eigenvalues of and , except for some zeros.

If , a corresponding decomposition is done where and switch roles:

where

The singular values are sorted in descending order.

For information about the method used in the SVD subroutine, see Wilkinson and Reinsch (1971).

The following example is taken from Wilkinson and Reinsch (1971):

a = {22  10   2   3   7,
     14   7  10   0   8,
     -1  13  -1 -11   3,
     -3  -2  13  -2   4,
      9   8   1  -2   4,
      9   1  -7   5  -1,
      2  -6   6   5   1,
      4   5   0  -2   2};
call svd(u, q, v, a);
print u, q, v;

/* check correctness of factors */
zero = ssq(a - u*diag(q)*v`);
reset fuzz;           /* print small numbers as zero */
print zero;

The matrix is rank-3 with exact singular values , , , , and . Because of the repeated singular values, the last two columns of the matrix are not uniquely determined. A valid result is shown in Figure 23.223:

The SVD routine performs most of its computations in the memory allocated for returning the singular value decomposition.