機器學習 —— 基礎整理（四）特徵提取之線性方法：主成分分析PCA、獨立成分分析ICA、線性判別分析LDA

     本文簡單整理了如下內容：

（一）維數災難

（二）特徵提取——線性方法

1. 主成分分析PCA

2. 獨立成分分析ICA

3. 線性判別分析LDA

 

（一）維數災難（Curse of dimensionality）

      維數災難就是說當樣本的維數增長時，若要保持與低維情形下相同的樣本密度，所須要的樣本數指數型增加。從下面的圖能夠直觀體會一下。當維度很大樣本數量少時，沒法經過它們學習到有價值的知識；因此須要降維，一方面在損失的信息量能夠接受的狀況下得到數據的低維表示，增長樣本的密度；另外一方面也能夠達到去噪的目的。

圖片來源：[1]

      有兩種辦法試圖克服這個問題：首先是特徵提取（Feature extraction），經過組合現有特徵來達到降維的目的；而後是特徵選擇，從現有的特徵裏選擇較小的一些來達到降維的目的。下面開始介紹特徵提取。

（二）特徵提取——線性方法

      首先約定一下記號：

      樣本矩陣 $X=(\textbf x_1,\textbf x_2,...,\textbf x_N)=({\textbf x^{(1)}}^{\top};{\textbf x^{(2)}}^{\top};...;{\textbf x^{(d)}}^{\top})\in\mathbb R^{d\times N}$ ，$X_{ji}=x_i^{(j)}$；

      $N$ 是訓練集的樣本數，每一個樣本都表示成 $\textbf x_i=(x_i^{(1)},x_i^{(2)},...,x_i^{(d)})\in\mathbb R^d$ 的列向量，樣本矩陣的每一列都是一個樣本；若是出現了 $\textbf x$ 這樣沒有上下標的記號就泛指任同樣本，至關於省略下標；

      $d$ 是特徵的維數，每維特徵都表示成 $\textbf x^{(j)}=(x_1^{(j)},x_2^{(j)},...,x_N^{(j)})\in\mathbb R^N$ 的列向量；若是出現了 $x^{(j)}$ 這樣的記號就泛指任同樣本的第 $j$ 維特徵，至關於省略下標；

      降維以後的樣本矩陣 $Z\in\mathbb R^{k\times N}$ ，$k$ 是降維以後的特徵維數，因此 $k\ll d$ ；每一個樣本被降維後都表示成 $\textbf z_i=(z_i^{(1)},z_i^{(2)},...,z_i^{(d)})\in\mathbb R^k$ 的列向量。

1、主成分分析（Principal component analysis，PCA）

      PCA也叫主元分析，可謂是很是經常使用的一種線性降維方式，好比在人臉識別中的「特徵臉」（Eigenfaces）。降維以後的每一個「新」特徵都被稱爲主成分。這是一種無監督的降維方法，沒有用到樣本的標記信息。

      1. 定義：線性投影

      下圖是一個二維的例子，圖示的兩個方向爲兩個主成分的基底，咱們能夠經過PCA將它降成一維，也就是隻保留一個基底。

圖片來源：[6]

      將原始數據 $X\in\mathbb R^{d\times N}$ 經過變換矩陣 $W=(\textbf w_1,\textbf w_2,...,\textbf w_k)\in\mathbb R^{d\times k}$ 投影到低維空間，獲得變換後的數據 $Z\in\mathbb R^{k\times N}$ 

$$\begin{aligned}Z&=W^{\top}X\\&=(\textbf z_1,\textbf z_2,...,\textbf z_N)\\&=({\textbf z^{(1)}}^{\top};{\textbf z^{(2)}}^{\top};...;{\textbf z^{(d)}}^{\top})\end{aligned}$$

樣本 $\textbf x_i$的通過映射後的獲得的「新」樣本 $\textbf z_i$ 爲

$$\textbf z_i=W^{\top}\textbf x_i$$

份量表示爲

$$z_i^{(j)}={\textbf w_j}^{\top}\textbf x_i,\quad j=1,2,...,d$$

