如何用矩阵运算实现多元高斯对数似然第三项的R高效计算

问题背景 #

在CrossValidated上看到多元高斯分布的对数似然公式：
$$
\begin{eqnarray}
\ell(\mu, \Sigma) &=& C - \frac{m}{2}\log|\Sigma|-\frac{1}{2} \sum_{i=1}^m \text{tr}\left[(\mathbf{x}^{(i)}-\mu)T \Sigma^{-1} (\mathbf{x}^{(i)}-\mu)\right]
\end{eqnarray}
$$
其中第三项可以去掉迹运算（因为标量的迹等于自身），已通过for循环实现计算，但希望用更简洁的矩阵运算优化，并理解线性代数等价性。

高效矩阵运算实现 #

以下是三种替代for循环的高效方法，结果与原循环完全一致：

library(MASS)

# 构造协方差矩阵
S = matrix(c(1.0, 0.2, 0.1, 0.35, 0.0,
             0.2, 1.0, 0.0, 0.4, 0.0,
             0.1, 0.0, 1.0, 0.0, 0.4,
             0.35, 0.4, 0.0, 1.0, 0.6,
             0.0, 0.0, 0.4, 0.6, 1.0), ncol = 5)
colnames(S) = c("Y1", "X1", "X2", "Z1" ,"Z2")
rownames(S) = colnames(S)

# 均值向量
mus = c(1, 2, 3, 4, 5); names(mus) = colnames(S)

# 生成观测数据
obs = mvrnorm(n = 200, mu = mus, Sigma = S)

# 原循环实现
llmat.term3 <- matrix(NA, 200,1)
for(n in 1:200) {
  llmat.term3[n] <- t(obs[n,]-mus) %*% solve(S) %*% (obs[n,]-mus)
}
sum_loop <- sum(llmat.term3)
sum_loop
#[1] 982.7356

# 方法1：点积等价矩阵运算
U <- obs - rep(mus, each = nrow(obs))  # 中心化观测矩阵（每行一个样本）
sum_matrix <- sum(U * solve(S, t(U)))
sum_matrix
#[1] 982.7356

# 方法2：矩阵对角线求和
sum_diag <- sum(diag(U %*% solve(S) %*% t(U)))
sum_diag
#[1] 982.7356

# 方法3：利用马氏距离函数
sum_mahalanobis <- sum(mahalanobis(obs, center = mus, cov = S))
sum_mahalanobis
#[1] 982.7356

线性代数等价性解释 #

1. 迹运算的冗余性 #

公式中的迹运算可以直接去掉：因为$(\mathbf{x}^{(i)}-\mu)T \Sigma^{-1} (\mathbf{x}^{{(i)}-\mu)$是一个**1×1的标量矩阵**，而标量的迹等于其自身，因此求和项可简化为$\sum_{i=1}}m (\mathbf{x}^{(i)}-\mu)T \Sigma^{-1} (\mathbf{x}^{(i)}-\mu)$。

2. 点积与矩阵乘法的等价性 #

对于单个中心化样本向量$\mathbf{u}^{(i)} = \mathbf{x}^{(i)} - \mu$，$\mathbf{u}^{(i)T} \Sigma^{-1} \mathbf{u}^{{(i)}$本质是$\mathbf{u}}{(i)}$与$\Sigma^{-1} \mathbf{u}^{(i)}$的点积：
$$
\mathbf{u}^{(i)T} \Sigma^{-1} \mathbf{u}^{(i)} = \sum_{j=1}^p u_j^{(i)} \cdot (\Sigma^{-1} \mathbf{u}^{(i)})j = \sum{j=1}^p \sum_{k=1}^p u_j^{(i)} (\Sigma^{-1})_{jk} u_k^{(i)}
$$
其中$p$为变量维度。因此所有样本的项之和，等于将中心化矩阵$U$与$\Sigma^{-1} U^T$转置后的矩阵对应元素相乘，再对所有元素求和，这就是sum(U * solve(S, t(U)))的底层逻辑。

3. 高效计算的原理 #

solve(S, t(U))直接求解线性方程组$S \cdot X = t(U)$，得到$X = \Sigma^{-1} t(U)$，比先计算$\Sigma^{-1}$再与$t(U)$相乘更高效，数值稳定性也更强，尤其适合高维度协方差矩阵场景。

内容的提问来源于stack exchange，提问作者IRTFM