答Regan Yuan──關於主成分分析的奇異值分解算法

網友Regan Yuan留言:

老师您好,首先对您以往的支持和耐心详细的讲解,表示由衷的敬意和感谢!再次虚心请教老师一个问题,对于采用singular value decomposition的principal components analysis算法,最近总是有些搞不清,能否提供一个具体的例子呢?比如下面这个问题:五个学生 (5 cases or observations) 的跳高,跳远,乒乓球三门具有相关性的体育成绩 (3 variables or dimensions) 的 5\times 3 矩阵,

\displaystyle  X=\begin{bmatrix}  6&5&8\\  7&7&6\\  8&7&7\\  7&8&6\\  9&9&8  \end{bmatrix}

用PCA进行降维度,具体解决方法如何?请明示,谢谢您!祝您开心每一天!

 
答曰:

為了相互參照,底下就你提供的數據矩陣說明兩種主成分分析的計算方法:樣本相關係數矩陣的譜分解 (spectrum decomposition) 和標準化後的數據矩陣的奇異值分解。詳細討論請見“主成分分析”和“主成分分析與奇異值分解”,在此不再贅述。

 
譜分解

令跳高、跳遠和乒乓球三項成績為變數 x_1,x_2,x_3。首先算出三個變數的樣本平均數 (見“樣本平均數、變異數和共變異數”)

\displaystyle  m_1=7.4,~~~m_2=7.2,~~~m_3=7

和樣本變異數

\displaystyle  s_1^2=1.3,~~~s_2^2=2.2,~~~s^2_3=1

故樣本標準差為 s_1=1.140s_2=1.483s_3=1。將數據矩陣 X=[x_{ij}] 予以標準化,算式為

\displaystyle  \tilde{x}_{ij}=\frac{x_{ij}-m_j}{s_j}

我們以 \tilde{X}=[\tilde{x}_{ij}] 表示標準化後的數據矩陣,如下:

\displaystyle  \tilde{X}=\left[\!\!\begin{array}{rrr}  -1.228&-1.483&1.000\\  -0.351&-0.135&-1.000\\  0.526&-0.135&0.000\\  -0.351&0.539&-1.000\\  1.403&1.214&1.000  \end{array}\!\!\right]

標準化後的數據集的樣本共變異數矩陣即為樣本相關矩陣 (見“相關係數”),故得 R=[r_{ij}]=\frac{1}{n-1}\tilde{X}^T\tilde{X},其中 n 為樣本大小,此例 n=5。計算得到

\displaystyle  R=\left[\!\!\begin{array}{crr}  1.000&0.828&0.219\\  0.828&1.000&-0.169\\  0.219&-0.169&1.000  \end{array}\!\!\right]

可知跳高和跳遠成績有頗強的正相關性 (r_{12}=0.828),跳高和乒乓球成績為正相關 (r_{13}=0.219),但跳遠和乒乓球成績為負相關 (r_{23}=-0.169)。相關矩陣 R 為一實對稱正定矩陣,其特徵值為 \lambda_1=1.829\lambda_2=1.081\lambda_3=0.090,對應的單範正交 (orthonormal) 特徵向量 (稱為主軸) 分別為

\displaystyle  \mathbf{w}_1=\begin{bmatrix}  0.712\\  0.701\\  0.046  \end{bmatrix},~~\mathbf{w}_2=\left[\!\!\begin{array}{r}  -0.172\\  0.237\\  -0.956  \end{array}\!\!\right],~~\mathbf{w}_3=\left[\!\!\begin{array}{r}  0.681\\  -0.673\\  -0.289  \end{array}\!\!\right]

注意,\text{trace}R=\lambda_1+\lambda_2+\lambda_3=3。令 D=\text{diag}(\lambda_1,\lambda_2,\lambda_3)W=\begin{bmatrix}  \mathbf{w}_1&\mathbf{w}_2&\mathbf{w}_3  \end{bmatrix} 為一正交矩陣 (orthogonal matrix)。主成分係數矩陣即為

