简介：这是一份用C语言实现SVD奇异值分解的代码包，适合需要掌握矩阵分解底层原理的算法学习者、C/C++开发者以及数据科学相关从业者。代码涵盖了特征值与特征向量计算、特征值排序、构造奇异值矩阵Σ以及标准正交矩阵U、V等核心步骤，并包含可执行文件，便于直接运行验证分解结果。压缩包共80个文件、大小约1.12MB，其中主要包含svd.cpp、main.cpp、svd.h等源码文件，以及Visual Studio工程配置、编译日志和说明文档，结构清晰，方便在Windows环境下打开编译或对照学习。已有1091人学习下载。通过这份代码，读者不仅能理解A=UΣV^T的完整实现流程，还能了解数值稳定性处理思路，并进一步迁移到PCA主成分分析、推荐系统协同过滤、图像压缩等典型应用场景中，是学习SVD算法不错的参考实现。
1. 奇异值分解的 C 代码落地：为什么这件事值得你亲自写一遍
提到奇异值分解，很多人的第一反应是 MATLAB 里一行svd(A)，或者 Python 里np.linalg.svd。可一旦场景换到嵌入式设备、实时信号处理，或者需要把某个算法固化到自研框架里，这些高层接口统统不适用。你需要的是一份能看懂、能改、能死磕的 C 实现。这篇文章要解决的问题很直接：在不依赖第三方数值库的前提下，用纯 C 把奇异值分解写出来，并且保证它在雅可比迭代下能收敛、能排序、能重构出原始矩阵。
这个方向适合谁？如果你正在做图像压缩、推荐系统里的矩阵降维、惯性导航里的姿态估计，或者只是想把《数值分析》课本上的奇异值分解真正跑起来，那这份代码就是你的起点。真正的工程价值不在于“能算出奇异值”，而在于理解迭代过程里每一个容差、每一次旋转换行，以及那些让你怀疑人生的收敛失败到底从哪来。接下来我会从算法选型开始，一路走到完整代码、验证方法和实际应用，最后把坑都摊开讲。
2. 奇异值分解的算法选型：双边雅可比和 Householder 到底该选谁
2.1 首先明确你要分解的是什么矩阵
在写任何一行代码之前，必须搞清楚一个基本问题：你要做的是方阵的还是一般矩阵的奇异值分解。对于方阵，我们可以直接做特征值分解；但对于m×n的矩形矩阵，标准做法是先把它缩小成方阵，再用对称特征值问题的思路处理。
这里提供一个稳定的化简路径：对于任意一个m×n矩阵 A，构造一个对称矩阵A^T·A。这个矩阵是n×n方阵，它的特征值就是 A 的奇异值的平方，特征向量就是右奇异向量 V。之后再用U = A·V·S^{-1}恢复出左奇异向量 U。这个做法的好处是——你只需要写一个对称矩阵的雅可比特征值分解，剩下的工作就是矩阵乘法和除法。坏处也很明显——A^T·A的条件数是原矩阵的平方，如果 A 本身就病态，这个中间步骤会放大误差，导致小奇异值精度严重丢失。
所以我在实际工程里有一条明确的边界线：如果矩阵的行数和列数都在几百以内，并且你能够容忍毫秒级的计算时间，直接用双边雅可比法处理A^T·A是简单可靠的路径；如果矩阵维度上千，或者奇异值动态范围超过 1e8，我宁可多写一个 Householder 双对角化。前者是这篇文章的主角，后者我会在最后一章给出迁移路径。
2.2 双边雅可比旋转的核心思想：逐对消除非对角元
双边雅可比法的思路非常朴素：通过一系列平面旋转矩阵，不断把对称矩阵的非对角元素“磨”小，直到所有非对角元的绝对值都小于某个容差。和普通雅可比迭代不同，双边雅可比每次旋转同时作用于左右两侧，也就是说「A ← J·A·J^T」，并且用旋转角度同时更新 U 和 V。
这里的关键参数就是旋转角度的计算。设当前要消去的元素是A[p][p]、A[q][q]和A[p][q]，旋转角度 θ 满足：
tan(2θ) = 2·A[p][q] / (A[q][q] - A[p][p])这个式子对应教科书的推导过程。在实际代码里，我们不直接用 tan，而是计算 c = cos(θ) 和 s = sin(θ)，然后构造旋转变换。如果你的矩阵是对称的，这个旋转会直接把A[p][q]位置清零（在浮点精度内）；但是后续的旋转又会把它重新填上一点非零值，这就是“迭代磨”的过程。
一个常见误区是：很多人以为一轮扫描就能把所有非对角元消干净。实际上雅可比迭代是多次扫描，每一轮扫描所有p < q的对，直到收敛。收敛判据和最大迭代次数是必须手动设置的，这也是后续调试时最容易出问题的地方。
2.3 为什么我最终选择了稳健的雅可比实现
对比一下主流路线，LAPACK 里用的dgesdd走的是 divide-and-conquer，速度快但代码量超过几千行，不适合教学和我们这种“要藏在单文件库里”的场景。GSL（某开源科学计算库）虽然有现成接口，但引入外部库违背了我们“纯 C 单文件”的假定。
双边雅可比实现的核心优势有三个。第一，逻辑极其透明，每迭代一轮你都能打印出中间矩阵观察收敛过程；第二，对于中小型稠密矩阵，它的精度是可控的，适当增加迭代次数完全可以把重构误差压到 1e-12 以下；第三，代码量小，不加注释约 150 行，便于移植到 DSP 或者 MCU 平台。
代价是收敛速度。当矩阵维度超过 200 之后，平方复杂度会让你等得不耐烦。所以再一次强调：这份代码是为“中小型矩阵、需要深度定制”的场景准备的，不是为大规模数值计算准备的。
3. 手写奇异值分解 C 代码：算法步骤、完整实现与关键参数
3.1 完整代码骨架：一个可运行的 C 语言实现
我这里给出一份可以直接编译运行的双边雅可比奇异值分解实现。它针对的是对称矩阵A^T·A的特征值问题，从而推导出奇异值和奇异向量。代码全部集中在三个函数里：一个做对称矩阵特征分解，一个组装 U 矩阵，一个做完整奇异值分解。
#include <stdio.h> #include <math.h> #define SVD_MAX_ITER 100 #define SVD_EPSILON 1e-10 // 对称矩阵的雅可比特征值分解 // a: 输入 n×n 对称矩阵，输出时被覆盖为对角矩阵 // v: 输出 n×n 特征向量矩阵（列向量是特征向量） // n: 矩阵维度 // 返回实际迭代轮数 int jacobi_eigen(double **a, double **v, int n) { for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { v[i][j] = (i == j) ? 1.0 : 0.0; } } int iter_count = 0; for (int iter = 0; iter < SVD_MAX_ITER; iter++) { // 寻找最大的非对角元素 double max_offdiag = 0.0; int p = 0, q = 1; for (int i = 0; i < n; i++) { for (int j = i + 1; j < n; j++) { if (fabs(a[i][j]) > max_offdiag) { max_offdiag = fabs(a[i][j]); p = i; q = j; } } } if (max_offdiag < SVD_EPSILON) { break; // 收敛 } // 计算旋转角度 double app = a[p][p], aqq = a[q][q], apq = a[p][q]; double theta = 0.5 * atan2(2.0 * apq, aqq - app); double c = cos(theta), s = sin(theta); // 更新 a 矩阵（对称旋转） for (int k = 0; k < n; k++) { double a_kp = a[k][p]; double a_kq = a[k][q]; a[k][p] = c * a_kp - s * a_kq; a[k][q] = s * a_kp + c * a_kq; a[p][k] = a[k][p]; // 保持对称 a[q][k] = a[k][q]; } // 更新特征向量矩阵 v for (int k = 0; k < n; k++) { double v_kp = v[k][p]; double v_kq = v[k][q]; v[k][p] = c * v_kp - s * v_kq; v[k][q] = s * v_kp + c * v_kq; } iter_count++; } return iter_count; } // 从对称矩阵 A^T·A 分解得到奇异值和左右奇异向量 // a: 输入 m×n 矩阵，会被覆盖 // u: 输出 m×n 左奇异向量矩阵 // s: 输出 n 维奇异值数组 // v: 输出 n×n 右奇异向量矩阵 // m, n: 输入矩阵维度 void svd_symmetric(double **a, double **u, double *s, double **v, int m, int n) { // 第一步：构造 A^T·A double **ata = (double **)malloc(n * sizeof(double *)); for (int i = 0; i < n; i++) { ata[i] = (double *)calloc(n, sizeof(double)); } for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { double sum = 0.0; for (int k = 0; k < m; k++) { sum += a[k][i] * a[k][j]; } ata[i][j] = sum; } } // 第二步：对 A^T·A 做特征分解 double **v_temp = (double **)malloc(n * sizeof(double *)); for (int i = 0; i < n; i++) { v_temp[i] = (double *)malloc(n * sizeof(double)); } jacobi_eigen(ata, v_temp, n); // 第三步：提取奇异值（特征值开方） for (int i = 0; i < n; i++) { s[i] = sqrt(fabs(ata[i][i])); for (int j = 0; j < n; j++) { v[i][j] = v_temp[j][i]; // 转置 } } // 第四步：恢复左奇异向量 U = A·V·S^{-1} for (int i = 0; i < m; i++) { for (int j = 0; j < n; j++) { double sum = 0.0; for (int k = 0; k < n; k++) { sum += a[i][k] * v[k][j]; } u[i][j] = (s[j] > SVD_EPSILON) ? sum / s[j] : 0.0; } } // 清理 for (int i = 0; i < n; i++) { free(ata[i]); free(v_temp[i]); } free(ata); free(v_temp); }3.2 代码逻辑说明：每一步到底在计算什么
这个实现里有几个关键点需要逐个讲清楚。第一步构造A^T·A，这是把矩形矩阵转换成方阵的必经桥段。注意我们是按行展开的，所以内层循环是a[k][i] * a[k][j]，这对应矩阵乘法的标准形式。如果用列优先存储，循环顺序需要调整，否则性能会差上好几个量级。
第二步调用jacobi_eigen函数，它的作用是把对称矩阵消成对角形。核心是那个atan2(2.0 * apq, aqq - app)—— 这里用atan2而不是atan，是为了自动处理app == aqq时角度为 45 度的边界情况，避免除零溢出。同时它把角度自动映射到正负 π 的范围，数值上比手动判断象限更稳。
第三步提取奇异值。这里有一个微妙的坑：理论上A^T·A的特征值都是非负的，但浮点计算会产生微小的负值，尤其当某个奇异值接近零时。所以我用sqrt(fabs(...))，并且这一步和后续的S^{-1}相互配合——当一个奇异值小于SVD_EPSILON时，直接把对应分量置零，不参与除法。这个处理方式叫“伪逆语义”，它保证了算法不会因为除以零而崩溃。
第四步恢复 U 矩阵用的是u = A·V·S^{-1}这个公式。这一步是 O(m·n^2) 的复杂度，所以对超大矩阵会偏慢；但在几百维之内完全无感。注意这里s[j] > SVD_EPSILON用的是绝对阈值，如果你的奇异值尺度本身就很大（比如 1e10 量级），需要把阈值改成相对阈值，我后面避坑章节会专门讲。
3.3 关键参数：三个最容易影响结果的数字
如果你只是“能跑起来”就满足，那这段代码已经够了。但工程上我们需要知道三个参数的调节逻辑，这也是面试和实际调试中的高频考点。
第一是SVD_MAX_ITER。它控制最大扫描轮数。100 轮对绝大多数中小矩阵都够了，因为雅可比迭代是平方收敛的，即每轮之后误差大约变成上一轮的平方。但如果你遇到高维的近奇异矩阵，100 轮可能不够，你会看到返回后非对角元还大得离谱。解决办法是把上限提到 500，同时加一个打印函数观察每轮的最大非对角元，判断是不是真的在下降。如果 500 轮都不降，基本可以断定矩阵病态到需要预处理。
第二是SVD_EPSILON。它同时承担两个职责：收敛阈值和零奇异值判定。1e-10 这个数值在双精度下比较稳妥，但如果你的数据本身噪声很大，设得太小没有意义。一个经验值是SVD_EPSILON = 1e-12 * max(m,n) * max(|A|)，这样可以适应性判断。用相对阈值而不是绝对阈值，能避开不同量级数据带来的误判。
第三是内存布局。我这里用的是double **指针数组形式，优点是写起来直观，缺点是多了一次间接寻址，缓存命中率差。如果你需要高性能，改成连续的一维数组double *a = (double *)calloc(m*n, sizeof(double))，索引用a[i*n + j]代替。这个改动很简单，但性能差距通常在 30% 以上。
4. 奇异值分解 C 代码的四个边界坑与排查经验
4.1 坑一：收敛后奇异值顺序混乱，导致 U 和 V 对不上
现象：分解完成后，打印奇异值s[0] > s[1]，但用U·S·V^T重构矩阵时误差巨大，甚至完全不同。
原因：雅可比迭代产出的特征值本来就是乱序的，因为每次消去非对角元时没有排序约束。特征值在矩阵对角上的位置取决于旋转顺序，所以最后得到的奇异值序列是任意排列的。如果你接下来的代码里假定奇异值降序排列，后续的截断和伪逆计算就会全部错乱。
解决：在拿到s数组后做一次冒泡排序或者插入排序，同时对v和u的列执行同步交换。注意这里必须同步操作三个数组，只排奇异值不排向量是最常见的翻车姿势。一个干净的实现是写一个sort_svd_result(u, s, v, n)函数，用下标索引数组记录排列顺序，然后一次性重排所有数据。
void sort_svd_result(double **u, double *s, double **v, int m, int n) { for (int i = 0; i < n - 1; i++) { int max_idx = i; for (int j = i + 1; j < n; j++) { if (s[j] > s[max_idx]) { max_idx = j; } } if (max_idx != i) { // 交换奇异值 double tmp_s = s[i]; s[i] = s[max_idx]; s[max_idx] = tmp_s; // 交换 U 的列 for (int k = 0; k < m; k++) { double tmp = u[k][i]; u[k][i] = u[k][max_idx]; u[k][max_idx] = tmp; } // 交换 V 的列 for (int k = 0; k < n; k++) { double tmp = v[k][i]; v[k][i] = v[k][max_idx]; v[k][max_idx] = tmp; } } } }4.2 坑二：A^T·A 的条件数爆炸，小奇异值精度完全丢失
现象：矩阵本身数值正常，但分解出来的最小奇异值比真值大了几个数量级，或者干脆是负的（被fabs掩藏了符号）。
原因：这是A^T·A路线的天然缺陷。矩阵乘法会把条件数从κ(A)变成κ(A)^2。如果 A 的条件数是 1e6，那么A^T·A的条件数就是 1e12，接近双精度浮点的极限。此时最小奇异值已经淹没在舍入误差里。
解决：这是算法的结构性短板，不是代码 bug。如果你对最小奇异值精度有硬性要求，就不要走这条路线，改用 Householder 双对角化加隐式 QR 迭代。如果只是做图像压缩里的截断操作，影响不大，因为截断只保留前几个大的奇异值。但做矩阵求逆或最小二乘时，小奇异值直接决定解的稳定性，这时候必须换路线。
4.3 坑三：非对角元判定阈值导致死循环或者提前退出
现象：设置SVD_EPSILON = 1e-10，结果矩阵维度稍微大一点（比如 50×50），算法跑了非常久；换一个小一点的矩阵，又很快跳出迭代，重构误差却不达标。
原因：阈值设置不是尺寸无关的。你的收敛判据是最大非对角元小于固定值，但矩阵元素的绝对量级会影响这个判断。如果矩阵元素本身是 1e6 量级，非对角元降到 1e-10 需要极端多的轮次；如果矩阵元素是 1e-8 量级，这个阈值又太松了。
解决：改用混合判据——同时要求最大非对角元小于SVD_EPSILON * max(|A_pq|)，并且小于一个绝对下限。同时在迭代函数里加入“最大迭代次数触发退出”的日志。我通常的做法是：每次迭代打印iter和max_offdiag，观察收敛曲线是否平滑下降。如果曲线在某个值上震荡，说明矩阵是近奇异的，这时候该做的是预处理而不是调阈值。
4.4 坑四：左奇异向量 U 的符号漂移，重构正确但语义错误
现象：U·S·V^T能完美重构 A，但拿 U 的某一列去做物理意义分析（比如姿态估计里的旋转轴），发现方向和预期差了 180 度。
原因：奇异值分解不是唯一的。你可以在任意一列 U 和对应列 V 上同时乘以 -1，等式依然成立。雅可比迭代的产物依赖于初始权矩阵和浮点舍入，符号是随机的。
解决：这不算 bug，但需要约定。在做姿态估计或子空间分析时，我会在分解完成后检查 U 的第一列或者某一指定列的符号，如果某个参考元素为负，就把该列和 V 的对应列一起取反。
5. 验证奇异值分解代码的三种手段：从单元测试到矩阵重构
5.1 手段一：用 A = U·S·V^T 重构误差来验证
这是最直接的验证方式。把分解得到的 U、S、V 重新乘起来，和原始矩阵比，计算最大绝对误差和 Frobenius 范数误差。对应的代码可以这么写：
double max_error = 0.0; double frob_error = 0.0; for (int i = 0; i < m; i++) { for (int j = 0; j < n; j++) { double recons = 0.0; for (int k = 0; k < n; k++) { recons += u[i][k] * s[k] * v[j][k]; } double err = fabs(recons - a_original[i][j]); if (err > max_error) max_error = err; frob_error += err * err; } } frob_error = sqrt(frob_error); printf("max error: %e, frob error: %e\n", max_error, frob_error);这里的关键是，重构时必须用v[j][k]而不是v[k][j]，取决于你设计矩阵时 V 的存储格式。我的实现里 V 是n×n，第 j 行第 k 列对应第 k 个右奇异向量的第 j 个分量，所以乘法索引是v[j][k]。如果不小心转置了，误差会大得离谱，这是非常容易犯的低级错误。
另一个细节是s[k]是标量数组，不需要展开成矩阵。这个重构代码本身就是最好的调试工具——误差一大，先用这个小程序定位，别急着改算法。
5.2 手段二：对比已知解析解验证正确性
拿一个已知奇异值和奇异向量的测试矩阵来验证。比如一个 2×2 对角矩阵diag(3, 1)，它的奇异值就是 3 和 1，左奇异向量和右奇异向量都是单位阵。如果代码输出不是这个结果，那基本上可以判定方位的索引有问题。
再比如一个秩为 1 的矩阵：A = u0·v0^T，其中 u0 和 v0 是任意单位向量。这个矩阵的奇异值应该是[1, 0]，U 的第一列是 u0，V 的第一列是 v0。这个用例特别擅长暴露零点处理的 bug——因为第二个奇异值为零，很多实现对零的处理是直接崩掉或者产生 NaN。
我在实际工作中会把这两组测试做成一个test_svd.c，每次改动代码就跑一遍，确保回归不破坏已有功能。这年头光靠“看着没问题”来交付代码，迟早要被线上事故教育。
5.3 手段三：和现成数值库做交叉验证
如果你手头有其他语言环境，可以用它们的结果当“标准答案”。比如在 Python 里用同样的矩阵做numpy.linalg.svd，然后把结果导出成文本，和你的 C 代码输出对比。对比时别比绝对值，要比相对误差，因为两个实现可能采用了不同的迭代路径，产生不同的舍入误差。
我一般写一个脚本读两个文件，输出每一列的奇异值相对误差和最大分量误差。如果相对误差在 1e-8 以内，基本可以确认实现是对的。如果超出这个范围，优先检查是否排序不一致或者符号漂移。这个方法做一次普通验证只需要 5 分钟，但能省下后面几天找 bug 的时间，是我的固定习惯。
6. 奇异值分解的进阶玩法：从压缩到伪逆，让这段 C 代码产生工程价值
6.1 低秩近似与矩阵压缩
拿到降序排列的奇异值之后，第一个最直接的应用就是低秩近似。对图像矩阵做奇异值分解，只保留前 k 个最大的奇异值，其余置零，再把三个矩阵乘回去，得到的矩阵就是原始矩阵的最优 k 秩近似。这个性质叫 Eckart-Young 定理，也是主成分分析和图像压缩的理论基石。
一个典型的调用流程是：读图像灰度矩阵 → 奇异值分解 → 置零小奇异值 → 重构矩阵。保留的奇异值数量决定了压缩比和失真程度。我见过最直观的实验是保留前 10% 的奇异值，重构图像视觉上几乎无损，文件体积却下降了 80% 以上。当然，实际工程里不会用这种纯 SVD 的方式做压缩（有 JPEG 等更高效的手段），但在理解数据能量分布时，这种实验非常有价值。
6.2 计算伪逆：让奇异值分解处理病态方程组
当你想解线性方程组Ax = b但 A 不可逆或者欠定时，可以用奇异值分解构造伪逆。公式是A^+ = V·S^{-1}·U^T，其中S^{-1}表示把每个非零奇异值取倒数，零奇异值保持为零。伪逆给出的解是“最小二乘意义下范数最小的解”，比普通高斯消元稳健得多。
使用场景包括：传感器标定里的过度约束问题、机器人运动学逆解、推荐系统里的矩阵补全。关键参数是“截断奇异值”——凡是小于阈值ε的奇异值直接忽略，避免病态分量放大噪声。这个阈值的选择直接影响解的稳定性，工程上常用的是ε = 1e-10 * σ_max，也就是最大奇异值的相对比例。
6.3 奇值分解与主成分分析的接口对接
如果你手里有某个做过零均值的样本矩阵，奇异值分解的右奇异向量 V 就是主成分方向，奇异值的平方正比于对应方向上的方差。按奇异值降序排列后的前几列 V，就是降维投影矩阵。
这个认知让奇异值分解从线性代数直接变成了数据分析工具。对于 C 代码落地，你只需要在分解得到 V 之后做一步列采样，就能立即输出降维结果，不需要额外写协方差矩阵的特征分解。这比单纯去调一个 PCA 库要快得多，也更容易和已有数据链路整合。
6.4 代码向量化与性能优化方向
当矩阵维度逐渐增大，纯标量版的雅可比迭代会越来越慢。优化的第一步是启用编译器自动向量化，比如在编译时加-O3和-march=native。第二步是改为连续内存布局并应用循环分块。第三步才是换成 Householder 双对角化。这三个层次逐步递进，每一步的改动量都不大，但收益明显。
另外，如果矩阵维度是固定的（比如姿态估计里恒为 3×3），完全可以硬编码展开循环，把迭代次数和旋转矩阵预先计算成常量表达式。这种特化能把单次分解时间从微秒级压到纳秒级，我在做惯性导航插件时就用了这个方案，效果非常明显。
做这套东西到现在，我的一个习惯是每个模块都自带一个自检函数，编译时加-DSELF_TEST开关就能跑内置用例。遇到奇异值不降序、重构误差超标或者收敛慢，先自检，再定位业务数据，这两步顺序颠倒会浪费时间。无论你是第一次接触奇异值分解，还是已经写过几版实现，记住一条：先把最小可验证用例跑通，再把算法塞进业务里替换旧的线性代数库。希望帮到你。
本文还有配套的精品资源，点击获取
