ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

NumPy 线性代数的具体实现

NumPy 线性代数的具体实现 前言先明确一个前提NumPy 不是 Python 标准库它是第三方库使用前需要pip install numpy。本机没有 Python 解释器也没有安装 NumPy所以本文的示例无法在本机运行验证只能逐行人工推演所有函数名、参数表和返回结构都以 NumPy 官方文档为准请你以自己的环境为准。另外NumPy 里的「矩阵」指的是二维的numpy.array对象官方已经不再推荐使用它的numpy.matrix子类即使在纯线性代数场景下也一样。第二个要纠正的误解是「线性代数函数都在numpy.linalg里」。不是矩阵乘法dot、matmul、inner、outer、kron都是顶层的numpy函数numpy.dot、numpy.matmul不是numpy.linalg的成员而求解、分解、特征值、行列式、秩、条件数这些才是numpy.linalg的强项。把这两类混为一谈是查文档时找不到函数的主要原因。第三个误解是「求逆就能解方程」。数学上确实可以写成x A⁻¹b但数值计算上不该这么做先求逆再相乘既慢又不稳正确做法是用numpy.linalg.solve。这一点会在第三节展开它是本文最实用的一条。本文按「乘法 → 解方程 → 分解 → 矩阵特征量」的顺序讲清numpy.linalg里各个函数在做什么、什么时候用哪个。一、dot/matmul/三者的关系先讲怎么乘。官方文档说得很明确运算符从 NumPy 1.10.0 引入在计算二维数组之间的矩阵乘法时它优于其他写法numpy.matmul就是的实现。也就是说A B与np.matmul(A, B)是一回事。numpy.dot也能做矩阵乘法但它和matmul并不是同义词差别在高维和标量上写法二维数组高维数组标量参与A B/np.matmul(A, B)矩阵乘法按「矩阵栈」处理并广播不允许标量np.dot(A, B)矩阵乘法对A的最后一轴与B的倒数第二轴求和允许退化为数乘所以一条实用规则是做矩阵乘法就用语义最清楚也不会在高维时踩到dot那种「按轴求和」的隐含行为。dot留给你明确需要它那套广义语义的地方。# 适用于 Python 3.8 且已安装 NumPy以官方文档为准import numpy as npA np.array([[1., 2.], [3., 4.]])B np.array([[5., 6.], [7., 8.]])print(A B) # 矩阵乘法print(np.matmul(A, B)) # 与上一行完全等价print(A * B) # 逐元素相乘不是矩阵乘法最后一行要单独强调*是逐元素相乘element-wise不是矩阵乘法。把*当成矩阵乘法是初学阶段最典型的一类 bug。顺带说一句 Python 2 的事Python 2.7 已于 2020 年 1 月 1 日停止维护新版 NumPy 也早已不支持 Python 2。网上老代码里的print语句、xrange、numpy.matrix的种种写法在 Python 3 的新版本里要么语法错误要么已不推荐不要照抄。二、解线性方程组solve与lstsq签名是numpy.linalg.solve(a, b)解的是a x b。官方文档给出的若干约束值得逐条记下a必须是方阵且满秩各行各列线性无关否则抛LinAlgError它内部走 LAPACK 的_gesv例程广播规则生效a可以是「堆叠」的多个方阵——这也是 NumPy 相对 SciPy 的一个优势numpy.linalg.solve能一次处理一批矩阵。还有一条版本差异要注意从 NumPy 2.0 起b只有恰好是一维ndim 1时才被当作形状(M,)的列向量其他情况一律按(M, K)的矩阵栈处理在此之前的版本若b.ndim等于a.ndim - 1就会按向量栈处理。跨版本写代码时这一点要留意。# 适用于 Python 3.8 且已安装 NumPy以官方文档为准import numpy as npA np.array([[3., 1.], [1., 2.]])b np.array([9., 8.])x np.linalg.solve(A, b)print(x) # 方程组的解print(np.allclose(A x, b)) # True验证一下如果a不是方阵或者虽然是方阵但奇异不满秩solve会直接报错——这时候不该去凑一个逆而应该换lstsq。numpy.linalg.lstsq(a, b, rcondNone)求的是最小二乘解方程可以是不定的、恰定的或超定的。如果a是方阵且满秩返回的解在舍入误差范围内就是精确解否则它最小化||b - a x||的欧氏二范数如果存在多个最小化解返回其中二范数最小的那个。它的返回值是一个四元组这一点经常被记错返回项形状含义x(N,)或(N, K)最小二乘解b是二维时解按列排residuals(1,)、(K,)或(0,)每列残差的平方和b是一维时形状为(1,)rankint矩阵a的秩s(min(M, N),)a的奇异值rcond是奇异值的截断比例小于「rcond乘以最大奇异值」的奇异值在判定秩时被当作零。它的默认值在 NumPy 2.0 变了以前默认是-1并用一次警告提示即将改变现在是None表示使用「机器精度乘以max(M, N)」显式传-1则是使用机器精度。这个改动会实实在在影响秩判定和返回的解所以显式写上rcondNone反而更稳。# 适用于 Python 3.8 且已安装 NumPy以官方文档为准import numpy as npA np.array([[0., 1.], [1., 1.], [2., 1.]])b np.array([1., 2., 3.])x, residuals, rank, s np.linalg.lstsq(A, b, rcondNone)print(x.shape) # (2,)print(residuals.shape) # (1,)print(rank) # 2三、为什么「能用solve就不要用inv再相乘」这是本文最值得记住的一条。先说清楚两种写法在数学上等价# 适用于 Python 3.8 且已安装 NumPy以官方文档为准import numpy as npA np.array([[3., 1.], [1., 2.]])b np.array([9., 8.])x1 np.linalg.solve(A, b) # 推荐x2 np.linalg.inv(A) b # 不推荐数值行为更差它们的结果在理想数学世界里一样在A可逆时但计算机用的是有限精度浮点两条路径的误差传播完全不同。原因有两层。第一层是步骤solve内部走的是分解后求解的路径LAPACK 的_gesv会做带主元的分解再回代一步到位而「先inv再相乘」相当于先完整构造出A⁻¹的每一个元素再做一次矩阵乘——显式求逆这个动作本身就会放大误差。第二层是代价显式求逆要做一次完整的矩阵分解加一次求逆之后还要再做一次乘法直接解只做一次分解加一次回代。矩阵越大差距越明显。还有一个容易被忽略的点如果A接近奇异显式求逆的结果可能完全没有意义——你会得到一组看起来很正常的数字但它们和真解差得很远而且没有任何报错。相比之下solve对真正奇异的输入会抛LinAlgError至少让你知道出事了。那numpy.linalg.inv(a)什么时候用——当你真的需要那个逆矩阵本身的时候比如后续要反复用它做别的计算、或者要把它交给另一个库。仅仅为了解一次方程组而求逆是典型的用错工具。numpy.linalg.pinv(a, rcond..., hermitian...)是伪逆Moore–Penrose 广义逆用于非方阵或奇异的矩阵。它的签名里除了rcond还有hermitian以及较新版本加入的rtol。具体参数的默认值和版本差异请以官方文档为准因为rcond/rtol的默认值在近几个版本调整过。四、分解与特征值det/eig/svd/qr这一组是「把矩阵拆开」的工具函数作用返回值要点numpy.linalg.det(a)行列式返回标量numpy.linalg.slogdet(a)符号与对数行列式返回(sign, logabsdet)numpy.linalg.eig(a)特征值与右特征向量一般方阵返回(w, v)numpy.linalg.eigh(a[, UPLO])对称厄米矩阵的特征值返回(w, v)只用下三角或上三角numpy.linalg.svd(a[, full_matrices, compute_uv, hermitian])奇异值分解默认返回(U, S, Vh)numpy.linalg.qr(a[, mode])QR 分解默认modereduced返回(Q, R)numpy.linalg.cholesky(a, /, *[, upper])Cholesky 分解要求对称正定几个容易记错的地方svd返回的第三个是Vh不是V。它已经转置准确说是共轭转置过了所以重构矩阵时应写U np.diag(S) Vh而不是再对Vh取转置。这一条错了通常不会报错只会得到一个形状对不上的结果或错误的重构。full_matrices决定U和Vh的形状。默认True时返回完整的方阵设成False时返回「经济型」分解U为(M, K)、Vh为(K, N)其中K min(M, N)。大矩阵上通常用False省内存。eig与eigh不能混用。一般矩阵用eig实对称或复厄米矩阵用eigh后者会利用对称性、更快也更稳并且返回的特征值默认是升序的。eigh还有UPLO参数决定读下三角L默认还是上三角U——只读一半意味着另一半被当成镜像输入不对称时结果会出乎意料。特征值的顺序不保证。eig返回的特征值没有规定顺序不要假设「第一个是最大的」。# 适用于 Python 3.8 且已安装 NumPy以官方文档为准import numpy as npA np.array([[0., 1.], [1., 1.], [2., 1.]])U, S, Vh np.linalg.svd(A, full_matricesFalse)print(U.shape, S.shape, Vh.shape) # (3, 2) (2,) (2,)Q, R np.linalg.qr(A)print(Q.shape, R.shape) # (3, 2) (2,)D np.array([[2., 0.], [0., 3.]])w, v np.linalg.eig(D)print(sorted(w)) # 特征值为 2 和 3顺序不保证所以先排序五、矩阵秩与条件数这两个量回答了「这个方程组好不好解」的问题做数值计算时经常要先用它们探路。numpy.linalg.matrix_rank(A[, tol, hermitian, rtol])用 SVD 求秩。为什么不用「数非零奇异值」这种朴素办法因为浮点误差会让本该为零的奇异值变成1e-16这种小量必须有一个容差。tol就是这个阈值不传时由最大奇异值和矩阵维度推算得出。numpy.linalg.cond(x[, p])求条件数。它衡量的是「输入的微小扰动会被放大多少倍」条件数大说明矩阵接近奇异解对输入极其敏感——b上一点点测量误差可能让x面目全非。看到大的条件数就该对结果保持怀疑而不是照单全收。# 适用于 Python 3.8 且已安装 NumPy以官方文档为准import numpy as npA np.array([[1., 2.], [2., 4.]]) # 两行线性相关矩阵奇异print(np.linalg.matrix_rank(A)) # 1print(np.linalg.cond(A)) # 很大的数1 除以 0 的极限情形是 inf顺便区分一下范数numpy.linalg.norm(x[, ord, axis, keepdims])同时管向量范数和矩阵范数靠ord和axis决定语义。如果你只想明确地要向量范数或矩阵范数官方提供了numpy.linalg.vector_norm和numpy.linalg.matrix_norm这两个更专一的函数语义比norm少一层歧义。常见坑点❌ 用np.linalg.inv(A) b解方程组认为它和solve等价。✅ 直接np.linalg.solve(A, b)显式求逆会放大误差、代价更高且接近奇异时得到的错误结果不会报错。❌ 把A * B当成矩阵乘法。✅*是逐元素相乘矩阵乘法用或np.matmul。❌ 认为numpy.dot和numpy.matmul完全等价在高维数组上随便换着用。✅ 二维时结果相同高维时dot按轴求和、matmul按矩阵栈广播做矩阵乘法优先用。❌ 调用np.linalg.lstsq只接一个返回值。✅ 它返回(x, residuals, rank, s)四元组residuals在b为一维时形状是(1,)秩不足时是空数组。❌ 以为lstsq的rcond默认值永远是-1。✅ 从 NumPy 2.0 起默认是None用机器精度乘以max(M, N)显式传-1才是机器精度跨版本时显式写出更稳。❌ 用svd返回的Vh时再对它取一次转置。✅ 第三个返回值已经是共轭转置后的Vh重构应写U np.diag(S) Vh。❌ 对实对称矩阵用np.linalg.eig并假设特征值按升序排列。✅ 对称厄米矩阵用eigh且无论用哪个都不要假设特征值的顺序。❌ 看到solve报了LinAlgError就去求伪逆硬算。✅ 报错说明a奇异或非方阵先检查矩阵秩与条件数确实需要最小二乘解时用lstsq。总结需求该用不该用矩阵乘法/np.matmul*那是逐元素乘解a x b方阵满秩np.linalg.solveinv再相乘非方阵 / 奇异矩阵求近似解np.linalg.lstsqsolve会直接报错真的需要逆矩阵本身np.linalg.inv/pinv仅为解方程而求逆对称矩阵特征值np.linalg.eighnp.linalg.eig奇异值分解np.linalg.svd注意Vh忘记第三个返回值已转置判断可解性matrix_rankcond直接算完就用把numpy.linalg用对核心不在于记住多少个函数名而在于两件事分清哪些是顶层的numpy函数、哪些属于numpy.linalg以及在数值路径上选对工具——解方程用solve、对称矩阵用eigh、可解性先看秩与条件数。做到这两点绝大多数「结果莫名其妙」的线性代数问题都会消失。参考numpy.linalg、numpy.matmul、numpy.linalg.solve、numpy.linalg.lstsq等的签名与版本说明以 NumPy 官方文档为准NumPy 为第三方库需pip install numpy本文代码未在本机运行仅作人工推演。
RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进