SVD分解原理及C++实现

SVD分解

1）特征值(EVD)变换

如果说一个向量v是方阵A的特征向量，将一定可以表示成下面的形式：

这时候λ 就被称为特征向量v对应的特征值，一个矩阵的一组特征向量是一组正交向量。特征值分解是将一个矩阵分解成下面的形式：若A 有n个特征值，则A矩阵

$P^{-1} AP = \sum$ 即$A P= P\sum$

由$A(\nu_1, \nu_2...\nu_i) =(\nu_1\lambda_1, \nu_2\lambda_2...\nu_i\lambda_i) = (\nu_1, \nu_2...\nu_i) \sum$ 得到$AW =W\sum$

即特征矩阵 $W$即为相似对角形的矩阵$P$ ，一个矩阵其实就是一个线性变换，因为一个矩阵乘以一个向量后得到的向量，其实就相当于将这个向量进行了线性变换。

当矩阵$A$ 是高维的情况下，那么这个矩阵就是高维空间下的一个线性变换，这个线性变化可能没法通过图片来表示，但是可以想象，这个变换也同样有很多的变换方向，我们通过特征值分解得到的前N个特征向量，那么就对应了这个矩阵最主要的N个变化方向。我们利用这前N个变化方向，就可以近似这个矩阵（变换）。也就是之前说的：提取这个矩阵最重要的特征。总结一下，特征值分解可以得到特征值与特征向量，特征值表示的是这个特征到底有多重要，而特征向量表示这个特征是什么 ，可以将每一个特征向量理解为一个线性的子空间，我们可以利用这些线性的子空间干很多的事情。不过，特征值分解也有很多的局限，比如说变换的矩阵必须是方阵。

2）SVD分解

奇异值分解是一个能适用于任意的矩阵的一种分解的方法，一般表示如下：
$A = U \sum V^T$

$A_{(n×m)} = U_{n×n} \sum_{n×m} V^T_{m×m}$ , $U$ 和$V$均为单位正交矩阵，即矩阵的行列向量线性无关，

上面的特征值分解的$A$矩阵是对称阵，根据EVD可以找到一个（超）矩形使得变换后还是（超）矩形，也即A可以将一组正交基映射到另一组正交基。SVD分解的精髓就是，对任意M*N的矩阵，找到一组正交基使得经过它变换后还是正交基。

现在假设存在M*N矩阵A，事实上，A矩阵将n维空间中的向量映射到k（k<=m）维空间中，k=Rank(A)。现在的目标就是：在n维空间中找一组正交基，使得经过A变换后还是正交的。假设已经找到这样一组正交基：

则A矩阵将这组基映射为：

如果要使他们两两正交，即

根据假设，存在

所以如果正交基v选择为A’A的特征向量的话，由于A’A是对称阵，v之间两两正交，那么

这样就找到了正交基使其映射后还是正交基了，现在，将映射后的正交基单位化：

因为

所以有

所以取单位向量

由此可得

当k < i <= m时，对u1，u2，…，uk进行扩展u(k+1),…,um，使得u1，u2，…，um为m维空间中的一组正交基，即

同样的，对v1，v2，…，vk进行扩展v(k+1),…,vn（这n-k个向量存在于A的零空间中，即Ax=0的解空间的基），使得v1，v2，…，vn为n维空间中的一组正交基，即

则可得到

继而可以得到A矩阵的奇异值分解：
$A = U \sum V^T$

$A_{(n×m)} = U_{n×n} \sum_{n×m} V^T_{m×m}$ , $U$ 和$V$均为单位正交矩阵，$\sum$ 为对角阵

现在可以来对A矩阵的映射过程进行分析了：如果在n维空间中找到一个（超）矩形，其边都落在A’A的特征向量的方向上，那么经过A变换后的形状仍然为（超）矩形！

vi为A’A的特征向量，称为A的右奇异向量，$u_i=Av_i$实际上为AA’的特征向量，称为A的左奇异向量。下面利用SVD证明满秩分解：

利用矩阵分块乘法展开得：

可以看到第二项为0，有

令

则$A=XY = U \sum V^T$，整个SVD的推导过程如上所述。

