并行计算奇异值分解--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的右奇异值，的特征值开方后就是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奇异值的个数），这也正是奇异值分解可以进行数据压缩的原因。

#include"matrix.h"
#include<cmath>
#include<map>
 
using namespace std;
 
const double THRESHOLD = 1E-8;
const int ITERATION = 30;   //迭代次数的上限
 
inline int sign(double number)
{
    if (number < 0)
        return -1;
    else
        return 1;
}
 
void rotate(Matrix < double >&matrix, int i, int j, bool & pass,
        Matrix < double >&J)
{
    double ele = matrix.get(i, j);
    if (fabs(ele) < THRESHOLD) {
        return;
    }
    pass = false;
    double ele1 = matrix.get(i, i);
    double ele2 = matrix.get(j, j);
    int size=matrix.getRows();
    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;
    Matrix < double >G(IdentityMatrix < double >(size, size));
    G.put(i, i, cos);
    G.put(i, j, -1 * sin);
    G.put(j, i, sin);
    G.put(j, j, cos);
    matrix = G.getTranspose() * matrix * G;
    J *= G;
}
 
void jacobi(Matrix < double >&matrix, int size, vector < double >&E,
        Matrix < double >&J)
{
 
    int iteration = ITERATION;
    while (iteration-- > 0) {
        bool pass = true;
        for (int i = 0; i < size; ++i) {
            for (int j = i+1; j < size; ++j) {
                rotate(matrix, i, j, pass, J);
            }
        }
        if (pass)   //当非对角元素全部变为0时迭代退出
            break;
    }
    cout << "迭代次数：" << ITERATION - iteration << endl;
    for (int i = 0; i < size; ++i) {
        E[i] = matrix.get(i, i);
        if (E[i] < THRESHOLD)
            E[i] = 0.0;
    }
}
 
void svd(Matrix < double >&A, Matrix < double >&U, Matrix < double >&V,
     vector < double >&E)
{
    int rows = A.