Alignment

Alignment

「對齊」。一堆數據,經過變換,盡量符合另一堆數據。已知甲堆數據、乙堆數據,求函數。

函數具有特殊限制,增進討論意義。以下討論相似變換。

Error

兩堆數據通常不一致。強硬地對齊,就會有「誤差」。

兩堆數據的誤差,有許多種衡量方式:

一、窮舉甲堆到乙堆的所有兩兩配對,累計距離。
二、窮舉甲堆每一筆數據,分別找到乙堆之中距離最近的那筆數據,累計距離。
三、窮舉乙堆每一筆數據,分別找到甲堆之中距離最近的那筆數據,累計距離。
四、二與三相加。
五、二與三聯集。
六、上述誤差,改為只累計前K短的距離。
七、上述誤差,只考慮互相鄰近的K筆數據。(部分銜接的效果)

隨隨便便就能發明一大堆誤差公式。重點在於計算時間長短、銜接密合程度。哪種誤差比較好?目前還沒有定論。

演算法(Principal Component Analysis)

兩堆數據各自計算PCA座標軸,然後對齊。

http://www.cs.tau.ac.il/~dcor/Graphics/cg-slides/svd_pca.pptx
http://www.cse.wustl.edu/~taoju/cse554/lectures/lect07_Alignment.pdf
http://perception.inrialpes.fr/~Horaud/Courses/pdf/Horaud_3DS_5.pdf

兩堆數據各自求得共變異矩陣,實施特徵分解。兩邊的特徵向量,依照特徵值大小一一對應。

最佳的位移量:數據中心的差距。最佳的旋轉量:特徵向量的夾角。最佳的縮放量:特徵值平方根的比例。

特徵向量不具方向性。指定每個特徵向量的方向,滿足右手座標系,才能求得旋轉矩陣,不含鏡射。

一種想法是窮舉各種可能的方向,從中找到旋轉角度最小者。

另一種想法是檢查determinant。det(E) = +1,是右手座標系。det(E) = -1,是左手座標系;只要再讓任意一個特徵向量顛倒方向,就是右手座標系。選擇最小的特徵值的特徵向量,讓誤差增加最少。

PCA座標軸容易因outlier而歪斜。是非常不精準的方法。

演算法(Procrustes Analysis)

當兩堆數據一樣多,並且預先知道數據對應關係,就有更精準的做法。

https://igl.ethz.ch/projects/ARAP/svd_rot.pdf
http://perception.inrialpes.fr/~Horaud/Courses/pdf/Horaud_3DS_6.pdf
    [ x1.x x2.x ... xN.x ]       [ y1.x y2.x ... yN.x ]
X = [ x1.y x2.y ... xN.y ]   Y = [ y1.y y2.y ... yN.y ]
    [ x1.z x2.z ... xN.z ]       [ y1.z y2.z ... yN.z ]
                                  _____________________ 
solve s R X + t = Y              | t: translate vector |
                            2    | s: scalar vector    |
minimize ‖ Y - (s R X + t) ‖     | R: rotation matrix  |
                                  ‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾ 
then t = Ȳ - s R X̄
      2      T         T          2        2
     s = ∑ Ŷi Ŷi ÷ ∑ X̂i X̂i = ‖ Ŷ ‖F ÷ ‖ X̂ ‖F
            T                  T       T
     R = V U    with C₃ₓ₃ = X̂ Ŷ = U Σ V 

令函數是位移、旋轉、縮放三個步驟。方便起見,把位移改為最終步驟。數據套用函數之後,令誤差越小越好。誤差訂為對應數據的距離的平方的總和。

最佳的位移量:兩堆數據各自的平均值相減,相減之前先讓第一個平均值套用函數。最佳的縮放量:兩堆數據各自的變異數相除。最佳的旋轉量:兩堆數據各自減去平均值(中心位移至原點)、除以變異數(大小縮放為相等),求得維度的共變異矩陣,實施奇異值分解,然後VUᵀ矩陣相乘,即是旋轉矩陣。

R是正規正交矩陣,功效是旋轉暨鏡射。大家通常不想鏡射。若要避免鏡射,則必須滿足右手座標系。