\displaystyle  Z=\tilde{X}W=\left[\!\!\begin{array}{rrr}  -1.868& -1.096& -0.127\\  -0.390&  0.984&  0.141\\   0.280& -0.122&  0.449\\   0.082&  1.144& -0.313\\   1.896& -0.910& -0.150  \end{array}\!\!\right]

因素負荷矩陣為

\displaystyle  F=WD^{\frac{1}{2}}=\left[\!\!\begin{array}{crr}  0.963&-0.178&0.202\\  0.948& 0.246&-0.202\\  0.062&-0.994&-0.088  \end{array}\!\!\right]

 
奇異值分解

設標準化後的數據矩陣 \tilde{X} 的「瘦」奇異值分解為 \tilde{X}=U\Sigma V^T (見“奇異值分解 (SVD)”),計算可得

\displaystyle\begin{aligned}  U&=\left[\!\!\begin{array}{rrr}  -0.691&-0.527&-0.212\\  -0.144&0.474&0.235\\  0.104&-0.059&0.751\\  0.031&0.550&-0.524\\  0.701&-0.438&-0.249  \end{array}\!\!\right]\\  \Sigma&=\begin{bmatrix}  2.705&&\\  &2.079&\\  &&0.598  \end{bmatrix}\\  V^T&=\left[\!\!\begin{array}{rrr}  0.712&0.701&0.046\\  -0.172&0.237&-0.956\\  0.681&-0.673&-0.289  \end{array}\!\!\right].\end{aligned}

