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

计算数学



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;      
}