能夠直觀看出，每一個新的特徵都是原先所有特徵的線性組合。

      「新」特徵 $\textbf z^{(1)}={\textbf w_1}^{\top}X$ 稱爲第一主成分，隨後是第二主成分 $\textbf z^{(2)}$，第三主成分……，只保留 $k$ 個主成分，構成樣本 $\textbf x_i$ 降維後的表示 $\textbf z_i$ 。

      2. 問題描述：矩陣分解

      將變換矩陣記成 $W=(\textbf w_1,\textbf w_2,...,\textbf w_k)$ ，這是低維空間的一組標準正交基，每一列 $\textbf w_i\in\mathbb R^d$ 與其餘列之間是互相正交的，而列自己的模爲1。換句話說，矩陣 $W$ 須要知足正交約束（不能叫正交矩陣，由於只有方陣才談得上是否是正交矩陣）：

$$W^{\top}W=I$$

      這樣看的話，PCA 其實就是正交約束下的矩陣分解問題：

$$X=WZ$$

$$\textbf x_i=W\textbf z_i$$

      3. 求解方法：協方差矩陣的特徵值分解

      如今的問題是，怎樣肯定出變換矩陣 $W$ 。PCA的基本作法是，在原始數據 $X$ （要去均值，也就是中心化）的協方差矩陣 $\varSigma_X$ 作特徵值分解（實對稱矩陣必定能找到一個正交矩陣使其對角化，這個正交矩陣正是由其特徵值對應的特徵向量所組成），將特徵值從大到小排序，其中最大的特徵值 $\lambda_1$ 所對應的特徵向量 $\textbf w_1$ 做用在樣本 $X$ 所獲得的「新」特徵 $\textbf z^{(1)}={\textbf w_1}^{\top}X$ 稱爲第一主成分，隨後是第二主成分 $\textbf z^{(2)}$，第三主成分……，只保留 $k$ 個主成分，構成樣本 $\textbf x_i$ 降維後的表示 $\textbf z_i$ 。$k$ 的值可經過下式來肯定：

$$\frac{\sum\limits_{i=1}^k\lambda_i}{\sum\limits_{i=1}^d\lambda_i}\geq 95\%$$

      4. PCA的motivation

      爲何 PCA 要對協方差矩陣 $\varSigma_X$ 作特徵值分解來獲得主成分？下面簡單說兩種角度。

      1) 最大投影方差

      假設訓練集樣本的每一維特徵都進行了去均值處理，使得 $\boldsymbol\mu (X)=\frac1N\sum_{i=1}^N\textbf x_i=\textbf 0\in\mathbb R^d$（須要注意，測試集的樣本在作去均值操做時去掉的是訓練集的均值，因此這個操做中獲得的訓練集的均值向量要保留）。

      直接求解

      考慮到降維到 $k$ 維的狀況：但願投影后的每一維都儘量地分開。

      因此咱們求取投影后的樣本的協方差矩陣：

$$\begin{aligned}\varSigma_Z&=\frac{1}{N-1}\sum_{i=1}^N(\textbf z_i - \boldsymbol\mu(Z))(\textbf z_i - \boldsymbol\mu(Z))^{\top}\\&=\frac{1}{N-1}\sum_{i=1}^N\textbf z_i\textbf z_i^{\top}\\&=\frac{1}{N-1}ZZ^{\top}\\&=\frac{1}{N-1}W^{\top}XX^{\top}W\end{aligned}$$

由於

$$\varSigma_X=\frac{1}{N-1}\sum_{i=1}^N\textbf x_i\textbf x_i^{\top}=\frac{1}{N-1}XX^{\top}$$

不難看出 $\varSigma_X$ 就是樣本去均值以後的協方差矩陣。因此有

$$\varSigma_Z=W^{\top}\varSigma_X W$$

      這個形式很眼熟：首先，這和高斯分佈經線性變換後的協方差矩陣形式同樣，即原始協方差矩陣的二次型；其次，若是不降維的話（也就是說 $W$ 是個正交方陣），單從代數的角度來看，這裏咱們須要求取的 $\varSigma_Z$ ，不正是把協方差矩陣 $\varSigma_X$ 對角化嗎？因此 $W$ 已經很顯然了，就是 $\varSigma_X$ 的特徵向量構成的矩陣。

      因爲 $\varSigma_Z$ 的對角線元素就表明了各維的方差，因此PCA的優化目標就是

