教科書上的LDA為什麼長這個樣子？

09-01

教科書上的LDA為什麼長這個樣子？

來自專欄 PaperWeekly13 人贊了文章

線性判別分析(Linear Discriminant Analysis, LDA)是一種有監督降維方法，有關機器學習的書上一定少不了對PCA和LDA這兩個演算法的介紹。LDA的標準建模形式是這樣的（這裡以兩類版本為例，文章會在幾個關鍵點上討論多類情況）：

$maxlimits_w quad frac{w^TS_Bw}{w^TS_Ww}$

其中， $S_B$ 是類間散布矩陣， $S_W$ 是類內散布矩陣， $w$ 是投影直線:

$S_B = (m_1 - m_2)(m_1 - m_2)^T, quad S_W = S_1 + S_2, quad w in R^d$

怎麼樣，一定非常熟悉吧，經典的LDA就是長這個樣子的。這個式子的目標也十分直觀：將兩類樣本投影到一條直線上，使得投影后的類間散布矩陣與類內散布矩陣的比值最大。

三個加粗的詞隱含著三個問題：

為什麼是類間散布矩陣 $(m_1-m_2)(m_1-m_2)^T$ 呢？直接均值之差 $m_1-m_2$ 不是更符合直覺嗎？這樣求出來的解和原來一樣嗎？
為什麼是投影到直線，而不是投影到超平面？PCA是把d維樣本投影到c維(c<d)，LDA為什麼不能也投影到c維，而是直接投影到1維呢？同樣地，在K類LDA中，為什麼書上寫的都是投影到K-1維，再高一點不行嗎？這是必然嗎？
為什麼是類間散布與類內散布的比值呢？差值不行嗎？

這篇文章就圍繞這三個問題展開。我們先回顧一下經典LDA的求解，然後順次講解分析這三個問題。

0. 回顧經典LDA：

原問題等價於這個形式：

$max limits_w quad w^TS_Bw \ qquad s.t. quad w^TS_Ww = 1$

然後就可以用拉格朗日乘子法了：

$L(w,lambda)=-w^TS_Bw+lambda(w^TS_Ww-1) \$

求導並令其為0：

$partial_w L(w,lambda) = -2S_Bw+2lambda S_Ww=0$

得到解：

$S_W^{-1}S_Bw=lambda w$

對矩陣 $S_W^{-1}S_B$ 進行特徵值分解就可以得到w。但是有更簡單的解法：

$S_W^{-1}(m_1-m_2)(m_1-m_2)^Tw=lambda w$

而其中 $(m_1-m_2)^Tw$ 是一個標量，所以 $S_W^{-1}(m_1-m_2)$ 和 $lambda w$ 共線，得到：

$w propto S_W^{-1}(m_1-m_2)$

求解完畢。

非常優雅，不愧是教科書級別的經典演算法，整個求解一氣呵成，標準的拉格朗日乘子法。但是求解中還是用到了一個小技巧： $(m_1-m_2)^Tw$ 是標量，從而可以免去特徵值分解的麻煩。那麼，我們能不能再貪心一點，找到一種連這個小技巧都不需要的求解方法呢？答案是可以，上面的問題(1)就是在做這個事情。

1. 類間散布&均值之差

我們不用類內散布矩陣了 $(m_1-m_2)(m_1-m_2)^T$ ，改用均值之差 $m_1-m_2$ 這個更符合直覺的東西：

$quad max limits_w quad w^T(m_1-m_2) \ s.t. quad w^TS_Ww = 1$

還是用拉格朗日乘子法：

$L(w,lambda)=-w^T(m_1-m_2)+lambda(w^TS_Ww-1)$

求導並令其為0：

$partial_w L(w,lambda) = -(m_1-m_2)+2lambda S_Ww=0$

得到解：

$w propto S_W^{-1}(m_1-m_2)$

怎麼樣，是不是求解更簡單了呢，不需要任何技巧，一步一步下來就好了。

為什麼說均值之差更符合直覺呢？大家想啊，LDA的目的是讓投影后的兩類之間離得更遠，類內離得更近。說到類內離得更近能想到的最直接的方法就是讓類內方差最小，這正是類內散布矩陣；說到類間離得更遠能想到的最直接的方法肯定是讓均值之差最大，而不是均值之差與自己的克羅內克積這個奇怪的東西最大。

那麼經典LDA為什麼會用類間散布矩陣呢？我個人認為是這樣的表達式看起來更加優雅：分子分母是齊次的，並且這個東西恰好就是廣義瑞利商：

$frac{x^TAx}{x^TBx}$

