Representation

Representation / Factorization

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

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

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

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

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

Matrix Representation / Matrix Factorization

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

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

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

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

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

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

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

備註

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

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

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

備註

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

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

Error(Loss)

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

Optimization

設定誤差之後,矩陣分解問題變成了最佳化問題!請見本站文件「Matrix Factorization」。

大家習慣令矩陣元素非負,符合實際情況,稱作「非負矩陣分解Nonnegative Matrix Factorization」。

最佳解缺乏實際意義。大家習慣修改誤差函數,獲得其他功效。

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

Principal Coordinates Analysis(Under Construction!)

Principal Coordinates Analysis

「主座標分析」。已知數據Y的兩兩距離M,求得數據Y。

兩兩距離  Mᵢⱼ = ‖pᵢ - pⱼ‖
每項平方  M⊙M = ‖pᵢ - pⱼ‖² = ‖pᵢ‖² + ‖pⱼ‖² - 2‖pᵢ ∙ pⱼ‖
行列中心化 -2YᵀY = C(M⊙M)Cᵀ = -2‖pᵢ ∙ pⱼ‖        C = (I - 𝟏𝟏ᵀ/n)
除以負二  YᵀY = -½C(M⊙M)Cᵀ = ‖pᵢ ∙ pⱼ‖

正確答案是Y = √ΛEᵀ。

Low-rank Principal Coordinates Analysis
(Metric Multidimensional Scaling)

已知數據X,求新數據Y,rank更小,盡量保留距離。

正確答案是XᵀX最大的特徵值暨特徵向量們,然後特徵值開根號。換句話說,X的奇異值暨奇異向量。

主座標分析:√XᵀX
主成分分析:√XXᵀ,然後算√XXᵀ⁻¹X

可以視作Low-rank Matrix Approximation的特例:XᵀX是對稱正定矩陣。

min ‖YᵀY - XᵀX‖²

延伸閱讀:採用Optimization觀點

https://users.ece.cmu.edu/~yuejiec/ece18898G_notes/ece18898g_nonconvex_lowrank_recovery.pdf

「純量對向量微分」請見本站文件「Linear Optimization」。

有人嘗試用「極值位於一次微分等於零的地方」解釋整件事:

min ‖uuᵀ - M‖²
solve 2uᵀ 2(uuᵀ - M) = 0   對u微分
solve uᵀM = uᵀuuᵀ
solve Mᵀu = ‖u‖²u          極值位於Mᵀ的奇異向量

然而上述計算步驟缺乏依據:

一、uᵀu微分推廣成uuᵀ微分。(變數替換v = uᵀ)
二、M從方陣推廣成矩陣。u從特徵向量推廣成奇異向量。
三、向量u推廣成矩陣U。

還必須追加定義,才能用微分解釋整件事。靠你了。

備忘

矩陣降秩應該放到其他網頁。可能要先講PCoA再講PCA。

矩陣降秩,最佳解是奇異值分解。

PCA和PCoA屬於矩陣降秩的特例:距離平方矩陣,經過推導簡化,成為點積矩陣Gram matrix。PCA是矩陣外積、PCoA是矩陣內積。恰是對稱正定矩陣,奇異值改成特徵值而且恆正,奇異向量改成特徵向量而且正規正交。最後從點積矩陣(二次方)還原成原本矩陣(一次方),於是特徵值另外再開根號。

所以說喔,介紹PCA和PCoA的時候,直接給出奇異值分解,其實是抄小路。平方和開根號通通都跳過了,沒有抓到真正精髓。

PCA和PCoA的點積矩陣,屬於二次方,因此可以使用二次型最佳化。矩陣降秩是一次方,因此無法使用二次型最佳化。硬是要用的話,可能必須寫成「先開根號再平方」,然而首先必須定義一下開根號的微分結果等於SVD云云。這個定義想必隱含Eckart-Young的證明原理。很不幸地,直至目前還沒看到有人定義這件事。

二次型最佳化,多向量版本是trace最佳化。PCA和PoCA其實都必須引入trace最佳化,不過這裡省略了。照理來說,trace最佳化必須隱含Eckart-Young的證明原理,但由於PCA、PCoA是正規正交矩陣,大家總是直接寫成subject to再用Lagrange multiplier硬幹,直接聲稱最佳解一定是正規正交、一定是特徵向量,絕口不提Eckart-Young。

