SVD 奇异值分解纯手工实现（C++）

6月10日，风雨大作。。。（好吧雨早停了）
一时兴起，想着好久没敲过CPP了，最近正好打算复习一下《统计学习方法》，本来想看看PCA来着，一打眼瞅到了SVD。。。SVD是一个比较基础的矩阵分解法，后边的PCA还会用到，所以干脆计划了晚上小小复习一下SVD，然后用久违的CPP实现一下，谁料想本以为挺早就能下班的，一直写到第二天下午。。。。太菜了！！！！但是实现完后还是挺开心的，因此记录一下这一天的成果。

关于SVD

个人以为SVD好用就是好用在条件弱，效果好，任意的实矩阵都能进行SVD，理论依据如下：

奇异值分解定理

若 A A A为一 m × n text{m} times{text{n}} m×n实矩阵， A ∈ R m × n Ain R^{text{m} times{text{n}}} A∈Rm×n，则 A A A的奇异值分解存在，且:
A = U Σ V T A=USigma V^T A=UΣVT
其中：
Σ = [ Σ 1 , Σ 2 ] = ( σ 1 σ 2 ⋱ σ r ⋱ σ n ) Sigma=[Sigma_1,Sigma_2] = begin{pmatrix} sigma_1\ & sigma_2 \ && ddots \ &&& sigma_r \ &&&& ddots\ &&&&&sigma_n\ end{pmatrix} Σ=[Σ1,Σ2]=⎝⎜⎜⎜⎜⎜⎜⎛σ1σ2⋱σr⋱σn⎠⎟⎟⎟⎟⎟⎟⎞
σ i sigma_i σi为 A T A A^TA ATA的特征值开根号，称为矩阵 A A A的奇异值。 r r r为非零特征值的个数， σ i ≥ σ i + 1 , i = 0 , 1 , . . . , n − 1 sigma_i ge sigma_{i+1},i=0,1,...,n-1 σi≥σi+1,i=0,1,...,n−1. σ j = 0 , j > r sigma_j=0,j>r σj=0,j>r 。
V = [ V 1 , V 2 ] = [ v 1 , v 2 , . . . , v r , . . . v n ] V=[V_1,V_2]=[v_1,v_2,...,v_r,...v_n] V=[V1,V2]=[v1,v2,...,vr,...vn]
V 1 V_1 V1为非零特征值对应的特征向量按特征值从大到小排列的列向量组， V 2 V_2 V2为特征值为零的特征向量组。
U = A V Σ − 1 = A [ V 1 , V 2 ] Σ U=AVSigma^{-1}=A[V_1,V_2]Sigma U=AVΣ−1=A[V1,V2]Σ
证明也为构造性证明，具体见《统计学习方法》。

紧奇异值分解

其实，通过一些矩阵运算，可以得到：
A = U Σ V T = U 1 Σ 1 V 1 T A=USigma V^T=U_1Sigma_1V_1^T A=UΣVT=U1Σ1V1T
U 1 ∈ R m × r U_1 in{R^{rm mtimes r}}, U1∈Rm×r, Σ 1 ∈ R r × r Sigma_1 in{R^{rm rtimes r}}, Σ1∈Rr×r, V 1 ∈ R n × r V_1in{rm R^{ntimes r}} V1∈Rn×r.
可以看到，原来矩阵 A ∈ R m × n Ain{rm R^{mtimes n}} A∈Rm×n,有 m ⋅ n mcdot n m⋅n个元素，而经过紧奇异值分解后，有 m ⋅ r + r ⋅ r + n ⋅ r mcdot r+rcdot r+ncdot r m⋅r+r⋅r+n⋅r，当 r r r较小时，数据占用内存大大减少，因此还可以用奇异值分解进行无损压缩矩阵。

截断奇异值分解

