clear all
echo on
!rm temp.txt
diary temp.txt
% enter a 'nice' matrix with orthogonal columns
orthog1=[ 4 -1  -1
          2  2  -1
          2  2   1
          4 -1   1
  ]
orthog1'*orthog1    % to test orthogonality of the columns
R=[ 1 3 2
    0 2 3
    0 0 2]              % upper triangular matrix
M=orthog1*R             % take linear combinations to destroy orthog.
M'*M                    % test that it is NOT orthogonal
keyboard
alpha=eye(3)            % matrix to store the multipliers of the projs
y1=M(:,1)        % first vector
u1=y1            % call it u1
y2=M(:,2)        % second vector
alpha(1,2)=  (y2'*u1)/(u1'*u1)    % calculate 'length' of proj of y2 on y1=u1
  %   find orthogonal vector to y1=u1
u2= y2- alpha(1,2)* u1         % y2 = u2 + alpha(1,2)*u1
disp('confirm u2, u1 are orthogonal')
u2'*u1        % confirm orthogonal
y3=M(:,3)      % third vector
alpha(1,3)=  (y3'*u1)/(u1'*u1)
alpha(2,3)=  (y3'*u2)/(u2'*u2)
u3=y3-alpha(1,3)*u1-alpha(2,3)*u2    % find third vector orthog to both u1,u2
disp('confirm u3 orthogonal to both u1 and u2')
u3'*u1        % confirm orthogonal
u3'*u2        % confirm orthogonal
U=[u1 u2 u3]
U'*U
disp('confirm that  M = U*alpha')
disp('original M')
M
disp('U*alpha')
U*alpha
disp('M-U*alpha')
M-U*alpha
disp('confirm span of orthogonal vectors')
rank([M U])
diary off
echo off
%%%%%%%%%%%%%% original M is:
%    4    10     3
%    2    10     8
%    2    10    12
%    4    10     7