奇異值分解提供了主成分係數矩陣 Z=U\Sigma,因素負荷矩陣 F=\frac{1}{\sqrt{n-1}}V\Sigma,主軸矩陣 W=V,以及特徵值矩陣 \text{diag}(\lambda_1,\lambda_2,\lambda_3)=\frac{1}{n-1}\Sigma^2

 
解析

  1. 需要保留幾個主成分?對角矩陣 \Sigma 的奇異值 (主對角元) \sigma_j 可推得樣本相關係數的特徵值 \lambda_j=\frac{\sigma_j^2}{n-1}j=1,2,3。特徵值 \lambda_j 代表第 j 個主成分的變異數 (即訊息量),\lambda_1=1.829\lambda_2=1.081\lambda_3=0.090 說明了包含三個變數的數據矩陣 \tilde{X} 僅須選取前二個主成分 (一般保留特徵值大於 1 的成分),故而保留 (\lambda_1+\lambda_2)/(\lambda_1+\lambda_2+\lambda_3)=2.91/3=97% 的變異量 (即經過標準化後的訊息量)。
  2. 降維後的數據資料有何意義?主成分係數矩陣 Z 的第一行 (column) 和第二行表示 5 位學生在第一主軸和第二主軸的座標,也就是第一主成分和第二主成分的含量。觀察發現第三行的數值相當小,故第三主成分可視為雜音。
  3. 主成分代表甚麼訊息?因素負荷矩陣 F(i,j)f_{ij} 揭露第 i 個變數和第 j 個主成分的相關係數。從 f_{11}=0.963f_{21}=0.948f_{31}=0.062 得知第一個主成分與跳高和跳遠兩項成績極為相關,但與乒乓球成績無關;從 f_{12}=-0.178f_{22}=0.246f_{32}=-0.994 得知第二個主成分與跳高和跳遠成績相關性不大,它反映的是乒乓球成績 (反轉 \mathbf{w}_2 方向即可改變 Z 的第二行正負號,並得到正相關係數)。因此,第一個主成分代表彈跳能力,第二個主成分代表肢體協調。
  4. 從所選取的主成分如何得到近似的數據矩陣?假設我們保留了前二個主成分,即刪除主成分係數矩陣 Z 的第三行。令

    \displaystyle  Z'=\left[\!\!\begin{array}{rrc}  -1.868& -1.096&0\\  -0.390&  0.984&0\\   0.280& -0.122&0\\   0.082&  1.144&0\\   1.896& -0.910&0  \end{array}\!\!\right]

  5. 因為 W^T=W^{-1},可得近似的標準化數據矩陣

    \displaystyle  \tilde{X}'=Z'W^T=Z'\begin{bmatrix}  \mathbf{w}_1^T\\  \mathbf{w}_2^T\\  \mathbf{w}_3^T  \end{bmatrix}=\begin{bmatrix}  -1.868\mathbf{w}_1^T-1.096\mathbf{w}_2^T\\  -0.390\mathbf{w}_1^T+0.984\mathbf{w}_2^T\\   0.280\mathbf{w}_1^T-0.122\mathbf{w}_2^T\\   0.082\mathbf{w}_1^T+1.144\mathbf{w}_2^T\\   1.896\mathbf{w}_1^T-0.910\mathbf{w}_2^T  \end{bmatrix}=\left[\!\!\begin{array}{rrr}  -1.142 & -1.569 & 0.962\\  -0.447 & -0.040 &-0.959\\   0.220 & 0.167  & 0.130\\  -0.138 & 0.329 & -1.090\\   1.506 & 1.113 & 0.957  \end{array}\!\!\right]

    再將近似的標準化數據矩陣 \tilde{X}'=[\tilde{x}'_{ij}] 轉換為原始尺規表達的近似數據矩陣 X'=[x'_{ij}]。套用公式 x'_{ij}=\tilde{x}'_{ij}s_j+m_j,或寫成矩陣運算:

    \displaystyle\begin{aligned}  X'&=\tilde{X}'\begin{bmatrix}  s_1&&\\  &s_2&\\  &&s_3  \end{bmatrix}+\begin{bmatrix}  m_1\mathbf{1}&m_2\mathbf{1}&m_3\mathbf{1}  \end{bmatrix}\\  &=\tilde{X}'\begin{bmatrix}  1.140&&\\  &1.483&\\  &&1  \end{bmatrix}+\begin{bmatrix}  7.4&7.2&7\\  7.4&7.2&7\\  7.4&7.2&7\\  7.4&7.2&7\\  7.4&7.2&7  \end{bmatrix}=\begin{bmatrix}  6.098 &4.873 &7.962\\  6.890 &7.141 &6.041\\  7.651 &7.448 &7.130\\  7.243 &7.688 &5.910\\  9.117 &8.851 &7.957  \end{bmatrix},\end{aligned}

    其中 \mathbf{1}=(1,1,1,1,1)^T。上式可解釋為伸縮與平移數據點 (見“數據矩陣的列與行”)。近似誤差為

    \displaystyle  X-X'=\left[\!\!\begin{array}{rrr}  -0.098  &0.127  &0.038\\   0.110  &-0.141 &-0.041\\   0.349  &-0.448 &-0.130\\  -0.243  &0.312  &0.090\\  -0.117  &0.149  &0.043  \end{array}\!\!\right]

廣告
本篇發表於 答讀者問, 應用之道 並標籤為 , , , 。將永久鏈結加入書籤。

