SVD奇异值分解（理论与C++实现）

前言

  奇异值分解(singular value decomposition，以下简称SVD)是线性代数中一种重要的矩阵分解。SVD将矩阵分解为奇异向量(singular vector)和奇异值(singular value)。SVD将矩阵$A$分解成三个矩阵的乘积
$$A=UD{V}^T$$
  设$A$是$m\times n$的矩阵，则$U$是一个$m\times m$的矩阵，$D$是一个$m\times n$的矩阵$V$是一个$n\times n$的矩阵。
  读者可访问以下链接获取相关代码
  代码链接

理论推导

部分代码实现

/**
 * @brief 实现奇异值分解
 * @note  基于Householder变换以及变星QR算法对一般实矩阵A进行奇异值分解
 * step1 用豪斯荷尔德变换将A约化为双对角线矩阵
 * step2 用变星的QR算法进行迭代，计算所有奇异值
 * steo3 对奇异值按非递增次序进行排列
*/
void Matrix::SVD( const Matrix & A, Matrix& U2, Matrix& W, Matrix& V) {
   
    Matrix U = A;
    int m = A.rows();
    int n = A.cols();
    U2 = Matrix(m, m);
    V = Matrix(n, n);

    NAFLOAT* w = (NAFLOAT*)malloc(n * sizeof(NAFLOAT));
    NAFLOAT* rv1 = (NAFLOAT*)malloc(n * sizeof(NAFLOAT));

    int32_t flag, i, its, j, jj, k, l, nm;
    NAFLOAT   anorm, c, f, g, h, s, scale, x, y, z;

    g = scale = anorm = 0.0;
    for (i = 0; i < n; i++) {
   
        l = i + 1;
        rv1[i] = scale * g;
        g = s = scale = 0.0;
        if (i < m) {
   
            for (k = i; k < m; k++) scale += fabs(U.m_val[k][i]);
            if (scale) {
   
                for (k = i; k < m; k++) {
   
                    U.m_val[k][i] /= scale;
                    s += U.m_val[k][i] * U.m_val[k][i];
                }
                f = U.m_val[i][i];
                g = -SIGN(sqrt(s), f);
                h = f * g - s;
                U.m_val[i][i] = f - g;
                for (j = l; j < n; j++) {
   
                    for (s