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

并行計(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)。一些鏈接

/Article/CDMD-10285-1012286387.htm

并行JACOBI方法求解矩陣奇異值的研究。本文呈現(xiàn)的代碼就是依據(jù)這篇論文寫出來的。/javanumerics/jama/

Jama包是用于基本線性代數(shù)運(yùn)算的java包，提供矩陣的cholesky分解、LUD分解、QR分解、奇異值分解，以及PCA中要用到的特征值分解，此外可以計(jì)算矩陣的乘除法、矩陣的范數(shù)和條件數(shù)、解線性方程組等。http://users.telenet.be/paul.larmuseau/SVD.htm

在線SVD運(yùn)算器。http://www.bluebit.gr/matrix-calculator/

bluebit在線矩陣運(yùn)算器，提供矩陣的各種運(yùn)算。/Projects/Matrix/

C++

Matrixlibrary提供矩陣的加減乘除、求行列式、LU分解、求逆、求轉(zhuǎn)置。本文的頭兩段程序就引用了這里面的matrix.h?；陔p邊Jacobi旋轉(zhuǎn)的奇異值分解算法V是A的右奇異向量，也是的特征向量;U是A的左奇異向量，也是的特征向量。特別地，當(dāng)A是對稱矩陣的時(shí)候，=，即U=V，U的列向量不僅是的特征向量，也是A的特征向量。這一點(diǎn)在主成分分析中會用到。對于正定的對稱矩陣，奇異值等于特征值，奇異向量等于特征向量。U、V都是正交矩陣，滿足矩陣的轉(zhuǎn)置即為矩陣的逆。雙邊Jacobi方法本來是用來求解對稱矩陣的特征值和特征向量的，由于就是對稱矩陣，求出的特征向量就求出了A的右奇異值，的特征值開方后就是A的奇異值。一個Jacobi旋轉(zhuǎn)矩陣J形如：它就是在一個單位矩陣的基礎(chǔ)上改變p行q行p列q列四個交點(diǎn)上的值，J矩陣是一個標(biāo)準(zhǔn)正交矩陣。當(dāng)我們要對H和p列和q列進(jìn)行正交變換時(shí)，就對H施行：在一次迭代過程當(dāng)中需要對A的任意兩列進(jìn)行一次正交變換，迭代多次等效于施行迭代的終止條件是為對角矩陣，即原矩陣H進(jìn)過多次的雙邊Jacobi旋轉(zhuǎn)后非對角線元素全部變成了0（實(shí)現(xiàn)計(jì)算當(dāng)中不可能真的變?yōu)?，只要小于一個很小的數(shù)就可以認(rèn)為是0了）。每一個J都是標(biāo)準(zhǔn)正交陣，所以也是標(biāo)準(zhǔn)正交陣，記為V，則是，則。從此式也可以看出對稱矩陣的左奇異向量和右奇異向量是相等的。V也是H的特征向量。當(dāng)特征值是0時(shí)，對應(yīng)的Ui和Vi不用求，我們只需要U和V的前r列就可以恢復(fù)矩陣A了（r是非0奇異值的個數(shù)），這也正是奇異值分解可以進(jìn)行數(shù)據(jù)壓縮的原因。+ViewCode基于單邊Jacobi旋轉(zhuǎn)的SVD算法相對于雙邊，單邊的計(jì)算量小，并且容易并行實(shí)現(xiàn)。

單邊Jacobi方法直接對原矩陣A進(jìn)行單邊正交旋轉(zhuǎn)，A可以是任意矩陣。

同樣每一輪的迭代都要使A的任意兩列都正交一次，迭代退出的條件是B的任意兩列都正交。svd.c#include"mpi.h"#include"matrix.h"#include<string.h>#include<stdlib.h>#include<math.h>

//gcc編譯的時(shí)候需要加-lm選項(xiàng)

#defineTHREASHOLD1e-8#defineITERATION20#define

ROW3

//每個計(jì)算節(jié)點(diǎn)上的矩陣塊有3行4列#defineCOL3#defineLINELEN5*COL

//矩陣文件每一行的長度

int

sign(double

number){

if(number<0)

return

-1;

else

return

1;}

int

myRank;

//計(jì)算結(jié)點(diǎn)的序號int

procSize;

//計(jì)算結(jié)點(diǎn)的數(shù)目MPI_Statusstatus;

//存儲狀態(tài)變量

/*從文件中讀取矩陣*/void

