Representation

Representation / Factorization

「重新表示」。數據套用函數,轉換成其他姿態,重新呈現。

「分解」。重新表示的反向操作,揭露數據原本姿態。

重新表示:已知原數據,求新數據、求函數。

分解:已知新數據,求原數據、求函數。

當函數擁有反函數,只要移項一下,重新表示、分解其實是同一個問題,沒有主從之分。

Matrix Representation / Matrix Factorization

當函數是矩陣,則有三種套用函數的方式。

一、矩陣乘在數據左邊:數據們實施線性變換。

二、矩陣乘在數據右邊:數據們的加權平均值。

三、矩陣相加:數據們加上誤差。

數據可以是一維,甚至多維。數據可以是一筆,甚至多筆。

一筆數據是一個直條。當矩陣乘在左邊,原數據與新數據從左往右依序一一對應。當矩陣乘在右邊,權重與新數據一一對應。當矩陣相加,原數據、新數據、誤差一一對應。

不熟悉線性代數的讀者,請見本站文件「Basis」,複習一下「矩陣切成直條」的觀點。

備註

representation通常是指「數據的等價呈現方式」,而不是這裡提到的「數據套用函數變成新數據」。

factorization通常是指「乘法的反向操作」,而不是這裡提到的「新數據分解成函數和原數據」。

這裡的定義源自機器學習領域,而該領域常常胡亂造詞。按理應當另造新詞,然而我並未發現更適切的詞彙,只好姑且沿用。

備註

矩陣乘法其實是函數複合。多個矩陣連乘,併成一個矩陣,這是「函數複合composition」。一個矩陣,拆成多個矩陣連乘,這是「函數分解decomposition」。函數複合的反向操作是函數分解。

儘管representation factorization外觀如同composition decomposition,但是其意義是數據、而非函數。本質略有差異。

Principal Component Analysis

Principal Component Analysis

Y = AX。已知數據X,求新數據Y、求轉換矩陣A。

YYᵀ = I。令新數據Y的維度是正規正交:相同維度的點積是1,相異維度的點積是0。既是單位向量、又互相垂直。

一個維度視作一個向量;兩個向量的點積,得到一個值;所有兩兩向量的點積,排列成矩陣。

推導過程

