song.yz@foxmail.com wechat: math-box

计算数学

修正的Gram-Schimdt正交化方法

令\[\left\{ {{{\bf{x}}_1}{\bf{,}}{{\bf{x}}_2}{\bf{,}} \cdots {\bf{,}}{{\bf{x}}_n}} \right\}\] 是p维向量子空间W的任意一组基,则子空间W的标准正交基可以通过Gram-Schimdt正交化构造, 即 $$ \begin{array}{l} {{\bf{p}}_1}{\bf{ = }}{{\bf{x}}_1}{\bf{,}}{{\bf{u}}_1}{\bf{ = }}\frac{{{{\bf{p}}_1}}}{{\left\| {\bf{p}} \right\|}}{\bf{ = }}\frac{{{{\bf{x}}_1}}}{{\left\| {{{\bf{x}}_1}} \right\|}}\\ {{\bf{p}}_k}{\bf{ = }}{{\bf{x}}_k}{\bf{ - }}\sum\limits_{i = 1}^{k = 1} {\left( {{{\bf{v}}_i}^H{{\bf{x}}_k}} \right){{\bf{u}}_i}{\bf{,}}{{\bf{u}}_k}{\bf{ = }}\frac{{{{\bf{p}}_k}}}{{\left\| {{{\bf{p}}_k}} \right\|}}} \end{array} $$ 对于超定的线性方程系数矩阵的QR分解可以通过Gram-Schimdt正交化方法来实现, 然而采用Gram-Schimdt正交化方法求解列正交矩阵Q时, 舍入误差较大,这在求解最小二乘法时候, 有时会不稳定。针对Gram-Schimdt正交化的缺点,下面给出修正的Gram-Schimdt正交化算法。 $$ \begin{array}{l} {R_{11}} = \left\| {{{\bf{a}}_1}} \right\|\\ {{\bf{q}}_1} = \frac{{{{\bf{q}}_1}}}{{{R_{11}}}} \end{array} $$ 对于k=2...n $$ \begin{array}{l} {R_{jk}} = {\bf{q}}_j^H{{\bf{a}}_k},j = 1, \cdots ,k - 1\\ {R_{kk}} = \left\| {{{\bf{a}}_k} - \sum\limits_{j = 1}^{k - 1} {{\bf{q}}{R_{jk}}} } \right\|\\ {\bf{q}} = \frac{{{{\bf{a}}_k} - \sum\limits_{j = 1}^{k - 1} {{\bf{q}}{R_{jk}}} }}{{{R_{kk}}}} \end{array} $$ 下面给出其C++版本的代码

void MGS(const MAT &A, MAT &Q, MAT &R)
 /*------------------------------------------ comment
 Author      : Song Yezhi
 Date        : 2015-7-28        
 ---------------------------------------------------
 Desciption  :  修正的Gram-Schmidt正交化方法
  *
 Post Script :
  *           分解结果为 Q是列正交的基向量,R为方阵
  *           注意分解结果与Householder变换是不同的
  *           关于分解的几何意义可以参看麻省工学院编写的代数教材
 parameters  :
  *     A----原矩阵
  *     Q----正交基组成的向量空间
  *     R----方阵
 Ref:   
        1. scientific computation in C#
        2. scientific computation in Fortran
           song.yz   Tsinghua Univ. Press ,2011 
 -------------------------------------------------*/
 {
     //获取A矩阵维数
     int M, N;
     M = A.size1();
     N = A.size2();    
     //为Q和R分配地址
     //Q = new MAT(M, N);
     //R = new MAT(N, N);
     VEC vtmp(M);
     double s=0.0;
     int i, j;
     for (i = 0; i < M; i++)
         s = s + A(i, 0) * A(i, 0);
     R(0, 0) = sqrt(s);
     for (i = 0; i < M; i++)
         Q(i, 0) = A(i, 0) / R(0, 0);
     for (int k = 1; k < N; k++)
     {
         for (j = 0; j <= k - 1; j++)
         {
             s=0;
             for (int p = 0; p < M; p++)
                 s = s + Q(p, j) * A(p, k);
             R(j, k) = s; 
         }
         for (i = 0; i < M; i++)
             vtmp(i) = A(i, k);
         for (j = 0; j <= k - 1; j++)
             for (i = 0; i < M; i++)
                 vtmp(i) = vtmp(i) - Q(i, j) * R(j, k);        
         R(k,k)=sqrt(vecdot(vtmp,vtmp));             
         for (i = 0; i < M; i++)
             Q(i, k) = vtmp(i) / R(k, k);
      } 
}