从零手写SVD：嵌入式C语言实现奇异值分解的完整实践

简介：SVD（奇异值分解）的C语言实现源码包，主要面向数值算法初学者、嵌入式开发者以及需要做矩阵运算的科研人员，帮助理解如何在C/C++环境下高效完成矩阵分解，并为PCA、协同过滤、图像压缩与降噪、文本分析等经典应用提供底层算法参考。资源包共80个文件、1.12MB，核心源码为svd.cpp、main.cpp与svd.h，完整覆盖从特征值/特征向量计算、奇异值排序到构造Σ、U、V矩阵的SVD主流程；同时附带Visual Studio工程配置（sln、vcxproj、dsp）、可执行exe、编译日志（tlog/log）和调试中间文件，方便直接打开、运行与单步跟踪。已有1087人学习。通过阅读代码可掌握数值稳定性处理与迭代算法（如Golub-Kahan、Lanczos）的具体实现思路，也可将分解结果与LAPACK等标准库对照验证，是学习矩阵分解和推荐系统底层原理的实用资料。 上个月在嵌入式端做一个多通道传感器数据预处理模块，需要反复做矩阵的奇异值分解（SVD），用来做PCA白化和数据去相关。板子是ARM架构的精简Linux，LAPACK装不上，OpenCV太重，GSL交叉编译后库体积也超标，最后决定直接用C语言从零实现SVD。这篇博文会完整记录我手写SVD的算法选型、核心代码、验证方式和踩坑经历，适合嵌入式开发者、算法移植工程师，以及想在C环境下把矩阵分解真正搞明白的同学参考。

1. 先说结论：三种实现路径的取舍，以及我为什么走到手写这一步

1.1 三条路的对比

很多朋友第一反应是：SVD不是有现成库吗？确实有，而且性能比我手写的版本好得多。但在实际工程里，选择哪条路取决于你的部署环境和约束条件，我把三种常见方案放在一起对比：

我在项目里最开始选的是GSL。交叉编译链都配好了，结果把库放进目标板一看，Flash和RAM的占用直接涨了一截，仅仅为了一个8×8矩阵的SVD，却要把整个GSL背在身上，这种资源浪费在工程上是说不过去的。后来换LAPACK，又发现它和BLAS的版本绑定关系复杂，交叉编译时各种符号缺失，折腾了一下午，最后还是放弃了。

1.2 手写之前必须想清楚的几个问题

• 你的矩阵是方阵还是非方阵？SVD算法需要明确支持m×n的一般矩阵，还是只要方阵就行。单边Jacobi法对m≥n的一般矩阵直接可算，但如果m<n，需要先转置处理。
• 精度要求是多高？工控场景里float往往不够，double是底线。如果你的应用里矩阵是病态的，还要考虑是否需要用更高精度的累加或预处理。
• 后续是否要长期维护？如果只是验证算法可行性，手写完全够用；如果是大型项目几年内持续演进，用成熟库能省掉不少维护成本。

就算最终打算调用现成库，我也不建议跳过这一步分析。搞清楚“我为什么不用LAPACK”，比“我会调用LAPACK”更能帮你形成可靠的技术判断力。

2. 从几何直觉到数学定义：SVD到底在做什么

2.1 矩阵也是一种变换

SVD说的是：任意m×n实矩阵A，都能被分解成A = UΣV^T。其中U和V都是正交矩阵，Σ是“对角矩阵”。从几何角度理解，任何线性变换都可以被拆成三步：先旋转一下，再沿坐标轴拉伸或压缩，最后再旋转一下。这个直觉很重要，因为后面所有SVD的应用——PCA、伪逆、最小二乘——本质上都是在利用这个“旋转-拉伸-旋转”的分解。

举个例子。一个摄像头把三维世界投影到二维平面，相机内参矩阵的SVD就能分解出等效焦距的方向和尺度。在传感器数据处理里，多通道信号之间的相关性也能通过SVD找到“主要方向”和“主要强度”。

2.2 U、Σ、V每一项的工程含义

• V的列：是A^T A的特征向量，代表数据在原始空间中的“主方向”。
• Σ的对角元：奇异值σ1 ≥ σ2 ≥ … ≥ σr ≥ 0，按降序排列时，前几个奇异值集中了绝大部分“能量”。
• U的列：是数据在这些主方向上的投影坐标，各列之间彼此正交。

对工程师来说，最直接的用处是PCA。很多人用协方差矩阵的特征分解做PCA，但换成SVD更稳妥，因为特征分解会先计算A^T A，数值上会把条件数平方，导致小奇异值误差被放大。SVD直接对A操作，数值稳定性更好。