{ Y = AX
{ YYᵀ = I
given X. find A and Y.
                       _____________________________________________ 
YYᵀ             = I   | XXᵀ is a real matrix.                       |
AX(AX)ᵀ         = I   | XXᵀ is a square matrix.                     |
AXXᵀAᵀ          = I   | XXᵀ is a symmetric matrix.                  |
A(XXᵀ)Aᵀ        = I   | XXᵀ is a positive semidefinite matrix.      |
A(EΛE⁻¹)Aᵀ      = I   | thus XXᵀ has real non-negative eigenvalues. |
A(E√Λ)(√ΛE⁻¹)Aᵀ = I   | thus XXᵀ has orthonormal eigenbasis.        |
A = √Λ⁻¹E⁻¹           | let XXᵀ = EΛE⁻¹ (E⁻¹ = Eᵀ)                  |
A = √Λ⁻¹Eᵀ             ‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾ 

XXᵀ是實數對稱半正定方陣,得以實施特徵分解。

E是特徵向量們(特徵基底),恰為正規正交矩陣,其轉置矩陣恰為反矩陣。而且都是實數。

Λ是特徵值們,恰為對角線矩陣。而且都是實數、都非負數。

欲求A,就想辦法讓A和Aᵀ能夠抵消EΛEᵀ,使之成為I。特徵基底求反矩陣,特徵值先開根號再求倒數,如此就湊出了A。

欲求Y,只消計算AX。

正確答案不只一種。特徵向量對調次序、顛倒方向,特徵值開根號時取正值、取負值,都是正確答案。

補充一下,正確答案的左邊乘上任意的正規正交矩陣,仍是正確答案。但由於這類答案並不實用,大家總是無視這類答案。

數據減去平均值

當X的中心(平均值)位於原點,則Y的中心也將位於原點。

大家總是預先將X減去平均值,將X的中心挪至原點。好處是:一、降低數值範圍,以減少浮點數運算誤差。二、具備直線擬合效果,詳見後面章節。

X減去平均值之後,XXᵀ和YYᵀ變成「維度的共變異矩陣」。

演算法(Eigendecomposition)

一、每筆數據減去其平均值。

二、求得維度之間(不是數據之間)的共變異矩陣。

三、共變異矩陣實施特徵分解。

X = {(1,2,3), (4,5,6), (5,0,4), (3,3,3), (7,5,9)}

    [ 1 4 5 3 7 ]   [ x₁ x₂ x₃ x₄ x₅ ]   [ |  |  |  |  |  ]
X = [ 2 5 0 3 5 ] = [ y₁ y₂ y₃ y₄ y₅ ] = [ p₁ p₂ p₃ p₄ p₅ ]
    [ 3 6 4 3 9 ]   [ z₁ z₂ z₃ z₄ z₅ ]   [ |  |  |  |  |  ]

    [ 4 4 4 4 4 ]
X̄ = [ 3 3 3 3 3 ]
    [ 5 5 5 5 5 ]

            [           ]   [ —— d₁ —— ]
X̂ = X - X̄ = [   3 × 5   ] = [ —— d₂ —— ]
            [           ]   [ —— d₃ —— ]

       [ d₁⋅d₁ d₁⋅d₂ d₁⋅d₃ ]   [       ]
X̂ X̂ᵀ = [ d₂⋅d₁ d₂⋅d₂ d₂⋅d₃ ] = [ 3 × 3 ] = E Λ Eᵀ
       [ d₃⋅d₁ d₃⋅d₁ d₃⋅d₃ ]   [       ]

    [       ]       [ λ₁ 0  0  ]
E = [ 3 × 3 ]   Λ = [ 0  λ₂ 0  ]
    [       ]       [ 0  0  λ₃ ]

              [ 1/√λ₁   0     0   ] [       ]
A = √Λ⁻¹ Eᵀ = [   0   1/√λ₂   0   ] [ 3 × 3 ]
              [   0     0   1/√λ₃ ] [       ]

推導過程

XXᵀ的特徵分解,正是X的奇異值分解的其中一個步驟。

XXᵀ的特徵值們開根號√Λ,正是X的奇異值們Σ。

XXᵀ的特徵向量們E,正是X的奇異向量們U。

X = UΣVᵀ
XXᵀ = UΣVᵀ(UΣVᵀ)ᵀ = UΣVᵀVΣᵀUᵀ = UΣΣᵀUᵀ = U(ΣΣᵀ)Uᵀ = EΛEᵀ
√Λ = Σ, E = U

以奇異值分解來詮釋PCA,簡潔直觀。X = UΣVᵀ,顯然Vᵀ的維度是正規正交,即為所求。令新數據Y = Vᵀ,然後將剩下的UΣ移項當作變換矩陣A = Σ⁻¹Uᵀ = √Λ⁻¹Eᵀ。

X = UΣVᵀ
UᵀX = ΣVᵀ     (U is orthonormal, thus U⁻¹ = Uᵀ)
Σ⁻¹UᵀX = Vᵀ
A = Σ⁻¹Uᵀ
Y = Vᵀ

演算法(Singular Value Decomposition)

一般來說,實數對稱半正定方陣XXᵀ的特徵分解,比普通矩陣X的奇異值分解來得快。可是當X是超大矩陣的情況下,預先計算XXᵀ相當費時,而奇異值分解不必計算XXᵀ,逆轉勝。

演算法(EM-PCA)(Direct Linear Transformation)

Generate random matrix A. Loop over until A converge to PCA basis:
  E-step:  X = A⁺Y  = (A Aᵀ)⁻¹ Aᵀ Y
  M-step:  A = Y X⁺ = X (XᵀX )⁻¹ Xᵀ 

明明是Fixed Point Iteration,結果被叫做EM Algorithm。

我不清楚這是否比奇異值分解來的快。

幾何意義

PCA找到了全新的垂直座標軸:原點是數據中心、座標軸是特徵基底、座標軸單位長度是特徵值平方根。

所有數據一齊位移、旋轉(與鏡射)、縮放(與鏡射)。

一、位移:數據中心位移至原點。

二、旋轉:以特徵基底,逆向旋轉數據。

三、縮放:各維度除以各特徵值平方根。

最後,各維度的平均數均為0,變異數均為1。

正確答案不只一種。特徵向量(連同特徵值)對調次序,效果是鏡射暨旋轉。特徵向量顛倒方向、特徵值變號,效果是鏡射。

正確答案可以包含鏡射,也可以不包含鏡射。如果討厭鏡射,就讓特徵基底成為右手座標系、讓特徵值平方根皆是正值。

補充一下,正確答案的左邊乘上任意的正規正交矩陣,也就是旋轉(與鏡射),仍是正確答案。但由於這類答案並不實用,大家總是無視這類答案。

讓特徵基底成為右手座標系

一、標準座標軸(單位矩陣I),是右手座標系。

二、右手座標系經過旋轉,仍是右手座標系。

三、偶數次鏡射,得視作旋轉。

四、右手座標系經過偶數次鏡射,仍是右手座標系。

五、右手座標系經過奇數次鏡射,若再增加一次鏡射,就是右手座標系。

determinant可以判斷鏡射次數。det(E) = +1是偶數次,是右手座標系。det(E) = -1是奇數次,是左手座標系;再增加一次鏡射,就是右手座標系,例如任取一個特徵向量顛倒方向(元素通通添上負號)。

幾何特性

一、所有數據投影到座標軸之後的變異數,第一座標軸最大,第二座標軸次大(須垂直於先前座標軸),……。

也就是說,座標軸的方向,就是數據最散開的方向。

二、所有數據到座標軸的距離平方和,第一座標軸最小,第二座標軸次小(須垂直於先前座標軸),……。此即直線擬合!

也就是說,座標軸的垂直方向,就是數據最聚合的方向。

座標軸總是穿越原點。當數據沒有預先減去平均值,則座標軸不會穿越數據中心,不過仍然具備前述性質。此即「必須穿越原點的直線擬合」!進一步改變原點,則可以製作「必須穿越任意點的直線擬合」。

數學證明

首先引入先備知識:Ax²求最大值、最小值。A是實數對稱半正定方陣。當x長度為1,答案是A的最大的、最小的特徵值暨特徵向量。

     T                 2
max x A x   subject ‖x‖ = 1
 x
     T            2
max x A x - λ (‖x‖ - 1)          Lagrange's multiplier
 x

∂     T            2
―― [ x A x - λ (‖x‖ - 1) ] = 0   derivative = 0
∂x

2 A x - 2 λ x = 0                expand (A is symmetric)

A x = λ x                        transposition

A 最大的特徵值暨特徵向量就是正解。

投影之後變異數最大、數據散開性證明:

                                 2
max ∑ ‖proj pᵢ - (1/n) ∑ proj pⱼ‖    subject to ‖v‖ = 1
 v  i    v             j   v

               2
max ∑ ‖proj pᵢ‖    subject to ‖v‖ = 1
 v  i    v

            2 
max ∑ ‖pᵢ∙v‖       subject to ‖v‖ = 1
 v  i

         T  2
max ∑ ‖pᵢ v‖       subject to ‖v‖ = 1
 v  i

      T  2
max ‖X v‖          subject to ‖v‖ = 1
 v

     T   T 
max v X X v        subject to ‖v‖ = 1
 v

     T    T
max v (X X ) v     subject to ‖v‖ = 1
 v

   T
X X 最大的特徵值暨特徵向量就是正解。

距離平方和最小、直線擬合、數據聚合性證明:

                    2
min ∑ ‖pᵢ - proj pᵢ‖     subject to ‖v‖ = 1
 v  i         v

                     2
min ∑ ‖pᵢ - (pᵢ∙v) v‖    subject to ‖v‖ = 1
 v  i

          2          2        2   2
min ∑ ‖pᵢ‖ - 2 ‖pᵢ∙v‖ + ‖pᵢ∙v‖ ‖v‖    subject to ‖v‖ = 1
 v  i

          2        2 
min ∑ ‖pᵢ‖ - ‖pᵢ∙v‖      subject to ‖v‖ = 1
 v  i

              2 
min ∑ - ‖pᵢ∙v‖           subject to ‖v‖ = 1
 v  i

            2 
max ∑ ‖pᵢ∙v‖             subject to ‖v‖ = 1
 v  i

   T
X X 最大的特徵值暨特徵向量就是正解。

Low-rank Principal Component Analysis

關於轉換矩陣A,先前討論方陣,現在討論矩陣。low-rank是指矩陣的rank比較小,外觀看起來是直條不足或者橫條不足。

轉換矩陣不再擁有反矩陣,但是仍擁有廣義反矩陣,重新表示、分解仍是同一個問題。為了清楚說明,以下將兩者皆列出。

重新表示版本:令新數據維度小於原數據維度,降低資訊量。

答案是A = √Λ⁻¹Eᵀ,保留大的、捨棄小的特徵值暨特徵向量。

分解版本:令原數據維度小於新數據維度,降低資訊量。

答案是A = E√Λ,保留大的、捨棄小的特徵值暨特徵向量。

如果式子再添加一個誤差項,稱作Factor Analysis。雞肋。

程式碼

幾何意義

PCA找到了一組新的垂直座標軸,而Low-rank PCA僅保留了特徵值較大的座標軸。

所有數據一齊位移、投影(與鏡射)、縮放(與鏡射)。

因為捨棄了一些座標軸,所以旋轉必須重新解讀為投影:數據投影至僅存的特徵基底。

p = RandomVariate[MultinormalDistribution[{0,0,0}, {{3,1,1},{1,2,1},{1,1,1}}], 12] / 3; ListPointPlot3D[p, PlotStyle -> PointSize[Large], BoxRatios->{1,1,1}]; p = {{-0.0561955,-0.00847422,0.453698},{-0.334564,-0.272408,-0.238724},{-0.567886,0.0607641,0.265588},{0.502573,-0.344719,-0.296151},{0.19767,0.450711,0.0407837},{-0.0795941,-0.316957,-0.129278},{0.388671,0.00273937,0.330277},{0.0718672,-0.0904262,0.213121},{0.0928513,-0.312574,0.213095},{0.0484546,0.251097,-0.522902},{0.0447417,0.575007,-0.0364518},{-0.308589,0.00523944,-0.293056}}; o = Mean[p]; p = p - Table[o, Length[p]]; e = Eigenvectors[N[Transpose[p] . p]]; normal = e[[3,All]]; q = p - (Dot[p, normal] * Table[normal, Length[p]]); l = Transpose[{p,q}]; pl = InfinitePlane[{0,0,0} , e[[1;;2,All]] ]; axis = {{{-2,0,0},{2,0,0}},{{0,-2,0},{0,2,0}},{{0,0,-2},{0,0,2}}}; Graphics3D[{Black, Thickness[0.015], Line[axis], Black, Specularity[White, 10], Sphere[p, 0.05], Thickness[0.025], CapForm["Butt"], RGBColor[255,192,0], Line[l], RGBColor[192,0,0], Opacity[0.5], EdgeForm[None], pl}, PlotRange -> {{-.6,.6},{-.6,.6},{-.6,.6}}, Boxed -> False] r = Transpose[{Dot[p, e[[1,All]]], Dot[p, e[[2,All]]]}]; Graphics[{Black, PointSize[0.05], Point[r]}, PlotRange -> {{-.6,.6},{-.6,.6},{-.6,.6}}, Boxed -> False]

幾何特性

除了先前提及的散開性、聚合性,又多了一個逼真性。

三、所有數據投影到座標軸空間的距離平方總和最小。

也就是說,座標軸所在位置,令數據投影之後的失真最少。

數學證明

首先引入先備知識:對角線矩陣D,套用正規正交矩陣Q,trace只會減少,或者不變(當正規正交矩陣Q是單位矩陣I)。

幾何觀點:向量們形成標準座標軸,長度不必是1。經過旋轉或翻轉,則會偏離標準座標軸,使得座標值減少,總和也減少。

max tr( Q D )

令 Q = I 以得到最大值

投影前後距離平方和最小、數據逼真性證明:

(Ǎ: remove some columns from square matrix A)

             T   2               T
min ‖ X - B̌ B̌ X ‖    subject to B̌ B̌ = I

                T  T         T
min tr( (X - B̌ B̌ X)  (X - B̌ B̌ X) )        trace of inner product

         T       T   T     T   T   T
min tr( X X - 2 X B̌ B̌ X + X B̌ B̌ B̌ B̌ X )   expand

         T     T   T                       T
min tr( X X - X B̌ B̌ X )                   B̌ B̌ = I

           T   T
min tr( - X B̌ B̌ X )                       remove constant

         T   T
max tr( X B̌ B̌ X )                         multiply -1

           T   T
max tr( B̌ B̌ X X )                         cyclic property of trace

           T     T
max tr( B̌ B̌ E D E )                       eigendecomposition

         T   T
max tr( E B̌ B̌ E D )                       cyclic property of trace

    T   T
令 E B̌ B̌ E = Ǐ 以得到最大值 Ď

令 B̌ = Ě 以得到最大值 Ď

   T
X X 較大的特徵值和特徵向量就是正解。

精髓

從不同角度看數據,時聚時散。不斷滾動數據,重新找一個視角,讓數據看起來最散開、最聚合。注意到,是聚散,不是寬窄。

延伸應用

whitening:實施PCA,位移、旋轉、縮放。將數據範圍大致調整成單位圓,方便比對數據。觀念類似於normalization。

orientation:實施PCA,位移、旋轉、略過縮放。功效是重設座標軸,擺正數據。

alignment:兩堆數據各自實施PCA,特徵值由大到小排序,特徵向量一一對應,得知數據對應關係。

embedding:實施Low-rank PCA,捨棄小的特徵值暨特徵向量,只做投影。功效是降維、壓縮,而且是誤差最小的方式。

Independent Component Analysis

Independent Component Analysis

Y = XA。已知數據X,求新數據Y、求轉換矩陣A。

YᵀY = I。令新數據Y是正規正交:相異數據的點積是0,相同數據的點積是1。

一筆數據視作一個向量;兩個向量的點積,得到一個值;所有兩兩向量的點積,排列成矩陣。

推導過程

{ Y = XA
{ YᵀY = I
given X. find A and Y.
                       _____________________________________________ 
YᵀY             = I   | XᵀX is a square matrix.                     |
(XA)ᵀXA         = I   | XᵀX is a symmetric matrix.                  |
AᵀXᵀXA          = I   | XᵀX is a real matrix.                       |
Aᵀ(XᵀX)A        = I   | XᵀX is a positive semidefinite matrix.      |
Aᵀ(EΛE⁻¹)A      = I   | thus XᵀX has real non-negative eigenvalues. |
Aᵀ(E√Λ)(√ΛE⁻¹)A = I   | thus XᵀX has orthonormal eigenbasis.        |
A = E √Λ⁻¹            | let XᵀX = EΛE⁻¹ (E⁻¹ = Eᵀ)                  |
                       ‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾ 

ICA也是預先將X減去平均值。X減去平均值之後,XᵀX和YᵀY變成「數據的共變異矩陣」。

Low-rank Independent Component Analysis

重新表示版本:令新數據數量小於原數據數量,降低資訊量。

答案是A = E√Λ⁻¹,保留大的、捨棄小的特徵值暨特徵向量。

分解版本:令原數據數量小於新數據數量,降低資訊量。

答案是A = √ΛEᵀ,保留大的、捨棄小的特徵值暨特徵向量。

PCA與ICA的特徵分解觀點

PCA:Y = AX,YYᵀ = I,XXᵀ = EΛEᵀ,A = √Λ⁻¹Eᵀ。

ICA:Y = XA,YᵀY = I,XᵀX = EΛEᵀ,A = E√Λ⁻¹。

PCA矩陣在左,維度獨立。ICA矩陣在右,數據獨立。

PCA與ICA的奇異值分解觀點

PCA:X = UΣVᵀ,(Σ⁻¹Uᵀ)X = Vᵀ,Y = Vᵀ,A = Σ⁻¹Uᵀ。

ICA:X = UΣVᵀ,X(VΣ⁻¹) = U,Y = U,A = VΣ⁻¹。

PCA是左邊移項。變換之後,橫條正規正交(維度獨立)。

ICA是右邊移項。變換之後,直條正規正交(數據獨立)。

推導了一大堆,精髓其實只有上面兩句話。

PCA與ICA互為轉置

原數據X、新數據Y、轉換矩陣A通通轉置之後實施PCA,等同於直接實施ICA。PCA和ICA對調亦如是。

Y  = A X   , Y Yᵀ = I     PCA
--------------------------------------------------
Yᵀ = AᵀXᵀ  , YᵀYᵀᵀ= I     PCA: transpose X Y A
Yᵀ = (XA)ᵀ , YᵀY  = I     PCA: transpose X Y A
Y  = X A   , YᵀY  = I     ICA

寫成數學公式就是PCA(Xᵀ)ᵀ = ICA(X)、PCA(X) = ICA(Xᵀ)ᵀ、PCA(X)ᵀ = ICA(Xᵀ),每個式子意義都相同。

當數據形成對稱矩陣,那麼ICA與PCA完全相同。然而一般不會遇到這種事,比中彩券還難。

備註:直式改成橫式

因為線性變換已經規定將矩陣乘在左邊,所以矩陣乘在右邊挺彆扭。於是數學家就想到,把左式右式一齊轉置,使得矩陣乘在左邊,方便討論;但是缺點是數據和矩陣須排列成橫條,更加彆扭。

幾乎所有ICA文獻都採用此形式。另外,除了規定新數據YᵀY = I,還額外規定原數據XᵀX = I,兩者又推導出轉換矩陣AᵀA = I。通通正規正交。

橫式的壞處是不符合線性變換。矩陣以直條為主,矩陣存放數據是以直條排列。橫式的好處是符合計算機資料結構。陣列以橫條為主,陣列儲存數據是以橫條排列。

Principal Component Pursuit(Under Construction!)

Principal Component Pursuit

Y = X + A。已知新數據Y,求數據X、求轉換矩陣A。

low-rank X。令原數據X的秩盡量小。

常見的、平淡無奇的地方。數據是其他數據的線性組合。宛如基因演算法,數據之間的組合,總是那幾套模式。

sparse A。令轉換矩陣A盡量稀疏。

罕見的、獨樹一格的地方。數據某些維度,只有自己有,而別人沒有。

為何可以這樣分解呢?目前沒有數學系統可以完整詮釋。

http://www.math.umn.edu/~lerman/Meetings/CPCP.pdf

Convex Relaxation

提升一個次方,方便最佳化。

     min ‖ X ‖₀ + α ‖ A ‖₀     (α ≥ 0)
---> min ‖ X ‖⁎ + α ‖ A ‖₁     (α ≥ 0)
rank -> nuclear norm
rank:         ∑ |λᵢ|⁰ count of non-zero eigenvalues
nuclear norm: ∑ |λᵢ|¹ sum of absolute value of eigenvalues
L₀ norm (sparsity) -> L₁ norm
L₀ norm: ∑|aᵢⱼ|⁰ count of non-zero elements
L₁ norm: ∑|aᵢⱼ|¹ sum of absolute value of elements

L₀ norm regularization。

L₀ norm和L₁ norm效果相同。

abs(x)、log(cosh(x))。

L₁ norm與L₂ norm差異。

備註

古代人的作法
假設 x 的出現機率呈 sigmoid function 或者 laplace distribution
因為獨立  所以所有 x 的出現機率是連乘積
用 maximum likelihood + gradient descent 來解

現代人的做法
把 sigmoid 改成 unit step
機率連乘積  取log變連加  這些過程通通精簡掉  直接變成 1-norm 最佳化
用 least squares (絕對值誤差改成平方誤差以利微分) + gradient descent 來解
感覺上
只要是 sigmoid / tanh
一定可以變成 unit step 之類的東西

如果我沒猜錯
1. 機率學所謂的的獨立事件 是稀疏(1-norm 最佳化)退化到零維的結果
2. 不含機率的稀疏(1-norm最佳化)
   是引入機率的稀疏 sigmoid/tanh/laplace + maximum likelihood 退化之後的結果
http://blog.sciencenet.cn/blog-261330-810156.html
http://www.zhihu.com/question/26602796/answer/33431062
http://blog.csdn.net/abcjennifer/article/details/7748833
PCA  min Frobenius norm >>> min sqrt( trace(AᵀA) )
     min 2-norm of singular values
PCS  min nuclear norm >>> min trace( sqrt(AᵀA) )
     min 1-norm of singular values

Sparse Coding

Sparse Principal Component Analysis

PCA可以改成最佳化問題:令新舊數據的距離平方總和盡量小,限制是新座標軸必須正規正交。

採用regularization,藉由調整參數,以控制距離遠近程度、正規正交程度。regularization所求得的新座標軸通常不是正規正交。然而座標軸歪斜一些也無妨,說不定能讓新舊數據距離平方和更小,找到更擬合數據的座標軸。

             T   2               T
min ‖ X - B̌ B̌ X ‖    subject to B̌ B̌ = I

             T   2            2
min ‖ X - B̌ B̌ X ‖  + α ( ‖ B̌ ‖ - I )     (α ≥ 0)

             T   2          2
min ‖ X - B̌ B̌ X ‖  + α ‖ B̌ ‖             (α ≥ 0)

乍看很理想,不過事情還沒有結束。先前提到,正確答案的左方乘上正規正交矩陣,仍是正確答案。由於答案眾多,於是再增加稀疏性限制:令新座標軸的L¹ norm總和盡量小。

             T   2           2
min ‖ X - B̌ B̌ X ‖₂  + α ‖ B̌ ‖₂  + β ‖ B̌ ‖₁    (α, β ≥ 0)

新座標軸B̌越靠近標準座標軸I,L¹ norm就越小。而B̌是正規正交矩陣,幾何意義是旋轉暨鏡射。也就是說,這是令旋轉角度盡量小、盡量不旋轉數據。

原始論文有提供專門的最佳化演算法,請讀者逕行參閱。

Reconstruction Independent Component Analysis

ICA也可以改成最佳化問題,數據轉置一下,仿效PCA。由於答案很多,於是再增加稀疏性限制:令數據投影到座標軸上的L¹ norm總和盡量小。

       T     T T 2           2         T T
min ‖ X - B̌ B̌ X ‖₂  + α ‖ B̌ ‖₂  + β ‖ B̌ X ‖₁   (α, β ≥ 0)

數據越靠近新座標軸B̌,L¹ norm就越小。也就是說,這是令座標軸盡量貼近數據、盡量讓數據外觀像個十字。

Sparse Coding(Dictionary Learning)

ICA分解版本,不限制正規正交。

演算法是K-SVD,請讀者自行查詢。

http://www.cnblogs.com/salan668/p/3555871.html

Non-negative Matrix Factorization

進一步限制數據非負數、權重非負數,符合真實世界常見情況。

例如圖片處理:數據是灰階值,非負數;權重是圖片合成比重,非負數。

演算法是Projected Gradient Descent,請讀者自行查詢。

https://zhuanlan.zhihu.com/p/22043930

Laplace Matrix(Under Construction!)

章節架構

基礎介紹:各點鄰邊和Lx、所有邊平方和xᵀLx。

推廣版本:邊乘上權重、高維度數據。

相關問題:Dirichlet Energy、Laplace Equation、Laplace Smoothing、Laplace Descent。

延伸應用:形變、平滑、合成、嵌入。

解Lx = b:公式解、不動點遞推法、共軛梯度法、生成樹。

Laplace Operator

「拉普拉斯運算子」。梯度的散度。以下簡稱梯散。

微積分和圖論的定義,相差一個負號。以下採用圖論定義。

微積分的梯散:鄰點減自身,通通加總。
圖論的梯散:自身減鄰點,通通加總。

Laplace Matrix

「拉普拉斯矩陣」。每一點都做拉普拉斯運算,寫成矩陣L。

點相減,做為邊。各點鄰邊和Lx,所有邊平方和xᵀLx。

Laplace Matrix推廣成加權版本

邊乘上權重。A變成權重,D變成權重和。

一點:0
兩點:w₁₂(x₁-x₂)²
三點:w₁₂(x₁-x₂)² + w₁₃(x₁-x₃)² + w₂₃(x₂-x₃)²
L = D - W  (diagonal of D is sum of W)
Lx = [sum Wij (xi - xj)]
xᵀLx = sum Wij (xi - xj)²

加權版本的L,不再是對稱矩陣!

如果你的任務是第一項min xᵀLx,那麼L可以等價地改成對稱矩陣:正向邊與反向邊,兩者權重取平均;原本矩陣與轉置矩陣,相加除以二。二次型函數圖形一模一樣,最佳化結果一模一樣。

如果你的任務是後三項,那麼L不能等價地改成對稱矩陣。

Wsym = (W + Wᵀ) / 2
Lsym = (L + Lᵀ) / 2

Laplace Matrix推廣成向量版本

一維數據們x推廣成多維數據們X。每個維度分開處理。

Lx推廣成LXᵀ。

LXᵀ = 0。梯散為零。求根。

LXᵀ = LYᵀ。保梯散變換。目前大家只求Y,沒人討論變換函數是什麼。

數學特性

  min xᵀLx                   Dirichlet energy
⇒ solve Lx = 0               Laplace equation
⇒ iterate x ≔ D⁻¹Ax         Laplace smoothing
⇔ iterate x ≔ x - D⁻¹Lx     Laplace descent

推導過程

一、當所有邊平方和最小,則各點鄰邊和皆零。

L是對稱正半定矩陣,xᵀLx的函數圖形是開口向上的拋物面,最小值位於一次微分等於零的地方。

min xᵀLx
solve (Lᵀ+L)x = 0   (L is positive semidefinite)
solve 2Lx = 0       (L is symmetric)
solve Lx = 0

二、各點鄰邊和皆零,宛如反覆實施「各點賦值為鄰點平均」。

L是對角優勢矩陣,解Lx = 0得以採用Jacobi Method,每回合恰是「各點賦值為鄰點平均」。

各點鄰點和Ax。各點鄰點平均D⁻¹Ax。

solve Lx = 0 by jacobi method
iterate x ≔ D⁻¹ [0 - (L-D)x]   (L is diagonally dominant matrix)
iterate x ≔ D⁻¹Ax

三、「各點賦值為鄰點平均」可以改寫成梯度下降法的形式。

橫條正規化的拉普拉斯矩陣D⁻¹L。

  D⁻¹Ax
= D⁻¹(D - L)x
= D⁻¹Dx - D⁻¹Lx
= x - D⁻¹Lx

Dirichlet Energy

「狄利克雷能量」。各處梯度長度平方和。等於xᵀLx。

正方形網格,梯度是右邊減自己、上面減自己。若將梯度值填在邊上,那麼各處梯度長度平方和,即是各點上邊右邊平方和,剛好每條邊各算一次,等於xᵀLx。

sum ‖ grad(xi) ‖² = xᵀLx = 1/2 sum sum (xi - xj)²

廣義的「狄利克雷能量」是將正方形網格推廣成任意圖。亦等於xᵀLx。

Laplace Equation

「拉普拉斯方程式」。處處梯散為零:Lx = 0。

解是x = 0或x = 1,缺乏討論意義。大家習慣追加限制條件,例如x部分已知。

經典應用是一坨東西的能量總和盡量小。

Poisson Equation

「帕松方程式」。處處梯散已知:Lx = b。

經典應用是兩坨東西的能量分布盡量相等。其中一坨已知。

圖論風格
min sum sum ‖(xi - xj) - (yi - yj)‖²   y is known
solve Lx = Ly

微積分風格
min ‖grad f(x) - grad g(x)‖²   g is known
min sum ‖f′(xi) - g′(xi)‖²
min sum sum ‖(f(xi) - f(xj)) - (g(xi) - g(xj))‖²
solve laplace f(x) = laplace g(x)

Laplace Smoothing

「拉普拉斯平滑」。每一點設定成鄰點的加權平均。注意到鄰點不含自己。

如果不想一步到位,可以添加一個小於1的倍率,只走一小步。

經典應用是一坨東西的能量總和盡量小。

Laplace Smoothing: Spectral Analysis

http://alice.loria.fr/publications/papers/2010/spectral_course/spectral_course.pdf

走一步的矩陣D⁻¹A,觀察特徵向量、特徵值。然並卵。

Laplace Matrix推廣成加權版本:加權平均座標系統

高維數據視作座標點。觀察「各點賦值為鄰點平均」,平均推廣成加權平均,權重是加權平均座標。

比方來說,鄰點恰好四點、鄰點形成矩形,就可以使用Bilinear Interpolation。

一般來說,鄰點數量不一、鄰點位置隨機,只能採用多點的加權平均座標。

Natural Neighbor Coordinates:
Mean Value Coordinates:遇到凹多邊形,避免權重為負值。
Green Coordinates:

W有很多種選擇。

Wij = 1                                          topology preserving
Wij = 1 / |xi - xj|²                             distance preserving???
Wij = (cot(a) + cot(b)) / 2     ab是邊ij對頂兩角 angle preserving
Wij = (tan(a/2) + tan(b/2)) / 2 ab是邊ij旁邊兩角 mean value coordinate
cotangent lapalcian
1. 三中垂線交於一點
2. 三個等腰三角形
3. alpha = a+b,邊ij對頂角是2a+2b,對半分是a+b = alpha
4. 底(|xi-xj|/2)  高cot(alpha)*(|xi-xj|/2)  相乘除以二就是半個等腰三角形面積
5. 所以總共除以8

tangent lapalcian (mean value coordinate)
(j-1,j,j+1)必須是凸的、alpha+beta < pi,不然會爆掉
折衷方案是theta=(j.i,j+1),把cot(alpha)改成tan(theta/2)
雖然角度不對,但是差不多啦...

面積對點微分  dA/dxi   grad(A)是對邊的垂直方向除以二
http://www.cs.jhu.edu/~misha/Fall09/5-mcf.pdf
http://brickisland.net/cs177/?p=309

如果mesh有boundary,位於boundary上面的點,位於boundary上面的邊,少半個三角形。

Laplace Descent【尚無正式名稱】

「拉普拉斯下降」。仿照梯度下降法,梯度改成了橫條正規化的拉普拉斯矩陣。

我不清楚它的性質。地球人才剛發現這件事:https://arxiv.org/abs/1806.06317

延伸應用

smoothing:數據是座標。能量盡量小。

deformation:數據是座標。先後能量分布盡量相等。

blending:數據是像素值。先後能量分布盡量相等。

embedding:數據是座標。先後能量分布盡量相等。

數學公式(Normal Equation)

追加條件限制之後,方陣變成矩陣,Linear Equation變成Linear Least Squares。

優點是容易追加條件限制,將新的方程式添入矩陣即可。

http://www-ui.is.s.u-tokyo.ac.jp/~takeo/research/rigid/index.html

演算法(Jacobi Method)

Laplace Matrix是對角優勢矩陣,得以使用不動點遞推法。Laplace Matrix構造極簡,不動點遞推法成為最佳首選。

Jacobi Method以舊值算新值,需要兩條陣列。Gauss-Seidel Method立即使用新值,只需要一條陣列。

Symmetric Gauss-Seidel Method修改了計算順序,原本是每回合從頭到尾,改成了頭尾來回。新值連貫,待命時間更少,收斂速度更快。此處的Symmetric,是指頭尾來回,不是指對稱矩陣。

演算法(Conjugate Gradient Method)

Laplace Matrix是對稱正定矩陣,得以使用共軛梯度法。一般矩陣,共軛梯度法才是最佳首選,然而Laplace Matrix構造極簡,導致共軛梯度法失去優勢。

不過還是有人採用Multigrid Preconditioned Conjugate Gradient Method,硬是追加外掛,成為了當前最快的演算法?我沒有仔細讀懂,有勞各位自立自強。

《Efficient Preconditioning of Laplacian Matrices for Computer Graphics》

演算法(KOSZ Algorithm)

Laplace Matrix是對稱對角優勢矩陣,得以使用專門的演算法。雖然時間複雜度極低,但是實際執行時間極高。僅有理論上的價值,沒有實務上的價值。好比Fibonacci Heap。

https://cs.uwaterloo.ca/~lapchi/cs270/notes/L27.pdf
http://appsrv.cse.cuhk.edu.hk/~chi/csc5160/notes/L13.pdf
https://arxiv.org/abs/1502.07888

Laplace Regularization(Under Construction!)

Principal Component Analysis:
Quadratic Form Optimization

xᵀAx是兩兩相乘、通通加總。

PCA可視作二次型最佳化。共變異數矩陣可視作Adjacency Matrix。

找到一個向量v(主成分),讓數據們投影到向量之後,變異數盡量大(盡量散開)。該向量即是最大特徵值的特徵向量。

maximize vᵀAv   A = XXᵀ is symmetric positive-definite

Laplace Equation:
Quadratic Form Optimization

xᵀLx是兩兩相減平方、通通加總。

Laplace Equation可視作二次型最佳化。矩陣是Laplace Matrix。

找到一個純量場x,讓梯度平方總和(Dirichlet Energy)盡量小。該純量場即是零向量,於是大家追加限制條件。

minimize xᵀLx   L is symmetric positive-definite

Laplace Regularization

最佳化的時候,追加目標函數xᵀLx,讓能量盡量小。值得一提的是,當數據是座標點,則是讓幾何形狀盡量平滑。最近幾年似乎很流行。比方說有人用於NMF。