trace最佳化當中,矩陣內積XᵀX、矩陣外積XXᵀ、矩陣相乘XY,最佳解可以用柯西不等式來證明:https://math.stackexchange.com/questions/754301/。矩陣相減X-Y,那就只好用Eckart-Young了?

Principal Component Analysis

Principal Component Analysis
採用Optimization觀點

「主成分分析」。大量數據之中,嘗試找到一條過原點直線,數據們投影到直線,數據們到原點的距離平方總和盡量大。這條直線稱做「主成分」,通常只有一條。

寫成數學式子,恰是二次型最佳化。

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ᵀ最大的特徵值的特徵向量。

argmax vᵀXXᵀv         subject to vᵀv = 1    距離平方=點積=矩陣轉置相乘
   ᵛ
argmax vᵀXXᵀv - λ(vᵀv - 1)                  拉格朗日乘數
   ᵛ
solve XXᵀv - λv = 0                         一次微分等於零的地方是極值

solve XXᵀv = λv                             特徵值暨特徵向量

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

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

數據減去平均值

大家總是讓主成分穿過數據中心。預先將數據減去平均值,將數據中心挪至原點;再實施主成分分析,就能讓主成分穿過數據中心。

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

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做。當遞推次數遠小於數據數量,可節省時間。

幾何特性

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

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

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

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

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

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

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

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

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

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

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

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

數學證明

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

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⁻¹。

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

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

正確答案的左邊乘上任意的正規正交矩陣,效果是旋轉與鏡射,仍是正確答案。但由於這類答案並不實用,大家總是無視這類答案。

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

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

