貝葉斯機器學習：最大熵及高斯分布

高斯分布，也被稱為正態分布，廣泛應用于連續型隨機變量分布的模型中。高斯分布可以從多個不同的角度來理解。例如，對于一個一元實值向量，使得熵取得最大值的是高斯分布。這個性質對于多元高斯分布也成立。當我們考慮多個隨機變量之和的時候，也會產生高斯分布。觀察式多元高斯分布的形式，考慮其中在指數位置上出現的二次型(x - mu)^T∑^{-1}(x - mu)。由于協方差矩陣∑是對稱矩陣，那么∑^{-1}也是對稱矩陣。我們假定∑是正定的，那么∑^{-1}也是正定的。于是，該二次型為x到mu的馬?距離（Mahalanobis distance）Delta的平方。當∑是單位陣時，就變成了歐氏距離。

高斯分布^[1]，也被稱為正態分布，廣泛應用于連續型隨機變量分布的模型中。對于一元變量$x$的情形。高斯分布可以寫成下列的形式：

\[\mathcal{N}(x\mid \mathcal{\mu}, \sigma^2) = \frac{1}{(2\pi\sigma^2)^{\frac{1}{2}}}\exp\left\{-\frac{1}{2\sigma^2}(x - \mu)^2\right\} \]

其中$\mu$是均值，$\sigma^2$是方差。對于$D$維向量$\boldsymbol{x}$，多元高斯分布的形式為：

\[\mathcal{N}(\boldsymbol{x}\mid \boldsymbol{\mu}, \boldsymbol{\Sigma}) = \frac{1}{(2\pi)^{\frac{D}{2}}}\frac{1}{(\mathrm{det}\boldsymbol{\Sigma})^{\frac{1}{2}}}\exp\left\{-\frac{1}{2}(\boldsymbol{x} - \boldsymbol{\mu})^T\boldsymbol{\Sigma}^{-1}(\boldsymbol{x} - \boldsymbol{\mu})\right\}\tag{1} \]

其中，$\boldsymbol{\mu}$是一個$D$維均值向量，$\boldsymbol{\Sigma}$是一個$D\times D$的協方差矩陣（covariance matrix），$\mathrm{det}\boldsymbol{\Sigma}$是$\boldsymbol{\Sigma}$的行列式。

1 高斯分布的物理意義

1.1 做為最大熵分布的高斯分布

高斯分布可以從多個不同的角度來理解。例如，對于一個一元實值向量，使得熵取得最大值的是高斯分布。這個性質對于多元高斯分布也成立。

熵（entropy） 的概念最早起源于物理學，是在熱力學平衡的背景中介紹的。后來，熵稱為描述統計力學中的無序程度的度量。在統計力學中，玻爾茲曼熵（Boltzmann's entropy）^[2][3]的定義為：

\[S \equiv k\ln\Omega \]

這里$\Omega$為系統宏觀態的重數，$k$為玻爾茲曼常量。

注關于統計力學的術語我們用一個例子^[4]來簡要做一下介紹。考慮將3個不同的球放進3個不同的箱子中，全部的$3^3 = 27$種可能放法如下圖所示：

這27種不同結果中的每一種都稱為微觀態（microstate）。通常，在統計力學中，為了知道系統的微觀態，我們必須清楚每個粒子的狀態，在我們這個例子中是每個球劃分到箱子中去的狀態。如果更一般地指定狀態——例如僅僅說各個箱子里有多少個球，我們稱它為宏觀態（macrostate）。當然，如果知道一個系統的微觀態（比如$\{ab\space\mid\space\space c\mid\space-\space\}$），那么我們肯定也能知道它的宏觀態（比如箱子中的球個數分別為2、1、0）。但反過來卻不行：知道箱子中的球個數分別為2、1、0并沒有告訴我們每個球的狀態，因為有3個微觀態都對應于這個宏觀態。對應于給定宏觀態的微觀態數量稱為該宏觀態的重數（multiplicity），在這種情況下為3。
更一般地，考慮將$N$個球放進$M$個箱子里，各個箱子里有$n_1, n_2,\cdots, n_M$個球時宏觀態的重數為我們在博客《概率論沉思錄：初等抽樣論》中提到的多項式系數：

\[\Omega(n_1, n_2, \cdots, n_M) = \frac{N!}{\prod_{i}^Mn_i!} \]
其中$\sum_{i=1}^Mn_i = N$。

對于將$N$個不同的球放進$M$個不同的箱子所形成的這樣一個系統，各個箱子里有$n_1, n_2,\cdots, n_M$個球時系統的熵為

\[S = k\ln\frac{N!}{\prod_{i}^Mn_i!} = k\ln N! - k\sum_{i=1}^M\ln n_i! \]

現在我們在此基礎上乘以一個縮放參數$\frac{1}{N}$，忽略掉常數$k$，并考慮$N\gg 1$，根據斯特林近似（Stirlings approximation） 我們有$\ln N!\approx N\ln N - N$，于是

\[\begin{aligned} \frac{1}{N}S &\approx \frac{1}{N}\left(N\ln N - N\right) - \frac{1}{N}\sum_{i=1}^M\left[n_i\ln n_i - n_i\right]\\ &= \ln N - \sum_{i=1}^M \frac{n_i}{N}\ln n_i\\ &= - \sum_{i=1}^M \left(\frac{n_i}{N}\right)\ln \left(\frac{n_i}{N}\right) \end{aligned} \]

注斯特林近似為：

\[ N!\approx \sqrt{2\pi N}\left(\frac{N}{e}\right)^N \]
當$N \gg 1$時，這個公式十分準確。該公式可以采用下列的直觀方式進行理解。$N!$是從$1$到$N$這$N$個因子的乘積，一個十分粗略的估計是把每個因子都替換成$N$，這就是$N!\approx N^N$。這是一個過高的估計，因為幾乎所有的因子都比$N$小，平均下來每個因子都大了大約$e$倍，也就是說

\[ N!\approx (\frac{N}{e})^N \]
這仍比$N!$差了一個大數字的因子，大約是$\sqrt{2\pi N}$。但若$N$是一個大數字，那么$N!$會是一個非常大的數字，所以這個（大數字）因子就可以被省略了。如果我們只關心$N!$的對數，通常上式就已經足夠準確：

