factoring
factoring【尚無正式名稱】
「分解」。擷取數據的特質,進而調整數據。
註:以下介紹的主題,屬於多變量統計分析multivariate statistics analysis領域的其中一部分。
principal coordinate analysis
Gram matrix
「點積矩陣」。儲存兩兩點積,元素Aᵢⱼ是i點與j點的點積。
⎡ p₀∙p₀ p₀∙p₁ p₀∙p₂ p₀∙p₃ ... ⎤ ⎢ p₁∙p₀ p₁∙p₁ p₁∙p₂ p₁∙p₃ ... ⎥ G = ⎢ p₂∙p₀ p₂∙p₁ p₂∙p₂ p₂∙p₃ ... ⎥ ⎢ p₃∙p₀ p₃∙p₁ p₃∙p₂ p₃∙p₃ ... ⎥ ⎣ : : : : ⎦
給定數據矩陣,求得點積矩陣:對稱半正定矩陣XᵀX。
⎡ | | ⎤ ⎡ ⎤⎡ | ⎤ ⎡ : ⎤ X = ⎢ p₀ p₁ .. ⎥ XᵀX = ⎢ —— pᵢ —— ⎥⎢ pⱼ ⎥ = ⎢.. pᵢ∙pⱼ ..⎥ ⎣ | | ⎦ ⎣ ⎦⎣ | ⎦ ⎣ : ⎦
給定點積矩陣,求得數據矩陣:答案無限多個,例如特徵分解的共軛分解X = √ΛEᵀ、Cholesky分解的共軛分解X = Lᵀ。
eigendecomposition: Cholesky decomposition: XᵀX = EΛEᵀ XᵀX = LLᵀ XᵀX = E√Λ√ΛEᵀ = (√ΛEᵀ)ᵀ√ΛEᵀ X = Lᵀ X = ⟌XᵀX = √ΛEᵀ
註:共軛分解的數學符號⟌,我自己瞎掰的。
distance matrix
「距離矩陣」。儲存兩兩距離,元素Aᵢⱼ是i點與j點的距離。
⎡ ‖p₀-p₀‖ ‖p₀-p₁‖ ‖p₀-p₂‖ ‖p₀-p₃‖ ... ⎤ ⎢ ‖p₁-p₀‖ ‖p₁-p₁‖ ‖p₁-p₂‖ ‖p₁-p₃‖ ... ⎥ M = ⎢ ‖p₂-p₀‖ ‖p₂-p₁‖ ‖p₂-p₂‖ ‖p₂-p₃‖ ... ⎥ ⎢ ‖p₃-p₀‖ ‖p₃-p₁‖ ‖p₃-p₂‖ ‖p₃-p₃‖ ... ⎥ ⎣ : : : : ⎦
給定數據矩陣,求得距離矩陣:利用全一矩陣。
每項平方 X⊙X 橫條平方和 (X⊙X)𝟏 = ‖pᵢ‖² 直條平方和 𝟏(X⊙X) = ‖pⱼ‖² 兩兩點積 XᵀX = (pᵢ ∙ pⱼ) 兩兩距離平方 (X⊙X)𝟏 + 𝟏(X⊙X) - 2(XᵀX) = M⊙M ‖pᵢ‖² + ‖pⱼ‖² - 2(pᵢ ∙ pⱼ) = ‖pᵢ - pⱼ‖² 每項開根號 √[(X⊙X)𝟏 + 𝟏(X⊙X) - 2(XᵀX)]ᵢⱼ = Mᵢⱼ
給定距離矩陣,求得數據矩陣:利用中心化矩陣。
兩兩距離 Mᵢⱼ = ‖pᵢ - pⱼ‖ 每項平方 M⊙M = ‖pᵢ - pⱼ‖² = ‖pᵢ‖² + ‖pⱼ‖² - 2(pᵢ ∙ pⱼ) 行列中心化 C(M⊙M)C = -2(pᵢ ∙ pⱼ) = -2XᵀX 除以負二 -½C(M⊙M)C = (pᵢ ∙ pⱼ) = XᵀX 共軛分解 ⟌-½C(M⊙M)C = X
證明:行列中心化
距離平方矩陣Mᵢⱼ² = ‖pᵢ‖² + ‖pⱼ‖² - 2(pᵢ ∙ pⱼ)。
試證:橫條中心化、直條中心化,可以消去平方項‖pᵢ‖²和‖pⱼ‖²,留下交叉項-2(pᵢ ∙ pⱼ)。
一、中心化可以消去平方項。
橫條擁有相同‖pᵢ‖²。各個橫條各自減去平均數,可以消去‖pᵢ‖²。橫條中心化消去橫條平方項。同理,直條中心化消去直條平方項。
⎡ ‖p₀‖²+‖p₀‖²-2(p₀∙p₀) ‖p₀‖²+‖p₁‖²-2(p₀∙p₁) ... ⎤ ⎢ ‖p₁‖²+‖p₀‖²-2(p₁∙p₀) ‖p₁‖²+‖p₁‖²-2(p₁∙p₁) ... ⎥ M⊙M = ⎢ ‖p₂‖²+‖p₀‖²-2(p₂∙p₀) ‖p₂‖²+‖p₁‖²-2(p₂∙p₁) ... ⎥ ⎢ ‖p₃‖²+‖p₀‖²-2(p₃∙p₀) ‖p₃‖²+‖p₁‖²-2(p₃∙p₁) ... ⎥ ⎣ : : ⎦ ^^^^^ ^^^^^
二、中心化也會改變交叉項。
但是我們需要交叉項。解決方法:假設數據已經中心化。
此時數據任一維度總和為零∑pⱼ = 0,使得橫條交叉項總和為零∑ⱼ(-2(pᵢ ∙ pⱼ)) = -2(pᵢ ∙ (∑pⱼ)) = 0。橫條中心化不改變橫條交叉項。同理,直條中心化不改變直條交叉項。
⎡ ‖p₀‖²+‖p₀‖²-2(p₀∙p₀) ‖p₀‖²+‖p₁‖²-2(p₀∙p₁) ... ⎤ ⎢ ‖p₁‖²+‖p₀‖²-2(p₁∙p₀) ‖p₁‖²+‖p₁‖²-2(p₁∙p₁) ... ⎥ M⊙M = ⎢ ‖p₂‖²+‖p₀‖²-2(p₂∙p₀) ‖p₂‖²+‖p₁‖²-2(p₂∙p₁) ... ⎥ ⎢ ‖p₃‖²+‖p₀‖²-2(p₃∙p₀) ‖p₃‖²+‖p₁‖²-2(p₃∙p₁) ... ⎥ ⎣ : : ⎦ ^^^^^^^^^ ^^^^^^^^^
三、橫條中心化、直條中心化,無論誰先做,結果都一樣。
principal coordinate analysis
「主座標分析」。已知距離矩陣,求得數據矩陣。
方才已介紹。換句話說吧。
一、距離矩陣M,每項平方,求得距離平方矩陣M⊙M。 二、距離平方矩陣M⊙M,橫條中心化,直條中心化,除以-2,求得點積矩陣XᵀX。 三、點積矩陣XᵀX,共軛分解,求得數據矩陣X。
low-rank principal coordinate analysis
「低秩主座標分析」。已知距離矩陣,求得數據矩陣。降低數據維度,盡量保持原本距離。
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}}; a = Flatten[Table[Table[p[[i]], Length[p]], {i, Length[p]}], 1]; b = Flatten[Table[p, Length[p]], 1]; l = Transpose[{a,b}]; Graphics3D[{Black, Specularity[White, 10], Opacity[0.2], Sphere[p, 0.05], Thickness[0.01], CapForm["Butt"], RGBColor[255,192,0], Opacity[0.9], Line[l]}, PlotRange -> {{-.6,.6},{-.6,.6},{-.6,.6}}, Boxed -> False] o = Mean[p]; p = p - Table[o, Length[p]]; G = p . Transpose[p]; e = Eigenvectors[N[G], 2]; v = Eigenvalues[N[G], 2] q = Transpose[DiagonalMatrix[v] . e]; Graphics[{Black, PointSize[0.05], Point[q]}, PlotRange -> {{-.6,.6},{-.6,.6},{-.6,.6}}, Boxed -> False]
修改最後一步。點積矩陣,保留前k大的特徵值們,其餘歸零。特徵分解的右半邊√ΛEᵀ,即是k維數據矩陣。
「low-rank matrix approximation」的最佳解是奇異值分解。點積矩陣是對稱半正定矩陣,於是奇異值分解等同特徵分解。
降維對象是點積矩陣。為何不是距離矩陣呢?這是一個謎,等你發表論文。先中心化再低秩,先低秩再中心化,兩者結果應該不同。
multidimensional scaling
「多維度縮放」。已知原始數據,求得降維數據,盡量保持原本距離。
求解步驟較少,不必求得距離矩陣、點積矩陣。
數據矩陣中心化。數據矩陣,保留前k大的奇異值們,其餘歸零。奇異值分解的右半邊ΣVᵀ,即是k維數據矩陣。
XᵀX即是點積矩陣。XᵀX的特徵向量,即是X的右奇異向量。XᵀX的特徵值開根號,即是X的奇異值。√ΛEᵀ = ΣVᵀ。
graph Laplacian analysis
楔子
PCoA:已知兩兩距離,求得數據。
GLA:盡量縮短兩兩距離,求得數據。
graph Laplacian analysis【尚無正式名稱】
「圖梯散分析」。製造大量數據,兩兩距離平方總和盡量小。
設計adjacency matrix,控制哪些距離需要納入計算。為了不讓數據集中於一點,必須追加限制條件,例如自訂部分數據座標。
adjacency matrix
「相鄰矩陣」。請見本站文件「adjacency matrix」。
兩點有關,邊為1。兩點無關,邊為0。
相鄰矩陣是對稱矩陣。
兩兩距離平方總和
兩兩距離平方總和。
sum ‖pᵢ - pⱼ‖² (i,j)
兩兩距離平方總和盡量小。
min sum ‖pᵢ - pⱼ‖² P (i,j)
引入adjacency matrix,只考慮部分的兩兩距離。
min sum Aᵢⱼ ‖pᵢ - pⱼ‖² P (i,j)
一個距離平方,可以拆解成各維距離平方。各維各自統計距離平方,各維各自最佳化,結果仍相同。多維問題等價於多個一維問題。
min sum Aᵢⱼ ‖xᵢ - xⱼ‖² + min sum Aᵢⱼ ‖yᵢ - yⱼ‖² x (i,j) y (i,j)
Laplacian matrix
「梯散矩陣」。請見本站文件「Laplacian matrix」。
自身減鄰點,通通加總。
點相減,作為邊。各點鄰邊和Lx。所有邊平方和xᵀLx。
梯散矩陣L=邊數矩陣D-相鄰矩陣A。
梯散矩陣是對稱半正定矩陣。
一維兩兩距離平方總和
一維兩兩距離平方總和,可以改寫成二次型xᵀLx。
點的數值,設定為座標。
sum Aᵢⱼ ‖xᵢ - xⱼ‖² = xᵀLx (i,j)
一維兩兩距離平方總和盡量小。
L是對稱半正定矩陣,xᵀLx的函數圖形是開口向上的拋物面,擁有最小值。
min sum Aᵢⱼ ‖xᵢ - xⱼ‖² = min xᵀLx x (i,j) x
高維兩兩距離平方總和
計算過程習慣分解成一維版本。數學式子習慣寫成高維版本。
高維兩兩距離平方總和,可以改寫成矩陣運算。
1D | KD ---------------|--------------------- x | X Lx | LXᵀ xᵀLx | XLXᵀ solve Lx = 0 | solve LXᵀ = 0 minimize xᵀLx | minimize tr(XLXᵀ)
演算法
一維兩兩距離平方總和盡量小。可以改寫成二次型最佳化。
該二次型稱為Dirichlet energy。
minimize xᵀLx
極值位於一次微分等於零的地方。可以改寫成一次方程組求解。
該一次方程組稱為Laplace equation。
solve Lx = 0
極值位於最小的特徵值的特徵向量。可以改寫成特徵分解。
其特徵值稱為Laplacian spectrum。
L = EΛEᵀ
演算法共三類:二次型最佳化、一次方程組求解、特徵分解。
演算法:二次型最佳化
目標函數:一維兩兩距離平方總和盡量小。
minimize f(x) where f(x) = sum Aᵢⱼ ‖xᵢ - xⱼ‖² = xᵀLx
追加限制條件:自訂部分數據座標。
minimize f(x) subject to x₀ = c₀, x₁ = c₁, ...
或者也可以追加函數。
minimize f(x) + λ g(x)
目標函數,分別對每筆數據偏微分,得到梯度。
⎡ ∂/∂x₀ sum Aᵢⱼ ‖xᵢ - xⱼ‖² ⎤ f′(x) = ⎢ ∂/∂x₁ sum Aᵢⱼ ‖xᵢ - xⱼ‖² ⎥ ⎣ : ⎦
目標函數,改寫成二次型,對向量微分,得到梯度。結果一樣。
f′(x) = (L + Lᵀ)x = 2Lx
演算法一覽。
最佳化:梯度下降法。 二次型最佳化:共軛梯度法。
演算法(gradient descent)
「梯度下降法」。請見本站文件「地面勘查類型」。
自訂初始位置。自訂步伐大小。係數2可以併入步伐大小。
xnext = x - step ⋅ 2Lx
改寫成不動點遞推法,以判斷是否收斂至最佳解。
xnext = (I - step ⋅ 2L)x
視作特徵向量演算法power iteration。
當step足夠小,則(I - step ⋅ 2L)是對稱正定矩陣,則收斂。
對稱正定矩陣,存在N個實數特徵向量,足以支配整個空間。任意初始向量(除了零向量)必定趨近最大的特徵值的特徵向量。
gradient descent:各座標同時移動。 coordinate descent:各座標輪流移動。
演算法(conjugate gradient method)
「共軛梯度法」。請見本站文件「quadratic optimization」。
實務上最快的演算法。N步收斂至最佳解。調整誤差容許範圍,提早結束演算法。
甚至追加外掛,採用multigrid preconditioned conjugate gradient method。我沒有仔細讀懂,有勞各位自立自強。
《Efficient Preconditioning of Laplacian Matrices for Computer Graphics》
principal component analysis
principal component analysis
「主成分分析」。大量數據之中,嘗試找到一條過原點直線,數據們投影到直線,數據們到原點的距離平方總和盡量大。這條直線稱作「主成分」,通常只有一條。
寫成數學式子,恰是二次型最佳化。
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ᵀ是對稱半正定矩陣,特徵向量互相垂直,特徵值皆是正數或零。二次型最佳化,駐點位於各個特徵向量,駐點高低順序就是特徵值大小順序。
主成分(極值):最大的特徵值的特徵向量。
主要與次要的主成分(駐點):各個特徵向量,而且互相垂直。
方便起見,一律叫做主成分。D維資料頂多有D個主成分。大家習慣找前幾大的主成分,其特徵值由大到小排序。
主成分可能一樣大,其特徵值一樣大,那麼這些主成分構成的平面、空間、……當中的任何一個向量,都是主成分。然而實務上不會遇到這種情況。
幾何特性
一、散開性:數據投影到主成分之後,長度平方總和,第一主成分最大,第二主成分次大(須垂直於先前主成分),……。
也就是說,主成分的方向,就是數據最散開的方向。
二、散開性:數據投影到主成分空間之後,長度平方總和最大。
也就是說,主成分所在位置,就是數據最散開的方向。
三、聚合性:數據投影到主成分之時,距離平方總和,第一主成分最小,第二主成分次小(須垂直於先前主成分),……。
也就是說,主成分的垂直方向,就是數據最聚合的方向。
直線擬合必然穿過數據中心,主成分必須穿過數據中心,因此主成分正是直線擬合!
四、聚合性:數據投影到主成分空間之時,距離平方總和最小。
也就是說,主成分所在位置,令數據投影之後的失真最少。
平面擬合必然穿過數據中心,主成分必須穿過數據中心,因此主成分正是平面擬合!
一二與三四是互補的,形成勾股定理。
已知斜邊平方(數據長度平方和),使得一股平方(投影後長度平方和)盡量大、另一股平方(投影時距離平方和)盡量小。
數學證明
散開性。投影之後,長度平方總和最大:
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ᵀ最大的特徵值暨特徵向量就是正解。
精髓
已知數據投影到X軸、Y軸、……的長度,求得數據投影到主成分的長度。長度平方總和盡量大。
從不同角度看數據,時聚時散。不斷滾動數據,重新找一個視角,讓數據看起來最散開。注意到,是聚散,不是寬窄。
演算法(NIPALS)(nonlinear iterative partial least squares)
演算法全名軂軂長,原因是作者的思路是一次迴歸、多變量最佳化、反矩陣,不是變異數、二次型最佳化、特徵向量。演算法縮寫多了個A,原因是多了個母音比較好念。
XXᵀ的特徵向量演算法。
找到主成分、消滅該維度,不斷遞迴下去。
找到主成分:共變異矩陣XXᵀ,最大的特徵值的特徵向量。
消滅該維度:全體數據垂直投影,沿著當前主成分方向。
power iteration求特徵向量,Hotelling deflation做垂直投影。
不預先計算XXᵀ。power iteration改成先乘Xᵀ、再乘X。Hotelling deflation改成對X做。當遞推次數遠小於數據數量,可節省時間。
principal component analysis
correlation matrix
「共相關矩陣」。儲存兩兩共相關數,元素Aᵢⱼ是i物件與j物件的共相關數。兩個向量的共相關數通常定義為點積。
⎡ x₀∙y₀ x₀∙y₁ x₀∙y₂ x₀∙y₃ ... ⎤ ⎢ x₁∙y₀ x₁∙y₁ x₁∙y₂ x₁∙y₃ ... ⎥ RXY = ⎢ x₂∙y₀ x₂∙y₁ x₂∙y₂ x₂∙y₃ ... ⎥ ⎢ x₃∙y₀ x₃∙y₁ x₃∙y₂ x₃∙y₃ ... ⎥ ⎣ : : : : ⎦
共相關矩陣有兩種版本:數據兩兩點積XᵀY、維度兩兩點積XYᵀ。大家習慣採用維度兩兩點積XYᵀ。數據數量相等,方可點積。
⎡ | | ⎤ ⎡ ⎤⎡ | ⎤ ⎡ : ⎤ X = ⎢ p₀ p₁ .. ⎥ XᵀY = ⎢ —— pᵢ —— ⎥⎢ qⱼ ⎥ = ⎢.. pᵢ∙qⱼ ..⎥ ⎣ | | ⎦ ⎣ ⎦⎣ | ⎦ ⎣ : ⎦ ⎡ —— d₀ —— ⎤ ⎡ ⎤⎡ | ⎤ ⎡ : ⎤ X = ⎢ —— d₁ —— ⎥ XYᵀ = ⎢ —— dᵢ —— ⎥⎢ eⱼ ⎥ = ⎢.. dᵢ∙eⱼ ..⎥ ⎣ : ⎦ ⎣ ⎦⎣ | ⎦ ⎣ : ⎦
共相關矩陣細分為互相關矩陣XYᵀ、自相關矩陣XXᵀ。
cross-correlation matrix:互相關矩陣。相異數據集的相關矩陣。XYᵀ。 autocorrelation matrix:自相關矩陣。相同數據集的相關矩陣。XXᵀ。
預先將數據中心化,則形成互變異矩陣(X-X̄)(Y-Ȳ)ᵀ、自變異矩陣(X-X̄)(X-X̄)ᵀ。
cross-covariance matrix:互變異矩陣。(X-X̄)(Y-Ȳ)ᵀ。 autocovariance matrix:自變異矩陣。(X-X̄)(X-X̄)ᵀ。
partial least squares regression
「偏最小平方迴歸」。兩堆數據的共相關數。
大家習慣預先將數據中心化。互變異矩陣的奇異值分解。
兩堆數據之中,分別找到一條直線穿過數據中心,數據們投影到直線,數據們到原點的距離,兩邊點積盡量大。
正確答案:最大的奇異值的左奇異向量,最大的奇異值的右奇異向量。點積結果是最大的奇異值,作為共相關數。
數據編號一旦改變,計算結果跟著改變。「偏最小平方迴歸」並不穩重成熟。
principal component analysis
「主成分分析」。兩堆相同數據的共相關數。
大家習慣預先將數據中心化。自變異矩陣的奇異值分解。
自變異矩陣是對稱半正定矩陣,奇異值分解等同特徵分解。
點積結果是最大的特徵值,即是距離平方總和,作為共相關數。
數據編號即使改變,計算結果依然不變!
精髓
GLA:已知點邊關係們M,求得一維數據x(點座標),讓共相關矩陣MMᵀ和xxᵀ的點積盡量小。
minimize xᵀLx subject to ‖x‖² = 1
xᵀLx = xᵀMMᵀx = (MMᵀ)∙(xxᵀ) = (d₁∙d₁)(x₁x₁) + (d₁∙d₂)(x₁x₂) + (d₁∙d₃)(x₁x₃) + ...
PCA:已知高維數據們X,求得一維向量v(主成分),讓共相關矩陣XXᵀ和vvᵀ的點積盡量大。
maximize vᵀXXᵀv subject to ‖v‖² = 1
vᵀCv = vᵀXXᵀv = (XXᵀ)∙(vvᵀ) = (d₁∙d₁)(v₁v₁) + (d₁∙d₂)(v₁v₂) + (d₁∙d₃)(v₁v₃) + ...
principal component analysis
楔子
GLA:數據到數據的距離平方總和盡量小。為免單調乏味,引入「adjacency matrix」以操作距離。最後將數據扯緊調勻。
PCA:數據到數據中心的距離平方總和盡量大。為免單調乏味,引入「投影到直線」以操作距離。最後將數據搓圓擺正。
principal component analysis
Y = AX。已知數據X,求新數據Y、求轉換矩陣A。
YYᵀ = I。令新數據Y的維度是正規正交:相同維度的點積是1,相異維度的點積是0。既是單位向量、又互相垂直。
一個維度視作一個向量;兩個向量的點積,得到一個值;所有兩兩向量的點積,排列成矩陣。
數據減去平均數
當X的中心(平均數)位於原點,則Y的中心也將位於原點。
大家總是預先將X減去平均數,將X的中心挪至原點。好處是:一、降低數值範圍,以減少浮點數運算誤差。二、具備直線擬合效果。三、這樣才是標準的主成分分析。
推導過程
⎰ 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 | XXᵀ has real non-negative eigenvalues. | A(E√Λ)(√ΛE⁻¹)Aᵀ = I | 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。
幾何意義
PCA找到了全新的垂直座標軸:原點是數據中心、座標軸是特徵基底、座標軸單位長度是特徵值平方根。
所有數據一齊位移、旋轉(與鏡射)、縮放(與鏡射)。
一、位移:數據中心位移至原點。
二、旋轉:以特徵基底,逆向旋轉數據。E⁻¹。
三、縮放:各維度除以各特徵值平方根。√Λ⁻¹。
正確答案不只一種。特徵向量(連同特徵值)對調次序,效果是鏡射暨旋轉。特徵向量顛倒方向、特徵值變號,效果是鏡射。
正確答案的左邊乘上任意的正規正交矩陣,效果是旋轉與鏡射,仍是正確答案。但由於這類答案並不實用,大家總是無視這類答案。
low-rank principal component 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]]; 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。或者計算時將XXᵀ改成XXᵀ/n。或者描述時將YYᵀ = I改成YYᵀ = nI。
orientation:實施PCA,位移、旋轉、略過縮放。功效是重設座標軸,擺正數據。
alignment:兩堆數據各自實施PCA,特徵值由大到小排序,特徵向量一一對應,得知數據對應關係。
embedding:實施low-rank PCA,捨棄小的特徵值暨特徵向量,只做投影。功效是降維、壓縮,而且是誤差最小的方式。
精髓
數據們實施仿射變換,挲圓調正,達成白化。
每個橫條的平均數都調成0,每個橫條的變異數都調成1/n,相異橫條的共相關數都調成0。
數學公式(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ᵀ
演算法(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。
我不清楚這是否比奇異值分解來的快。
Laplacian principal component analysis
graph Laplacian analysis【尚無正式名稱】
數據兩兩距離平方總和盡量小。
然而正解缺乏討論意義:Y的元素通通相等。
minimize tr(YLYᵀ) given L. find Y.
改為討論歪解:次小的特徵值的特徵向量們,作為Y的橫條。
由於特徵向量互相垂直,追加限制條件YYᵀ = I。
⎰ minimize tr(YLYᵀ) ⎱ YYᵀ = I given L. find Y.
principal component analysis
數據到原點的距離平方總和盡量大。
改寫成數據的線性變換。
⎰ Y = AX ⎱ YYᵀ = I given X. find A and Y.
Laplacian principal component analysis
Microsoft的專利。pen-pineapple-apple-pen。
⎧ minimize tr(YLYᵀ) ⎧ minimize tr(AXLXᵀAᵀ) ⎨ Y = AX ⎨ Y = AX ⎩ YYᵀ = I ---> ⎩ AXXᵀAᵀ = I given X and L. find A and Y. given X and L. find A and Y.
兩種思路,答案好像不太一樣,我沒有仔細研究。
一、GLA視作主角:
首先實施PCA,使得XXᵀ = YYᵀ = AAᵀ = I。
然後XLXᵀ次小的特徵值的特徵向量們,作為A的橫條。再以A和X求得Y。
二、PCA視作主角:
正確答案的左邊乘上任意的正規正交矩陣,仍是正確答案。追加GLA,從中取得更對味的答案──所有點對距離盡量均勻。
PCA當中,正確答案們的點對距離完全一致,追加GLA沒有任何意義,多此一舉。low-rank PCA當中,追加GLA就有意義。
independent component analysis
楔子
PCA:數據實施仿射變換,成為不相關數據。
ICA:數據實施仿射變換,成為獨立數據。
independent component analysis
「獨立成分分析」。數據實施仿射變換,成為獨立數據。
然而,獨立相當嚴苛。一般數據無法透過仿射變換成為獨立數據。於是數學家想到兩種折衷方案。
一、特殊數據:給定數據是實施了某種仿射變換的獨立數據。
這種特殊數據,可以透過仿射變換(反向變換),成為獨立數據(恢復原貌)。
二、獨立程度:無法完全獨立,那就盡量獨立。設計指標,評估獨立程度,化作最佳化問題。
常見指標mutual information、cross-cumulant,並未充分運用獨立的特性,既失準、又繁冗。等你發明更好的指標。
原數據是實施了某種仿射變換的獨立數據。
獨立數據的座標軸,反向變換,得到原數據的「獨立成分」。
先算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,:));