det(R) = +1,是右手座標系。det(R) = -1,是左手座標系;只要再讓任意一個向量顛倒方向,就是右手座標系。選擇最小的特徵值所對應的向量,讓誤差增加最少。

演算法(Iterative Closest Point)

一、初始化:
 甲、實施PCA,求得對齊函數。
 乙、甲堆套用函數。讓甲堆靠近乙堆。
二、重複此步驟,直到最近點不變為止:
 甲、甲堆每一點各自找到在乙堆的最近點,建立對應關係。
   無視其餘點,實施PA,求得對齊函數。
 乙、甲堆套用函數。讓甲堆更加靠近乙堆。

尋找最近點,簡易的方式是窮舉法,進階的方式是乙堆套用分割區域的資料結構,諸如Uniform Grid、KD-Tree。請見本站文件「Region」。

演算法(RANSAC)

一、甲堆隨機抓幾點,乙堆隨機抓幾點,點數一樣多,推定它們依序一一對應。
  以這些點實施Procrustes Analysis,求得對齊函數。
二、計算「甲堆套用函數」與「乙堆」的誤差大小。
三、重複一二,取誤差最小者。

絕聖棄智,民利百倍。

Manifold Alignment

https://en.wikipedia.org/wiki/Manifold_alignment

Embedding(Under Construction!)

Embedding(Dimensionality Reduction)

「嵌入」。更換數據所在空間。例如從三維換成二維。例如從二維換成三維。例如從球面換成平面。例如從曲面換成環面。

嵌入方式有許多種。例如投影就是其中一種嵌入方式。例如訂立座標系統並且實施座標轉換,也是其中一種嵌入方式。

嵌入時,我們通常希望儘量保留數據的原始特性。特性可以是數據的聚散程度、數據之間的距離遠近。計算學家針對各種特性,設計各種演算法。

這裡先不談太廣,僅討論最簡單的情況:數據從N維空間換成M維空間。細分為三種情況:維度變小、維度不變、維度變大。

後兩者缺乏討論意義──數據不變、數據補零不就好了。因此我們僅討論維度變小。此時「嵌入」的意義就變成了:找到一個函數,讓一堆數據經過函數變換,減少數據維度。

演算法(Principal Component Analysis)

嵌入時,採用垂直投影。

實施Low-rank PCA,保留大的、捨棄小的特徵值暨特徵向量,只做投影、略過位移、略過縮放。

演算法(Metric Multidimensional Scaling)

嵌入時,盡量保留所有點對距離。

首先計算原數據X的兩兩距離,做為新數據Y的兩兩距離。而兩兩距離可以寫成「距離矩陣」:Mᵢⱼ = ‖pᵢ - pⱼ‖。

問題變成:已知距離矩陣M、求得數據Y。

Mᵢⱼ² = ‖qᵢ - qⱼ‖² = ‖qᵢ‖² + ‖qⱼ‖² - 2‖qᵢ ∙ qⱼ‖,橫條擁有相同‖qᵢ‖²,直條擁有相同‖qⱼ‖²。嘗試消去它們:每個橫條減去橫條總平均,然後每個直條減去直條總平均。不失一般性,假設原點位於數據中心,qᵢ總和是0,如此便剩下-2‖qᵢ ∙ qⱼ‖這一項。全部元素再除以-2,得到‖qᵢ ∙ qⱼ‖。最終形成了矩陣YᵀY,好算多了!簡言之:距離平方矩陣Mᵢⱼ²,行列中心化,除以-2,變成矩陣YᵀY。

問題變成:給定YᵀY,求得Y。

仿照PCA的手法,YᵀY實施特徵分解,得到答案Y = √ΛEᵀ。保留大的、捨棄小的特徵值和特徵向量,以符合新維度,同時也讓誤差最小。

YᵀY = EΛEᵀ = E√Λ√ΛEᵀ = (E√Λ)(√ΛEᵀ) = (√ΛEᵀ)ᵀ(√ΛEᵀ)

還有另一種手法,YᵀY實施Cholesky分解(對稱矩陣的LU分解),得到答案Y = Lᵀ。然而降維時誤差大,沒人這樣做。

