Gram Schmidt

Suppose we have a set of linearly independent vectors $V = \{v_1, v_2, ..., v_n\}$ and we want to find a new set of orthonormal vectors $Q = \{q_1, q_2, ..., q_n\}$ that span the same space. We can use the Gram Schmidt (GS) process, commonly implemented in three ways:

• Initially let $Q = \{\emptyset\}$
• Pivot: Find the vector $v_i$ not in $Q$ with the largest $\ell^2-$norm

• Classical: Orthogonalize the next vector $v_i$ against all vectors previously added to $Q$. Then normalize $v$ and add it to $Q$.
• Modified: Normalize $v$ and add it to $Q$. Then orthogonalize all vectors not in $Q$ against $v$ .
• Complete: When adding $v$, do Classical twice, or Modified twice, or both once.

By orthogonalize v against u, I mean to subtract the $u$-component from $v$ ($v - $proj$_uv$). Each method is more numerically stable than the last, and the stability of complete Gram Schmidt is on par with the Householder reflections or Givens rotations. A good example is if you let your vectors be the columns of matrix $A_{i,j} = 1/(j+k^2) \;\; ; \;\; j = 1,2,...,100 \;\; , \;\; k = 1,2,...,50$. The maximum dot product of any two vectors in $Q$ for classical = $0.62$, modified $= 0.02$, and complete $= 3.47$e-16, found using the Matlab code below.

GS can be used to find the $QR$ decomposition of a matrix, $Q$ unitary and $R$ upper-triangular, which lets you solve the linear system $Ax=b$ more efficiently than directly computing $A^{-1}$.

    function [Q] = GramSchmidt(Q, type)
       
        EM = 2^-53;
        n = size(Q,2);

        for i = 1:n

            % PIVOT
            qdotq = sum(Q(:,i:n).^2,1);
            j = i-1 + find(qdotq == max(qdotq),1);
            Q(:,[i,j]) = Q(:,[j,i]);
            if norm(Q(:,i)) < EM, break, end

            % CLASSIC AND COMPLETE
            if ~strcmp(type,'modified')
                for j = 1:i-1
                    Q(:,i) = Q(:,i) - dot(Q(:,i),Q(:,j))*Q(:,j);
                end
                Q(:,i) = Q(:,i) / norm(Q(:,i));
            end

            % MODIFIED AND COMPLETE
            if ~strcmp(type,'classical')
                Q(:,i) = Q(:,i) / norm(Q(:,i));
                for j = i+1:n
                    Q(:,j) = Q(:,j) - dot(Q(:,j),Q(:,i))*Q(:,i);
                end
            end
        end
    end

___________________________________________________________________________________