并行计算奇异值分解--Jacobi旋转
鉴于矩阵的奇异值分解SVD在工程领域的广泛应用（如数据压缩、噪声去除、数值分析等等，包括在NLP领域的潜在语义索引LSI核心操作也是SVD），今天就详细介绍一种SVD的实现方法--Jacobi旋转法。跟其它SVD算法相比，Jacobi法精度高，虽然速度慢，但容易并行实现。
一些链接
http://cdmd.cnki.com.cn/Article/CDMD-10285-1012286387.htm 并行JACOBI方法求解矩阵奇异值的研究。本文呈现的代码就是依据这篇论文写出来的。
http://math.nist.gov/javanumerics/jama/ Jama包是用于基本线性代数运算的java包，提供矩阵的cholesky分解、LUD分解、QR分解、奇异值分解，以及PCA中要用到的特征值分解，此外可以计算矩阵的乘除法、矩阵的范数和条件数、解线性方程组等。
http://users.telenet.be/paul.larmuseau/SVD.htm 在线SVD运算器。
http://www.bluebit.gr/matrix-calculator/ bluebit在线矩阵运算器，提供矩阵的各种运算。
http://www.drque.net/Projects/Matrix/ C++ Matrix library提供矩阵的加减乘除、求行列式、LU分解、求逆、求转置。本文的头两段程序就引用了这里面的matrix.h。
基于双边Jacobi旋转的奇异值分解算法
V是A的右奇异向量，也是的特征向量;
U是A的左奇异向量，也是的特征向量。
特别地，当A是对称矩阵的时候，=，即U=V，U的列向量不仅是的特征向量，也是A
的特征向量。这一点在主成分分析中会用到。
对于正定的对称矩阵，奇异值等于特征值，奇异向量等于特征向量。 U、V都是正交矩阵，满足矩阵的转置即为矩阵的逆。
双边Jacobi方法本来是用来求解对称矩阵的特征值和特征向量的，由于就是对称矩阵，求出的特
征向量就求出了A的右奇异值，一个Jacobi旋转矩阵J形如：
的特征值开方后就是A的奇异值。
它就是在一个单位矩阵的基础上改变p行q行p列q列四个交点上的值，J矩阵是一个标准正交矩阵。
当我们要对H和p列和q列进行正交变换时，就对H施行：
在一次迭代过程当中需要对A的任意两列进行一次正交变换，迭代多次等效于施行
迭代的终止条件是为对角矩阵，即原矩阵H进过多次的双边
Jacobi旋转后非对角线元素全部变成了0（实现计算当中不可能真的变为0，只要小于一个很小的数就可以认为是0了）。
每一个J都是标准正交阵，所以也是标准正交阵，记为V，则是
，则
征向量。
。从此式也可以看出对称矩阵的左奇异向量和右奇异向量是相等的。V也是H的特
当特征值是0时，对应的Ui和Vi不用求，我们只需要U和V的前r列就可以恢复矩阵A了（r是非0奇异值的
个数），这也正是奇异值分解可以进行数据压缩的原因。 + View Code
基于单边Jacobi旋转的SVD算法
相对于双边，单边的计算量小，并且容易并行实现。
单边Jacobi方法直接对原矩阵A进行单边正交旋转，A可以是任意矩阵。
同样每一轮的迭代都要使A的任意两列都正交一次，迭代退出的条件是B的任意两列都正交。
单边Jacobi旋转有这么一个性质：旋转前若反之亦然。
，则旋转后依然是；
把矩阵B每列的模长提取出来，+ View Code
，把记为V，则。
基于单边Jacobi旋转的SVD并行算法
单边Jacobi之于双边Jacobi的一个不同就是每次旋转只影响到矩阵A的两列，因经很多列对的正交旋转变换是可以并行执行的。
基本思想是将A按列分块，每一块分配到一个计算节点上，块内任意两列之间进行单边Jacobi正交变换，所有的计算节点可以同时进行。问题就是矩阵块内部列与列之间进行正交变换是不够的，我们需要所有矩阵块上的任意两列之间都进行正交变换，这就需要计算节点之间交换矩阵的列。本文采用奇偶序列的方法，具体可参考文章 http://cdmd.cnki.com.cn/Article/CDMD-10285-1012286387.htm。
由于这次我用的是C版的MPI（MPI也有C++和Fortan版的），所以上面代码用到的C++版的matrix.h就不能用了，需要自己写一个C版的matrix.h。 matrix.h
#ifndef _MATRIX_H #define _MATRIX_H #include