$$\max_W\quad\text{tr}(W^{\top}\varSigma_X W)$$

約束是 $W^{\top}W=I$ ，而後用拉格朗日乘數法求解。但我以爲這樣不直觀，接下來從降維到一維入手，逐步求解主成分。

      逐步求解

      不妨先實現一個小目標：咱們須要降維到一維，樣本點在投影到一維以後，可以儘量地「分開」——也就是投影后的方差最大。這時變換矩陣 $W$ 退化爲一個向量 $\textbf w_1$ ，樣本 $\textbf x_i$ 投影的結果爲 $\textbf z_i =\textbf w_1^{\top}\textbf x_i\in\mathbb R$。

      投影后的方差爲：

$$\begin{aligned}&\frac{1}{N-1}\sum_{i=1}^N\textbf z_i\textbf z_i^{\top}\\=&\frac{1}{N-1}\sum_{i=1}^N(\textbf w_1^{\top}\textbf x_i)(\textbf w_1^{\top}\textbf x_i)^{\top}\\=&\textbf w_1^{\top}(\frac{1}{N-1}\sum_{i=1}^N\textbf x_i\textbf x_i^{\top})\textbf w_1\\=&\textbf w_1^{\top}\varSigma_X\textbf w_1\end{aligned}$$

      爲了實現上面降維到一維後方差最大的這個小目標，咱們須要求解以下的約束最優化問題：

$$\max_{\textbf w_1}\quad\textbf w_1^{\top}\varSigma_X\textbf w_1$$

$$\text{s.t.}\quad\quad\|\textbf w_1\|=1$$

使用拉格朗日乘數法，那麼問題的解就是拉格朗日函數 $L$ 的偏導數等於零這個方程的解：

$$L=\textbf w_1^{\top}\varSigma_X\textbf w_1-\lambda(\textbf w_1^{\top}\textbf w_1-1)$$

$$\frac{\partial L}{\partial \textbf w_1}=2\varSigma_X\textbf w_1-2\lambda\textbf w_1=0$$

$$\varSigma_X\textbf w_1=\lambda\textbf w_1$$

因此問題的解須要是 $\varSigma_X$ 的特徵向量。優化目標 $\textbf w_1^{\top}\varSigma_X\textbf w_1=\textbf w_1^{\top}\lambda\textbf w_1=\lambda$ ，若要使它最大化，便要使 $\lambda$ 最大化。因此問題的解是 $\varSigma_X$ 的最大特徵值所對應的特徵向量。

      如今已經證實了，第一主成分怎樣得來。那麼考慮到 $k$ 個主成分的狀況：如何證實，最大的 $k$ 個特徵值對應的特徵向量所組成的矩陣 $W$ ，知足投影后的各維方差都儘量大？

      固然是數學概括法：如今 $k=1$ 時成立（概括基礎），咱們假設 $k=m$ 時成立，只要論證出 $k=m+1$ 時仍成立，那結論就是成立的。

      如今已知 $\textbf w_1,\textbf w_2,...,\textbf w_m$ 是一組在新空間的維度爲 $m$ 下知足投影方差最大的基底。$\textbf w_{m+1}$ 要知足的條件有：

      (1) 模爲1，$\|\textbf w_{m+1}\|=1$；(2) 與 $\textbf w_1,\textbf w_2,...,\textbf w_m$ 都正交，$\textbf w_{m+1}^{\top}\textbf w_j=0$ ；(3) $\textbf w_{m+1}^{\top}X$ 的方差最大

      寫成約束最優化的形式，就是

$$\max_W\quad \textbf w_{m+1}^{\top}\varSigma_X\textbf w_{m+1}$$

$$ \begin{aligned} \text{s.t.}\quad & \|\textbf w_{m+1}\|=1 \\& \textbf w_{m+1}^{\top}\textbf w_j=0,\quad j=1,2,...,m \end{aligned} $$

使用拉格朗日乘數，可得

$$L=\textbf w_{m+1}^{\top}\varSigma_X\textbf w_{m+1}-\lambda (\textbf w_{m+1}^{\top}\textbf w_{m+1}-1)-\sum_{j=1}^m\eta_j\textbf w_{m+1}^{\top}\textbf w_j$$