雖然經典LDA求解沒有上面這個方法直接，但是問題表述更加規範，所用到的技巧也非常簡單，不會給是使用者帶來困擾，所以LDA最終採用的就是類間散布矩陣了吧。

2. 直線&超平面

上面那個問題只算是小打小鬧，沒有太大的意義，但是這個問題就很有意義了：LDA為什麼直接投影到直線（一維），而不能像PCA一樣投影到超平面（多維）呢？

我們試試不就完了hhh。

假設將樣本投影到 $W = (w_1, w_2, ..., w_c)in R^{d imes c}$ 上，其中每一個 $w_i$ 都是經典LDA中的 $w$ 。也就相當於我們不是把樣本投影到一條直線上，而是投影到c條直線上，也就相當於投影到了超平面上。投影后的樣本坐標為：

$W^Tx = left[ egin{matrix} w_1^Tx \ w_2^Tx \ . \ . \ w_c^Tx end{matrix} ight] in R^c$

所以樣本的投影過程就是：

$xRightarrow W^Tx$

那麼，均值的投影過程也是這樣：

$m_1Rightarrow W^Tm_1 \ m_2Rightarrow W^Tm_2$

投影后的均值之差的二範數：

$egin{equation} egin{aligned} ||W^T(m_1-m_2)||_2^2 &=(m_1-m_2)^TWW^T(m_1-m_2) \ &=tr[W^T(m_1-m_2)(m_1-m_2)^TW] \ &=tr[W^TS_BW] end{aligned} end{equation}$

為什麼不用第一行的向量內積而偏要用第二行的跡運算呢？因為這樣可以拼湊出類間散布 $(m_1-m_2)(m_1-m_2)^T$ 來，和經典LDA保持形式的一致。

回顧一下經典LDA的形式：

$max limits_w quad w^TS_Bw \ qquad s.t. quad w^TS_Ww = 1$

現在我們有了 $tr[W^TS_BW]$ ，還缺個約束，類比一下就可以得到 $tr[W^TS_WW]$ 了：

$max limits_W quad tr[W^TS_BW] \ qquad s.t. quad tr[W^TS_WW] = 1$

實際上，約束也可以選擇成 $W^TS_WW = I_c$ ，這兩個約束實際上都是在限制W的解空間，得出來的解是等價的，這兩個約束有什麼區別呢？我只發現了一點：

回想 $W = (w_1, w_2, ..., w_c)in R^{d imes c}$ 是c條投影直線，為了確保向這c條直線投影能等價於向c維子空間投影，我們需要保證c條直線是線性無關的，即 $rank(W) =c$ 。看一下約束：

$W^TS_WW = I_c$

右邊是秩為c的單位矩陣，因為矩陣乘積的秩不大於每一個矩陣的秩，所以左邊這三個矩陣的秩都要不小於c，因此我們得到了 $rank(W) = c$ 。也就是說， $W^TS_WW = I_c$ 能夠保證我們在向一個c維子空間投影。而約束 $tr[W^TS_WW] = 1$ 中沒有顯式地表達這一點。我對矩陣的理解還不夠深入，不知道 $tr[W^TS_WW] = 1$ 是否也隱含了對秩的約束，所以為了保險起見，我選擇了帶有顯式秩約束的 $W^TS_WW = I_c$ ，這樣就得到了我們的高維投影版LDA：

$max limits_W quad tr[W^TS_BW] \ quad s.t. quad W^TS_WW = I_c \$

下面來求解這個問題。還是拉格朗日乘子法：

$L(W,Lambda) = -tr[W^TS_BW] + tr[Lambda^T(W^TS_WW-I_c)]$

求導並令其為0：

$partial_W L(W,Lambda) = -2S_BW + 2S_WWLambda = 0$

得到了：

$S_BW = S_WWLambda$

在大部分情況下，一些協方差矩陣的和 $S_W$ 是可逆的。即使不可逆，上面這個也可以用廣義特徵值問題的方法來求解，但是這裡方便起見我們認為 $S_W$ 可逆：

$S_W^{-1}S_BW = WLambda$

我們只要對 $S_W^{-1}S_B$ 進行特徵值分解，就可以得到d個特徵向量了，挑出最大特徵值對應的c個特徵向量來組成W，我們不就得到向c維子空間的投影了嗎？！

真的是這樣嗎？

不是這樣的。誠然，我們可以選出c個特徵向量，但是其中只有1個特徵向量真正是我們想要的，另外c-1個沒有意義。

觀察：

$S_W^{-1}S_B = S_W^{-1}(m_1-m_2)(m_1-m_2)^T$

發現了嗎？等式右邊的 $(m_1-m_2)$ 是一個向量，換句話說，是一個秩為1的矩陣。那麼，這個乘積的秩也不能大於1，並且它不是0矩陣，所以

