《并行计算奇异值分解.docx》由会员分享，可在线阅读，更多相关《并行计算奇异值分解.docx（14页珍藏版）》请在三一文库上搜索。
1、并行计算奇异值分解-Jacobi旋转鉴于矩阵的奇异值分解SVD在工程领域的广泛应用（如数据压缩、噪声去除、数值分析等等，包括在 NLP领域的潜在 语义索引LSI核心操作也是SVD）,今天就详细介绍一种SVD的实现方法-Jacobi旋转法。跟其它SVD算法相比, Jacobi法精度高，虽然速度慢，但容易并行实现。一些链接并彳T JACOBI方法求解矩阵奇异值的研究。本文呈现的代码就是依据这篇论文写出来的。Jama包是用于基本线性代数运算的java包，提供矩阵的cholesky分解、LUD分解、QR分解、奇异值分解，以及 PCA中要用到的特征值分解，此外可以计算矩阵的乘除法、矩阵的范数和条件数、解
2、线性方程组等。在线SVD运算器。bluebit在线矩阵运算器，提供矩阵的各种运算。C+ Matrix library提供矩阵的加减乘除、求行列式、LU分解、求逆、求转置。本文的头两段程序就引用了这里面的 matrix.h 。基于双边Jacobi旋转的奇异值分解算法V是A的右奇异向量，也是的特征向量;U是A的左奇异向量，也是AA的特征向量。特别地，当A是对称矩阵的时候,AS A4,即U=V , U的列向量不仅是的特征向量，也是A的特征向量。这一点在 主成分分析中会用至上 对于正定的对称矩阵，奇异值等于特征值，奇异向量等于特征向量。U、V都是正交矩阵，满足矩阵的转置即为矩阵的逆。双边Jacobi方
3、法本来是用来求解对称矩阵的特征值和特征向量的，由于就是对称矩阵，求出的特征向量就求出了 A的右奇异值,的特征值开方后就是A的奇异值。一个Jacobi旋转矩阵J形如:l -sinO：cosO它就是在一个单位矩阵的基础上改变 p行q行p列q列四个交点上的值，J矩阵是一个标准正交矩阵。H=AiA当我们要对H和p列和q列进行正交变换时，就对H施行:在一次迭代过程当中需要对A的任意两列进行一次正交变换，迭代多次等效于施行迭代的终止条件是为对角矩阵，即原矩阵H进过多次的双边Jacobi旋转后非对角线元素全部变成了 0（实现计算当中不可能真的变为0 ,只要小于一个很小的数就可以认为是 0 了）每一个J都是标
4、准正交阵，所以也是标准正交阵，记为V也是H的特从此式也可以看出对称矩阵的左奇异向量和右奇异向量是相等的征向量。当特征值是0时，对应的Ui和Vi不用求，我们只需要U和V的前r列就可以恢复矩阵A 了（r是非0奇异值的 个数），这也正是奇异值分解可以进行数据压缩的原因+ View Code基于单边Jacobi 旋转的SVD 算法相对于双边，单边的计算量小，并且容易并行实现。单边Jacobi方法直接对原矩阵A进行单边正交旋转，A可以是任意矩阵同样每一轮的迭代都要使 A的任意两列都正交一次，迭代退出的条件是 B的任意两列都正交。单边Jacobi旋转有这么一个性质：旋转前若则旋转后依然是II端11叫I反之
5、亦然。把矩阵B每列的模长提取出来,/蚱，把小以人记为V,则+ View Code基于单边Jacobi旋转的SVD并行算法单边Jacobi之于双边Jacobi的一个不同就是每次旋转只影响到矩阵 A的两列，因经很多列对的正交旋转变换是可以并 行执行的。基本思想是将A按列分块，每一块分配到一个计算节点上，块内任意两列之间进行单边Jacobi正交变换，所有的计算节点可以同时进行。问题就是矩阵块内部列与列之间进行正交变换是不够的，我们需要所有矩阵块上的任意两列之间都 进行正交变换，这就需要计算节点之间交换矩阵的列。本文采用奇偶序列的方法，具体可参考文章。由于这次我用的是C版的MPI (MPI也有C+和F
6、ortan版的)，所以上面代码用到的 C+版的matrix.h 就不能用 了，需要自己写一个C版的matrix.h 。matrix.h#ifndef _MATRIX_H#define _MATRIX_H#include#include#include/初始化一个二维矩阵double* getMatrix(int rows,int columns)double *rect=(double*)calloc(rows,sizeof(double*);int i;for(i=0;irows;+i)recti=(double*)calloc(columns,sizeof(double); return
7、rect;/返回一个单位矩阵double* getIndentityMatrix(introws)double* IM=getMatrix(rows,rows);int i;for(i=0;irows;+i)IMii=1.0;return IM; /返回一个矩阵的副本double* copyMatrix(double* matrix,introws,int columns)double* rect=getMatrix(rows,columns);int i,j;for(i=0;irows;+i)for(j=0;jcolumns;+j) rectij=matrixij； return rect;
8、/从一个一维矩阵得到一个二维矩阵void getFromArray(double* matrix,int rows,int columns,double *arr) int i,j,k=0;for(i=0;irows;+i)for(j=0;jcolumns;+j) matrixij=arrk+; /打印二维矩阵void printMatrix(double* matrix,introws,int columns)int i,j;for(i=0;irows;+i)for(j=0;jcolumns;+j) printf(%-10付,matrixij); printf(n); /释放二维矩阵void
9、 freeMatrix(double* matrix,introws)int i;for(i=0;irows;+i) free(matrixi);free(matrix);/获取二维矩阵的某一行double* getRow(double *matrix,int rows,int columns,int index) assert(indexrows);double *rect=(double*)calloc(columns,sizeof(double); int i;for(i=0;icolumns;+i)recti=matrixindexi; return rect;/获取二维矩阵的某一列d
10、ouble* getColumn(double *matrix,int rows,int columns,int index) assert(indexcolumns);double *rect=(double*)calloc(rows,sizeof(double); int i;for(i=0;irows;+i)recti=matrixiindex; return rect;/设置二维矩阵的某一列void setColumn(double *matrix,int rows,int columns,int index,double *arr) assert(indexcolumns);int
11、i;for(i=0;irows;+i) matrixiindex=arri;/交换矩阵的某两列void exchangeColumn(double *matrix,int rows,int columns,int i,int j) assert(icolumns);assert(jcolumns); int row;for(row=0;rowrows;+row) double tmp=matrixrowi; matrixrowi=matrixrowj; matrixrowj=tmp;/得到矩阵的转置double* getTranspose(double *matrix,int rows,int
12、 columns) double *rect=getMatrix(columns,rows);int i,j;for(i=0;icolumns;+i)for(j=0;jrows;+j) rectij=matrixji;return rect;/计算两向量内积double vectorProduct(double *vector1,double *vector2,int len) double rect=0.0; int i; for(i=0;ilen;+i)rect+=vector1i*vector2i; return rect;/两个矩阵相乘double*matrixProduct(doub
13、le *matrix1,int rows1,int columns1,double *matrix2,intcolurdouble *rect=getMatrix(rows1,columns2);int i,j;for(i=0;irows1;+i)for(j=0;jcolumns2;+j)double *vec1=getRow(matrix1,rows1,columns1,i);double *vec2=getColumn(matrix2,columns1,columns2,j);rectij=vectorProduct(vec1,vec2,columns1);free(vec1); free
14、(vec2); return rect;/得到某一列元素的平方和double getColumnNorm(double* matrix,int rows,int columns,int index) assert(indexcolumns);double* vector=getColumn(matrix,rows,columns,index);double norm=vectorProduct(vector,vector,rows); free(vector); return norm;/打印向量void printVector(double* vector,intlen)int i;for(
15、i=0;ilen;+i)printf(%-15.8ft,vectori);printf(n);#endifsvd.c#includempi.h#includematrix.h#include#include#include/gcc 编译的时候需要加-Im选项#define THREASHOLD 1e-8#define ITERATION 20#define ROW 3 /每个计算节点上的矩阵块有 3行4列#define COL 3#define LINELEN 5*COL/矩阵文件每一行的长度int sign(double number) if(number0)return -1;elsere
16、turn 1;int myRank;/计算结点的序号int procSize; /计算结点的数目MPI_Status status;/ 存储状态变量/*从文件中读取矩阵*/void readFrom *matrix,int row,int col,char* file)FILE *fp;int len=col*10;char *buf=(char*)calloc(len,sizeof(char);if(fp=fopenfile,r)=NULL)perror(fopen);printf(%snfile);exit(1);int i,j;for(i=0;irow;+i)if(fgets(buf,l
17、en,fp)=NULL)fprintf(stderr,文件的行数小于矩阵需要的行数n);exit；char *seg=strtok(buf,t);double ele=atof(seg);matrixi0=ele;for(j=1;jcol;+j)if(seg=strtok(NULL,t)=NULL)fprintf(stderr,文件的列数小于矩阵需要的列数n);exit； ele=atof(seg); matrixij=ele; memset(buf,0x00,len);free(buf);fclose(fp);/*把矩阵写入文件*/void writeTo *matrix,int rows,
18、int columns,char* file)FILE *fp;if(fp=fopen(file,w)=NULL) perror(fopen);exit(1);fprintf(fp,%dt%dn,rows,columns);/ 文件首行记录矩阵的行数和int i,j; for(i=0;irows;+i)for(j=0;jcolumns;+j)fprintf(fp,%-10ft,matrixij); fprintf(fp,n);fclose(fp);/*把向量写入文件*/void vectorTo *vector,int len,char* file)FILE *fp;if(fp=fopen(f
19、ile,w)=NULL) perror(fopen);exit(1);int i;for(i=0;ilen;+i)fprintf(fp,%-10付”,vectori); fclose(fp);/*两个向量进行单边Jacobi正交变换*/void orthogonalVector(double *Ci,double *Cj,int len1,double *Vi,double *Vj,int double ele=vectorProduct(Ci,Cj,len1);if(fabs(ele)THREASHOLD) return;/如果两列已经正交，不需要进行变换，则返回 true*pass=0;d
20、ouble ele1=vectorProduct(Ci,Ci,len1);double ele2=vectorProduct(Cj,Cj,len1);double tao=(ele1-ele2)/(2*ele);double tan=sign(tao)/(fabs(tao)+sqrt(1+pow(tao,2);double cos=1/sqrt(1+pow(tan,2);double sin=cos*tan;int row;for(row=0;rowlen1;+row)double var1=Cirow*cos+Cjrow*sin;double var2=Cjrow*cos-Cirow*sin
21、;Cirow=var1;Cjrow=var2;for(row=0;rowlen2;+row)double var1=Virow*cos+Vjrow*sin;double var2=Vjrow*cos-Virow*sin;Virow=var1;Vjrow=var2;/*矩阵的两列进行单边Jacobi正交变换。V是方阵，行/列数为columns*/void orthogonal(double *matrix,int rows,int columns,int i,int j,int *pass,dou assert(ij);double* Ci=getColumn(matrix,rows,colum
22、ns,i);double* Cj=getColumn(matrix,rows,columns,j);double* Vi=getColumn(V,columns,columns,i);double* Vj=getColumn(V,columns,columns,j);orthogonalVector(Ci,Cj,rows,Vi,Vj,columns,pass);int row;for(row=0;rowrows;+row) matrixrowi=Cirow; matrixrowj=Cjrow;for(row=0;rowcolumns;+row) Vrowi=Virow;Vrowj=Vjrow;
23、free(Ci);free(Cj)； free(Vi); free(Vj);void normalize(double *A,int rows,int columns)double *sigular=(double*)calloc(columns,sizeof(double);int i,j;for(i=0;icolumns;+i)double *vector=getColumn(A,rows,columns,i);double norm=sqrt(vectorProduct(vector,vector,rows); sigulari=norm;charoutFileS7=S,X,.,m,a,
24、t,0;outFileS1=0+myRank;vectorTo);double *U=getMatrix(rows,columns);for(j=0;jcolumns;+j)if(sigularj=0)for(i=0;irows;+i) U皿=0; elsefor(i=0;irows;+i) Uij=Aij/sigularj;char outFileU7=U,X,.,m,a,t,0;outFileU1=0+myRank;writeTo);free(sigular);freeMatrix(U,rows);void main(int argc, char*argv口)MPI_Init(&argc,
25、&argv);MPI_Comm_rank(MPI_COMM_WORLD,&myRank);MPI_Comm_size(MPI_COMM_WORLD,&procSize);assert(myRank0)/*奇数次迭代后矩阵按列范数从大到小排列；偶数次迭代后矩阵按列范数从小到int pass=1;int allpass=0;/*每个计算节点上相邻两列进行单边 Jacobi变换*/int i;for(i=1;i0) recv=1;for(+j;jCOL;j+=2)/首列需要和上一个计算结点的最后一/不是第1个计算节点/需要从上一个计算节点接收最力/第j列与j-1列进行单边Jacobiorthogon
26、al(A,ROW,COL,j-1,j,&pass,V,totalColumn);exchangeColumn(A,ROW,COL,j-1,j);/exchangeColumn(V,totalColumn,COL,j-1,j);if(j=COL)/ 最后一列剩下了,if(myRankprocSize-1) / send=1;无法配对，它需要发送给下 如果不是最后1个计算节大 /需要把最后一列发给if(send)/把最后一列发给下一个计算节点double *lastColumnA=getColumn(A,ROW,COL,COL-1);double *lastColumnV=getColumn(V,
27、totalColumn,COL,COL-MPIfree(lastColumnA);free(lastColumnV);if(recv)/从上一个计算节点接收最后一列double* preColumnA=(double*)calloc(ROW,sizeof(double)double* preColumnV=(double*)calloc(totalColumn,sizeof/本行首列与上一个计算节点末列进行正交变换double* firstColumnA=getColumn(A,ROW,COL,0);double* firstColumnV=getColumn(V,totalColumn,CO
28、L,0);orthogonalVector(preColumnA,firstColumnA,ROW,preColum/ 把 preColumn 留下setColumn(A,ROW,COL,0,preColumnA);setColumn(V,totalColumn,COL,0,preColumnV);/把巾rstColumn 发送给上一个计算结点MPI_Send(firstColumnV,totalColumn,MPI_DOUBLE,myRank-1 free(preColumnA);free(preColumnV);free(firstColumnA);free(firstColumnV);
29、if(send)/把最后一列发给下一个计算节点后，下一个计算节点做完了需要接收 double* nextColumnA=(double*)calloc(ROW,sizeof(double double* nextColumnV=(double*)calloc(totalColumn,sizeo MPI_Recv(nextColumnA,ROW,MPI_DOUBLE,myRank+1,49,M MPI_Recv(nextColumnV,totalColumn,MPI_DOUBLE,myRank+ /才亚收到的那一列赋给本块的最后一列一 setColumn(A,ROW,COL,COL-1,next
30、ColumnA);setColumn(V,totalColumn,COL,COL-1,nextColumnV); free(nextColumnA); free(nextColumnV); MPI_Barrier(MPI_COMM_WORLD);/各个计算节点都完成一次迭MPI_Reduce(&pass,&allpass,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD); MPI_Bcast(&allpass,1,MPI_INT,0,MPI_COMM_WORLD);/ 把if(allpass=procSize) break;/退出迭代 if(myRank=0) printf(迭代次数:%dn,ITERATION-iteration-1); char outFileV7=V,X,.,m,a,t,0; outFileV1=0+myRank; writeTo); normalize(A,ROW,COL); freeMatrix(A,ROW); freeMatrix(V,totalColumn); MPI_Finalize();