$$\frac{\partial L}{\partial \textbf w_1}=2\varSigma_X\textbf w_{m+1}-2\lambda\textbf w_{m+1}-\sum_{j=1}^m\eta_j\textbf w_j=0$$

將等式兩邊依次右乘 $\textbf w_j, j=1,2,...,m$ ，就能夠依次獲得 $\eta_j=0, j=1,2,...,m$，因此有

$$\varSigma_X\textbf w_{m+1}=\lambda\textbf w_{m+1}$$

這樣就論證了 $\textbf w_{m+1}$ 就是協方差矩陣的第 $m+1$ 大的特徵值所對應的特徵向量。

      2) 最小均方重建偏差（mean square reconstruction error，MSRE）

      咱們知道，若是要從 $Z$ 再重建到原先數據所在的空間中，須要作的變換是 $WZ$ 。最小均方重建偏差的優化目標爲 $\min_W\quad\|X-WW^{\top}X\|^2$ ，約束爲 $W^{\top}W=I$ 。能夠推導出，優化目標等價於 $\max_W\quad\text{tr}(W^{\top}\varSigma_X W)$ 。也就是說和最大投影方差的優化目標一致。

      3) 除此以外...

      PCA還有另外若干種理解角度，如高斯隨機採樣。從這種角度理解能夠參考[2]，此外還介紹了PCA在指數族分佈上的推廣。

      5. PCA 和 SVD 的關係

      參考知乎的這個回答。

      關於SVD，曾經在介紹 LSA 模型時提到過：

$$X=U\varSigma V^{\top}$$

那麼就能夠獲得

$$XX^{\top}=U\varSigma \varSigma^{\top}U^{\top}$$

      因此對比一下

$$\varSigma_X=W\varSigma_Z W^{\top}$$

就知道 U 和 W 的地位是同樣的。

      這反映出了奇異值與特徵值的關係：矩陣 $X$ 的奇異值對應於 $XX^{\top}$ 的特徵值的平方根，$X$ 的左奇異向量對應於 $XX^{\top}$ 的特徵向量。

      6. 其餘說明

      對於測試集，中心化操做中所減去的均值是訓練集的；降維操做中的矩陣使用的是訓練集的 $W$ 。機器學習方法中全部的預處理、調超參數等操做都不能混入任何測試集的信息，不妨想象：測試集就一個樣本。

      變換矩陣 $W$ 能夠經過一個單隱層的AutoEncoder來求解，隱層神經元個數就是所謂的 $k$ ，這樣就能夠用神經網絡那一套來訓練了。這招是從 [1] 看的，恍然大悟的感受，閱讀經典的專著總會有意想不到的收穫。

      在實際的工業生產環境下，每每要面對數據量很大的狀況，須要在線計算協方差矩陣，使得主成分隨着新數據的到來而實時更新。方法是Oja’s Rule。

      降維效果好壞的比較，應該仍是從具體任務上的分類效果孰優孰劣來比較的。從低維的圖裏能夠看出來，PCA仍是比較適合於樣本服從高斯分佈，因此有些時候PCA降維後未必有很好的效果。

      下面這個圖是我本科畢設的實驗，虛線的是PCA以前的，可見在這個特徵空間下PCA降維對每一個類的效果都有提高。

2、獨立成分分析（Independent component analysis，ICA）

      ICA相比於PCA，其追求的效果是不同的：ICA尋找的是最能使數據的相互獨立的方向，而PCA僅要求方向是不相關的。咱們知道，獨立能夠推出不相關，反之則不能夠，而高斯分佈的狀況下獨立等價於不相關。所以ICA須要數據的高階統計量，PCA則只須要二階統計量。

 

圖片來源：[6]

      1. 背景

      考慮盲信號分離（Blind signal separation）的問題：設有 $d$ 個獨立的標量信號源發出聲音，其在時刻 $t$ 發出的聲音可表示爲 $\textbf s_t=(s_t^{(1)},s_t^{(2)},...,s_t^{(d)})^{\top}\in\mathbb R^d$ 。一樣地，有 $d$ 個觀測器在進行採樣，其在時刻 $t$ 記錄的信號可表示爲：$\textbf x_t\in\mathbb R^d$ 。認爲兩者知足下式，其中矩陣 $A\in\mathbb R^{d\times d}$ 被稱爲mixing matrix，反映信道衰減參數：

