本篇主要回答下面这个问题：

Faq

当维度 $d$ 很大时，如何高效计算 $\det(\Sigma)$ 和 $\Sigma^{-1}$ ？如何高效计算矩阵-向量乘法？

这个问题其实是两部分，而且第二句（矩阵-向量乘法）才是关键——它在暗示我们：我们根本不需要真的算出 $\Sigma^{-1}$ .

答案：用矩阵分解代替求逆。

对对称正定的 $\Sigma$ 做 Cholesky 分解 $\Sigma=LL^T$ ，于是

$\det(\Sigma)=(\prod_i L_{ii})^2$ ，计算成本极小；
$\Sigma^{-1}v$ 通过两次三角求解（前代 + 回代）得到，每次 $O(d^2)$ ，不用显式求逆。

总代价：一次分解 $O(d^3)$ ，之后每次查询 $O(d^2)$ 。

1. 如何高效计算 $\det(\Sigma)$ 和 $\Sigma^{-1}$ #

1.1 “直接求逆” 的坏处#

朴素做法有三宗罪：

做法	代价	问题
显式算 $\Sigma^{-1}$ （如 Gauss–Jordan）	$O(d^3)$	数值不稳定；多花 $O(d^2)$ 内存
直接展开算 $\det(\Sigma)$	比 $O(d^3)$ 还糟	大 $d$ 下上溢/下溢
然后矩阵-向量乘	每次 $O(d^2)$	但其实分解后能更省

1.1.1 为什么显式求逆是 $O(d^3)$ ？#

我们可以数”扫过矩阵的次数”：消到第 $k$ 个主元时，要把主元下方约 $d$ 行的该列清零，每次行变换要更新整行（约 $d$ 个数），故单个主元 $O(d^2)$ ；共 $d$ 个主元：

\underbrace{d}_{\text{主元个数}}\times\underbrace{O(d^2)}_{\text{单主元工作量}}=O(d^3)

更精确地，待消的块逐次缩小， $\sum_{k=1}^{d}k^2\approx\tfrac13 d^3$ ，加上回代共约 $\tfrac23 d^3$ ——这正是 Gauss 消元 / LU 是 $\tfrac23 d^3$ 的来源。

关于避免直接求逆，最核心的原因在于：在实际问题里，我们几乎从不需要 $\Sigma^{-1}$ 本身，我们需要的只是它能作用于一个向量：

\Sigma^{-1}v,\qquad \text{以及}\qquad \det(\Sigma)\ (\text{或}\ \log\det\Sigma).

也就是说，所有用到 $\Sigma^{-1}$ 的公式，都会以上面两种形式出现。我们以多元高斯对数密度为例：

\log p(x)=-\frac d2\log(2\pi)-\frac12\underbrace{\log\det\Sigma}_{\text{标量}}-\frac12\underbrace{(x-\mu)^T\Sigma^{-1}(x-\mu)}_{\text{一次"咬向量"}}

拆开看第二项：先算 $u=\Sigma^{-1}(x-\mu)$ ，再算 $(x-\mu)^T u$ （一个内积）。整个过程只用了”作用一次”的结果，从来没有需要逆矩阵本身。

于是就可以跳过求逆，改用分解（Cholesky）：

\Sigma = LL^T \;\Longrightarrow\; \Sigma^{-1}v \text{ 变成解 } Ly=v \text{ 再解 } L^Tx=y

下面我们将详细介绍 Cholesky 分解。

1.2 Cholesky 分解#

Cholesky Decomposition

设 $A\in\mathbb{R}^{d\times d}$ 对称正定（symmetric positive definite, SPD），即

A=A^T,\qquad \mathbf{x}^TA\mathbf{x}>0\quad \text{对所有 } \mathbf{0}\neq\mathbf{x}\in\mathbb{R}^d

则存在唯一的下三角矩阵 $L\in\mathbb{R}^{d\times d}$ ，其对角元全为正（ $L_{ii}>0$ ），使得

A=LL^T

$L$ 称为 $A$ 的 Cholesky 因子（Cholesky factor）。