\[\ln N! \approx N\ln N - N \]

取極限$N\rightarrow \infty$，并保持比值$\frac{n_i}{N}$固定，我們可以得到

\[-\lim_{N\rightarrow \infty}\sum_{i=1}^M \left(\frac{n_i}{N}\right)\ln \left(\frac{n_i}{N}\right) = -\sum_{i=1}^M p_i\ln p_i \]

這里$p_i = \lim_{N\rightarrow \infty}(\frac{n_i}{N})$為給定系統宏觀態時，一個球被分配到第$i$個箱子的概率。對于一個球，它被分配到的箱子可以表述成一個離散隨機變量$X$（$X$一共有$M$個狀態），且$p(X = x_i) = p_i\space (i = 1, \cdots, M)$。這樣，隨機變量$X$的熵就為

\[H[p] = -\sum_{i=1}^M p(x_i)\ln p(x_i) \]

如果分布$p(x_i)$在幾個值周圍有尖銳的峰值，熵就會相對較低。如果分布$p(x_i)$相對平衡地跨過許多值，那么熵就會相對較高。下圖展示了兩個概率分布在30個箱子上的直方圖：

可以看到，熵值越大，$H$越寬（最大的熵值產生于均勻分布，此時的熵值為$H = -\ln(\frac{1}{30})\approx 3.40$）。

由于$0\leqslant p_i\leqslant 1$，因此熵是非負的（由于當$x\rightarrow 0$時，$x\ln x\rightarrow 0$，我們約定$0\ln 0 = 0$）。當$p_i = 1$且所有其他的$p_{j\neq i}=0$時，熵取得最小值$0$。在$\sum_i p(x_i) = 1$（概率歸一化）的約束下，使用拉格朗日乘數法可以找到熵的最大值。因此，我們要最大化

\[\tilde{H} = -\sum_{i=1}^M p(x_i)\ln p(x_i) + \lambda\left(\sum_{i=1}^M p(x_i) - 1\right) \]

令$\frac{\mathrmw0obha2h00(\tilde{H})}{\mathrmw0obha2h00p(x_i)}=0$，解得$p(x_i) = e^{\lambda - 1}$，然后代回約束條件$\sum_i p(x_i)$，得到$\lambda = 1 - \ln M$，于是最優解為$p(x_i) = \frac{1}{M}$（此時所有的$p(x_i)$都相等），此時熵取得最大值$\ln M$。

我們可以把熵的定義擴展到連續隨機變量$X$的概率分布$p(x)$，方法如下。首先把連續隨機變量$X$的取值范圍切分成寬度為$\Delta$的小段區間。根據均值定理（mean value theorem）^[5]，對于每個這樣的小段區間$\Delta$，一定在區間中存在一個值$x_i$使得

\[\int_{i\Delta}^{(i + 1)\Delta}p(x)\mathrmw0obha2h00x = p(x_i)\Delta \]

現在我們可以這樣量化連續隨機變量$X$：只要$X$落在第$i$個區間中，我們就把$X$賦值為$x_i$。因此觀察到值$x_i$的概率為$p(x_i)\Delta$。這樣就變成了離散的分布，這種情況下熵的形式為：

\[H_{\Delta} = -\sum_{i=1}\left(p(x_i)\Delta\right) \ln \left(p(x_i)\Delta\right) = -\sum_{i=1}p(x_i)\Delta \ln p(x_i) - \ln \Delta \]

推導時我們使用了$\sum_ip(x_i)\Delta$，這是因為$\sum_{i}p(x_i)\Delta = \sum_i\int_{i\Delta}^{(i + 1)\Delta}p(x)\mathrmw0obha2h00x = 1$。我們現在省略上式右側的第二項$-\ln \Delta$，然后考慮極限$\Delta\rightarrow 0$，此時有：

\[\lim_{\Delta\rightarrow 0}\left\{-\sum_{i=1}p(x_i)\Delta \ln p(x_i)\right\} = -\int p(x)\ln p(x) \mathrmw0obha2h00x \]

其中，右側的量被稱為微分熵（differential entropy）。我們看到，熵的離散形式與連續形式的差是$\ln\Delta$，這在極限$\Delta\rightarrow 0$的情形下發散。這反映出一個事實：具體化一個連續隨機變量需要大量的比特位。和離散情形類似，我們將微分熵記為關于概率密度函數$p(x)$的函數：

\[H[p] = -\int p(x)\ln p(x) \mathrmw0obha2h00x \]

在離散分布的情況下，我們看到最大熵對應于變量所有可能狀態的均勻分布。現在讓我們考慮連續隨機變量的最大熵。為了讓這個最大值有一個合理的定義，除了保留歸一化的約束之外，還有必要限制$p(x)$的均值（一階矩）和方差（二階矩）。之所以需要限制其方差，是因為當方差增大時，熵也會無限制地增加，因此除非我們給定固定的方差$\sigma^2$，否則尋找哪一個分布有最大熵這個問題是沒有意義的。之所以需要限制其均值，是因為在不改變熵的條件下一個分布可以被隨意地改變，因此我們需要再加一個均值為$\mu$的約束以獲得一個唯一的解^[6]。綜上，我們最大化微分熵的時候要施加下面三個約束：

\[\begin{aligned} \int p(x) \mathrmw0obha2h00x &= 1\\ \int xp(x) \mathrmw0obha2h00x &= \mu\\ \int (x - \mu)^2p(x) \mathrmw0obha2h00x &= \sigma^2 \end{aligned} \]

那么這個問題的拉格朗日泛函如下：

\[J(p) = H[p] + \lambda_1 (\int p(x) \mathrmw0obha2h00x - 1) + \lambda_2 (\int xp(x) \mathrmw0obha2h00x - \mu) + \lambda_3 (\int (x - \mu)^2p(x) \mathrmw0obha2h00x - \sigma^2) \]

根據變分法，可以得到泛函導數為：

\[\frac{\delta J}{\delta p(x)} = -(\ln p(x) + 1) + \lambda_1 + \lambda_2 x + \lambda_3(x - \mu)^2 \]

令$\frac{\delta J}{\delta p(x)} = 0$（對$\forall x$），解得$p(x) = \exp\{-1 + \lambda_1 + \lambda_2x + \lambda_3(x - \mu)^2\}$。

