% svdgames.m
% some applications of the SVD and the pseudoinverse
clear

% The first application is simple

A=randn(4,2);
b=randn(4,1);
x=A\b;		% Least squares solution to the overdetermined system

[U,S,V]=svd(A);		% A=U*S*V', U,V unitary, S diagonal
Spinv=(S*diag(1./diag(S)./diag(S)))';	% The pseudoinverse to a diagonal matrix
Apinv=V*Spinv*U';	% Apinv is the LEFT inverse of A
norm(x-Apinv*b,'fro')	% we get the same answer
pause

% Now, let us see a stupid example

A=[1 1 1 ; 1 1 1];
b=[2 ; 4];
x=A\b		% Not liked by matlab, also wrong solution
[U,S,V]=svd(A);
Spinv=S'/S(1,1)/S(1,1);	% invert the non-zero singular values
Apinv=V*Spinv*U';
x=Apinv*b	% Produces the minimum norm, least
		% squares solution to the stupid example
		% namely x1=x2=x3=1
pause

% Now, let us see some noise reduction

N=70;
vec=1:N;
A=vec'*vec;		% A has rank 1
noise=5*randn(N,N);
data=A+noise;
[U,S,V]=svd(data);
Strunc=zeros(N,N);
Strunc(1,1)=S(1,1);	% As A has rank 1 we only keep the largest singular value
Aapprox=U*Strunc*V';
norm(A-data,'fro')
norm(A-Aapprox,'fro')	% An improvement