$A_{(n×m)} = U_{n×n} \sum_{n×m} V^T_{m×m}$，对矩阵$\sum$ 来说，奇异值$\sigma$减少的特征快，在很多情况下，前1%或者10%的奇异值之和就占了整个矩阵奇异值和的99%，所以我们可以用前K大的奇异值来近似描述矩阵，称为部分奇异值：

$A_{(n×m)} = U_{n×r} \sum_{r×r} V^T_{r×m}$ ，因此，矩阵A 的巨大数据就可以存储在三个小矩阵中，实现了数据的压缩和降维，奇异值的计算是一个$O(N^3)$的算法，计算复杂度高。

3）SVD与PCA

PCA的全部工作简单点说，就是对原始的空间中顺序地找一组相互正交的坐标轴，第一个轴是使得方差最大的，第二个轴是在与第一个轴正交的平面中使得方差最大的，第三个轴是在与第1、2个轴正交的平面中方差最大的，这样假设在N维空间中，我们可以找到N个这样的坐标轴，我们取前r个去近似这个空间，这样就从一个N维的空间压缩到r维的空间了，但是我们选择的r个坐标轴能够使得空间的压缩使得数据的损失最小。

还是假设我们矩阵每一行表示一个样本，每一列表示一个feature，用矩阵的语言来表示，将一个n * m的矩阵A的进行坐标轴的变化，P就是一个变换的矩阵从一个N维的空间变换到另一个N维的空间，在空间中就会进行一些类似于旋转、拉伸的变化。

$A_{n×m} P_{m×m} = \hat A_{n×m}$

（1）若要将矩阵$A_{n×m}$实现降维，设n为样本维度，m为样本数，要实现数据从n维降为k 维，则根据SVD分解，和PCA转换的矩阵表示：

$P_{k×n}A_{n×m} = \hat A_{k×m}$ , $\hat A_{k×m} = P_{k×n} A_{n×m}$

P 即为PCA原理中的k个特征向量$Wi$ (n×1) 组成的特征矩阵$W_{(n×k)}$d的转置矩阵$W_{(k×n)}^T$ ，

SVD分解$A_{(n×m)} = U_{n×k} \sum_{k×k} V^T_{k×m}$ ， 等式两边同时左乘一个$U_{k×n}^T$ 矩阵，

$U_{k×n}^TA_{(n×m)} = \sum_{k×k} V^T_{k×m} = \hat A_{k×m}$，即转换矩阵$P_{k×n} =W_{(k×n)}^T =U_{k×n}^T$ ，

（2）若要将矩阵$A_{n×m}$中的一些相似样本进行数据压缩，降维，实现原始数据从m维降为r 维，与上述转换类似，

$A_{n×m} P_{m×r} = \hat A_{n×r}$

由SVD分解 $A_{(n×m)} = U_{n×r} \sum_{r×r} V^T_{r×m}$ ，等式两边同时右乘一个$V_{m×r}$ 矩阵，

$A_{(n×m)} V_{m×r} = U_{n×r} \sum_{r×r} =\hat A_{n×r}$ ，即转换矩阵$P_{m×r} = V_{m×r}$ ，

由上述转化可以看出，其实PCA几乎可以说是对SVD的一个包装，如果我们实现了SVD，那也就实现了PCA了，而且更好的地方是，有了SVD，我们就可以得到两个方向的PCA，如果我们对$A^ T A$进行特征值的分解，只能得到一个方向的PCA。

4）SVD分解矩阵的代码实现

const int MAX_ITER=100000;      //迭代器迭代次数
const double eps=0.0000001;      //允许的最大误差

double get_norm(double *x, int n)   //求sqrt(x[1]*x[1] + x[2]*x[2] +......+ x[n]*x[n])，                                     //求向量x的模长|x|；
{
    double r=0;
    for(int i=0;i<n;i++)
        r+=x[i]*x[i];
    return sqrt(r);
}


double normalize(double *x, int n)   //当求的的向量x的模小于eps最小值时，舍去；
{                        //当x的模大于eps时，保留；并将x[i]/|x|，将向量x归一化为单位向量e；
    double r=get_norm(x,n);        
    if(r<eps)
    return 0;
    for(int i=0;i<n;i++)
        x[i]/=r;
    return r;                 //返回向量x的模长，模长不能太小，
}