$$\textbf x_t=A\textbf s_t$$

      顯然，有多少個採樣時刻，就能夠理解爲有多少個樣本；而信號源的個數能夠理解爲特徵的維數。ICA的目標就是從 $\textbf x$ 中提取出 $d$ 個獨立成分，也就是找到矩陣unmixing matrix $W$ ：

$$\textbf s_t=W\textbf x_t,\quad W=A^{-1}$$ 

圖片來源：[1]

      將矩陣 $W$ 記爲 $W=({\textbf w_1}^{\top};{\textbf w_2}^{\top};...;{\textbf w_d}^{\top})$ ，也就是它的第 $j$ 行是 ${\textbf w_j}^{\top}$ ，那麼 $s_i^{(j)}={\textbf w_j}^{\top}\textbf x_i$ 。（這裏的 $W$ 相比於PCA推導中的 $W$ 差一個轉置）

      2. 求解

      在沒有其餘先驗知識的狀況下，因爲上面式子中的只有觀測信號是已知的，故沒法求解：

      首先是源信號的幅值不肯定，致使矩陣 $W$ 隨幅值的變化而變化；其次是源信號的順序不肯定，即便是徹底相同的兩組源信號而僅僅是順序不一樣，由於當順序更改了以後致使矩陣 $W$ 的行排列隨之變化，因此沒法肯定出惟一的矩陣 $W$ 。

      再有，源信號不能是高斯分佈的。考慮2維狀況，$\textbf s\sim N(\textbf 0,I)$ 。因爲觀測信號 $\textbf x_t=A\textbf s_t$ ，根據高斯分佈的線性不變性，$\textbf x\sim N(\textbf 0,AIA^{\top})=N(\textbf 0,AA^{\top})$ 。現考慮另一個mixing matrix $A'=AR$ ，其中 $R$ 是正交矩陣。則這時的觀測信號 $\textbf x_t'=A'\textbf s_t$ ，且$\textbf x'\sim N(\textbf 0,A'IA'^{\top})=N(\textbf 0,ARR^{\top}A^{\top})=N(\textbf 0,AA^{\top})$ ，不一樣的mixing matrix卻獲得了相同的觀測信號。可見當源信號服從高斯分佈時 $W$ 不能惟一肯定。從圖示能夠看出，高斯分佈的線性不變性致使ICA失效。

 

圖片來源：[1]，圖中的下標對應於本文的上標

      如今考慮ICA的求解。以前說過，$d$ 個源信號是相互獨立的（且沒有噪聲），因此源信號的密度函數能夠表示爲

$$p_{\textbf s}(\textbf s)=\prod_{j=1}^dp_s(s^{(j)})$$

      觀測信號 $\textbf s_t$ 和觀測信號 $\textbf x_t$ 的關係是$\textbf x_t=A\textbf s_t$，它們的機率密度函數有以下關係： 

$$p_{\textbf x}(\textbf x)=\frac{p_{\textbf s}(\textbf s)}{|A|}=p_{\textbf s}(\textbf s)|W|$$

能夠獲得下式

$$p_{\textbf x}(\textbf x)=|W|\prod_{j=1}^dp_s({\textbf w_j}^{\top}\textbf x)$$

      如今須要作的是指定 $p_s(\cdot)$ 。在沒有任何先驗知識的狀況下，能夠指定

$$\begin{aligned}p_s(\cdot)&=\sigma '(\cdot)\\&=(\dfrac{1}{1+\exp(-\cdot)})'\\&=\dfrac{\exp(-\cdot)}{(1+\exp(-\cdot))^2}\\&=\sigma(\cdot)(1-\sigma(\cdot))\end{aligned}$$

這至關於分佈函數是logistic函數，密度函數天然就是它的導函數。這是一個合理的指定，由於在不少問題上都work well；若是有先驗知識固然能夠指定爲其餘的形式。

      給定 $N$ 個時刻的觀測值，使用極大似然估計，獲得似然函數爲 $L(W)=\prod_{i=1}^Np_{\textbf x}(\textbf x_i)$，進一步獲得對數似然函數爲

