机器学习线性代数Python实战：矩阵运算、SVD与PCA降维

简介：面向机器学习入门者与高校课程学习者的线性代数配套代码包，围绕张雨萌《机器学习线性代数基础（Python语言描述）》组织，适合阅读教材时同步动手验证抽象概念，也可供授课教师选用练习素材。压缩包共53个文件、约2.48MB，主体是50个py示例脚本，按第1至第7章分目录编号，覆盖Numpy科学计算、矩阵运算、线性方程组与特征值分解等章节内容；另有1份docx库说明、1个csv数据文件及1张png测试图，用于环境入门、数据读取和结果对照。已有187人学习下载，具备一定参考价值。读者可按章节运行代码、对照书中推导复现实验，理解关键算法的实现思路，并修改参数扩展练习，提升用代码表达线性代数运算与机器学习预处理的能力；清晰编号也便于查漏补缺，适合课程实验与自学复盘。

1. 从一次矩阵维度对不上说起

前阵子帮一个做推荐系统的朋友看代码，他卡在一个特别典型的场景：特征矩阵X是(10000, 128)，想用 PCA 压到 32 维，sklearn一行就出来了，但被问“为什么是 32”“这 32 维到底保留了多少信息”的时候答不上来。再往前翻，从梯度下降的矩阵求导到 SVD 的低秩近似，全是靠调包绕过去的。这份「机器学习线性代数基础（Python语言描述）」配套代码，恰好就是填这个坑的东西——它不摆定理，把向量、矩阵、行列式、特征分解、SVD、PCA 一条线用 Python 从零写一遍，每段代码都能单独跑。适合两类人：刚上手机器学习被numpy形状和轴搞晕的新手，以及写了几年业务代码、想回头把数学直觉补上的工程老手。

2. 拆包与目录结构：从解压到跑通第一个脚本

2.1 目录怎么分层，先搞清楚每个文件对应哪一知识点

教材配套代码有个通病：作者按书里章节顺序排目录，读者按知识点回查。这类包常见的组织方式是按章编号分目录，每个目录里放.py脚本或.ipynb笔记本，文件名基本就是知识点名字。我在本地整理成下面这样，好处是后面改代码、写对照实验时不用再翻书目录。

实际拿到的压缩包如果编号不一致，不用纠结，按文件名里的关键词认就行，vector、matrix、eigen、svd、pca这几个词一定在。真正的坑在于：很多包把同一份功能既提供了手写实现，又提供了numpy一行版，命名上只差一个后缀。建议先把每个目录下的文件列表拉出来看一眼，搞清楚哪几个是“教学版手写”，哪几个是“标准库对照版”，后面调试的时候才不会被两个结果打架搞懵。

# 只看结构，不解压全部内容，先摸清有多少个知识点目录 unzip -l machine-learning-linear-algebra-python-code.rar | head -50 # 解压到指定目录，避免中文路径带来的编码问题 unzip machine-learning-linear-algebra-python-code.rar -d ./la-code cd ./la-code && find . -maxdepth 2 -name "*.py" | sort

第一条命令只列清单，用来看压缩包内层级是否多套了一层目录（很常见，解压后变成包名/包名/...）。第二条解压时指定英文目录，是因为部分脚本在读写图片、CSV 时用的是相对路径，中文路径在个别环境上会因为编码问题直接抛FileNotFoundError。find那行是给自己建索引，心里有个数：这个包到底覆盖了多少个知识点。

2.2 环境依赖：Python 版本与 numpy / matplotlib 的最小组合

这套代码的核心依赖非常薄，主体是numpy，画图部分用matplotlib，SVD 图像那节可能还要Pillow读图。没有深度学习框架，没有 GPU 依赖，一台普通笔记本足够。麻烦在于版本：书里的代码大概率是按numpy1.x 写的，而 2.x 删掉了一批老的类型别名（np.float、np.int、np.float_之类），直接跑会一片AttributeError。

# 建一个独立环境，别污染主环境 python -m venv .venv source .venv/bin/activate # Windows 用 .venv\Scripts\activate # 保守起步：先装 1.x 稳定线，跑通再考虑升级 pip install "numpy>=1.24,<2.0" matplotlib pillow jupyter # 装完自检，确认版本和关键函数都在 python -c "import numpy as np; print(np.__version__); print(np.linalg.svd(np.eye(3))[1])"

