function ShowLSSolvers()
% Solving the LS problem three ways
m = 5;
n = 3;
A = randn(m,n);
b = randn(m,1);
clc
disp('   Normal Equations      via A = QR         via A = USV''')
disp('-------------------------------------------------------------')
[xLSNE,rhoLS]  = LS_normaleqn(A,b);
[xLSQR,rhoLS]  = LS_qr(A,b);
[xLSSVD,rhoLS] = LS_svd(A,b);
format long
disp([xLSNE xLSQR xLSSVD])





  function [xLS,rhoLS] = LS_normaleqn(A,b)
% A is mxn with full column rank and b is mx1.
% xLS (nx1) minimizes || Ax - b || and rhoLS = || A*xLS - b ||.
% Uses the method of normal equations.
% min ||Ax - b|| is equivalent to solving A'Ax = A'b
C = A'*A;
d = A'*b;
G = chol(C,'lower');
xLS = G'\(G\d);
rhoLS = norm(A*xLS-b);

  function [xLS,rhoLS] = LS_qr(A,b)
% A is mxn with full column rank and b is mx1.
% xLS (nx1) minimizes || Ax - b || and rhoLS = || A*xLS - b ||.
% Uses the QR factorization.
% || Ax - b ||^2 = ||Q'*A - Q'b ||^2 = || Rx - d ||^2 
%                = || R1*x - d(1:n)||^2 + || d(n+1:m) ||^2
[Q1,R1] = qr(A,0);
xLS = R1\(Q1'*b);
rhoLS = norm(A*xLS-b);

  function [xLS,rhoLS] = LS_svd(A,b)
% A is mxn with full column rank and b is mx1.
% xLS (nx1) minimizes || Ax - b || and rhoLS = || A*xLS - b ||.
% Uses the singular value decomposition.
% || Ax - b ||^2 = || U'Ax - U'b ||^2 = || SVx - U'b ||^2 = || Sy - d ||^2
%                = || S(1:n,1:n)y - d(1:n) ||^2 + || d(n+1:m) ||^2.
[U,S,V] = svd(A,0);
xLS = V*((U'*b)./diag(S));
rhoLS = norm(A*xLS-b);