YᵀY = LLᵀ

演算法(Random Projection)

https://en.wikipedia.org/wiki/Random_projection
https://en.wikipedia.org/wiki/Johnson–Lindenstrauss_lemma
https://homes.cs.washington.edu/~thickstn/docs/fast_jlt.pdf

線性函數的本質是旋轉、鏡射、縮放,數據的相對位置不變。隨便找個線性函數做為嵌入函數,都能大致保留數據的相對位置,但是會大幅改變數據的相對距離。

而隨機導致分散,分散導致正交,正交導致距離不變。

絕聖棄智,民利百倍。

Manifold Embedding(Manifold Learning)

更換數據所在空間,從N維空間的流形換成M維空間。

一、選擇鄰居:全部數據、Fixed-radius Near Neighbor、K-Nearest Neighbor。
二、套用核函數:無、Gaussian Kernel、Heat Kernel。
三、保留特性:距離、梯散、內插函數。

核函數採用鐘形函數,以對付流形:超過一定距離,鄰居關係迅速減弱。

各點邊數總和Wᵢⱼ = 1 / ‖xᵢ - xⱼ‖²。

演算法清單:http://www.cs.cmu.edu/~liuy/distlearn.htm

演算法(Isomap)

嵌入時,盡量保留所有點對距離。

一、K-Nearest Neighbor,盡量讓距離較近的數據成為鄰居。
二、All-Pairs Shortest Paths,求出所有點對距離。
三、Metric Multidimensional Scaling,盡量保留所有點對距離。

演算法(Laplacian Eigenmap)

嵌入時,盡量讓相鄰數據距離均衡。

一、K-Nearest Neighbor,盡量讓距離較近的數據成為鄰居。
二、Laplacain Smoothing,盡量讓相鄰數據距離均衡。

Laplacian Smoothing改寫成Dirichlet Energy Minimization。

minimize yᵀLy

各點邊數不一致、梯散不等比、距離不均衡。因此實施正規化。

正解是最小的特徵值暨特徵向量,缺乏討論意義。改取次小的特徵值的特徵向量們。M個特徵向量構成M維數據。v變回y。

minimize yᵀ√D⁻¹L√D⁻¹y
minimize vᵀLv           (v = √D⁻¹y)

令特徵向量長度為一,化作廣義特徵值。

minimize yᵀ√D⁻¹L√D⁻¹y          subject to yᵀy = 1 
minimize vᵀLv                  subject to vᵀDv = 1   (v = √D⁻¹y)
minimize vᵀLv - λ (vᵀDv - 1)
solve Lv - λDv = 0
solve Lv = λDv

延伸閱讀:Laplacian Eigenmap原始版本

原始論文邏輯不通、思維跳躍。我試著腦補一下。

原始論文採取了其他解:令平均值是零、變異數是一。

minimize yᵀ√D⁻¹L√D⁻¹y   subject to yᵀ𝟏 = 0 and yᵀy = 1 
minimize vᵀLv           subject to vᵀ√D𝟏 = 0 and vᵀDv = 1

投影梯度法,求得最佳解位置,只能得到一維數據,無法得到高維數據(所有維度數值相同實在太醜)。折衷方式是特徵向量。

其中一個限制條件當作沒看到,另一個限制條件併入目標函數。化作廣義特徵值,取次小的特徵值的特徵向量們。M個特徵向量構成M維數據。

minimize vᵀLv                  subject to vᵀ√D𝟏 = 0 and vᵀDv = 1
minimize vᵀLv - λ (vᵀDv - 1)   subject to vᵀ√D𝟏 = 0
solve Lv - λDv = 0             subject to vᵀ√D𝟏 = 0
solve Lv = λDv                 subject to vᵀ√D𝟏 = 0

原始論文聲稱限制條件是vᵀD𝟏 = 0與vᵀDv = 1。反正一樣。

solve Lv = λDv                 subject to vᵀD𝟏 = 0

v直接當作答案,不必變回y。v誤植為y。

solve Ly = λDy                 subject to yᵀD𝟏 = 0

演算法(Linear Laplacian Eigenmap)

