13.3. Finding Orthogonal Basis Vectors

We sometimes need to find a set of orthogonal basis vectors for the columns of a matrix. Three such needs are for vector projections (Alternate Projection Equation), finding least squares solutions to rectangular systems of equations (Left-Divide Operator), and finding the eigenvalues of a matrix (Finding Eigenvalues and Eigenvectors, Iterative QR Algorithm).

One of three algorithms is 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 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 rearranging it into what is known as the modified Gram–Schmidt process (MGS) [GOLUB13]. MATLAB does not include functions for either the CGS or MGS algorithms.

QR decomposition is the second algorithm for finding orthogonal basis vectors (QR Decomposition). Like Gram–Schmidt, QR is a quick, direct calculation, and returns functionally the same results. However, it is more accurate than either Gram–Schmidt algorithm because it uses orthogonal factoring methods [DEMMEL97].

The third algorithm for finding orthogonal basis vectors is from the SVD. It is found from an iterative algorithm as described in Calculating the SVD with Unitary Transformations, and is thus slower to determine either Gram–Schmidt or QR. MATLAB provides a function called orth, which returns the \(\bf{U}\) matrix from the economy singular value decomposition, [U, S]= svd(A,’econ’) (Projection and the Economy SVD). The SVD-based algorithm yields a different set of vectors than either Gram–Schmidt or QR.

13.3.1. Gram–Schmidt

Although MATLAB does not provide a function that implements the Gram–Schmidt process, we will briefly describe it here as a reference for those who 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 is orthogonal to the projection vector. Here, we want to find the vector representing the error from the projection.

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

\[\begin{split}\mathbf{A} = \begin{bmatrix} \vertbar{} & \vertbar{} & {} & \vertbar{} \\ \bm{a}_1 & \bm{a}_2 & \cdots{} & \bm{a}_n \\ \vertbar{} & \vertbar{} & {} & \vertbar{} \end{bmatrix}\end{split}\]

We will first find a matrix, \(\bf{B}\), who’s columns are orthogonal and are formed from projections of the \(\bf{A}\) columns. 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{split}\begin{aligned} \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}_n} {\bm{b}_1^T\,\bm{b}_1}\,\bm{b}_1 - \frac{\bm{b}_2^T\, \bm{a}_n}{\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{aligned}\end{split}\]

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

13.3.2. Implementation of Classic Gram–Schmidt

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 algorithm in the mod_gram_schmidt function is adapted from a code segment in Golub and Van Loan’s MATRIX Computations text [GOLUB13]. In this implementation, an upper triangular \(\bf{R}\) matrix is also returned. As with the \(\bf{Q}\) and \(\bf{R}\) matrices returned from MATLAB’s qr function, \(\bf{Q}\) and \(\bf{R}\) are factors of \(\bf{A}\), such that \(\mathbf{A} = \mathbf{Q\,R}\). For all but poorly conditioned matrices, the \(\bf{Q}\) and \(\bf{R}\) matrices from QR are nearly the same as found from the modified Gram–Schmidt process.

function [Q, R] = mod_gram_schmidt(A)
% MOD_GRAM_SCMIDT - Modified Gram-Schmidt Process
%   Variation 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