Representation

Representation / Factorization

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

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

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

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

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

Matrix Representation / Matrix Factorization

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

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

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

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

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

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

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

備註

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

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

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

備註

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

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

Error(Loss)

Y = AX。誤差定為‖Y - AX‖²,矩陣元素平方和。

Optimization

設定誤差之後,矩陣分解問題變成了最佳化問題!

令矩陣元素非負,稱作「非負矩陣分解Nonnegative Matrix Factorization」。導致最佳化過程收斂,導致多項式時間演算法。

矩陣分解的非負限制、一次規劃的第一象限限制、二次規劃的半正定限制,概念皆類似。

非負矩陣分解沒有任何用處,最佳解沒有任何意義。大家習慣修改誤差函數、追加目標函數,獲得其他功效。

Kernel NMF:修改誤差函數,改成kernel,讓數據形狀改變。
Sparse NMF:修改誤差函數,將L₂ norm改成L₀ norm,讓數據形狀盡量稀疏。
Regularized NMF:追加目標函數。例如Laplace Regularization:
                 追加xᵀLx,讓數據形狀盡量平滑。

演算法(Projected Gradient Method)

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

方法一:投影梯度法,max用來投影到限制條件:非負數。
方法二:座標下降法,W和H視作兩個維度。
方法三:修改一下步伐大小、有的沒的,催出速度。

Principal Component Analysis

先備知識

大量數據之中,嘗試找到一條過原點直線,數據們投影到直線,數據們到原點的距離平方總和盡量大。這樣的直線通常只有一條。

寫成數學式子,恰是二次型最佳化。最大值位於XXᵀ最大的特徵值的特徵向量。

argmax ∑ ‖proj pᵢ‖²   subject to ‖v‖² = 1   距離平方總和
   ᵛ   ⁱ    ᵛ
argmax ∑ ‖pᵢ∙v‖²      subject to ‖v‖² = 1   投影=點積(直線向量長度為一)
   ᵛ   ⁱ
argmax ∑ ‖pᵢᵀv‖²      subject to ‖v‖² = 1   點積=矩陣轉置相乘
   ᵛ   ⁱ
argmax ‖Xᵀv‖²         subject to ‖v‖² = 1   數據們合併成一個大矩陣
   ᵛ
argmax (Xᵀv)ᵀ(Xᵀv)    subject to ‖v‖² = 1   距離平方=點積=矩陣轉置相乘
   ᵛ
argmax (vᵀX)(Xᵀv)     subject to ‖v‖² = 1   轉置公式
   ᵛ
argmax vᵀ(XXᵀ)v       subject to ‖v‖² = 1   二次型最佳化
   ᵛ

XXᵀ是維度的兩兩點積,稱作「維度的共相關矩陣」。

數據投影到X軸、Y軸、……,數據的長度的數列,兩兩點積,特徵向量是主成分。數據投影到主成分,數據的長度的數列,平方總和最大。

Principal Component Analysis
採用Optimization觀點

「主成分分析」。大量數據之中,嘗試找到一條直線穿過數據中心,數據們投影到直線,數據們到數據中心的距離平方的平均盡量大(變異數盡量大)。這條直線稱做「主成分」,通常只有一條。

數據中心挪至原點(數據減去平均值),再實施先備知識,即是「主成分分析」。總和換成了平均,不會影響最佳化結果。

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

主成分互相垂直

XXᵀ是對稱半正定矩陣,特徵向量互相垂直,特徵值皆是正數或零。二次型最佳化,駐點位於各個特徵向量,駐點高低順序就是特徵值大小順序。

主成分(極值):最大的特徵值的特徵向量。

主要與次要的主成分(駐點):各個特徵向量,而且互相垂直。

方便起見,一律叫做主成分。D維資料頂多有D個主成分。大家習慣找前幾大的主成分,其特徵值由大到小排序。

