矩阵分析之Givens Reduction

1.Givens Rotation介绍

2.Givens Reduction算法原理

(1).获取待约减的矩阵AA: (m−by−n)(m-by-n)

⎡⎣⎢⎢⎢a11a21⋯am1a12a22⋯am2a13a23⋯am3⋯⋯⋯⋯a1na2n⋯amn⎤⎦⎥⎥⎥

\left [\begin{array}{}a_{11}& a_{12}& a_{13}&\cdots &a_{1n}\\a_{21}& a_{22}& a_{23}&\cdots &a_{2n}\\\cdots & \cdots &\cdots &\cdots &\cdots \\a_{m1}& a_{m2}& a_{m3}&\cdots&a_{mn}\\\end{array} \right]
(2)初始化矩阵P为单位矩阵
(3) k=1:n k=1 :n 循环下面操作：
利用旋转矩阵是的矩阵 AA中的第kk列的 kk个元素下面的所有元素为0.则j=k+1j=k+1取到 mm
(3).a计算矩阵Pk=P∗AP_k=P*A;
(3).b计算 c和sc和s

c=Ak[k][k]sqrt(Ak[k][k]2)+Ak[j][k]2

c=\frac{A_k[k][k]}{sqrt({A_k[k][k]}^2)+{A_k[j][k]}^2}

s=Ak[j][k]sqrt(Ak[k][k]2)+Ak[j][k]2

s=\frac{A_k[j][k]}{sqrt({A_k[k][k]}^2)+{A_k[j][k]}^2}
(3).c 令 P¯¯¯=I\overline P=I；然后用 cc去替代矩阵P¯¯¯\overline P中的第 kk行第 kk列的元素；然后用 ss去替代矩阵P¯¯¯\overline P中的第 kk行第 jj列的元素;然后用 −s-s去替代矩阵 P¯¯¯\overline P中的第 jj行第 jj列的元素;然后用 cc去替代矩阵P¯¯¯\overline P中的第 jj行第 jj列的元素;
(3).d,更新矩阵 PP，

P=P¯¯¯∗P

P=\overline P*P
(3).e,更新矩阵 AkA_k，

Ak=P∗A

A_k=P*A
(4).循环结束则可得到矩阵 PP和T=P∗AT=P*A

3.具体实例

给定矩阵A<script type="math/tex" id="MathJax-Element-3101">A</script>,求其Givens Reduction?
实现过程如下：

