并行計(jì)算奇異值分解

1、并行計(jì)算奇異值分解-jacobi旋轉(zhuǎn)鑒于矩陣的奇異值分解svd在工程領(lǐng)域的廣泛應(yīng)用（如數(shù)據(jù)壓縮、噪聲去除、數(shù)值分析等等，包括在 nlp領(lǐng)域的潛在 語義索引lsi核心操作也是svd）,今天就詳細(xì)介紹一種svd的實(shí)現(xiàn)方法-jacobi旋轉(zhuǎn)法。跟其它svd算法相比, jacobi法精度高，雖然速度慢，但容易并行實(shí)現(xiàn)。一些鏈接并彳t jacobi方法求解矩陣奇異值的研究。本文呈現(xiàn)的代碼就是依據(jù)這篇論文寫出來的。jama包是用于基本線性代數(shù)運(yùn)算的java包，提供矩陣的cholesky分解、lud分解、qr分解、奇異值分解，以及 pca中要用到的特征值分解，此外可以計(jì)算矩陣的乘除法、矩陣的范數(shù)和條件數(shù)、解

2、線性方程組等。在線svd運(yùn)算器。bluebit在線矩陣運(yùn)算器，提供矩陣的各種運(yùn)算。c+ matrix library提供矩陣的加減乘除、求行列式、lu分解、求逆、求轉(zhuǎn)置。本文的頭兩段程序就引用了這里面的 matrix.h ?；陔p邊jacobi旋轉(zhuǎn)的奇異值分解算法v是a的右奇異向量，也是的特征向量;u是a的左奇異向量，也是aa的特征向量。特別地，當(dāng)a是對(duì)稱矩陣的時(shí)候,as a4,即u=v , u的列向量不僅是的特征向量，也是a的特征向量。這一點(diǎn)在 主成分分析中會(huì)用至上 對(duì)于正定的對(duì)稱矩陣，奇異值等于特征值，奇異向量等于特征向量。u、v都是正交矩陣，滿足矩陣的轉(zhuǎn)置即為矩陣的逆。雙邊jacobi方

3、法本來是用來求解對(duì)稱矩陣的特征值和特征向量的，由于就是對(duì)稱矩陣，求出的特征向量就求出了 a的右奇異值,的特征值開方后就是a的奇異值。一個(gè)jacobi旋轉(zhuǎn)矩陣j形如:l -sino：coso它就是在一個(gè)單位矩陣的基礎(chǔ)上改變 p行q行p列q列四個(gè)交點(diǎn)上的值，j矩陣是一個(gè)標(biāo)準(zhǔn)正交矩陣。h=aia當(dāng)我們要對(duì)h和p列和q列進(jìn)行正交變換時(shí)，就對(duì)h施行:在一次迭代過程當(dāng)中需要對(duì)a的任意兩列進(jìn)行一次正交變換，迭代多次等效于施行迭代的終止條件是為對(duì)角矩陣，即原矩陣h進(jìn)過多次的雙邊jacobi旋轉(zhuǎn)后非對(duì)角線元素全部變成了 0（實(shí)現(xiàn)計(jì)算當(dāng)中不可能真的變?yōu)? ,只要小于一個(gè)很小的數(shù)就可以認(rèn)為是 0 了）每一個(gè)j都是標(biāo)

4、準(zhǔn)正交陣，所以也是標(biāo)準(zhǔn)正交陣，記為v也是h的特從此式也可以看出對(duì)稱矩陣的左奇異向量和右奇異向量是相等的征向量。當(dāng)特征值是0時(shí)，對(duì)應(yīng)的ui和vi不用求，我們只需要u和v的前r列就可以恢復(fù)矩陣a 了（r是非0奇異值的 個(gè)數(shù)），這也正是奇異值分解可以進(jìn)行數(shù)據(jù)壓縮的原因+ view code基于單邊jacobi 旋轉(zhuǎn)的svd 算法相對(duì)于雙邊，單邊的計(jì)算量小，并且容易并行實(shí)現(xiàn)。單邊jacobi方法直接對(duì)原矩陣a進(jìn)行單邊正交旋轉(zhuǎn)，a可以是任意矩陣同樣每一輪的迭代都要使 a的任意兩列都正交一次，迭代退出的條件是 b的任意兩列都正交。單邊jacobi旋轉(zhuǎn)有這么一個(gè)性質(zhì)：旋轉(zhuǎn)前若則旋轉(zhuǎn)后依然是ii端11叫i反之

