矩阵从会考试到会用：运算、数值稳定性与代码实战

矩阵这两个字很多人是在大二那门线性代数课上第一次见到的然后在大三的数据分析课上第二次见到第三次见到是在工作第一年的代码评审里——同事甩过来一行np.linalg.solve(A, b)控制台甩回来一个LinAlgError: Singular matrix你盯着屏幕想了半分钟才意识到问题不在代码在你对矩阵这个东西本身的理解。矩阵线性代数说白了就是三样东西一张按行列排好的数表、一台把向量搬来搬去的机器、一套解方程组的工具箱。它的用途远比课本上写的广泛图像是矩阵推荐系统的评分表是矩阵神经网络的权重是矩阵结构力学里的刚度矩阵、电路里的节点导纳矩阵、统计里的协方差矩阵全是它。这篇文章想干的事很具体——把矩阵从会考试推进到会用从手算流程讲到代码实现从形状约定讲到条件数把中间那些课本上不写、但工程上天天坑人的细节摊开讲。适合刚学完线性代数想找落地感觉的学生也适合工作中偶尔要碰数值计算、但一看到矩阵就心里发虚的开发者。1. 先把矩阵看明白一张表、一台机器、一套方程组1.1 两个视角数据表和线性变换初学者最容易犯的错是把矩阵只当成二维数组。这个理解不算错但只覆盖了一半。我一般建议同时建立两个视角视场景切换。第一个视角是数据表。一个 m 行 n 列的矩阵就是 m 条记录、n 个特征。比如一张用户评分表行是用户、列是商品元素是打的分空缺表示没看过。这个视角下你关心的是形状、稀疏性、每一列的量纲操作是取子块、转置、按列求均值。做数据分析的人九成时间用的是这个视角。第二个视角是线性变换。把矩阵看成一个函数喂进去一个向量吐出来另一个向量。一个 2×2 的矩阵作用在平面上可以把圆压成椭圆、把正方形旋转一个角度、把图形沿某条线做镜像。这个视角下你关心的是特征值哪些方向只被拉伸不被转向、行列式面积缩放了多少倍、秩维度塌缩了多少。做算法和图形的人靠这个视角建立直觉。两个视角其实是同一件事的两面关键在于矩阵的列就是基向量的像。一个 2×2 矩阵[[a, b], [c, d]]第一列(a, c)是 x 轴单位向量被搬到的位置第二列(b, d)是 y 轴单位向量被搬到的位置。这条结论一旦记住后面看任何矩阵都不再是一堆数字而是一组落脚点。我自己判断一个矩阵大概长什么样第一件事就是看它的列列的长度大致反映各方向的拉伸幅度列之间的夹角反映扭曲程度。1.2 为什么工程里到处都是矩阵矩阵之所以统治了数值计算核心原因只有一个它把批量线性组合这件事写成了紧凑、可并行、可分析的符号。拿图像来说一张 1024×1024 的灰度图本质就是一个 1024 行 1024 列的矩阵每个元素是 0 到 255 的亮度。做一次模糊就是把这个矩阵和一个小尺寸的核矩阵做卷积做一次边缘检测同样是换一个核。你以为在调滤镜其实在跑矩阵运算。拿拟合来说你有一百个观测点想拟合一条直线 y kx b课本上要列正规方程组工程上就是构造一个 100×2 的矩阵 A第一列是 x第二列全是 1然后解A^T A · [k, b]^T A^T y。换成拟合二次曲线就是把 A 扩成三列代码一行不改。拿推荐系统来说用户-商品评分矩阵通常极其稀疏要做的是矩阵分解——把它拆成两个低秩矩阵的乘积用已知评分去预测未知评分。核心假设是用户的偏好维度远少于商品数量落到数学上就是这个矩阵近似低秩。还有一类更隐蔽的用法用矩阵表达约束。线性规划里所有的约束都是Ax ≤ b的形式电路分析里基尔霍夫定律直接写成矩阵方程结构力学里上千个节点的受力平衡是一个巨大的稀疏线性系统。这些场景的共同点是变量多、方程多、但每个方程只涉及少数变量所以矩阵是稀疏的求解器要专门优化。你要是拿个稠密矩阵去解内存先炸。1.3 形状约定与命名m×n 到底谁在前谁在后形状是矩阵世界里最容易出错、也最不该出错的地方。约定很简单m×n 表示 m 行 n 列先行后列。一个 3×2 的矩阵有 3 行 2 列共 6 个元素。我建议养成三个习惯。第一写索引时想清楚是从 0 开始还是从 1 开始。数学教材用 a₁₁ 表示第一行第一列代码里是A[0, 0]。两套体系混着用调试时能把人逼疯。我的做法是数学推导全部用 1 起始一进代码立刻切换到 0 起始中间不来回横跳。第二每个矩阵都标注形状。推导时在草稿纸上写A(3×2) · B(2×4) C(3×4)养成习惯以后形状不匹配的乘法你一眼就能看出来不用等编译器报错。这个习惯在推导多层网络的反向传播时尤其救命。第三区分行向量和列向量。课本里向量默认是列向量写成 n×1。代码里numpy的一维数组既不是行也不是列做乘法时会自动按广播规则处理很容易埋雷。我现在的做法是只要涉及矩阵乘法一律显式写成二维x.reshape(-1, 1)或者x.reshape(1, -1)宁可在形状上多打几个字也不要让广播机制替我猜。记号含义典型形状容易混淆的点A(m×n)m 行 n 列矩阵3×2先记行数再记列数A^T转置行列互换2×3转置后形状反过来了x(n×1)列向量3×1一维数组在代码里不算列向量I_nn 阶单位矩阵3×3主对角线为 1其余为 00零矩阵任意形状不同不能相加2. 核心运算的规则与手算流程2.1 加法、数乘、转置看着简单错得最多矩阵加法要求两个矩阵形状完全一致对应位置相加。这一条谁都知道但实际错误往往出在广播上两个形状是 (3,1) 和 (3,) 的数组在你脑子里都是三个数的向量程序却不这么认为结果要么报错要么悄悄算出一个 (3,3) 的矩阵。这类错误在生产环境里排查起来非常费劲因为不报错。数乘是每个元素都乘同一个标量没有形状限制。要注意的是它的几何含义数乘改变的是所有方向的长度正数不翻转方向负数会把图形整体翻转 180 度。转置是把第 i 行第 j 列的元素搬到第 j 行第 i 列形状从 m×n 变成 n×m。转置的性质里最常用的是(AB)^T B^T A^T——注意顺序反了。这条性质在推导向量内积、协方差矩阵、梯度公式时反复出现顺序写反是新手高频错误。我见过有人在推导梯度的时候把顺序搞错最后整个模型的符号是反的训练完全不收敛查了两天才发现问题出在一行转置上。提示转置的反序规律值得单独记住因为它背后是先做 B 变换再做 A 变换反过来看就要先撤销 A 再撤销 B和穿脱衣服的顺序一个道理。2.2 矩阵乘法为什么是行乘列这是线性代数里最反直觉的一条规则也是最多人死记硬背、最不理解的一条。规则本身是C AB 中C 的第 i 行第 j 列元素等于 A 的第 i 行与 B 的第 j 列对应相乘再求和。形状上A 必须是 m×pB 必须是 p×n中间的 p 必须对上结果 C 是 m×n。为什么非得这么定用一个 2×2 的例子算给你看。设 A 是变换先做某件事B 是变换再做某件事那么 AB 应该表示先 B 后 A。关键是A [[1, 2], B [[0, 1], [3, 4]] [1, 0]]B 把 x 轴单位向量 (1,0) 送到 (0,1)把 y 轴单位向量 (0,1) 送到 (1,0)也就是交换两个坐标。那么 AB 的第一列应该是B 把 x 轴送到的位置再被 A 作用B 把 (1,0) 送到 (0,1)A 作用在 (0,1) 上得到 (2, 4)。所以你算出来的 AB 第一列是 (2,4)。手工按行乘列验证第一行第一列 1×0 2×1 2 ✓第二行第一列 3×0 4×1 4 ✓也就是说行乘列不是谁拍脑袋规定的它是复合变换的坐标计算结果。只要接受矩阵的列是基向量的像这一条行乘列就是唯一合理的算法。反过来如果按对应元素相乘定义得到的就不是复合变换只是个逐点缩放几何上没有意义。手算时我推荐一个校验技巧先算结果的形状再逐行逐列填。比如 A(2×3) 乘 B(3×2) 得到 C(2×2)一共四个数一个一个填。填完以后用一个随机向量去验算A (B x)和(A B) x是否相等两边都是同一结果说明乘法没算错。这个验算成本极低能拦住绝大部分手误。还要强调一条矩阵乘法一般不可交换AB 和 BA 通常不相等甚至形状都可能对不上。这不是数学家的恶作剧是先旋转再缩放和先缩放再旋转本来就不同。2.3 行列式与逆矩阵能不能求逆一句话判断行列式的几何含义非常直白它表示矩阵对应的线性变换对体积的缩放倍数。二维里是面积倍数三维里是体积倍数。行列式为正说明变换保持定向左右手不变为负说明翻转了定向等于零说明至少有一个维度被压扁了——原本的平面被压成一条线信息丢失不可逆转。2×2 的行列式最好记det([[a,b],[c,d]]) ad - bc。3×3 用萨吕斯法则或者按行展开实际工作中手算 3×3 已经很少了我一般直接上代码只在推导时用符号算。有一点要额外注意行列式可能非常大或非常小一个 20×20 的随机矩阵行列式可能大到超出浮点数表示范围也可能小到直接下溢成 0。所以工程里判断是否接近奇异看行列式是不靠谱的应该看条件数。逆矩阵只有方阵才谈得上而且必须行列式不为零满秩。2×2 的逆有现成公式A [[a, b], A^-1 1/(ad-bc) · [[ d, -b], [c, d]] [-c, a]]记忆要点主对角线元素互换副对角线元素取负再除以行列式。三阶以上就得用高斯-约当消元手算量大基本交给代码。这里有个经典坑不要为了解方程去求逆。方程 Ax b 的正确解法是消元LU 分解直接解而x A^-1 b在数值上更慢、误差更大。行列式接近零时求逆出来的结果会放得极大稍微一点输入扰动输出就完全不靠谱。我在一次小规模回归计算里就吃过这个亏矩阵 A^T A 的条件数到了 1e10求逆的结果里有一个系数是 1.2e-3 但符号是反的看起来还挺正常直到用另一组数据复现才发现模型一直在乱跳。后来改成np.linalg.lstsq一步解决这个坑从此记忆深刻。3. 手算到代码用 Python 把矩阵跑通3.1 环境准备与最小可用示例我用的是一套最省事的组合Python 3.10 以上加 NumPy交互环境随意。安装只有一条命令pip install numpy下面这段代码把矩阵创建、形状查看、乘法、转置、行列式、求逆、解方程一次跑完可以直接复制到一个脚本里验证import numpy as np A np.array([[2.0, 1.0], [1.0, 3.0]]) b np.array([5.0, 6.0]) print(shape:, A.shape) print(T:, A.T) print(det:, np.linalg.det(A)) print(cond:, np.linalg.cond(A)) x np.linalg.solve(A, b) print(solution:, x) print(residual:, A x - b)输出应当看到行列式是 5.0解是[1.8, 1.4]残差接近零。手算验证2×1.8 1×1.4 51×1.8 3×1.4 6两边都对得上。这段代码里我特意加了cond它返回条件数。经验上条件数在 1 到 100 之间解是可信的到 1e6 就开始丢失有效数字到 1e12 以上双精度浮点基本给不出有意义的结果。这个数比行列式有用得多建议养成习惯解方程之前先看一眼。3.2 解线性方程组的三条路同一个 Ax b代码里有三条常见路径选错会付出代价。第一条是np.linalg.solve(A, b)内部走 LU 分解。要求 A 是方阵且满秩速度最快、数值最稳能选就选它。第二条是np.linalg.inv(A) b。能用但没必要。它先老老实实求一遍逆再做一次乘法运算量大约是前者的两倍以上误差还可能更大。除非你确实需要这个逆矩阵本身比如后面还要反复用否则别走这条路。第三条是np.linalg.lstsq(A, b, rcondNone)最小二乘 / 伪逆。A 不方、不满秩都能跑代价是速度慢。实际工作中我有一半以上的时间在用这个因为现实里的方程组经常是超定的方程比未知数多本来就无精确解只能求近似解。方法适用条件相对速度典型场景solve方阵、满秩最快已知精确模型求解inv后乘方阵、满秩较慢需要反复使用逆矩阵lstsq任意形状最慢但稳拟合、超定方程、秩亏pinv任意形状慢需要伪逆本身还有一条经验A^T A 这一步能不做就不做。很多人解最小二乘的本能反应是构造正规方程(A^T A) x A^T b再solve数学上没错但 A^T A 会把条件数平方A 的条件数是 1e4 时A^T A 就到 1e8误差直接放大。直接用lstsq内部走 SVD或者 QR 分解稳定性好一个档位。3.3 特征值分解和 SVD 的实操特征值分解解决的是这个矩阵在哪些方向上只拉伸不转向。2×2 矩阵 A [[2,1],[1,3]]特征多项式是 λ² - 5λ 5 0解得 λ (5 ± √5)/2约等于 3.618 和 1.382。代码一行vals, vecs np.linalg.eigh(A) # A 对称时用 eigh更快更稳 print(vals) # [1.382, 3.618] print(vecs)注意这里用eigh而不是eig因为 A 是对称的。对称矩阵的特征值必为实数特征向量正交eigh就是利用这个性质做的优化速度快、结果稳。如果你的矩阵不对称只能用eig返回的可能是复数这一点要提前有心理准备。奇异值分解SVD把任意矩阵拆成A U Σ V^T对形状没有任何要求是矩阵分解里的万金油。奇异值就是 Σ 对角线上的非负数按从大到小排列。它的含义是A 把一组正交方向按比例 σ₁, σ₂... 拉伸到另一组正交方向。前几个奇异值占比越大说明这个矩阵的有效维度越低。U, s, Vt np.linalg.svd(A, full_matricesFalse) print(s) # 奇异值序列 rank_est (s 1e-8 * s[0]).sum() print(numerical rank:, rank_est) print(energy ratio:, (s**2).cumsum() / (s**2).sum())最后一行算的是累计能量占比[0.9, 1.0]之类的输出告诉你前 k 个奇异值能解释多少信息这是做降维和压缩时最重要的一个参考量。3.4 数值稳定性条件数是怎么坑人的数值稳定性这个概念课本上讲得少工程上遇到的概率却极高。核心结论有三条。第一浮点数不精确。0.1 在二进制里是无限循环小数存成双精度就带了一点误差一万次累加以后误差能到 1e-12 量级。这不是 bug是物理事实所有应对手段都是在和它共处。第二减去两个相近的数会丢有效位。1.0000001 - 1.0000000结果只剩一位有效数字。矩阵消元过程中主元很小的时候就会出现这种情况所以实际实现要做部分主元选取每一步挑当前列绝对值最大的元素当主元交换行使之就位。NumPy 的solve内部已经做了自己手写高斯消元时千万别省这一步。第三病态问题无解。一个条件数 1e15 的系统输入里最后一位的舍入误差会被放大到和真值同样的量级这时候无论用什么算法都救不回来只能改模型——加正则化、减变量、换数据。加一点轻微的正则化A^T A λI往往能把条件数压下去好几个数量级代价是引入一点偏差。这是工程上非常常用的折中我在做传感器标定时就靠这一招把不稳定的系数按住了。lam 1e-6 A_reg A.T A lam * np.eye(A.shape[1])4. 实操项目用矩阵做最小二乘拟合与图像压缩4.1 最小二乘拟合从公式到可复现代码这个例子我拿来当矩阵入门实战很多次了因为它把形状、转置、求解、残差全串了一遍。目标是拟合一条直线 y kx b。给定五个观测点 (1,2.1)、(2,3.9)、(3,6.2)、(4,7.8)、(5,10.1)肉眼能看出大概有斜率但点不在一条线上需要最小二乘。构造矩阵A 是 5×2第一列是 x第二列全 1b 是 5×1。要解的是 A · [k, b]^T ≈ y。代码import numpy as np x np.array([1, 2, 3, 4, 5], dtypefloat) y np.array([2.1, 3.9, 6.2, 7.8, 10.1]) A np.column_stack([x, np.ones_like(x)]) print(A shape:, A.shape) # (5, 2) coef, residuals, rank, sv np.linalg.lstsq(A, y, rcondNone) k, b0 coef print(fy {k:.4f} x {b0:.4f}) pred A coef resid y - pred print(SSE:, (resid**2).sum()) print(RMSE:, np.sqrt((resid**2).mean()))跑出来斜率大概 1.98 左右截距接近 0.1 出头SSE 在 0.1 量级。这里有几个值得说的细节。为什么用column_stack而不是hstackcolumn_stack把一维数组当列拼hstack会把它们首尾相接成一维。形状错了后面全错这是新手最容易卡住的地方。rcondNone是什么意思控制判定奇异值多小算零的阈值。新版 NumPy 建议显式传None走机器精度的默认规则。不传会收到警告虽然不影响结果但日志里一堆警告很难看。残差的单位SSE 是平方和量纲是 y 的平方不好直观理解RMSE 开根号后和 y 同量纲可以直接说平均偏差零点几。我汇报结果时一律用 RMSE不用 SSE。想扩成二次拟合只需把 A 改成三列np.column_stack([x**2, x, np.ones_like(x)])其余代码一行不改。这就是矩阵写法的好处——模型阶数变了代码结构不变。4.2 图像压缩SVD 去掉不重要的方向这个实验能让你直观看到低秩近似的威力。思路是一张灰度图是一个矩阵做 SVD 得到U Σ V^T只保留前 k 个奇异值用U[:, :k] diag(s[:k]) Vt[:k, :]重建就得到一个压缩后的近似图。import numpy as np from PIL import Image img np.array(Image.open(sample_gray.png).convert(L), dtypefloat) U, s, Vt np.linalg.svd(img, full_matricesFalse) print(original shape:, img.shape, singular values:, s.size) for k in (5, 20, 50, 100): approx (U[:, :k] * s[:k]) Vt[:k, :] err np.linalg.norm(img - approx) / np.linalg.norm(img) keep (img.size and (U[:, :k].size s[:k].size Vt[:k, :].size)) / img.size print(fk{k:3d} 相对误差{err:.4f} 存储占比{keep:.2%})一张 512×512 的图原始存储 262144 个数。取 k50需要存 512×50 50 50×512 51712 个数只有原来的 20% 左右而肉眼看上去几乎没差别。相对误差一般在 0.05 上下取决于图像内容的复杂度。这里有两个实操心得。第一U[:, :k] * s[:k]用的是广播乘法把每一列乘以对应的奇异值等价于U[:, :k] np.diag(s[:k])但省掉了构造对角矩阵的开销速度和内存都更友好。第二存储占比要按数字个数算不能按文件字节算因为 PNG 本身有压缩你拿 SVD 重建后的数组再存成 PNG文件大小未必更小。这个实验的意义在于理解低秩结构不是真的做图片压缩工具。顺带说一句SVD 在这个场景里比特征值分解合适得多因为图像矩阵不是方阵甚至长宽比可能很怪特征值分解根本用不了。4.3 结果验证与指标记录做完上面两个实验我建议养成一个习惯把关键指标记成一张表别只看一眼输出就完事。表格能帮你发现参数稍微一改结果就崩的情况。实验关键参数观察指标我常用的合理区间直线拟合阶数 1RMSE与数据噪声同量级二次拟合阶数 2RMSE 是否明显下降下降不足 10% 就别升阶图像压缩k相对误差 / 存储占比k 取到奇异值总数的 10% 左右线性求解条件数残差范数残差应远小于解的范数用这张表回头看很多看起来正常的结果会立刻暴露问题。比如二次拟合的 RMSE 只比直线低 0.3%那说明二次项完全是噪声拟合属于典型过拟合不该采纳。再比如某个线性系统的残差范数居然是解范数的十分之一那基本可以断定矩阵病态解不可信。5. 常见问题与排查手册5.1 手算阶段的排错流程手算矩阵题错的地方高度集中。我总结了一套从后往前的检查顺序比从头重算快得多。第一步查形状。结果的形状对不对A 是 2×3、B 是 3×3结果必须是 2×3。形状错了后面全白算。第二步查量纲级别的异常。结果里出现一个特别大或者特别小的数通常意味着某个元素抄错了尤其是符号。第三步用特例验证。令某个变量为 0 或者 1看结果是否符合直觉。行列式算出来是 0那就在草稿上确认一下是不是某两行成比例了。第四步反向验算。求出逆矩阵就乘回去看是不是单位矩阵解出方程组就代回原方程做特征分解就验证A v λ v。反向验算的成本通常只有重算的三分之一命中率却高得多。注意反向验算一定要用原始数据不要用你中间算出来的那些数。中间量错了正着验和反着验会一起错验不出来。5.2 代码报错的对照表跑矩阵代码时报的错来来回回就那么几个。整理成表遇到直接查。报错信息根本原因处理方式shapes not aligned内维不匹配如 (3,2)(3,2)检查形状必要时转置或 reshapeSingular matrix矩阵不满秩行列式为 0改用lstsq或加正则项Last 2 dimensions must be square传给det/inv的不是方阵确认是否该用svd或pinv结果形状是 (3,3) 而非 (3,)一维数组触发广播显式 reshape 成 (3,1) 或 (1,3)结果全是nan除零、溢出、或输入本身含 nan检查cond检查分母特征值出现复数矩阵不对称用了eig正常现象若要实数需确认矩阵对称并用eigh我最常踩的是最后一行往上一行那个一维数组广播。w x里如果 w 是 (3,)、x 是 (3,)得到的是一个标量内积如果 x 是 (3,1)得到的是 (3,1)。同一个表达式形状差一个维度含义完全不同而且很多时候不报错。所以我的代码里现在有一条硬规矩凡是要做矩阵乘法的变量一律保持二维不做一维数组。5.3 秩亏、病态和那些课本不写的事最后聊几个真正让人头疼的场景。秩亏。矩阵的行或列线性相关有效维度不足。比如你在做回归两个特征其实是同一个量的两种度量那设计矩阵就秩亏正规方程有无穷多解。这时候solve会报错lstsq会给一个最小范数解。但更该做的是回数据层面把冗余特征删掉一个。工具能兜底但兜底的结果通常不是你想要的。病态。矩阵满秩但接近奇异条件数很大。表现是换一组采集的数据解出来的系数符号都变了。这类问题靠换算法解决不了只能靠正则化、减特征、或者增加样本量。我一般的判断标准是条件数超过 1e8 就停下来看看数据不硬算。稀疏性被忽略。实际工程里的矩阵经常 99% 是零比如图的邻接矩阵、有限元刚度矩阵。拿稠密存储去算内存和时间都浪费几个数量级。Python 里该用scipy.sparsescipy.sparse.linalg.spsolve解稀疏系统。这一点课本完全不提但工作中一旦遇到就会明白差距有多大。整数类型溢出。np.array([1,2,3])默认是 int64做矩阵乘法时如果数值大可能溢出成负数。我踩过一次一个规模不大的整数矩阵相乘结果里出现了一个 -9223372036854775808看上去像乱码实际就是溢出。养成习惯做数值计算建数组时就写dtypefloat。重复求逆。循环里对同一个矩阵反复求逆或反复分解性能会非常惨。能分解一次复用就复用LU 分解和 Cholesky 分解都支持这样用。做完 Cholesky 分解后续每次求解的代价就降到 O(n²)而不是每次 O(n³)。我个人在实际操作中的体会是矩阵这东西的难点从来不在公式本身而在形状、数值、语义这三件事上。形状决定能不能算数值决定算出来可不可信语义决定算出来该不该用。课本把九成功夫花在公式推导上剩下的三成要自己补多写几行带print(shape)的代码多看一眼条件数多把结果代回原始数据验一遍。这三件小事坚持下来矩阵就不再是一门考完就忘的课而是一个能随时拿起来用的工具。