Implementation of Gram Schmidt in Matlab

 

http://faculty.salina.k-state.edu/tim/DAT/linAlg/appendixLA/QR.html

Implementation of Gram Schmidt

function U = gram_schmidt(V)
% GRAM_SCHMIDT - Classic Gram-Schmidt Process
%   Input vectors should be the columns of input matrix.
%   Output = unitary, orthogonal vectors in columns of matrix

n = size(V, 2); % number of columns

U(:,1) = V(:,1);

% find next orthogonal column
for ii = 2:n
    b = V(:,ii);
    c = b;
    for k = 1:ii-1
        a = U(:,k);
        c = c - (a*a'*b)/(a'*a);
    end
    U(:,ii) = c;
end

% now normalize
for ii = 1:n
    u = U(:,ii);
    U(:,ii) = u/norm(u);
end

Gram Schmidt and Orthogonal (QR) Factoring

The Gram Schmidt algorithm has two primary applications. In some physics applications one might have a set basis vectors which are not orthonormal. The Gram Schmidt algorithm can be used to find a new set of basis vectors that are orthonormal. The other application is a matrix factoring resulting in one matrix, , with orthonormal columns and another matrix, , which is upper triangular.

The discussion here focuses on QR factoring, but can be equally be applied to finding orthonormal basis vectors from the columns of a matrix.

Note

 

What is defined here is the classic Gram – Schmidt process. MATLAB has two functions that implement a modified version of the classic Gram – Schmidt process. The function orth( ) finds orthogonal basis vectors, while the function qr( ) factors a matrix into an orthogonal matrix and an upper triangular matrix. The results of both functions are different from what results from a direct implementation of the classic Gram – Schmidt process. According to the MATLAB documentation, the classic Gram – Schmidt process is numerically unstable. But for our introduction, we will stick with the classic algorithm.

6.11.5.1. The Gram Schmidt Algorithm

The columns of  are first formed from vector projections (see Projections Onto a Line), and then made unit length. Recall that when a vector is projected onto another vector, the vector representing the error between the projection and original vector are perpendicular to each other. Here it is the vector representing the error from projection is what we want.

Let matrix  be formed from column vectors.

We will first find a matrix, , who’s columns are orthogonal and are formed from projections of the columns of . For each column after the first, we subtract the projection of the column onto the previous columns. Thus we are subtracting away the projection leaving a vector that is orthogonal to the previous column vector.

The columns of  are then unit vectors of the columns of .

Then  is found from  and . Because  has orthonormal columns, . Because the matrix multiplication involves dot products of orthogonal vectors, matrix  will be upper triangular.

MATLAB has a command, [Q, R] = qr(A) that finds the factors  and .

a1 = [7 1 1]';
a2 = [6 10 3]';
a3 = [4 5 8]';
A = [a1 a2 a3];
b1 = a1;
b2 = a2 - (b1'*a2)/(b1'*b1)*b1;
b3 = a3 - (b1'*a3)/(b1'*b1)*b1 - (b2'*a3)/(b2'*b2)*b2;

q1 = b1/norm(b1);
q2 = b2/norm(b2);
q3 = b3/norm(b3);
Q = [q1 q2 q3];
R = Q'*A;
z = [0 0 0];
figure, axis([-2 8 -2 11 0 9]), daspect([1 1 1])
grid on
hold on
arrow3(z, a1', 'b')
arrow3(z, a2', 'b')
text(a2(1)+0.3, a2(2)-0.3, a2(3)-0.3, 'a2')
arrow3(z, a3', 'b')
text(a3(1)+0.3, a3(2)-0.3, a3(3)-0.3, 'a3')
arrow3(z, b1', 'r')
text(b1(1)+0.3, b1(2)-0.3, b1(3)-0.3, 'a1 == b1')
arrow3(z, b2', 'r')
text(b2(1)+0.8, b2(2)-0.3, b2(3)-0.3, 'b2')
arrow3(z, b3', 'r')
text(b3(1)+0.3, b3(2)-0.3, b3(3)-0.3, 'b3')
hold off
xlabel('x')
ylabel('y')
zlabel('z')
title('Gram-Schmidt Orthogonal Transform')

After Gram Schmidt, the columns of  and  are orthogonal vectors.