13.3. Finding Orthogonal Basis Vectors¶

We sometimes need to find a set of orthogonal basis vectors for the columns of a matrix. One such need is for vector projections (see Alternate Projection Equation).

One of three algorithms are typically used to find orthogonal basis vectors. A well known algorithm is the classic Gram–Schmidt process (CGS). It was proposed by Laplace (1749 - 1827) and later refined by Gram (1850 - 1916) and Schmidt (1876 - 1959). Gram–Schmidt is often taught in linear algebra courses, usually along with or shortly after the study of vector projections. However, the classic Gram–Schmidt process is known to be numerically unstable when the matrix is poorly conditioned with nearly dependent columns resulting in a loss of orthogonality [DEMMEL97]. Interestingly, the algorithm can be made more reliable by rearranged it into what is know as the modified Gram–Schmidt process (MGS) [GOLUB13]. MATLAB does not include functions for either the CGS or MGS algorithms. Of course CGS is not included because of its lack of numerical stability and the MGS is not essential because MATLAB’s qr function is faster and returns functionally the same results.

Some implementations of CGS and MGS return both an orthogonal matrix, $\bf{Q}$ , whose columns form basis vectors for the columns of $\bf{A}$ , and an upper triangular matrix, $\bf{R}$ , such that $\mathbf{A} = \mathbf{Q\,R}$ . The $\bf{Q}$ and $\bf{R}$ factors of $\bf{A}$ and the same as, or sometimes the negative of, the matrices returned by MATLAB’s qr function. The QR factorization described in QR Decomposition is the second algorithm for finding orthogonal basis vectors. Although QR returns the same functional results as MGS, it is implemented using a different algorithm and is considerably faster.

The third algorithm for finding orthogonal basis vectors is from the SVD. MATLAB provides a function called orth, which returns the $\bf{U}$ matrix from the economy singular value decomposition, [U, S] = svd(A,'econ'). For more details, see Projection and the Economy SVD. The SVD based algorithm yields a different set of vectors than either Gram–Schmidt or QR.

13.3.1. The Gram–Schmidt Algorithm¶

Although MATLAB does not provide a function that implements the Gram–Schmidt process, we will briefly describe it here as a reference for those that encounter it in other linear algebra studies.

The input to the Gram–Schmidt process is a matrix, $\bf{A}$ , and the output is a matrix, $\bf{Q}$ , whose columns are an orthogonal basis of $\bf{A}$ .

The columns of $\bf{Q}$ 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 the original vector are orthogonal to each other. Here we want to find the vector representing the error from projection.

Let matrix $\bf{A}$ be formed from column vectors.

$\mathbf{A} = \mat{ \vertbar{} \vertbar{} {} \vertbar{}; \bm{a}_1 \bm{a}_2 \cdots{} \bm{a}_n; \vertbar{} \vertbar{} {} \vertbar{}}$

We will first find a matrix, $\bf{B}$ , who’s columns are orthogonal and are formed from projections of the columns of $\bf{A}$ . 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 all of the previous column vectors.

$\begin{array}{rl} \bm{b}_1 &= \bm{a}_1 \\ \bm{b}_2 &= \bm{a}_2 - \frac{\bm{b}_1^T\, \bm{a}_2} {\bm{b}_1^T\,\bm{b}_1}\,\bm{b}_1\\ \bm{b}_3 &= \bm{a}_3 - \frac{\bm{b}_1^T\, \bm{a}_3} {\bm{b}_1^T\,\bm{b}_1}\,\bm{b}_1 - \frac{\bm{b}_2^T\, \bm{a}_3}{\bm{b}_2^T\,\bm{b}_2}\,\bm{b}_2 \\ &\vdots{} \\ \bm{b}_n &= \bm{a}_n - \frac{\bm{b}_1^T\,\bm{a}_3} {\bm{b}_1^T\,\bm{b}_1}\,\bm{b}_1 - \frac{\bm{b}_2^T\, \bm{a}_3}{\bm{b}_2^T\,\bm{b}_2}\,\bm{b}_2 - \cdots - \frac{\bm{b}_{n-1}^T\, \bm{a}_n} {\bm{b}_{n-1}^T\,\bm{b}_{n-1}}\,\bm{b}_{n-1} \\ \end{array}$

The columns of $\bf{Q}$ are the columns of $\bf{B}$ scaled to be unit vectors.

$\bm{q}_i = \frac{\bm{b}_i}{\norm{\bm{b}_i}}$

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];
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')

Fig. 13.3 After Gram–Schmidt, the columns of $\bf{B}$ and $\bf{Q}$ are orthogonal vectors.

13.3.2. Implementation of Classic Gram–Schmidt¶

Here is an implementation of the classic Gram–Schmidt (CGS) algorithm. Note that CGS is considered to be numerically unstable for some matrices, so it is for educational purposed only. Although, the CGS, MGS, and QR algorithms will return the same $\bf{Q}$ matrix for well conditioned input matrices.

function Q = gram_schmidt(A)
% GRAM_SCHMIDT - Classic Gram-Schmidt Process (CGS)
%   Input - A matrix. The algorithm operates on the columns.
%   Output - unitary matrix - columns are basis for A.

[m, n] = size(A);
Q = zeros(m, n);
Q(:,1) = A(:,1)/norm(A(:,1));    % the first column

% find next orthogonal column and normalize the column
for ii = 2:n
    b = A(:,ii);
    c = b;
    for k = 1:ii-1
        a = Q(:,k);
        c = c - (a*a'*b);
    end
    % now normalize
    Q(:,ii) = c/norm(c);
end

13.3.3. Implementation of Modified Gram–Schmidt¶

The modified Gram–Schmidt (MGS) algorithm improves, but does not completely correct, the stability problems of the CGS process. The core MGS algorithm in this function is from a code segment in Golub and Van Loan’s MATRIX Computations text [GOLUB13].

function [Q, R] = mod_gram_schmidt(A)
% MOD_GRAM_SCMIDT - Modified Gram-Schmidt Process
%   Variantion of the Gram-Schmidt process with
%   improved numerical stability.
%  This function is slower, but functionally
%  equivalent to MATLAB's qr function.
%
%  A - Input matrix
%  Q - Matrix - Unitary basis vectors of columns of A
%  R - Upper triangular matrix, A = Q*R

    [m, n] = size(A);
    Q = zeros(m, n);
    R = zeros(m, n);
    for k = 1:n
        R(k,k) = norm(A(1:m,k));    % start with normalized
        Q(1:m,k) = A(1:m,k)/R(k,k); % column vector
        for j = k + 1:n
            R(k,j) = Q(1:m,k)'*A(1:m, j);
            A(1:m,j) = A(1:m,j) - Q(1:m,k)*R(k,j);
        end
    end

13.3. Finding Orthogonal Basis Vectors¶

13.3.1. The Gram–Schmidt Algorithm¶

13.3.2. Implementation of Classic Gram–Schmidt¶

13.3.3. Implementation of Modified Gram–Schmidt¶

Table of Contents

Previous topic

Next topic

This Page