SVD分解的并行实现

在之前的文章中，我对SVD进行了大致的了解，它在互联网高端领域中有广泛的应用，至于它的一些详细应

用以后再进一步学习。现在主要的问题是如何有效地实现SVD分解，接下来我会先用两种方法来实现SVD分

解，即基于双边Jacobi旋转的SVD和基于单边Jacobi旋转的SVD。

这两种方法重点参考中国知网上的一篇论文：《并行JACOBI方法求解矩阵奇异值的研究》

代码：

#include <string.h>
#include <stdio.h>
#include <math.h>
#include <map>
 
#include "matrix.h"
 
using namespace std;
 
const double EPS = 1e-8;
const int ITER_NUM = 30;
 
int sign(double x)
{
    if(x < 0) return -1;
    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) < EPS) return;
 
    pass = false;
    double ele1 = matrix.get(i, i);
    double ele2 = matrix.get(j, j);
    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 size = matrix.getRows();
    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 iter_num = ITER_NUM;
    while(iter_num-- > 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) break;
    }
    cout << "迭代次数：" << ITER_NUM - iter_num << endl;
 
    for(int i = 0; i < size; i++)
    {
        E[i] = matrix.get(i, i);
        if(E[i] < EPS)
            E[i] = 0.0;
    }
}
 
void svd(Matrix<double> &A, Matrix<double> &U, Matrix<double> &V, vector<double> &E)
{
    int rows = A.getRows();
    int columns = A.getColumns();
    assert(rows <= columns);
    assert(U.getRows() == rows);
    assert(U.getColumns() == rows);
    assert(V.getRows() == columns);
    assert(V.getColumns() == columns);
    assert(E.size() == columns);
 
    Matrix<double> B = A.getTranspose() * A;
    Matrix<double> J(IdentityMatrix<double>(columns, columns));
    vector<double> S(columns);
    jacobi(B, columns, S, J);
    for(int i = 0; i < S.size(); i++)
        S[i] = sqrt(S[i]);
 
    multimap<double, int> eigen;
    for(int i = 0; i < S.size(); i++)
        eigen.insert(make_pair(S[i], i));
    multimap<double, int>::const_iterator iter = --eigen.end();
    int num_eig = 0;
    for(int i = 0; i < columns; i++, iter--)
    {
        int index = iter->second;
        E[i] = S[index];
        if(E[i] > EPS)
            num_eig++;
        for(int row = 0; row < columns; row++)
            V.put(row, i, J.get(row, index));
    }
 
    assert(num_eig <= rows);
    for(int i = 0; i < num_eig; i++)
    {
        Matrix<double> vi = V.getColumn(i);
        double sigma = E[i];
        Matrix<double> ui(rows, 1);
        ui = A * vi;
        for(int j = 0; j < rows; j++)
            U.put(j, i, ui.get(j, 0) / sigma);
    }
}
 
int main()
{
    const int ROW = 4;
    const int COL = 5;
    assert(ROW <= COL);
    double arr[ROW * COL] =
    {1, 0, 0, 0, 2,
     0, 0, 3, 0, 0,
     0, 0, 0, 0, 0,
     0, 4, 0, 0, 0
     };
    Matrix<double> A(ROW, COL);
    A = arr;
 
    Matrix<double> U(ROW, ROW);
    Matrix<double> V(COL, COL);
    vector<double> E(COL);
 
    svd(A, U, V, E);
    Matrix<double> Sigma(ROW, COL);
    for(int i = 0; i < ROW; i++)
        Sigma.put(i, i, E[i]);
 
    cout << "U = " << endl << U;
    cout << "Sigma = " << endl << Sigma;
    cout << "V^T = " << endl << V.getTranspose();
 
    Matrix<double> AA = U * Sigma * V.getTranspose();
    cout << "Result AA = " << endl << AA;
 
    return 0;
}