奇异值分解(SVD) C++代码

/**
函数原型:
bool svd(vector > A, int K, std::vector > &U, std::vector &S, std::vector > &V);
输入矩阵A,分解矩阵的秩K
输出U,S,V
本函数将A分解为U diag(S) V'
S[i],U[i],V[i]是A的第i大奇异值，及其对应的左歧义向量和右奇异向量
S,U,V的size由K指定
K是需要分解的rank，0本程序采用的是最基本幂迭代算法，在linux g++下编译通过
**/
#include
#include
#include
#include
#include
#include
#include
using namespace std;
const int MAX_ITER=100000;
const double eps=0.0000001;
double get_norm(double *x, int n){
double r=0;
for(int i=0;i r+=x[i]*x[i];
return sqrt(r);
}
double normalize(double *x, int n){
double r=get_norm(x,n);
if(r return 0;
for(int i=0;i x[i]/=r;
return r;
}
inline double product(double*a, double *b,int n){
double r=0;
for(int i=0;i r+=a[i]*b[i];
return r;
}
void orth(double *a, double *b, int n){//|a|=1
double r=product(a,b,n);
for(int i=0;i b[i]-=r*a[i];

}
bool svd(vector > A, int K, std::vector > &U, std::vector &S, std::vector > &V){
int M=A.size();
int N=A[0].size();
U.clear();
V.clear();
S.clear();
S.resize(K,0);
U.resize(K);
for(int i=0;i U[i].resize(M,0);
V.resize(K);
for(int i=0;i V[i].resize(N,0);

srand(time(0));
double *left_vector=new double[M];
double *next_left_vector=new double[M];
double *right_vector=new double[N];
double *next_right_vector=new double[N];
while(1){
for(int i=0;i left_vector[i]= (float)rand() / RAND_MAX;
if(normalize(left_vector, M)>eps)
break;
}
int col=0;
for(int col=0;col double diff=1;
double r=-1;
for(int iter=0;diff>=eps && iter memset(next_left_vector,0,sizeof(double)*M);
memset(next_right_vector,0,sizeof(double)*N);
for(int i=0;i for(int j=0;j next_right_vector[j]+=left_vector[i]*A[i][j];
r=normalize(next_right_vector,N);
if(r for(int i=0;i orth(&V[i][0],next_right_vector,N);
normalize(next_right_vector,N);
for(int i=0;i for(int j=0;j next_left_vector[i]+=next_right_vector[j]*A[i][j];
r=normalize(next_left_vector,M);
if(r for(int i=0;i orth(&U[i][0],next_left_vector,M);
normalize(next_left_vector,M);
diff=0;
for(int i=0;i double d

=next_left_vector[i]-left_vector[i];
diff+=d*d;
}
memcpy(left_vector,next_left_vector,sizeof(double)*M);
memcpy(right_vector,next_right_vector,sizeof(double)*N);
}
if(r>=eps){
S[col]=r;
memcpy((char *)&U[col][0],left_vector,sizeof(double)*M);
memcpy((char *)&V[col][0],right_vector,sizeof(double)*N);
}else
break;
}
delete [] next_left_vector;
delete [] next_right_vector;
delete [] left_vector;
delete [] right_vector;
return true;
}
void print(vector > &A){
for(int i=0;i for(int j=0;j cout< }
cout< }
}
int main(){
int m=10;
int n=5;
srand(time(0));
vector > A;
A.resize(m);

for(int i=0;i A[i].resize(n);
for(int j=0;j A[i][j]=(float)rand()/RAND_MAX;
}
print(A);
cout< vector > U;
vector S;
vector > V;
svd(A,2,U,S,V);
cout<<"U="< print(U);
cout< cout<<"S="< for(int i=0;i cout< }
cout< cout<<"V="< print(V);
return 0;
}