在 MATLAB 里判断一个矩阵是否正定看起来就是一行chol(A)或者all(eig(A)0)的事但真正落到工程里我吃过的亏、踩过的坑还真不少。做优化算法要把 Hessian 矩阵拿来判断收敛性做统计得确认样本协方差矩阵是不是满秩正定的做有限元计算要验证刚度矩阵性质——正定性几乎贯穿了所有数值计算场景。很多教程只会告诉你“特征值大于零就是正定”可实际跑起来问题就来了矩阵不对称怎么算特征值小到 1e-15 到底算正还是算负为什么chol明明报错了数学上它却应该是正定的这篇文章我从定义、判据、MATLAB 例程到常见坑点一次讲清楚尽量让新手看完能直接拿去用也让写过不少代码的同行能补充一些细节认知。1. 判断矩阵正定的数学基础我们得先回到定义本身去理解这几个判据否则用起chol、eig这些函数时你只知道“能判断”不知道“为什么它这么判断”遇到例外情况很容易懵。1.1 正定矩阵的定义与几何意义一个 n 阶实对称矩阵 A 被称为正定矩阵指的是对于任意非零向量 x ∈ R^n二次型 x^T A x 始终大于 0即x^T A x 0∀ x ≠ 0我最早学这个定义的时候不太有感觉后来把它和二次曲线、等高线联系起来才觉得形象。你可以把 x^T A x 想象成一个曲面在向量 x 方向的“高度”正定就意味着这个曲面在所有方向上都开口向上原点是个全局唯一的极小点。这个性质放到优化里太重要了——如果目标函数的 Hessian 矩阵正定二阶必要条件就保证你站在一个极小点附近算法往哪个方向走都“有底”。放到统计里协方差矩阵正定意味着数据在每个方向上都存在真实的方差没有被完全共线的冗余变量弄退化放到有限元里刚度矩阵正定意味着系统在约束下是稳定的没有刚体位移或者机构失稳。所以我常说判断矩阵正定不是做数学题是在给整个数值算法的稳定性做体检。这里有一个非常容易被忽略的点为什么定义里默认说“实对称矩阵”如果你拿到一个非对称矩阵 B也想去算 x^T B x那实际上可以把 B 拆成对称部分和反对称部分B (B B^T)/2 (B - B^T)/2后一项是反对称矩阵任意反对称矩阵 K 都满足 x^T K x 0因为每一项 k_ij · x_i · x_j 和它的对称镜像恰好抵消了。所以整个二次型 x^T B x 完全由 B 的对称部分 (BB^T)/2 决定。这个结论我在后文“避坑”部分还会反复用到因为实际工程拿来的矩阵经常不是对称的但你用chol之前必须先把这个问题处理掉否则结果可能根本没有意义。1.2 三大等价判据特征值、顺序主子式、Cholesky判断正定数学上我有三个认知路径它们互相等价但在计算机上的表现却完全不同。第一是特征值判据。对于实对称矩阵 A它可以正交对角化即存在正交矩阵 Q 使得 Q^T A Q diag(λ_1, ..., λ_n)。于是 x^T A x Σ λ_i y_i^2要让这个值对所有非零 x 恒正必须每个 λ_i 都大于 0。这个方法最直观缺点是要计算特征分解矩阵大一点时间开销就上来了。第二是顺序主子式判据也叫 Sylvester 判据。它的表述是A 正定当且仅当所有顺序主子式 det(A_k) 都大于 0其中 A_k 是 A 的前 k 行前 k 列组成的子矩阵。这个判据证明起来可以用归纳法和 Schur 补核心逻辑就是“每次消掉前 k-1 个变量后剩下的末维系数仍然为正”。缺点也很明显它需要用行列式而行列式在计算机里是个数值上很脆弱的量大矩阵算起来又慢又不稳定。第三是 Cholesky 分解判据。实对称正定矩阵可以唯一分解为 A R^T R其中 R 是一个对角元全为正的上三角矩阵。这个分解可以看作是高斯消去法里“把每个主元再开一次平方根”的矩阵版本主元大于零开平方根才有意义只要某个主元变成 0 或者负数平方根就取不出来矩阵就一定不是正定。MATLAB 里的chol函数就是干这件事的。除此之外还有一个主元判据对称矩阵正定当且仅当不进行行列交换的高斯消去过程中所有主元都大于 0。这个判据其实跟 Cholesky 是同一件事的不同表述数学上等价数值上也同样依赖主元是否“足够正”。1.3 理论等价计算机上却要分情况选择三个判据在精确算术下等价但浮点数世界里完全不是一回事。你直接用det去判断顺序主子式小矩阵还行矩阵一超过 5 阶行列式的数值范围可能变得极其夸张接近 0 的行列式加上浮点误差非常容易误判。特征值分解本质上是迭代法要反复做 QR 变换、Householder 约化所以它最稳健但也最慢。Cholesky 分解是直接法只需 O(n^3) 的加减乘除和开方速度快得多但遇到接近奇异的矩阵时会因为主元接近 0 而失败这也算是一种“宁可错杀也不放过”的保守判断。所以工程上的选择逻辑很清晰只想知道是/否用chol想知道矩阵离正定有多远、特征值分布长什么样用eig给学生讲概念才用顺序主子式。2. MATLAB 实现正定性判断的三种方法这一节是核心实操内容。我分别把chol、eig、顺序主子式三种方法在 MATLAB 里怎么写、输出怎么解读、各自的坑在哪都过一遍。2.1 方法一基于 Cholesky 分解的 chol 方案chol是 MATLAB 中判断正定最常用、也是我自己的默认首选。基础语法有两种。一种是直接调用A [4 2; 2 3]; R chol(A);如果 A 是正定的R 是一个上三角矩阵满足 R*R A也就是 Cholesky 因子。如果 A 不是正定直接这样写会报错中断脚本。所以日常我更习惯用第二种带 p 的写法A [4 2; 2 3]; [R, p] chol(A); if p 0 disp(矩阵正定); disp(R); else disp([矩阵不正定分解失败位置在第 , num2str(p), 行]); disp(R); end这里的[R, p] chol(A)有两个关键点。第一当 p 0 时表示分解成功可以确认矩阵正定当 p 0 时表示分解在第 p 步停止前 p-1 个对角线元素都做完了但从第 p 行开始主元不再大于 0。第二即使失败R 中已经成功求得的部分依然会返回这很有用——你可以从 R 的规模大致看出矩阵在哪个维度上开始“撑不住了”。我拿一个最经典的半正定矩阵演示A [1 1; 1 1]。它的特征值是 0 和 2行列式为 0是半正定但不是正定。执行[R, p] chol(A)后p 会等于 2。原因很简单第一步主元是 sqrt(1) 1R(1,1) 1第二步计算对角线元素 R(2,2) 时要对 A(2,2) - R(1,2)^2 开平方这个值是 1 - 1 0平方根开不出来只能在 p 2 处停住。这个例子我建议所有初学者都亲手跑一遍它能让你真正理解“Cholesky 前进一步背后就是一次高斯消去主元”这件事。实际运算时还要注意chol默认只用矩阵的上三角部分做计算假定下三角是它的对称副本。这背后的逻辑是既然理论上 Cholesky 只适用于对称矩阵那就不再看下三角省点计算量。可如果你的 A 在数值上并不对称chol算的其实是一个“假装对称”后的矩阵这就会造成误导。所以我在正式代码里一律先做对称化后面第 3 章的例程你会看到这个习惯。2.2 方法二基于特征值的 eig 方案特征值判据适合两种场景一是你想知道特征值的具体分布而不仅仅是一个是/否的结论二是你要同时判断半正定、负定、不定这些情况。核心代码也很短A [4 2; 2 3]; d eig(A); n size(A, 1); tol n * eps(max(abs(d))); if all(d tol) disp(矩阵正定); elseif all(d -tol) disp(矩阵半正定); elseif all(d -tol) disp(矩阵负定); else disp(矩阵不定); end注意我加了容差 tol这是和“裸写all(d 0)”最本质的区别。因为浮点数计算出来的特征值不可能是精确的 0一个理论上为 0 的特征值在计算机里可能变成 1e-18 或者 -1e-18。如果你用all(d 0)去判断一个数学上半正定的矩阵到你这里就会被判成“不定”运气好是负定反正不正定但它其实只是缺了那么一丝精度。容差的常用近似取法是 n 乘以最大特征值对应的浮点精度eps(max(abs(d)))这背后是“误差随矩阵维度和量级增长”的经验公式。它不是一个严格定理但工程上足够好使。再说eig的代价。特征值的计算不是直接法它要通过反复迭代把矩阵约化成对角形时间复杂度和chol差不多也是 O(n^3)但常数项大得多。一个 1000 阶的矩阵chol可能几毫秒就完事eig就要花掉几十上百毫秒。所以如果只是单纯的“是不是正定”我会用chol如果矩阵本身是要做奇异值分析、主成分分析、谱聚类这类需要特征谱的任务我才会顺手用eig的结果顺带判断正定性这样不浪费计算量。eig对对称性的容忍度其实比chol更宽松它默认对非对称矩阵也照样算特征值只是返回的特征值可能是复数。如果d里出现复数all(d tol)就会毫不犹豫地返回 false因为 MATLAB 里复数和实数比大小会直接报错。这其实是个变相提醒你的矩阵可能不是对称的。我建议判断之前先用issymmetric(A)看一眼或者干脆执行对称化。2.3 方法三基于顺序主子式的 Sylvester 判据顺序主子式判据在数学上是完整的但它真的只适合在教学和小矩阵演示中使用。我用 MATLAB 实现一个最直白的版本A [4 2; 2 3]; n size(A, 1); flag true; th n * eps; for k 1:n if det(A(1:k, 1:k)) th flag false; break; end end if flag disp(所有顺序主子式均大于0矩阵正定); else disp([矩阵不正定第一个不满足的主子式阶数: , num2str(k)]); end这段代码的原理不难依次取 A 的前 1 阶、前 2 阶……前 n 阶子矩阵求行列式任何一个不大于 0 就立刻终止。比如 A [4 2; 2 3]一阶主子式 det(4) 4二阶主子式 det(A) 43 - 22 8都为正所以正定。为什么我不推荐实际用它做工程判断第一det的数值稳定性在浮点世界里非常让人头疼。行列式本质上是一个连乘量n 阶矩阵的行列式数值范围可能是天文数字也可能小到超出浮点表示范围。即使理论上为正算出来也可能因为溢出或者舍入误差变成 0 甚至负数。第二每求一个顺序主子式都要对那个子矩阵做一次新的分解相当于从 1 阶到 n 阶一共做了 n 次矩阵分解复杂度急剧上升。曾经我拿一个 50 阶矩阵跑顺序主子式判据卡了很长时间换成chol一瞬间就出了结果。所以这个方法我现在的定位是“讲课时演示逻辑用”生产级代码永远不要这么写。2.4 三种方法对比与使用建议我把三种方法的核心特征整理成一张表方便你直接对照选择判断方法核心函数复杂度数值稳健性主要适用场景Cholesky 分解cholO(n^3)常数小速度快中等接近奇异时可能误判多数工程判断判断仅正定/不正定特征值eigO(n^3)常数大较慢较好配套容差后很稳定需要特征谱、半正定/负定/不定分类顺序主子式detO(n^3)·n 次分解慢差浮点下容易误判教学演示、极小规模矩阵我的选择逻辑可以总结成一句话如果你的问题只问“这个矩阵正不正定”并且你接下来还想解一个 Ax b 或者算一个二次型最小化那就直接chol因为它分解出来的 R 还能继续用来解方程一石二鸟如果你在怀疑矩阵为什么不正定、想知道是不是有某个方向上的方差几乎为 0那就eig看最小特征值到底小到什么程度比任何标志位都直观顺序主子式就当成一个“理论存在”的参照物吧看看就好。3. 可直接运行的完整例程与实战演练光讲函数用法还不够我把判断逻辑封装成一个可复用的函数并给出两个实际工程项目里会遇到的调用场景。这段代码我直接放到自己的工具库里反复用遇到问题也更容易定位。3.1 封装一个可复用的 isPositiveDefinite 函数我写这个函数时做了几个我认为很重要的设计。第一是默认方法选chol因为工程场景大多数只需要一个快速结论第二是支持eig和minor顺序主子式两种切换用于诊断和数据探查第三是整个函数内部先做对称化处理避免矩阵不对称带来的隐性错误第四是输出结构里带上对称化后的矩阵和特征值如果算过这样调试时能多看一眼。function [isPos, res] isPositiveDefinite(A, method, tol) % 判断实数方阵是否正定 % 输入: % A - 实数方阵 % method - chol默认| eig | minor % tol - 判断阈值默认根据数值精度自动计算 % 输出: % isPos - 逻辑值true表示正定 % res - 结构体含对称化矩阵、方法、可选特征值等 if nargin 2 || isempty(method) method chol; end if nargin 3 || isempty(tol) tol []; end % 方阵检查 if ~ismatrix(A) || size(A, 1) ~ size(A, 2) error(输入必须是方阵); end % 强制对称化并检查原始矩阵的对称偏差 Asym (A A.) / 2; skewPart norm(A - Asym, fro); if skewPart 1e-10 warning(原始矩阵不对称对称部分偏差范数为 %.4e将使用对称部分判断, skewPart); end n size(A, 1); isPos false; res struct(); res.Asym Asym; res.method method; switch method case chol [~, p] chol(Asym); isPos (p 0); res.info p; case eig d eig(Asym); if isempty(tol) tol n * eps(max(abs(d))); end isPos all(d tol); res.eigvals d; res.tol tol; case minor % 顺序主子式判据仅适合教学演示 if isempty(tol) tol n * eps; end isPos true; for k 1:n if det(Asym(1:k, 1:k)) tol isPos false; res.failStep k; break; end end otherwise error(未知方法: %s可选 chol/eig/minor, method); end res.flag isPos; end函数不长但每一块都有明确意图。比如skewPart那一项我曾经在处理一批传感器数据生成的协方差矩阵时被一个不对称的输入坑惨了协方差矩阵理论上是绝对对称的但采集数据在某些通道上有毫秒级的对齐误差导致算出来的矩阵出现 1e-8 量级的非对称扰动直接让chol在某个维度上崩掉。加了对称化之后再配合这个警告问题一下就暴露了。3.2 实战一判断样本协方差矩阵是否正定做统计的人最常遇到的一个问题就是协方差矩阵是否正定。如果样本量大于变量数而且变量之间没有严格线性关系样本协方差矩阵理论上是正定的但实际数据里常见缺失值插补、多重共线性、数值噪声都可能让情况偏离理论。我们构造一个三维数据人为设置三个不同的标准差和旋转相关结构然后计算它的样本协方差矩阵并判断rng(42); n 1000; X randn(n, 3) * [3 0 0; 0.5 1.5 0; 0.3 0.4 0.8]; % 样本协方差矩阵分母用 n-1 得到无偏估计 C X * X / (n - 1); % 三种方法依次判断 [isPos1, res1] isPositiveDefinite(C, chol); [isPos2, res2] isPositiveDefinite(C, eig); [isPos3, res3] isPositiveDefinite(C, minor); fprintf(chol 判断结果: %d\n, isPos1); fprintf(eig 判断结果: %d最小特征值: %.6f\n, isPos2, min(res2.eigvals)); fprintf(minor判断结果: %d\n, isPos3);不出意外三次判断都会返回 1。因为矩阵的行列式大于 0特征值全为正顺序主子式也都为正。这个例子主要是演示函数的标准用法和输出结构。实际跑下来你会发现minor方法在小矩阵上也能正常工作但如果你把示例里的 X 维度加到 100 维同时让样本数缩到和维度接近minor就会开始出现某些顺序主子式因为数值误差变得非常接近 0 却又不完全等 0 的情况。这就是我前面说的数值问题顺手可以体验一下。3.3 实战二接近奇异的病态矩阵与正则化处理更真实的问题来自一个接近奇异的矩阵。比如你在做岭回归或者高斯过程时核矩阵本身应该正定但由于样本里有几乎重复的点核矩阵的最小特征值会小到 1e-12 以下chol会直接失败而eig告诉你“它在数值意义上已经半正定了”。我构造一个极端的例子两个高度相关的变量生成协方差矩阵相关性系数 0.999999rho 0.999999; S2 [1 rho; rho 1]; % 理论上是正定但矩阵病态严重 condVal cond(S2); [isPosC, resC] isPositiveDefinite(S2, chol); [isPosE, resE] isPositiveDefinite(S2, eig); fprintf(条件数: %.3e\n, condVal); fprintf(chol 判断: %d\n, isPosC); fprintf(eig 判断: %d最小特征值: %.4e\n, isPosE, min(resE.eigvals));这个矩阵的两个特征值分别是 1 rho ≈ 1.999999 和 1 - rho ≈ 1.000000000001e-6最小的那个是正的所以在数学上是正定。但由于chol要对这个接近 0 的主元开平方根实际计算中可能因为舍入误差把它当成 0 处理导致 p 0。这就是理论正定和数值正定之间的灰色地带。碰到这种问题工程上通常用正则化处理给对角线上加一个很小的正值让矩阵严格离开奇异区域alpha 1e-8; S2_reg S2 alpha * eye(2); [isPosC2, resC2] isPositiveDefinite(S2_reg, chol); fprintf(正则化后 chol 判断: %d\n, isPosC2);加了对角扰动之后最小特征值变成约 1e-6 1e-8明显为正chol也就成功了。要注意这个操作实际判断的已经不是原来那个矩阵而是它附近的良态矩阵。这跟岭回归里的“岭参数”是同样的思想牺牲一点点无偏性换取数值稳定性和唯一解。我会建议任何做矩阵正定判断的人都养成一个习惯——当chol失败但理论上矩阵应当正定时不要急着下结论说矩阵不正定先检查条件数、最小特征值再考虑加对角扰动这个操作。4. 常见问题与避坑技巧实录每次写这类数值判断代码总能在细枝末节上栽跟头。下面这些问题都是我在实际项目里真实遇到过的不是从教科书上抄来的。4.1 矩阵不对称到底应该怎么办这是最普遍的一个坑。很多教材默认“正定矩阵本身是实对称矩阵”但工程里拿来的矩阵哪会那么干净。我前面已经讲过二次型只取决于矩阵的对称部分所以对任意方阵 A判断正定的严谨做法是先把 A 替换成 Asym (A A) / 2再用chol或eig去判断。为什么要这样做因为如果直接用chol(A)MATLAB 默认只读上三角忽略下三角那么下三角里那些可能写错的数据就直接被无视了你判断的根本不是你手头那个矩阵。举个例子假设你的矩阵是 [4 2; 1 3]看着差不多对称但严格意义上它不是对称矩阵。用chol算的时候MATLAB 会把上三角 [4,2; 0,3] 镜像成 [4 2; 2 3] 来处理而不是用你下三角那个 1。如果下三角的 1 是真实数据那这个结果就完全错了。所以我在封装函数时坚持先检查非对称偏差再强制对称化这条规则让我的判断代码在整个部门里几乎从不背锅。4.2 浮点误差特征值 1e-16 到底算正数还是零浮点世界里根本没有精确的 0只有“在误差范围内可以当作 0”的数。一个理论半正定的矩阵eig算出来最小特征值可能是 -1e-16直接all(d 0)就返回 false另一个接近奇异的正定矩阵最小特征值可能是 1e-17倒是all(d 0)返回 true 了但这个“正定”毫无意义因为任何一点输入扰动都可能让它翻负。正确做法是设置一个容差常用的经验公式是tol n * eps(max(abs(d)))。这个值约为 n 乘以最大特征量级的浮点间隔相当于“在这个矩阵本身的数值精度范围内认为小于它的量都是噪声”。我见过不少论文复现代码里直接写if min(eig(A)) 0结果就是复现时出现天壤之别。如果你也想严谨一点我建议至少把我上面的经验公式用上并且在论文或者项目文档里写明这个容差取值别人也能复现你的结论。4.3 半正定矩阵的判断与处理半正定矩阵在数据领域太常见了样本数小于变量数的协方差矩阵必然只有部分方向有方差这种矩阵一定是半正定而不是正定。使用chol判断它结果是 p 0函数只会告诉你“不正定”但不会告诉你“它是半正定”。所以如果你想区分半正定和非正定必须靠eig看特征值只要所有特征值在容差范围内都大于等于 0它就是半正定。处理半正定矩阵也有两种典型需求。如果你只需要“判断”那就用eig方法加容差判断。如果你下一步要做 Cholesky、求逆、解方程半正定矩阵会让这些操作全部崩掉。实用做法是给对角加一个极小正则项 A alpha*eye(n)alpha 根据矩阵量级选比如 1e-8 到 1e-12。虽然严格来说你算的是一个新矩阵但它与原矩阵在数值上非常接近在统计建模里这等价于岭回归、正则化最小二乘是可以接受的折中。4.4 大矩阵与稀疏矩阵的判断效率当矩阵规模变大——比如有限元模型的刚度矩阵能到十万阶——再调用完整的eig(A)基本属于灾难行为因为特征值分解会把一个稀疏矩阵填充成稠密矩阵内存和时间都扛不住。这种时候chol依然很争气因为 MATLAB 的chol原生支持稀疏矩阵分解过程中能尽量维持稀疏结构。如果你的问题是“这个大稀疏矩阵正定吗”最直接的做法就是[~, p] chol(A)p 为 0 就是正定。如果它不满足正定可你又想知道最小特征值大概是多少可以考虑用eigs(A, 2, smallestreal)求最小的几个特征值避免全量特征分解。但要小心eigs对于特征值非常接近 0 的矩阵收敛很慢有时求出来的结果也不稳定。所以我的习惯是稀疏大矩阵先上chol看看是否通过通过就万事大吉不通过再退到eigs去诊断而不是一开始就死磕特征值。4.5 常见问题速查表现象可能原因推荐处理方式chol(A)报错提示矩阵必须正定矩阵确实不正定或数值上接近半正定用[~,p]chol(A)避免报错再用eig检查最小特征值chol判断结果是 false但理论上应该正定矩阵被污染存在 1e-8 量级非对称扰动或共线变量对称化后重测检查cond(A)必要时加对角正则项eig(A)返回复数特征值A 不是对称矩阵先执行Asym(AA.)/2再对Asym判断特征值里有 -1e-16理论 0 的浮点残差设置tol n*eps(max(abs(d)))大于 -tol 就按 0 处理顺序主子式都大于 0 但chol失败det数值误差或矩阵练习病态以chol和eig结果为准不要迷信det大稀疏矩阵判断正定太慢用了eig而没用chol直接用[~,p] chol(A)判断诊断时用eigs求最小特征值结尾这套判断方法和例程我在多个项目里反复验证过最大的体会是多数“正定判断失败”其实是数据质量问题不是矩阵本身数学性质的问题。拿到一个理论应当正定的矩阵却被chol拒绝时先别急着怀疑算法和 MATLAB先检查矩阵是否对称、是否含噪声、是否病态。我的习惯动作是先跑一遍issymmetric(A)看一眼min(eig(A))和cond(A)这三个检查能解决掉九成以上的异常。另一个非常实用的技巧是把isPositiveDefinite这种函数封装好之后丢进自己的工具库以后每次遇到协方差矩阵、Hessian 矩阵、刚度矩阵的判断需求都先统一过一遍函数输出结构统一调试时候真的能省下大把时间。矩阵正定判断本身不复杂但它在数值算法里的地位就像地基里的钢筋——平时你注意不到它一旦它出了问题整栋楼都得跟着遭殃。