% Instead of solving the normal equations, we can use singular value decomposition 

% i.e. we will seek a 'solution' x to Ex=y that satisfies the
% overdetermined system (M > N) in a least squares sense.

% The SVD factorization of E is E = U*S*V' and is obtained simply with:

disp(' ')
disp('Solve the least squares fit by Singular Value Decomposition')
disp(' ')

[U,S,V] = svd(E,0);

% You can verify the orthonormality property of U and V: U'U=I

% We need the inverse of the diagonal matrix S, but don't need to do inv(S)
% because we know we can obtain the inverse of a diagonal matrix simply by
% taking the reciprocal of the diagonal elements

Si = 1./S;

% But this leaves us with Inf on the off diagonals. Find these values and
% replace them with zeros.

Si(find(~isfinite(Si))) = 0;

% Now the solution is

xs = V*Si*U'*y

% Compare: the columns should be the same:

[x xs]

% Now obtain the same solution simply by the Matlab matrix left divide
% operation. For an overdetermined system (M > N), this obtains the
% solution that satifies the system equations with the least mean squared
% error 

xmld = E\y;

% Compare: the columns should be the same

[x xs xmld]