★ FEATURED ARTICLE 前阵子帮一位做数值模拟的朋友排查一个特别隐蔽的问题他用 NumPy 矩阵行列式来判断某个 3x3 变换矩阵是否可逆代码里写if det 0结果跑出来一个-1.8e-16判断永远不成立后面求逆直接炸出一串inf。这种坑在纯计算脚本里顶多算个小小的事故可一旦你要把结果画出来给人看、还要在汇报里解释这个矩阵为什么把平面压扁成了一条线那 NumPy 矩阵行列式计算与可视化这两件事就绕不开了——算得准不准和看得清不清是同一个问题的两面。这篇东西想聊的就是这条链路从np.linalg.det到底在底层干了什么到行列式的几何意义怎么用二维平行四边形和三维平行六面体画出来再到参数扫描时怎么用热力图一眼定位奇异区最后是批量计算、大矩阵溢出、整数精确计算这些真刀真枪的问题。适合谁看写过几行 NumPy、但被浮点精度坑过的同学要给学生讲线性代数又缺一张好图的老师以及做图形学、机器学习特征分析、结构力学小工具需要判断矩阵退没退化的工程人员。数学基础不用太深会矩阵乘法就够剩下的我用代码和图慢慢补齐。1. 先把行列式这件事想清楚再动手1.1 从一次翻车的代码评审说起那次评审的原始代码大概长这样一个旋转变换加一点缩放矩阵元素全是小数det算出来是2.9999999999999996而不是3。朋友的第一反应是NumPy 算错了第二反应是把判断条件从 0改成 1e-12然后就把 bug 藏起来了。我后来让他打印np.linalg.cond条件数已经到1e16量级这说明矩阵离奇异只差一点点改阈值只是把错误推迟到下一次。这件事给我两个教训也构成了后面所有内容的前提。第一行列式在数学上是精确的代数对象在计算机里是一个有误差的浮点数任何跟 0 的比较都必须带上下文第二行列式本身不擅长表达离奇异有多远它是个带量纲的乘积型指标元素量级一变就失去可比性真正稳定的可逆性判据是条件数和奇异值。所以动手写代码之前先在心里摆好三件事的位置行列式的数学定义用于理解、NumPy 的实现路径用于知道误差从哪来、可视化手段用于向别人解释。三者缺一块写出来的东西都容易在某个地方漏气。1.2 行列式的几何直觉它就是面积缩放倍率线性代数教材里那句行列式表示线性变换对体积的缩放倍数第一次看很难有体感。换个说法把一个 2x2 矩阵看成两根列向量它们从原点伸出去张开一个平行四边形这个平行四边形的面积就是行列式的绝对值符号则告诉你这两根向量是右手系还是左手系也就是定向翻没翻。三维里换成三根列向量张成的平行六面体体积等于行列式的绝对值。这个直觉极其好用因为它把一堆抽象性质瞬间变得显而易见。行列式为 0说明这个体被压扁了二维里压成一条线段面积为 0三维里压成一个平面或一条线体积为 0这种变换不可逆因为多个不同的输入点被映射到了同一个输出点信息丢了。行列式为负说明变换把整个空间翻了个面就像把一张纸翻过来看。行列式的绝对值大于 1空间被拉伸小于 1被压缩等于 1体积守恒——这也是为什么纯旋转矩阵的行列式恒为 1它只是转不变形。再补一个常用性质后面写代码时天天用det(A B) det(A) * det(B)。这意味着如果你先做一次变换再做一次变换总的体积缩放倍率就是两者相乘。反过来det(A^-1) 1 / det(A)所以行列式越小逆矩阵的元素整体越大数值上越不稳定。这条性质直接解释了为什么接近奇异比完全奇异更危险完全奇异求逆会直接报错或者给inf你一眼就发现了接近奇异求逆会给你一堆看起来人模人样、实际误差爆炸的数字这才是最阴的。1.3 为什么是 NumPy而不是自己写循环有人会问行列式公式不就那点东西吗2x2 是ad - bc3x3 是六个乘积相加减为什么还要用库小矩阵确实可以手写而且手写有时更快后面第 4 节会给数据。但一旦矩阵维度到 5 以上按定义展开的复杂度是阶乘级的n10 就要算三百六十多万项纯属自残。NumPy 底下挂的是 LAPACK用的是 LU 分解路线复杂度 O(n³)而且经过几十年的数值优化缓存友好、向量化、分块都做透了。更关键的是生态衔接。你算完行列式多半要画图画图要用 matplotlib而 matplotlib 吃的就是 NumPy 数组。这条链路是通的np.array造矩阵 →np.linalg.det算值 →matplotlib.patches.Polygon画形状 →annotate标数值中间不需要任何格式转换。如果换成手写 C 扩展或者别的数学库你还得处理数据搬运和类型映射光调试这些就够喝一壶。依据我自己的习惯日常的组合是numpy负责计算matplotlib负责静态图需要交互时用plotly或者ipywidgets挂个滑块实时改矩阵元素。这套组合的学习成本很低但能把抽象概念讲得特别直观尤其适合给别人做科普或者给自己做检查。2. NumPy 里的 det 到底是怎么算出来的2.1 打开黑箱LU 分解与误差的来源np.linalg.det的底层路径可以概括成一句话对矩阵做带部分主元选取的 LU 分解A P·L·U其中 P 是置换矩阵L 是单位下三角U 是上三角然后利用行列式的乘法性质det(A) det(P)·det(L)·det(U)。det(L)因为对角线全是 1恒等于 1det(U)就是对角线元素连乘det(P)是 ±1取决于做了奇数次还是偶数次行交换。所以整个计算退化成做一次 LU取 U 的对角线连乘再根据换行次数定符号。这个路径解释了误差从哪来。连乘会累积相对误差每做一次乘法相对误差大致按机器精度量级往上叠。float64 的机器精度约2.2e-16n 阶矩阵要连乘 n 次最坏情况下相对误差的量级在n·eps附近一般远小于这个上界但绝对不可能为 0。这就是为什么你算det(np.eye(3))能得到漂亮的1.0而算一个元素全是小数的 3x3 矩阵结果可能是2.9999999999999996。它不是 bug是浮点运算的本性。还有一个细节值得知道主元选取会引入行交换这会影响符号。NumPy 内部帮你处理了但如果哪天你自己实现 LU 去验证忘了跟踪交换次数算出来的行列式符号会莫名其妙反掉然后你会以为是数学错了。我当年干过这事对着一个 3x3 矩阵验算了半小时最后发现是换行没记符号。2.2 det、slogdet、cond 三个函数各自的分工这三个函数经常被混用其实定位完全不同我列个表说清楚。函数返回值适用场景主要风险np.linalg.det(A)单个浮点标量或复数中小矩阵、量级正常、只关心数值大矩阵上溢/下溢接近 0 时精度不可靠np.linalg.slogdet(A)(sign, logabsdet)二元组大矩阵、行列式量级跨度极大、需要连乘比较仍需自己判断 sign 是否为 0np.linalg.cond(A)条件数正实数判断可逆性、评估数值稳定性大矩阵计算开销高于 detnp.linalg.matrix_rank(A)整数秩判断是否满秩、奇异检测依赖容差参数需要自己设定slogdet是我最推荐给工程场景的一个。它返回符号和绝对值的自然对数你拿sign * np.exp(logabsdet)就能还原行列式。好处是大矩阵也不会溢出假设一个 500x500 的矩阵元素量级在 2 附近行列式的真实值可能在2^500这个量级float64最大才约1.8e3082^500早就飞出去了det直接返回inf但log|det|大概是500·ln2 ≈ 346舒舒服服在范围内。这个差别在实际项目里非常致命尤其是你在做高斯过程、贝叶斯模型或者任何涉及多元概率密度计算的活儿密度函数里那项1/sqrt(det(Σ))一旦溢出整个似然就变成nan而且你很难从报错信息里看出是行列式炸的。cond的价值在于它没有量纲。行列式的数值大小受矩阵元素量级影响极大把整个矩阵乘 10n 阶矩阵的行列式就乘10^n但矩阵的可逆性一点没变。条件数则不受这个影响它只反映最大奇异值比最小奇异值大多少接近 1 说明矩阵性质很好到1e15就基本可以认为在 float64 下不可逆了。所以判断这个矩阵能不能求逆我的习惯是先看cond再看rankdet只用来做展示和几何解释。2.3 那些容易忽略的参数和返回值陷阱说几个我自己踩过或者见别人踩过的点都很琐碎但每一个都能让人卡半天。第一个整数矩阵传入det后返回的是浮点数。np.linalg.det(np.array([[1,2],[3,4]]))返回-2.0000000000000004这样的浮点值不是整数-2。如果你后面拿它做整除、取模、或者当字典 key行为会跟预期不同。想要精确整数结果得走符号计算路线第 4.3 节会讲。第二个非方阵直接抛异常。np.linalg.LinAlgError: Last 2 dimensions of the array must be square。这个好理解但错误信息里last 2 dimensions提示了一个重要信息det支持批量输入只要最后两维是方阵就行。这个特性在第 4 节会大用特用。第三个空矩阵的行为。np.linalg.det(np.empty((0,0)))返回1.0。这是数学上的约定零阶行列式等于 1逻辑上说得通但第一次遇到会愣一下。如果你在写递归或者动态规划用 0x0 当边界条件记得这个 1.0 是故意的。第四个复数矩阵返回复数。det对复数矩阵照样工作返回值是complex128。判断它是否为 0 的时候不能直接用 0得看模长。而且复数的接近 0更需要小心因为实部和虚部各自都可能带误差。第五个float32 的精度陷阱。如果你的数组是从某个深度学习框架里拿出来的float32直接丢给det误差会被放大得非常明显小矩阵就可能出现符号错误。稳妥的做法是进det之前先.astype(np.float64)多占一点内存换来十位以上有效数字的差距完全值得。提示在写任何涉及行列式是否为 0的判断之前先问自己三个问题——矩阵元素的量级是多少矩阵维度是多少我真正想判断的是能否求逆还是体积是否为零三个答案会直接决定你该用det、slogdet还是cond。3. 把行列式画出来从平行四边形到热力图3.1 二维平行四边形面积是最直观的教具二维情形是最适合入门的因为一切都看得见。给定一个 2x2 矩阵它的两列就是两根向量从原点张开的平行四边形面积等于|det|。画图的关键是三件事画原始单位方格作参照、画变换后的平行四边形、标出数值。set_aspect(equal)这行不能省否则坐标轴比例不一致你画出来的平行四边形视觉面积和det对不上会误导观众。import numpy as np import matplotlib.pyplot as plt from matplotlib.patches import Polygon A np.array([[1.5, 0.8], [0.3, 1.2]]) det np.linalg.det(A) u, v A[:, 0], A[:, 1] fig, ax plt.subplots(figsize(6, 6)) # 原始单位方格虚线作参照 unit Polygon([[0, 0], [1, 0], [1, 1], [0, 1]], closedTrue, fillFalse, ls--, ecgray, lw1) ax.add_patch(unit) # 变换后的平行四边形 para Polygon([[0, 0], u, u v, v], closedTrue, alpha0.35, fctab:blue, ectab:blue, lw1.5) ax.add_patch(para) # 两根列向量 ax.annotate(, xyu, xytext(0, 0), arrowpropsdict(arrowstyle-, colortab:red, lw2)) ax.annotate(, xyv, xytext(0, 0), arrowpropsdict(arrowstyle-, colortab:green, lw2)) ax.axhline(0, ck, lw0.5) ax.axvline(0, ck, lw0.5) ax.set_xlim(-2.5, 3.0) ax.set_ylim(-2.5, 3.0) ax.set_aspect(equal) ax.set_title(fdet {det:.4f} |area| {abs(det):.4f}) plt.tight_layout() plt.show()这段代码跑出来的效果很直接浅蓝色的平行四边形面积就是det的绝对值红色和绿色箭头分别对应矩阵的两列。把A改成[[1, 1], [1, 1]]平行四边形立刻退化成一个点加两条重合的线段面积归零det也归零那种压扁的感觉一下就出来了。我习惯在这个图上再加两个小变化对教学特别有用。一个是在det为负的时候把平行四边形换成另一种颜色并且给个翻转的标注让人意识到定向变了另一个是把单位方格的四个顶点也做一次变换用散点画上去这样能看出方格内部的点被搬到了哪里。3.2 三维平行六面体与六个面三维的画法在思路上完全一致只是平行六面体没有现成的Polygon可用得自己拼六个四边形面。三根列向量a, b, c从原点出发八个顶点分别是0、a、b、c、ab、ac、bc、abc六个面就是这些顶点的组合。用Poly3DCollection一次性塞进去透明度调到 0.35 左右既能看出体积感又不至于互相遮挡。from mpl_toolkits.mplot3d.art3d import Poly3DCollection A np.array([[1.2, 0.4, 0.1], [0.0, 1.5, 0.6], [0.2, 0.0, 0.9]]) det np.linalg.det(A) a, b, c A[:, 0], A[:, 1], A[:, 2] o np.zeros(3) faces [ [o, a, a b, b], [o, b, b c, c], [o, a, a c, c], [a, a b, a b c, a c], [b, a b, a b c, b c], [c, a c, a b c, b c], ] fig plt.figure(figsize(7, 6)) ax fig.add_subplot(111, projection3d) ax.add_collection3d( Poly3DCollection(faces, alpha0.3, facecolortab:blue, edgecolork, lw0.6) ) for vec, col in zip([a, b, c], [tab:red, tab:green, tab:orange]): ax.quiver(0, 0, 0, *vec, colorcol, lw2, arrow_length_ratio0.1) ax.set_xlim(0, 2.2) ax.set_ylim(0, 2.2) ax.set_zlim(0, 2.2) ax.set_xlabel(x) ax.set_ylabel(y) ax.set_zlabel(z) ax.set_title(fVolume |det| {abs(det):.4f}) plt.show()三维图有个老生常谈的问题默认视角经常把某个面挡得严严实实看起来像一片糊。解决办法是调整ax.view_init(elev22, azim35)这两个参数多试几次找到能看清三条棱的角度。另一个技巧是把六个面的透明度设成不同值让靠近观察者的面更不透明远处的更透人眼会自动补出立体结构。这些不是数学问题是纯粹的视觉工程但做得好不好直接决定你的图能不能用。3.3 让线性变换动起来插值动画静态图能说明一个矩阵的效果动画能说明一族矩阵的变化趋势。做法很朴素在单位矩阵I和目标矩阵A之间做线性插值M(t) (1-t)I tAt 从 0 走到 1每一帧把M(t)作用到一组网格线上同时实时更新det的显示。你会看到网格先被拉伸、旋转然后在某个 t 上突然塌陷——那个瞬间就是行列式过零的时刻。from matplotlib.animation import FuncAnimation A np.array([[1.8, 1.0], [0.5, 0.6]]) frames, n_line, n_pt 60, 8, 80 # 预先生成单位方格内的网格线每条线用一串采样点表示 s np.linspace(-1, 1, n_pt) segs [] for c in np.linspace(-1, 1, n_line): segs.append(np.stack([s, np.full_like(s, c)])) # 水平线 segs.append(np.stack([np.full_like(s, c), s])) # 竖直线 segs np.array(segs) # (2*n_line, n_pt, 2) fig, ax plt.subplots(figsize(5.5, 5.5)) ax.set_xlim(-2.5, 2.5) ax.set_ylim(-2.5, 2.5) ax.set_aspect(equal) lines [ax.plot([], [], lw0.8, colorsteelblue)[0] for _ in segs] title ax.set_title() def update(i): t i / (frames - 1) M (1 - t) * np.eye(2) t * A d np.linalg.det(M) transformed segs M.T # 一次性变换所有点 for ln, pts in zip(lines, transformed): ln.set_data(pts[:, 0], pts[:, 1]) title.set_text(ft {t:.2f} det {d:.3f}) return lines [title] anim FuncAnimation(fig, update, framesframes, interval60, blitFalse) plt.show() # 需要保存时anim.save(transform.gif, fps20)这里有个性能细节值得说不要在内层循环里对每个点单独做矩阵乘法而是把所有采样点堆成一个(K, 2)数组一次segs M.T全搞定。n_line8、n_pt80的时候一共一千多个点逐点乘和批量乘的差距不明显但如果你想画 50 条网格线、每条 500 个点逐点乘会让帧率掉到个位数批量乘依然流畅。这个习惯在后续做参数扫描时更重要。blitFalse是我保守的选择。blitTrue理论上更快但一旦你的图里有set_title、colorbar或者坐标轴标签变化很容易出现文字不刷新的鬼影。为了少踩坑我一律关掉除非帧率实在撑不住。3.4 参数扫描热力图一眼找到奇异面如果矩阵的形式是带参数的比如[[1, x], [y, 1]]那它的行列式是1 - xy等于 0 的地方构成双曲线xy 1。把x、y在[-3, 3]上扫一遍用等高线图把det画出来再把det 0那条线加粗描出来奇异区域的位置一目了然。这种图在调参时非常好用你能看到安全区和危险区的分布形状而不是盲目试数。x np.linspace(-3, 3, 401) y np.linspace(-3, 3, 401) X, Y np.meshgrid(x, y) Z 1 - X * Y # det([[1, x],[y, 1]]) fig, ax plt.subplots(figsize(6.2, 5.2)) cf ax.contourf(X, Y, Z, levels40, cmapRdBu_r) ax.contour(X, Y, Z, levels[0], colorsk, linewidths2.2) ax.contour(X, Y, Z, levels[1], colorsk, linestyles:, linewidths1.2) fig.colorbar(cf, axax, labeldet) ax.set_xlabel(x) ax.set_ylabel(y) ax.set_aspect(equal) ax.set_title(det 1 - x*y黑实线为 det0) plt.show()这类图的读法需要一点训练。实线det 0两侧颜色对比最剧烈说明行列式在这附近变化快也意味着矩阵在这个区域对参数极其敏感虚线是det 1代表面积守恒的那条线红蓝分界说明符号翻转变换的定向在这里改变。我常用这张图给做控制系统的人解释为什么参数不能随便取边界值比讲一堆特征值理论有效得多。如果想换成三维曲面看直接把contourf换成plot_surface再加个contour投影到底面效果更震撼但读数值不如二维等高线准。我的建议是二维等高线用来定量读数三维曲面用来做展示两者别混用。4. 批量、精度与性能实战里的硬骨头4.1 批量计算能堆叠就别写循环np.linalg.det支持最后两维是方阵的高维数组这意味着你可以一次性算十万个 2x2 矩阵的行列式而不需要 Python 层的 for 循环。这个特性在有大量小矩阵的场景里是救命稻草比如图形学里每个三角面片的变换、机器学习里每个样本的协方差矩阵、或者参数扫描时每个网格点对应的矩阵。rng np.random.default_rng(0) N 200_000 A rng.standard_normal((N, 2, 2)) # 方式一批量调用推荐 d1 np.linalg.det(A) # 方式二显式公式避免了 LAPACK 调用开销 a, b A[:, 0, 0], A[:, 0, 1] c, d A[:, 1, 0], A[:, 1, 1] d2 a * d - b * c print(np.allclose(d1, d2)) # True我在这台机器上实测的耗时大致是这个量级批量det走 LAPACK200k 个 2x2 大概 40 到 60 毫秒显式公式只要 3 到 5 毫秒快一个数量级。原因不复杂det虽然向量化了但每个矩阵还是要走一遍完整的 LU 流程包括主元选取和符号处理而 2x2 的公式就是四次索引加两次乘一次减几乎没有额外开销。所以经验是维度 ≤ 3 且调用量大时手写公式维度更大或者追求代码简洁时用det。3x3 的显式公式稍微长一点但依然值得写def det3_batch(A): a, b, c A[:, 0, 0], A[:, 0, 1], A[:, 0, 2] d, e, f A[:, 1, 0], A[:, 1, 1], A[:, 1, 2] g, h, i A[:, 2, 0], A[:, 2, 1], A[:, 2, 2] return (a * (e * i - f * h) - b * (d * i - f * g) c * (d * h - e * g))不要小看这点优化。我做过一个粒子模拟的小工具每帧要处理三万多个 3x3 变换用np.linalg.det的时候帧率卡在 22 左右换成显式公式后直接上到 55。数值结果完全一致在1e-12容差内代码也只是多了五行。4.2 大矩阵防溢出slogdet 的正确用法前面提过slogdet这里给一个具体的量化例子。取一个 300x300 的对角矩阵对角线元素全为 3真实行列式是3^300这个数有多大3^300 ≈ 10^143还在 float64 范围内上限约1.8e308。但如果元素是 10 呢10^300就是1e300勉强在范围内再大一点就inf了。而如果矩阵是 400 阶、元素量级 55^400 ≈ 10^280也悬。真实的协方差矩阵维度经常上千元素量级不大但连乘起来很容易越界。n 800 A np.eye(n) * 2.0 np.random.default_rng(1).standard_normal((n, n)) * 0.01 sign, logabsdet np.linalg.slogdet(A) print(fsign {sign}, log|det| {logabsdet:.4f}) # 需要真实值时注意可能溢出只在确实需要时才还原 try: det_val sign * np.exp(logabsdet) print(fdet {det_val:.4e}) except OverflowError: print(数值过大保留对数形式)这里有一个我踩过的细节np.exp对超大值会返回inf并给一个运行时警告而不是抛OverflowError。所以try/except其实是抓不住的得用np.errstate配合np.isinf来检查。更稳妥的写法是用np.exp(np.clip(logabsdet, -700, 700))把指数限幅在安全区间宁可牺牲精度也不要拿到inf污染下游计算。另一个常见需求是比较两个矩阵的行列式谁大。直接比较det会溢出比较logabsdet就稳了如果两个符号不同那比较的其实是带符号的大小得先判断 sign 再比绝对值。这种逻辑写起来有点啰嗦但一旦写进工具函数后面复用起来很舒服。4.3 整数矩阵想要精确结果就得换工具NumPy 的det一律走浮点哪怕输入是int64。这对纯理论计算或者需要精确判等的场景是个问题。比如判断一个整数矩阵是否奇异浮点算出来1e-17你说它是 0 还是不是 0数学上它可能就是个非零的整数。第一种替代方案是用fractions.Fraction配合分数版高斯消元全程精确有理运算。小矩阵10 阶以内完全够用代码也就二十来行。from fractions import Fraction def det_exact(mat): n len(mat) M [[Fraction(x) for x in row] for row in mat] det Fraction(1) for i in range(n): piv next((r for r in range(i, n) if M[r][i] ! 0), None) if piv is None: return Fraction(0) if piv ! i: M[i], M[piv] M[piv], M[i] det -det det * M[i][i] for r in range(i 1, n): factor M[r][i] / M[i][i] if factor 0: continue for c in range(i, n): M[r][c] - factor * M[i][c] return det print(det_exact([[1, 2], [3, 4]])) # -2 print(det_exact([[1, 2], [2, 4]])) # 0第二种方案是引入符号计算库直接调用它的行列式函数。优点是功能全、支持符号变量缺点是有额外依赖、大矩阵慢得离谱、返回值类型不是 NumPy 数组和后面的画图流程要额外转换。我的判断标准很简单如果只是要个精确整数答案用分数版如果要把符号变量也带进去比如推导含参矩阵的行列式表达式那再上符号库。第三种方案适合固定小维度把 2x2、3x3 的整数行列式公式写出来用 Python 原生大整数算永远不会溢出、永远精确。a*d - b*c这种表达式在整数上是封闭的结果一定是精确整数。这个方法适合嵌入式或者需要极致确定性的场景。注意不要用round(np.linalg.det(A))来恢复整数结果。矩阵维度稍大一点浮点误差就可能大到1e-3以上四舍五入会给你一个错误的整数而且错得很隐蔽。整数精确计算没有捷径只能换算法。5. 常见问题与排查速查表5.1 精度与判零类问题这类问题占了实际求助的一大半。典型现象是明明数学上是 0代码算出来不是 0或者两个应该相等的行列式算出来差了一点。根因无非三种浮点累积误差、元素量级差异过大导致的小值被吃掉、以及用了float32。排查顺序我一般这样走先打印矩阵的条件数如果cond 1e12那这个矩阵本身就病态任何浮点算法都救不了得从建模层面改再看元素量级如果最大元素和最小元素差十几个数量级先做归一化再算最后检查 dtype确保进det之前是float64。判断近似为零时绝对阈值1e-12只在元素量级为 1 左右时成立。更稳的做法是相对判据abs(det) tol * np.prod(np.linalg.norm(A, axis0))也就是拿列范数的乘积当尺度参考。这个尺度在几何上就是如果这些列正交体积该有多大用它做分母判据才有量纲一致性。5.2 奇异检测与可逆性判断很多人用det 0判断不可逆这在教学里没问题在工程里会出事。正确的优先级是cond优先rank次之det只做辅助说明。对于方阵cond无穷大意味着严格奇异有限但极大意味着数值上不可逆。np.linalg.cond默认用 2-范数计算代价是 SVD对 1000 阶矩阵大概几百毫秒能接受。如果只想要个快速筛查可以先用np.linalg.det看个大概再用cond确认两段式判断能省不少时间。matrix_rank的问题在于它内部也用 SVD并且有一个容差参数默认是max(M, N) * eps * sigma_max。这个默认值在大多数情况下合理但对于元素量级极端的矩阵可能不合适需要手动传tol。我遇到过一个小工具输入矩阵元素全是1e-8量级matrix_rank判定秩为 0实际上满秩就是因为容差相对元素量级太大了。归一化之后一切正常。5.3 可视化环节的坑画图这块的问题不影响数值正确性但影响你能否说服别人所以也值得列一列。第一个set_aspect(equal)忘了加二维平行四边形看起来是歪的面积比例全乱。这个错误极其常见因为默认的坐标轴比例会把图自动拉伸以填满画布。第二个三维图的遮挡。Poly3DCollection的绘制顺序依赖深度透明面叠在一起的时候会出现该在前面的跑到后面的错觉。缓解办法是减少透明度层次、调高alpha的对比度、或者干脆只画面不画边。第三个中文标签显示成方框。这是字体配置问题需要在matplotlib.rcParams里指定一个支持中文的字体族并且把axes.unicode_minus设为False否则负号也会出问题。这两个设置建议放在脚本最开头统一做一次。第四个动画保存失败。anim.save依赖外部编码器机器上没有的话会报一个含糊的错误。开发阶段先用plt.show()看效果确认逻辑对了再折腾保存。第五个热力图的levels设太少会让等高线变形成诡异的块状设太多又会让渲染变慢。40 到 60 之间是个舒服的区间contourf和contour用同一组levels能保证颜色和线条对齐。现象可能原因处理办法det返回inf或nan大矩阵连乘溢出、矩阵含nan改用slogdet先np.isnan(A).any()检查行列式应该为 0 却得到小值浮点误差、矩阵病态改用cond判断归一化后重算批量计算报维度错误最后两维不是方阵检查A.shape[-2:]是否相等符号与理论值相反手写 LU 未跟踪行交换记录交换次数奇数次取负结果显示位数为 1数组 dtype 是float32转float64二维图形状看起来不对缺set_aspect(equal)补上该调用动画卡顿逐点矩阵乘法批量segs M.T中文显示为方块字体未配置设置字体族并关闭 unicode_minus6. 我自己的几个使用习惯最后分享几个用久了之后固化下来的做法都是小细节但省了很多来回。第一个习惯是先建一个矩阵工厂函数。不管是教学还是调试我总会准备一个make_mat(kind, **kw)之类的小函数能快速产出旋转、缩放、剪切、近奇异、严格奇异这几类矩阵。这样画图脚本里就不用手写一堆np.array([[...]])改参数也方便。尤其是画对比图的时候从同一个工厂拿两个只差一点的矩阵视觉效果最能说明问题。第二个习惯是把算和画彻底分开。计算部分写成纯函数输入矩阵输出数值不做任何绘图副作用绘图部分只吃数值不碰np.linalg。这样当数值出问题时可以单独在 REPL 里验证计算函数不受图形环境干扰反过来当图不好看时也不用担心改动影响数值。这个分离看起来是常识但在赶进度的时候特别容易被破坏一旦混在一起调试成本会翻倍。第三个习惯是给每张图都标上数值。人眼对面积的判断精度远不如对数字的判断一张没有标注det的平行四边形图观众只能看出大概是个平行四边形看不出它到底是 3.2 还是 0.3。把数值直接写进标题或者图内注释图的解释力立刻上一个档次。第四个习惯是奇异判定永远留一条退路。任何依赖矩阵可逆的代码路径我都会额外准备一个分支条件数超阈值时走伪逆或者降维而不是硬求逆然后等着inf冒出来。这条习惯来自一次线上事故当时一个协方差矩阵在某个数据分布下接近奇异求逆出来的结果后面被当成正常值用了三个小时最后是一堆完全不合常理的业务指标才让人发现不对。第五个习惯大矩阵的行列式只在对数域里用。除了展示给非技术同事看的时候需要还原成真实值其他所有内部计算我都停在logabsdet上。概率密度、模型比较、高斯过程的边际似然这些场景全都能在对数域里完成而且数值稳定得多。
阅读完成 · 觉得有帮助？