第一段的venv是硬性建议，不是洁癖。线性代数代码里大量用float64，如果你的主环境里混了某个要求numpy2.x 的包，回退版本会把别的项目搞崩，血泪经验。第二段把numpy锁在 1.x，是为了让书里的老写法原样能跑；等你把代码全过了一遍、心里有数了，再统一升级到 2.x，那时候哪一行报错你都知道该改成什么。最后那段自检跑通，说明numpy的线性代数子模块没问题，可以往下走了。

提示：如果你的机器上已经装了 2.x 又不想降级，把脚本里所有np.float换成float、np.int换成int、np.bool换成bool，绝大多数报错就消失了。这属于机械替换，用编辑器的全局搜索一次搞定。

2.3 两种运行方式：脚本适合看结果，notebook 适合看过程

这个包里的代码基本是双形态：.py脚本和.ipynb笔记本。脚本适合“我要一个确定的结果”，notebook 适合“我要一步一步看中间变量长什么样”。学线性代数这件事上，notebook 的优势非常明显——你可以把矩阵A单独跑一行，看看它的形状、秩、特征值，再决定下一步怎么改。

# quick_check.py —— 放在包根目录，用来验证环境和包内模块都能导入 import numpy as np import matplotlib matplotlib.use("Agg") # 无界面环境必须，否则 plt.show() 会卡住 import matplotlib.pyplot as plt A = np.array([[2.0, 1.0], [1.0, 3.0]]) print("shape :", A.shape) print("rank :", np.linalg.matrix_rank(A)) print("eigval:", np.linalg.eigvalsh(A)) # 对称矩阵用 eigvalsh，结果必为实数 plt.imshow(A, cmap="viridis") plt.colorbar() plt.savefig("check.png", dpi=120) print("figure saved -> check.png")

这段代码有三个作用。matplotlib.use("Agg")把后端切到非交互模式，在服务器、容器或者 SSH 环境里跑脚本时，不加这一行经常出现“程序不报错但也永远不结束”的玄学现象。matrix_rank用来确认矩阵是否满秩，后面解线性方程组、求逆之前先看一眼秩，能省掉大量LinAlgError。eigvalsh的h表示 Hermitian，对实对称矩阵来说比通用的eigvals更快也更稳，返回的一定是实数且已排序——这一点在对照书上的特征值结果时特别重要，因为eig返回的顺序是任意的，直接和书上比会以为算错了。

跑通这个脚本，环境这一关就算过了。接下来才是真正值钱的部分：把每个知识点的手写实现拆开看。

3. 核心模块拆解：向量、矩阵与手写运算层

3.1 向量运算：点积、范数与余弦相似度为什么总被写错

向量这一层看着最简单，实际上是翻车最集中的地方。numpy里*是逐元素乘、@是矩阵乘、np.dot在一维数组上是内积在二维数组上是矩阵乘，三种语义混在一个符号体系里，新手极容易写出一行不报错但结果错的代码。配套代码在向量这节通常会给出“手写版 + 标准库版”两套，对照着看能把这层语义彻底理清。

import numpy as np def cosine_similarity(a, b, eps=1e-12): """余弦相似度 = 点积 / 模长乘积，返回 [-1, 1] 之间的浮点数""" a = np.asarray(a, dtype=np.float64) # 统一升到 float64，避免整型数组溢出 b = np.asarray(b, dtype=np.float64) if a.shape != b.shape: raise ValueError(f"维度不一致: {a.shape} vs {b.shape}") na = np.linalg.norm(a) # 默认 2-范数，即欧氏长度 nb = np.linalg.norm(b) if na < eps or nb < eps: # 零向量没有方向，夹角无定义 return 0.0 return float(np.dot(a, b) / (na * nb)) v1 = np.array([1, 2, 3]) v2 = np.array([2, 4, 6]) # 与 v1 同方向、长度翻倍 print(cosine_similarity(v1, v2)) # 期望 1.0

先说dtype=np.float64。如果传进来的是np.array([1, 2, 3])这种整型数组，np.dot的结果还是整型，两个大整数相乘在 32 位平台上会静默溢出，得到一个负数，而代码不会报任何错。转成float64是最省事的规避方式，代价是内存翻倍，对向量这个量级完全无所谓。

