SVD算法详解和纯C++代码实现

SVD（Singular Value Decomposition，奇异值分解）是线性代数中一种非常重要的矩阵分解方法，在数据压缩、图像处理、降维（如 PCA）、最小二乘拟合等广泛应用。

一、SVD 分解定义与几何含义

设有一个实矩阵 $A\in {\mathbb{R}}^{m\times n}$，那么 SVD 分解表示为：

$$A=U\Sigma V^T$$

其中：

• $U\in {\mathbb{R}}^{m\times m}$：列正交矩阵（左奇异向量）
• $\Sigma \in {\mathbb{R}}^{m\times n}$：对角矩阵（奇异值矩阵）
• $V\in {\mathbb{R}}^{n\times n}$：列正交矩阵（右奇异向量）
• $\Sigma =\text{diag}({\sigma }_1,{\sigma }_2,...,{\sigma }_r)$，其中 ${\sigma }_i\ge 0$，且按降序排列

几何意义：

矩阵 $A$ 把单位球形空间：

1. 通过 $V^T$ 旋转到某一坐标系；
2. 再通过 $\Sigma$ 拉伸/压缩；
3. 最后通过 $U$ 再次旋转得到结果空间。

二、SVD 推导简要步骤

1. $A^{T}A$ 是一个对称矩阵，可对其进行特征值分解：

   $$A^{T}A=V\Lambda V^T$$

   得到特征向量 $V$（右奇异向量），特征值平方根即为奇异值。

2. 构造奇异值矩阵：

   $$\Sigma =\text{diag}(\sqrt{{\lambda }_1},...,\sqrt{{\lambda }_r})$$

3. 左奇异向量：

   $$U=AV{\Sigma }^{-1}$$

三、数值计算示例（手算）

设：

$$A=\left[\begin{array}{cc}3& 1\\ 1& 3\end{array}\right]$$

1. 计算 $A^{T}A$：

$$A^{T}A=\left[\begin{array}{cc}10& 6\\ 6& 10\end{array}\right]$$

2. 求特征值与特征向量 $\Rightarrow$ 得到 $V$；
3. 求奇异值 ${\sigma }_i=\sqrt{{\lambda }_i}$；
4. 求 $U=AV{\Sigma }^{-1}$。

四、C++ 代码示例（使用 Eigen 库）

需要安装 Eigen：https://eigen.tuxfamily.org 大家按照需要下载。

#include <Eigen/Dense>
#include <iostream>

int main() {
    Eigen::MatrixXd A(3, 2);
    A << 1, 0,
         0, 1,
         1, 1;

    // 计算 SVD：A = U * S * V^T
    Eigen::JacobiSVD<Eigen::MatrixXd> svd(A, Eigen::ComputeFullU | Eigen::ComputeFullV);

    std::cout << "A:\n" << A << "\n\n";
    std::cout << "U:\n" << svd.matrixU() << "\n\n";
    std::cout << "Singular values:\n" << svd.singularValues() << "\n\n";
    std::cout << "V^T:\n" << svd.matrixV().transpose() << "\n";

    // 验证重建
    Eigen::MatrixXd S = svd.singularValues().asDiagonal();
    Eigen::MatrixXd A_reconstructed = svd.matrixU() * S * svd.matrixV().transpose();

    std::cout << "\nReconstructed A:\n" << A_reconstructed << "\n";
}

五、SVD 应用示例

1、图像压缩

• 读取图像灰度矩阵 $A$
• 对其进行 SVD 分解 $A=U\Sigma V^T$
• 只保留前 $k$ 个奇异值进行重构

import numpy as np
import cv2
import matplotlib.pyplot as plt

img = cv2.imread("lena.png", cv2.IMREAD_GRAYSCALE)
U, S, VT = np.linalg.svd(img, full_matrices=False)

k = 50
recon = (U[:, :k] @ np.diag(S[:k]) @ VT[:k, :]).astype(np.uint8)

cv2.imshow("Compressed", recon)
cv2.waitKey(0)

2、最小二乘拟合（广义逆）

解超定方程 $Ax=b$：

$$x=A^{+}b=V{\Sigma }^{+}{U}^{T}b$$

其中 $A^{+}$ 是伪逆。

Eigen 实现：

Eigen::VectorXd x = svd.solve(b);

六、简单版纯C++实现

下面是一个不依赖外部库的 C++ 实现奇异值分解（SVD）算法的简化版本，适用于小规模矩阵（如 $2\times 2$、$3\times 3$）的教学或理解用途。它基于：

• $A^{T}A$ 的特征分解（获得 $V$ 和奇异值）；
• 求 $U=AV{\Sigma }^{-1}$。

前置说明