2.3 SVD的应用远不止PCA

推荐系统里的潜在语义分析（LSA）也是SVD的典型应用：把词项-文档矩阵做奇异值分解，去掉小奇异值对应的分量，就相当于把高维稀疏的文本向量压缩成低维稠密的语义向量。这个思想后来演变成了主题模型里很多方法的数学基础。还有最小二乘里求伪逆，统一公式是A⁺ = VΣ⁺U^T，有了SVD就能一并搞定。

3. 单边Jacobi旋转法：代码量最小且精度最高的选择

3.1 为什么不用Golub-Kahan算法

Golub-Kahan是LAPACK dgesvd_背后的经典算法：先把矩阵双对角化，再用隐式QR迭代求奇异值。这个算法性能极佳，但实现细节非常多。双对角化阶段要处理Householder变换的符号选择，隐式QR迭代阶段要设计位移策略和收敛判据。我评估了一下，完整手写一遍至少要上千行，调试周期很长，而且稍不留神就会在某些矩阵上不收敛。

相比之下，单边Jacobi法的实现简练很多，核心操作只有一个：选两列，做个旋转，让它们正交。这个操作本身不复杂，每一步都看得见摸得着，调试起来非常直观。对于中小规模的矩阵（比如n ≤ 20），精度甚至可以超过Golub-Kahan，因为它本质上是在做正交相似变换，不会引入额外的数值损失。

3.2 单边Jacobi的数学推导

算法的目标是：对A右乘一系列旋转矩阵J，得到B = A·J，让B的列两两正交。当B的列两两正交时，可以写B = UΣ，其中U是B各列归一化后拼成的矩阵，Σ是对角线上放各列二范数的对角阵。因为J是正交矩阵的乘积，所以V = J仍然是正交矩阵，于是A = B·J^T = UΣV^T。SVD就得到了。

问题变成：怎么选旋转矩阵，让一对列正交？对列(p, q)，令：

• alpha = ‖a_p‖²（第p列的范数平方）
• beta = ‖a_q‖²（第q列的范数平方）
• gamma = a_p·a_q（两列内积）

我们要找旋转角θ，使得旋转后的第p列和第q列内积为零。旋转公式是：

• a_p' = c · a_p - s · a_q
• a_q' = s · a_p + c · a_q

其中c = cos θ，s = sin θ。代入内积为零的条件，可以解出θ。但在实际代码里，我不会直接用atan2去算θ，而是用下面的数值稳定公式。

3.3 数值稳定的旋转角计算

直接调用atan2(2·gamma, alpha - beta)再算cos、sin，代码看上去简单，但某些角度下精度损失较大，而且多次调用三角函数开销不小。更常用的做法是先算一个中间变量zeta：

• zeta = (beta - alpha) / (2 · gamma)
• t = sign(zeta) / (|zeta| + sqrt(1 + zeta²))
• c = 1 / sqrt(1 + t²)
• s = c · t

这个公式来自解二次方程时取较小根的技巧，可以避免zeta趋于无穷时带来的数值精度问题。第一次看到这个写法可能觉得绕，但它在double精度下非常稳，我建议直接沿用。

4. C语言实现：完整代码与关键细节

4.1 内存布局和接口约定

这个实现采用行主序存储，接口传裸指针double*，不额外封装结构体。好处是零额外内存开销，调用方在嵌入式环境下可以自由选择静态数组或内存池。函数只支持m ≥ n，如果m < n，直接返回-1。注意这是thin SVD：U是m×n，不是m×m。对大多数应用（PCA、伪逆、最小二乘）来说thin SVD已经够了。

4.2 完整C代码