注函數$f$的函數被稱為泛函（functional）$J(f)$。正如許多情況下對一個函數求關于以向量的各元素為變量的偏導數一樣，我們可以使用泛函導數（functional derivative），即在任意特定的$x$值，對一個泛函$J(f)$求關于函數$f(x)$的導數，這也被稱為變分導數（variational derivative）。泛函$J$的關于函數$f$在點$x$處的泛函導數被記作$\frac{\delta J}{\delta f(x)}$。

對于可微分函數$f(x)$以及帶有連續導數的可微分函數$g(y, x)$（其中$y=f(x)$），有下列結論：

\[\frac{\delta}{\delta f(x)} \int g(f(x), x)\mathrmw0obha2h00x = \frac{\partial}{\partial y} g(f(x), x) \]
為了使上述等式更加直觀，我們可以把$f(x)$看做是一個有著無窮不可數多個元素的向量。在這里，這種關系式中描述的對于特定的$x$，關于$f(x)$的泛函導數可以類比為對于特定的下標$i$，關于向量$\boldsymbol{\theta}\in \mathbb{R}^n$的第$i$個元素的偏導數：

\[\frac{\partial}{\partial \theta_i}\sum_j g(\theta_j, j) = \frac{\partial}{\partial \theta_i} g(\theta_i, i) \]
為了關于一個向量優化某個函數，我們可以求出這個函數關于這個向量的梯度，然后找到使這個梯度中每一個元素都為0的點。類似地，我們可以通過尋找一個函數使得泛函導數在每個點上都等于0，從而來優化一個泛函。

為了滿足所有的約束，我們可以令$\lambda_1 = 1 - \ln \sigma \sqrt{2\pi}$，$\lambda_2 = 0$，$\lambda_3 = - \frac{1}{2\sigma^2}$，從而得到

\[p(x) = \frac{1}{(2\pi \sigma^2)^{\frac{1}{2}}}\exp\left\{ - \frac{(x - \mu)^2}{2\sigma^2}\right\} \]

因此最大化微分熵的分布是高斯分布$\mathcal{N}(x\mid \mathcal{\mu}, \sigma^2)$。這也是當我們不知道真實的分布時，往往使用高斯分布的一個原因。因為高斯分布擁有最大的熵，我們通過這個假定來保證了最小可能量的結構。

注若$X\in [a, b]$，無其它約束條件，則此時最大熵分布就是該區間上的均勻分布，這里可以和$X$為離散隨機變量的情形聯系起來。

如果我們求高斯分布的微分熵，我們會得到：

\[\begin{aligned} H(p) &= -\int p(x) \ln p(x) \mathrmw0obha2h00x\\ &= -\int \frac{1}{(2\pi \sigma^2)^{\frac{1}{2}}}\exp\left\{ - \frac{(x - \mu)^2}{2\sigma^2}\right\}\left[\ln \frac{1}{(2\pi \sigma^2)^{\frac{1}{2}}} - \frac{(x - \mu)^2}{2\sigma^2}\right]\\ &= -\ln \frac{1}{(2\pi \sigma^2)^{\frac{1}{2}}} + \frac{1}{2\sigma^2} \int \underbrace{\frac{1}{(2\pi \sigma^2)^{\frac{1}{2}}} \exp\left\{ - \frac{(x - \mu)^2}{2\sigma^2}\right\} (x - \mu)^2 \mathrmw0obha2h00x}_{\sigma^2}\\ &= \frac{1}{2}\left\{\ln (2\pi\sigma^2) + 1\right\} \end{aligned} \]

我們可以看到熵隨著分布寬度（即$\sigma^2$）的增加而增加。這個結果也表明，與離散熵不同，微分熵可以為負，因為對于上式而言，當$\sigma^2 < \frac{1}{2\pi e}$時，$H(p) < 0$。

1.2 做為隨機變量和分布的高斯分布

當我們考慮多個隨機變量之和的時候，也會產生高斯分布。根據林德伯格所證明的中心極限定理（central limit theorem）（為我們在博客《概率論沉思錄：初等假設檢驗》中提到過的伯努利試驗的棣莫弗-拉普拉斯極限定理的推廣），設$\{X_k\}_{k=1}^N$是相互獨立且具有共同分布的隨機變量序列，假定$\mu=\mathbb{E}[X_k]$和$\sigma^2 = \mathrm{Var}[X_k]$都存在，并令$S_N = X_1 + \cdots + X_N$，則當$N\rightarrow \infty$時，有$S_N \rightarrow \mathcal{N}(N\mu, \sigma\sqrt{N})$（依分布）。

例如，考慮$N$個隨機變量$X_1, \cdots, X_N$，每一個都是區間$[0, 1]$上的均勻分布，然后考慮$N$個隨機變量的均值$\frac{S_N}{N} = \frac{1}{N}(X_1 + \cdots + X_N)$的分布。對于大的$N$，這個分布趨向于高斯分布，如下圖所示：

在實際應用中，隨著$N$的增加，分布會很迅速收斂為高斯分布。

2 高斯分布的性質

2.1 高斯分布的幾何性質

觀察式$(1)$中的多元高斯分布的形式，考慮其中在指數位置上出現的二次型

\[(\boldsymbol{x} - \boldsymbol{\mu})^T\boldsymbol{\Sigma}^{-1}(\boldsymbol{x} - \boldsymbol{\mu})\tag{2} \]

由于協方差矩陣$\mathbf{\Sigma}$是對稱矩陣，那么$\mathbf{\Sigma}^{-1}$也是對稱矩陣（注意$\mathbf{\Sigma}^{-1}$的第$j, i$個元素為$(-1)^{i + j}\left[\mathrm{det}\mathbf{\Sigma}_{ij} / \mathrm{det}\mathbf{\Sigma}\right]$，其中$\mathbf{\Sigma}_{ij}$是從$\mathbf{\Sigma}$中除去第$i$行和第$j$列后得到的矩陣，而對于對稱矩陣$\mathbf{\Sigma}$有$\mathrm{det}\mathbf{\Sigma}_{ij} = \mathrm{det}\mathbf{\Sigma}_{ji}$）。我們假定$\mathbf{\Sigma}$是正定的，那么$\mathbf{\Sigma}^{-1}$也是正定的（后面我們會證明）。于是，式$(2)$為$\boldsymbol{x}$到$\boldsymbol{\mu}$的馬?距離（Mahalanobis distance）$\Delta$的平方。當$\boldsymbol{\Sigma}$是單位陣時，就變成了歐氏距離。