主成分可能一樣大,其特徵值一樣大,那麼這些主成分構成的平面、空間、……當中的任何一個向量,都是主成分。然而實務上不會遇到這種情況。

演算法(NIPALS)(Nonlinear Iterative Partial Least Squares)

http://www.statistics4u.com/fundstat_eng/dd_nipals_algo.html

演算法全名軂軂長,原因是作者的思路是線性迴歸、多變量最佳化、反矩陣,不是變異數、二次型最佳化、特徵向量。演算法縮寫多了個A,原因是多了個母音比較好念。

找到主成分、消滅該維度,不斷遞迴下去。

找到主成分:共變異矩陣XXᵀ,最大的特徵值的特徵向量。

消滅該維度:全體數據垂直投影,沿著當前主成分方向。

Power Iteration求特徵向量,Hotelling Deflation做垂直投影。

不預先計算XXᵀ。Power Iteration改成先乘Xᵀ、再乘X。Hotelling Deflation改成對X做。當遞推次數遠小於數據數量,可節省時間。

幾何特性

一、散開性:數據投影到主成分之後,長度平方總和,第一主成分最大,第二主成分次大(須垂直於先前主成分),……。

也就是說,主成分的方向,就是數據最散開的方向。

二、散開性:數據投影到主成分空間之後,長度平方總和最大。

也就是說,主成分所在位置,就是數據最散開的方向。

三、聚合性:數據投影到主成分之時,距離平方總和,第一主成分最小,第二主成分次小(須垂直於先前主成分),……。

也就是說,主成分的垂直方向,就是數據最聚合的方向。

直線擬合必然穿過數據中心,主成分必須穿過數據中心,因此PCA正是直線擬合!

四、聚合性:數據投影到主成分空間之時,距離平方總和最小。

也就是說,主成分所在位置,令數據投影之後的失真最少。

平面擬合必然穿過數據中心,主成分必須穿過數據中心,因此PCA正是平面擬合!

一二與三四是互補的,形成勾股定理。

已知斜邊平方(數據長度平方和),使得一股平方(投影後長度平方和)盡量大、另一股平方(投影時距離平方和)盡量小。

數學證明

散開性。投影之後,長度平方總和最大:

max ∑ ‖proj pᵢ - (1/n) ∑ proj pⱼ‖²      subject to ‖v‖² = 1
 ᵛ  ⁱ    ᵛ             ʲ   ᵛ
max ∑ ‖proj pᵢ‖²                        subject to ‖v‖² = 1
 ᵛ  ⁱ    ᵛ
max ∑ ‖pᵢ∙v‖²                           subject to ‖v‖² = 1
 ᵛ  ⁱ
max ∑ ‖pᵢᵀv‖²                           subject to ‖v‖² = 1
 ᵛ  ⁱ
max ‖Xᵀv‖²                              subject to ‖v‖² = 1
 ᵛ
max (vᵀX)(Xᵀv)                          subject to ‖v‖² = 1
 ᵛ
max vᵀ(XXᵀ)v                            subject to ‖v‖² = 1
 ᵛ
XXᵀ最大的特徵值暨特徵向量就是正解。

聚合性。投影之時,距離平方總和最小:

min ∑ ‖pᵢ - proj pᵢ‖²                   subject to ‖v‖² = 1
 ᵛ  ⁱ         ᵛ
min ∑ ‖pᵢ - (pᵢ∙v) v‖²                  subject to ‖v‖² = 1
 ᵛ  ⁱ
min ∑ ‖pᵢ‖² - 2 ‖pᵢ∙v‖² + ‖pᵢ∙v‖²‖v‖²   subject to ‖v‖² = 1
 ᵛ  ⁱ
min ∑ ‖pᵢ‖² - ‖pᵢ∙v‖²                   subject to ‖v‖² = 1
 ᵛ  ⁱ
min ∑ - ‖pᵢ∙v‖²                         subject to ‖v‖² = 1
 ᵛ  ⁱ