31 Responses to 答Regan Yuan──關於主成分分析的奇異值分解算法

  1. Regan Yuan 說道:

    老师您好,再次感谢,昨天用matlab软件 中的princomp实验了一下,发现结果略有不同,不知道什么原因,此外我用了这个文章《A tutorial onprincipal components analysis》中的代码based on SVD:http://arxiv.org/pdf/1404.1100v1.pdf,发现princomp和这个文章的两个代码的运行结果一样,请教是不是这里的算法和您所提出的算法,是不是设置上有所不同,我还是有些困惑。请老师赐教。谢谢 附录代码: function [signals,PC,V] = pca2(data)   % PCA2: Perform PCA using SVD.   % data – MxN matrix of input data   % (M dimensions, N trials)   % signals – MxN matrix of projected data   % PC – each column is a PC   % V – Mx1 matrix of variances   [M,N] = size(data);   % subtract off the mean for each dimension   mn = mean(data,2);   data = data – repmat(mn,1,N);   % construct the matrix Y   Y = data’ / sqrt(N-1);   % SVD does it all   [u,S,PC] = svd(Y);   % calculate the variances   S = diag(S);   V = S .* S;   % project the original data   signals = PC’ * data;

    • ccjou 說道:

      PCA演算程序的輸入數據矩陣可能僅去除平均數(稱為demean)或者標準化(demean後再除以樣本標準差),前者得到樣本共變異數矩陣,後者產生樣本相關矩陣。你給出的文獻採用樣本共變異數矩陣,本文使用的是相關矩陣,兩者的PCA結果不同。如果變數的變異數差異很大,通常我們會選用樣本相關矩陣,詳細討論見下文的最後一段:
      https://ccjou.wordpress.com/2013/04/15/%E4%B8%BB%E6%88%90%E5%88%86%E5%88%86%E6%9E%90/

      Matlab 文件說明:
      princomp centers X by subtracting off column means, but does not rescale the columns of X. To perform principal components analysis with standardized variables, that is, based on correlations, use princomp(zscore(X)).

  2. Regan Yuan 說道:

    老师您好,请问上次的问题,是否收到?能否给我指教,谢谢

  3. wang 說道:

    老師您好,針對文忠寫道:『對應的標準正交 (orthonormal) 特徵向量 (稱為主軸) 分別為W1 W2 W3』 怎麼度應的呢?謝謝

  4. Wang 說道:

    老師您好,我想針對本題目稍作修改,請您幫我看看!
    原始舉陣X(5位學員,3門具有相關性的體育成績)若將題目解釋成(5種果茶,3種原料比例)
    3種原料比例,我設z1、z2、z3(因此第一杯果茶分別由6*z1+5*z2+8*z3所組成;第二杯果茶分別由7*z1+7*z2+6*z3所組成…以此類推!)

    若我用主成分分析X
    (a)使用樣本共變異數矩陣,僅去除平均值。(矩陣分別減去3種原料的平均值)
    ※關於僅去除平均值(口與或許能解釋成讓原始舉陣轉換座標以平均值為基準)但如此一來原始舉陣便會改變(各比例關係已和原始舉陣無關)
    Xmean=[7.4 7.2 7]
    原始舉陣X=[-1.4 -2.2 1;-0.4 -0.2 -1;0.6 -0.2 0;-0.4 0.8 -1;1.6 1.8 1](各比例關係已和原始舉陣無關)
    若還是硬著頭皮做分析,得W1=[0.5865 0.8096 -0.0251]、W2=[0.3181 -0.2017 0.9264]、W3=[0.7449 -0.5513 -0.3758]
    W1約占71.6%、W2約占25.3%、W3約占3.06%
    ※※這跟原始舉陣直接作分析,得到的結果〝竟然〞一樣???這該如何解釋呢?這樣做的目的是?
    明明分析時原比例6:5:8扣除平均值後變成-1.4 -2.2 1,不了解這中間的關係!

    (b)使用樣本相關矩陣
    ※算法如同老師的參數一樣
    得W1=[0.7117 0.7010 0.0457]、W2=[0.1716 -0.2365 0.9564]、W3=[0.6812 -0.6728 -0.2886]
    W1約占60.9826%、W2約占36.0340%、W3約占2.9834%

    最後想再問,若直接用數據來看W1+W2=97%(那麼說明了W1和W2這2種線性比例關係能夠完整表達原始矩陣的97%訊息)

    W1=[0.7117 0.7010 0.0457]皆是正相關,每次調配依照此比例=某一種果茶
    W2=[0.1716 -0.2365 0.9564]有負相關,這點就不知道怎麼去做調配了,理論上依這線性比例,也會是一種果茶

    而用近式解求得出(W1*8.1404)+(W2*7.4973)=原矩陣6*z1+5*z2+8*z3

    還請老師指導以上觀念有誤的地方,謝謝!

  5. Wang 說道:

    s1跑出來和SVD分解後的轉置V 一樣

    • Wang 說道:

      以本題為例
      SVD 分解後的V transpose
      矩陣=[0.712 0.701 0.046;-0.172 0.237 -0.956;0.681 -0.673 -0.289]

      但有後正負號卻會剛好相反
      譬如PCA後矩陣=[0.712 0.701 0.046;0.172 -0.237 0.956;-0.681 0.673 0.289]

      • ccjou 說道:

        但 V 不就是由特徵向量所組成的矩陣嗎?V就是上文的W。

      • Wang 說道:

        data=[6 5 8 ; 7 7 6 ; 8 7 7 ; 7 8 6 ;9 9 8]
        [s1 t1 l1]=princomp(data)
        s1 =
        0.5865 0.3181 0.7449
        0.8096 -0.2017 -0.5513
        -0.0251 0.9264 -0.3758
        t1 =
        -2.6272 0.9248 -0.2059
        -0.3714 -1.0132 0.1881
        0.1900 0.2312 0.5572
        0.4382 -1.2150 -0.3631
        2.3705 1.0722 -0.1763
        l1 =
        3.2219
        1.1403
        0.1378

        dataMean =
        7.4000 7.2000 7.0000
        B = data-repmat(dataMean,[n 1])
        B =
        -1.4000 -2.2000 1.0000
        -0.4000 -0.2000 -1.0000
        0.6000 -0.2000 0
        -0.4000 0.8000 -1.0000
        1.6000 1.8000 1.0000
        [s2 t2 l2]=princomp(B)
        s2 =
        0.5865 0.3181 0.7449
        0.8096 -0.2017 -0.5513
        -0.0251 0.9264 -0.3758
        t2 =
        -2.6272 0.9248 -0.2059
        -0.3714 -1.0132 0.1881
        0.1900 0.2312 0.5572
        0.4382 -1.2150 -0.3631
        2.3705 1.0722 -0.1763
        l2 =
        3.2219
        1.1403
        0.1378

  6. Wang 說道:

    資質愚鈍,容許我再度提問!
    上述老師解釋的w1、w2、w3(分別對應特徵值123);100*(特徵)/特徵1+2+3=X%
    很多書上也幾乎只寫道這個地方(當好對應Princomp指令),
    才會想說(若6成=[0.712 0.701 0.046],可以依這種線性組合的比例來合成該有多好,但負相關卻不知如何調配!!!呵呵~!)
    上述這些稱特徵向量(主軸)…書本解釋為:主成分
    那麼Z矩陣又稱為主成分係數,有點混亂了!
    Z矩陣=標準化後乘上特徵值=?
    若只保留前2個成分係數
    想問一下若『單獨』來看成分1及成分2,這2個各自代表的又是什麼呢?

    • ccjou 說道:

      上文解析3說明了如何從因素負載矩陣解讀主成分的意義,對應的特徵值表示該主成分的強度。第一個主成分表現彈跳能力(跳高與跳遠),第二個主成分表現肢體協調(乒乓球)。特徵向量各元之所以出現負值是因為我們以平均量作為參考基準,如果用果茶當例子,以中等(平均)甜度的果茶當作基準,則負值表示減少糖量,正值表示增加糖量。

  7. 龙一 說道:

    老师解析3用因素负载说明了主成分的强度,文中算至特征向量主轴时也说明了正负相关,如何解释这两者差异?而最期待获得的解析是Wang所询问的如何同时考虑所有变数,正如Wang第一次提问,原矩阵为5X3,主成分应有3个3X1所组成,依特征值决定其重要性,方可得知需选成个数,负值理应是取决标准化是否低于样本平均,令x_1柑橘x_2薄荷x_3茉莉,如何诠释成分1至3完整调配的比例,并非单独考虑,以100c.c果茶解析成分关系,w_1’=[0.712 0.701 0.046]=48.8c. c柑橘+48c. c薄荷+3.2c. c茉莉=成分1,w_2’=[-0.172 0.237 -0.956]=?还是需先加回样本平均?

    • ccjou 說道:

      我不是很清楚你的問題。一個一個來吧。首先,如何解释这两者差异?哪兩者?至於你的第二個問題,請先閱讀解析4後再重寫過問題。

    • ccjou 說道:

      因素負載矩陣 F 由伸縮特徵向量矩陣 W 的行 (column) 而得,F=WD^{1/2}=\begin{bmatrix} \sqrt{\lambda_1}\mathbf{w}_1&\sqrt{\lambda_2}\mathbf{w}_2&\sqrt{\lambda_3}\mathbf{w}_3\end{bmatrix}

      • 龙一 說道:

        成分1成分2成分3各是什么型态,而并非概略得知,w_1’=[0.712 0.701 0.046]各变数正负相关性,以及成分1说明了x_1柑橘0.712与x_2薄荷0.701影响程度较为相关,5个样本对应的3种配方变量5X3,而今告知成分1可表达61%原矩阵讯息,换言之原矩阵61%皆占有成分1的成分,如今却仅以解释在成分1内3变量各自对于成分1的重要性

        • ccjou 說道:

          我仍不明白你的問題。這樣吧,我寫出配方公式給你參考。符號如上文所述,數據矩陣代表五種果茶A,B,C,D,E的柑橘、薄荷、茉莉含量,果茶A為 (x_1,x_2,x_3)=(6,5,8),標準化如下:

          (\tilde{x}_1,\tilde{x}_2,\tilde{x}_3)=(\frac{x_1-m_1}{s_1},\frac{x_2-m_2}{s_2},\frac{x_3-m_3}{s_3})=(-1.228,-1.483,1)

          我想利用主成分簡化果茶調製過程,我知道第一、二、三主成分的變異量分別為61%、36%、3%,故可忽略第三主成分,也就是果茶口味不同主要決定於第一和第二主成分的添加量。果茶A的主成分係數為 (z_1,z_2,z_3)=(-1.868,-1.096,-0.127),於是設簡化配方為 (z_1,z_2,0)=(-1.868,-1.096,0)。首先得到近似的標準化數據

          (\tilde{x}'_1,\tilde{x}'_2,\tilde{x}'_3)=z_1\mathbf{w}_1^T+z_2\mathbf{w}_2^T=(-1.868)(0.712,0.701,0.046)+(-1.096)(-0.172,0.237,-0.956)=(-1.142,-1.569,0.962)

          最後可得簡化配方的柑橘、薄荷、茉莉含量

          (x'_1,x'_2,x'_3)=(s_1\tilde{x}'_1+m_1,s_2\tilde{x}'_2+m_2,s_3\tilde{x}'_3+m_3)=(6.098,4.873,7.962)\sim(6,5,8)

  8. 洪嘉隆 說道:

    老師你好!我想請問我將數據取共變異數矩陣然後取特徵值,或是直接用matlab的函數PCA我得到的值,正負號所代表的意思是甚麼?謝謝老師

    • ccjou 說道:

      共變異數矩陣是半正定,特徵值必不為負。我不知道你說的PCA得到的值是甚麼變數的值?

發表迴響

在下方填入你的資料或按右方圖示以社群網站登入:

WordPress.com Logo

您的留言將使用 WordPress.com 帳號。 登出 / 變更 )

Twitter picture

您的留言將使用 Twitter 帳號。 登出 / 變更 )

Facebook照片

您的留言將使用 Facebook 帳號。 登出 / 變更 )

Google+ photo

您的留言將使用 Google+ 帳號。 登出 / 變更 )

連結到 %s