用householder方法进行矩阵的QR分解

【算法原理】
[Householder阵]
（1） 设a Rn， = ||a||2，通常取 与a1同号，记H=I-2vvT，（v= ）,
则Ha= - e1. H=I -2vvT称为Householder阵。
（2） 更一般地，对a=(a1,a2,…am,am+1,…,an)T,记 =  ，可求出H，使
Ha=(a1,a2,…am,  ,0,…,0)T。
为此，先在Rn-m中求 使 满足
     =(am+1,…,an)T=(- ,0,…,0,0)T，
再作H= ,则Ha= (a1,a2,…am,am+1,…,an)T =( a1,a2,…am,- ,0,…,0,0)T

[用Householder方法求矩阵的QR分解]
记A=(aij)n*n，由1可知，存在H1=I -2v1v1T，使
H1(a11,a21,…,an1)T=(a11(1),0,…,0)T,
于是            H1A= 
又由1知，存在H2= ,使 ，于是
H1A=  = 
类似地依次进行n-1次，得出
Hn-1Hn-2…H1A= 。
记R=Hn-1Hn-2…H1A，Q=HnHn-1…H1,得A=Q*R


【函数表】
函数名    作用
Matix_Multiply（double A[n][n],double B[n][n],double C[n][n]）    计算A*B 并将结果返回到C中(计算H*A和Q*H)
Matix_copy(double A[n][n],double B[n][n])    用于矩阵的拷贝，在进行矩阵乘法时，需要一个过渡的变量矩阵，采用此函数进行矩阵的复制。
Householder_trans    Householder变换函数。
Matrix_print    矩阵打印函数。
Matrix_input    矩阵输入函数。
【源代码】
用vc++6.0编译通过
//利用Householder方法，进行矩阵的QR分解。                          
//作者：Kevin Mauchly [房智勇]                                     
//Date:2004/11/11
#include 
#include 
#define n 100 //maxN
void Matrix_print(double A[n][n])
{
    for (int i=0;i    {
        for (int j=0;j            printf("%8.4lf\t",A[i][j]);
        printf("\n");
    }
}
double Matrix_norm(double a[n]) 
{
    double d=0;
    for (int i=0;i        d+=a[i]*a[i];
    return sqrt(d);
}

void Matrix_multiply(double A[n][n],double B[n][n],double C[n][n])
{

    for (int i=0;i        for (int j=0;j        {
            C[i][j]=0;
            for (int t=0;t                C[i][j]+=A[i][t]*B[t][j];
        }

}
void Matrix_copy(double A[n][n],double B[n][n])
{
    for(int i=0;i        for(int j=0;j        A[i][j]=B[i][j];
}

void householder_trans(double A[n][n],int k,double Q[n][n])
{
    double a[n];
    for(int i=0;i        a[i]=0;
    for(i=n-k;i        a[i]=A[i][n-k];  
    a[n-k]+=Matrix_norm(a);
    double d=Matrix_norm(a);
    for(i=0;i        a[i]=a[i]/d;  
    double H[n][n];
    for(i=0;i    {
        for(int j=0;j            H[i][j]=-2*a[i]*a[j];
        H[i][i]++;
    }  //μ?H=I-2vvT
    double temp[n][n];
    Matrix_multiply(H,A,temp);
    Matrix_copy(A,temp);
    Matrix_multiply(Q,H,temp);
    Matrix_copy(Q,temp);

}
Matrix_input(double A[n][n])
{
   For(int i=0 ;i      For(int j=0;j        {
           Printf(“a[%d][%d]=”,i,j);
           Scanf(“%lf”,&a[i][j]);
}
}
void main()
{
    double Q[n][n];
    double A[n][n];
     int n;
     Matrix_input(A)
    printf("A: \n");
    Matrix_print(A);
    for (int i=0;i    {
        for (int j=0;j        Q[i][j]=0;
        Q[i][i]=1;
    }

    for (i=n;i>=2;i--) householder_trans(A,i,Q);
    printf("R: \n");
    Matrix_print(A);
    printf("Q: \n");
    Matrix_print(Q);
}
【DEBUG】
经检验，该算法结果可靠，可行。没有发现数值计算中的偏差。