Laplacian Eigenmap改良版。嵌入時,必須是線性變換。

線性變換的優點:數據相對位置大致相同。

LA追加線性變換,效果等同PCA暨LA。

一、K-Nearest Neighbor,盡量讓距離較近的數據成為鄰居。
二、Linear Transformation & Laplacian Smoothing,
  數據實施某種線性變換,盡量讓相鄰數據距離均衡。

二次型多維版本,代入Y = AX,事情迎刃而解。

求得X√D⁻¹L√D⁻¹Xᵀ,取次小的特徵值的特徵向量們,M個特徵向量做為A的橫條,再以A和X求得Y。

minimize Y√D⁻¹L√D⁻¹Yᵀ      subject to Y = AX
minimize AX√D⁻¹L√D⁻¹XᵀAᵀ

令特徵向量長度為一、互相垂直,化作廣義特徵值。

minimize Y√D⁻¹L√D⁻¹Yᵀ      subject to Y = AX and YYᵀ = I
minimize AX√D⁻¹L√D⁻¹XᵀAᵀ   subject to AXXᵀAᵀ = I
minimize AX√D⁻¹L√D⁻¹XᵀAᵀ - λ (AXXᵀAᵀ - I)
solve X√D⁻¹L√D⁻¹Xᵀa - λXXᵀa = 0
solve (X√D⁻¹L√D⁻¹Xᵀ)a = λ(XXᵀ)a

原始論文令V = X√D⁻¹,式子稍微清爽。

solve VLVᵀa = λVDVᵀa   (V = X√D⁻¹)

原始論文繼承了Laplacian Eigenmap的錯誤,V誤植為X。

solve XLXᵀa = λXDXᵀa

演算法(Self-organizing Map)

https://medium.com/@CCstruggled/d212534160e0

嵌入時,盡量讓相鄰數據距離縮短、互相靠近?

找到最相像的節點,更新鄰居,有如注水。

似乎可以解釋成Local Laplacian Smoothing?

演算法(Poisson Embedding)【尚待確認】

嵌入時,盡量保持梯度相同(梯散相同)。

一、K-Nearest Neighbor,盡量讓距離較近的數據成為鄰居。
二、Poisson Equation,盡量保持梯散相同。

我不知道這要怎麼解。我還沒有看人家解過。

minimize (Xᵀ-Yᵀ)L(Xᵀ-Yᵀ)
solve LXᵀ = LYᵀ
minimize ‖M - LY‖²
minimize MᵀM - MᵀLY - (LY)ᵀM + (LY)ᵀLY

調整Laplace Matrix權重,調整梯散大小。

演算法(Locally Linear Embedding)

https://www.cnblogs.com/xbinworld/archive/2012/07/09/LLE.html

Poisson Equation改良版。嵌入時,盡量保持梯散為零。

一、K-Nearest Neighbor,盡量讓距離較近的數據成為鄰居。
二、Laplace Equation,針對原數據,令梯散為零,求得權重。
三、Laplace Equation,針對新數據,盡量保持梯散相同。

數據減去相鄰數據的加權平均值,就是梯散!

minimize ‖Y - YWᵀ‖²   subject to YYᵀ = I
minimize YMYᵀ         subject to YYᵀ = I   where M = (I-W)ᵀ(I-W)
solve My = λy
minimize ‖LYᵀ‖²       subject to YYᵀ = I   where L = I - W
minimize YLᵀLYᵀ       subject to YYᵀ = I
solve LᵀLy = λy

演算法(Linear Locally Linear Embedding)

Locally Linear Embedding改良版。嵌入時,必須是線性變換。

solve (XMXᵀ)a = λ(XXᵀ)a     where M = (I-W)ᵀ(I-W)

演算法(Stochastic Neighbor Embedding)

http://www.datakit.cn/blog/2017/02/05/t_sne_full.html

嵌入時,盡量保持分布相同。

原數據、新數據,分別求得高斯混和模型?令兩個分布的距離盡量小,距離採用Kullback-Leibler Divergence。

核函數是常態分布,即是此演算法,簡稱SNE。

核函數是學生t分布,效果更好,簡稱t-SNE。