$$l(W)=\sum_{i=1}^N(\sum_{j=1}^d\log \sigma'({\textbf w_j}^{\top}\textbf x_i)+\log |W|)$$

      爲了求極大值，須要用梯度上升法。根據公式 $\dfrac{\partial |W|}{\partial W}=|W|(W^{-1})^{\top}$ 以及 $\log\sigma'(\cdot)=1-2\sigma(\cdot)$，可求得對於一個樣本的梯度爲

$$\frac{\partial l_{1}}{\partial W}=(\textbf 1-2\sigma({\textbf w_j}^{\top}\textbf x_i))\textbf x_i^{\top}+(W^{-1})^{\top}$$

這裏若是logistic函數的自變量爲矩陣，就對逐個元素作運算，函數值是階數相同的矩陣。這樣就獲得了優化目標對於一個樣本的梯度，進而能夠用梯度上升法迭代更新最大值。

      附上sklearn中對PCA和ICA的介紹：2.5. Decomposing signals in components (matrix factorization problems)

 

3、線性判別分析（Linear Discriminant Analysis，LDA）

      LDA也稱爲Fisher判別分析，是從更利於分類的角度來降維，利用到了訓練樣本的類別標記，追求的是最可以分開各個類別數據的投影方法。從下圖能夠看出，相比於直接用PCA降維，若是將樣本投影到下圖所示的直線上則會更利於分類。

      LDA在NLP裏一般指的是 latent Dirichlet allocation（隱狄利克雷分配），是一種主題模型的簡稱。不要搞混了。

      嚴格說的話，LDA還要求了各個類別的協方差矩陣相等，這樣能夠從貝葉斯決策論的角度推導出LDA分類器的判別函數 $g_i(\textbf x)$ 爲線性函數（Fisher判別分析則沒有協方差矩陣相同的要求）。

 

圖片來源：[3]

      若想使得投影后的結果有利於分類，須要從兩方面的需求考慮：一方面，投影以後相同類別的樣本之間要儘量近；另外一方面，投影以後各個類別之間要儘量遠。

      1. 二分類的狀況

      首先考慮二類問題，將 $d$ 維樣本投影到一條方向由 $\textbf w$ 肯定的的直線上。如今的問題就是找出 $\textbf w$ 。

      一個容易想到的思路是投影后各個類的均值儘量遠。對於第 $i$ 類樣本，投影前的均值爲 $\boldsymbol\mu_i=\frac{1}{n_i}\sum_{\textbf x\in\omega_i}\textbf x$ ，投影后的均值爲 $\widetilde{\boldsymbol\mu}_i=\textbf w^{\top}\boldsymbol\mu_i$ ，而後使 $|\widetilde\mu_1-\widetilde\mu_2|^2$ 最大化（由於降維後是一維，因此均值向量就是個標量，用了細體字母）。但這樣實際上是不行的，從下圖一看便知：

圖片來源：[7]

      上述作法的癥結在於，同類樣本之間顯然太遠了，它只考慮了第二個需求而沒有考慮第一個需求。因此，須要用總類內散度來歸一化上述目標。首先定義以下幾個概念：

      類內散度矩陣（從定義中可看出與協方差矩陣相差一個常數倍）：

$$S_i=\sum_{\textbf x\in\omega_i}(\textbf x-\boldsymbol\mu_i)(\textbf x-\boldsymbol\mu_i)^{\top}$$

      總類內散度矩陣就是全部類別的類內散度矩陣之和：

$$S_W=\sum_{i=1}^cS_i$$

      類間散度矩陣（式中 $\boldsymbol\mu=\frac1n\sum_{\textbf x}\textbf x$ 爲所有數據的均值向量）：

      (1) 類別數等於2

$$S_B=(\boldsymbol\mu_1-\boldsymbol\mu_2)(\boldsymbol\mu_1-\boldsymbol\mu_2)^{\top}$$

      (2) 類別數大於2

$$S_B=\sum_{i=1}^cn_i(\boldsymbol\mu_i-\boldsymbol\mu)(\boldsymbol\mu_i-\boldsymbol\mu)^{\top}$$

      整體散度矩陣

$$S_T=S_W+S_B=\frac1n\sum_{\textbf x}(\textbf x-\boldsymbol\mu)(\textbf x-\boldsymbol\mu)^{\top}$$

      通過變換矩陣 $W$ 後，類內散度矩陣和類間散度矩陣變爲

$$\widetilde S_W=W^{\top}S_WW$$

$$\widetilde S_B=W^{\top}S_BW$$

（從形式上能夠看出，與高斯分佈通過線性變換後的新協方差矩陣的形式是同樣的）

      因此對於二類的狀況，優化目標能夠寫爲

$$J=\frac{|\widetilde\mu_1-\widetilde\mu_2|^2}{\widetilde s_1+\widetilde s_2}=\frac{\textbf w^{\top}S_B\textbf w}{\textbf w^{\top}S_W\textbf w}$$

該式是兩個散度矩陣 $S_B$ 、$S_W$ 的廣義瑞利商。

      1) 普通解法：特徵值分解

      直接令該式的偏導數等於零，能夠獲得 $S_W^{-1}S_B\textbf w=J\textbf w$ ，因此能夠看做是矩陣 $S_W^{-1}S_B$ 的特徵值分解問題。另外，因爲 $\textbf w$ 的幅值不影響問題的解，因此也能夠看做下述約束最優化問題：

