这篇本来是想和PCA写在一起的，写了一上午终于写完了Household变换，先把结果贴一下，首先说下Household的基本目的：就是把一个向量A=[a1,a2,a3.....an]转化成另一个有同样维度的向量B，其中B和A拥有相同的模，||A||==||B||。而B可以根据你自己的需要设计成任意模式比如B(b1,b2,0,0...0)或者其他形式。
而求矩阵特征值要做的目的是把一个矩阵转换成上三角矩阵，如下图，即把A转换成Q*R而R的形式如下图所示
过程：A为原始矩阵。
step1:
step2:
step3，step4:同step2，最终结果为一个上三角矩阵。
最后我们得到Q=H4*H3*H2*H1；R=最后的上三角矩阵，也等于Q*A；
关键问题是如何构建矩阵H
H的构建就是上述的Household 变换：
构造H的方法：
C语言程序：
//
// main.c
// Matrix
//
// Created by 谭升 on 14/12/16.
// Copyright (c) 2014年 谭升. All rights reserved.
//
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define DataType double
struct intM_{
int rowNum;
int colNum;
DataType *data;
};
typedef struct intM_ Mat;
void FreeMatrix(Mat *src){
free(src->data);
}
void MakeMatrix(Mat *src){
DataType *data=(DataType *)malloc(sizeof(DataType)*src->rowNum*src->colNum);
if(data==NULL){
printf("malloc memory fail\n");
exit(0);
}
src->data=data;
}
void ZeroMat(Mat *src){
for(int i=0;i<src->rowNum*src->colNum;i++)
src->data[i]=0;
}
void InitMatrix(Mat *src,int rowNum,int colNum){
src->colNum=colNum;
src->rowNum=rowNum;
MakeMatrix(src);
ZeroMat(src);
}
void MatrixCopy(Mat *src,Mat *dst){
if(src->rowNum!=dst->rowNum||
src->colNum!=dst->colNum){
printf("MatrixCopy wrong input !\n");
exit(0);
}
for(int i=0;i<src->rowNum*src->colNum;i++)
dst->data[i]=src->data[i];
}
void MatrixMulty(Mat * src1,Mat *src2,Mat *dst){
if((src1->colNum!=src2->rowNum)||
(dst->colNum!=src2->colNum) ||
(dst->rowNum!=src1->rowNum)){
printf("wrong parameter for matrix multy\n");
exit(0);
}
Mat temp;
InitMatrix(&temp, dst->rowNum, dst->colNum);
int rowNum=dst->rowNum;
int colNum=dst->colNum;
int step=src1->colNum;
for(int i=0;i<rowNum;i++)
for(int j=0;j<colNum;j++)
for(int k=0;k<step;k++){
temp.data[i*colNum+j]+=src1->data[i*step+k]*src2->data[k*colNum+j];
}
MatrixCopy(&temp, dst);
FreeMatrix(&temp);
}
void MatrixTranspose(Mat *src,Mat *dst){
Mat temp;
InitMatrix(&temp, dst->rowNum, dst->colNum);
if(src->rowNum!=dst->colNum||
src->colNum!=dst->rowNum){
printf("Transpose wrong parameter\n");
exit(0);
}
for(int i=0;i<src->rowNum;i++)
for(int j=0;j<src->colNum;j++){
temp.data[j*dst->colNum+i]=src->data[i*src->colNum+j];
}
MatrixCopy(&temp, dst);
FreeMatrix(&temp);
}
void MatrixEye(Mat * src){
int row=src->rowNum;
int col=src->colNum;
ZeroMat(src);
for(int i=0;i<row;i++)
src->data[i*col+i]=1;
}
void MatrixSub(Mat *src1,Mat *src2,Mat *dst){
if(src1->rowNum!=src2->rowNum||
src1->colNum!=src2->colNum){
printf("MatrixSub wrong input !\n");
exit(0);
}
for(int i=0;i<src1->rowNum*src1->colNum;i++)
dst->data[i]=src1->data[i]-src2->data[i];
}
void MatrixAdd(Mat *src1,Mat *src2,Mat *dst){
if(src1->rowNum!=src2->rowNum||
src1->colNum!=src2->colNum){
printf("MatrixAdd wrong input !\n");
exit(0);
}
for(int i=0;i<src1->rowNum*src1->colNum;i++)
dst->data[i]=src1->data[i]+src2->data[i];
}
void MatrixMulNum(Mat *src,double num,Mat *dst){
for(int i=0;i<src->rowNum*src->colNum;i++)
dst->data[i]=src->data[i]*num;
}
double VectorNorm1(Mat *src){
if(src->rowNum!=1&&src->colNum!=1){
printf("vector's Norm shoule input a vector\n");
exit(0);
}
Mat res;
Mat vt;
InitMatrix(&res, 1, 1);
InitMatrix(&vt, src->colNum, src->rowNum);
MatrixTranspose(src, &vt);
if(src->colNum==1)
MatrixMulty(&vt, src, &res);
if(src->rowNum==1)
MatrixMulty( src,&vt, &res);
DataType returnresult=res.data[0];
FreeMatrix(&res);
FreeMatrix(&vt);
return sqrt(returnresult);
}
double VectorNorm2(Mat *src){
if(src->rowNum!=1&&src->colNum!=1){
printf("vector's Norm shoule input a vector\n");
exit(0);
}
Mat res;
Mat vt;
InitMatrix(&res, 1, 1);
InitMatrix(&vt, src->colNum, src->rowNum);
MatrixTranspose(src, &vt);
if(src->colNum==1)
MatrixMulty(&vt, src, &res);
if(src->rowNum==1)
MatrixMulty( src,&vt, &res);
DataType returnresult=res.data[0];
FreeMatrix(&res);
FreeMatrix(&vt);
return returnresult;
}
void House_U(Mat * src,int Trans_col,Mat *u){
if(u->rowNum!=src->rowNum||u->colNum!=1){
printf("wrong input u\n");
exit(0);
}
Mat trans_vector;
InitMatrix(&trans_vector, src->rowNum-Trans_col, 1);
for(int i=Trans_col;i<src->rowNum;i++){
trans_vector.data[i-Trans_col]=src->data[i*src->colNum+Trans_col];
}
double xn=src->data[Trans_col*src->colNum+Trans_col];
DataType deta=0.0;
if(xn<0)
deta=-1.0*VectorNorm1(&trans_vector);
else if(xn>=0)
deta=VectorNorm1(&trans_vector);
for(int i=0;i<src->rowNum;i++){
if(i<Trans_col)
u->data[i]=src->data[i*src->colNum+Trans_col];
else if(i==Trans_col)
u->data[i]=deta;
else
u->data[i]=0;
}
for(int i=0;i<src->rowNum;i++){
u->data[i]=src->data[i*src->colNum+Trans_col]-u->data[i];
}
FreeMatrix(&trans_vector);
}
void House_P(Mat * u,Mat *P){
Mat E;
InitMatrix(&E, P->rowNum, P->colNum);
MatrixEye(&E);
Mat normu;
InitMatrix(&normu, 1, 1);
double unorm_1=2.0/VectorNorm2(u);
Mat u_tran;
InitMatrix(&u_tran, u->colNum, u->rowNum);
MatrixTranspose(u, &u_tran);
MatrixMulty(u, &u_tran, P);
MatrixMulNum(P, unorm_1, P);
MatrixSub(&E, P, P);
FreeMatrix(&E);
FreeMatrix(&normu);
FreeMatrix(&u_tran);
}
void House_nstep(Mat *src,Mat *pn,int Trans_col){
Mat pn_temp;
InitMatrix(&pn_temp, pn->rowNum,pn->colNum);
Mat u;
InitMatrix(&u, src->rowNum, 1);
House_U(src,Trans_col,&u);
House_P(&u,&pn_temp);
MatrixCopy(&pn_temp, pn);
FreeMatrix(&u);
FreeMatrix(&pn_temp);
}
void HouseholderTrans(Mat * src,Mat *dst_Q,Mat *dst_R){
int row=src->rowNum;
int col=src->colNum;
Mat p_temp;
Mat q_temp;
Mat src_temp;
InitMatrix(&p_temp, row, col);
InitMatrix(&q_temp, row, col);
InitMatrix(&src_temp, row, col);
MatrixCopy(src, &src_temp);
for(int i=0;i<col-1;i++){
House_nstep(&src_temp, &p_temp, i);
if(i==0)
MatrixCopy(&p_temp, &q_temp);
else
MatrixMulty(&p_temp, &q_temp, &q_temp);
MatrixMulty(&p_temp, &src_temp, &src_temp);
}
MatrixCopy(&q_temp,dst_Q);
MatrixTranspose(dst_Q, dst_Q);
MatrixCopy(&src_temp,dst_R);
//MatrixMulty(dst_Q, src, dst_R);
FreeMatrix(&q_temp);
FreeMatrix(&src_temp);
FreeMatrix(&p_temp);
}
int main() {
Mat A;
A.rowNum=4;
A.colNum=4;
DataType a[16]={1,3,3,4,
3,1,2,3,
3,2,1,5,
4,3,5,1,
};
A.data=a;
Mat Q;
InitMatrix(&Q, 4, 4);
Mat R;
InitMatrix(&R, 4, 4);
HouseholderTrans(&A, &Q, &R);
//MatrixMulty(&A, &A, &A);
printf("A=[\n");
for(int i=0;i<A.rowNum;i++){
for(int j=0;j<A.colNum;j++)
printf("%10lf",A.data[i*A.colNum+j]);
printf("\n");
}
printf("]\nQ=[\n");
for(int i=0;i<Q.rowNum;i++){
for(int j=0;j<Q.colNum;j++)
printf("%10lf",Q.data[i*Q.colNum+j]);
printf("\n");
}
printf("]\nR=[\n");
for(int i=0;i<R.rowNum;i++){
for(int j=0;j<R.colNum;j++)
printf("%10lf",R.data[i*R.colNum+j]);
printf("\n");
}
printf("]\n");
//FreeMatrix(&C);
}
结果：
