% This file contains several different implementations of algorithms for
% computing a QR factorization via the Gram-Schmidt process. The purpose is
% merely to illustrate concepts covered in class. These codes are NOT
% efficient, are not particularly numerically stable, and are written to
% handle only the case where the matrix has full column rank. (In
% particular, they cannot handle matrices with more columns than rows.)
%
% P.G. Martinsson, University of Texas at Austin, September 2024

function main_GS

rng(0)
m = 5;
n = 3;
A = randn(m,n);

[Q,R] = QR_GS(A);
%[Q,R] = QR_GS_proj(A);
%[Q,R] = QR_GS_blocked(A);
%[Q,R] = QR_MGS(A);
%[Q,R] = QR_MGS_blocked(A);
%[Q,R] = QR_via_triangular_transforms(A);

fprintf(1,'|A - Q*R|             = %10.3e\n',max(max(abs(A - Q*R))))
fprintf(1,'|trans(Q)*Q - eye(n)| = %10.3e\n',max(max(abs(Q'*Q - eye(n)))))

keyboard
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function computes a QR factorization via classical Gram - Schmidt.
% Note that the Q matrix overwrites the A matrix.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Q,R] = QR_GS(A)

n = size(A,2);
R = zeros(n,n);

for j = 1:n  
  %%% Compute the j'th column of A and subtract away the components of
  %%% A(:,j) that can be spanned with the existing basis vectors.
  for i = 1:(j-1)
    R(i,j) = A(:,i)'*A(:,j);  % Note that A(:,i) = Q(:,i) here!
    A(:,j) = A(:,j)- A(:,i)*R(i,j);
  end
  %%% Compute R(j,j).
  R(j,j) = norm(A(:,j));
  %%% Compute Q(:,j) and overwrite it on top of A(:,j):
  A(:,j) = (1/R(j,j))*A(:,j); 
end

Q = A;

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function computes a QR factorization via classical Gram - Schmidt.
% Note that the Q matrix overwrites the A matrix.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Q,R] = QR_GS_proj(A)

fprintf(1,'Warning: This code only demonstrates the math. It is VERY slow.\n')

m = size(A,1);
n = size(A,2);
R = zeros(n,n);

for j = 1:n
  for i = 1:(j-1)
    R(i,j) = A(:,i)'*A(:,j);
    A(:,j) = (eye(m) - A(:,i)*A(:,i)')*A(:,j);   % Super inefficient!
  end
  R(j,j) = norm(A(:,j));
  A(:,j) = (1/R(j,j))*A(:,j);
end

Q = A;

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Classical Gram-Schmidt, in blocked form.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Q,R] = QR_GS_blocked(A)

n = size(A,2);
R = zeros(n,n);

for j = 1:n
  % Compute the j'th column of R.
  R(1:(j-1),j) = A(:,1:(j-1))'*A(:,j);
  % Project away from span of Q(:,1:(j-1)):
  A(:,j)       = A(:,j) - A(:,1:(j-1))*R(1:(j-1),j);
  % Compute R(j,j):
  R(j,j)       = norm(A(:,j));
  % Compute Q(:,j) and overwrite it on the j'th column of A:
  A(:,j)       = (1/R(j,j))*A(:,j);
end

Q = A;

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Modified Gram-Schmidt
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Q,R] = QR_MGS(A)

n = size(A,2);
R = zeros(n,n);

for j = 1:n
  R(j,j) = norm(A(:,j));
  A(:,j) = (1/R(j,j))*A(:,j); % We overwrite A(:,j) by Q(:,j).
  %%% Now project the remaining columns away from the span of Q(:,j)
  %%% In the process, compute the j'th row of R.
  for i = (j+1):n
    R(j,i) = A(:,j)'*A(:,i);         % Recall: A(:,j) holds Q(:,j)!
    A(:,i) = A(:,i) - A(:,j)*R(j,i); 
  end
end

Q = A;

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Modified Gram-Schmidt, written in blocked form.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Q,R] = QR_MGS_blocked(A)

n = size(A,2);
R = zeros(n,n);

for j = 1:n
  R(j,j)       = norm(A(:,j));
  A(:,j)       = (1/R(j,j))*A(:,j);
  R(j,(j+1):n) = A(:,j)'*A(:,(j+1):n);
  A(:,(j+1):n) = A(:,(j+1):n) - A(:,j)*R(j,(j+1):n);
end

Q = A;

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Drive A to orthonormal form via triangular transformations.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Q,R] = QR_via_triangular_transforms(A)

fprintf(1,'Warning: This code only demonstrates the math. It is VERY slow.\n')

n      = size(A,2);
R      = eye(n,n);
A_orig = A;

for j = 1:n
  %%% Build the triangular transformation matrix Rloc.
  Rloc        = eye(n);
  rjj         = norm(A(:,j));
  qj          = (1/rjj)*A(:,j);
  Rloc(j,j:n) = (1/rjj)*[1,-qj'*A(:,(j+1):n)];
  %%% Apply Rloc to A and to R to advance the factorization one step.
  A           = A*Rloc;
  R           = inv(Rloc)*R;
  %fprintf(1,'Step %3d: Error = %10.3e\n',j,max(max(abs(A_orig - A*R))))
end
  
Q = A;

return