再说eps这个参数。零向量的模长是 0，直接做除法会得到nan或者inf，而且这个nan会沿着后续计算一路传染，最后你在损失函数里看到一个nan完全不知道从哪来的——这种问题排查起来就是黑匣子。用eps兜一下，返回 0 表示“方向不可比”，语义上也说得通。

最后是float()强转。np.dot返回的是np.float64标量，它在大多数场景下和 Pythonfloat行为一致，但在json.dumps、某些序列化框架、以及pandas的某些操作里会出问题。养成在函数出口转成原生类型的习惯，能省掉不少莫名其妙的TypeError。

3.2 矩阵乘法：手写三重循环与@的差别不只是快慢

矩阵乘法是这个包里最值得动手写一遍的模块。不是为了在生产里用（生产当然用@），而是为了把“行乘列”这个动作刻进肌肉记忆——后面看注意力机制里的QK^T、看线性层里的xW^T，脑子里能立刻浮现出形状怎么变。

import numpy as np def matmul_naive(A, B): """手写矩阵乘法，严格按 (m,n) x (n,p) -> (m,p) 定义""" A = np.asarray(A, dtype=np.float64) B = np.asarray(B, dtype=np.float64) if A.ndim != 2 or B.ndim != 2: raise ValueError("只接受二维矩阵") if A.shape[1] != B.shape[0]: raise ValueError(f"内维不匹配: {A.shape} x {B.shape}") m, n = A.shape p = B.shape[1] C = np.zeros((m, p), dtype=np.float64) for i in range(m): for j in range(p): s = 0.0 for k in range(n): s += A[i, k] * B[k, j] C[i, j] = s return C A = np.arange(6).reshape(2, 3) B = np.arange(12).reshape(3, 4) print(np.allclose(matmul_naive(A, B), A @ B)) # True

函数前两行的形状校验是重点。A.shape[1] != B.shape[0]这个判断，看起来是废话，但真实项目里维度错误往往发生在三层嵌套之后，报错信息只告诉你“shapes not aligned”，你还得手工回溯每一层的形状。自己写的函数里把两侧形状打进异常信息，排查速度差一个数量级。

np.zeros((m, p), dtype=np.float64)而不是先建空列表再np.array，是因为列表追加再转换的过程会额外分配一次内存，矩阵大一点就能感觉到差别；更重要的是dtype明确，不会因为初始值推断出整型。

三重循环的复杂度是O(m·n·p)，1000×1000 的矩阵乘一次大概要几秒，而A @ B走的是 BLAS，毫秒级。这个差距不是让你去优化循环，而是让你明白：为什么线性代数库要用 Fortran 写、为什么要分块、为什么要考虑缓存局部性。理解了这个量级差，你才会在写模型时下意识避免“在 Python 循环里做矩阵运算”。

比乘法本身更容易踩的是广播。想给矩阵每行减去该行的均值，写成X - X.mean(axis=0)是对的；写成X - X.mean(axis=1)就会因为形状(m,n)和(m,)在最后一维对不上而报错，或者在你调换过来之后得到一个完全错误的结果。判断规则很简单：广播从尾部维度开始对齐，(m, n)和(n,)能广播，(m, n)和(m,)不能，除非你显式加一维变成(m, 1)。

3.3 行列式与高斯消元：解线性方程组的正确姿势

行列式这节在书里通常用来判断矩阵是否可逆，配套代码会给出递归展开的实现和np.linalg.det的对照。递归版只适合教学，O(n!)的复杂度在 10 阶以上就没法用了，实际计算都是先做 LU 分解再取对角线乘积。真正有工程价值的是紧接着的高斯消元——它把“解方程组”这件事从公式变成了可执行流程。

import numpy as np def gauss_solve(A, b, tol=1e-12): """带部分主元的高斯消元，求解 Ax = b""" A = np.array(A, dtype=np.float64, copy=True) # copy=True，绝不能改调用方的数组 b = np.array(b, dtype=np.float64, copy=True) n = b.shape[0] if A.shape != (n, n): raise ValueError(f"系数矩阵形状应为 ({n},{n})，实际 {A.shape}") for k in range(n): # 部分主元：在当前列剩余行里挑绝对值最大的，交换到第 k 行 p = k + int(np.argmax(np.abs(A[k:, k]))) if abs(A[p, k]) < tol: raise ValueError(f"第 {k} 列主元接近 0，矩阵奇异或严重病态") if p != k: A[[k, p]] = A[[p, k]] b[[k, p]] = b[[p, k]] # 消元：把第 k 列下方全部清零 for i in range(k + 1, n): factor = A[i, k] / A[k, k] A[i, k:] -= factor * A[k, k:] b[i] -= factor * b[k] # 回代 x = np.zeros(n, dtype=np.float64) for i in range(n - 1, -1, -1): x[i] = (b[i] - A[i, i + 1:] @ x[i + 1:]) / A[i, i] return x A = np.array([[2., 1., -1.], [-3., -1., 2.], [-2., 1., 2.]]) b = np.array([8., -11., -3.]) print(gauss_solve(A, b)) # [ 2. 3. -1.] print(np.allclose(gauss_solve(A, b), np.linalg.solve(A, b)))