透过现象看本质，根据前面奇异值分解的公式，我们可以推导出，每一个矩阵都可以分解成如下形式：
A = σ 1 ⋅ A 1 + σ 2 ⋅ A 2 + . . . + σ r ⋅ A r A=sigma_1cdot A_1+sigma_2cdot A_2+...+sigma_rcdot A_r A=σ1⋅A1+σ2⋅A2+...+σr⋅Ar
这里，暂且不管 A i , i = 1 , 2 , . . . , r A_i,i=1,2,...,r Ai,i=1,2,...,r是什么（可以推，太麻烦了不好写），可以很明显的看出 A A A可以由 σ i sigma_i σi为权值通过矩阵叠加组成，因此，当某些 σ i sigma_i σi很小时，对组成 A A A的贡献很小，可以看成 A A A的噪声，将其去除，这样就得到了截断奇异值分解，其矩阵形式如下:
A ≈ U k Σ k V k T Aapprox U_kSigma_kV_k^T A≈UkΣkVkT
k k k大小等于保留的奇异值个数，由于奇异值矩阵 Σ Sigma Σ是按从大到小排列的，因此 U k U_k Uk表示取 U U U的前 k k k列， Σ k Sigma_k Σk表示取 Σ Sigma Σ的前 k k k行 k k k列所组成的方阵， V k V_k Vk表示取 V V V的前 k k k列。
根据上述思想，我们可以对矩阵通过截断奇异值分解进行有损压缩，去噪等。

代码代码构架

为了实现SVD，根据上述公式可以得出，我们需要用到的工具有：

盛放矩阵的数据结构；（在此我选用了vector容器）

vector> A = {{1,0,0,0},{0,0,0,4},{0,3,0,0},{0,0,0,0},{2,0,0,0}};

矩阵乘法运算；

// 矩阵乘法
template
vector> matrix_multiply(vector> const arrL, vector> const arrR)
{
    int rowL = arrL.size();// 左矩阵行数
    int colL = arrL[0].size();// 左矩阵列数
    int rowR = arrR.size();// 右矩阵行数
    int colR = arrR[0].size();// 右矩阵列数
    // 判断是否能够相乘
    if(colL != rowR)
    {
        throw "left matrix's row not should equal with right matrix!";
    }
    // initialize result matrix 
    vector> res(rowL);
    for(int i=0; i
        res[i].resize(colR);
    }
    // compute matrix multiplication
    for(int i=0; i
        for(int j=0; j
            for(int k=0; k
                res[i][j] += arrL[i][k]*arrR[k][j];
            }
        }
    }

    return res;

}

矩阵转置；

// 矩阵转置
template
vector> transpose(vector> const arr)
{
    int row = arr.size();
    int col = arr[0].size();
    // initialize transpose matrix col*row
    vector> trans(col);
    for(int i=0;i
        trans[i].resize(row);
    }
    // fill elements
    for(int i=0; i
        for(int j=0;j
            trans[i][j] = arr[j][i];
        }
    }
    return trans;
}

实对称矩阵特征值特征向量求解（选用Jacobi迭代法）

