修正的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);
}
}