5、亦然。把矩陣b每列的模長提取出來,/蚱，把小以人記為v,則+ view code基于單邊jacobi旋轉(zhuǎn)的svd并行算法單邊jacobi之于雙邊jacobi的一個(gè)不同就是每次旋轉(zhuǎn)只影響到矩陣 a的兩列，因經(jīng)很多列對(duì)的正交旋轉(zhuǎn)變換是可以并 行執(zhí)行的?；舅枷胧菍按列分塊，每一塊分配到一個(gè)計(jì)算節(jié)點(diǎn)上，塊內(nèi)任意兩列之間進(jìn)行單邊jacobi正交變換，所有的計(jì)算節(jié)點(diǎn)可以同時(shí)進(jìn)行。問題就是矩陣塊內(nèi)部列與列之間進(jìn)行正交變換是不夠的，我們需要所有矩陣塊上的任意兩列之間都 進(jìn)行正交變換，這就需要計(jì)算節(jié)點(diǎn)之間交換矩陣的列。本文采用奇偶序列的方法，具體可參考文章。由于這次我用的是c版的mpi (mpi也有c+和f

6、ortan版的)，所以上面代碼用到的 c+版的matrix.h 就不能用 了，需要自己寫一個(gè)c版的matrix.h 。matrix.h#ifndef _matrix_h#define _matrix_h#include#include#include/初始化一個(gè)二維矩陣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;/返回一個(gè)單位矩陣double* getindentitymatrix(introws)double* im=getmatrix(rows,rows);int i;for(i=0;irows;+i)imii=1.0;return im; /返回一個(gè)矩陣的副本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、/從一個(gè)一維矩陣得到一個(gè)二維矩陣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;/設(shè)置二維矩陣的某一列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;/得到矩陣的轉(zhuǎn)置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;/計(jì)算兩向量內(nèi)積double vectorproduct(double *vector1,double *vector2,int len) double rect=0.0; int i; for(i=0;ilen;+i)rect+=vector1i*vector2i; return rect;/兩個(gè)矩陣相乘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 編譯的時(shí)候需要加-im選項(xiàng)#define threashold 1e-8#define iteration 20#define row 3 /每個(gè)計(jì)算節(jié)點(diǎn)上的矩陣塊有 3行4列#define col 3#define linelen 5*col/矩陣文件每一行的長度int sign(double number) if(number0)return -1;elsere

16、turn 1;int myrank;/計(jì)算結(jié)點(diǎn)的序號(hào)int procsize; /計(jì)算結(jié)點(diǎn)的數(shù)目mpi_status status;/ 存儲(chǔ)狀態(tài)變量/*從文件中讀取矩陣*/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,文件的行數(shù)小于矩陣需要的行數(shù)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,文件的列數(shù)小于矩陣需要的列數(shù)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);/ 文件首行記錄矩陣的行數(shù)和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);/*兩個(gè)向量進(jìn)行單邊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;/如果兩列已經(jīng)正交，不需要進(jìn)行變換，則返回 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;/*矩陣的兩列進(jìn)行單邊jacobi正交變換。v是方陣，行/列數(shù)為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)/*奇數(shù)次迭代后矩陣按列范數(shù)從大到小排列；偶數(shù)次迭代后矩陣按列范數(shù)從小到int pass=1;int allpass=0;/*每個(gè)計(jì)算節(jié)點(diǎn)上相鄰兩列進(jìn)行單邊 jacobi變換*/int i;for(i=1;i0) recv=1;for(+j;jcol;j+=2)/首列需要和上一個(gè)計(jì)算結(jié)點(diǎn)的最后一/不是第1個(gè)計(jì)算節(jié)點(diǎn)/需要從上一個(gè)計(jì)算節(jié)點(diǎn)接收最力/第j列與j-1列進(jìn)行單邊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;無法配對(duì)，它需要發(fā)送給下 如果不是最后1個(gè)計(jì)算節(jié)大 /需要把最后一列發(fā)給if(send)/把最后一列發(fā)給下一個(gè)計(jì)算節(jié)點(diǎn)double *lastcolumna=getcolumn(a,row,col,col-1);double *lastcolumnv=getcolumn(v,

27、totalcolumn,col,col-mpifree(lastcolumna);free(lastcolumnv);if(recv)/從上一個(gè)計(jì)算節(jié)點(diǎn)接收最后一列double* precolumna=(double*)calloc(row,sizeof(double)double* precolumnv=(double*)calloc(totalcolumn,sizeof/本行首列與上一個(gè)計(jì)算節(jié)點(diǎn)末列進(jìn)行正交變換double* firstcolumna=getcolumn(a,row,col,0);double* firstcolumnv=getcolumn(v,totalcolumn,col,0);orthogonalvector(precolumna,firstcolumna,row,precolum/ 把 precolumn 留下setcolumn(a,row,col,0,precolumna);setcolumn(v,totalcolumn,col,0,precolumnv);/把巾rstcolumn