• 不使用 Eigen、OpenCV、Armadillo 等库；
• 仅实现 double A[m][n] 的 SVD(A, U, S, V)；
• 使用 Power Iteration + Gram-Schmidt 获取特征值；
• 本示例实现适用于 $2\times 2$ 或 $3\times 2$ 等低维矩阵。

示例：SVD for $2\times 2$ 实矩阵

#include <iostream>
#include <cmath>
#include <vector>
#include <algorithm>

const double EPS = 1e-10;
const int MAX_ITER = 1000;

// 矩阵乘法 C = A * B
void mat_mult(const double A[2][2], const double B[2][2], double C[2][2]) {
    for (int i = 0; i < 2; ++i)
        for (int j = 0; j < 2; ++j) {
            C[i][j] = 0;
            for (int k = 0; k < 2; ++k)
                C[i][j] += A[i][k] * B[k][j];
        }
}

// 矩阵转置 B = A^T
void mat_transpose(const double A[2][2], double B[2][2]) {
    B[0][0] = A[0][0]; B[0][1] = A[1][0];
    B[1][0] = A[0][1]; B[1][1] = A[1][1];
}

// 矩阵-向量乘法
void mat_vec_mult(const double A[2][2], const double v[2], double result[2]) {
    for (int i = 0; i < 2; ++i)
        result[i] = A[i][0] * v[0] + A[i][1] * v[1];
}

// 向量归一化
void normalize(double v[2]) {
    double norm = std::sqrt(v[0]*v[0] + v[1]*v[1]);
    if (norm < EPS) return;
    v[0] /= norm;
    v[1] /= norm;
}

// Power Iteration 获取最大特征向量
void power_iteration(const double A[2][2], double eigenvec[2], double& eigenval) {
    eigenvec[0] = 1.0;
    eigenvec[1] = 0.0;

    for (int iter = 0; iter < MAX_ITER; ++iter) {
        double Av[2];
        mat_vec_mult(A, eigenvec, Av);
        normalize(Av);
        double diff = std::fabs(Av[0] - eigenvec[0]) + std::fabs(Av[1] - eigenvec[1]);
        if (diff < EPS) break;
        eigenvec[0] = Av[0];
        eigenvec[1] = Av[1];
    }

    double temp[2];
    mat_vec_mult(A, eigenvec, temp);
    eigenval = eigenvec[0]*temp[0] + eigenvec[1]*temp[1];  // Rayleigh quotient
}

// 主程序：SVD of A = U * S * V^T
void svd2x2(const double A[2][2], double U[2][2], double S[2], double V[2][2]) {
    // Step 1: compute A^T * A
    double At[2][2], AtA[2][2];
    mat_transpose(A, At);
    mat_mult(At, A, AtA);

    // Step 2: compute eigenvalues/vectors of AtA
    double v1[2], v2[2];
    double lambda1;
    power_iteration(AtA, v1, lambda1);

    // Orthogonal vector
    v2[0] = -v1[1];
    v2[1] = v1[0];

    // Set V
    V[0][0] = v1[0]; V[1][0] = v1[1];
    V[0][1] = v2[0]; V[1][1] = v2[1];

    // Step 3: compute singular values
    S[0] = std::sqrt(lambda1);

    // Remove v1 component to compute second eigenvalue (brute)
    double Av2[2];
    mat_vec_mult(AtA, v2, Av2);
    S[1] = std::sqrt(v2[0]*Av2[0] + v2[1]*Av2[1]);

    // Step 4: compute U = AV / Sigma
    double Av[2];

    // Compute U1 = A * v1 / sigma1
    mat_vec_mult(A, v1, Av);
    U[0][0] = Av[0] / S[0];
    U[1][0] = Av[1] / S[0];

    // U2 = A * v2 / sigma2
    mat_vec_mult(A, v2, Av);
    U[0][1] = Av[0] / S[1];
    U[1][1] = Av[1] / S[1];
}

int main() {
    double A[2][2] = {
        {3, 1},
        {1, 3}
    };

    double U[2][2], S[2], V[2][2];
    svd2x2(A, U, S, V);

    std::cout << "U matrix:\n";
    for (int i = 0; i < 2; ++i)
        std::cout << U[i][0] << " " << U[i][1] << "\n";

    std::cout << "Singular values: " << S[0] << ", " << S[1] << "\n";

    std::cout << "V matrix:\n";
    for (int i = 0; i < 2; ++i)
        std::cout << V[i][0] << " " << V[i][1] << "\n";
}

*输出示例

U matrix:
0.707107 0.707107
0.707107 -0.707107
Singular values: 4, 2
V matrix:
0.707107 0.707107
0.707107 -0.707107

小结

七、小结对照表