max ∑ ‖pᵢ∙v‖²                           subject to ‖v‖² = 1
 ᵛ  ⁱ
XXᵀ最大的特徵值暨特徵向量就是正解。

精髓

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

Principal Component Analysis

Principal Component Analysis
採用Representation觀點

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ᵀ變成「維度的共變異矩陣」。

幾何意義

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Low-rank Principal Component Analysis

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

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

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

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

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

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

幾何意義

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]]; e1 = InfiniteLine[o, e[[1,All]]]; e2 = InfiniteLine[o, e[[2,All]]]; 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] l = Transpose[{Table[o, Length[p]],q}]; Graphics3D[{Yellow, Specularity[Black, 10], Sphere[q, 0.05], Thickness[0.025], CapForm["Butt"], RGBColor[255,192,0], Line[l], RGBColor[192,80,77], e1, RGBColor[0,176,80], e2, 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]

延伸應用

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

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

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

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

數學公式(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̂ᵀ = [ 3 × 3 ] = [ d₂⋅d₁ d₂⋅d₂ d₂⋅d₃ ] = 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/√λ₃ ] [       ]

數學公式(Singular Value Decomposition)

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

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ᵀ的維度正規正交,即為所求。令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ᵀ

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

數學公式(Normalization)【尚無正式名稱】

矩陣正規化,得到PCA的其中一解!一句話就講完了。

數學與計算學傾向稱作「正規化Normalization」。
統計學傾向稱作「標準化Standardization」或「白化Whitening」。

「共變異矩陣的平方根倒數」,恰是轉換矩陣A乘上正規正交矩陣E,仍是正確答案。

數據X乘上「共變異矩陣的平方根倒數」,恰是PCA乘上正規正交矩陣E,仍是正確答案。

 XXᵀ    = EΛE⁻¹
√XXᵀ    = E√ΛE⁻¹
√XXᵀ⁻¹  = E√Λ⁻¹E⁻¹  = EA  := AwhitenXXᵀ⁻¹X = E√Λ⁻¹E⁻¹X = EAX := AwhitenX

「共變異矩陣的平方根倒數」,宛如除以長度。

數據X乘上「共變異矩陣的平方根倒數」,宛如正規化。

 1       1
———— := ———— = √XXᵀ⁻¹  := Awhiten
‖Xᵀ‖    √XXᵀ
 X       X
———— := ———— = √XXᵀ⁻¹X := AwhitenX
‖Xᵀ‖    √XXᵀ

註:向量x的長度‖x‖定義為√xᵀx。矩陣X的長度‖X‖定義為√XᵀX。
  此處的定義有別於一般的定義。我查不到專有名詞。

矩陣的平方根倒數,目前沒有高速演算法。等你幫忙翻身!

數學公式(Direct Linear Transformation)
演算法(EM-PCA)

E-step:  X = A⁺Y  = (A Aᵀ)⁻¹ Aᵀ Y
M-step:  A = Y X⁺ = X (XᵀX )⁻¹ Xᵀ

一開始建立隨機矩陣A,不斷遞推,直到收斂。

明明就是Block Coordinate Descent,或者說是Fixed Point Iteration,結果卻是EM Algorithm。

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

Principal Component Analysis(Under Construction!)

先備知識

請見本站文件「Floating Tuple」。

數據們視做一個多變量隨機變數。實施線性變換。

獨立數據實施線性變換。

independent  各截面相同(僅倍率不同)
集體位移 ---> 仍然獨立   (改變原點) (原點挪至數據中心)
集體縮放 ---> 仍然獨立
集體旋轉 ---> 失效  
集體歪斜 ---> 失效     (X,X+Y)

不相關數據實施線性變換。

uncorrelated   E[XY] = 0
集體位移 ---> 失效      E[(X+a)(Y+b)] = E[XY] + aE[Y] + bE[X] + ab != 0
集體縮放 ---> 仍然不相關   E[(sX)(tY)] = stE[XY] = 0
集體旋轉 ---> 仍然不相關   E[(RX)(RY)] = RE[XY] = 0
集體歪斜 ---> 失效      E[(X)(X+Y)] = E[XX] + E[XY] != 0

Principal Component Analysis
採用統計學觀點

數據實施線性變換,成為不相關數據。

YYᵀ = I,相異維度的點積等於零,正是E[YᵢYⱼ] = 0。

Independent Component Analysis(Under Construction!)

導讀

PCA專注於不相關。數據實施線性變換,成為不相關數據。

不相關條件寬鬆,任意數據皆可透過線性變換成為不相關數據。

ICA專注於獨立。數據實施線性變換,成為獨立數據。

獨立條件嚴苛,一般數據無法透過線性變換成為獨立數據。

既然不能完全獨立,那就盡量獨立吧。然而目前沒有任何指標可以評估獨立程度,也就沒有演算法了。等你發明。

統計學家改為討論特殊的數據:經過線性變換的獨立數據。如此就不必討論獨立了,有名無實。演算法相當多。

Independent Component Analysis

「獨立成分分析」。數據實施線性變換,成為獨立數據。

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

X是特殊的數據:某種獨立數據,實施某種線性變換。

Y是獨立數據。

新數據的座標軸,經過反向線性變換,得到原數據的「獨立成分」。

計算過程:

一、實施PCA,得到不相關數據。推定獨立成分互相垂直。

二、旋轉數據,將獨立成分擺正。自行定義何謂獨立成分。

一、X實施PCA:A = √Λ⁻¹Eᵀ以及Y = AX。
  獨立成分變得互相垂直。
二、Y實施垂直ICA:A = Vᵀ以及Z = AY。
  消滅該維度,遞迴下去,得到互相垂直的獨立成分。
三、兩個矩陣相乘,即是X的ICA:A = A⟂ICAAPCA以及Z
  反矩陣A⁻¹ = E√ΛV,推定為獨立成分。

如果討厭鏡射,記得讓特徵基底成為右手座標系。

延伸閱讀:交叉數據

交叉數據不是獨立數據的線性變換,不符合獨立成分分析的前提。不知為何,大家習慣拿交叉分布的數據當作範例,甚至大家認為獨立成分分析專門處理交叉數據。積非成是的都市傳說。

聲音訊號不是交叉數據。不知為何,大家習慣將多道聲音訊號畫成交叉數據。真的是非常奇葩。

演算法(FOBI)(Fourth-order Blind Identification)

三步驟:PCA、數據個別乘上自身長度、PCA。

首先求得不相關數據。不相關數據個別乘上自身長度,得到新數據X'。X'X'ᵀ最大的特徵值的特徵向量,推定為獨立成分。

數據預先乘上自身長度的共變異矩陣X'X'ᵀ,用來模仿四次動差矩陣X⊗Xᵀ⊗X⊗Xᵀ,用來模仿峰度kurtosis,用來模仿常態分布不相似度non-gaussianity,用來當作獨立成分。

演算法(FastICA)

作者本人專著《Independent Component Analysis》

f(p) = ‖ E[G(p)] - E[G(q)] ‖²
E[]:期望值。
G():常態分布函數,添加負號。(添加負號是為了稍微加速演算法)
p:給定的一維分布隨機變數。(給定的一維數據們)
q:一維常態分布隨機變數。(無法化作一維數據,無法計算。啊就作者唬爛。)

目標函數‖ E[G(p)] - E[G(q)] ‖²,用來模仿負熵negentropy,用來模仿常態分布不相似度non-gaussianity,用來當作獨立成分。

然而根據推導過程,事實並非如此。目標函數其實是E[G(p)],用來模仿聚散程度,用來當作獨立成分。

數據們靠攏原點,E[G(p)]很小;數據們遠離原點,E[G(p)]很大。其餘性質不明。

推導過程取自專著(8.25)至(8.43),精簡如下:

optimize E[G(vᵀx)]                           數據x投影到獨立成分v:vᵀx
                                             盡量靠攏或遠離原點

v_next = v - F′(v) / F″(v)                   牛頓法求極值
       = v - E[x G′(vᵀx)] / E[xxᵀ G″(vᵀx)]   最小值或最大值

v_next = E[xxᵀ G″(vᵀx)] v - E[x G′(vᵀx)]     右式乘上分母,以消去分母
                                             避免矩陣除法,然而答案不對

v_next = E[G″(vᵀx)] v - E[x G′(vᵀx)]         數據經過PCA處理之後xxᵀ= I
                                             專著當中,此式算錯了,變號了

找到獨立成分、消滅該維度,不斷遞迴下去。

找到獨立成分:牛頓法。目標函數求極值。

消滅該維度:投影牛頓法。投影到先前獨立成分們的垂直方向。

數據實施投影,正確但是慢了點;數據不投影,獨立成分時時投影,不正確但是快了些。當遞推次數遠小於數據數量,可節省時間。

Laplace Matrix(Under Construction!)

導讀

二次型最佳化,請見本站文件「Quadratic Form Optimization」。

PCA專注於數據到數據中心的距離平方總和。為免單調乏味,引入「投影到直線」的想法,操作距離。

Laplace Equation專注於數據到數據的距離平方總和。為免單調乏味,引入「Adjacency Matrix」的想法,操作距離。

距離總和可以視作二次型。PCA與Laplace Equation可以視作二次型最佳化。

章節架構

基礎介紹:各點鄰邊和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。

無向邊,A必須是對稱矩陣。有向邊,A不必是對稱矩陣,但是所有邊平方和不再是xᵀLx,必須另行處置。簡單起見,本篇文章只討論無向邊。

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

Laplace Matrix推廣成向量版本

一維數據們x推廣成多維數據們X。

每個維度視作一維數據,分開套用Laplace Matrix。

Lx推廣成LXᵀ。xᵀLx推廣成XLXᵀ

四項任務

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

推導過程

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

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

L是對稱矩陣,兩個任務等價;L不是對稱矩陣,只能從左到右,無法從右到左。

所有邊平方和xᵀLx,各點鄰邊和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 / Poisson Equation

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

經典應用是一坨東西的勢力均衡。

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

常見形式是兩邊梯散相同Lx = Ly,y已知,或者Ly已知。

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

圖論風格
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,走小步。

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

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

經典應用是讓一坨東西變得均勻。

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【尚無正式名稱】

「拉普拉斯下降法」。仿照梯度下降法,梯度改成了橫條正規化的拉普拉斯矩陣。地球人才剛發現,正體不明。

《Laplacian Smoothing Gradient Descent》

延伸應用

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

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

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

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

將任務化作Dirichlet Energy

這四項任務,大家習慣化作第一項任務min xᵀLx。

優點:容易追加最佳化目標、追加限制條件。缺點:L是對稱矩陣才能等價化作第一項任務;幸好大家習慣討論對稱矩陣。

最小值是0,位於最小的特徵值0的特徵向量1。換句話說,正解是x = 1k,1是向量[1,1,1,...]ᵀ,k是任意倍率。換句話說,正解是各點數值均等。

正解缺乏討論意義,於是改為討論歪解。兩種做法:

一、追加限制條件:已知x一部分數值,令其不均等。

二、消滅維度x = 1k:取次小的特徵值的特徵向量。

將任務化作Laplace Equation

這四項任務,亦得化作第二項任務solve Lx = 0。

優點:Lx = b求解,擁有高速演算法。以下介紹這些演算法。

數學公式(Inverse Matrix)

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

求反矩陣。Normal Equation、QR、SVD。

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

演算法(Jacobi Method)

https://math.stackexchange.com/questions/1727204/

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
https://sites.google.com/a/yale.edu/laplacian/

Kernel Representation(Under Construction!)

Kernel Representation

數據預先實施變換。Kernel是數據變換之後的距離函數。

Kernel Principal Component Analysis

maximize vᵀXXᵀv   subject to vᵀv = 1
maximize vᵀXXᵀv - λ(vᵀv - 1)              Lagrange multiplier
solve XXᵀv = λv                           derivative = 0
solve Φ(X)Φ(X)ᵀw = λw                     transformation
solve Φ(X)Φ(X)ᵀΦ(X)α = λΦ(X)α             kernel trick
solve Φ(X)ᵀΦ(X)Φ(X)ᵀΦ(X)α = λΦ(X)ᵀΦ(X)α
solve K²α = λKα   let K = Φ(X)ᵀΦ(X)
solve Kα = λα
solve K'α = λα    let K'= (I - eeᵀ/n)K(I - eeᵀ/n)  centering

Kernel Laplace Smoothing

why heat kernel? because d/dt = laplace(x) is heat kernel?

Gram Matrix / Squared Distance Matrix

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

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

Isomap / Locally Linear Embedding

C改成A?A改成D?

Correlation Matrix

共相關數矩陣C。

兩個數字a b的共相關數是兩數相乘ab。

一堆數字x1 x2 x3的共相關數們是兩兩相乘x1x2, x1x3, x2x3。

一堆高維度數據x1 x2 x3的共相關數們是兩兩點積x1∙x2, x1∙x3, x2∙x3。

如果數字預先減去平均值,就是共變異數。

延伸閱讀:Correlation Matrix與Adjacency Matrix不相容

共相關矩陣換成相鄰矩陣會發生什麼呢?A用來設定哪些數據配對互相相關。

通常無解。比方說,三點兩邊:兩邊端點之共相關數是1,缺邊端點之共相關數是0,各點找不到合適的數值。ab = ac = 1 and bc = 0無解。

Principal Component Analysis

maximize xᵀCx   subject to ‖x‖² = 1
maximize vᵀCv   C = XXᵀ is symmetric positive-definite

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

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

xᵀAx是兩兩相乘、通通加總。因此yᵀXXᵀy = (d1∙d2)(y1y2) + (d1∙d3)(y1y3) + (d2∙d3)(y2y3),高維數據們X、主成分向量y,維度的共相關數矩陣的點積盡量大!

給定高維數據,找到一維數據,共相關數越像越好。

找到一維數據,維度的共相關數總和越大越好(不是平方和)。

Sparse Representation(Under Construction!)

Sparse Representation

數據盡量貼近座標軸。

Convex Relaxation

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

     min ‖ L ‖₀ + α ‖ S ‖₀     (α ≥ 0)
---> min ‖ L ‖⁎ + α ‖ S ‖₁     (α ≥ 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差異。

Basis Pursuit

Sparse Principal Component Analysis

PCA追加稀疏性。其實應該叫做ICA吧。

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̌是正規正交矩陣,幾何意義是旋轉暨鏡射。也就是說,這是令旋轉角度盡量小、盡量不旋轉數據。

Sparse Coding(Dictionary Learning)

NMF追加稀疏性。不限制正規正交。

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

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

http://blog.sciencenet.cn/blog-261330-810156.html
http://www.zhihu.com/question/26602796/answer/33431062
http://blog.csdn.net/abcjennifer/article/details/7748833

compressed sensing = underdetermined system sparse coding

https://dsp.stackexchange.com/questions/44257/
https://zhuanlan.zhihu.com/p/22445302
https://www.cnblogs.com/AndyJee/category/579870.html

Principal Component Pursuit

X = L + S。

low-rank L。令秩盡量小。

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

sparse S。令矩陣盡量稀疏。

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

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

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

Robust Principal Component Analysis

X = L + S + E。

https://users.ece.cmu.edu/~yuejiec/teaching.html
https://sites.google.com/site/thierrybouwmans/
https://sites.google.com/site/robustdlam/

備註

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

備註

古代人的作法
假設 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 退化之後的結果