//提前声明后续用到的argsort函数，功能类似于numpy的那个
template 
vector argsort(const vector& array);
// 实对称矩阵特征值特征向量
// param: arr   :input array
// param: E     :eigen vectors
// param: e     :eigen values
template
void eigen(vector> arr, vector> &E, vector &e)
{
    //vector> arr = arr_ori;
    int n = arr.size();// size of matrix
    int row = 0;// row index of max
    int col = 0;// col index of max
    int iter_max_num = 10000;//迭代总次数
    int iter_num = 0;
    double eps = 1e-40;//误差
    double max = eps;// 非对角元素最大值
    // 初始化特征向量矩阵为单位阵,初始化特征值
    E.resize(n);
    e.resize(n);
    for(int i=0; i
        E[i].resize(n,0);
        E[i][i] = 1;
    }

    while(iter_num=eps)
    {
        max = fabs(arr[0][1]);
        row = 0;
        col = 1;
        // find max value and index
        for(int i=0;i
            for(int j=0;j
                if(i!=j && fabs(arr[i][j])>max){
                    max = fabs(arr[i][j]);
                    row = i;
                    col = j;
                }
            }
        }
        double theta = 0.5*atan2(2 * arr[row][col] , (arr[row][row] - arr[col][col]));
        //update arr
        double aii = arr[row][row];
        double ajj = arr[col][col];
        double aij = arr[row][col];
        double sin_theta = sin(theta);
        double cos_theta = cos(theta);
        double sin_2theta = sin(2 * theta);
        double cos_2theta = cos(2 * theta);
        arr[row][row] = aii*cos_theta*cos_theta + ajj*sin_theta*sin_theta + aij*sin_2theta;//Sii'
        arr[col][col] = aii*sin_theta*sin_theta + ajj*cos_theta*cos_theta - aij*sin_2theta;//Sjj'
        arr[row][col] = 0.5*(ajj - aii)*sin_2theta + aij*cos_2theta;//Sij'
        arr[col][row] = arr[row][col];//Sji'
        for (int k = 0; k < n; k++)
        {
            if (k != row && k != col)
            {
                double arowk = arr[row][k];
                double acolk = arr[col][k];
                arr[row][k] = arowk * cos_theta + acolk * sin_theta;
                arr[k][row] = arr[row][k];
                arr[col][k] = acolk * cos_theta - arowk * sin_theta;
                arr[k][col] = arr[col][k];
            }
        }
        // update E
        double Eki;
        double Ekj;
        for(int k=0; k
            Eki = E[k][row];
            Ekj = E[k][col];
            E[k][row] = Eki*cos_theta + Ekj*sin_theta;
            E[k][col] = Ekj*cos_theta - Eki*sin_theta;
        }
        iter_num++;
    }
    //update e
    for(int i=0;i
        e[i] = arr[i][i];
    }
    // sort E by e
    vector sort_index;
    sort_index = argsort(e);
    // initialize E_sorted, e_sorted
    vector> E_sorted(n);
    for(int i=0;i
        E_sorted[i].resize(n);
    }
    vector e_sorted(n);
    for(int i=0;i
        e_sorted[i] = e[sort_index[i]];
        for(int j=0;j
            E_sorted[i][j] = E[i][sort_index[j]];
        }
    }
    E = E_sorted;
    e = e_sorted;
    //delete &E_sorted, &e_sorted;
    cout<<"max element is: "< 
在实现完上述工具后，即可实现SVD类，其构造如下： 
//
//########################SVD##########################
//
//params:
//      --arr   :input matrix m*n
//      --U     :left matrix m*r , r <= rank(arr)
//      --S     :medium matrix r*r 
//      --V     :right matrix n*r
class SVD
{
    public:
        vector> U,S,V,ATA,A;
        int n,m,r;
        vector> E; //特征向量矩阵
        vector e; // 特征值向量
        SVD(vector> arr);
        void tight_svd();//紧奇异值分解
        void truncated_svd(int);//截断奇异值分解
    
};
SVD::SVD(vector> arr){
    m = arr.size();
    n = arr[0].size();
    A = arr;
    ATA = matrix_multiply(transpose(A),A);
    // 计算ATA特征值特征向量
    eigen(ATA,E,e);
}
void SVD::tight_svd(){
    r = 0;
    // 确定秩
    for(int i=0;i
        if(e[i]>1e-10){
            r++;
        }
        else break;
    }
    //确定V
    V = E;
    for(int i=0; i
        V[i].resize(r);
    }
    //确定S
    S.resize(r);
    for(int i=0;i
        S[i].resize(r);
        S[i][i] = sqrt(e[i]);
    }
    //确定U
    vector> Sinv = S;
    for(int i=0;i
        Sinv[i][i] = 1/S[i][i];
    }
    U = matrix_multiply(matrix_multiply(A,V),Sinv);
}
void SVD::truncated_svd(int rr){
    r = rr;
    //确定V
    V = E;
    for(int i=0; i
        V[i].resize(r);
    }
    //确定S
    S.resize(r);
    for(int i=0;i
        S[i].resize(r);
        S[i][i] = sqrt(e[i]);
    }
    //确定U
    vector> Sinv = S;
    for(int i=0;i
        Sinv[i][i] = 1/S[i][i];
    }
    U = matrix_multiply(matrix_multiply(A,V),Sinv);
}

 
实验测试 
以《统计机器学习》P277页例15.2为例
  
代码如下： 
int main()
{
    vector> A = {{1,0,0,0},{0,0,0,4},{0,3,0,0},{0,0,0,0},{2,0,0,0}};
    SVD svd(A);
    // tight SVD
    svd.tight_svd();
    cout< 
测试结果
 
 完整代码见我的GitHub仓库：https://github.com/wjtgoo/SVD-CPP（万水千山总是情，给个星星行不行/(ㄒoㄒ)/~~）

SVD 奇异值分解纯手工实现（C++）

C/C++/C#相关栏目本月热门文章