第一次认真啃多元正态分布是在做一个多通道传感器异常检测项目的时候。当时我手里有一批十几维的特征数据想用马氏距离判定哪个点偏得离谱结果直接把协方差矩阵求逆求炸了——矩阵接近奇异逆矩阵里的元素大得离谱算出来的距离全是噪声正常的样本被标成异常真正的异常反而藏起来了。那一刻我才意识到一元正态里那个简简单单的 σ²搬到多维之后变成了一件多么麻烦的事它不再是一个数而是一整个协方差矩阵里面既藏着维度之间的全部相关性也藏着数值计算中的大部分坑。这篇文章想聊的就是多元正态分布Multivariate normal distribution到底是怎么回事它在工程里怎么用以及它在代码里会以什么姿势坑你。内容会覆盖协方差矩阵的几何含义、密度函数的来历、条件分布和边缘分布这两个最值钱的性质、仿射变换下的封闭性、参数估计的闭式解以及白化、采样、数值稳定性这些真正落到代码上的环节。不管你是刚学概率论的学生还是已经在做卡尔曼滤波、高斯过程、混合高斯聚类、风险建模的工程师下面这些应该都能对得上号。1. 从一元正态到多元正态协方差矩阵才是真正的主角1.1 一元正态那两个参数到了高维为什么就不够用一元正态 N(μ, σ²) 只有两个参数均值决定位置方差决定胖瘦。你画一条钟形曲线它是对称的、单峰的形状完全由这两个数固定下来。很多人第一次接触多元正态会很自然地想那就每个维度各自来一个一元正态不就行了这个想法对应的模型是各维度独立也就是协方差矩阵是个对角矩阵非对角元素全是零。它确实是多元正态的一个特例但它丢掉的东西恰恰是实际数据里最值钱的部分——维度之间的相关性。举个具体例子。假设你在做一个传感器阵列两个通道测的是同一个物理量只是位置略有不同。这两个通道的读数几乎总是同涨同跌你单独看每个通道的分布都是标准的一元正态可你一旦把它们放到同一张散点图上会发现所有点几乎落在一条斜线上而不是铺满整个矩形区域。这个斜就是相关性一元分布看不出来对角协方差矩阵也表达不出来。要描述它你必须引入协方差这个量Cov(x₁, x₂) E[(x₁-μ₁)(x₂-μ₂)]它刻画两个维度一起偏离均值的倾向。当协方差为正两个维度同向偏离为负反向偏离为零才叫不相关在联合正态假设下就等于独立。所以多元正态的完整参数是均值向量 μk 维加上协方差矩阵 Σk×k 对称正定。为什么必须是正定而不是半正定这是一个经常被忽略的细节。正定意味着任意非零方向的方差都严格大于零概率密度才在每一点都有定义、可归一化。如果 Σ 只是半正定比如有零特征值那么数据实际上被压在一个低维子空间里这种情况叫退化分布密度函数在通常意义下不存在要用广义函数或者降维处理。后面讲数值计算的时候我会专门说实践中遇到的协方差矩阵求逆炸了多半就是这个问题。1.2 协方差矩阵的几何解读椭球、主轴和方向把 Σ 做特征值分解Σ QΛQ⁻¹因为 Σ 对称Q 是正交矩阵列向量是一组标准正交基Λ 是对角矩阵对角元是特征值 λ₁ ≥ λ₂ ≥ ... ≥ λ 0。这个分解是理解多元正态几何形状的钥匙。密度函数里的等值面 (x-μ)ᵀΣ⁻¹(x-μ) c² 是一个椭球。这个椭球的中心在 μ主轴方向由 Q 的列向量给出半轴长度是 c√λᵢ。也就是说特征值大的方向对应椭球胖的方向数据在那个方向散布得开特征值小的方向对应椭球瘦的方向数据很快收敛到均值附近。特征值越小说明那个方向的方差越小数据越集中越接近退化。举个二维的例子帮你建立直觉。如果 Σ [[4, 0], [0, 1]]Q 就是单位矩阵Λ 对角是 4 和 1椭球的长轴沿 x₁ 方向长度 2√4短轴沿 x₂ 方向长度 1数据呈横向拉长的椭圆。如果 Σ [[2, 1.8], [1.8, 2]]特征值大约是 3.8 和 0.2特征向量大约在 (1,1)/√2 和 (1,-1)/√2 方向。这个形状是一个沿 45 度方向拉得很长的细椭圆对应两个维度高度正相关。特征值 0.2 很小说明垂直于这个方向的数据几乎不动你在这个方向上稍微偏一点都会得到巨大的马氏距离。理解这一点非常重要因为它直接解释了为什么马氏距离在异常检测里这么灵敏它把数据空间按方差归一化然后在方差小的方向上给予极端的惩罚。一个在原始欧氏距离下看起来很小的偏移如果恰好发生在特征值很小的方向上马氏距离会大得惊人。1.3 独立、不相关、正定三个经常被混为一谈的概念这里要专门拎出来讲三个概念因为它们经常被混着用。独立是说两个随机变量的联合分布等于各自分布之积它是个很强的条件。不相关是说协方差为零它只描述线性关系。联合正态假设下不相关等价于独立这句话是成立的但前提是联合正态。如果不是联合正态不相关推不出独立这在很多教材习题里都是经典陷阱。还有一个容易搞混的是协方差为零和协方差矩阵是对角矩阵。前者是单个元素的性质后者是整个矩阵的性质。对角协方差矩阵意味着所有维度两两不相关那么在正态假设下就是完全独立联合密度可以写成各维度一元密度的乘积。这也是为什么很多算法比如朴素贝叶斯会假设特征独立——一旦独立协方差矩阵变成对角求逆变成逐元素取倒数计算量从 O(k³) 降到 O(k)省事得多。但代价是丢掉了所有维度间的相关性信息模型能力下降。至于正定我在前面提过。判断一个协方差矩阵是不是正定最实用的方法不是算特征值而是试图做 Cholesky 分解Σ LLᵀ。Cholesky 能成功就是正定报错或者对角元出现非正数就不是。这个方法比直接算特征值快而且顺手给你分解结果后面求逆、求行列式、采样都能接着用。我在代码里检查协方差矩阵时基本都这么干。2. 密度函数那个归一化常数到底是干嘛的2.1 二次型式与马氏距离多元正态的密度函数长这样f(x) (2π)^(-k/2) |Σ|^(-1/2) exp( -1/2 (x-μ)ᵀ Σ⁻¹ (x-μ) )指数里那个 (x-μ)ᵀ Σ⁻¹ (x-μ) 就是马氏距离的平方记作 d²(x, μ)。它和欧氏距离的差别就在中间那个 Σ⁻¹。欧氏距离是 (x-μ)ᵀ(x-μ)每个方向的偏离权重都一样马氏距离用 Σ⁻¹ 给每个方向加权方差小的方向权重1/λᵢ大方差大的方向权重小。你可以把它理解为先白化再做欧氏距离。我在实际项目里踩过一个坑用马氏距离做异常检测时如果协方差矩阵是从带异常点的数据里估出来的异常点本身会把某些方向的方差拉大导致这些方向的马氏距离被稀释真正的异常反而不容易被检出。这时候通常要用稳健估计比如最小协方差行列式MCD或者基于分位数的重估。这个坑我后面会专门展开。密度函数里那个 exp 前面的常数是 (2π)^(-k/2) |Σ|^(-1/2)它的唯一作用是把整个函数在 ℝᵏ 上的积分归一化成 1。很多人写代码的时候直接忽略常数项只比较指数部分这在做分类、做似然比的时候是合理的因为常数对所有样本一样不影响比较结果但如果你要做概率密度回归、要算真实的概率值就不能省略。2.2 归一化常数的来龙去脉为什么常数项里会有 |Σ|^(-1/2)这个可以通过变量代换推导。做变换 z Σ^(-1/2)(x-μ)那么 x μ Σ^(1/2) zdx |Σ^(1/2)| dz |Σ|^(1/2) dz。在 z 空间里二次型变成 zᵀz就是标准正态的模平方。于是积分 ∫ exp(-1/2 (x-μ)ᵀΣ⁻¹(x-μ)) dx |Σ|^(1/2) ∫ exp(-1/2 zz) dz |Σ|^(1/2) (2π)^(k/2)。所以为了让密度积分为 1系数必须乘上 (2π)^(-k/2) |Σ|^(-1/2)。这里 |Σ|^(1/2) 就是雅可比行列式是变量代换带来的体积缩放因子。理解来源之后你在代码里看到 log|Σ| 就不会觉得突兀了——对数似然里那一项 -1/2 log|Σ| 正是来自这个归一化常数它是体积罚项惩罚那些把椭球撑得太大的协方差参数。举个具体的数值感受一下。如果 k 10Σ 是单位矩阵|Σ| 1常数是 (2π)^(-5) ≈ 0.00117。密度最大处在原点值约 0.00117。你得意识到高维空间里密度函数的绝对数值都很小所以直接比较密度值的大小没什么意义通常都要转成对数。这也是为什么我强烈建议所有涉及高维正态的计算都走对数域后面会详细讲。2.3 什么时候密度函数写不出来我在 1.1 里提过退化的情形这里展开说。如果 Σ 奇异行列式为零|Σ|^(-1/2) 就没有意义除以零密度函数在通常意义下不存在。这不是数学上的吹毛求疵实践中非常常见。两个典型场景。第一个是数据维度大于样本数比如你有 100 个样本、500 个特征样本协方差矩阵的秩最多是 99必然是奇异的。第二个是特征之间存在精确的线性关系比如你把某三个变量的和作为第四个变量录进去协方差矩阵的相关列线性相关也奇异。处理方式有几种。最简单的是降维用 PCA 把数据投到主成分空间只保留非零特征值对应的方向。稍微复杂但更优雅的是加一点对角扰动Σ εIε 取个 1e-6 到 1e-3 之间保证正定这叫 ridge 正则或 Tikhonov 正则本质上就是往协方差矩阵里加一点各向同性的方差。如果数据有团结构也可以用因子模型或者图模型来约束协方差结构减少参数个数。我个人的经验是先诊断再处理。诊断就是看样本协方差矩阵的秩和条件数最大特征值除以最小特征值条件数超过 10^6 基本就可以认为数值上不可靠了。处理上如果样本量够n 远大于 k优先考虑用无偏估计加轻微正则如果样本量不够果断降维硬扛只会得到一堆数值垃圾。3. 条件分布与边缘分布这块才是每天真正要用到的3.1 分块记号与两个分量的分解多元正态最值钱的性质就是边缘分布还是正态条件分布也是正态。这两个性质让它在贝叶斯推断、滤波、回归里无处不在。先把记号立清楚。把随机向量分成两块x [x_a; x_b]对应的均值 μ [μ_a; μ_b]协方差矩阵分块为Σ [[Σ_aa, Σ_ab], [Σ_ba, Σ_bb]]其中 Σ_aa 是 x_a 的自协方差Σ_bb 是 x_b 的自协方差Σ_ab 是交叉协方差Σ_ba Σ_ab。分块的大小按你的需求切比如 x_a 是待推断的量x_b 是观测到的量。边缘分布有多简单x_a 的边缘分布就是 N(μ_a, Σ_aa)你只要把对应的行和列切出来就行其他部分直接扔掉。没有任何计算没有积分要算。这就是为什么多元正态在工程里这么友好——边缘化是一次切片概念成本几乎为零。条件分布稍微复杂一点但也是闭式的x_a | x_b ~ N( μ_a Σ_ab Σ_bb⁻¹ (x_b - μ_b), Σ_aa - Σ_ab Σ_bb⁻¹ Σ_ba )条件均值是 μ_a Σ_ab Σ_bb⁻¹ (x_b - μ_b)条件协方差是 Σ_aa - Σ_ab Σ_bb⁻¹ Σ_ba。后者叫舒尔补Schur complement它在数值线性代数里是个反复出现的角色。3.2 条件均值就是线性回归卡尔曼滤波的脊梁盯着条件均值那个式子看两秒钟μ_a Σ_ab Σ_bb⁻¹ (x_b - μ_b)。把 (x_b - μ_b) 记作 Δ这就是一个关于 Δ 的线性函数斜率矩阵是 Σ_ab Σ_bb⁻¹截距是 μ_a。如果你把 x_a 当作要预测的目标x_b 当作已知的特征那么这个条件均值就是最优线性预测。这其实就是线性回归回归系数等于 Σ_ab Σ_bb⁻¹和最小二乘的闭式解一模一样。换句话说在联合正态的假设下给定 x_b 预测 x_a这个问题的最优解是线性的而且系数有闭式表达不需要迭代优化。卡尔曼滤波整套理论就建立在这个上面。状态转移和观测方程都是线性的噪声是高斯的那么状态和观测的联合分布就是多元正态滤波的每一步做的就是计算给定观测之后状态的均值和协方差也就是上面那个条件分布的公式。只不过在卡尔曼滤波的记号里Σ_ab Σ_bb⁻¹ 这个斜率矩阵被叫做卡尔曼增益 KΣ_aa - Σ_ab Σ_bb⁻¹ Σ_ba 被叫做后验协方差。很多人学卡尔曼滤波时被那堆公式绕晕其实退一步看它就是在反复地做多元正态的条件分布计算。我当初就是靠这个视角才把卡尔曼滤波真正搞明白的。3.3 舒尔补为什么不能直接拿工程公式硬套公式好看落到实现上有讲究。条件均值里需要 Σ_bb⁻¹ (x_b - μ_b)注意这里是解一个线性方程组不是显式求逆。对是解而不是求逆——Σ_bb⁻¹ (x_b - μ_b) 应该用 solve(Σ_bb, x_b - μ_b) 来算而不是 inv(Σ_bb) (x_b - μ_b)。两者的数学结果相同数值上后者误差更大、速度更慢。这是我在代码评审里反复纠正新人的一点。import numpy as np # 推荐解线性方程组 delta x_b - mu_b adj np.linalg.solve(Sigma_bb, delta) cond_mean mu_a Sigma_ab adj # 不推荐显式求逆 # cond_mean mu_a Sigma_ab np.linalg.inv(Sigma_bb) delta条件协方差里的 Σ_aa - Σ_ab Σ_bb⁻¹ Σ_ba 也没必要真的去求 Σ_bb⁻¹。可以用 Cholesky 分解 Σ_bb LLᵀ然后算 Σ_ab L⁻ᵀ L⁻¹ Σ_ba。具体写法是先解 L Y Σ_ba前代再解 L Z Y回代最后拼出 Σ_ab Σ_bb⁻¹ Σ_ba。这个技巧在统计和信号处理里叫用求解法求舒尔补能显著提升数值稳定性。还有一个坑是条件协方差的正定性。理论上它是正定的因为它是原协方差在子空间上的限制但数值上如果 Σ_bb 病态算出来的 Σ_aa - Σ_ab Σ_bb⁻¹ Σ_ba 可能出现小的负特征值导致后续的 Cholesky 或者采样直接挂掉。处理办法是做对称化(A Aᵀ)/2加轻微的对角扰动或者干脆用平方根滤波square-root filter的形式直接在平方根因子上做更新避免形成完整的协方差矩阵。做惯导、做组合导航的朋友对这套应该不陌生。4. 线性变换、白化与采样工程实现里最常写的那几行代码4.1 仿射变换的封闭性多元正态对线性变换是封闭的这是个非常强的性质值得单独说。如果 x ~ N(μ, Σ)A 是一个 m×k 的矩阵b 是 m 维向量那么 y Ax b ~ N(Aμ b, AΣAᵀ)。推导不难用特征函数或者矩母函数都可以。这个性质的实用价值在于它让你可以用线性工具来构造、变换、简化正态分布不用反复算积分。降维投影是它白化也是它卡尔曼滤波里的状态转移还是它。只要变换是线性的严格说是仿射的带偏移结果一定是正态协方差按 AΣAᵀ 走。一个常被问到的推论两个联合正态的随机变量的和还是正态X Y ~ N(μ_x μ_y, Σ_x Σ_y Cov(x,y) Cov(y,x))。注意这里必须带上互协方差项如果 x 和 y 独立互协方差为零退化成 Σ_x Σ_y。我见过不少人做随机变量求和时直接把协方差相加忘了互协方差项结果方差算小了一半。做金融风险里组合方差、做测量误差传播的时候这个错误代价很大。4.2 白化把相关变不相关的标准操作白化的目标是把一个协方差为 Σ 的多元正态变成协方差为单位矩阵的标准正态。做法是找一个矩阵 W使得 WΣWᵀ I。这样的 W 不唯一常见的有三种。第一种是基于 Cholesky 分解的。令 Σ LLᵀ取 W L⁻¹那么 WΣWᵀ L⁻¹LLᵀL⁻ I。这种白化叫下三角白化好处是 W 是下三角的求解快是唯一确定的一种。第二种是基于特征分解的W Λ^(-1/2) Qᵀ。这种白化得到的坐标轴就是协方差椭球的主轴物理意义更清楚但计算成本高一些要算特征分解。白化之后的样本在原始空间里是沿主轴分开放置的。第三种是带子空间的白化W Λ_r^(-1/2) Q_rᵀ只保留 r 个非零特征值对应的方向把数据从 k 维压到 r 维。这是在数据奇异时最自然的处理方式PCA 白化通常指的就是这一种。白化的实际用途很广。数据预处理里它能让不同特征的量纲和相关性统一让后续的算法距离度量、梯度下降更稳定。在生成模型里采样常常是先采样标准正态再白化逆变换。在信号处理里白化是预白化滤波的基础目的是把有色噪声变成白噪声方便后续匹配滤波。4.3 采样从标准正态到目标正态采样这个操作看起来简单——x μ Lzz 是 k 维标准正态L 是 Σ 的 Cholesky 因子——但细节不少。首先为什么用 Cholesky 而不是特征分解Cholesky 是 O(k³/3)特征分解大概要多几倍的成本虽然都是 O(k³)而且 Cholesky 因子是下三角在需要重复采样的时候特别方便因为 L 算一次就能一直复用。只有当 Σ 非常接近奇异Cholesky 会失败的时候才退回到特征分解用 Σ^(1/2) QΛ^(1/2) 来做。特征分解还能天然地做截断把负的或极小的特征值置零保证半正定。其次采样的数值稳定性有个容易忽略的点如果你要采一大批样本不要每次重新做 Cholesky。Σ 不变的话L 算一次就够了。我见过有人在一个百万次循环里每次都算 np.linalg.cholesky(cov)运行时间直接爆炸。正确的做法是把 L 提到循环外面。import numpy as np k 50 mu np.zeros(k) Sigma np.eye(k) 0.1 * np.random.randn(k, k) Sigma Sigma Sigma.T # 构造一个正定矩阵 # 分解一次重复使用 L np.linalg.cholesky(Sigma) n_samples 100000 z np.random.randn(n_samples, k) samples mu z L.T # 每行是一个样本再就是随机数种子的管理。做可复现实验的时候采样器最好有独立的随机数生成器不要和全局的 np.random 混用否则实验的可重复性会变得很糟。我自己项目里习惯每个模块用独立的 Generator 对象用固定种子初始化这样跑两遍结果完全一致。还有一个纯工程问题k 很大的时候比如 k 1000Cholesky 本身也可能因为浮点误差而失败。这时候通常用带 jitter 的版本对 Σ 的对角线加一点正数或者用 LDL 分解代替 CholeskyLDL 不需要开方对半正定矩阵也适用。5. 参数估计极大似然那套闭式解背后的事5.1 均值和协方差的极大似然估计给定 n 个独立同分布的样本 x₁, ..., xₙ每个都是 k 维假设来自 N(μ, Σ)极大似然估计是有闭式解的μ̂ (1/n) Σ xᵢΣ̂ (1/n) Σ (xᵢ - μ̂)(x - μ̂)ᵀ第一个就是样本均值第二个是样本协方差注意分母是 n 而不是 n-1。这个推导过程用对数似然后求导就能得到不算难。有意思的是μ̂ 和 Σ̂ 是分开估计的不像有些分布那样耦合这也是多元正态好用的原因之一。但我要强调一点协方差的极大似然估计是有偏的。在无偏性检验里E[Σ̂] (n-1)/n Σ所以当 n 比较小的时候偏差不可忽略。举个例子n 10、k 3 的时候偏差因子是 0.9主对角线的方差被系统性低估了 10%。无偏估计就是乘上 n/(n-1)Σ̂_unbiased (1/(n-1)) Σ (xᵢ - μ̂)(xᵢ - μ̂)这就是为什么 numpy 的 np.cov 默认用 ddof1用 n-1 做分母。5.2 到底用 n 还是 n-1这是一个取舍问题用 n 还是 n-1这个问题没有绝对答案取决于你接下来要干什么。从统计推断的角度如果你要构造无偏估计并且关心估计量的期望用 n-1。如果你要做极大似然并且关心的是让似然最大化用 n因为极大似然是给出一组参数让观测到的数据在给定模型下的概率最大不追求无偏。这两种目标本身就不一样。从工程实践的角度当 n 远远大于 k 的时候比如 n 100k两者差异可以忽略用哪个都行。真正有影响的是小样本、高维的场景这时用 n-1 会好一点但更本质的问题其实是样本量根本不够加 1 个偏差修正救不了你。还有一个更隐蔽的坑如果你先估计 μ̂再用同一个数据估计 Σ̂两者不是独立的虽然最大似然意义下它们是分开的最优解但有限样本下样本均值和样本协方差是负相关的这个影响很小在极端高维下会被放大。这种时候可以考虑用留一法或者交叉验证来估计或者对 μ 和 Σ 都用更稳健的估计方法。5.3 高维陷阱样本协方差矩阵什么时候不可信当特征维度 k 接近或超过样本量 n 的时候样本协方差矩阵一定不可信。具体来说如果 k ≥ n样本协方差矩阵必然奇异秩最多 n-1求逆直接失败。如果 k 接近 n比如 k 0.5n虽然理论上可逆但条件数非常大最大特征值受少数几个方向主导最小特征值接近零求逆得到的矩阵数值误差巨大。这时候用马氏距离做异常检测结果基本是乱七八糟的。即使 k 远小于 n如果数据本身存在强共线性比如特征之间有近似线性关系条件数也会飙高。高维场景下标准解法有几种。收缩估计shrinkage把样本协方差往目标矩阵通常是单位矩阵或者对角矩阵方向拉Σ̂_shrink (1-α) Σ̂ α Tα 用 Ledoit-Wolf 方法选。因子模型假设协方差由少数几个公共因子驱动参数个数大幅减少。稀疏逆协方差估计graphical lasso直接估计 Σ⁻¹ 并施加稀疏约束避免显式求逆。低秩加对角分解把 Σ 写成 UUᵀ D 的形式U 是低秩的载荷矩阵D 是对角矩阵这是因子分析的标准形式。我自己的经验是先算条件数看数据是不是高维小样本。如果是别硬上选一个合适的结构化模型。收缩估计是最省事的实现简单超参数有理论指导因子模型解释性强但参数估计麻烦稀疏逆协方差的调参成本高。具体选哪个看你的下游任务和你有多少调参预算。6. 实际用起来会踩的坑6.1 马氏距离做异常检测的三个陷阱马氏距离是多元正态自然的离群程度度量d²(x) (x-μ)ᵀΣ⁻¹(x-μ)理论上服从自由度为 k 的卡方分布。所以一个很自然的做法是算 d²和卡方分布的分位数比如 chi2.ppf(0.99, dfk)比超了就判异常。第一个陷阱是 μ 和 Σ 的估计被异常点污染。前面提过样本均值和协方差对离群点都不稳健。一个离均值很远的点会把均值拉过去把自己的距离降低同时会把某些方向的方差拉大稀释其他异常。推荐用稳健估计最小协方差行列式MCD是最常用的方法它在 sklearn 里叫 MinCovDet实现上就是找一个子集让这个子集的协方差矩阵的行列式最小再基于它估计位置和散布。经验上当污染比例小于 20% 时 MCD 很可靠。第二个陷阱是掩盖和淹没效应。掩盖是多个异常点互相掩护让彼此的马氏距离看起来正常淹没是一个异常点把协方差撑大导致其他真异常被误判为正常。这两个效应在经典的多元离群点检测文献里讨论很多处理方法是重加权迭代先用一个宽松阈值筛出疑似正常点用它们重新估参数再算距离。这个迭代重复两三轮基本就收敛了。第三个陷阱是阈值的选择。如果数据本身不严格正态卡方分位数给出的阈值就会偏。轻尾的数据用卡方阈值会误报太多重尾数据会漏报。实践中我会用一个基于数据本身的稳健分位数比如把 d² 排序后取 97.5% 分位作为阈值再结合业务含义去调整。别迷信理论阈值实际数据很少完全满足多元正态假设。6.2 混合高斯与 EM 里的数值下溢混合高斯模型GMM是多元正态的组合版它假设数据由若干个多元正态按一定权重混合而成用 EM 算法估计参数。在实现 EM 的时候有一个绕不过去的数值问题似然值下溢。多变量密度函数在高维下取值很小前面举例 k10、单位协方差时峰值才 0.00117如果维度更高密度值可能小到 10^-30 以下直接相乘就下溢到 0 了。这时候必须在对数域里算log-sum-exp 技巧。对每一个样本计算它属于每个高斯分量的对数似然然后在 log 域做加权求和再用 log-sum-exp 转换成概率。import numpy as np def log_sum_exp(log_weights): M np.max(log_weights) return M np.log(np.sum(np.exp(log_weights - M)))这个技巧保证不会上溢也不会下溢是所有涉及混合模型和对数似然求和的地方必须用的。我审过的代码里超过一半的 EM 实现对这个问题处理得不好要么日志里出现 nan要么在小数据集上看起来正常数据集一放大就崩。还有一个常被忽略的问题是协方差的数值退化。EM 迭代过程中某个分量的协方差可能变得越来越小数据点高度集中最终接近奇异行列式为零log|Σ| 变成负无穷。解决办法是给协方差加一个最小正则项Σ εIε 取 1e-6 左右。或者在每次迭代后检查协方差的最小特征值如果低于阈值就往单位矩阵方向拉一点。这两个方法在 sklearn 的 GaussianMixture 里都有实现对应 reg_covar 参数。我做实验的时候默认把 reg_covar 设成 1e-6关键项目上调到 1e-4代价是稍微降低模型的拟合紧致度但换来数值稳定。6.3 高维退化带来的连锁反应最后聊一个比较隐蔽但很致命的问题高维多元正态的集中现象。随着维度 k 增加多元正态的样本会越来越集中在一个半径约 √k 的薄壳上。同时任意两个样本之间的距离会趋于相等——这是维度诅咒的一个具体表现。这意味着在高维下马氏距离的区分度在下降欧氏距离基本失效很多基于距离的算法都会退化。对多元正态来说这个现象的数学含义是当 k 很大时密度函数的大部分质量集中在椭球的表壳上而不是中心附近。这和低维直觉相反低维下最大密度在中心也会让人以为概率质量集中在中心是需要习惯的一件事。另一个连锁反应是协方差矩阵的估计噪声。即使 n 和 k 的比例还不错高维下样本协方差的特征值谱也会扭曲最大特征值被高估最小特征值被低估。这叫 Marchenko-Pastur 现象是随机矩阵理论里的经典结果。如果你直接拿高维样本协方差做 PCA会得到一堆虚假的主成分做聚类会得到虚假的簇结构做异常检测会把随机噪声当成异常。处理这个问题的工具箱也比较成熟了。特征值谱的校正可以用 Marchenko-Pastur 分布去拟合噪声部分然后做截断或者收缩。随机矩阵去噪Random Matrix Theory based denoising适用于信噪比中等偏低的场景。低秩修正把大的特征值对应的方向保留小的统一收缩到一个共同值保证协方差矩阵的条件数可控。我个人的判断标准是如果 k 50 或者 k/n 0.1就要警惕这些问题了如果 k 200基本必须上结构化估计或者降维。这个阈值不是绝对的取决于数据的信噪比和相关性结构但作为一个经验起点够用了。最后一点个人体会。多元正态这个东西初学的时候会觉得它只是一元正态的推广公式一堆、看着吓人真正用起来之后才发现它的价值恰恰在于那些让公式变复杂的性质——封闭性、条件分布、边缘分布。你把这些性质吃透会在很多地方突然看懂卡尔曼滤波是条件分布高斯过程是无限维正态线性回归在正态假设下是条件均值GMM 是正态的凸组合PCA 是正态椭球的坐标变换。它们不是五个不同的知识点而是同一个东西的五张面孔。如果要我给一条最实用的建议把 Cholesky 分解当成你的默认工具。求逆、求行列式、采样、白化、算马氏距离、做条件分布它都能用而且比显式求逆稳得多。我做过很多次对比同一份数据、同一个模型把 inv 换成 solve、把 det 换成对角元的对数求和数值误差能差好几个数量级。这些细节不会写在教科书里但它们在真实项目里决定了你的代码是跑得通还是跑得对。