/*****************************************************
Description:Matrix_GivensReduction
Author:Robert.Tianyi
Date:2016.11.3
******************************************************/
#include<iostream>
#include<iomanip>
#include<cmath>
#include<cstdlib>
using namespace std;
float * MatrixMulMatrix(float *Pone,float *Rone,int m,int n);
float  VectorMUL(float a[],float b[],int len);
void Matrix_GivensReduction(float *B,int m,int n);
int main(){float A[3][3]={{0,-20,-14},{3,27,-4},{4,11,-2}};float B[9];cout<<"原始矩阵A："<<endl;for(int i=0;i<3;i++){for(int j=0;j<3;j++){cout<<setw(6)<<A[i][j]<<" ";B[i*3+j]=A[i][j];}cout<<endl;}Matrix_GivensReduction(B,3,3);
}
void Matrix_GivensReduction(float *B,int m,int n){cout<<"进行Matrix_GivensReduction:"<<endl; float A[m][n];for(int i=0;i<m;i++)for(int j=0;j<n;j++)A[i][j]=B[i*n+j];cout<<"矩阵A:"<<endl;for(int i=0;i<m;i++){for(int j=0;j<n;j++){cout<<setw(6)<<A[i][j]<<" ";    }cout<<endl;}float P[m][m];for(int i=0;i<m;i++){for(int j=0;j<m;j++){P[i][j]=0;if(i==j)P[i][j]=1;}}float Ak[m][n];for(int k=0;k<n;k++){for(int j=k+1;j<m;j++){//让第K列的第K个元素后面的元素进行Givens rotation //计算矩阵P*A 等于zAk矩阵 float *Pone=new float[m*m]; float *Aone=new float[m*n];int t=0;for(int i=0;i<m;i++){for(int j=0;j<m;j++){Pone[t++]=P[i][j];}}t=0;for(int i=0;i<m;i++){for(int j=0;j<n;j++){Aone[t++]=A[i][j];}}float *PmulA=new float [m*n];t=0;PmulA=MatrixMulMatrix(Pone,Aone,m,n);for(int i=0;i<m;i++){for(int j=0;j<n;j++){Ak[i][j]=PmulA[t++];}} float c,s;c= Ak[k][k]/sqrt((Ak[k][k])*(Ak[k][k])+(Ak[j][k])*(Ak[j][k]));s= Ak[j][k]/sqrt((Ak[k][k])*(Ak[k][k])+(Ak[j][k])*(Ak[j][k]));float temp_P[m][m];for(int i=0;i<m;i++){for(int j=0;j<m;j++){temp_P[i][j]=0;if(i==j)temp_P[i][j]=1;}}temp_P[k][k]=c;temp_P[k][j]=s;temp_P[j][k]=-s;temp_P[j][j]=c;float *temp_Pone=new float[m*m]; t=0;for(int i=0;i<m;i++){for(int j=0;j<m;j++){temp_Pone[t++]=temp_P[i][j];}}//计算矩阵temp_P乘以P float *temp_PmulP=new float [m*m];temp_PmulP=MatrixMulMatrix(temp_Pone,Pone,m,m);//更新矩阵Pfor(int i=0;i<m;i++){for(int j=0;j<m;j++){P[i][j]=temp_PmulP[i*m+j];}} }   }cout<<"矩阵P：\n"<<endl;for(int i=0;i<m;i++){for(int j=0;j<m;j++){cout<<setw(8)<<P[i][j];}cout<<endl;} //最后更新矩阵Ak=P*A; float *Pone=new float[m*m]; float *Aone=new float[m*n];int t=0;for(int i=0;i<m;i++){for(int j=0;j<m;j++){Pone[t++]=P[i][j];}}t=0;for(int i=0;i<m;i++){for(int j=0;j<n;j++){Aone[t++]=A[i][j];}}float *PmulA=new float [m*n];t=0;PmulA=MatrixMulMatrix(Pone,Aone,m,n);for(int i=0;i<m;i++){for(int j=0;j<n;j++){Ak[i][j]=PmulA[t++];if(Ak[i][j]>-0.00001&&Ak[i][j]<0.00001)Ak[i][j]=0;}} cout<<"矩阵T=P*A：\n"<<endl;for(int i=0;i<m;i++){for(int j=0;j<n;j++){cout<<setw(5)<<Ak[i][j];}cout<<endl;}
}
float * MatrixMulMatrix(float *Pone,float *Aone,int m,int n){float *PmulA=new float [m*n];float result[m][n];float (*P)[m];P=(float (*)[m])malloc(sizeof(float)*(m)*(m));float (*A)[n];A=(float (*)[n])malloc(sizeof(float)*(m)*(n));int k=0;for(int i=0;i<m;i++){for(int j=0;j<m;j++){P[i][j]=Pone[k++];}}k=0;for(int i=0;i<m;i++){for(int j=0;j<n;j++){A[i][j]=Aone[k++];}}for(int i=0;i<m;i++){//读取P的第i行 float *Pvector=new float[m];for(int k=0;k<m;k++){Pvector[k]=P[i][k]; }for(int j=0;j<n;j++){//读取A第J列 float *Avector=new float[m];for(int k=0;k<m;k++){Avector[k]=A[k][j]; }result[i][j]=VectorMUL(Pvector,Avector,m);PmulA[i*m+j]=result[i][j];}   }return PmulA;
}
/*向量内积*/
float  VectorMUL(float a[],float b[],int len){float *p;float result=0;for(int i=0;i<len;i++)result+=a[i]*b[i];return result;
}

运行结果如下：