copy=True这一行值得单独说。np.array(A)默认就会复制，但如果写成np.asarray(A)，当输入已经是float64数组时它会直接返回原对象，后面的原地消元就会把调用方的矩阵改掉。这种 bug 特别隐蔽：你调用完求解函数，回头再看A，发现它变成了上三角矩阵，而你还以为是自己哪里算错了。加上显式copy=True，语义清晰，也不会因为numpy版本行为差异翻车。

部分主元那几行是数值稳定性的核心。如果不交换行，遇到主元非常小的情况（比如1e-18），消元时除以它会把误差放大若干数量级，最后解出来的x可能连数量级都不对。tol=1e-12是判断“这个主元算不算 0”的阈值，太小会把病态矩阵当成可解，太大又会把本来能解的正常矩阵拒掉。实际项目里更稳妥的做法是先算条件数np.linalg.cond(A)，量级超过1e10就基本别指望直接解了，该上最小二乘或者正则化。

回代那行A[i, i+1:] @ x[i+1:]用的是切片加内积，比再写一层循环清爽，也顺便复用了前面讲的矩阵运算语义。整个函数是O(n^3)，和np.linalg.solve同阶，但常数大得多，教学价值大于实用价值——真要解方程，一律用np.linalg.solve，它内部走的是 LAPACK 的分块 LU。

4. 特征分解、SVD 与 PCA：从线性变换到降维

4.1 特征值与特征向量：幂迭代为什么能逼近主特征值

特征分解是机器学习里出现频率最高的线性代数工具之一，PCA、谱聚类、图神经网络、马尔可夫链稳态分布，背后都是同一个问题：找一个向量，使得矩阵作用在它上面只改变长度不改变方向。配套代码在这节一般会给出np.linalg.eig的用法，但真正让人有感觉的是手写一次幂迭代——它用最朴素的方式说明了“为什么最大特征值能被迭代出来”。

import numpy as np def power_iteration(A, num_iter=1000, tol=1e-10, seed=0): """幂迭代求主特征值和对应的单位特征向量（要求 |λ1| > |λ2|）""" A = np.asarray(A, dtype=np.float64) n = A.shape[0] rng = np.random.default_rng(seed) v = rng.standard_normal(n) v /= np.linalg.norm(v) # 初始向量归一化，避免数值爆掉 lam = 0.0 for it in range(num_iter): Av = A @ v lam = float(v @ Av) # Rayleigh 商，比单纯取模长更准 norm = np.linalg.norm(Av) if norm == 0: raise ValueError("矩阵把非零向量映到零向量，主特征值为 0") v_new = Av / norm # 每步归一化，否则向量长度指数增长 if v_new @ v < 0: # 方向对齐，消除正负号抖动 v_new = -v_new if np.linalg.norm(v_new - v) < tol: v = v_new break v = v_new return lam, v A = np.array([[4., 1.], [2., 3.]]) lam, v = power_iteration(A) print("lambda:", lam) # 5.0 print("residual:", np.linalg.norm(A @ v - lam * v)) # 应接近 0

rng.standard_normal而不是np.random.rand，是因为前者给出的是标准正态分布、均值 0、有正有负，能避免初始向量恰好与主特征向量正交这种极小概率但确实存在的翻车情况。固定seed是为了结果可复现，调试时能稳定复现问题。

v /= np.linalg.norm(v)每轮都做一次归一化，这是幂迭代能收敛的前提。如果不归一化，每次乘A向量长度都会乘以约λ1，几十轮之后就溢出成inf；归一化之后，长度信息被剥离，只剩下方向在迭代中逐渐对齐主特征向量，而λ1通过 Rayleigh 商v^T A v单独提取出来。