延伸閱讀:Stochastic Neighbor Embedding等價觀點

嵌入時,盡量保持內插函數相同。

數據添加一個新維度做為函數值,一律是1。原數據X、新數據Y,分別實施Inverse Distance Weighting Interpolation(原始論文的分母少了一項)。令兩個內插函數的距離盡量小。距離是N個數據的函數值,函數值取log相減,越大越好。

核函數是指數函數,超過一定距離,距離迅速暴增。

Unsupervised Manifold Embedding(Clustering)

數據沒有類別。嵌入時,盡量同類相聚。

演算法(Spectral Clustering)

一、K-Nearest Neighbor,盡量讓距離較近的數據成為鄰居。
二、Laplacian Smoothing,盡量讓相鄰數據距離均衡。
三、K-Means Clustering,盡量讓相鄰數據歸類於同一組。

前兩步驟正是Laplacian Eigenmap:Laplace Matrix取次小的M個特徵向量,數據從N維降為M維。

原始論文聲稱此方法等價於Maximum Cut近似解。我想應該是哪裡搞錯了。Maximum Cut當中,向量元素必須是0或1;而特徵向量元素通常不是0或1。

Supervised Manifold Embedding(Classification)

數據具有類別。嵌入時,盡量同類相聚。

演算法(Linear Discriminative Analysis)

援引PCA的概念。

尋找一條直線,數據投影到直線之後,各類數據的平均值的變異數盡量大(異類盡量散開),各類數據各自的變異數盡量小(同類盡量聚集)。

演算法(Neighborhood Component Analysis)

https://en.wikipedia.org/wiki/Neighbourhood_components_analysis

援引Stochastic Neighbor Embedding的概念。

尋找一個線性變換,數據實施變換之後,各類數據各自求得高斯混和模型,分布總和越大越好?

鄰居:K-Nearest Neighbor。核函數:Gaussian kernel。

Deformation(Under Construction!)

Deformation

「形變」。改變數據造型。

點邊面。

平面圖、球面圖,則可以討論面。

laplace operator 點是主角。

延伸應用

1. Graph Drawing: Stress Majorization。無向圖畫成圖片。
2. Fluid Simulation: Position-based Dynamics。調整水滴位置。
3. Structural Optimization: Topology Optimization。調整模型結構。
4. Molecular Modeling: Energy Minimization。調整化學結構。

無限制

min sum ‖(yi - yj) - (xi - xj)‖²

化作Laplace Matrix。各座標軸分開計算。

答案是所有點皆相等,縮成一點,缺乏討論意義。因此必須追加限制條件,實施約束最佳化。

精確位置限制

挑選一些點,設定明確座標。

解的一部分變成了定值,演算法只能使用Normal Equation系列、Jacobi Iteration系列,無法使用Conjugated Gradient Method系列。

相對位置限制、位置限制

min sum sum wij ‖xi - xj - vij‖² + sum ui ‖xi - pi‖²
min sum (Dx-v)ᵀW(Dx-v) + (x-p)ᵀU(x-p)
W and U are diagonal matrix
D is the discrete derivative operator
solve Lx = b such that L = DᵀWD + U and b = DᵀWv + Up

化作對稱對角優勢矩陣。

嚴格的說法是化作M-Matrix:Laplace Matrix的對角線元素各自加上任意正數或零。不過M-Matrix缺乏實務上的高速演算法。

《Efficient Preconditioning of Laplacian Matrices for Computer Graphics》

相對距離限制

min sum (‖xi - xj‖ - dij)²

公式解:Stress Majorization。演算法是Metric Multidimensional Scaling。

數值模擬:Euler Method。每個點分別求得加速度、更新速度、更新座標。請見本站文件「Particle Simulation」。

相對應力限制

min sum kij (‖xi - xj‖ - dij)²
k是彈性係數 d是彈簧長度

公式解:Matrix Completion。奇異值減去倍率。

數值模擬:Euler Method。

位置變化限制、速度變化限制、能量總和限制

min sum ‖xi' - xi‖² + ‖vi' - vi‖²  subject to c(x) = 0
x x'是先後位置,v v'是先後速度。c(x)是能量。

