Centering

Centering

「中心化」。數據平均值挪至原點。數據減去平均值。

每個維度分開處理。

中心化可以寫成矩陣運算,然而非常複雜,充滿人類的惡意。

All-one Matrix

「全一矩陣」。元素全是1的矩陣。

    [ 1 1 1 ]
𝟏 = [ 1 1 1 ]
    [ 1 1 1 ]

可以用來計算總和、平均值、離差。

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

data:               one:
                        [ 1 1 1 1 1 ]
    [ 1 4 5 3 7 ]       [ 1 1 1 1 1 ]
X = [ 2 5 0 3 5 ]   𝟏 = [ 1 1 1 1 1 ]
    [ 3 6 4 3 9 ]       [ 1 1 1 1 1 ]
                        [ 1 1 1 1 1 ]

sum:
       [ 20 20 20 20 20 ]
Xsum = [ 15 15 15 15 15 ] = X𝟏
       [ 25 25 25 25 25 ]

mean:
    [ 4 4 4 4 4 ]
X̄ = [ 3 3 3 3 3 ] = (X𝟏)/n = X(𝟏/n)
    [ 5 5 5 5 5 ]

deviation:
    [ 1-4 4-4 5-4 3-4 7-4 ]
X̂ = [ 2-3 5-3 0-3 3-3 5-3 ] = X - X(𝟏/n) = X(I - 𝟏/n) = XC
    [ 3-5 6-5 4-5 3-5 9-5 ]

中心化的結果,就是離差。

Centering Matrix

「中心化矩陣」。方便計算離差的矩陣。

              [ 1 0 0 ]     [ 1 1 1 ]
C = I - 𝟏/n = [ 0 1 0 ] - ⅓ [ 1 1 1 ]
              [ 0 0 1 ]     [ 1 1 1 ]
數據平均值 X(𝟏n×n/n)
維度平均值 (𝟏d×d/d)X
X的各個直條,減去各自平均值:X - X(𝟏/n) = X(I - 𝟏/n) = XC
X的各個橫條,減去各自平均值:X - (𝟏/d)X = (I - 𝟏/d)X = CX

全一矩陣、中心化矩陣沒有反矩陣。

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(Under Construction!)

Graph Laplacian Analysis【尚無正式名稱】

PCoA:已知兩兩距離,求得數據。

GLA:令兩兩距離盡量小,求得數據。

「圖梯散分析」。憑空產生大量數據,兩兩距離平方總和盡量小。設計Adjacency Matrix,控制哪些距離需要納入計算。為了不讓數據集中於一點,必須追加限制條件,例如自訂部分數據座標。

Adjacency Matrix

「相鄰矩陣」。請見本站文件「Adjacency Matrix」。

兩點有關,邊為1。兩點無關,邊為0。

相鄰矩陣是對稱矩陣。

兩兩距離平方總和

兩兩距離平方總和。

 sum ‖xᵢ - xⱼ‖²
(i,j)

兩兩距離平方總和盡量小。

min sum ‖xᵢ - xⱼ‖²
 x (i,j)

引入Adjacency Matrix,只考慮部分的兩兩距離。

min sum Aᵢⱼ ‖xᵢ - xⱼ‖²
 x (i,j)

一個距離平方,可以拆解成各維距離平方。各維各自統計距離平方,各維各自最佳化,結果仍相同。多維問題等價於多個一維問題。

Laplacian Matrix

「梯散矩陣」。請見本站文件「Laplacian Matrix」。

梯散矩陣L=邊數矩陣D-相鄰矩陣A。

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

梯散矩陣是對稱半正定矩陣。

一維兩兩距離平方總和