更详细的介绍参见 Cholesky Decomposition.

协方差矩阵 $\Sigma$ 对称正定，正好命中 Cholesky 的适用条件：

\Sigma=LL^T,\qquad L\text{ 为下三角}

下面我们用一个例子来展示：

我们的目标是求解方程

$\Sigma x = v$

其中

\Sigma=\begin{bmatrix}4&2\\2&3\end{bmatrix} \;\Longrightarrow\; L=\begin{bmatrix}2&0\\1&\sqrt2\end{bmatrix}, \qquad LL^T=\begin{bmatrix}2&0\\1&\sqrt2\end{bmatrix}\begin{bmatrix}2&1\\0&\sqrt2\end{bmatrix}=\begin{bmatrix}4&2\\2&3\end{bmatrix}\;\checkmark

$v=(1,0)$

我们有

$\Sigma x = v \quad\Longleftrightarrow\quad LL^Tx=v$

设 $y=L^Tx$ ，于是要解 $\Sigma x = v$ ，等价于解：

\begin{cases} Ly = v &\text{（前代）}\\ L^T x = y &\text{（回代）} \end{cases}

第一步前代： $Ly=v$ 展开成

\begin{align} \begin{bmatrix}2&0\\1&\sqrt2\end{bmatrix}\begin{bmatrix}y_1\\y_2\end{bmatrix}&= \begin{bmatrix}1\\0\end{bmatrix}\\ &\Downarrow \\ &\begin{cases} 2y_1 = 1\\ 1\cdot y_1 + \sqrt2\,y_2 = 0 \end{cases} \end{align}

我们可以轻松解出来

\begin{cases} y_1 = 0.5\\ y_2 = -0.3536 \end{cases}

第二步回代： $L^Tx=y$ 展开（ $L^T=\begin{bmatrix}2&1\\0&\sqrt2\end{bmatrix}$ ）

\begin{cases} 2x_1 + 1\cdot x_2 = 0.5\\ \sqrt2\,x_2 = -0.3536 \end{cases}

然后算出 $x_1=0.375,x_2=-0.25$ .

简单验算一下：

\Sigma^{-1}=\frac{1}{8}\begin{bmatrix}3&-2\\-2&4\end{bmatrix}, \qquad \Sigma^{-1}\begin{bmatrix}1\\0\end{bmatrix}=\begin{bmatrix}0.375\\-0.25\end{bmatrix} \;\checkmark

可以看到，我们从头到尾没有算过 $\Sigma^{-1}$ ，但拿到了 $\Sigma^{-1}v$ .

对于 $\det$ 的计算也是一样的。我们知道，三角矩阵的行列式 = 对角元之积（这是三角矩阵的好性质之一）。再由 $\det(LL^T)=\det L\cdot\det L$ ：

\det(\Sigma)=\Big(\prod_i L_{ii}\Big)^2

$L$ 已经在第 1 步算好了，所以取行列式只是把 $d$ 个数乘一乘： $O(d)$ 。

在实际代码中取

\log

$\det$ 这个数会爆炸或消失。举个例子： $d=1000$ 、每个 $L_{ii}=2$ ，那么

\det(\Sigma)=2^{2000}\approx10^{602}

而计算机的双精度浮点数最大只能表示到约 $1.8\times10^{308}$ ——直接溢出了，变成 inf（反过来， $L_{ii}$ 小于 1 时会一路下溢成 0，信息全丢）。但取对数：

\log\det(\Sigma)=2\times1000\times\ln 2\approx1386.3

所以，在实际代码里（比如算高斯密度时）永远取 $\log$ 再往下算，这也是 “机器学习的 log-likelihood 为什么到处是 log” 的一个技术原因。

1.3 回到高斯密度：整套流程#

把上面拼起来，多元正态密度的计算变成：

\log p_X(x)=-\frac d2\log(2\pi)-\underbrace{\frac12\log\det\Sigma}_{\text{由 }L\text{ 得来}}-\underbrace{\frac12(x-\mu)^T\Sigma^{-1}(x-\mu)}_{\text{由两次三角求解得来}}