注給定$D\times D$的對稱正定矩陣$\mathbf{A}$，則從點$\boldsymbol{x}$到原點的馬氏距離$\Delta$由正定二次型$0 < \Delta^2 = \boldsymbol{x}^T\mathbf{A}\boldsymbol{x}$（$\boldsymbol{x} \neq \boldsymbol{0}$）確定；反過來，一個正定二次型$\boldsymbol{x}^T\mathbf{A}\boldsymbol{x} > 0$（$\boldsymbol{x} \neq \boldsymbol{0}$）可以解釋為點$\boldsymbol{x}$到原點的馬氏距離$\Delta$的平方^[7]。此外，從點$\boldsymbol{x}$到任意一固定點$\boldsymbol{\mu}$的馬氏距離$\Delta$的平方由公式$(\boldsymbol{x} - \boldsymbol{\mu})^T\mathbf{A}(\boldsymbol{x} - \boldsymbol{\mu})$給出。

馬氏距離可以理解為在新的基底（也即$\mathbf{A}$的規范正交特征向量組成的基底）下兩點之間的距離。例如，設$D=2$，則到原點的馬氏距離為常數$c$的點$\boldsymbol{x} = (x_1, x_2)^T$滿足$\boldsymbol{x}^T\mathbf{A}\boldsymbol{x} = c^2$。對$\mathbf{A}$進行譜分解得：

\[\mathbf{A} = \lambda_1 \boldsymbol{e}_1\boldsymbol{e}_1^T + \lambda_2\boldsymbol{e}_2\boldsymbol{e}_2^T \]
所以

\[\begin{aligned} \boldsymbol{x}^T\mathbf{A}\boldsymbol{x} &= \boldsymbol{x}^T(\lambda_1 \boldsymbol{e}_1\boldsymbol{e}_1^T + \lambda_2\boldsymbol{e}_2\boldsymbol{e}_2^T)\boldsymbol{x}\\ & = \lambda_1 (\boldsymbol{e}_1^T\boldsymbol{x})^2 + \lambda_2 (\boldsymbol{e}_2^T\boldsymbol{x})^2 \end{aligned} \]
令$y_1 = \boldsymbol{e}_1^T\boldsymbol{x}, y_2 = \boldsymbol{e}_2^T\boldsymbol{x}$，則$\boldsymbol{y} = (y_1, y_2)^T$可視為點$\boldsymbol{x}=(x_1, x_2)^T$在新的基底$\boldsymbol{e}_1, \boldsymbol{e}_2$下的坐標，基變換關系為$\left(\begin{matrix} y_1\\ y_2 \end{matrix}\right) = \left(\begin{array}{c:c}\boldsymbol{e}_1 & \boldsymbol{e}_2\end{array}\right)^T\left(\begin{matrix} x_1\\ x_2 \end{matrix}\right) = \left(\begin{matrix} \boldsymbol{e}_1^T\\ \hdashline \boldsymbol{e}_2^T \end{matrix}\right)\left(\begin{matrix} x_1\\ x_2 \end{matrix}\right)$。此時，$\boldsymbol{y}=(y_1, y_2)^T$在以原點為中心的橢圓上，其方程為：

\[ \begin{aligned} \lambda_1 y_1^2 + \lambda_2 y_2^2 = c^2 \end{aligned} \]
該橢圓的軸即為做為新基底的規范正交特征向量$\boldsymbol{e}_1, \boldsymbol{e}_2$。我們可以發現將$\boldsymbol{x} = c\lambda_1^{-1/2}\boldsymbol{e}_1$代入橢圓方程得到$\boldsymbol{x}^T\mathbf{A}\boldsymbol{x} = \lambda_1 (c\lambda_1^{-1/2}\boldsymbol{e}_1^T\boldsymbol{e}_1)^2=c^2$，可見橢圓沿$\boldsymbol{e}_1$方向的半軸長為$c\lambda_1^{-1/2}$；同理，$\boldsymbol{x} = c\lambda_2^{-1/2}\boldsymbol{e}_2$也給出了橢圓沿$\boldsymbol{e}_2$方向的半軸長$c\lambda_2^{-1/2}$。到原點的馬氏距離為常數$c$的點的位置如下圖所示（$D = 2, \lambda_1 < \lambda_2$）：

如果$D > 2$，到原點的馬氏距離為常數$c = \sqrt{\boldsymbol{x}^T\mathbf{A}\boldsymbol{x}}$的點$\boldsymbol{x} = (x_1, \cdots, x_D)$在橢球面$\lambda_1 (\boldsymbol{e_1}^T\boldsymbol{x})^2 + \cdots + \lambda_D (\boldsymbol{e}_D^T\boldsymbol{x}^T)^2 = c^2$上，其軸由$\mathbf{A}$的規范正交特征向量給出。沿$\boldsymbol{e}_i$方向的半軸長為$\frac{c}{\sqrt{\lambda_i}}, i = 1, 2, \cdots, p$，其中$\lambda_1, \cdots, \lambda_D$是$\mathbf{A}$的特征值。

下圖直觀地說明了馬氏距離相比歐氏距離具有優越性的情況：

該圖描述了重心（樣本均值）在點$Q$的一組點。$P$到$Q$的歐氏距離大于$O$到$Q$的歐氏距離。然而，$P$點卻比原點$O$更像是屬于這一組點內的點。如果我們用馬氏距離來度量距離，則$P$距離$Q$就會比$O$距離$Q$要近了。

至于馬氏距離中的矩陣$\mathbf{A}$如何確定則和樣本各維度的標準差有關。馬氏距離相當于將樣本各維度除以樣本標準差之后（也即“標準化”之后），再采用歐氏距離的公式進行計算，感興趣的讀者可以閱讀《Applied Multivariate Statistical Analysis》1.5節。