一維兩兩距離平方總和,可以改寫成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 trace(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
Gradient Descent:各座標同時移動。
Coordinate Descent:各座標輪流移動。

演算法(Conjugate Gradient Method)

「共軛梯度法」。請見本站文件「Quadratic Optimization」。

實務上最快的演算法。調整誤差容許範圍,提早結束演算法。

甚至追加外掛,採用Multigrid Preconditioned Conjugate Gradient Method。我沒有仔細讀懂,有勞各位自立自強。

《Efficient Preconditioning of Laplacian Matrices for Computer Graphics》

演算法:一次方程組求解

目標函數,分別對每筆數據偏微分,使之為零。

{ ∂/∂x₀ sum Aᵢⱼ ‖xᵢ - xⱼ‖² = 0
{ ∂/∂x₁ sum Aᵢⱼ ‖xᵢ - xⱼ‖² = 0
{            :               :

目標函數,改寫成二次型,對向量微分,使之為零。結果一樣。

solve 2Lx = 0

係數2可以消去。等號兩邊同乘1/2。

solve Lx = 0

追加限制條件:自訂部分數據座標。

solve Lx = 0   subject to x₀ = c₀, x₁ = c₁, ...

已知數移項併入常數項,矩陣尺寸縮小,形成新的方程組。

solve [ L₀₀ L₀₁ ] [ c ] = [ 0 ]
      [ L₁₀ L₁₁ ] [ x ]   [ 0 ]
solve L₀₀ c + L₀₁ x = 0
solve L₀₁ x = - L₀₀ c

演算法一覽。

矩陣求解:高斯消去法。
對稱半正定矩陣求解:Cholesky Decomposition反向代入。
對角優勢矩陣求解:鬆弛法。(宛如梯度下降法)
對稱對角優勢矩陣求解:生成樹。(僅有理論上的價值,沒有實務上的價值。)
微分方程式求解:歐拉法。

演算法(Jacobi Iteration)

「雅可比法」。請見本站文件「Rexalation」。

雅可比法宛如梯度下降法的改良版本,step⋅2L改成了D⁻¹L。

xnext = D⁻¹ [0 - (L-D)x]   Jacobi Iteration
      = D⁻¹Ax
      = D⁻¹(D - L)x
      = D⁻¹Dx - D⁻¹Lx
      = x - D⁻¹Lx          Normalized Gradient Descent
Jacobi Iteration:以舊值算新值,需要兩條陣列。
Gauss-Seidel Iteration:立即使用新值,只需要一條陣列。
Symmetric Gauss-Seidel Iteration:
修改了計算順序,原本是每回合從頭到尾,改成了頭尾往返。
新值連貫,待命時間更少,收斂速度更快。
此處的Symmetric,是指頭尾往返,不是指對稱矩陣。
SOR:新值與舊值,調整比重,浮報新值份量,漏報舊值份量。

演算法(Spielman-Teng Algorithm)

運用生成樹切割矩陣、決定計算順序。我沒有仔細讀懂,有勞各位自立自強。

時間複雜度極低,實際執行時間極高。僅有理論上的價值,沒有實務上的價值。

Spielman-Teng Algorithm:最初版本。
KOSZ Algorithm:時間複雜度更低。
https://sites.google.com/a/yale.edu/laplacian/
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

《Nearly-Linear Time Algorithms for Preconditioning and Solving Symmetric, Diagonally Dominant Linear Systems》

演算法(Euler Method)

「歐拉法」。請見本站文件「Differential Equation」。

Laplace Equation改成Diffusion Equation。Lx = 0改成-Lx = d/dt x。

令d/dt x = 0,得到原式Lx = 0。令x隨著時間流逝逐漸收斂至定值,使得「-Lx = d/dt x的收斂結果」等於「Lx = 0的解」。

引入時間變數,將問題化作微分方程式求解,採用歐拉法。

歐拉法宛如梯度下降法的改良版本,step⋅2L改成了Δt⋅L。

Δt宛如步伐大小,需要自訂大小。

Explicit Euler Method:取舊值。
(xnext - x) / Δt = -Lx
xnext = x - Δt⋅Lx
xnext = (I - Δt⋅L)x

Implicit Euler Method:取新值。
(xnext - x) / Δt = -Lxnext
xnext = (I + Δt⋅L)⁻¹x

再改寫成特徵向量演算法Power Iteration,以判斷收斂。

對稱正定矩陣,存在N個實數特徵向量,足以支配整個空間。任意初始向量(不含零向量)必定趨近最大的特徵值的特徵向量。

Explicit Euler Method:(I - Δt⋅L)當Δt太大就不是對稱正定矩陣。不一定收斂。
Implicit Euler Method:(I + Δt⋅L)⁻¹是對稱正定矩陣。必定收斂。

不想煩惱步伐大小、想要保證收斂,代價是預先計算反矩陣。

與其以反矩陣(I + Δt⋅L)⁻¹遞推,不如直接以反矩陣L⁻¹求解。

演算法:特徵分解

L是對稱半正定矩陣,xᵀLx的最小值是0,位於最小的特徵值0的特徵向量[1,1,1,...]ᵀ。

最後求得數據x = [k,k,k,...]ᵀ,k是任意數。數據座標均等,集中於一點。

最小值缺乏討論意義,改為討論次小值。

梯散矩陣,保留次k大的特徵值們,其餘歸零。特徵分解L = EΛEᵀ右半邊√ΛEᵀ,即是k維數據矩陣。

另外,梯散矩陣等於關聯矩陣的轉置的內積。L = MMᵀ。

關聯矩陣,保留次k大的奇異值們,其餘歸零。奇異值分解M = UΣVᵀ左半邊UΣ,即是k維數據矩陣。

Laplacian Smoothing / Diffusion

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

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

xnext = D⁻¹Ax

「擴散」。每一點移動一步,根據自身與鄰點的差異平均值。自訂擴散率,調整移動距離。

各點鄰邊和Lx。各點鄰邊平均D⁻¹Lx。

xnext = x - ratio ⋅ D⁻¹Lx;

當擴散率是1,擴散等價於梯散平滑化。

擴散讓相鄰距離逐漸均勻、讓兩兩距離總和逐漸變小。

minimize xᵀLx    	Dirichlet energy minimization
solve Lx = 0     	Laplace equation
xnext = D⁻¹Ax    	Laplacian smoothing
xnext = x - D⁻¹Lx	diffusion

梯度下降法、鬆弛法、歐拉法,皆是擴散。各點擴散率不同,各點擴散次序不同。

GLA的精髓,就是擴散。

權重操作

Adjacency Matrix可以引入權重。各個兩兩距離,調整比重。

從擴散的觀點來看,權重有著幾何意義。

一、正規化:各邊除以鄰邊總和(或者數量,當權重為1)。

擴散時,剛好走到鄰點的加權平均值,避免鄰點較多而走較遠。

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

二、加權平均座標:各種加權平均座標都可以當作權重。

擴散時,調整行進方向。

三、擴散無效化:給定數據,求得權重,使得處處梯散為零。

這樣的權重,性質與用途皆不明。等你發論文。

類似概念出現於Locally Linear Embedding。

given x, find A, such that Lx = 0.

延伸應用

drawing:數據是點座標。盡量均勻。

blending:數據是像素值。先後梯度盡量相等。

deformation:數據是點座標。先後形狀盡量相等。

embedding:數據是點座標。先後距離盡量相等。

Whitening(Under Construction!)

Whitening

「白化」。數據相關矩陣縮放為一。數據整體除以相關矩陣。

讓相同維度的相關數是1,相異維度的相關數是0。換句話說,共相關矩陣是I。

【註:統計學家定義成變異數,也就是說,數據預先中心化。】

白化結果無限多種。任意一種白化結果,經過正規正交矩陣(旋轉暨鏡射),仍是白化結果。

因此白化方式有許多種,其中一種是正規化。

W = √XXᵀ⁻¹。

延伸閱讀:數學符號

由於是新名詞新觀念,我只好自己瞎掰數學符號了。

X̄ = X(𝟏/n)  	mean
X̂ = X - X̄ = XC	deviation (centerized)
‖X‖ = √XᵀX	absolute value
⟪X⟫ = X / ‖X‖	unitary   (normalized = standardized and decorrelated)

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ᵀ的特徵向量演算法。

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

找到主成分:共變異矩陣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

Correlation Matrix

「相關矩陣」。矩陣儲存兩兩相關數,矩陣元素Aᵢⱼ是i維與j維的點積。

維度的兩兩共相關數(X橫條配Y橫條),得到共相關矩陣XYᵀ。

分成數據版本和特徵版本?【尚待確認】

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

相關矩陣細分為自相關矩陣、互相關矩陣。

Autocorrelation:自相關矩陣。相同數據集的相關矩陣。XXᵀ。
Cross-correlation Matrix:互相關矩陣。相異數據集的相關矩陣。XYᵀ。

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

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

直線擬合(主成分分析)是自相關矩陣XXᵀ。

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 Principal Component Analysis(Under Construction!)

備忘

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

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

Graph Laplacian Analysis【尚無正式名稱】

新數據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和GLA合體!

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

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

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

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

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

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

延伸閱讀:矩陣降秩

PCoA屬於Low-rank Matrix Approximation的特例:XᵀX是對稱正定矩陣。

min ‖YᵀY - XᵀX‖²

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

PCA和PCoA屬於矩陣降秩的特例: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了?

Normalization(Under Construction!)

Normalization

一個向量的正規化

一個矩陣的正規化

一堆數據的標準化

統計學當中,數據的概觀模樣,有兩個重要指標:

平均值:數據的中心。數據總和,除以數據數量。
變異數:數據的疏密。數據減去平均值,平方,總和。

數據可以標準化,總共兩個步驟。

中心化:挪動位置至原點0。減去平均值。
標準化:縮放範圍至長度1。除以變異數的平方根。
平均值是零xᵀ1⃗ = 0、變異數是一‖x‖² = xᵀx = 1。

一堆數據的標準化【尚無正式名稱】

中心化、白化。

因為沒有正式名稱,所以大家稱作白化。

Correlation(Under Construction!)

Correlation

「相關」。兩堆數據的類似程度。

大家習慣先做中心化、白化,校正位置、範圍,再來計算兩堆數據的類似程度。

cosine similarity = cross correlation
a dot b / |a||b| = a/|a| dot b/|b|
= cov(a/|a|,b/|b|)
= cov(a/sqrt(var(a)),b/sqrt(var(a)))
= cov(a,b) / sqrt(var(a)) sqrt(var(b))
= σab / σa σb

matrix version use outer product!
var(X) = XXᵀ    squared length
cov(X,Y) = XYᵀ

XYᵀ / |X||Y| = X/|X| (Y/|Y|)ᵀ
= cov(X,Y) / sqrt(var(a)) sqrt(var(b))

normalization   make squared length equals 1.
standardization make variance       equals 1.

名稱有許多種:

Partial Least Squares Regression  沒有白化
Orthonormalized PLS Regression    X白化
Canonical Correlation Analysis    X與Y都白化

解讀方式有許多種:

一、加權平均值:
  Xa與Yb,令相關係數(餘弦相似度)越大越好。
  Cov是點積、Var開根號是長度。
二、直線擬合:
  首先轉置一下。X投影到直線a,Y投影到直線b,令相關係數越大越好。
三、主成分分析:
  PCA(Xᵀ)的變換矩陣(係數)、PCA(Yᵀ)的變換矩陣(係數),實施多變量迴歸。

數據的共變異矩陣XᵀY的奇異值分解,奇異值越大越好。

演算法有奇異值分解或者NIPALS。

https://www.ijcai.org/Proceedings/09/Papers/207.pdf
http://users.cecs.anu.edu.au/~kee/pls.pdf
https://stats.stackexchange.com/questions/23353/

最佳化形式:

                                Cov[aᵀX, bᵀY]           
max Corr[aᵀX, bᵀY] = max ———————————————————————————————
a,b                  a,b  sqrt(Var[aᵀX]) sqrt(Var[bᵀY]) 
                                   aᵀXYᵀb           
                   = max ———————————————————————————
                     a,b  sqrt(aᵀXXᵀa) sqrt(bᵀYYᵀb) 

直線可以推廣成平面、空間、……。

max tr(AᵀXYᵀB)   subject to AᵀXXᵀA = I, BᵀYYᵀB = I
A,B

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(主成分向量)。兩者的維度的共相關矩陣,XXᵀ和vvᵀ,點積盡量大。

Graph 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是座標,GLA是圖。你能銜接兩者嗎?