矩阵的QR分解

这篇本来是想和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);
    
}
阅读更多
换一批

没有更多推荐了,返回首页