一、標準座標軸(單位矩陣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,位移、旋轉、縮放,再乘以√n。將數據範圍大致調整成單位圓,方便比對數據。

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ᵀ
let A = Σ⁻¹Uᵀ, let Y = Vᵀ

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

數學公式(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

Principal Component Analysis
採用統計學觀點

先備知識「Floating Tuple」。

數據們視做一個多變量隨機變數。出現機率均等。

主成分穿越平均值。數據投影到主成分之後,變異數盡量大。

數據們實施仿射變換,每個維度的平均值都變成0,每個維度的變異數都變成1/n,相異維度的共相關數都變成0(不相關)。

YYᵀ = I等同於E[YᵢYᵢ] = 1/n暨E[YᵢYⱼ] = 0。E[YᵢYᵢ] = 1/n是變異數等於1/n(數據中心挪至原點)。E[YᵢYⱼ] = 0是不相關。

正確答案的左邊乘上任意的正規正交矩陣,仍是正確答案。換句話說,每個維度變異數等於1/n、不相關,旋轉之後,仍是如此。

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

矩陣正規化,得到PCA的其中一解!

「共變異矩陣的平方根倒數」,恰是轉換矩陣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。
  此處的定義有別於一般的定義。我查不到專有名詞。

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

Whitening

「白化」。三個願望,一次滿足。

註:三個願望目前沒有廣為人知的正式定義。

centering:中心化。每個維度的平均值都變成零。
standardization:標準化。每個維度的變異數都變成一。(甚至中心化。)
decorrelation:不相關化。相異維度的共變異數(共相關係數)都變成零。

白化可以視作正規化的推廣版本。

normalization:正規化。將每筆數據的範圍,調整成單位長度。
whitening:白化。將多筆數據的範圍,大致調整成單位圓。

白化有許多種演算法,其中一種演算法是PCA:PCA計算結果再乘以√n,讓變異數是1。或者計算時將XXᵀ改成XXᵀ/n。或者描述時將YYᵀ = I改成YYᵀ = nI。

Independent Component Analysis

導讀

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

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

然而,獨立的條件相當嚴苛。一般數據,無法透過仿射變換,成為獨立數據。於是數學家想到兩種折衷方案。

特殊數據

既然一般數據不行,那就討論特殊數據吧:實施了某種仿射變換的獨立數據。

這種特殊數據,可以透過仿射變換(反向變換),成為獨立數據(恢復原貌)。

獨立程度

既然無法完全獨立,那就盡量獨立吧:自行設計指標,評估獨立程度,化作最佳化問題。

常見指標mutual information、cross-cumulant,並未充分運用獨立的特性,既失準、又繁冗。等你發明更好的指標。

Independent Component Analysis

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

原數據是實施了某種仿射變換的獨立數據。

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

先算PCA,再算⟂ICA,即得ICA。

先用PCA扳正獨立成分,再用⟂ICA擺正獨立成分。

一、PCA:得到不相關數據。推定獨立成分互相垂直。
 甲、獨立導致不相關。不相關未必獨立。
 乙、原數據的性質:獨立數據(挪至數據中心時,亦是不相關數據),實施仿射變換。
   PCA的功能:任意數據,實施仿射變換,成為不相關數據。
 丙、原數據實施PCA,成為不相關數據,
   推定為獨立數據,推定獨立成分互相垂直。
   【註:很有可能不是推定、而是確定。我不會證明,也查無證明。】
 丁、然而,PCA所得到的不相關數據,實施任意旋轉,仍是不相關數據。
   必須找到正確方位。
二、⟂ICA:得到獨立數據。得到互相垂直的獨立成分。
 甲、找到最佳向量,做為獨立成分:數據投影到該向量,看起來最獨立。
 乙、自行設計指標,評估獨立程度,化作最佳化問題。
 丙、消滅該維度,遞迴下去,得到互相垂直的獨立成分。
 丁、獨立成分的反矩陣,就是⟂ICA變換矩陣,得到獨立數據。
三、PCA變換矩陣、⟂ICA變換矩陣,依序複合(相乘),即是ICA變換矩陣。
 甲、ICA變換矩陣的反矩陣,推定為原數據的獨立成分。
〇、X減去平均值。
一、PCA(X'):  X' = A₁X    (A₁ = √Λ⁻¹E⁻¹)
二、⟂ICA(X"): Y = A₂X'    (A₂ = F⁻¹)
三、ICA(X):   Y = A₂A₁X   (A = A₂A₁ = F⁻¹√Λ⁻¹E⁻¹)
  A⁻¹ = E√ΛF推定為獨立成分。
  如果討厭鏡射,記得讓E和F成為右手座標系。

據我所知,目前所有的演算法,都是基於此手法。儘管邏輯不夠嚴謹,但是我們也沒有更好的方法了。面對亂七八糟的數據,只好用亂七八糟的方法。

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

PCA:求得不相關數據。

⟂ICA:不相關數據X,個別乘上自身長度X'。共變異矩陣X'X'ᵀ的特徵向量,推定為互相垂直的獨立成分。

共變異矩陣X'X'ᵀ,用來模仿四次動差矩陣X⊗Xᵀ⊗X⊗Xᵀ,用來模仿峰度kurtosis,用來評估獨立程度。

演算法如同兩次PCA。第二次時,數據個別乘上自身長度。

演算法(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
                                             專著當中,此式算錯了,變號了

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

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

消滅該維度:獨立成分垂直投影。沿著先前獨立成分方向。。

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

延伸閱讀:交叉數據

交叉數據,扳正擺正的時候,最接近獨立數據。大家經常拿交叉數據當作ICA的範例。

交叉數據讓人產生誤解,讓人誤以為ICA專門處理交叉數據。其實只是恰好能夠處理而已。

X = [5 1 2 3 4 1 2 3 4 -1 -2 -3 -4 -1 -2 -3 -4; 5 1 2 3 4 -1 -2 -3 -4 1 2 3 4 -1 -2 -3 -4]; X(2, :) = X(2, :) / 2; [d, n] = size(X); X = X - repmat(mean(X,2),1,n); [E, D] = eig(X * X'); Xw = inv(sqrt(D)) * E' * X; [E2, D2] = eig((repmat(sum(Xw.*Xw,1),d,1).*Xw)*Xw'); Y = E2' * Xw; A = E * sqrt(D) * E2; T = inv(sqrt(D)) * E'; scatter(X(1,:), X(2,:)); hold on scatter(E(1,:), E(2,:), 150, 'filled'); hold off axis([-5 5 -5 5]) scatter(Xw(1,:), Xw(2,:));

Laplacian Analysis(Under Construction!)

Laplacian Analysis【尚無正式名稱】
採用Optimization觀點

「梯散分析」。憑空產生大量數據,兩兩距離平方總和盡量小。設計Adjacency Matrix,控制哪些距離要納入計算。如果拿不定主意,那麼可以運用K-nearest Neighbor建立Adjacency Matrix。

為了不讓數據集中於一點,大家習慣追加限制條件。

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

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

距離總和是二次型。PCA與LA是二次型最佳化。

章節架構

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

四項任務:Dirichlet Energy Minimization、Laplace Equation、Laplacian Smoothing、Laplacian Descent。

推廣版本:多維數據、權重、正規化。

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

Laplace Operator

「拉普拉斯運算子」。處處求得梯度的散度。

微積分的梯散:連續版本。鄰點減自身,通通加總。

圖論的梯散:離散版本。自身減鄰點,通通加總。

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

Laplace Matrix

「拉普拉斯矩陣」。拉普拉斯運算子,寫成矩陣L。

拉普拉斯矩陣L=邊數矩陣D-相鄰矩陣A。自身減鄰點,通通加總。

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

數學式子。

L = D - A   let Dᵢᵢ = sumⱼ Aᵢⱼ
Lx = [sum A₀ⱼ (x₀ - xⱼ) , sum A₁ⱼ (x₁ - xⱼ) , ...]ᵀ
xᵀLx = ½ sum sum Aᵢⱼ (xᵢ - xⱼ)²

無向圖是對稱矩陣。有向圖不是對稱矩陣,所有邊平方和不再是xᵀLx,必須另行處置。簡單起見,以下只討論無向圖。

四項任務

主線任務。通通等價。

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

支線任務。通通等價。簡單起見,以下只討論前兩個。

minimize (x-y)ᵀL(x-y)
solve Lx = Ly
iterate x ≔ D⁻¹(Ax + b)         let b = Ly
iterate x ≔ x - D⁻¹(Lx - b)     let b = Ly

推導過程

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

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

L是對稱矩陣,兩個任務等價;L不是對稱矩陣,只能從左到右,無法從右到左。幸好大家習慣討論對稱矩陣。

所有邊平方和xᵀLx,各點鄰邊和Lx。

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

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

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

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

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

三、反覆實施「各點賦值為鄰點平均」⬌梯散下降法。

梯度下降法被改成了梯散下降法。大家正在研究當中。

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

iterate x ≔ D⁻¹Ax
iterate x ≔ D⁻¹(D - L)x
iterate x ≔ D⁻¹Dx - D⁻¹Lx
iterate x ≔ x - D⁻¹Lx

Dirichlet Energy

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

sum ‖ grad xᵢ ‖² = xᵀLx
                 = ½ sum sum Aᵢⱼ (xᵢ - xⱼ)²

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

圖論的狄利克雷能量:任意圖。亦等於xᵀLx。

經典應用是一坨東西盡量平緩。

「對應梯度差異平方和」。尚無正式名稱。等於(x-y)ᵀL(x-y)。

sum ‖ grad xᵢ - grad yᵢ ‖² = (x-y)ᵀL(x-y)
                           = ½ sum sum Aᵢⱼ [(xᵢ - yᵢ) - (xⱼ - yⱼ)]²
                           = ½ sum sum Aᵢⱼ [(xᵢ - xⱼ) - (yᵢ - yⱼ)]² 

若將梯度值填在邊上,那麼對應梯度差異平方和,即是各邊差異平方和,剛好每條邊各算一次,等於(x-y)ᵀL(x-y)。

經典應用是兩坨東西起伏相仿。其中一坨已知。

Laplace Equation / Poisson Equation

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

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

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

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

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

Laplacian Smoothing

「梯散平滑化」。各點賦值為鄰點平均。

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

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

Laplacian Descent【尚無正式名稱】

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

Laplacian Smoothing仿照不動點遞推法,Laplacian Descent仿照梯度下降法。最佳化當中,兩個演算法等價。

Laplace Matrix引入多維數據

一維數據推廣成多維數據,每個維度互不干涉。多組一維數據,各組各自套用Laplace Matrix。

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

solve Lx = 0推廣成solve LXᵀ = 0。

solve Lx = Ly推廣成solve LXᵀ = LYᵀ。

minimize xᵀLx推廣成minimize trace(XLXᵀ)。

minimize (x-y)ᵀL(x-y)推廣成minimize trace((X-Y)L(X-Y)ᵀ)。

Laplace Matrix引入權重

邊乘上權重。各邊權重W。

L = D - W   let Dᵢᵢ = sumⱼ Wᵢⱼ
Lx = [sum W₀ⱼ (x₀ - xⱼ) , sum W₁ⱼ (x₁ - xⱼ) , ...]ᵀ
xᵀLx = ½ sum sum Wᵢⱼ (xᵢ - xⱼ)²

Laplace Matrix引入權重:特殊幾何意義

數據可以是各種東西,例如像素值、點座標。

當數據是點座標,設定特殊權重,XLXᵀ獲得特殊幾何意義。

點對距離總和 Wᵢⱼ = 1
表面積(外心)Wᵢⱼ = (cot(γ₋) + cot(γ₊)) / 4
表面積(內心)Wᵢⱼ = ((tan(α₋/2)⋅tan(β₋/2)) / (tan(α₋/2)+tan(β₋/2)) / 2
           + ((tan(α₊/2)⋅tan(β₊/2)) / (tan(α₊/2)+tan(β₊/2)) / 2
體積     Wᵢⱼ = ?

α₋ = ∠j₋ij, β₋ = ∠j₋ji, γ₋ = ∠ij₋j
α₊ = ∠j₊ij, β₊ = ∠j₊ji, γ₊ = ∠ij₊j

表面積,Adjacency Matrix必須是多面體,每個面都是三角形。

表面積外心版本(中垂線):每個面切成三塊,外心到三頂點。每條邊的左右兩塊三角形相加。缺點:直角導致外心落在邊界,面積為零;鈍角導致外心落在外部,面積為負值。

表面積內心版本(角平分線):同前。缺點:公式太長。

1. 三中垂線交於一點
2. 三個等腰三角形
3. alpha = a+b,邊ij對頂角是2a+2b,對半分是a+b = alpha
4. 底(|xi-xj|/2)  高cot(alpha)*(|xi-xj|/2)
5. 底乘以高除以二就是半個等腰三角形面積
6. 所以總共除以8

體積:我查不到資料。有誰可以幫忙算一下?

minimize trace((X-Y)L(X-Y)ᵀ)使得兩坨東西的對應邊的幾何意義的差異平方和盡量小,宛如處處保X變換。

Laplace Matrix引入權重:加權平均座標

當數據是點座標,觀察Laplacian Smoothing,「各點賦值為『鄰點的加權平均值』」變成了「各點移至『鄰點的加權平均座標』」。

各種加權平均座標都可以當作權重!比方說,表面積計算費時,大家習慣改成Mean Value Coordinate,效果差不多。

Laplace Matrix引入權重:處處梯散為零

存在一種權重,使得處處梯散為零。

性質與用途不明。類似概念出現於Locally Linear Embedding。

Laplace Matrix引入正規化

各邊預先除以鄰點數(引入權重後則是鄰邊總和),就是正規化,可以避免鄰點太多而走遠。W正規化等同L正規化。

橫條正規化(非對稱矩陣)   D⁻¹L
橫條暨直條正規化(對稱矩陣) √D⁻¹L√D⁻¹

延伸應用

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

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

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

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

將任務化作Dirichlet Energy Minimization

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

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

正解缺乏討論意義,於是改為討論歪解。三種經典歪解:

一、消滅維度𝟏:取次小的特徵值的特徵向量。
二、追加限制條件:已知x一部分數值,令其不均等。
三、追加限制條件:平均值是零xᵀ𝟏 = 0、變異數是一‖x‖² = xᵀx = 1。

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

將任務化作Laplace Equation

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

追加限制條件,方陣變成矩陣。

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

數學公式(Inverse)

Laplace Matrix是方陣,但是不存在反矩陣,只存在虛擬反矩陣。

Linear Equation變成Linear Least Squares,數學公式有Normal Equation、QR Decomposition、SVD。

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

演算法(Jacobi Iteration)

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

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

Symmetric Gauss-Seidel Iteration修改了計算順序,原本是每回合從頭到尾,改成了頭尾往返。新值連貫,待命時間更少,收斂速度更快。此處的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/

Laplacian Analysis(Under Construction!)

Laplacian Analysis【尚無正式名稱】
採用Representation觀點

新數據Y處處梯散為零。

然而正解缺乏討論意義:Y的元素通通相等。

LYᵀ = 0
given L. find Y.

改為討論歪解:次小的特徵值的特徵向量們,做為Y的橫條。

因應歪解,solve LYᵀ = 0改寫成minimize trace(YLYᵀ)。

因應特徵向量互相垂直,追加限制條件YYᵀ = I。

{ YYᵀ = I
{ minimize trace(YLYᵀ)
given L. find Y.

Laplacian Principal Component Analysis

數據不是憑空產生,而是已知數據的線性變換。

{ Y = AX                           { Y = AX
{ YYᵀ = I                   --->   { AXXᵀAᵀ = I
{ minimize trace(YLYᵀ)             { minimize trace(AXLXᵀAᵀ)
given X and L. find A and Y.       given X and L. find A and Y.

結果恰是PCA和LA合體!

兩種計算方式,答案好像不太一樣。大家採用第二種。

一、PCA視作主角。取XXᵀ的最大的特徵值的特徵向量們。

正確答案的左邊乘上任意的正規正交矩陣,仍是正確答案。追加LA,從中取得更對味的答案:令所有點對距離盡量均勻。

PCA當中,正確答案們的點對距離完全一致,追加LA沒有任何意義,多此一舉。Low-rank PCA當中,追加LA就有意義。

二、LA視作主角。取XLXᵀ的最小的特徵值的特徵向量們,做為A的橫條。再以A和X求得Y。

因應特徵向量互相垂直,數據必須預先做PCA,使得XXᵀ = I、AXXᵀAᵀ = I、AAᵀ = I。

Laplacian Analysis(Under Construction!)

Correlation

兩個數字a b的共相關數:兩數相乘ab。
一堆數字x₁ x₂ x₃的共相關矩陣:兩兩相乘x₁x₂ x₁x₃ x₂x₃。
一堆高維度數據x₁ x₂ x₃的共相關矩陣:兩兩點積x₁∙x₂ x₁∙x₃ x₂∙x₃。

數據預先減去平均值,就是共變異數、共變異矩陣。

共相關矩陣、共變異矩陣,一定是對稱半正定矩陣。

Principal Component Analysis
採用Correlation觀點

Rank-one Measurement
http://www.princeton.edu/~yc5/ele538b_sparsity/lectures/matrix_recovery.pdf

given yi = (Axi)², solve x
yi = (aiᵀx)² = aiᵀxxᵀai = dot(aiaiᵀ, xxᵀ)

let X = xxᵀ
solve positive semidefinite matrix X
subject to yi = dot(aiaiᵀ,X) and rank(X) = 1
maximize vᵀXXᵀv   subject to ‖v‖² = 1
vᵀSv = vᵀXXᵀv = (d₁∙d₂)(v₁v₂) + (d₁∙d₃)(v₁v₃) + (d₂∙d₃)(v₂v₃)

已知高維數據們X,求得一維數據v(主成分向量)。兩者的維度的共相關矩陣,點積盡量大。

Laplacian Analysis
採用Correlation觀點

minimize xᵀLx   subject to ‖x‖² = 1
xᵀLx = xᵀMMᵀx = (d₁∙d₂)(x₁x₂) + (d₁∙d₃)(x₁x₃) + (d₂∙d₃)(x₂x₃)

已知點邊關係們M,求得一維數據x(各點的權重)。兩者的維度的共相關矩陣,點積盡量小。

PCA是座標,LA是圖。你能銜接兩者嗎?

Kernel Representation(Under Construction!)

Kernel Representation

數據預先實施變換。

Kernel Principal Component Analysis

PCA:

maximize vᵀXXᵀv   subject to vᵀv = 1
maximize vᵀXXᵀv - λ(vᵀv - 1)           Lagrange multiplier
solve XXᵀv = λv                        derivative = 0

數據預先實施變換的PCA:

maximize vᵀΦ(X)Φ(X)ᵀv     st vᵀv = 1   transformation: Φ(X)
maximize vᵀΦ(X)Φ(X)ᵀv - λ(vᵀv - 1)     Lagrange multiplier
solve Φ(X)Φ(X)ᵀv = λv                  derivative = 0
solve Φ(X)Φ(X)ᵀΦ(X)w = λΦ(X)w          kernel trick: v = Φ(X)w
solve Φ(X)ᵀΦ(X)w = λw                  drop Φ(X)
solve K(X)w = λw                       kernel: K(X) = Φ(X)ᵀΦ(X)

變換函數、核函數不是線性函數,不能寫成矩陣乘法ΦX與KX,只能寫成函數求值Φ(X)與K(X)。至於Φ(X)與K(X)都是矩陣。

不設計變換函數Φ(X),設計核函數K(X),精簡計算步驟。

設計核函數K(X),求得最大的特徵值暨特徵向量。w不必還原成v,意思到了就好,精簡計算步驟。

數據減去平均值

中心化可以寫成矩陣C = (I - 𝟏𝟏ᵀ/n),稱作「中心化矩陣」。𝟏是向量,元素皆1。n是數據數量。X減去平均值,可以寫成矩陣相乘CX。

數據實施變換之後,數據中心挪至原點。

maximize vᵀΦ(X)Φ(X)ᵀv     st vᵀv = 1   transformation: Φ(X)
maximize vᵀΦ'(X)Φ'(X)ᵀv   st vᵀv = 1   centering: Φ'(X) = CΦ(X)
maximize vᵀΦ'(X)Φ'(X)ᵀv - λ(vᵀv - 1)   Lagrange multiplier
solve Φ'(X)Φ'(X)ᵀv = λv                derivative = 0
solve Φ'(X)Φ'(X)ᵀΦ'(X)w = λΦ'(X)w      kernel trick: v = Φ'(X)w
solve Φ'(X)ᵀΦ'(X)w = λw                drop Φ'(X)
solve K'(X)w = λw                      kernel: K'(X) = Φ'(X)ᵀΦ'(X)

仍是Kernel PCA,但是很難設計核函數。必須滿足中心化。

有些文獻聲稱:新核函數等於原核函數的橫條直條中心化。我想應該是搞錯了。

K(X) = Φ(X)ᵀΦ(X)
K'(X) = Φ'(X)ᵀΦ'(X) = Φ(X)ᵀCᵀCΦ(X) = Φ(X)ᵀCΦ(X)
K'(X) ≠ CᵀΦ(X)ᵀΦ(X)C = CK(X)C

交換律:核函數、中心化

核函數是距離平方函數,不受位移影響(兩向量同減一向量,距離不變),也不受中心化影響(兩向量同減平均值,距離不變),故有交換律(先中心化再核函數=先核函數再中心化)。

K(CX) = CK(X)
Φ(CX)ᵀΦ(CX) = CΦ(X)ᵀΦ(X)

Kernel Laplacian Analysis【尚無正式名稱】

LA:

minimize vᵀLv   subject to vᵀv = 1
minimize vᵀLv - λ(vᵀv - 1)           Lagrange multiplier
solve Lv = λv                        derivative = 0

核函數硬湊上去:

minimize vᵀΦ(X)Φ(X)ᵀv   st vᵀv = 1   L = Φ(X)Φ(X)ᵀ
solve K(X)w = λw                     kernel: K(X) = Φ(X)ᵀΦ(X)

運用K-nearest Neighbor或其他方法建立Incidence Matrix,做為Φ(X)。

Gram Matrix / Squared Distance Matrix

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

xᵀDx是兩兩一對、各自平方、通通加總。

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

K = Φ(X)ᵀΦ(X)是兩兩點積,Gram Matrix。

K = Φ(X)ᵀΦ(X)習慣設計成兩兩相減平方,Squared Distance Matrix。

距離平方矩陣,行列中心化,得到點積矩陣的-2倍。

CDCᵀ = (I - 𝟏𝟏ᵀ/n)D(I - 𝟏𝟏ᵀ/n)ᵀ = (I - 𝟏𝟏ᵀ/n)D(I - 𝟏𝟏ᵀ/n)

Sparse Representation(Under Construction!)

Sparse Representation

數據維度盡量少,座標數值盡量都是零。距離函數從L₂ norm改成L₀ norm。

Sparse Principal Component Analysis

PCA追加稀疏性。

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

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

min ‖ X - B̌B̌ᵀX ‖²                      subject to B̌ᵀB̌ = I
min ‖ X - B̌B̌ᵀX ‖² + α ( ‖ B̌ ‖² - I )   (α ≥ 0)
min ‖ X - B̌B̌ᵀX ‖² + α ‖ B̌ ‖²           (α ≥ 0)

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

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

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

Laplacian Sparse Coding

Sparse Coding追加Dirichlet Energy Minimization。

Separable Representation(Under Construction!)

Separable Representation

https://arxiv.org/abs/1401.5226

《Algorithmic Aspects of Machine Learning》

Separable Matrix:每個直條j,皆存在橫條i,該橫條僅(i,j)處不為零。

Separable NMF:好像有多項式時間演算法。

數據盡量貼近座標軸、遠離座標軸?ICA和non-gaussianity?

一側為1、一側為0,得到最大割?

max xᵀLx + sum gauss(xi)   subject to ‖x‖² = 1