收敛速度取决于|λ2 / λ1|这个比值。比值是 0.9 的话，每轮误差只缩小 10%，可能要几百轮；比值是 0.1 的话，十几轮就够了。这也是为什么 PCA 之前要做中心化——数据中心化之后协方差矩阵的最大特征值和其他特征值差距通常更明显，收敛更快。

注意：幂迭代只能求绝对值最大的那个特征值。如果矩阵有-5和5两个特征值，绝对值相同，迭代会在两个方向之间来回震荡永远不收敛。想求最小特征值，对A的逆做幂迭代；想求中间的，得用位移策略或者直接上np.linalg.eig。

4.2 SVD 与图像压缩：截断到多少维才不算糟蹋数据

奇异值分解是线性代数里最“万能”的工具：任意m×n矩阵都能拆成U Σ V^T，而且这个分解在低秩近似意义下是最优的——保留前k个奇异值，得到的A_k是所有秩为k的矩阵里与A最接近的那个（Frobenius 范数意义下）。图像压缩是把这个性质可视化最直观的方式。

import numpy as np from PIL import Image def svd_compress(path, ks=(10, 50, 120)): img = Image.open(path).convert("L") # 转灰度，单通道矩阵，避免三通道分别处理 A = np.asarray(img, dtype=np.float64) m, n = A.shape U, s, Vt = np.linalg.svd(A, full_matrices=False) # 经济型 SVD，省内存 total_energy = (s ** 2).sum() for k in ks: Ak = (U[:, :k] * s[:k]) @ Vt[:k, :] # 秩-k 近似 ratio = (s[:k] ** 2).sum() / total_energy # 保留的能量占比 out = np.clip(Ak, 0, 255).astype(np.uint8) # 先截断再转 uint8 Image.fromarray(out).save(f"out_k{k}.png") print(f"k={k:3d} energy={ratio:.4f} m*k+k+k*n={m*k + k + k*n}") return U, s, Vt svd_compress("sample.jpg")

full_matrices=False这个参数非常关键。默认的True会返回完整的m×m的U和n×n的Vt，一张 4000×3000 的图，U就是 1600 万个float64，128MB 起步，加上中间过程很容易把内存打满。经济型只返回m×min(m,n)和min(m,n)×n，省下的内存是数量级的差别。

(U[:, :k] * s[:k]) @ Vt[:k, :]这一行里的*是逐元素乘，s[:k]是一维数组，会自动广播到U[:, :k]的每一列上，等价于U[:, :k] @ np.diag(s[:k]) @ Vt[:k, :]但快得多也省内存。

np.clip(Ak, 0, 255)是防翻车的关键一步。Ak是float64，理论上重建结果应该在 0-255 之间，但因为浮点误差，边界像素可能变成-0.0001或255.0001。直接astype(np.uint8)会做截断式转换，负数会绕回成很大的正数（比如-1变成255），画面上就是一片刺眼的噪点。先clip再转，是图像处理里必须养成的习惯。

energy那个比例是回答“保留多少维合适”的量化依据。通常保留 90%~95% 的能量，视觉上就基本看不出损失；k从 10 加到 50，能量可能从 0.6 跳到 0.95，但从 50 加到 120 只多 0.03，性价比断崖式下跌。这也是 PCA 里n_components该怎么设的直观答案。

4.3 PCA 两条路线：协方差特征分解 vs 直接 SVD

PCA 本质上就是“中心化之后做 SVD”，但教材里通常先讲协方差矩阵的特征分解，因为这样能和第 4.1 节的幂迭代串起来。配套代码如果两条路线都给了，一定要对照着跑一遍，能看清一个经常被忽略的事实：两者的数学结果一致，数值稳定性和计算成本差别很大。