一次 Cholesky $O(d^3)$ ，之后可重复使用；
每个样本： $\log\det$ 是 $O(d)$ ，马氏距离 $O(d^2)$ 。

两个阶段

上面给出了两种不同的时间复杂度，对应计算过程中的两个阶段：

阶段	依赖什么	要做几次	成本
① 准备：对 $\Sigma$ 做 Cholesky，得到 $L$	只依赖分布（ $\Sigma$ ）	一次	$O(d^3)$
② 查询：给一个数据点 $x$ ，算 $\log p(x)$	每个样本各不相同	$N$ 次（ $N$ 个样本）	每次 $O(d^2)$

单个样本的公式是

\log p(x)=-\frac d2\log(2\pi)-\frac12\underbrace{\log\det\Sigma}_{\text{常数项}}-\frac12\underbrace{(x-\mu)^T\Sigma^{-1}(x-\mu)}_{\text{马氏距离（Mahalanobis distance）}}

第一项 $-\frac d2\log(2\pi)$ ：纯常数， $O(1)$ 。
第二项 $\log\det\Sigma=2\sum_i\log L_{ii}$ ：把 $d$ 个对角元取对数再加起来， $O(d)$ . 注意它根本不含 $x$ ——对所有样本都一样，所以严格说算一次就够了（之后缓存成一个数，每次查表 $O(1)$ ）。
第三项（马氏距离） 才是真正”每个样本都要重算”的部分：

(x-\mu)^T\Sigma^{-1}(x-\mu)

怎么算？不碰 $\Sigma^{-1}$ ，而是解两次三角方程：

\text{\textcircled{1}令 } b=x-\mu \;\to\; \text{\textcircled{2}解 } Ly=b \;(O(d^2)) \;\to\; \text{\textcircled{3}解 } L^Tu=y\;(O(d^2)) \;\to\; \text{\textcircled{4}算 } b^Tu\;(O(d))

合计 $O(d^2)$ 。每来一个新样本，就重复一次这个 $O(d^2)$ 的流程。

1.4 当 $d$ 极大、连 $O(d^3)$ 都受不了时#

这是真正的”大维度”对策，和我们后面将要介绍的 SVD / 低秩相关：

(a) 低秩 + 对角结构（Woodbury 恒等式）

若

\Sigma = D + UU^T,\qquad D\text{ 对角},\ U\in\mathbb R^{d\times k},\ k\ll d

则用 Woodbury 矩阵恒等式：

\Sigma^{-1}=D^{-1}-D^{-1}U\big(I+U^T D^{-1}U\big)^{-1}U^T D^{-1}

\log\det\Sigma=\sum_i\log D_{ii}+\log\det\!\big(I+U^T D^{-1}U\big)

所有昂贵的运算都只作用在小矩阵 $I+U^TD^{-1}U$ （ $k\times k$ ）上，代价降到 $O(dk^2)$ 。

(b) 稀疏 / 带状 / 层次结构

若 $\Sigma^{-1}$ （精度矩阵）稀疏（如高斯随机场、稀疏高斯过程），就用稀疏 Cholesky / 稀疏 LU，代价 $\approx O(d^{1.5})$ 甚至 $O(d)$ 。

(c) 随机化方法

用随机化 SVD 求低秩近似——这正是 Lec2 结尾预告的”randomized SVD”。

(d) 不要 $\det$ 本身，只要它的对数 —— 之前已说明。

2. 如何高效计算矩阵-向量乘法#

2.1 稠密一般矩阵： $O(d^2)$ 是下界#

对一般稠密 $A\in\mathbb R^{d\times d}$ ， $y=Av$ 需要

y_i=\sum_{j=1}^d A_{ij}v_j

你必须碰到矩阵的每一个元素，所以 $\text{FLOPs}=O(d^2)$ 是渐近最优的，无法再降。

\text{FLOPs}

FLOPs = Floating-point Operations（浮点运算次数），即做了多少次的浮点数加减乘除。

