并行计算奇异值分解

并行计算奇异值分解--Jacobi旋转

鉴于矩阵的奇异值分解SVD在工程领域的广泛应用（如数据压缩、噪声去除、数值分析等等，包括在NLP领域的潜在语义索引LSI核心操作也是SVD），今天就详细介绍一种SVD的实现方法--Jacobi旋转法。跟其它SVD算法相比，Jacobi法精度高，虽然速度慢，但容易并行实现。

一些链接

/Article/CDMD-10285-1012286387.htm并行JACOBI方法求解矩阵奇异值的研究。本文呈现的代码就是依据这篇论文写出来的。

/javanumerics/jama/ Jama包是用于基本线性代数运算的java包，提供矩阵的cholesky 分解、LUD分解、QR分解、奇异值分解，以及PCA中要用到的特征值分解，此外可以计算矩阵的乘除法、矩阵的范数和条件数、解线性方程组等。

http://users.telenet.be/rmuseau/SVD.htm在线SVD运算器。

http://www.bluebit.gr/matrix-calculator/ bluebit在线矩阵运算器，提供矩阵的各种运算。

/Projects/Matrix/C++ Matrix library提供矩阵的加减乘除、求行列式、LU分解、求逆、求转置。本文的头两段程序就引用了这里面的matrix.h。

基于双边Jacobi旋转的奇异值分解算法

V是A的右奇异向量，也是的特征向量;

U是A的左奇异向量，也是的特征向量。

特别地，当A是对称矩阵的时候，=，即U=V，U的列向量不仅是的特征向量，也是A 的特征向量。这一点在主成分分析中会用到。

对于正定的对称矩阵，奇异值等于特征值，奇异向量等于特征向量。

U、V都是正交矩阵，满足矩阵的转置即为矩阵的逆。

双边Jacobi方法本来是用来求解对称矩阵的特征值和特征向量的，由于就是对称矩阵，求出的特征向量就求出了A的右奇异值，的特征值开方后就是A的奇异值。

一个Jacobi旋转矩阵J形如：

它就是在一个单位矩阵的基础上改变p行q行p列q列四个交点上的值，J矩阵是一个标准正交矩阵。

当我们要对H和p列和q列进行正交变换时，就对H施行：

在一次迭代过程当中需要对A的任意两列进行一次正交变换，迭代多次等效于施行

迭代的终止条件是为对角矩阵，即原矩阵H进过多次的双边Jacobi旋转后非对角线元素全部变成了0（实现计算当中不可能真的变为0，只要小于一个很小的数就可以认为是0了）。

每一个J都是标准正交阵，所以也是标准正交阵，记为V，则是

，则。从此式也可以看出对称矩阵的左奇异向量和右奇异向量是相等的。V也是H的特征向量。

当特征值是0时，对应的U i和V i不用求，我们只需要U和V的前r列就可以恢复矩阵A了（r是非0奇异值的个数），这也正是奇异值分解可以进行数据压缩的原因。

+ View Code

基于单边Jacobi旋转的SVD算法

相对于双边，单边的计算量小，并且容易并行实现。

单边Jacobi方法直接对原矩阵A进行单边正交旋转，A可以是任意矩阵。

同样每一轮的迭代都要使A的任意两列都正交一次，迭代退出的条件是B的任意两列都正交。

单边Jacobi旋转有这么一个性质：旋转前若，则旋转后依然是；反之亦然。

把矩阵B每列的模长提取出来，，把记为V，则。+ View Code

基于单边Jacobi旋转的SVD并行算法

单边Jacobi之于双边Jacobi的一个不同就是每次旋转只影响到矩阵A的两列，因经很多列对的正交旋转变换是可以并行执行的。

基本思想是将A按列分块，每一块分配到一个计算节点上，块内任意两列之间进行单边Jacobi正交变换，所有的计算节点可以同时进行。问题就是矩阵块内部列与列之间进行正交变换是不够的，我们需要所有矩阵块上的任意两列之间都进行正交变换，这就需要计算节点之间交换矩阵的列。本文采用奇偶序列的方法，具体可参考文

章/Article/CDMD-10285-1012286387.htm。

由于这次我用的是C版的MPI（MPI也有C++和Fortan版的），所以上面代码用到的C++版的matrix.h就不能用了，需要自己写一个C版的matrix.h。

matrix.h

svd.c