奇异值分解 SVD C++代码

来源：互联网发布：域名不以封疆编辑：程序博客网时间：2024/06/04 22:02
/*本程序在linux g++下编译通过bool svd(vector<vector<double> > A, int K, vector<vector<double> > &U, vector<double> &S, vector<vector<double> > &V);A: 输入待分解矩阵K: 输入，取前K大奇异值及奇异向量U[0],U[1],...,U[K-1]: 前K大奇异值对应的左奇异向量S[0],S[1],...,S[K-1]: 前K大奇异值 S[0]>=S[1]>=...>=S[K-1]V[0],V[1],...,V[K-1]: 前K大奇异值对应的右奇异向量*/#include <cmath>#include <iostream>#include <iomanip>#include <cstdlib>#include <cstring>#include <fstream>#include <vector>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<n;i++)        r+=x[i]*x[i];    return sqrt(r);}double normalize(double *x, int n){    double r=get_norm(x,n);    if(r<eps)        return 0;    for(int i=0;i<n;i++)        x[i]/=r;    return r;}inline double product(double*a, double *b,int n){    double r=0;    for(int i=0;i<n;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<n;i++)        b[i]-=r*a[i];    }bool svd(vector<vector<double> > A, int K, vector<vector<double> > &U, vector<double> &S, vector<vector<double> > &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<K;i++)        U[i].resize(M,0);    V.resize(K);    for(int i=0;i<K;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];    int col=0;    for(int col=0;col<K;col++){        double diff=1;        double r=-1;        while(1){            for(int i=0;i<M;i++)                left_vector[i]= (float)rand() / RAND_MAX;            if(normalize(left_vector, M)>eps)                break;        }        for(int iter=0;diff>=eps && iter<MAX_ITER;iter++){            memset(next_left_vector,0,sizeof(double)*M);            memset(next_right_vector,0,sizeof(double)*N);            for(int i=0;i<M;i++)                for(int j=0;j<N;j++)                    next_right_vector[j]+=left_vector[i]*A[i][j];            r=normalize(next_right_vector,N);            if(r<eps) break;            for(int i=0;i<col;i++)                orth(&V[i][0],next_right_vector,N);            normalize(next_right_vector,N);            for(int i=0;i<M;i++)                for(int j=0;j<N;j++)                    next_left_vector[i]+=next_right_vector[j]*A[i][j];            r=normalize(next_left_vector,M);            if(r<eps) break;            for(int i=0;i<col;i++)                orth(&U[i][0],next_left_vector,M);            normalize(next_left_vector,M);            diff=0;            for(int i=0;i<M;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{            cout<<r<<endl;            break;        }    }    delete [] next_left_vector;    delete [] next_right_vector;    delete [] left_vector;    delete [] right_vector;    return true;}void print(vector<vector<double> > &A){}int main(){    int m=1000;    int n=800;    int k=5;    //分解一个10*8的矩阵A，求其前5个奇异值和奇异向量    srand(time(0));    vector<vector<double> > A;    A.resize(m);        for(int i=0;i<m;i++){        A[i].resize(n);        for(int j=0;j<n;j++)            A[i][j]=(float)rand()/RAND_MAX-0.5;    }        cout<<"A="<<endl;    for(int i=0;i<A.size();i++){        for(int j=0;j<A[i].size();j++){            cout<<setw(12)<<A[i][j]<<' ';        }        cout<<endl;    }    cout<<endl;    vector<vector<double> > U;    vector<double> S;    vector<vector<double> > V;    svd(A,k,U,S,V);    cout<<"U="<<endl;    for(int i=0;i<U[0].size();i++){        for(int j=0;j<U.size();j++){            cout<<setw(12)<<U[j][i]<<' ';        }        cout<<endl;    }    cout<<endl;    cout<<"S="<<endl;    for(int i=0;i<S.size();i++){        cout<<setw(7)<<S[i]<<' ';    }    cout<<endl;    cout<<"V="<<endl;    for(int i=0;i<V[0].size();i++){        for(int j=0;j<V.size();j++){            cout<<setw(12)<<V[j][i]<<' ';        }        cout<<endl;    }    return 0;}
0 0