readFromFile(double

**matrix,int

row,int

col,char*file){

FILE

*fp;

int

len=col*10;

char

*buf=(char*)calloc(len,sizeof(char));

if((fp=fopen(file,"r"))==NULL){

perror("fopen");

printf("%s\n",file);

exit(1);

}

int

i,j;

for(i=0;i<row;++i){

if(fgets(buf,len,fp)==NULL){

fprintf(stderr,"文件的行數(shù)小于矩陣需要的行數(shù)\n");

exit(1);

}

char

*seg=strtok(buf,"\t");

double

ele=atof(seg);

matrix[i][0]=ele;

for(j=1;j<col;++j){

if((seg=strtok(NULL,"\t"))==NULL){

fprintf(stderr,"文件的列數(shù)小于矩陣需要的列數(shù)\n");

exit(1);

}

ele=atof(seg);

matrix[i][j]=ele;

}

memset(buf,0x00,len);

}

free(buf);

fclose(fp);}

/*把矩陣寫入文件*/void

writeToFile(double

**matrix,int

rows,int

columns,char*file){

FILE

*fp;

if((fp=fopen(file,"w"))==NULL){

perror("fopen");

exit(1);

}

fprintf(fp,"%d\t%d\n",rows,columns);

//文件首行記錄矩陣的行數(shù)和列數(shù)

int

i,j;

for(i=0;i<rows;++i){

for(j=0;j<columns;++j){

fprintf(fp,"%-10f\t",matrix[i][j]);

}

fprintf(fp,"\n");

}

fclose(fp);}

/*把向量寫入文件*/void

vectorToFile(double

*vector,int

len,char*file){

FILE

*fp;

if((fp=fopen(file,"w"))==NULL){

perror("fopen");

exit(1);

}

int

i;

for(i=0;i<len;++i){

fprintf(fp,"%-10f\t",vector[i]);

}

fclose(fp);}

/*兩個向量進(jìn)行單邊Jacobi正交變換*/void

orthogonalVector(double

*Ci,double

*Cj,int

len1,double

*Vi,double

*Vj,int

len2,int

*pass){

double

ele=vectorProduct(Ci,Cj,len1);

if(fabs(ele)<THREASHOLD)

return;

//如果兩列已經(jīng)正交，不需要進(jìn)行變換，則返回true

*pass=0;

double

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;row<len1;++row){

double

var1=Ci[row]*cos+Cj[row]*sin;

double

var2=Cj[row]*cos-Ci[row]*sin;

Ci[row]=var1;

Cj[row]=var2;

}

for(row=0;row<len2;++row){

double

var1=Vi[row]*cos+Vj[row]*sin;

double

var2=Vj[row]*cos-Vi[row]*sin;

Vi[row]=var1;

Vj[row]=var2;

}}

/*矩陣的兩列進(jìn)行單邊Jacobi正交變換。V是方陣，行/列數(shù)為columns*/void

orthogonal(double

**matrix,int

rows,int

columns,int

i,int

j,int

*pass,double

**V){

assert(i<j);

double*Ci=getColumn(matrix,rows,columns,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;row<rows;++row){

matrix[row][i]=Ci[row];

matrix[row][j]=Cj[row];

}

for(row=0;row<columns;++row){

V[row][i]=Vi[row];

V[row][j]=Vj[row];

}

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;i<columns;++i){

double

*vector=getColumn(A,rows,columns,i);

double

norm=sqrt(vectorProduct(vector,vector,rows));

sigular[i]=norm;

}

char

outFileS[7]={'S','X','.','m','a','t','\0'};

outFileS[1]='0'+myRank;

vectorToFile(sigular,columns,outFileS);

double

**U=getMatrix(rows,columns);

for(j=0;j<columns;++j){

if(sigular[j]==0)

for(i=0;i<rows;++i)

U[i][j]=0;

else

for(i=0;i<rows;++i)

U[i][j]=A[i][j]/sigular[j];

}

char

outFileU[7]={'U','X','.','m','a','t','\0'};

outFileU[1]='0'+myRank;

writeToFile(U,rows,columns,outFileU);

free(sigular);

freeMatrix(U,rows);}

void

