SVD分解的并行实现

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

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

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

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

代码：

[cpp]  view plain  copy  

1. #include <string.h>  
2. #include <stdio.h>  
3. #include <math.h>  
4. #include <map>  
5.    
6. #include "matrix.h"  
7.    
8. using namespace std;  
9.    
10. const double EPS = 1e-8;  
11. const int ITER_NUM = 30;  
12.    
13. int sign(double x)  
14. {  
15.     if(x < 0) return -1;  
16.     return 1;  
17. }  
18.    
19. void rotate(Matrix<double> &matrix, int i, int j, bool &pass, Matrix<double> &J)  
20. {  
21.     double ele = matrix.get(i, j);  
22.     if(fabs(ele) < EPS) return;  
23.    
24.     pass = false;  
25.     double ele1 = matrix.get(i, i);  
26.     double ele2 = matrix.get(j, j);  
27.     double tao = (ele1 - ele2) / (2 * ele);  
28.     double tan = sign(tao) / (fabs(tao) + sqrt(1 + pow(tao, 2)));  
29.     double cos = 1 / sqrt(1 + pow(tan, 2));  
30.     double sin = cos * tan;  
31.    
32.     int size = matrix.getRows();  
33.     Matrix<double>G(IdentityMatrix<double>(size, size));  
34.     G.put(i, i, cos);  
35.     G.put(i, j, -1 * sin);  
36.     G.put(j, i, sin);  
37.     G.put(j, j, cos);  
38.    
39.     matrix = G.getTranspose() * matrix * G;  
40.     J *= G;  
41. }  
42.    
43. void jacobi(Matrix<double> &matrix, int size, vector<double> &E, Matrix<double> &J)  
44. {  
45.     int iter_num = ITER_NUM;  
46.     while(iter_num-- > 0)  
47.     {  
48.         bool pass = true;  
49.         for(int i = 0; i < size; i++)  
50.         {  
51.             for(int j = i + 1; j < size; j++)  
52.                 rotate(matrix, i, j, pass, J);  
53.         }  
54.         if(pass) break;  
55.     }  
56.     cout << "迭代次数：" << ITER_NUM - iter_num << endl;  
57.    
58.     for(int i = 0; i < size; i++)  
59.     {  
60.         E[i] = matrix.get(i, i);  
61.         if(E[i] < EPS)  
62.             E[i] = 0.0;  
63.     }  
64. }  
65.    
66. void svd(Matrix<double> &A, Matrix<double> &U, Matrix<double> &V, vector<double> &E)  
67. {  
68.     int rows = A.getRows();  
69.     int columns = A.getColumns();  
70.     assert(rows <= columns);  
71.     assert(U.getRows() == rows);  
72.     assert(U.getColumns() == rows);  
73.     assert(V.getRows() == columns);  
74.     assert(V.getColumns() == columns);  
75.     assert(E.size() == columns);  
76.    
77.     Matrix<double> B = A.getTranspose() * A;  
78.     Matrix<double> J(IdentityMatrix<double>(columns, columns));  
79.     vector<double> S(columns);  
80.     jacobi(B, columns, S, J);  
81.     for(int i = 0; i < S.size(); i++)  
82.         S[i] = sqrt(S[i]);  
83.    
84.     multimap<double, int> eigen;  
85.     for(int i = 0; i < S.size(); i++)  
86.         eigen.insert(make_pair(S[i], i));  
87.     multimap<double, int>::const_iterator iter = --eigen.end();  
88.     int num_eig = 0;  
89.     for(int i = 0; i < columns; i++, iter--)  
90.     {  
91.         int index = iter->second;  
92.         E[i] = S[index];  
93.         if(E[i] > EPS)  
94.             num_eig++;  
95.         for(int row = 0; row < columns; row++)  
96.             V.put(row, i, J.get(row, index));  
97.     }  
98.    
99.     assert(num_eig <= rows);  
100.     for(int i = 0; i < num_eig; i++)  
101.     {  
102.         Matrix<double> vi = V.getColumn(i);  
103.         double sigma = E[i];  
104.         Matrix<double> ui(rows, 1);  
105.         ui = A * vi;  
106.         for(int j = 0; j < rows; j++)  
107.             U.put(j, i, ui.get(j, 0) / sigma);  
108.     }  
109. }  
110.    
111. int main()  
112. {  
113.     const int ROW = 4;  
114.     const int COL = 5;  
115.     assert(ROW <= COL);  
116.     double arr[ROW * COL] =  
117.     {1, 0, 0, 0, 2,  
118.      0, 0, 3, 0, 0,  
119.      0, 0, 0, 0, 0,  
120.      0, 4, 0, 0, 0  
121.      };  
122.     Matrix<double> A(ROW, COL);  
123.     A = arr;  
124.    
125.     Matrix<double> U(ROW, ROW);  
126.     Matrix<double> V(COL, COL);  
127.     vector<double> E(COL);  
128.    
129.     svd(A, U, V, E);  
130.     Matrix<double> Sigma(ROW, COL);  
131.     for(int i = 0; i < ROW; i++)  
132.         Sigma.put(i, i, E[i]);  
133.    
134.     cout << "U = " << endl << U;  
135.     cout << "Sigma = " << endl << Sigma;  
136.     cout << "V^T = " << endl << V.getTranspose();  
137.    
138.     Matrix<double> AA = U * Sigma * V.getTranspose();  
139.     cout << "Result AA = " << endl << AA;  
140.    
141.     return 0;  
142. }  

供参考文章：http://www.cnblogs.com/zhangchaoyang/articles/2575948.html