生活随笔
收集整理的這篇文章主要介紹了
SVD分解的并行实现
小編覺得挺不錯的,現在分享給大家,幫大家做個參考.
在之前的文章中,我對SVD進行了大致的了解,它在互聯網高端領域中有廣泛的應用,至于它的一些詳細應
用以后再進一步學習。現在主要的問題是如何有效地實現SVD分解,接下來我會先用兩種方法來實現SVD分
解,即基于雙邊Jacobi旋轉的SVD和基于單邊Jacobi旋轉的SVD。
?
這兩種方法重點參考中國知網上的一篇論文:《并行JACOBI方法求解矩陣奇異值的研究》
?
代碼:
[cpp]?view plaincopy
#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;?? }??
?
供參考文章:http://www.cnblogs.com/zhangchaoyang/articles/2575948.html
總結
以上是生活随笔為你收集整理的SVD分解的并行实现的全部內容,希望文章能夠幫你解決所遇到的問題。
如果覺得生活随笔網站內容還不錯,歡迎將生活随笔推薦給好友。