★ FEATURED ARTICLE 1. 动机从一张“糊掉”的照片到老老实实手写SVD入行头两年我在一次数据清洗任务里遇到了一个非常尴尬的问题手里是一个超高维的特征矩阵正交归一化做完之后还要跑一遍矩阵分解做数据压缩结果系统直接把内存打满了。同事随口丢了一句“用SVD截断一下”我盯着那三个字母发呆。后来我把奇异值分解SVD从头推导、从零实现、又在真实场景里完整用了一遍才明白为什么说SVD是机器学习里最不该跳过的基础设施。我现在在做的事情就是把这套完整的“奇异值分解SVD与实现”笔记整理出来。内容既包含手写幂迭代的核心思想也包括使用成熟库完成截断SVD、图像压缩、以及与主成分分析PCA的换算关系。适合刚学完线性代数觉得矩阵很抽象的人也适合想做推荐系统、降维、去噪、压缩但总在“调库”的人。看完以后你至少能做到三件事第一说出U、Σ、V这三个矩阵到底各自在干什么第二遇到千万级维度的矩阵时知道为什么库函数也会卡死第三亲手写一段能压缩图像的代码直观感受奇异值个数如何影响重建质量。很多人第一次接触SVD都是在一个公式里被吓退的A UΣVᵀ三个矩阵相乘而已。但真正让我“脑子通电”的是一个实际问题把图片变成矩阵再把它拆开发现绝大部分信息被集中在少数几个奇异值上丢掉后面几百个奇异值图片肉眼几乎看不出变化。那一刻我才意识到SVD不是考试题而是一台智能信息压缩机。2. 奇异值分解的核心直觉矩阵不是表格是线性变换2.1 矩阵拆解成三层的理解方式不用管复杂的推理我们先把SVD当作一次“拆解”。对任意一个形状为 m×n 的实矩阵 ASVD把它写成A(m×n) U(m×m) · Σ(m×n) · Vᵀ(n×n)但按工程习惯我们通常只用“瘦身版”U取前n列Σ取n×n对角阵Vᵀ取n×n。这样处理之后Σ从上到下依次放着奇异值它们都是非负的通常我们从大到小排序。怎么理解这三块我用一个不太严谨但很好用的生活类比U是“输出方向”的坐标系V是“输入方向”的坐标系Σ是连接两个坐标系的“缩放力度”。矩阵乘向量时系统是先借用V那边换到主轴再用Σ逐轴拉伸/压缩最后用U换回原坐标。所以SVD的本质是在揭示“一个线性变换最核心的几个拉伸方向”。2.2 奇异值大小直接告诉你信息藏在哪奇异值不是均匀出现的。真实数据里的矩阵前几个奇异值往往非常大后面的迅速衰减。这个现象背后的意义极其重要它说明矩阵的有效秩远远低于它的实际行数或列数。也就是说一张1024×1024的图名义上有一百万个数字但真正承载视觉轮廓的可能只有前50个奇异值。剩下的数字更像是噪声、纹理细节或完全可忽略的冗余。机器学习里做降维本质上就是保留这前k个奇异值丢掉尾巴。2.3 低秩近似公式千万不要强行逐元素复制我将每次实际操作中都会用到的低秩近似公式写在下面A ≈ U_k · Σ_k · V_kᵀ这里的下标k表示只取前k列、前k个奇异值和前k行。这个近似不是简单把矩阵里某些“看起来不重要的元素”设为零而是从全局角度寻找一个低秩矩阵使它和原始矩阵的差异尽可能小。数学上有理有据现实中效果也惊人图像压缩、协同过滤、主题模型全都建立在这一条式子上。初次上手的人容易产生一个误区认为SVD得到的U、V正交向量没什么区别。但实际代码里U是列正交的V是行正交的而Σ排列顺序直接影响你截断时的取舍。调试代码时我经常要先打印出s[:10]看一下奇异值的分布再决定保留个数。直接拍脑袋定k十有八九会出错。3. 动手前的准备环境、工具选型与整体设计思路3.1 依赖选型别一上来就只装库我做SVD实践时minimal环境一般是这三样Python 3.10、NumPy、SciPy。NumPy里有个linalg.svdSciPy里有个sparse.linalg.svds这两个接口我反复用但它们的定位完全不同。如果你打算只装NumPy那么对于形状为1000×800的矩阵跑完整SVD完全没问题如果矩阵到了20000×15000完整SVD的时间成本和内存占用就会很肉疼此时应该用SciPy的截断版本或随机化SVD。3.2 功能模块拆分从输入到可视化一次跑通我不会把一个SVD项目写成一大坨脚本。以“图像压缩”这个目标为例我把整个实现拆成了五块模块职责核心函数建议数据读取模块读取图片、转灰度、归一化load_gray_img(path)分解模块得到U、s、Vtnp.linalg.svd(...)截断重建模块按指定k重建矩阵reconstruct(U_k, s_k, Vt_k)评估模块计算相对误差、峰值信噪比PSNRcal_error(...)可视化模块对比不同k下的重建结果matplotlib拆完之后每个模块都能单独测试调试时不会把所有问题混在一起。更重要的是这样的结构方便你替换“分解模块”今天用NumPy明天换成SciPy或手写版本后面章节我会展示完整的多版本实现方案。3.3 数据准备的三个检查点确认矩阵元素不是整数int8如果是图像数据读进来后务必用.astype(np.float64)转换否则SVD计算过程中的除法会累积误差甚至溢出。确认矩阵中没有NaN值SVD对缺失值极其敏感有一两个NaN就能导致结果整体无法使用。在推荐系统场景里缺失评分通常会先用均值填充再分解。确认矩阵不是“行秩严重不足”如果矩阵里大量行完全相同奇异值会出现大量精确的0截断后解析容易但不要误以为分解失败。4. 核心实现从零手写SVD到调用成熟库函数4.1 手写幂迭代一口气挤出最大奇异值网上很多“手写SVD”教程直接照抄数值计算书籍里的双重雅可比算法代码巨长反而让人看不懂。我选了一个真正能讲清楚原理的方法幂迭代。它的思路很朴素随机初始化一个向量反复应用矩阵A及其转置Aᵀ让向量在主奇异向量的方向上不断膨胀最后归一化得到主奇异向量。下面是我自己整理的精简版本只负责求第一个奇异值和对应向量import numpy as np def power_iteration(A, num_iters500, tol1e-8): m, n A.shape v np.random.default_rng(42).standard_normal(n) v / np.linalg.norm(v) sigma 0.0 for _ in range(num_iters): u A v new_sigma np.linalg.norm(u) if new_sigma 1e-12: break u u / new_sigma v_new A.T u v_new_norm np.linalg.norm(v_new) if v_new_norm 1e-12: break v_new v_new / v_new_norm # 由于奇异向量只有方向有意义正负方向都需要判停 if np.linalg.norm(v_new - v) tol or np.linalg.norm(v_new v) tol: v v_new sigma new_sigma break sigma new_sigma v v_new return sigma, v这段代码最容易被忽略的地方有两个回归归一化。每次把它拆掉奇异值就是本来要更新的。拆完后如果直接除不会稳定但如果你先除n个再除n个还是会因为浮点误差导致震荡。第二个是判停条件里必须考虑符号因为奇异向量v和-v代表同一个方向不能因为符号翻转就判为不收敛。4.2 补全多个奇异值幂迭代加压缩收缩数据一维的情况谁都会实际工程里我们往往要前k个奇异值。经典做法是算出最大奇异值和向量后把这部分贡献从矩阵中扣除再对剩余矩阵继续调用幂迭代。这种压缩法在代码里只要一行A_work A_work - sigma * np.outer(u, v)但请注意这个操作会引入累积误差k越大后面的奇异值误差越大。我自己在小矩阵上测试时k取10以内效果还能接受再往后误差就开始肉眼可见地增大。所以工程上一旦要取几十个奇异值我的建议是直接交给np.linalg.svd手写版本用于理解原理和教学演示就够了。4.3 复现幂迭代过程中最常见的坑我特别想强调一下初始向量的影响。幂迭代收敛速度取决于第二大奇异值和最大奇异值的比值σ₂/σ₁。这个比值越接近1收敛越慢。遇到两个奇异值几乎相等的情况比如矩阵的两条主轴强度一致幂迭代会有很大概率收敛到其中一个而不是稳定输出完整解集。解决方法是多尝试几个不同的随机种子或者改用带正交化处理的迭代方法。真正常见的场景不是等我拍脑袋而是数据量极大时矩阵的差异其实被压缩得很小往往收敛很慢此时我会直接上随机化算法。4.4 成熟库函数的正确打开方式手写版本理解原理之后我们的决胜武器还是标准库。最常用的写法有两种一种我总结到下表需求推荐函数关键点中小矩阵完整分解np.linalg.svd(A, full_matricesFalse)输出U, s, Vts为一维数组不是对角阵超大稀疏矩阵只取前k个奇异值scipy.sparse.linalg.svds(A, k10)返回U, s, Vt但奇异值默认升序排序需要稳定符号方向手动统一符号比如让U每列第一个元素为正SVD结果不是唯一的正负号可以整体翻转我最常犯的错误是把svds的结果直接当成降序。SciPy里svds返回的奇异值是升序的记住了否则你去截断时会把最小的几个留下来。这是一个非常容易复现的坑我在一次项目里踩完后几乎每次查看s[0]之前都会强制np.sort(s)[::-1]。5. 实战用SVD压缩图像亲手感受“有效低秩”5.1 为什么拿图像开刀图像是人类知觉系统最容易感知“信息损失”的对象。把RGB图片转成灰度矩阵以后你会看到一个数值范围在0到255之间的二维矩阵SVD的每个步骤都直观可见。图像压缩还能直接看出截断维度k与重建画面质量的关系比推荐系统里那些看不见的评价指标更有说服力。5.2 完整图像压缩实现我和你们分享一段我在实验里反复使用的脚本。我这里用一张名称为demo_photo.jpg的灰度图做演示。图片不是关键你也可以用任意一张自己的照片。import numpy as np import matplotlib.pyplot as plt import matplotlib.image as mpimg def load_gray_img(path): img mpimg.imread(path).astype(np.float64) if img.ndim 3: # 处理RGB或RGBA图 img img[..., :3].mean(axis-1) return img def compress_svd(img, k): U, s, Vt np.linalg.svd(img, full_matricesFalse) U_k U[:, :k] s_k s[:k] Vt_k Vt[:k, :] return U_k np.diag(s_k) Vt_k, s if __name__ __main__: A load_gray_img(demo_photo.jpg) print(原始矩阵形状, A.shape) for k in [5, 20, 60, 120]: A_rebuild, s compress_svd(A, k) # 计算相对误差 error np.linalg.norm(A - A_rebuild) / np.linalg.norm(A) print(fk{k:3d}, 相对误差{error:.4f}) # 画对比图 ks [5, 20, 60, 120] fig, axes plt.subplots(1, 4, figsize(16, 5)) for ax, k in zip(axes, ks): A_rebuild, _ compress_svd(A, k) ax.imshow(A_rebuild, cmapgray) ax.set_title(fk{k}) ax.axis(off) plt.show()运行这段代码你会看到一组非常直观的结果k5时图基本只有“一团雾”k20能看出轮廓k60已经接近原图k120肉眼几乎分辨不出区别。而事实上原始矩阵可能有几百甚至上千的奇异值。5.3 奇异值能量占比的参数选取方法“k到底选多少”这是新手最常见的疑问。我的做法是计算“前k个奇异值的能量占比”近似公式用奇异值的平方和来算def energy_ratio(s, k): return np.sum(s[:k] ** 2) / np.sum(s ** 2) # 例子找到达到90%能量的最小k s np.linalg.svd(A, compute_uvFalse) ratio_cumsum np.cumsum(s ** 2) / np.sum(s ** 2) r np.argmax(ratio_cumsum 0.9) 1 print(保留90%能量所需的k, r)能量占比达到0.9到0.95通常就是一个性价比很高的压缩点。但不要迷信固定的阈值。有一次我用一个人像照片测试90%能量对应的k非常小人脸轮廓依然清晰换到一张布满细密纹理的叶子照片后同样的90%能量重建出来的图就发糊。所以场景才是最终决定k的主体。想让肉眼满意直接固定几个k值对比着看。5.4 存储与内存优化方案别直接重建大矩阵如果你是为了“压缩存储”完整重建U_k np.diag(s_k) Vt_k只适合可视化。真正部署时应该直接保存U_k、s_k、Vt_k三个小矩阵。假如原图是1024×1024k60原本需要存储一百万个浮点数现在只存储1024×60 60 60×1024也就是大约12万个浮点数压缩了将近八成。这个思想延伸到机器学习里的模型压缩和数据降维都是同一个套路。我通常在压缩前把图片转成float16保存成npz文件内存和磁盘同时收益。6. 常见问题与排障技巧实录6.1 SVD计算告警“不收敛”怎么办当使用scipy.sparse.linalg.svds时它本质上是迭代算法如果矩阵维度过大且分布非常怪异可能出现不收敛或收敛到局部的情况。我的第一反应是增加迭代次数并设置tol更大一些比如maxiter2000。如果还是不收敛就要检查矩阵是否包含全零行或全零列先去掉这些数据再去跑分解往往立刻见效。6.2 奇异向量符号随机翻转np.linalg.svd每次返回的U、V列向量本身不是唯一确定的。如果矩阵存在近似重复的奇异值向量符号会在不同运行中随机变化。这在做增量更新或需要“对齐”模型时会造成困扰。我的解决方式很粗暴在得到U之后让每一列第一个绝对值最大的元素强制为正。这个规则你写在数据处理函数里后面所有下游模块都能保持统一。6.3 完整SVD和截断SVD的内存量差异很多新人在处理百万行矩阵时兴致勃勃写U, s, Vt np.linalg.svd(A)然后内存爆掉。核心原因是完整SVD会尝试返回U和V这两个可能非常大且稠密的矩阵。此时务必改成from scipy.sparse.linalg import svds # 只要前20个奇异值 U, s, Vt svds(A, k20)但我们仍然必须记得svds返回的奇异值默认升序所以使用前最好排序。如果矩阵是稀疏且稀疏度很高可以先用scipy.sparse.csr_matrix格式存储再传给svds内存消耗可以显著降低。6.4 SVD和PCA的等价与迁移最后一条避坑知识协方差矩阵的特征分解和居中后的SVD有明确的对应关系。我们在做PCA时先对特征矩阵X按列中心化再对X做SVD得到的右奇异向量V就是PCA的主方向。而左奇异向量U与标准化后的主成分得分成正比。不少机器学习框架在实现PCA时底层根本不做协方差矩阵求逆而是直接调SVD。所以你学会SVD等于顺手把PCA的原理也掌握了一大半。但要注意SVD本身不要求中心化PCA必须先中心化直接对原始矩阵做SVD结果不等同于对原始特征做PCA。6.5 奇异值骤降为0怎么排查真实业务矩阵里常出现完全重复的行/列或者多重共线性字段。此时奇异值会出现一个或多个精确的0。这种情况并不代表SVD失败反而是矩阵秩不足的信号。我会用它来诊断数据质量奇异值中出现大量0说明数据冗余度极高可以先做特征去重再进入模型训练。7. 踩坑之后我对SVD的看法如果只让我说一条个人体会一定要在真实数据上“看”一次奇异值分布。你盯着代码里的s数组看五分钟比看十篇SVD推导更管用。我尝到的甜头之一是在一次多模态特征融合项目中直接用奇异值能量占比找出了一批冗余特征把维度从两千多压到200下游模型训练时间肉眼可见缩短。还有一个小技巧也是我现在每回写SVD相关代码都会做的把SVD过程封装成一个统一函数内部自动处理符号一致性、奇异值排序、非零值过滤接口只暴露矩阵和k。别人接手你的代码时就不用再踩scipy.sparse.linalg.svds升序那个坑了。最后建议你今晚就找一张自己的照片跑一遍上文第5节的脚本分别用k5、20、60、100看一遍重建效果。亲眼看到奇异值排名越靠前的信息越“值钱”你对SVD的理解才算真正落地生根。
阅读完成 · 觉得有帮助？