在上面的例子中，计算某个 $y_i$ 需要 $d$ 次乘法 + $d$ 次加法，一共有 $d$ 个 $y_i$ ，所以总运算量约为 $d\times 2d=2d^2$ ，用大 $O$ 记就是 $O(d^2)$ .

那”高效”从哪来？硬件层面：

矩阵-向量乘法算术强度低（每读一个矩阵元素只做一次乘法+一次加法），是内存带宽瓶颈；
优化手段：cache 分块（blocking）、向量化、利用对称性（ $\Sigma$ 对称 → 只存一半）；
批量做多个向量：把 $v_1,\ldots,v_N$ 拼成矩阵，变成 Level 3 矩阵-矩阵乘法，算术强度大幅提升，性能天差地别。

2.2 利用结构：真正能把指数降下来的地方#

矩阵结构	朴素	结构化后	手段
一般稠密	$O(d^2)$	$O(d^2)$	BLAS、分块、向量化
对称	$O(d^2)$	$\tfrac12 O(d^2)$	只存/算一半
低秩 $A=UV^T$	$O(d^2)$	$O(dk)$	先算 $U^Tv$ ，再乘 $U$
带状（带宽 $w$ ）	$O(d^2)$	$O(dw)$	只遍历带宽
稀疏（ $nnz$ 非零）	$O(d^2)$	$O(nnz)$	稀疏存储 + SpMV
Toeplitz/循环	$O(d^2)$	$O(d\log d)$	用 FFT 卷积
层次/H-矩阵	$O(d^2)$	$O(d\log d)$ 或 $O(d)$	低秩分块

回到高斯： $\Sigma^{-1}v$ 之所以能”高效”，本质就是把 $\Sigma^{-1}$ 看成”三角分解的作用”而不是一个显式矩阵——用两次 $O(d^2)$ 的三角求解，取代一次 $O(d^3)$ 求逆 + 一次 $O(d^2)$ 乘法。

2.3 近似计算#

若允许近似：

随机化投影 / Johnson–Lindenstrauss：把 $d$ 维投影到低维再做；
随机化 SVD（后续内容）：用少量矩阵-向量乘逼近低秩结构；
采样/草图（sketching）： $O(d)$ 或 $O(d\log d)$ 。

音乐

音乐

1. 如何高效计算 $\det(\Sigma)$ 和 $\Sigma^{-1}$ #

1.1 “直接求逆” 的坏处#

1.1.1 为什么显式求逆是 $O(d^3)$ ？#

1.2 Cholesky 分解#

1.3 回到高斯密度：整套流程#

1.4 当 $d$ 极大、连 $O(d^3)$ 都受不了时#

2. 如何高效计算矩阵-向量乘法#

2.1 稠密一般矩阵： $O(d^2)$ 是下界#

2.2 利用结构：真正能把指数降下来的地方#

2.3 近似计算#

文章分享

评论区

音乐

音乐

如何高效计算高维行列式与逆矩阵

1. 如何高效计算 det⁡(Σ)\det(\Sigma)det(Σ) 和 Σ−1\Sigma^{-1}Σ−1#

1.1 “直接求逆” 的坏处#

1.1.1 为什么显式求逆是 O(d3)O(d^3)O(d3)？#

1.2 Cholesky 分解#

1.3 回到高斯密度：整套流程#

1.4 当 ddd 极大、连 O(d3)O(d^3)O(d3) 都受不了时#

2. 如何高效计算矩阵-向量乘法#

2.1 稠密一般矩阵：O(d2)O(d^2)O(d2) 是下界#

2.2 利用结构：真正能把指数降下来的地方#

2.3 近似计算#

文章分享

评论区

音乐

1. 如何高效计算 $\det(\Sigma)$ 和 $\Sigma^{-1}$ #

1.1.1 为什么显式求逆是 $O(d^3)$ ？#

1.4 当 $d$ 极大、连 $O(d^3)$ 都受不了时#

2.1 稠密一般矩阵： $O(d^2)$ 是下界#