Classical Gram-Schmidt and Modified Gram-Schmidt

---

一、Classical Gram-Schmidt (CGS)

For vecotr set { x 1 , x 2 , ⋯   , x m } \{\boldsymbol{x}_{1},\boldsymbol{x}_{2},\cdots,\boldsymbol{x}_{m}\} {x1​,x2​,⋯,xm​},

There are two versions of the above algorithm implementation in matlab. (The first one is realized by myself, and the second one is modified by my supervisor)

• CGS Algorithm implementation (v1)

function [Q] = my_cgs(A)
% A is a n-by-n matrix
% Q is column orthonormal matrix. If matrix A is nonsigular, the size Q will
% be n-by-n.
V = zeros(size(A));
Q = zeros(size(A));
r = zeros(size(A));
V(:,1) = A(:,1);
Q(:,1) = V(:,1)/sqrt(V(:,1)'*V(:,1));
for j=2:size(A,2)
    proj_sum = zeros(size(A,1),1);
    V(:,j)=A(:,j);
    for k=1:j-1
        r(j,k)=V(:,j)'*Q(:,k);
        proj_sum = proj_sum  + r(j,k)*Q(:,k);
    end
    V(:,j)=V(:,j)-proj_sum;
    if V(:,j)'*V(:,j) == 0
        break
    end
    Q(:,j)=V(:,j)/sqrt(V(:,j)'*V(:,j));
end

• CGS Algorithm implementation (v2)

function [Q] = my_cgs_v2(A)
% A is a n-by-n matrix
% Q is column orthonormal matrix. If matrix A is nonsigular, the size Q will
% be n-by-n.
%V = zeros(size(A));
Q = zeros(size(A));
Q(:,1) = A(:,1)/sqrt(A(:,1)'*A(:,1));
for j=2:size(A,2)
    v = A(:,j);
    r = v'*Q(:,1:j-1); %
    q_j=v-Q(:,1:j-1)*r';
    if q_j'*q_j == 0
        break
    end
    Q(:,j)=q_j/norm(q_j);
end

二、Modified Gram-Schmidt (MGS)

For vecotr set { x 1 , x 2 , ⋯   , x m } \{\boldsymbol{x}_{1},\boldsymbol{x}_{2},\cdots,\boldsymbol{x}_{m}\} {x1​,x2​,⋯,xm​},

The above MGS algorithm is also implemented in matlab by myself.

• MGS Algorithm implementation

function [Q] = my_mgs(A)
% A is a n-by-n matrix
% Q is column orthonormal matrix. If matrix A is nonsigular, the size Q will
% be n-by-n.
V = zeros(size(A));
Q = zeros(size(A));
r = zeros(size(A));
V(:,1) = A(:,1);
Q(:,1) = V(:,1)/sqrt(V(:,1)'*V(:,1));

for j=2:size(A,2)
    %proj_sum = zeros(size(A,1),1);
    V(:,j)=A(:,j);
    for k=1:j-1
        r(j,k)=V(:,j)'*Q(:,k);
        %proj_sum = proj_sum  + r(j,k)*Q(:,k);
        V(:,j)=V(:,j)-r(j,k)*Q(:,k);
    end
    %V(:,j)=V(:,j)-proj_sum;
    if V(:,j)'*V(:,j) == 0
        break
    end
    Q(:,j)=V(:,j)/sqrt(V(:,j)'*V(:,j));
end

三、Differences in Error Accumulation and Elimination

If an error is made in computing ${\mathbf{q}}_2$, so that $({\mathbf{q}}_1,{\mathbf{q}}_2)=\delta$ is small, but not-zero.

• In CGS, this error will not be corrected for in any of the computations that follow.
  Such as, $${\hat{\mathbf{q}}}_3={\mathbf{x}}_3-({\mathbf{x}}_3,{\mathbf{q}}_1){\mathbf{q}}_1-({\mathbf{x}}_3,{\mathbf{q}}_2){\mathbf{q}}_2$$
  then
  $$({\mathbf{q}}_1,\widehat{{\mathbf{q}}_3})=({\mathbf{q}}_1,{\mathbf{x}}_3)-({\mathbf{x}}_3,{\mathbf{q}}_1)-({\mathbf{x}}_3,{\mathbf{q}}_2)\delta =-({\mathbf{x}}_3,{\mathbf{q}}_2)\delta$$
  $$({\mathbf{q}}_2,\widehat{{\mathbf{q}}_3})=({\mathbf{q}}_2,{\mathbf{x}}_3)-({\mathbf{x}}_3,{\mathbf{q}}_1)\delta -({\mathbf{x}}_3,{\mathbf{q}}_2)=-({\mathbf{x}}_3,{\mathbf{q}}_1)\delta$$
  that is, $\widehat{{\mathbf{q}}_3}$ is not orthogonal to ${\mathbf{q}}_2$ and ${\mathbf{q}}_1$.

• In MGS, this error can be eliminated in ${\mathbf{q}}_2^T{\mathbf{q}}_3$ and the follow. As,
  $$\begin{array}{rl}& {\widehat{{\mathbf{q}}_3}}^{(0)}={\mathbf{x}}_3;\\ & {\widehat{{\mathbf{q}}_3}}^{(1)}={\widehat{{\mathbf{q}}_3}}^{(0)}-({\mathbf{q}}_1,{\widehat{{\mathbf{q}}_3}}^{(0)}){\mathbf{q}}_1;\\ & \widehat{{\mathbf{q}}_3}={\widehat{{\mathbf{q}}_3}}^{(1)}-({\mathbf{q}}_2,{\widehat{{\mathbf{q}}_3}}^{(1)}){\mathbf{q}}_2;\end{array}$$

  • then for ${\mathbf{q}}_2$,
    $$({\mathbf{q}}_2,\widehat{{\mathbf{q}}_3})=({\mathbf{q}}_2,{\widehat{{\mathbf{q}}_3}}^{(1)})-({\mathbf{q}}_2,{\widehat{{\mathbf{q}}_3}}^{(1)})=0$$
    so, $\widehat{{\mathbf{q}}_3}$ perserves orthogonality to ${\mathbf{q}}_2$.
  • then for ${\mathbf{q}}_1$,
    $$\begin{array}{ll}({\mathbf{q}}_1,\widehat{{\mathbf{q}}_3})& =({\mathbf{q}}_1,{\widehat{{\mathbf{q}}_3}}^{(1)})-({\mathbf{q}}_2,{\widehat{{\mathbf{q}}_3}}^{(1)})\delta \\ & =-(({\mathbf{q}}_2,{\widehat{{\mathbf{q}}_3}}^{(0)})-({\mathbf{q}}_1,{\widehat{{\mathbf{q}}_3}}^{(0)})\delta )\delta \\ & =-({\mathbf{q}}_2,{\widehat{{\mathbf{q}}_3}}^{(0)})\delta +({\mathbf{q}}_1,{\widehat{{\mathbf{q}}_3}}^{(0)}){\delta }^2\end{array}$$
    since, $\delta$ is very small, ${\delta }^2$ is much smaller, then the error in $({\mathbf{q}}_1,\widehat{{\mathbf{q}}_3})$ is not worse than in CGS, but MGS has eliminated the errors in $({\mathbf{q}}_2,\widehat{{\mathbf{q}}_3})$.

---

References:
[1] 谷同祥. 迭代方法和预处理技术(上).
[2] 潘建瑜. 矩阵计算讲义.
[3] The Modified Gram-Schmidt Procedure. [https://www.math.uci.edu/~ttrogdon/105A/html/Lecture23.html]