import numpy as np def pca_via_eig(X, k): """路线一：中心化 -> 协方差矩阵 -> 特征分解 -> 投影""" X = np.asarray(X, dtype=np.float64) Xc = X - X.mean(axis=0) # 必须中心化，否则第一主成分被均值方向主导 C = (Xc.T @ Xc) / (X.shape[0] - 1) # dxd 协方差矩阵，d 是特征维数 w, V = np.linalg.eigh(C) # 对称矩阵用 eigh，结果实且升序 idx = np.argsort(w)[::-1][:k] # 取最大的 k 个 return Xc @ V[:, idx], w[idx] def pca_via_svd(X, k): """路线二：中心化 -> 直接对数据矩阵做 SVD -> 取前 k 列""" X = np.asarray(X, dtype=np.float64) Xc = X - X.mean(axis=0) U, s, Vt = np.linalg.svd(Xc, full_matrices=False) return U[:, :k] * s[:k], (s[:k] ** 2) / (X.shape[0] - 1) rng = np.random.default_rng(42) X = rng.standard_normal((500, 20)) @ rng.standard_normal((20, 20)) Z1, ev1 = pca_via_eig(X, 3) Z2, ev2 = pca_via_svd(X, 3) print("投影一致:", np.allclose(np.abs(Z1), np.abs(Z2), atol=1e-8)) print("解释方差一致:", np.allclose(ev1, ev2))

路线一的成本是O(n·d²)构造协方差加上O(d³)做特征分解，d是特征维数。路线二的成本是O(m·n·min(m,n))，取决于样本数m和特征数n的相对大小。经验规则是：特征数远小于样本数（比如d=50, n=10000）时，协方差路线更划算；特征数很大（比如图像，d=10000）时，直接 SVD 更稳——因为构造d×d协方差矩阵本身就是O(n·d²)，d=10000时那个矩阵是 800MB。

eigh和eig的选择也值得说一句。协方差矩阵一定是对称半正定的，eigh专门针对这种情况，内部走的是更稳定的算法，返回的特征值一定是实数且升序排列。用eig也能算，但可能返回极小的虚部，而且顺序是乱的，你得自己排序，多一步就多一个出错的机会。

np.allclose(np.abs(Z1), np.abs(Z2))里套了个abs，是因为特征向量的符号本身就不唯一——v和-v都是合法的特征向量，两条路线算出来的方向可能正好相反。这不是 bug，是数学性质。判断结果对不对，看投影后的绝对值或者看解释方差，别去比原始符号，否则会陷入“明明数学一样为什么结果不等”的自我怀疑。

5. 避坑与排查：五个真实翻车现场

5.1AttributeError: module 'numpy' has no attribute 'float'

现象：脚本一跑就报这个错，位置在某个看起来完全正常的dtype=np.float或者np.float(3.14)上。

原因：numpy2.0 移除了np.float、np.int、np.bool、np.object等一批从 Python 内置类型别名过来的符号，它们从 1.20 开始就被标记为废弃。书里的代码写在更早的版本上，自然中招。

解决：全局搜索替换，np.float→float、np.int→int、np.bool→bool。注意别误伤np.float64、np.float32这些正规的类型名——它们是合法的，替换的时候用正则加词边界更稳妥。如果不想改代码，就把环境降回 1.x 并锁死版本。

5.2LinAlgError: Singular matrix出现在求逆那一步

现象：np.linalg.inv(A)直接抛异常，或者更气人的是——不抛异常，但求出来的逆矩阵乘回去不等于单位阵，误差大得离谱。

原因：前者是矩阵行列式严格为 0（比如某两行线性相关）；后者是矩阵“接近奇异”，行列式不是 0 但极小，条件数量级到了1e15，浮点精度根本撑不住求逆这个过程。

解决：求逆之前先看np.linalg.cond(A)，超过1e10就要警惕。绝大多数场景根本不需要显式求逆——解方程用np.linalg.solve，超定或欠定用np.linalg.lstsq，做正则化用(A^T A + λI)^{-1} A^T b的形式。显式求逆是数值线性代数里最不推荐的操作之一，能不用就不用。

5.3 特征向量和书上对不上，符号相反、顺序也乱

现象：手算出来的特征向量是[1, 1]，代码跑出来是[-0.707, -0.707]；或者书上是先λ=5再λ=2，代码跑出来顺序反了。

原因：特征向量的方向本身就不唯一，v和-v、v和2v都是同一个特征值对应的特征向量，不同的库、不同的算法会给出不同的归一化方式和符号。特征值的返回顺序在eig里也是任意的，取决于内部算法，不保证按大小排。

解决：比较时先归一化再取绝对值，用np.allclose(np.abs(v1), np.abs(v2))；需要固定顺序就自己np.argsort排一下；需要固定符号就约定“第一个非零分量为正”，加一行判断翻转。这是数学性质，不是代码错了，别在这上面耗时间。

5.4 图像压缩后整张图发灰、发黑或者出现雪花噪点