#include <stdio.h> #include <stdlib.h> #include <math.h> #include <string.h> #define SVD_EPS 1e-15 #define SVD_MAX_SWEEP 100 // Compute inner product of column p and column q in mat (row-major m x n) static double col_dot(const double *mat, int rows, int cols, int p, int q) { double sum = 0.0; for (int i = 0; i < rows; i++) { sum += mat[i * cols + p] * mat[i * cols + q]; } return sum; } // Compute norm of column j static double col_norm(const double *mat, int rows, int cols, int j) { return sqrt(col_dot(mat, rows, cols, j, j)); } // Single-sided Jacobi SVD. // A: m x n row-major matrix, m >= n // Output: // U: m x n, columns are left singular vectors // S: n, singular values in descending order // V: n x n, rows? no, V is n x n row-major, stored as V^T semantics // The caller can reconstruct A as U * diag(S) * V^T int svd_jacobi(const double *A, int m, int n, double *U, double *S, double *V) { if (m < n) { return -1; // only support m >= n } // U = A, V = I memcpy(U, A, m * n * sizeof(double)); memset(V, 0, n * n * sizeof(double)); for (int i = 0; i < n; i++) { V[i * n + i] = 1.0; } // Jacobi sweep: make columns of U orthogonal for (int sweep = 0; sweep < SVD_MAX_SWEEP; sweep++) { double max_off = 0.0; for (int p = 0; p < n - 1; p++) { for (int q = p + 1; q < n; q++) { double alpha = col_dot(U, m, n, p, p); double beta = col_dot(U, m, n, q, q); double gamma = col_dot(U, m, n, p, q); if (alpha == 0.0 || beta == 0.0) { continue; // zero column, nothing to rotate } double off = fabs(gamma) / sqrt(alpha * beta); if (off > max_off) { max_off = off; } if (off < SVD_EPS) { continue; // already orthogonal enough } // Numerically stable Jacobi rotation angle double zeta = (beta - alpha) / (2.0 * gamma); double t = (zeta >= 0.0 ? 1.0 : -1.0) / (fabs(zeta) + sqrt(1.0 + zeta * zeta)); double c = 1.0 / sqrt(1.0 + t * t); double s = c * t; // Rotate columns p, q of U for (int i = 0; i < m; i++) { double up = U[i * n + p]; double uq = U[i * n + q]; U[i * n + p] = c * up - s * uq; U[i * n + q] = s * up + c * uq; } // Rotate columns p, q of V for (int i = 0; i < n; i++) { double vp = V[i * n + p]; double vq = V[i * n + q]; V[i * n + p] = c * vp - s * vq; V[i * n + q] = s * vp + c * vq; } } } if (max_off < SVD_EPS) { break; } } // Extract singular values, normalize U columns for (int j = 0; j < n; j++) { double norm = col_norm(U, m, n, j); S[j] = norm; if (norm > 1e-300) { for (int i = 0; i < m; i++) { U[i * n + j] /= norm; } } } // Sort singular values descending, sync columns of U and V 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) { continue; } double tmp = S[i]; S[i] = S[max_idx]; S[max_idx] = tmp; for (int r = 0; r < m; r++) { tmp = U[r * n + i]; U[r * n + i] = U[r * n + max_idx]; U[r * n + max_idx] = tmp; } for (int r = 0; r < n; r++) { tmp = V[r * n + i]; V[r * n + i] = V[r * n + max_idx]; V[r * n + max_idx] = tmp; } } return 0; }

4.3 代码里几个容易忽略的细节

这个代码看起来不算长，但有几个细节是调试时容易踩坑的地方。

第一，V的初始化必须是单位矩阵。后续每次旋转都会累积在V上，V最后就是总的右乘矩阵J。如果初始化为零矩阵，后面怎么转都白搭。

第二，为什么用max_off判断收敛？因为一轮扫描中可能会有多个列对需要旋转，只要其中有一对还不够正交，就得继续下一轮。我用网格扫描的方式遍历所有列对，一轮结束后统计最大的非正交度量，小于阈值就终止。

第三，排序动作必须在最后统一做。奇异值降序排列是SVD的通用约定，排序时U的列、V的列、S的值三处必须同步交换，漏掉任何一处，重构出来的矩阵就会错。

5. 实测与排坑：从“能算”到“算得对”

5.1 验证正确性的黄金标准：重构误差

代码写完，第一步不是看奇异值是否和Python算出的结果一致，而是直接做重构。对随机矩阵A，用分解出的U、S、V重新计算A_rec = U·diag(S)·V^T，然后看最大元素误差。如果误差在1e-14量级，基本可以确定实现正确。

下面这段测试代码我放在了main函数里：

