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 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

改寫成不動點遞推法,以判斷是否收斂至最佳解。

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》

演算法:一次方程組求解

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

{ ∂/∂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反向代入。
對角優勢矩陣求解:不動點遞推法、鬆弛法。(宛如梯度下降法)
對稱對角優勢矩陣求解:生成樹。(僅有理論上的價值,沒有實務上的價值。)
微分方程式求解:歐拉法。

演算法(Richardson Iteration)

「不動點遞推法」。請見本站文件「Root Finding」。

梯度下降法與不動點遞推法的差異,-step⋅2L改成了+α⋅L。

xnext = x + α ⋅ (Lx - 0)
      = x + α ⋅ Lx

演算法(Jacobi Iteration)

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

梯度下降法與雅可比法的差異,step⋅2L改成了D⁻¹L。

xnext = D⁻¹ [0 - (L-D)x]
      = D⁻¹Ax
      = D⁻¹(D - L)x
      = D⁻¹Dx - D⁻¹Lx
      = x - D⁻¹Lx
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。

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

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

Explicit Euler Method:(I - Δt⋅L)當Δt足夠小才是對稱正定矩陣。
Implicit Euler Method:(I + Δt⋅L)⁻¹是對稱正定矩陣。必定收斂。

與其以反矩陣(I + Δt⋅L)⁻¹遞推,不如直接以反矩陣L₀₁⁻¹求解。就當作是學個想法吧。

Laplacian Smoothing / Diffusion

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

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

xnext = D⁻¹Ax

亦得解釋成每一點減去鄰邊平均數。邊是兩點差距。

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

xnext = x - D⁻¹Lx;

「擴散」。每一點移動一步,根據鄰邊平均數多寡。自訂擴散率,調整移動量。當擴散率是1,擴散即是梯散平滑化。

xnext = x - rate ⋅ D⁻¹Lx;

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

minimize xᵀLx         	Dirichlet Energy Minimization
solve Lx = 0          	Laplace Equation
xnext = x + α ⋅ Lx    	Richardson Iteration
xnext = x - D⁻¹Lx     	Jacobi Iteration
xnext = x - step ⋅ 2Lx	Gradient Descent
xnext = x - Δt⋅Lx     	Explicit Euler Method

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

GLA的精髓,就是擴散!儘管最快的演算法不是擴散,而是共軛梯度法。

引入權重

一、權重:Adjacency Matrix可以引入非負權重。

兩兩距離,調整比重。權重大,引力強,擴散後,距離近。

L = D - A   where A is weighted adjacency matrix
                  D is weight sum

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

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

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

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

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

類似概念出現於Locally Linear Embedding。

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

演算法:特徵分解

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

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

最小值缺乏討論意義,改為討論次小值。又為了避免各維座標均等,改為討論前d個次小值,d是維度。

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

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

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

【尚無正式名稱】

給定原數據,產生新數據,兩兩相對位置盡量相同。

設計Adjacency Matrix,控制哪些相對位置需要納入計算。為了讓新數據有所改變,必須追加限制條件,例如自訂部分數據座標。

兩兩向量差異平方總和盡量小。

min sum ‖(pᵢ - pⱼ) - (p̕ᵢ - p̕ⱼ)‖²
 P̕ (i,j)

引入Adjacency Matrix,只考慮部分的兩兩向量。

min sum Aᵢⱼ ‖(pᵢ - pⱼ) - (p̕ᵢ - p̕ⱼ)‖²
 P̕ (i,j)

多維問題等價於多個一維問題。

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

一維兩兩向量差異平方總和,可以改寫成二次函數。

  sum Aᵢⱼ ‖(xᵢ - xⱼ) - (x̕ᵢ - x̕ⱼ)‖²
= sum Aᵢⱼ ‖(xᵢ - x̕ᵢ) - (xⱼ - x̕ⱼ)‖²
= (x-x̕)ᵀL(x-x̕)
= xᵀLx - 2xᵀLx̕ + x̕ᵀLx̕

一維兩兩向量差異平方總和盡量小,可以改寫成二次最佳化。

  min sum Aᵢⱼ ‖(xᵢ - xⱼ) - (x̕ᵢ - x̕ⱼ)‖²
= min sum Aᵢⱼ ‖(xᵢ - x̕ᵢ) - (xⱼ - x̕ⱼ)‖²
= min (x-x̕)ᵀL(x-x̕)
= min xᵀLx - 2xᵀLx̕ + x̕ᵀLx̕
= min x̕ᵀLx̕ - 2bᵀx̕ + c     where b = Lx, c = xᵀLx

極值位於一次微分等於零的地方。可以改寫成一次方程組求解。

該一次方程組稱為Poisson Equation。

solve Lx̕ = Lx

二次最佳化,無法使用特徵分解求得極值。

經典應用

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

embedding:數據是點座標。降低維度,盡量均勻攤平。

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

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

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)

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

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

幾何意義

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 trace(YLYᵀ)
given L. find Y.

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

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

{ minimize trace(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 trace(YLYᵀ)             { minimize trace(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,:));