如果令多元高斯分布中的馬氏距離的平方$\Delta^2 = (\boldsymbol{x} - \boldsymbol{\mu})^T\boldsymbol{\Sigma}^{-1}(\boldsymbol{x} - \boldsymbol{\mu})$為常數$c^2$，那么多元高斯分布的概率密度也為常數，這也就意味著到均值$\boldsymbol{\mu}$的馬氏距離相等的點擁有相同的賦概。我們前面提到過，這些點在一個橢球面上，我們將其稱為輪廓線（contours）：

\[\begin{aligned} \quad\text{常數概率密度輪廓線} &= \{\boldsymbol{x}\mid (\boldsymbol{x} - \boldsymbol{\mu})^T \mathbf{\Sigma}^{-1} (\boldsymbol{x} - \boldsymbol{\mu}) = c^2\}\\ &= \text{中心在}\boldsymbol{\mu}的橢球面 \end{aligned} \]

為了確定橢球面的軸方向和半軸長，我們對協方差矩陣的逆矩陣$\mathbf{\Sigma}^{-1}$進行譜分解：

\[\mathbf{\Sigma}^{-1} = \sum_{i=1}^D \frac{1}{\lambda_i}\boldsymbol{e}_i\boldsymbol{e}_i^T\tag{3} \]

其中$\lambda_1, \cdots, \lambda_D$為$\mathbf{\Sigma}$的特征值，$\boldsymbol{e}_1, \cdots, \boldsymbol{e}_D$為與之相伴的規范正交特征向量。

注關于$\mathbf{\Sigma}^{-1}$，我們有下列結論：
$\quad$ 結論若$\mathbf{\Sigma}$是正定的，其逆為$\mathbf{\Sigma}^{-1}$，則$$
\mathbf{\Sigma}\boldsymbol{e} = \lambda \boldsymbol{e}\Rightarrow\mathbf{\Sigma}^{-1}\boldsymbol{e} = \left(\frac{1}{\lambda}\right)\boldsymbol{e}$$

所以$\mathbf{\Sigma}$的特征值-特征向量對$(\lambda, \boldsymbol{e})$對應于$\mathbf{\Sigma}^{-1}$的特征值-特征向量對$(\frac{1}{\lambda}, \boldsymbol{e})$。且$\mathbf{\Sigma}^{-1}$也是正定的。

$\quad$ 證明對于正定的$\mathbf{\Sigma}$以及一個特征向量$\boldsymbol{e}\neq \boldsymbol{0}$，我們有$0 < \boldsymbol{e}^T\mathbf{\Sigma}\boldsymbol{e} = \boldsymbol{e}^T(\mathbf{\Sigma}\boldsymbol{e}) = \boldsymbol{e}^T(\lambda\boldsymbol{e}) = \lambda\boldsymbol{e}^T\boldsymbol{e} = \lambda$。而且$\boldsymbol{e} = \mathbf{\Sigma}^{-1}(\mathbf{\Sigma}\boldsymbol{e}) = \mathbf{\Sigma}^{-1}(\lambda\boldsymbol{e}) = \lambda\mathbf{\Sigma}^{-1}\boldsymbol{e}$，用$\lambda > 0$除，得到$\mathbf{\Sigma}^{-1}\boldsymbol{e} = \left(\frac{1}{\lambda}\right)\boldsymbol{e}$。于是，$(\frac{1}{\lambda}, \boldsymbol{e})$是$\mathbf{\Sigma}^{-1}$的一對特征值-特征向量。因此，我們對任意$D$維向量$\boldsymbol{x}\neq \boldsymbol{0}$有

\[\begin{aligned} \boldsymbol{x}^T\mathbf{\Sigma}^{-1}\boldsymbol{x} &= \boldsymbol{x}^T\left(\sum_{i=1}^D \frac{1}{\lambda_i}\boldsymbol{e}_i\boldsymbol{e}_i^T\right)\boldsymbol{x}\\ &= \sum_{i=1}^D (\frac{1}{\lambda_i}) (\boldsymbol{e}_i^T\boldsymbol{x})^2 > 0 \end{aligned} \]
由此得出$\mathbf{\Sigma}^{-1}$是正定的。

將式$(3)$的$\mathbf{\Sigma}^{-1}$的譜分解結果代入式$(2)$所示的二次型中，有

\[\begin{aligned} c^2 &= (\boldsymbol{x} - \boldsymbol{\mu})^T \mathbf{\Sigma}^{-1} (\boldsymbol{x} - \boldsymbol{\mu})\\ &= (\boldsymbol{x} - \boldsymbol{\mu})^T \left(\sum_{i=1}^D \frac{1}{\lambda_i}\boldsymbol{e}_i\boldsymbol{e}_i^T\right) (\boldsymbol{x} - \boldsymbol{\mu})\\ &= \sum_{i=1}^D \frac{\left(\boldsymbol{e}_i^T(\boldsymbol{x} - \boldsymbol{\mu})\right)^2}{\lambda_i}\\ \end{aligned} \]

令$y_i = \boldsymbol{e}_i^T(\boldsymbol{x} - \boldsymbol{\mu})$，則$\boldsymbol{y}$可視為點$\boldsymbol{x} - \boldsymbol{\mu}$在新的基底$\{\boldsymbol{e}_i\}$下的坐標，基變換關系為$\left(\begin{matrix} y_1\\ \vdots\\ y_D \end{matrix}\right) = \left(\begin{array}{c:c:c}\boldsymbol{e}_1& \cdots &\boldsymbol{e}_D\end{array}\right)^T(\boldsymbol{x} - \boldsymbol{\mu}) = \left(\begin{matrix} \boldsymbol{e}_1^T\\ \hdashline \vdots\\ \hdashline \boldsymbol{e}_D^T \end{matrix}\right)(\boldsymbol{x} - \boldsymbol{\mu})$。此時，$\boldsymbol{y}=(y_1, \cdots, y_D)^T$在以$\boldsymbol{\mu}$為中心的橢球面上，其方程為：

\[ \begin{aligned} \sum_{i=1}^D \frac{y_i^2}{\lambda_i} = c^2 \end{aligned} \]

