% Math 32a Fall 2003 Richard Palais Brandeis Univ.
% Discussion of Project 1
% Here is how I would go about developing the solution
% to Project 1. I would not start with an M-File, but
% first try to develop commands to implement the various
% steps of the algorithm in a text file such as this,
% and test them out in the mMatlab Command Window, and
% then when the pieces of my code seemed to be working
% properly, I would put them in a function M-file and
% add error checking. (It might even pay to have an
% intermediate step consisting of a script M-File, but
% for a simple project such as Project 1, that is not
% necessary.
% Recall, the goal is to apply Gram-Schmidt algorithm to
% the rows of a matrix M.
% Here again is a statement of the algorithm:
% 1) First normalize the initial row.
% 2) Then, in order, project each row (the current_row)
% onto the (already orthonormalized) earlier rows.
% 3) Subtract this projection from the current row to
% get the projection normal to earlier rows, and then
% 4) normalize the result.
% OK, here goes:
% To normalize the initial row of M, namely M[1,:].
% Compute its norm
norm1 = sqrt(M(1,:) * M(1,:)');
% We could also have written this as
% norm1 = sqrt( dot(M(1,:) , M(1,:)));
% or as:
% norm1 = norm(M(1,:));
% using Matlabs built-in functions dot and norm.
% Ignore possibility that norm1 might be zero for now
% and divide the first row by its norm
M(1,:) = M(1,:)/norm1;
% This completes step 1) above.
% Now start the loop that will orthonormalize later rows.
% Recall that the number of rows of M is size(M,1)
% so we want to iterate over the list of rows from
% row number 2 to row number size(M,1). In matlab-speak this is:
for row_num = 2:size(M,1)
current_row = M(row_num,:);
% project current_row of M onto each earlier row and add the projections
% to get the projection of current_row on span of earlier rows
proj = 0;
for j = 1:row_num - 1
proj = proj + dot(current_row,M(j,:)) * M(j,:);
end
% Now subtract proj from the current row to get the projection of the
% current row orthogonal to earlier rows.
u = current_row - proj;
norm1 = sqrt(u*u');
M(row_num,:) = u/norm1;
end;
%
%
% Set up a test matrix.
tst = [1 2 3 4;2 3 4 1;3 4 1 2; 4 1 2 3]
% set M equal to tst
M = tst
% Then copy the above code and test it on M
%
%
% next we add the boilerplate necessary to make our early
% code into an M-File, and also add error checking.
%
% First, here is the Help Comment:
% The function GramSchmidt takes as its single input
% argument and m by n matrix M of real numbers and returns
% an error message if the rows of M are linearly dependent,
% but if they are independent then it returns a matrix
% with orthonormal rows having the same size as M and such
% that the first k rows of Q have the same span as th first
% k rows of M.
% Programmer: Richard S. Palais
function X = GramSchmidt(M)
norm1 = sqrt(M(1,:) * M(1,:)'); % norm of the first row;
if (norm1 <0.000000001)
Q = 'First row is zero. Cannot continue.';
disp(Q);
return
else
M(1,:) = M(1,:)/norm1;
end % if
% Start loop to orthonormalize later rows.
for row_num = 2:size(M,1)
current_row = M(row_num,:);
% project current_row of M onto each earlier row and add the projections
% to get the projection of current_row on span of earlier rows
proj = 0;
for j = 1:row_num - 1
proj = proj + dot(current_row,M(j,:)) * M(j,:);
end % for j = 1:row_num - 1
% Now subtract proj from the current row to get the projection of the
% current row orthogonal to earlier rows.
u = current_row - proj;
norm1 = sqrt(u*u');
if (norm1 <0.000000001)
Q = 'Rows are linearly dependent, cannot continue.';
disp(Q);
return
else
M(row_num,:) = u/norm1;
end %if (norm1 <0.000000001)
end; % for row_num = 2:size(M,1)
X = M;
% END OF GramSchmidt
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Here is a somewhat selection of a few more solutions developed by
% various class members.
% First, David Diamondstone's version
% Call as N = GramSchmidt(M).
% M should be a rectangular
% matrix with linearly independent rows,
% and then N will be a matrix of the same size with orthonormal rows
% and the linear spans of the first k rows of M and N are the same.
% Programmer: David Diamondstone
function N=GramSchmidt(M)
r=size(M,1); c=size(M,2);
for i=2:r,
for j=1:i-1,
M(i:r:r*c)=M(i:r:r*c)-(M(i:r:r*c)*M(j:r:r*c)')/(M(j:r:r*c)*M(j:r:r*c)')*M(j:r:r*c);
if M(i:r:r*c)*2==M(i:r:r*c), 'The input matrix was not linearly independant.',
Failure_Row=i,
end;
end;
end;
for i=1:r,
M(i:r:r*c)=M(i:r:r*c)/sqrt(M(i:r:r*c)*M(i:r:r*c)');
end;
N=M;
% Comments on David's M-File.
% Note that David first orthogonalizes
% all the rows, and only then normalizes them.
% This has a number of advantages! Can you see them?
% Pros: It is VERY short and succinct.
% Cons: It is a little hard to understand.
% The error handling is less than prfect, and
% in particular, it fails to handle
% the case that the firt row is zero.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Here is Lacra Bintu's M-File:
% GramSchmidt takes a matrix with linearly independent rows and returns another matrix
% whose rows are orthonormal to each other.
%If the rows of the input matrix are not linearly independent, the GramSchmidt will give
% an error message.
% Programmer: Lacramioara Bintu
function F=GramSchmidt(M)
n=size(M,1); m=size(M,2); %matrix dimensions
if M(1,:)*M(1,:)'==0
F='The rows of your matrix are not linearly independent! Orthonormalization is impossible!';
else
F(1,:)=M(1,:)/((M(1,:)*M(1,:)')^(1/2)); %normalize the first row
for i=2:n
vi=0;
for k=1:i-1
vi=vi+(M(i,:)*F(k,:)')*F(k,:); %project the k'th row onto the space of the first k-1 orthonormalized vectors
end
ui=M(i,:)-vi; %subtract the projection from the original i'th row vector
u=(ui*ui')^(1/2); %find the norm of the difference
if (u<10^(-10)*mean(mean(abs(M))))
% if the norm is zero give error message
F='The rows of your matrix are not linearly independent! Orthonormalization is impossible!';
break %stop going through the loop
else
F(i,:)=ui/u; %normalize the difference obtained above
end
end
end
% Comments on Lacra's M-File
% Pros: Quite short, but well-commented and so reasonably easy to understand.
% Appears to work correctly for matrices with linearly independent rows,
% and error handling is good.
% The scaling in the line if (u<10^(-10)*mean(mean(abs(M))))
% is a very nice touch, since it correctly makes the test depend
% on the relative round_off error rather than absolute round-off error.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Chandni Rajan Valiathan's M-File
%This function takes as input, a rectangular matrix M of size m x n whose
%rows are linearly independant. It transforms the matrix to one whose rows
%are orthonormal and whose first K rows span the same subspace as that
%spanned by the first K rows of M. If it gets, as input, a matrix that has
%linearly dependant rows, it gives an error and no matrix. You can call the
%function as follows >> GramSchmidt (A) (where A is the matrix whose rows
%you want to be orthonormalized.)
%Programmer: Chandni Rajan
%Valiathan
%Math 32a Fall 2003
%Project 1
function F= GramSchmidt(M)
m= size (M, 1); %no. of rows
v= M(1,:);
normv= sqrt(v*v');
if (normv == 0) %if the first row is the zero vector, then the matrix is
%linearly dependant so display an error and return to the prompt
F= 'Input matrix has linearly dependant rows, so impossible to orthonormalize';
return;
else
F(1,:)= (1/normv)*v; %otherwise normalize it and replace the first row of M with
%the new orthonormalized vector.
end;
for i=2:m,
w=0; %reset w which stores the orthogonal projection
vi= M(i, :); %vi= ith row of vectors
for j=1:(i-1),
w= w + (vi*(F(j, :))')*F(j, :);
end; %orthogonal projection of v onto the subspace W of V
u= vi - w;
normu= sqrt(u*u'); %n=norm of u
if (normu < 1e-5)
F= 'Input matrix has linearly dependant rows, so impossible to orthonormalize';
return;
%vectors not linearly dependant so error displayed
else
F(i,:)= (1/normu)* u;
end;
end;
% Comments: Clear and well-documented. Works correctly and has good error handling.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Anny Jones M-File is very carefully commentes and
% works correctly with both "good" and "bad" matrices
% GramSchmidt
% This program takes as input a rectangular matrix (M) of arbitrary size m x n,
% and assuming that m rows are linearly independent, transforms (M) into another
% m x n matrix (Q) whose rows are orthonormal, vis-a-vis the Gram-Schmidt
% Algorithm. The subspace spanned by the first k rows of (Q) is equal
% to the subspace spanned by the first k rows of (M).
% Programmer: M. Anny Jones
% M = input('matrix M = ');
function GramSchmidt(M)
% First I will determine the size of input matrix (M):
m = size(M,1); n = size(M,2);
% I will duplicate matrix (M) and overlay each row as it is calculated
Q = M;
% Next, I will define Q1, the first element of output matrix (Q):
w = M( 1, 1:n );
% Since we know that the norm of vector w is the square root of the inner product
% , which can be written w*w' when w is a row vector:
normw = (w*w')^0.5;
if normw == 0
disp('The first vector of (M) is zero. Please correct.')
return
end
Q( 1, 1:n ) = w./normw;
% Here I have overlaid my first vector (Q1, if you will) over (M)
% Once m rows have been replaced, matrix Q will be the solution.
% Now, using Q1, we can iterate the process to find the final elements Qi:
for i = 2:m,
w = M( i, 1:n );
s = 0;
for j = 1:i-1,
v = (w*Q( j, 1:n)')*Q( j, 1:n);
s = s + v;
end
% Here, s is the projection of Mi onto Wi-1, ie. the projection of the
% ith row of M onto the subspace spanned by the first (i-1) rows of M (or Q).
% I have used a second 'for' to represent the summation over j=1 to j=i-1
u = w - s;
if (u*u' < 0.000001)
disp('Program cannot be completed! Please check that the')
disp('rows of the entered matrix are linearly independent.')
return
else
normu = (u*u')^0.5;
Q( i, 1:n ) = u./normu;
end
end
% Now I have found each subsequent row of Q and the result is displayed.
Q
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%