Why Modified Gram-Schmidt Is More Stable Than Classical – Numerical Linear Algebra

algorithmsgram-schmidtnumerical linear algebraorthogonality

This may be an old question, and there are certainly some related posts which I will mention below. However, there seems no clear answer to me yet. The question is: is there an intuitive way to explain why the modified Gram-Schmidt (MGS) process for doing QR factorization of a matrix $A\in\mathbb{C} ^{m\times n}$ gives a $Q$ matrix that is "more orthogonal" than that from the classical Gram-Schmidt (CGS) process? By "intuitive", I hope the explanation can be related to the procedural difference between MGS and CGS in a transparent way.

In Trefethen's Numerical Linear Algebra, the distinction between the CGS and MGS is as below:

At the $j$th step, both GS processes compute $q_j$ as
$$ q_j=\frac{P_j a_j }{\|| P_j a_j \|| } $$
while for CGS,
$$ P_j=I-Q_{j-1}Q_{j-1}^* $$
but for MGS,
$$ P_j=(I-q_{j-1}q_{j-1}^* )…(I-q_2q_2^* )(I-q_1q_1^* ) $$

Trefethen does not discuss why this procedural difference leads to the better numerical stability of MGS.

@AlgebraicPavel has given quantitative bounds here on the orthogonality factors: $\||I-Q^* Q\||\leq O(\epsilon \kappa(A))$ for MGS, while $\||I-Q^* Q\||\leq O(\epsilon \kappa^2(A))$ for CGS. These results are quantitative enough. However, as mentioned above, I would like a more intuitive reasoning for how this comes out.

@Ian said here that:

"Classical Gram-Schmidt, in which you subtract off the projections of the (k+1)th vector onto the first k vectors, is quite unstable, especially in high dimensions, because you essentially ensure that your new vector is orthogonal to the input vector in question but fail to ensure that the vectors you get at the end of the process are orthogonal to each other. Combine that with the fact that you can wind up subtracting nearly equal numbers and you get a bad situation."

This does sound like an intuitive and qualitative explanation for the problem of CGS. However, going into the detail, I don't feel comfortable about this line of reasoning. Specifically, saying that the "new vector is orthogonal to the input vector in question" does not seem to agree with what CGS is doing. For both CGS and MGS, the new vector ($a_j$) is subtracted in an attempt to make it orthogonal to the existing $q_i, i=1,…,j-1$. It might not be proper to call these $q_i$ "input vector", and this does not address the main procedural difference between MGS and CGS.

In this post, the $4\times 3$ Lauchli matrix is used as an example to demo the different results between MGS and CGS. Though there is still no intuitive explanation to the question either, I notice that for this Lauchli example, the result that $q_3^{CGS}$ fails to be orthogonal to $q_2^{CGS}$ is because the $r_{23}^{CGS}$ is wrongly computed, with a relative error of 100%. However, I can't figure out why the MGS procedure can alleviate this problem significantly.

I greatly appreciate any comments.

Best Answer

In both CGS and MGS, the orthogonalization step of subtracting off projections onto the columns of $Q$ that have already been computed introduces errors to due to finite-precision arithmetic. Each column $\mathbf{q}_i$ of $Q$ therefore has some error component in the direction of previously-computed columns $\{\mathbf{q}_1….\mathbf{q}_{i-1}\}$. The error accumulates for increasing column number $i$, which is an inherent weakness in both algorithms.

In CGS, the orthogonalization of a column $n$ against column $\mathbf{q}_{i}$ ($i<n$) is performed by projecting the original column of $A$ (call this $\mathbf{a}_n$) onto $\mathbf{q}_{i}$ and subtracting. $$ \begin{split} \mathbf{p}_{n} &\equiv \mathbf{a_n} - \sum_{i=1}^{n-1}(\mathbf{q_i^T}\cdot \mathbf{a_n})\mathbf{q_i} \\ \mathbf{q}_{n} &= \frac{\mathbf{p}_{n}}{\|\mathbf{p}_{n}\|} \end{split} $$ In MGS, on the other hand, the components along each $\mathbf{q}_i$ are immediately subtracted out of rest of the columns to the right of column $i$ as soon as the $\mathbf{q}_i$ are computed. Therefore the orthogonalization of column $n$ against $\mathbf{q}_{i}$ is not performed by projecting $\mathbf{q}_{i}$ against the original column of $A$ as it is in CGS, but rather against a vector obtained by subtracting from that column of $A$ the components in span($\mathbf{q}_1….\mathbf{q}_{i-1}$). This is important because of the error components of $\mathbf{q}_i$, which span $\{\mathbf{q}_1….\mathbf{q}_{i-1}\}$.