該橢球面的軸即為做為新基底的規范正交特征向量$\boldsymbol{e}_1, \cdots, \boldsymbol{e}_D$。我們可以發現將$\boldsymbol{x}_i = \boldsymbol{\mu} + c\lambda_i^{1/2}\boldsymbol{e}_i$代入橢圓方程得到$(\boldsymbol{x} - \boldsymbol{\mu})^T \mathbf{\Sigma}^{-1} (\boldsymbol{x} - \boldsymbol{\mu}) = \frac{\left(c\lambda_i^{1/2}\boldsymbol{e}_i^T\boldsymbol{e}_i\right)^2}{\lambda_i} = c^2$，可見橢圓沿$\boldsymbol{e}_i$方向的半軸長為$c\lambda_i^{1/2}$（$i = 1, \cdots, D$）。在$D=2$的情況下，該橢球面為中心為$\boldsymbol{\mu}$的橢圓，如下圖所示（$\lambda_2 < \lambda_1$）：

現在我們考慮多元高斯分布是否是歸一化的。我們之前使用了變量替換$\boldsymbol{y} = T(\boldsymbol{x}) = \left(\begin{matrix} \boldsymbol{e}_1^T\\ \hdashline \vdots\\ \hdashline \boldsymbol{e}_D^T \end{matrix}\right)(\boldsymbol{x} - \boldsymbol{\mu})$，于是在新的基底下的積分公式可以經由變量替換^[8][9]表示為

\[\int \mathcal{N}(\boldsymbol{y}\mid \boldsymbol{0}, \boldsymbol{\Sigma})\mathrmw0obha2h00\boldsymbol{y} = \int \mathcal{N}\left(T(\boldsymbol{x})\mid \boldsymbol{0}, \boldsymbol{\Sigma}\right)|\mathrm{det}\mathbf{J}|\mathrmw0obha2h00\boldsymbol{x} \]

這里$\mathbf{J}$為仿射變換$T$的Jacobian矩陣：

\[\mathbf{J} = \left(\begin{matrix} &\frac{\partial T_1(\boldsymbol{x})}{\partial x_1} &\cdots &\frac{\partial T_1(\boldsymbol{x})}{\partial x_D}\\ &\vdots & &\vdots\\ &\frac{\partial T_D(\boldsymbol{x})}{\partial x_1} &\cdots &\frac{\partial T_D(\boldsymbol{x})}{\partial x_D} \end{matrix}\right) = \left(\begin{matrix} &e_{11} &\cdots &e_{1D}\\ &\vdots & &\vdots\\ &e_{D1} &\cdots &e_{DD} \end{matrix}\right) = \left(\begin{matrix} \boldsymbol{e}_1^T\\ \hdashline \vdots\\ \hdashline \boldsymbol{e}_D^T \end{matrix}\right) \]

設$\mathbf{U} = \left(\begin{array}{c:c:c}\boldsymbol{e}_1& \cdots &\boldsymbol{e}_D\end{array}\right)$（此時$\mathbf{J} = \left(\begin{matrix} \boldsymbol{e}_1^T\\ \hdashline \vdots\\ \hdashline \boldsymbol{e}_D^T \end{matrix}\right) = \mathbf{U}^T$），由于新的基底$\boldsymbol{e}_1, \cdots, \boldsymbol{e}_D$是規范正交的，因此$\mathbf{U}$是正交矩陣，此時我們有：

\[|\mathrm{det}\mathbf{J}| = |\mathrm{det}\mathbf{U}^T| = \sqrt{\left(\mathrm{det}\mathbf{U}^T\right)^2} = \sqrt{\mathrm{det}\mathbf{U}^T\mathrm{det}\mathbf{U}} = \sqrt{\mathrm{det}\left(\mathbf{U}^T\mathbf{U}\right)} = \sqrt{\mathrm{det}\mathbf{I}} = 1 \]

而我們發現又多元高斯分布$\mathcal{N}\left(\boldsymbol{y}\mid \boldsymbol{0}, \boldsymbol{\Sigma}\right)$可以分解成$D$個獨立一元高斯分布$\mathcal{N}(y_i\mid 0, \lambda_i)$的乘積：

\[\begin{aligned} \mathcal{N}\left(\boldsymbol{y}\mid \boldsymbol{\mu}, \boldsymbol{\Sigma}\right) &= \frac{1}{(2\pi)^{\frac{D}{2}}}\frac{1}{(\mathrm{det}\boldsymbol{\Sigma})^{\frac{1}{2}}}\exp\left\{-\frac{1}{2}\sum_{i=1}^D \frac{y_i^2}{\lambda_i}\right\}\\ &= \prod_{i=1}^D\frac{1}{(2\pi\lambda_i)^{\frac{1}{2}}}\exp\left\{-\frac{y_i^2}{2\lambda_i}\right\}\\ &= \prod_{i=1}^D\mathcal{N}(y_i\mid 0, \lambda_i) \end{aligned} \]

（其中我們用到了結論$\mathrm{det}\mathbf{\Sigma} = \prod_i\lambda_i$）因此

\[\int \mathcal{N}(\boldsymbol{y}\mid \boldsymbol{0}, \boldsymbol{\Sigma})\mathrmw0obha2h00\boldsymbol{y} = \prod_{i=1}^D \int \mathcal{N}(y_i\mid 0, \lambda_i)\mathrmw0obha2h00y_i = 1 \]

于是我們有

\[\int \mathcal{N}\left(\boldsymbol{x}\mid \boldsymbol{\mu}, \boldsymbol{\Sigma}\right)\mathrmw0obha2h00\boldsymbol{x} = \int \mathcal{N}\left(T(\boldsymbol{x})\mid \boldsymbol{0}, \boldsymbol{\Sigma}\right)\mathrmw0obha2h00\boldsymbol{x} = \int \mathcal{N}(\boldsymbol{y}\mid \boldsymbol{0}, \boldsymbol{\Sigma})\mathrmw0obha2h00\boldsymbol{y} = 1 \]

至此我們證明了多元高斯分布是歸一化的。

2.2 高斯分布的一階矩和二階矩

現在我們考察多元高斯分布的一階矩和二階矩，這可以提供參數$\boldsymbol{\mu}$和$\mathbf{\Sigma}$的描述。多元高斯分布下$\boldsymbol{x}$的期望為：