$$\min_{\textbf w}\quad -\textbf w^{\top}S_B\textbf w$$

$$\text{s.t.}\quad \textbf w^{\top}S_W\textbf w=1$$

使用拉格朗日乘數法可求得該問題的解就是下式的解：

$$S_W^{-1}S_B\textbf w=\lambda\textbf w$$

      2) 另外一種解法：典範變量

      在這個問題中，其實不必作特徵值分解：

      根據 $S_B=(\boldsymbol\mu_1-\boldsymbol\mu_2)(\boldsymbol\mu_1-\boldsymbol\mu_2)^{\top}$ ，獲得 $S_B\textbf w=(\boldsymbol\mu_1-\boldsymbol\mu_2)(\boldsymbol\mu_1-\boldsymbol\mu_2)^{\top}\textbf w$ ，由於後兩個因子的乘積是標量（記做 $a$ ），因此 $S_B\textbf w$ 老是位於 $\boldsymbol\mu_1-\boldsymbol\mu_2$ 的方向上，進而有 $\lambda\textbf w=aS_W^{-1}(\boldsymbol\mu_1-\boldsymbol\mu_2)$ 。因此能夠馬上寫出問題的解爲

$$\textbf w=S_W^{-1}(\boldsymbol\mu_1-\boldsymbol\mu_2)$$

這個解有時被稱爲典範變量（Canonical variable）。

      2. 多分類的狀況

      對於多分類問題，設類別個數是 $c$ ，從分類的角度講是從 $d$ 維空間向 $c-1$ 維空間投影，從降維的角度來講的話是從 $d$ 維空間向 $k$ 維空間投影（ $k\leq c-1$ ）。記矩陣 $W=(\textbf w_1,...,\textbf w_{c-1})\in\mathbb R^{d\times(c-1)}$ 。優化目標一般是

$$J=\frac{|\widetilde S_B|}{|\widetilde S_W|}=\frac{\text{tr}(W^{\top}S_BW)}{\text{tr}(W^{\top}S_WW)}$$

其解的方程以下，最優矩陣 $W$ 的列向量是下列等式中最大特徵值對應的特徵向量：

$$S_W^{-1}S_B\textbf w_i=\lambda_i\textbf w_i$$

由於 $S_B$ 是 $c$ 個秩爲1或0的矩陣的和，其中只有 $c-1$ 個是獨立的，因此 $S_B$ 的秩不會超過 $c-1$ ，也就是說非零特徵值至多隻有 $c-1$ 個。

 

 

參考資料：

[1] 《模式分類》及slides

[2] 66天寫的邏輯迴歸

[3] Dimension Reduction: A Guided Tour

[4]［非線性方法推薦看這篇］【機器學習算法系列之三】簡述多種降維算法

[5] CS229 Lecture notes11：ICA

[6] Assessment of multidimensional functional neuroimaging data model by statistical resampling

[7] A Tutorial on Data Reduction Linear Discriminant Analysis

[8] Independent Component Analysis: A Tutorial