现象：SVD 重建后保存的图片，要么整体偏暗，要么边缘出现一片刺眼的白点黑点，和原图完全不像。

原因：float64数组直接astype(np.uint8)是截断转换，超出[0, 255]的值会回绕——-1变成255，256变成0。重建结果因为浮点误差在边界上越界一点点，画面上就是一整片噪点。另一种情况是原图读进来时已经被归一化到[0, 1]，你却按[0, 255]去截断，整张图自然接近全黑。

解决：转换前先np.clip(A, 0, 255)，再astype(np.uint8)。同时确认输入数组的取值范围——PIL读进来是uint8的[0, 255]，matplotlib.image.imread读 PNG 可能是[0, 1]的 float，两者混用必然出问题。归一化数据要保存成图片，先乘回 255。

5.5 matplotlib 不出图、中文变方块、脚本跑完卡住不退

现象：plt.show()之后什么都没发生；图上的中文标题全变成小方块；在服务器上跑脚本，日志打完最后一行就挂在那里不动了。

原因：前两个是后端和字体配置问题，默认后端在无显示环境下无法创建窗口，默认字体里没有中文字形。第三个是交互式后端在等你关窗口，而你在 SSH 里根本关不了。

解决：脚本开头固定matplotlib.use("Agg")，把渲染结果直接写文件，用plt.savefig()代替plt.show()。中文显示加两行plt.rcParams["font.sans-serif"] = ["SimHei"]和plt.rcParams["axes.unicode_minus"] = False（后者解决负号变方块）。如果只是想在 notebook 里看图，用%matplotlib inline，别用notebook后端。

提示：这五类问题占了配套代码跑不动的九成以上。遇到报错先按“版本 → 形状 → 数值范围 → 后端”这个顺序排查，比从头读代码快得多。

6. 二次开发：把手写实现接进真实工程的三条改法

教学代码和工程代码之间隔着一层“防御性”。书里的函数假定输入永远合法、矩阵永远可逆、数据永远干净，真实项目里这三条一条都不成立。我一般按下面三步改造，改完之后这套代码就能当成一个轻量的数值工具箱用。

第一步，抽成独立包并加形状断言。把散落在各章的cosine_similarity、gauss_solve、power_iteration收进一个la_lite/目录，每个函数入口加assert或显式raise，把形状、dtype、NaN检查前置。这一步不增加任何计算量，但能把后面 80% 的调试时间省下来。第二步，用标准库做回归测试，核心就一句np.allclose(handwritten(...), np.linalg.xxx(...))，容差按数据类型定，float64用1e-10，涉及迭代的用1e-6。

import numpy as np import pytest from la_lite.vector import cosine_similarity from la_lite.solve import gauss_solve def test_cosine_against_numpy(): rng = np.random.default_rng(0) for _ in range(50): a, b = rng.standard_normal(8), rng.standard_normal(8) ref = a @ b / (np.linalg.norm(a) * np.linalg.norm(b)) assert cosine_similarity(a, b) == pytest.approx(ref, abs=1e-12) def test_gauss_matches_numpy(): rng = np.random.default_rng(1) for _ in range(20): A = rng.standard_normal((6, 6)) + np.eye(6) * 3 # 加对角占优，降低病态概率 b = rng.standard_normal(6) assert np.allclose(gauss_solve(A, b), np.linalg.solve(A, b), atol=1e-9)

+ np.eye(6) * 3这一步是故意的，随机矩阵有一定概率接近奇异，测试用例本身就不该依赖“运气好”。pytest.approx用于标量比较，np.allclose用于数组，两者别混用。这套测试跑通，说明你的手写实现至少在数值上是可信的，后面替换成库函数也有对照基准。

第三步，数值稳定性上做取舍。工程里默认用solve不用inv，默认用lstsq处理可能秩亏的方程，默认用svd而不是特征分解处理非对称矩阵。手写版留着做教学和兜底，主路径走 LAPACK。有一个具体的量化门槛我一直在用：np.linalg.cond(A) > 1e12就直接放弃直接解，改用带正则的最小二乘，宁可引入一点偏差，也不要一个误差放大到面目全非的解。

从那以后我每次把一个教学性质的数值代码往工程里搬，都会先做一件事：拿一组已知答案的输入跑一遍，把输出和标准库结果并排放，差值超过预期就停下来查，不往下写。这个习惯让我少踩了很多“结果不对但代码不报错”的坑。希望帮到你。

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