计算数学
Householder变换
对于给定的向量X与Y,二者的2范数相同，则可以找到正交变换H， 使得HX=Y。其变换就是镜像变换或称为Householder变换。取
\[{\bf{u = }}\frac{{{\bf{x - y}}}}{{{{\left\| {{\bf{x - y}}} \right\|}_2}}}\] 令 \[{\bf{H = I - 2u}}{{\bf{u}}^T}\] 可以验证H为正交矩阵，保持向量2范数不变。 对于单个向量X，欲用镜像变换使之 \[{\bf{Hx = w}} = \left\| {\bf{x}} \right\|{{\bf{e}}_1},{{\bf{e}}_1} = {\left( {1,0, \cdots } \right)^T}\] 可令 \[{\bf{u = w - x}}\] 继而可以构造变换矩阵 \[{\bf{H = I - }}2\frac{{{\bf{u}}{{\bf{u}}^T}}}{{\left\| {\bf{u}} \right\|_2^2}}\] 对于矩阵A做QR分解，第一次正交变换后，形如下 \[{{\bf{H}}_1}{\bf{A}} = \left[ {\begin{array}{*{20}{c}} *&*&*& \cdots &*\\ 0&*&*& \cdots &*\\ 0&*&*& \cdots &*\\ 0&*&*& \cdots &*\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0&*&*& \cdots &* \end{array}} \right]\] 第二次变换形如下 \[{{\bf{H}}_2}{{\bf{H}}_1}{\bf{A}} = \left[ {\begin{array}{*{20}{c}} *&*&*& \cdots &*\\ 0&*&*& \cdots &*\\ 0&0&*& \cdots &*\\ 0&0&*& \cdots &*\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0&0&*& \cdots &* \end{array}} \right]\] 如此循环到矩阵的最后一列，至此变换后的矩阵已经成为上三角矩阵。 由数值代数理论可知，如果R对角元素正负号都选正号或者负号，则QR分解是唯一的。计算H时候，并不需要存储各次的变换结果， 而是逐步矩阵累成而得，变换到矩阵的最后一列，新的矩阵已经是H。 上三角阵R也不需要再执行 ，当变换到最后一列时已经自动形成上三角阵， 该矩阵即为R 。
下面给出C++版本源代码，其中用到矩阵类，见附录部分。
void householder(const MAT &A,MAT &Q,MAT &R)
/*------------------------------------------------------
Versions and Changes :
v1.0---- 2015-09-02
householder 变换
REF：《C#科学计算讲义》
--------------------------------------------------------
Input Parameters :
A---要分解的矩阵
Output Parameters :
Q---正交矩阵
R---分解后的上三角矩阵
Notes :
--------------------------------------------------------
Author : 宋叶志
Copyrigt(C) : Shanghai Astronomical Observatory, CAS
(All rights reserved) 2015
-------------------------------------------------------*/
{
int M=A.size1();
int N=A.size2();
MAT H0(M,M);
MAT H1(M,M);
MAT H2(M,M);
MAT A1(M,N);
MAT A2(M,N);
VEC u(M);
int i,j;
matcopy( A1, A);
//copy matrix A to A1
for(i=0;i< M;i++)
for(j=0;j< M;j++)
H1(i,j)=0.0;
for(i=0;i< M;i++)
H1(i,i)=1.0;
int k;
for(k=0;k< N;k++)
// k is index of all col
{
//set the H to eye matrix
for(i=0;i< M;i++)
for(j=0;j< M;j++)
H0(i,j)=0.0;
for(j=0;j< M;j++)
H0(j,j)=1.0;
double s=0.0;
for(i=k;i< M;i++)
s=s+A1(i,k)*A1(i,k);
s=sqrt(s);
for(i=0;i< M;i++)
u(i)=0.0;
if (A(k,k)>=0.0)
u(k)=A1(k,k)+s;
else
u(k)=A1(k,k)-s;
for(i=k+1;i< M;i++)
u(i)=A1(i,k);
double du;
du=vecdot(u,u);
for(i=k;i< M;i++)
for(j=k;j< M;j++)
{
H0(i,j)=-2.0*u(i)*u(j)/du;
if(i==j)
H0(i,j)=1.0+H0(i,j);
}
matmul(A2,H0,A1);
matcopy(A1,A2);
//在 C#科学计算讲义中（song.yz）
//直接 matmul(H1,H1,H0)
//这里不行，因为 matmul中使用 const &
matmul(Q,H1,H0);
matcopy(H1,Q);
}
matcopy(Q,H1);
matcopy(R,A1);
//------set the lower part of matrix R to zero
// in fact ,it is close to zero after above process
// that means it is no necessary to set the part to zero actually
for(i=0;i< N;i++)
for(j=i+1;j< M;j++)
R(j,i)=0.0;
}