function FactorizationReview()

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% QR factorization
 
A = randn(8,5); b = randn(8,1); c = randn(5,1);
[Q,R] = qr(A,0);
% xLS =         % minimizes ||Ax - b|| assuming rank(A) = 5;
xLS = R\(Q'*b);

% v =           % a vector that is orthogonal to A(:,j),j=1:4
v = Q(:,5);     % By comparing columns in A = QR we see that the 
                % range of A(:,1:4) is the range of Q(:,1:4)
[Q,R] = qr(A);
% Z =           % the columns of this matrix span the null space of A'
Z = Q(:,6:8);   % A' = R'Q' so A'Q = R'. Look at the last three columns

% w =           % Solves A'w = c
y = R(1:5,1:5)'\c; w = Q(:,1:5)*y;    % R'(Q'w) = c

[Q,R,P] = qr(A);
% xB = approximately minimizes ||Ax - b|| assuming that A is close to having rank
% 3 and xB has three nonzero components.
btilde = Q(:,1:3)'*b; z = R(1:3,1:3)\btilde; xB = P*[z;zeros(2,1)];
      

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SVD
A = randn(8,4); A = [A A*randn(4,1)]; b = randn(8,1);  c = randn(5,1);
[U,S,V] = svd(A,0);
% xLS =         % minimizes ||Ax - b|| assuming rank(A) = 4;
xLS = V(:,1:4)*((U(:,1:4)'*b)./diag(S(1:4,1:4)));

[U,S,V] = svd(A);
% Z =           % the columns of this matrix span the null space of A'
Z = U(:,5:8);   %    A = USV', A' = VS'U' A'U = VS'  look at the last three columns

% w =           % Solves A'w = c
                % V(S'(U'w) = c
f = V'*c; g = [f./(diag(S(1:5,1:5)));zeros(3,1)]; w = U*g
       
% Q =           % A matrix of eigenvectors for A'*A
Q = V;          % A'A = (USV')'(USV') = V(S'S)V', i.e., V'*(A'A)V = S'S = diagonal

% Z =           % A matrix of eigenvectors for A*A'
Z = U;          % AA' = (USV')(USV')' = U(SS')U' i.e., U'(AA')U = SS' = diagonal

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Schur decomposition
A = randn(5,5); A = A'*A;  b = randn(5,1);
[Q,D] = eig(A);

%  f =        % The vector of eigenvalues for Q(:,1:2)'*A*Q(:,1:2)
f = diag(D(1:2,1:2))  % D(1:2,1:2) = Q(:,1:2)'A*Q(:,1:2)

%  x =        % the dominant eigenvector of inv(A)
[m,k] = min(abs(diag(D))); x = Q(:,k) ; %Ax = mu*x  inv(A)x = (1/mu)x
%  x =        % the dominant eigenvector of inv(A-7I)
[m,k] = min(abs(diag(D)-7)); x = Q(:,k) ; %Ax = mu*x  inv(A-7I)x = (1/(mu-7))x

%  x =        % Solves (A+ 5I)(A+2I)x = b
d = diag(D); w = Q'*b, y = w./(d+5); z = y./(d+2); x = Q*z;% Q(D+5I)(D+2I)Q'x = b

%  d =        % the vector of eigenvalues for A^10
d = diag(D).^10;   % Ax = mu*x  (A^10)x = (mu^10)x


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\
% Tridiagonalization
A = randn(500,500); A = A'*A; b = randn(500,1);
[U,T] = hess(A);

% x =    % Solves (A+5I)x = b
y = U'*b; z = (T+5*eye(500,500))\y; x = U*z; % (A+5I)x = U(T+5I)U'x = b
% Note that T+5I is pos def and tridiagonal so could use band cholesky.

% v =    % The vector A*U(:,500) 
alfa = T(500,500); beta = T(499,1); 
v = alfa*U(:,500)+beta*U(:,499);% AU  = UT, so A*U(:,500) = U*T(:,500)