More precisely, in MGS the orthogonalization of column $n$ against $\mathbf{q}_{i}$ is performed by subtracting the component of $\mathbf{q}_{i}$ from the vector $\mathbf{v}_n^{i-1}$, where $\mathbf{v}_n^0\equiv \mathbf{a}_n$ and $\mathbf{v}_n^i$ ($0<i<n$) is defined as $$ \begin{split} \mathbf{v}_n^{i}&\equiv \mathbf{v}_n^{i-1} -(\mathbf{q}_{i}^T\cdot \mathbf{v}_n^{i-1})\mathbf{q}_{i}, \\ \mathbf{q}_n &= \frac{\mathbf{v}_n^{n-1}}{\|\mathbf{v}_n^{n-1}\|} \end{split} $$ Notice the difference in the projection factors in parenthesis in the above expression, $(\mathbf{q}_{i}^T\cdot \mathbf{v}_n^{i-1})$, and the corresponding one for CGS, ($\mathbf{q_i^T}\cdot \mathbf{a_n}$). The vector $\mathbf{q}_i$ has error components in span($\mathbf{q}_1….\mathbf{q}_{i-1}$) that will be introduce error into this projection factor. Whereas the vector $\mathbf{a}_n$ can in general have large components in span($\mathbf{q}_1….\mathbf{q}_{i-1}$), the vector $\mathbf{v}_n^{i-1}$ has only error components in span($\mathbf{q}_1….\mathbf{q}_{i-1}$) because in computing $\mathbf{v}_n^{i-1}$ those components of $\mathbf{a}_n$ in span($\mathbf{q}_1….\mathbf{q}_{i-1}$) have already been subtracted off. As a result, the error in this multiplicative factor due to the imperfect orthogonality between $\mathbf{q}_i$ and $\{\mathbf{q}_1...\mathbf{q}_{i-1}\}$ is much smaller in MGS than it is in CGS.

Because of the much smaller error in this projection factor, the MGS introduces less orthogonalization error at each subtraction step than does CGS.

Related Solutions

Algorithm for orthogonalizing polynomials with specific inner product

I can't tell from what you said what you implemented, this is an example called [Chebfun]. 1. It'd be better if you posted your code. Note this is some code for the Legendgre polynomials with pseudo-code.

x = chebfun('x',[-1 1]);
w = exp(pi*x);
N = 5;
P = OrthPoly(w,N);

    function P = OrthPoly(w,N)
        if isnumeric(w), w = chebfun(w,[-1 1]); end
        d = w.ends;                     % the domain
        x = chebfun('x',d);             % linear chebfun
        P = chebfun(1./sqrt(sum(w)),d); % the constant (normalised)
        for k = 1:N;
            xk = x.*P(:,k);
            P(:,k+1) = xk;
            for j = 1:k       % Subtract out the components
                C = sum(w.*xk.*P(:,j));
                P(:,k+1) = P(:,k+1) - C*P(:,j);
            end
            P(:,k+1) = P(:,k+1)./sqrt(sum(w.*P(:,k+1).^2)); % normalise
        end
    end

It's actually the lanczos iteration..

Note this is from Trefethen..

Connect the normal equation to the Gram-Schmidt method of finding an orthogonal vector

You might be having trouble showing the equivalence because you need additional assumptions for the projection matrix to be $QQ^T$ in the first place: it’s actually $Q(Q^TQ)^{-1}Q^T$. With the additional assumption that the columns of $Q$ are unit vectors, this reduces to $QQ^T$. If we remove this restriction, however, we’re no longer guaranteed that $Q^TQ$ is invertible: we need to add the condition that the columns of $Q$ are linearly independent. (Pairwise orthogonality of the columns isn’t enough since the zero vector is orthogonal to everything.) Equivalently, $Q$ must have full column rank.

Setting $Q=\begin{bmatrix}v_1&v_2&\cdots&v_n\end{bmatrix}$, we have by orthogonality of the columns $$Q^TQ = \begin{bmatrix}\langle v_1,v_1\rangle&0&\cdots&0\\0&\langle v_2,v_2\rangle&\cdots&0\\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&\langle v_n,v_n\rangle\end{bmatrix}$$ with all of the diagonal entries nonzero. We can then expand directly: $$Q(Q^TQ)^{-1}Q^Tb = Q(Q^TQ)^{-1}\begin{bmatrix}\langle v_1,b\rangle \\ \langle v_2,b\rangle \\ \vdots \\ \langle v_n,b\rangle\end{bmatrix} = Q\begin{bmatrix}{\langle v_1,b\rangle \over \langle v_1,v_1\rangle} \\ {\langle v_2,b\rangle \over \langle v_2,v_2\rangle} \\ \vdots \\ {\langle v_n,b\rangle \over \langle v_n,v_n\rangle} \end{bmatrix} = \sum_{j=1}^n {\langle v_j,b\rangle \over \langle v_j,v_j\rangle} v_j.$$

Best Answer

Related Solutions

Algorithm for orthogonalizing polynomials with specific inner product

Connect the normal equation to the Gram-Schmidt method of finding an orthogonal vector

Related Question