svd分解 c语言实现

/*

本程序在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=0i<ni++)

        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=0i<ni++)

        x[i]/=r

    return r

}

inline double product(double*a, double *b,int n){

    double r=0

    for(int i=0i<ni++)

        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=0i<ni++)

        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=0i<Ki++)

        U[i].resize(M,0)

    V.resize(K)

    for(int i=0i<Ki++)

        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=0col<Kcol++){

        double diff=1

        double r=-1

        while(1){

            for(int i=0i<Mi++)

                left_vector[i]= (float)rand() / RAND_MAX

            if(normalize(left_vector, M)>eps)

                break

        }

        for(int iter=0diff>=eps &&iter<MAX_ITERiter++){

            memset(next_left_vector,0,sizeof(double)*M)

            memset(next_right_vector,0,sizeof(double)*N)

            for(int i=0i<Mi++)

                for(int j=0j<Nj++)

                    next_right_vector[j]+=left_vector[i]*A[i][j]

            r=normalize(next_right_vector,N)

            if(r<eps) break

            for(int i=0i<coli++)

                orth(&V[i][0],next_right_vector,N)

            normalize(next_right_vector,N)

            for(int i=0i<Mi++)

                for(int j=0j<Nj++)

                    next_left_vector[i]+=next_right_vector[j]*A[i][j]

            r=normalize(next_left_vector,M)

            if(r<eps) break

            for(int i=0i<coli++)

                orth(&U[i][0],next_left_vector,M)

            normalize(next_left_vector,M)

            diff=0

            for(int i=0i<Mi++){

                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=10

    int n=8

    int k=5

    //分解一个10*8的矩阵A，求其前5个奇异值和奇异向量

    srand(time(0))

    vector<vector<double>>A

    A.resize(m)

    for(int i=0i<mi++){

        A[i].resize(n)

        for(int j=0j<nj++)

            A[i][j]=(float)rand()/RAND_MAX-0.5

    }

    cout<<"A="<<endl

    for(int i=0i<A.size()i++){

        for(int j=0j<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=0i<U[0].size()i++){

        for(int j=0j<U.size()j++){

            cout<<setw(12)<<U[j][i]<<' '

        }

        cout<<endl

    }

    cout<<endl

    cout<<"S="<<endl

    for(int i=0i<S.size()i++){

        cout<<setw(7)<<S[i]<<' '

    }

    cout<<endl

    cout<<"V="<<endl

    for(int i=0i<V[0].size()i++){

        for(int j=0j<V.size()j++){

            cout<<setw(12)<<V[j][i]<<' '

        }

        cout<<endl

    }

    return 0

}

/** C++ function for SVD

函数原型:

bool svd(vector<vector<double>>A, int K, std::vector<std::vector<double>>&U, std::vector<double>&S, std::vector<std::vector<double>>&V)

其中

A是输入矩阵，假设A的维数是m*n，那么本函数将A分解为U diag(S) V'

其中U是m*K的列正交的矩阵. V是n*K的列正交矩阵，S是K维向量。K由第二个参数指定。

U的第i列是A的第i大奇异值对应的左歧义向量,S[i]=A的第 i大奇异值,V的第i列是A的第i大奇异值对应的右歧义响亮.

K是需要分解的rank，0<K<=min(m,n)

本程序采用的是最基本幂迭代算法，在linux g++下编译通过

**/

#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=0i<ni++)

        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=0i<ni++)

        x[i]/=r

    return r

}

inline double product(double*a, double *b,int n){

    double r=0

    for(int i=0i<ni++)

        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=0i<ni++)

        b[i]-=r*a[i]

}

bool svd(vector<vector<double>>A, int K, std::vector<std::vector<double>>&U, std::vector<double>&S, std::vector<std::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=0i<Ki++)

        U[i].resize(M,0)

    V.resize(K)

    for(int i=0i<Ki++)

        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=0i<Mi++)

            left_vector[i]= (float)rand() / RAND_MAX

        if(normalize(left_vector, M)>eps)

            break

    }

    int col=0

    for(int col=0col<Kcol++){

        double diff=1

        double r=-1

        for(int iter=0diff>=eps &&iter<MAX_ITERiter++){

            memset(next_left_vector,0,sizeof(double)*M)

            memset(next_right_vector,0,sizeof(double)*N)

            for(int i=0i<Mi++)

                for(int j=0j<Nj++)

                    next_right_vector[j]+=left_vector[i]*A[i][j]

            r=normalize(next_right_vector,N)

            if(r<eps) break

            for(int i=0i<coli++)

                orth(&V[i][0],next_right_vector,N)

            normalize(next_right_vector,N)

            for(int i=0i<Mi++)

                for(int j=0j<Nj++)

                    next_left_vector[i]+=next_right_vector[j]*A[i][j]

            r=normalize(next_left_vector,M)

            if(r<eps) break

            for(int i=0i<coli++)

                orth(&U[i][0],next_left_vector,M)

            normalize(next_left_vector,M)

            diff=0

            for(int i=0i<Mi++){

                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<vector<double>>&A){

    for(int i=0i<A.size()i++){

        for(int j=0j<A[i].size()j++){

            cout<<setprecision(3)<<A[i][j]<<' '

        }

        cout<<endl

    }

}

int main(){

    int m=10

    int n=5

    srand(time(0))

    vector<vector<double>>A

    A.resize(m)

    for(int i=0i<mi++){

        A[i].resize(n)

        for(int j=0j<nj++)

            A[i][j]=(float)rand()/RAND_MAX

    }

    print(A)

    cout<<endl

    vector<vector<double>>U

    vector<double>S

    vector<vector<double>>V

    svd(A,2,U,S,V)

    cout<<"U="<<endl

    print(U)

    cout<<endl

    cout<<"S="<<endl

    for(int i=0i<S.size()i++){

        cout<<S[i]<<' '

    }

    cout<<endl

    cout<<"V="<<endl

    print(V)

    return 0

}

X= U*S*V'

这里的V是转置的

所以跟你的C程序结果是转置

其次正负号可以随便分配的，例如你把U的符号全取负号，那么V也是全取负号

然后就是精度问题，你在matlab 先输入format long

就能看到更多位数了，大小基本是一样的

所以其实结果是等价的，说白了就是S按特征值大小排序的话是唯一的，U,V不唯一

不明白可追问