答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]

This entry was posted in 答讀者問, 應用之道 and tagged , , , . Bookmark the permalink.

31 則回應給 答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