svd-decomposition.html

This is an example of a singular value decomposition. A = U E V^t

> with( LinearAlgebra ):
A := Matrix( [[ -2, 1, 2], [6, 6, 3]] );

Compute the singular values, which are the square roots of the spectrum( A*A ). Note that A*A is hermitian,
so if A is real, all of the eigenvalues will be positive reals as well.

> ( singVals, eVecs ) := Eigenvectors( HermitianTranspose( A ) . A ):
Map( x->sqrt(x), singVals ):

( cDim, rDim ) := Dimension( eVecs ):
for i from 1 to Dimension( singVals ) do
   for j from i to Dimension( singVals ) do
       if ( singVals[ i ] < singVals[ j ] ) then
           tempVal := singVals[ i ];
           singVals[ i ] := singVals[ j ];
           singVals[ j ] := tempVal;
           tempVec := Column( eVecs, i );
           eVecs[ 1..rDim, i ] := Column( eVecs, j );
           eVecs[ 1..rDim, j ] := tempVec;
       end if;
   end do;
end do;

( singVals, eVecs );

We now normalize all of the eigenvectors. Since A*A is hermitian, we know that all the eigenvectors are
already orthogonal. These normalized eigenvectors give us the matrix V.

> V := Matrix( eVecs ):
for i from 1 to cDim do
V[ 1..rDim, i ] := ( 1 / sqrt( Transpose( Column( V, i ) ) . Column( V, i ) ) ) . Column( V, i );
end do:
V;

The E matrix is built from the nonzero singular values on the diagonal. The overall shape is the same as A.

> E := Matrix( Dimension( A ) ):
for i from 1 to Dimension( singVals ) do
   if ( singVals[ i ] <> 0 ) then
       E[ i, i ] := singVals[ i ];
   end if:
end do:
E;

Now I'm trying to build U

> U := ( A . Column( V, 1 ) ) / singVals[ 1 ]:
for i from 2 to Dimension( singVals ) do
   if ( singVals[ i ] <> 0 ) then
       U := < U | (( A . Column( V, i ) ) / singVals[ i ] ) >;
   end if:
end do;
U;

That's it. A = U E V^t

> Equal( A, U . E . Transpose( V ) );

All this work could easily have been achieved using Maple's built-in functionality.
Note that while the built in function can spit out the singular values, it doesn't build a pretty

Sigma (E) matrix like I did above.

> ( mU, mVt ) := SingularValues( A, output=['U', 'Vt'] );
Equal( A, mU . E . mVt );