inline double product(double*a, double *b,int n)  //求的向量a和向量b的内积，即(a[1]*b[1] + a[2]*b[2] +......+ a[n]*b[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 ，向量a为单位向量e，将向量b与单位
{                                         // 向量进行正交，
    double r=product(a,b,n);        //向量b与单位向量e的内积，即为向量b在单位向量e上的投影值
    for(int i=0;i<n;i++)
        b[i]-=r*a[i];                   //b'[i] = b[i] - (e, b)*b[i] ,得到与向量a正交的向                                         //量b' 即（a, b')的内积为0；
}

SVD函数实现如下：使用幂迭代法对矩阵进行SVD分解：

//函数将A分解为U diag(S) V'；S[i],U[i],V[i]是A的第i大奇异值，及其对应的左歧义向量和右奇异向量
//S,U,V的size由K指定,K是需要分解的rank，0<K<=min(m,n)

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(M×N) ，A = U S (V^T)； U(M×K)  S(K×K)  V(N×K)
    int M=A.size();
    int N=A[0].size();
    U.clear();
    V.clear();
    S.clear();

    S.resize(K,0);      //S为对角矩阵，只需保存对角线上的奇异值即可，故S取向量保存，节约空间；                         //S保存K个元素即可   

    U.resize(K);              // U矩阵为K个列向量；每个列向量有M个元素；
    for(int i=0;i<K;i++)
    {
        U[i].resize(M,0);
    }

    V.resize(K);            // V矩阵也为K个列向量；每个列向量有N个元素；
    for(int i=0;i<K;i++)
    {
        V[i].resize(N,0);
    }

    srand(time(0));            //给一个随机种子数,

    double *left_vector=new double[M];          //动态分配内存,产生left_vector和                                                           //next_left_vector向量用于迭代运算
    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<M;i++)
            left_vector[i]= (float)rand() / RAND_MAX;    //随机生成一个向量vector（即α），                                                       //共有M个元素，用来求解左奇异矩阵 
        if(normalize(left_vector, M)>eps)               //当向量的模长大于eps时，生成结束；                                                        //生成一个模长合适的迭代向量
        break;
    }

    int col=0;
    for(int col=0;col<K;col++)    //按列计算左右奇异矩阵,迭代
    {
        double diff=1;
        double r=-1;

        for(int iter=0; diff>=eps && iter<MAX_ITER; iter++)     //迭代器iter的迭代次数为10000次
        {
            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];   //向量α(1×M)×A(M×N)矩阵，得到右迭代向量β_next(1×N)

            r=normalize(next_right_vector,N);                 //单位化向量β_next

            if(r<eps) break;                                  //若β_next 模长太小，则退出内层循环
            for(int i=0;i<col;i++)
                orth(&V[i][0],next_right_vector,N);              //与右矩阵V正交，正交化得到β'_next
            normalize(next_right_vector,N);                      //单位化β'_next



            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]; //矩阵A(M×N)×(β'_next)^T(N×1)向量，得到左迭代向量α_next(1×M)

            r=normalize(next_left_vector,M);                       //单位化向量α_next

            if(r<eps) break;                                       //若α_next 模长太小，则退出内层循环
            for(int i=0;i<col;i++)                           
                orth(&U[i][0],next_left_vector,M);                 //与左矩阵U正交，正交化得到α'_next
            normalize(next_left_vector,M);                         //单位化α'_next

            diff=0;

            for(int i=0;i<M;i++)
            {
                double d=next_left_vector[i]-left_vector[i];      //计算前后两个迭代向量的2范数的平方，即两个向量差距的平方
                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;                                                  //若向量的模长大于最小值eps，则取为矩阵A的奇异值

            memcpy((char *)&U[col][0],left_vector,sizeof(double)*M);   //按列拷贝左迭代向量到左奇异矩阵U中

            memcpy((char *)&V[col][0],right_vector,sizeof(double)*N);  //按列拷贝右迭代向量到右奇异矩阵V中

        }
        else break;
    }

    delete [] next_left_vector;
    delete [] next_right_vector;
    delete [] left_vector;
    delete [] right_vector;

    return true;
}