$rank(S_W^{-1}S_B ) = 1$

秩為1的矩陣只有1個非零特徵值，也只有1個非零特徵值對應的特徵向量w。

可能有人會問了，那不是還有零特徵值對應的特徵向量嗎，用它們不行嗎？

不行。來看一下目標函數：

$max limits_W quad tr[W^TS_BW]$

我們剛才得到的最優性條件：

$S_BW = S_WWLambda$

所以目標函數為：

$egin{equation} egin{aligned} max limits_W quad tr[W^TS_BW] &= max limits_W quad tr[W^TS_WWLambda] \ &= max limits_W quad tr[I_cLambda] \ &= max limits_W quad tr[Lambda] \ &= max limits_W quad lambda_1 + lambda_2 + ... + lambda_d end{aligned} end{equation}$

而我們的W只能保證 $lambda_1, lambda_2, ..., lambda_d$ 中的一個非零，無論我們怎麼選取剩下的c-1個w，目標函數也不會再增大了，因為唯一一個非零特徵值對應的特徵向量已經被選走了。

所以，兩類LDA只能向一條直線投影。

這裡順便解釋一下K類LDA為什麼只能投影到K-1維，其實道理是一樣的。K類LDA的類間散布矩陣是：

$S_B=sum_{k=1}^KN_k(m_k-m)(m_k-m)^T$

可以看出， $S_B$ 是K個秩一矩陣 $(m_k-m)(m_k-m)^T$ 的和（因為 $m_k-m$ 是秩一的向量），所以它的秩最大為K。並且 $Nm = N_1m_1 + N_2m_2 + ... + N_km_k$ ，所以這K項中有一項可以被線性表出。所以， $S_B$ 的秩最大為K-1。也即：

$rank(S_W^{-1}B_W) = K-1$

只有K-1個非零特徵值。所以，K類LDA最高只能投影到K-1維。

咦？剛才第三個問題怎麼提的來著，可不可以不用比值用差值，用差值的話會不會解決這個投影維數的限制呢？

3. 比值&差值

經典LDA的目標函數是投影后的類間散布與類內散布的比值，我們很自然地就會想，為什麼非得用比值呢，差值有什麼不妥嗎？

再試試不就完了hhh。

注意，這一節我們不用向量w，而使用矩陣W來討論，這也就意味著我們實際上在同時討論二類LDA和多類LDA，只要把 $S_B$ 和 $S_W$ 換成對應的就好了。

$max limits_W quad tr[W^TS_BW] - alpha tr[W^TS_WW]$

注意到可以通過放縮W來得到任意大的目標函數，所以我們要對W的規模進行限制，同時也進行秩限制：

$W^TW=I_c$

也就得到了差值版的LDA：

$max limits_W quad tr[W^TS_BW] - alpha tr[W^TS_WW]\ s.t. quad W^TW=I_c qquad qquad qquad qquad$

依然拉格朗日乘子法：

$L(W, Lambda) = -tr[W^TS_BW]+alpha tr[W^TS_WW] + tr[Lambda^T(W^TW-I_c)]$

求導並令其為0：

$partial_W L(W,Lambda) = -2S_BW + 2alpha S_WW+2WLambda = 0$

得到了：

$S_BW =alpha S_WW+WLambda$

由 $W^TW=I_c$ ，有：

$S_BW =alpha S_WW+WLambda W^TW$

可以得到：

$S_BW=(alpha S_W+WLambda W^T)W$

若括弧內的東西可逆，則上式可以寫為：

$(alpha S_W+WLambda W^T)^{-1}S_BW=W$

注意到， $S_B$ 的秩不大於K-1，說明等號左邊的秩不大於K-1，那麼等號右邊的秩也不大於K-1，即

$rank(W) <= K-1$

所以我們還是會遇到秩不足，無法求出K-1個以上的非零特徵值和對應的特徵向量。這樣還不夠，我們還需要證明的一點是，新的目標函數在零特徵值對應的特徵向量下依然不會增加。

目標函數（稍加變形）為：

$max limits_W quad tr[W^T(S_B-alpha S_W)W]$

再利用剛才我們得到的最優性條件（稍加變形）：

$(S_B -alpha S_W)W=WLambda$

所以目標函數為：

$egin{equation} egin{aligned} max limits_W quad tr[W^T(S_B-alpha S_W)W] &= max limits_W quad tr[W^TWLambda] \ &= max limits_W quad tr[I_cLambda] \ &= max limits_W quad tr[Lambda] \ &= max limits_W quad lambda_1 + lambda_2 + ... + lambda_d end{aligned} end{equation}$