《FEPR: Fast Energy Projection for Real-Time Simulation of Deformable Objects》

二次函數有著高速演算法:Newton-Lagrange Method與Projected Gradient Method兩者合體。

乘子法-->微分方程組--->手動微分--->係數和乘數取出來變向量--->解大矩陣
qk+1 = argmin 1/2 ‖q - qₖ‖²D  subject to c(qₖ) + ∇c(qₖ)ᵀ(q - qₖ) = 0
         q
solve λk+1: [∇c(qₖ)ᵀ D⁻¹ ∇c(qₖ)] λk+1 = c(qₖ)
solve qk+1: qk+1 = qk - D⁻¹ ∇c(qk) λk+1

剛體變換限制、相似變換限制

min sum ‖(yi - yj) - T(xi - xj)‖²

Alignment + Deformation。

剛體變換是位移、旋轉。相似變換是位移、整體縮放、旋轉。

主角是點、主角是面,兩者式子不同。此處主角是點。

演算法是Block Coordinate Descent:固定R求得y,固定y求得R,兩者不斷輪流。換句話說,Procrustes Analysis與Laplace Equation不斷輪流。

Laplace Matrix是定值、是對稱矩陣,可以預先做Cholesky Decomposition,以便快速求解。

重心是定值,可以預先處理位移,將重心移到原點。

註:此演算法有人稱作Global-Local Method。

《As-Rigid-As-Possible Surface Modeling》,剛體,座標下降法,主角是點。

《Least Squares Conformal Maps for Automatic Texture Atlas Generation》,相似,Normal Equation,主角是面。

Manifold Deformation

請見本站文件「Model」。

edge difference sum = vertex laplace sum。

二維:整個面每一點都形變,必須計算面積,形成cot。大家將cot視作laplace matrix的權重,雖然很奇怪但是很簡潔。

三維:角椎對面積微分,得到cot。角錐底部不平坦。

[PP93] Computing Discrete Minimal Surfaces and Their Conjugates

piecewise linear deformation:

min sum Area_t ‖ xi' - Tt(xi) ‖²
     t

min sum sum 1/4 cot_ij ‖ (xj' - xi') - Tt(xj - xi) ‖²
     t   3

min sum ‖ Laplace_cot(xi') - Ti(Laplace_cot(xi)) ‖²
     i

Correlation(Under Construction!)

Correlation

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

https://stats.stackexchange.com/questions/23353/

演算法(Canonical Correlation Analysis)

https://www.ijcai.org/Proceedings/09/Papers/207.pdf

數據範圍太廣,難以解析。把數據投影到直線上面,方便解析。

替數據X找一條直線a,替數據Y找一條直線b。X投影到a,Y投影到b,兩堆數據的相關係數越大越好。

相關係數即是餘弦相似度:Cov是點積、Var開根號是長度。

兩堆數據必須一樣多。點對點相乘,而非兩兩交叉相乘。

                                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

演算法(Partial Least Squares Regression)

http://users.cecs.anu.edu.au/~kee/pls.pdf

模仿PCA。維度的共相關矩陣XYᵀ投影到直線w,變異數盡量大。

max Var[wᵀXYᵀ] = max wᵀXYᵀYXᵀw
 w                w

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

max tr(WᵀXYᵀYXᵀW)   subject to WᵀXXᵀW = I
 W

Power Iteration與Hotelling Deflation求得特徵值暨特徵向量。

Ranking(Under Construction!)

Ranking(Link Analysis)

「排名」。有向圖的節點排名。

演算法(PageRank)

演算法(Hyperlink-Induced Topic Search)

演算法(SpringRank)

min sum kij (‖xi - xj‖ - dij)²
k是彈性係數 d是彈簧長度
kij = dij = 1

Recommendation(Under Construction!)

Recommendation

「推薦」。相關程度排名。

點積是相關程度。收集兩兩相關程度,形成矩陣;分解矩陣,推敲出原始事物。

topic model = collaborative filtering
TF-IDF
word2vector
matrix completion
factorization machine
tensor decomposition http://people.csail.mit.edu/moitra/docs/Tensors.pdf
Laplacian Assignment
SimRank