%
disp(' ');
disp(' SVD.m  ver 1.0  January 8, 2010 ');
disp(' by Tom Irvine ');
disp(' ');
disp(' Singular Value Decomposition of a Matrix:  ');
disp(' A=USV^T ');
disp(' ');
disp(' V is the orthonomal input ');
disp(' U is the orthonomal output ');
disp(' S is the gain controls ');
disp(' ');
%
%  References:
%  http://web.mit.edu/be.400/www/SVD/Singular_Value_Decomposition.htm
%  http://www.cse.unr.edu/~bebis/MathMethods/SVD/lecture.pdf
%  http://www.mathworks.com/access/helpdesk/help/techdoc/ref/svd.html
%  http://www.ling.ohio-state.edu/~kbaker/pubs/Singular_Value_Decomposition_Tutorial.pdf
%  http://www.miislita.com/information-retrieval-tutorial/singular-value-de
%  composition-fast-track-tutorial.pdf
%  http://www.cs.iastate.edu/~cs577/handouts/svd.pdf
%  http://www.ee.ic.ac.uk/hp/staff/dmb/matrix/decomp.html
%  http://www.utdallas.edu/~herve/Abdi-SVD2007-pretty.pdf
%
%  Note that Matlab already has a built-in function:
%  [U,S,V] = svd(X)
%
%  The matrix S is unique but the matrices U and V are not,
%  but U and V must be a matched pair.
%
%
%
%
%  This script is meant to demonstate the underlying calculation.
%
clear A;
clear W1;
clear W2;
clear S;
clear U;
clear V;
clear yhat;
clear temp;
%
A = input(' Enter matrix name ');
%
sz=size(A);
nr=sz(1);
nc=sz(2);
%
iflag=0;
if(nr<nc)
    A=A';
    iflag=1;
    temp=nr;
    nr=nc;
    nc=temp;
end
%
W1=zeros(nr,nr);
W2=zeros(nc,nc);
S=zeros(nc,nc);
U=zeros(nr,nc);
V=zeros(nc,nc);
%
W1=A*A';
W2=A'*A;
%
% Calculate V
%
    [MS2,Eig2]=eig(W2)
%
    for(i=1:nc)
        j=nc+1-i;
        V(:,i)=MS2(:,j);
    end
%
    for(i=1:nc)
        k=nc+1-i;
        S(i,i)=sqrt(Eig2(k,k));
    end
%
% Gram-Schmidt orthonormalization
%
%%    clear x;
%%    clear yhat;
%%    x=V;
%%    n=nc;
%%    yhat=zeros(n,n);
%%    yhat(:,1)=x(:,1)/norm(x(:,1));
%
%%    for i=2:n
%%        for j=i:n
%%            x(:,j)=x(:,j)-(yhat(:,i-1)'*x(:,j))*yhat(:,i-1);
%%        end
%%        yhat(:,i)=x(:,i)/norm(x(:,i));
%%    end    
%%    yhat
%%    V=yhat;
%
    U=A/(S*V');  %  because needs to be matched pair in terms of polarities

if(iflag==1)
    clear temp;
    temp=U;
    U=V';
    V=temp';
    A=A';
    U=U';
    V=V';
end
%
A
W1
W2
U
S
V
U*S*V'
error=A-(U*S*V')