content

• Preface
• Orthogonal Transformation and Least Squares
• Cholesky Decomposition
• The Square Root Free Cholesky Algorithm
• The Householder Transformation
• Modified Gram–Schmidt orthogonalization
• Solving nolinear Equation
• Newton-Raphson Method
• mulit-root iterative improvement
• Muhammad method
• numerical Optimization
• Successive parabolic interpolation
• Newton’s Method
• Quasi-Newton Methods(DFP)
• Quasi-Newton Methods(BFGS)
• Interpolation
• Lagrange interpolation
• Newton method
• Hermit interpolation
• data fitting
• Polynomial fitting ephemeris
• Chebyshev fitting ephemeris
• ODE
• Symplectic Geometric Algorithms for Hamiltonian Systems
• Precise integration method
• Appendix：Matrix and Linear algebra class design
• C++ matrix class design
• C# matrix class design
• golang's slice for linear algebra

修正的Gram-Schimdt正交化方法

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