householder.c

Householder法求解最小二乘问题。可以避免常规方法遇到奇异矩阵(即行列式|A|接近零)时误差太大的问题。本方法的精度非常高。

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

void main()
{ 
	int Householder(double a[],int m,int n,double b[],double q[]);
	
	int i,j,m,n;
    static double a[30]={-544.20625,1,
      -990.1204,1,
    -1528.59745,1,
      -2169.914,1,
    -2909.32685,1,
     -3731.9989,1,
    -4620.26925,1,
    -5573.62405,1,
     -6610.6716,1,
    -7821.08335,1,
      -9217.586,1,
   -10790.71085,1,
    -12650.2696,1};
    static double b[13]={406.3421,
		485.4862,
		591.4679,
		691.1652,
		787.6605,
		857.6836,
		918.8571,
		987.8525,
		1086.2426,
		1334.5809,
		1458.4244,
		1687.8253,
		2031.2922};
	static double q[30];

    m=13;
	n=2;
    i=Householder(a,m,n,b,q);
    if(i!=0)
    { 
		for(i=0;i<n;i++)
			printf("x(%d)=%e\n",i,b[i]);
        printf("\n");
        printf("MAT Q IS:\n");
        for(i=0;i<m;i++)
			for(j=0;j<n;j++)
			{
				printf("%e  ",q[i*m+j]);
				printf("\n");
			}
        printf("\n");
        printf("MAT R IS:\n");
        for(i=0;i<m;i++)
	        for(j=0;j<m;j++)
			{
				printf("%e  ",a[i*m+j]);
				printf("\n");
			}
    }
	return;
}
/* 
求解AX=B的豪斯荷尔德(Householder)变换法

输入： 
  A是m*n的矩阵 
  B是m*1的向量 

返回值: 
  q是m*m的矩阵，存放A的QR分解式中的正交矩阵Q 
  A中返回矩阵R 
  B的前n个分量中返回最小二乘解X
*/
int Householder(double a[],int m,int n,double b[],double q[])
{ 
	int qr(double a[],int m,int n,double q[]);
	
	int i,j;
    double d,*c;
    
    c=malloc(n*sizeof(double));

	i=qr(a,m,n,q);
    if(i==0) 
	{ 
		free(c);
		return(0);
	}
    for(i=0;i<=n-1;i++)
    { 
		d=0.0;
        for(j=0;j<=m-1;j++)
			d=d+q[j*m+i]*b[j];
        c[i]=d;
    }
    b[n-1]=c[n-1]/a[n*n-1];
    for(i=n-2;i>=0;i--)
    { 
		d=0.0;
        for(j=i+1;j<=n-1;j++)
			d=d+a[i*n+j]*b[j];
        b[i]=(c[i]-d)/a[i*n+i];
    }
    free(c);
	return(1);
}

/* 矩阵的 QR 分解 */
int qr(double a[],int m,int n,double q[])
{ 
	int i,j,k,l,nn,p,jj;
    double u,alpha,w,t;
    
	if(m<n)
    { 
		printf("fail\n");
		return(0);
	}
    for(i=0;i<=m-1;i++)
		for(j=0;j<=m-1;j++)
		{ 
			l=i*m+j;
			q[l]=0.0;
			if(i==j) q[l]=1.0;
		}
    nn=n;
    if(m==n) nn=m-1;
    for(k=0;k<=nn-1;k++)
    { 
		u=0.0;
		l=k*n+k;
        for(i=k;i<=m-1;i++)
        { 
			w=fabs(a[i*n+k]);
            if(w>u) u=w;
        }
        alpha=0.0;
        for(i=k;i<=m-1;i++)
        { 
			t=a[i*n+k]/u;
			alpha=alpha+t*t;
		}
        if(a[l]>0.0) u=-u;
		
		alpha=u*sqrt(alpha);
        if(fabs(alpha)+1.0==1.0)
        { 
			printf("fail\n");
			return(0);
		}
        u=sqrt(2.0*alpha*(alpha-a[l]));
        if((u+1.0)!=1.0)
        { 
			a[l]=(a[l]-alpha)/u;
            for(i=k+1;i<=m-1;i++)
            { 
				p=i*n+k;
				a[p]=a[p]/u;
			}
            for(j=0;j<=m-1;j++)
            { 
				t=0.0;
                for(jj=k;jj<=m-1;jj++)
					t=t+a[jj*n+k]*q[jj*m+j];
                for(i=k;i<=m-1;i++)
                { 
					p=i*m+j;
					q[p]=q[p]-2.0*t*a[i*n+k];
				}
            }
            for(j=k+1;j<=n-1;j++)
            { 
				t=0.0;
                for(jj=k;jj<=m-1;jj++)
					t=t+a[jj*n+k]*a[jj*n+j];
                for(i=k;i<=m-1;i++)
                { 
					p=i*n+j;
					a[p]=a[p]-2.0*t*a[i*n+k];
				}
            }
            a[l]=alpha;
            for(i=k+1;i<=m-1;i++)
				a[i*n+k]=0.0;
        }
    }
    for(i=0;i<=m-2;i++)
		for(j=i+1;j<=m-1;j++)
		{ 
			p=i*m+j;
			l=j*m+i;
			t=q[p];
			q[p]=q[l];
			q[l]=t;
		}
    return(1);
}