% play with stability of QR. lec 5 Math 116. Barnett 4/8/14

% Gram-Schmidt (needs clgs.m and mgs.m). See NLA Ch. 9, etc

m = 200; % matrix size
[U,~] = qr(randn(m,m)); [V,~] = qr(randn(m,m)); % make random orthog U,V
A = U * diag(1.3.^-(1:m)) * V';   % build matrix with decaying sing vals

[Q,R] = clgs(A);
figure; semilogy(diag(R),'+');      % r_ii only decays to 1e-8
figure; imagesc(log10(abs(Q'*Q))); colorbar   % not orthog! so not a QR decomp
norm(A - Q*R) / norm(A)  % is backward stable, not relevant since Q not orthog

[Q,R] = mgs(A);
figure; semilogy(diag(R),'+');       % r_ii goes down to 1e-16 as should
figure; imagesc(log10(abs(Q'*Q))); colorbar   % still not orthog, but better
norm(A - Q*R) / norm(A)  % is backward stable, not relevant since Q not orthog

% The professional way to do it:

[Q,R] = qr(A);                    % via Householder reflections
figure; semilogy(diag(R),'+');            % r_ii goes down to 1e-16 as should
figure; imagesc(log10(abs(Q'*Q))); colorbar   % orthog to e_mach
norm(A - Q*R) / norm(A)  % is backward stable
% good