\[\begin{aligned} \mathbb{E}[\boldsymbol{x}] &= \frac{1}{(2\pi)^{\frac{D}{2}}}\frac{1}{(\mathrm{det}\boldsymbol{\Sigma})^{\frac{1}{2}}} \int \exp\left\{-\frac{1}{2}(\boldsymbol{x} - \boldsymbol{\mu})^T\boldsymbol{\Sigma}^{-1}(\boldsymbol{x} - \boldsymbol{\mu})\right\} \boldsymbol{x} \mathrmw0obha2h00\boldsymbol{x}\\ &= \frac{1}{(2\pi)^{\frac{D}{2}}}\frac{1}{(\mathrm{det}\boldsymbol{\Sigma})^{\frac{1}{2}}} \int \exp\left\{-\frac{1}{2}\boldsymbol{z}^T\boldsymbol{\Sigma}^{-1}\boldsymbol{z}\right\} (\boldsymbol{z} + \boldsymbol{\mu}) \mathrmw0obha2h00\boldsymbol{z}\\ &= \frac{1}{(2\pi)^{\frac{D}{2}}}\frac{1}{(\mathrm{det}\boldsymbol{\Sigma})^{\frac{1}{2}}} \int \exp\left\{-\frac{1}{2}\boldsymbol{z}^T\boldsymbol{\Sigma}^{-1}\boldsymbol{z}\right\} \boldsymbol{z} \mathrmw0obha2h00\boldsymbol{z} + \boldsymbol{\mu} \end{aligned} \]

其中在第二個等式中我們使用了$\boldsymbol{z} = \boldsymbol{x} - \boldsymbol{\mu}$進行了變量替換，在第三個等式中我們利用了$\int \mathcal{N}(\boldsymbol{z}\mid \boldsymbol{0}, \mathbf{\Sigma})\mathrmw0obha2h00\boldsymbol{z} = 1$。現在我們來考慮第一個積分項，我們注意到指數位置是$\boldsymbol{z}$的偶函數，而$\boldsymbol{z}$是奇函數，因此被積函數是奇函數^[10]，又由于積分區域關于原點對稱，因此第一個積分項為0。于是我們有

\[\mathbb{E}[\boldsymbol{x}] = \boldsymbol{\mu} \]

因此我們把$\boldsymbol{\mu}$稱為多元高斯分布的均值。

現在我們考慮多元高斯分布的二階矩。在一元變量的情形下，二階矩由$\mathbb{E}[x^2]$給出。對于多元高斯分布，有$D^2$個由$\mathbb{E}[x_ix_j]$給出的二階矩，可以聚集在一起組成矩陣$\mathbb{E}[\boldsymbol{x}\boldsymbol{x}^T]$。這個矩陣可以表示為：

\[ \begin{aligned} \mathbb{E}[\boldsymbol{x}\boldsymbol{x}^T] &= \frac{1}{(2\pi)^{\frac{D}{2}}}\frac{1}{(\mathrm{det}\boldsymbol{\Sigma})^{\frac{1}{2}}} \int \exp\left\{-\frac{1}{2}(\boldsymbol{x} - \boldsymbol{\mu})^T\boldsymbol{\Sigma}^{-1}(\boldsymbol{x} - \boldsymbol{\mu})\right\} \boldsymbol{x}\boldsymbol{x}^T \mathrmw0obha2h00\boldsymbol{x}\\ &= \frac{1}{(2\pi)^{\frac{D}{2}}}\frac{1}{(\mathrm{det}\boldsymbol{\Sigma})^{\frac{1}{2}}} \int \exp\left\{-\frac{1}{2}\boldsymbol{z}^T\boldsymbol{\Sigma}^{-1}\boldsymbol{z}\right\} (\boldsymbol{z} + \boldsymbol{\mu})(\boldsymbol{z} + \boldsymbol{\mu})^T \mathrmw0obha2h00\boldsymbol{z}\\ \end{aligned} \]

其中在第二個等式中我們再次使用了$\boldsymbol{z} = \boldsymbol{x} - \boldsymbol{\mu}$來進行變量替換，涉及到$\boldsymbol{z}\boldsymbol{\mu}^T$和$\boldsymbol{\mu}\boldsymbol{z}^T$的交叉項將再次變為0，而$\boldsymbol{\mu}\boldsymbol{\mu}^T$也可以拿出。于是我們有

\[\mathbb{E}[\boldsymbol{x}\boldsymbol{x}^T] = \frac{1}{(2\pi)^{\frac{D}{2}}}\frac{1}{(\mathrm{det}\boldsymbol{\Sigma})^{\frac{1}{2}}} \int \exp\left\{-\frac{1}{2}\boldsymbol{z}^T\boldsymbol{\Sigma}^{-1}\boldsymbol{z}\right\} \boldsymbol{z}\boldsymbol{z}^T \mathrmw0obha2h00\boldsymbol{z} + \boldsymbol{\mu}\boldsymbol{\mu}^T \]

接下來考慮第一個積分項，我們采用和之前類似的做法，對$\mathbf{\Sigma}^{-1}$進行譜分解得$\mathbf{\Sigma}^{-1} = \sum_{i} \frac{1}{\lambda_i}\boldsymbol{e}_i\boldsymbol{e}_i^T$，并使用變量替換令$y_i = \boldsymbol{e}_i^T\boldsymbol{z}$（基變換關系為$\left(\begin{matrix} y_1\\ \vdots\\ y_D \end{matrix}\right) = \left(\begin{matrix} \boldsymbol{e}_1^T\\ \hdashline \vdots\\ \hdashline \boldsymbol{e}_D^T \end{matrix}\right)\boldsymbol{z}$，因此$\boldsymbol{z} = \left(\begin{array}{c:c:c}\boldsymbol{e}_1& \cdots &\boldsymbol{e}_D\end{array}\right)\boldsymbol{y} = \sum_iy_i\boldsymbol{e}_i$），于是有