int main(void) { const int m = 5, n = 3; double A[m * n] = { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 2.0, 3.0, 1.0, 5.0, 4.0, 2.0 }; double U[m * n], V[n * n], S[n]; if (svd_jacobi(A, m, n, U, S, V) != 0) { fprintf(stderr, "SVD failed: m must be >= n\n"); return 1; } printf("singular values:\n"); for (int i = 0; i < n; i++) { printf(" s[%d] = %.15e\n", i, S[i]); } double max_err = 0.0; for (int i = 0; i < m; i++) { for (int j = 0; j < n; j++) { double rec = 0.0; for (int k = 0; k < n; k++) { // V[j][k] is V^T[k][j] rec += U[i * n + k] * S[k] * V[j * n + k]; } double err = fabs(A[i * n + j] - rec); if (err > max_err) { max_err = err; } } } printf("max reconstruction error = %e\n", max_err); return 0; }

随机测试矩阵是验证算法正确性的最佳选择，因为随机矩阵几乎必然是列满秩的，不会出现零列或线性相关这种“幸运”情况。如果随机矩阵都能通过重构检验，那实现正确性的可信度就很高了。

5.2 病态矩阵、秩亏矩阵和边界情况

随机矩阵测试通过之后，我开始加边界情况。最典型的是Hilbert矩阵（H_ij = 1/(i+j+1)），它的条件数随着维数指数增长，是公认的病态矩阵。用6×6的Hilbert矩阵实测，Jacobi法依然能保持不错的重构精度，但奇异值很小的那几个分量，误差会稍微大一点。这是数值上不可避免的，不必太焦虑。

更需要注意的是秩亏矩阵。比如某列本身就是零向量，那么计算alpha或beta时会得到0.0，直接continue。真正麻烦的是排序完成后，零奇异值对应的列不能参与归一化，否则0/0会得到NaN。代码里用if (norm > 1e-300)做了保护，这个阈值不是随手写的，而是double能表示的最小非零规格化数的量级附近。

5.3 一次真实排坑记录

我调试时遇到过一个诡异现象：sweep=20时，大部分矩阵的重构误差已经到1e-14，但个别矩阵的误差始终停在1e-10，怎么都降不下来。定位后发现，问题出在阈值上。我一开始把SVD_EPS设成1e-12，想省几轮迭代，结果列对没有充分正交，旋转就提前终止了。把SVD_EPS改回1e-15后，同样的矩阵误差立刻降到1e-14量级。

这说明两个问题：第一，阈值必须和double精度匹配，不是越大越好；第二，迭代轮数不够时，即使max_off已经很小，也不代表所有列对都严格正交了。后来我在工程版里加了一个保护：当sweep达到上限但max_off仍然较大时，返回-2提示迭代不收敛。

另一个坑是m < n的情况。我最初的实现没有检查维度，直接对A做SVD，结果拿到的U和V维度怎么都对不上。后来明确约定：函数只支持m ≥ n，m < n时由调用方对A^T求SVD，再把U和V交换回来。只要在注释里写清楚，这个限制完全不影响使用。

6. 工程化要点：API设计、内存策略和后续方向

6.1 对外接口应该怎么设计

实际项目里，我不建议让调用方直接面对这个几百行的实现。更合理的做法是封装成一个稳定的接口，返回int错误码，0表示正常，-1表示维度错误，-2表示迭代不收敛。内存方面，U、S、V都由调用方预留，库内部不使用malloc。这一点在嵌入式环境里几乎是必须的，因为目标板上可能没有标准堆，或者需要统一走内存池。

接口设计还有一个细节：SVD的输出符号不是唯一的。同一个矩阵，把U的某一列和V的对应列同时取反，重构结果完全不变。如果你的上层算法对符号有严格要求，需要在封装层统一符号约定，比如约定U的第一行非零元素为正。

6.2 嵌入式场景下的优化思路

如果后续性能不够，可以考虑这几个方向：

• 单精度float版本。很多工控场景不需要double级别的精度，float版本体积更小、速度更快，但SVD_EPS要放宽到1e-6左右。
• 固定迭代次数。对n ≤ 8的小矩阵，可以固定sweep=10次，不必每次都判断max_off，减少分支开销。
• 并行化。Jacobi每轮扫描里，多个列对（p,q）的旋转在数学上是近似独立的，可以在多核平台用OpenMP并行，但要注意同一个列不能同时被两个线程更新。

我在目前的板子上测过，8×8矩阵double精度下，100次SVD总耗时不到2毫秒，对于数据预处理管线的要求完全够用。

6.3 最后一点个人体会

写SVD的整个过程比直接调库多花了一周时间，但收益也是明显的：我现在对“矩阵分解”这件事的理解，远比只会调包时深入。如果项目周期紧张，直接用成熟库才是正解；但如果是做算法移植、嵌入式调优，或者单纯想把线代基础打牢，手写一遍的收获是不可替代的。这个实现已经在我板子上稳定跑了一个多月，后续我准备在此基础上封装伪逆和最小二乘求解器，等做完了再继续分享。

本文还有配套的精品资源，点击获取