main(int

argc,char

*argv[]){

MPI_Init(&argc,&argv);

MPI_Comm_rank(MPI_COMM_WORLD,&myRank);

MPI_Comm_size(MPI_COMM_WORLD,&procSize);

assert(myRank<10);

int

totalColumn=COL*procSize;

//算出原矩陣總共有多少列

/*準(zhǔn)備矩陣塊A和V*/

char

matrixFile[11]={'b','l','o','c','k','X','.','m','a','t','\0'};

matrixFile[5]='0'+myRank;

double

**A=getMatrix(ROW,COL);

readFromFile(A,ROW,COL,matrixFile);

double

**V=getMatrix(totalColumn,COL);

int

j;

for(j=0;j<COL;++j){

V[COL*myRank+j][j]=1.0;

}

/*開始進(jìn)行單邊Jacobi旋轉(zhuǎn)迭代*/

int

iteration=ITERATION;

while(iteration-->0){

/*奇數(shù)次迭代后矩陣按列范數(shù)從大到小排列；偶數(shù)次迭代后矩陣按列范數(shù)從小到大排列*/

int

pass=1;

int

allpass=0;

/*每個計(jì)算節(jié)點(diǎn)上相鄰兩列進(jìn)行單邊Jacobi變換*/

int

i;

for(i=1;i<=totalColumn;++i){

//原矩陣有幾列就需要做幾輪的交換

int

j;

int

send=0,recv=0;

//send記錄是否要把本矩陣塊的最后一列發(fā)送給下一個計(jì)算結(jié)點(diǎn)；recv記錄是否要從上一個計(jì)算結(jié)點(diǎn)接收一列數(shù)據(jù)

int

mod1=i%2;

//余數(shù)為0時(shí)是奇序，否則為偶序

int

mod2=(myRank*COL)%2;

//判斷本塊的第1列(首列)是否為原矩陣的第奇數(shù)列，為0則是，為1則不是

if(mod1^mod2){

//融合了奇序和偶序的情況

j=0;

//首列可以和第二列進(jìn)行正交變換

}

else{

j=1;

//首列需要和上一個計(jì)算結(jié)點(diǎn)的最后一列進(jìn)行正交變換

if(myRank>0){

//不是第1個計(jì)算節(jié)點(diǎn)

recv=1;

//需要從上一個計(jì)算節(jié)點(diǎn)接收最后一列

}

}

for(++j;j<COL;j+=2){

//第j列與j-1列進(jìn)行單邊Jacobi正交變換

orthogonal(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){

//最后一列剩下了，無法配對，它需要發(fā)送給下一個計(jì)算節(jié)點(diǎn)

if(myRank<procSize-1){

//如果不是最后1個計(jì)算節(jié)點(diǎn)

send=1;

//需要把最后一列發(fā)給下一個計(jì)算節(jié)點(diǎn)

}

}

if(send){

//把最后一列發(fā)給下一個計(jì)算節(jié)點(diǎn)

double

*lastColumnA=getColumn(A,ROW,COL,COL-1);

double

*lastColumnV=getColumn(V,totalColumn,COL,COL-1);

MPI_Send(lastColumnA,ROW,MPI_DOUBLE,myRank+1,59,MPI_COMM_WORLD);

MPI_Send(lastColumnV,totalColumn,MPI_DOUBLE,myRank+1,60,MPI_COMM_WORLD);

free(lastColumnA);

free(lastColumnV);

}

if(recv){

//從上一個計(jì)算節(jié)點(diǎn)接收最后一列

double*preColumnA=(double*)calloc(ROW,sizeof(double));

double*preColumnV=(double*)calloc(totalColumn,sizeof(double));

MPI_Recv(preColumnA,ROW,MPI_DOUBLE,myRank-1,59,MPI_COMM_WORLD,&status);

MPI_Recv(preColumnV,totalColumn,MPI_DOUBLE,myRank-1,60,MPI_COMM_WORLD,&status);

//本行首列與上一個計(jì)算節(jié)點(diǎn)末列進(jìn)行正交變換

double*firstColumnA=getColumn(A,ROW,COL,0);

double*firstColumnV=getColumn(V,totalColumn,COL,0);

orthogonalVector(preColumnA,firstColumnA,ROW,preColumnV,firstColumnV,totalColumn,&pass);

//把preColumn留下

setColumn(A,ROW,COL,0,preColumnA);

setColumn(V,totalColumn,COL,0,preColumnV);

//把firstColumn發(fā)送給上一個計(jì)算結(jié)點(diǎn)

MPI_Send(firstColumnA,ROW,MPI_DOUBLE,myRank-1,49,MPI_COMM_WORLD);

MPI_Send(firstColumnV,totalColumn,MPI_DOUBLE,myRank-1,50,MPI_COMM_WORLD);

free(preColumnA);

free(preColumnV);

free(firstColumnA);

free(firstColumnV);

}

if(send){

//把最后一列發(fā)給下一個計(jì)算節(jié)點(diǎn)后，下一個計(jì)算節(jié)點(diǎn)做完了正交變換，又把一列發(fā)送回來了，現(xiàn)在需要接收

double*nextColumnA=(double*)calloc(ROW,sizeof(double));

double*nextColumnV=(doub