\[\begin{aligned} &\frac{1}{(2\pi)^{\frac{D}{2}}}\frac{1}{(\mathrm{det}\boldsymbol{\Sigma})^{\frac{1}{2}}} \int \exp\left\{-\frac{1}{2}\boldsymbol{z}^T\boldsymbol{\Sigma}^{-1}\boldsymbol{z}\right\} \boldsymbol{z}\boldsymbol{z}^T \mathrmw0obha2h00\boldsymbol{z}\\ &= \frac{1}{(2\pi)^{\frac{D}{2}}}\frac{1}{(\mathrm{det}\boldsymbol{\Sigma})^{\frac{1}{2}}} \sum_{i=1}^D\sum_{j=1}^D\boldsymbol{e}_i\boldsymbol{e}_j^T\int \exp\left\{-\frac{1}{2}\sum_{k=1}^D \frac{y_k^2}{\lambda_k}\right\} y_iy_j \mathrmw0obha2h00\boldsymbol{y}\\ &= \frac{1}{(2\pi)^{\frac{D}{2}}}\frac{1}{(\mathrm{det}\boldsymbol{\Sigma})^{\frac{1}{2}}} \sum_{i=1}^D\boldsymbol{e}_i\boldsymbol{e}_i^T\int \exp\left\{-\frac{1}{2}\sum_{k=1}^D \frac{y_k^2}{\lambda_k}\right\} y_i^2 \mathrmw0obha2h00\boldsymbol{y}\\ &= \sum_{i=1}^D\boldsymbol{e}_i\boldsymbol{e}_i^T\lambda_i\\ &= \mathbf{\Sigma} \end{aligned} \]

其中第二個等式是由于當$i\neq j$時，被積函數是奇函數，導致積分項為0；第三個等式是由于$\mathbb{E}[y_i^2] = \int \mathcal{N}(y_i\mid 0, \lambda_i)y_i^2\mathrmw0obha2h00y_i=\lambda_i + 0^2 = \lambda_i$；最后一個等式是根據$\mathbf{\Sigma}$的譜分解得到。

這樣，我們就得到了

\[\mathbb{E}[\boldsymbol{x}\boldsymbol{x}^T] = \mathbf{\Sigma} + \boldsymbol{\mu}\boldsymbol{\mu}^T \]

對于一元隨機變量的方差，為了定義方差，我們在取二階矩之前會減掉均值。類似地，對于多元變量的情形，把均值減掉同樣很方便。這給出了隨機變量$\boldsymbol{x}$的協方差（covariance） 的定義：

\[\mathrm{Cov}[\boldsymbol{x}] = \mathbb{E}\left[(\boldsymbol{x} - \mathbb{E}[\boldsymbol{x}])(\boldsymbol{x} - \mathbb{E}[\boldsymbol{x}])^T\right] \]

對于多元高斯分布這一特例，我們有

\[\begin{aligned} \mathrm{Cov}[\boldsymbol{x}] &= \mathbb{E}\left[\boldsymbol{x}\boldsymbol{x}^T - \boldsymbol{x}\boldsymbol{\mu}^T - \boldsymbol{\mu}\boldsymbol{x}^T + \boldsymbol{\mu}\boldsymbol{\mu}^T\right]\\ &= \mathbb{E}[\boldsymbol{x}\boldsymbol{x}^T] - \boldsymbol{\mu}\boldsymbol{\mu}^T\\ &= \mathbf{\Sigma} \end{aligned} \]

由于參數$\mathbf{\Sigma}$控制了多元高斯分布下$\boldsymbol{x}$的協方差，因此它被稱為協方差矩陣。

雖然式$(1)$定義的高斯分布$\mathcal{N}(\boldsymbol{x}\mid \boldsymbol{\mu}, \boldsymbol{\Sigma})$被廣泛用作概率密度模型，但是它有著一些巨大的局限性。考慮分布中自由參數的數量。一個通常的對稱協方差矩陣$\boldsymbol{\Sigma}$有$1 + 2 + \cdots + D = \frac{D(D + 1)}{2}$個獨立參數，$\boldsymbol{\mu}$中有另外$D$個獨立參數，因此總計有$\frac{D(D + 1)}{2} + D = \frac{D(D + 3)}{2}$個獨立參數。對于大的$D$值，參數的總數隨著$D$以平方的方式增長，導致對大的矩陣進行操作（如求逆）的計算變得不可行。解決這個問題的一種方式是使用協方差矩陣的限制形式。如果我們考慮對角的協方差矩陣，即$\mathbf{\Sigma} = \mathrm{diag}(\sigma_i^2)$，那么在概率密度模型中，我們就有總數$2D$個獨立參數。此時常數概率密度輪廓線為與軸對齊的橢球。我們可以進一步地把協方差矩陣限制成正比于單位矩陣，也即$\mathbf{\Sigma} = \sigma^2 \mathbf{I}$，此時它被稱為各向同性（isotropic）的協方差。這使得模型有$D + 1$個獨立的參數，并且常數概率密度輪廓線為球面。下圖展示了通常的協方差矩陣、對角的協方差矩陣以及各向同性協方差矩陣的概率密度輪廓線（$D=2$）：

高斯分布的另一個局限性是它本質上是單峰的（即只有一個最大值），因此不能很好地近似多峰分布。不過，相當多的多峰分布可以使用混合高斯分布來描述（參見博客《統計學習：EM算法及其在高斯混合模型(GMM)中的應用》）。

參考

[1] Bishop C M, Nasrabadi N M. Pattern recognition and machine learning[M]. New York: springer, 2006.
[2] Schroeder D V. An introduction to thermal physics[M]. Oxford University Press, 2020.
[3] 《維基百科：熵 (統計物理學)》
[4] Feller W. An introduction to probability theory and its applications, Volume 1[M]. John Wiley & Sons, 1991.
[5] Weisstein E W. CRC concise encyclopedia of mathematics[M]. Chapman and Hall/CRC, 2002.
[6] Bengio Y, Goodfellow I, Courville A. Deep learning[M]. Cambridge, MA, USA: MIT press, 2017.
[7] Johnson R A, Wichern D W. Applied multivariate statistical analysis[J]. 2002.
[8] Rudin W. Principles of mathematical analysis[M]. New York: McGraw-hill, 1964.
[9] Axler S. Linear algebra done right[M]. springer publication, 2015.
[10] 維基百科：《奇函數與偶函數》)

posted @ 2025-01-23 23:12 orion-orion 閱讀(882) 評論(0) 收藏舉報

刷新頁面返回頂部

Orion's Blog

學習理論、在線優化、數據挖掘