Question

Matrix Transformation

当函数是矩阵，经典的变换有缩放、旋转、平移等等。请参考本站文件“ 2D Transformation ”。

Representation / Factorization

“重新表示”。数据套用函数，转换成其他姿态，重新呈现。

“分解”。重新表示的反向操作，揭露数据原本姿态。

重新表示：已知原数据，求新数据、求函数。

分解：已知新数据，求原数据、求函数。

当函数拥有反函数，只要移项一下，重新表示、分解其实是同一个问题，没有主从之分。

Matrix Representation / Matrix Factorization

当函数是矩阵，则有三种套用函数的方式。

一、矩阵乘在数据左边：数据们实施线性变换。

二、矩阵乘在数据右边：数据们的加权平均值。

三、矩阵相加：数据们加上误差。

数据可以是一维，甚至多维。数据可以是一笔，甚至多笔。

一笔数据是一个直条。矩阵乘在左边，原数据与新数据从左往右依序一一对应。矩阵乘在右边，权重与新数据一一对应。矩阵相加，原数据、新数据、误差一一对应。

不熟悉线性代数的读者，请参考本站文件“ Basis ”，複习一下“矩阵切成直条”的观点。

备注

representation 通常是指“数据套用函数变成新数据”，所以矩阵乘在左边是没问题的。但是矩阵乘在右边呢？姑且算是吧。

factorization 通常是指“乘法的反向操作”，而不是这裡提到的“新数据分解成函数和原数据”。按理应当另造一词，然而我并未发现更适切的词彙，只好姑且使用 factorization。

备注

矩阵乘法其实是函数複合。多个矩阵连乘，併成一个矩阵，这是“函数複合 composition”。一个矩阵，拆成多个矩阵连乘，这是“函数分解 decomposition”。函数複合的反向操作是函数分解。

儘管 representation factorization 外观如同 composition decomposition，但是其意义是数据、而非函数。本质略有差异。

Principal Component Analysis

Principal Component Analysis

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 semi-definite.            |
A(EDE⁻¹)Aᵀ      = I   | thus XXᵀ has real non-negative eigenvalues. |
A(E√D)(√DE⁻¹)Aᵀ = I   | thus XXᵀ has orthonormal eigenbasis.        |
A = √D⁻¹E⁻¹           | let XXᵀ = EDE⁻¹ (E⁻¹ = Eᵀ)                  |
A = √D⁻¹Eᵀ             ‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾ 

XXᵀ是实数对称正半定方阵，得以实施特徵分解。

E 是特徵向量们（特徵基底），恰为正规正交矩阵，其转置矩阵恰为反矩阵。而且都是实数。

D 是特徵值们，恰为对角线矩阵。而且都是实数、都非负数。

欲求 A，就想办法让 A 和 Aᵀ能够抵消 EDEᵀ，使之成为 I。特徵基底求反矩阵，特徵值先开根号再求倒数，如此就凑出了 A。

欲求 Y，只消计算 AX。

正确答案不只一种。特徵向量对调次序、颠倒方向，特徵值开根号时取正值、取负值，都是正确答案。

补充一下，正确答案的左边乘上任意的正规正交矩阵，仍是正确答案。但由于这类答案并不实用，大家总是无视这类答案。

数据减去平均值

当 X 的中心（平均值）位于原点，则 Y 的中心也将位于原点。

大家总是预先将 X 减去平均值，将 X 的中心挪至原点。好处是：一、降低数值范围，以减少浮点数运算误差。二、具备直线拟合效果，详见后面章节。

X 减去平均值之后，XXᵀ和 YYᵀ变成“维度的共变异矩阵”。

演算法（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̂ᵀ = [ d₂⋅d₁ d₂⋅d₂ d₂⋅d₃ ] = [ 3 × 3 ] = E D Eᵀ
       [ d₃⋅d₁ d₃⋅d₁ d₃⋅d₃ ]   [       ]

    [       ]       [ λ₁ 0  0  ]
E = [ 3 × 3 ]   D = [ 0  λ₂ 0  ]
    [       ]       [ 0  0  λ₃ ]

              [ 1/√λ₁   0     0   ] [       ]
A = √D⁻¹ Eᵀ = [   0   1/√λ₂   0   ] [ 3 × 3 ]
              [   0     0   1/√λ₃ ] [       ]

演算法（Singular Value Decomposition）

X 实施奇异值分解，恰好可以取代 XXᵀ实施特徵分解。

X 的奇异值、基底，恰好可以兜出 XXᵀ的特徵值、特徵基底。

X = UΣVᵀ
XXᵀ = UΣVᵀ(UΣVᵀ)ᵀ = UΣVᵀVΣᵀUᵀ = UΣΣᵀUᵀ = U(ΣΣᵀ)Uᵀ = EDEᵀ
ΣΣᵀ = D, U = E

一般来说，实数对称正半定方阵 XXᵀ的特徵分解，比普通矩阵 X 的奇异值分解来得快。可是当 X 是超大矩阵的情况下，预先计算 XXᵀ相当费时，而奇异值分解不必计算 XXᵀ，逆转胜。

演算法（EM-PCA）（Direct Linear Transformation）

Generate random matrix A. Loop over until A converge to PCA basis:
  E-step:  X = A⁺Y  = (A Aᵀ)⁻¹ Aᵀ Y
  M-step:  A = Y X⁺ = X (XᵀX )⁻¹ Xᵀ 

明明是 Fixed Point Iteration，结果被叫做 EM Algorithm。

我不清楚这是否比奇异值分解来的快。

几何意义

PCA 找到了全新的垂直座标轴：原点是数据中心、座标轴是特徵基底、座标轴单位长度是特徵值平方根。

所有数据一齐位移、旋转（与镜射）、缩放（与镜射）。

一、位移：数据中心位移至原点。

二、旋转：以特徵基底，逆向旋转数据。

三、缩放：各维度除以各特徵值平方根。

最后，各维度的平均数均为 0，变异数均为 1。

正确答案不只一种。特徵向量（连同特徵值）对调次序，效果是镜射暨旋转。特徵向量颠倒方向、特徵值变号，效果是镜射。

正确答案可以包含镜射，也可以不包含镜射。如果讨厌镜射，就让特徵基底成为右手座标系、让特徵值平方根皆是正值。

补充一下，正确答案的左边乘上任意的正规正交矩阵，也就是旋转（与镜射），仍是正确答案。但由于这类答案并不实用，大家总是无视这类答案。

让特徵基底成为右手座标系

一、标准座标轴（单位矩阵 I），是右手座标系。

二、右手座标系经过旋转，仍是右手座标系。

三、偶数次镜射，得视作旋转。

四、右手座标系经过偶数次镜射，仍是右手座标系。

五、右手座标系经过奇数次镜射，若再增加一次镜射，就是右手座标系。

determinant 可以判断镜射次数。det(E) = +1 是偶数次，det(E) = -1 是奇数次。当 det(E) = -1，则再增加一次镜射，就是右手座标系。例如任取一个特徵向量颠倒方向（元素通通添上负号）。

几何特性

一、所有数据投影到座标轴之后的变异数，第一座标轴最大，第二座标轴次大（须垂直于先前座标轴），……。

也就是说，座标轴的方向，就是数据最散开的方向。

二、所有数据到座标轴的距离平方和，第一座标轴最小，第二座标轴次小（须垂直于先前座标轴），……。此即直线拟合！

也就是说，座标轴的垂直方向，就是数据最聚合的方向。

座标轴总是穿越原点。当数据没有预先减去平均值，则座标轴不会穿越数据中心，不过仍然具备前述性质。此即“必须穿越原点的直线拟合”！进一步改变原点，则可以制作“必须穿越任意点的直线拟合”。

数学证明

首先引入先备知识：Ax²求最大值、最小值。A 是实数对称正半定方阵。当 x 长度为 1，答案是 A 的最大的、最小的特徵值暨特徵向量。

     T                 2
max x A x   subject ‖x‖ = 1
 x
     T            2
max x A x - λ (‖x‖ - 1)          Lagrange's multiplier
 x

∂     T            2
―― [ x A x - λ (‖x‖ - 1) ] = 0   derivative = 0
∂x

2 A x - 2 λ x = 0                expand (A is symmetric)

A x = λ x                        transposition

A 最大的特徵值暨特徵向量就是正解。

投影之后变异数最大、数据散开性证明：

                                 2
max ∑ ‖proj pᵢ - (1/n) ∑ proj pⱼ‖    subject to ‖v‖ = 1
 v  i    v             j   v

               2
max ∑ ‖proj pᵢ‖    subject to ‖v‖ = 1
 v  i    v

            2 
max ∑ ‖pᵢ∙v‖       subject to ‖v‖ = 1
 v  i

         T  2
max ∑ ‖pᵢ v‖       subject to ‖v‖ = 1
 v  i

      T  2
max ‖X v‖          subject to ‖v‖ = 1
 v

     T   T 
max v X X v        subject to ‖v‖ = 1
 v

     T    T
max v (X X ) v     subject to ‖v‖ = 1
 v

   T
X X 最大的特徵值暨特徵向量就是正解。

距离平方和最小、直线拟合、数据聚合性证明：

                    2
min ∑ ‖pᵢ - proj pᵢ‖     subject to ‖v‖ = 1
 v  i         v

                     2
min ∑ ‖pᵢ - (pᵢ∙v) v‖    subject to ‖v‖ = 1
 v  i

          2          2        2   2
min ∑ ‖pᵢ‖ - 2 ‖pᵢ∙v‖ + ‖pᵢ∙v‖ ‖v‖    subject to ‖v‖ = 1
 v  i

          2        2 
min ∑ ‖pᵢ‖ - ‖pᵢ∙v‖      subject to ‖v‖ = 1
 v  i

              2 
min ∑ - ‖pᵢ∙v‖           subject to ‖v‖ = 1
 v  i

            2 
max ∑ ‖pᵢ∙v‖             subject to ‖v‖ = 1
 v  i

   T
X X 最大的特徵值暨特徵向量就是正解。

Low-rank Principal Component Analysis

关于转换矩阵 A，先前讨论方阵，现在讨论矩阵。low-rank 是指矩阵的 rank 比较小，外观看起来是直条不足或者横条不足。

转换矩阵不再拥有反矩阵，但是仍拥有广义反矩阵，重新表示、分解仍是同一个问题。为了清楚说明，以下将两者皆列出。

重新表示版本：令新数据维度小于原数据维度，降低资讯量。

答案是 A = √D⁻¹Eᵀ，保留大的、捨弃小的特徵值暨特徵向量。

分解版本：令原数据维度小于新数据维度，降低资讯量。

答案是 A = E√D，保留大的、捨弃小的特徵值暨特徵向量。

如果式子再添加一个误差项，称作 Factor 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]]; 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] 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]

几何特性

除了先前提及的散开性、聚合性，又多了一个逼真性。

三、所有数据投影到座标轴空间的距离平方总和最小。

也就是说，座标轴所在位置，令数据投影之后的失真最少。

数学证明

首先引入先备知识：对角线矩阵 D，套用正规正交矩阵 Q，trace 只会减少，或者不变（当正规正交矩阵 Q 是单位矩阵 I）。

几何观点：向量们形成标准座标轴，长度不必是 1。经过旋转或翻转，则会偏离标准座标轴，使得座标值减少，总和也减少。

max tr( Q D )

令 Q = I 以得到最大值

投影前后距离平方和最小、数据逼真性证明：

(Ǎ: remove some columns from square matrix A)

             T   2               T
min ‖ X - B̌ B̌ X ‖    subject to B̌ B̌ = I

                T  T         T
min tr( (X - B̌ B̌ X)  (X - B̌ B̌ X) )        trace of inner product

         T       T   T     T   T   T
min tr( X X - 2 X B̌ B̌ X + X B̌ B̌ B̌ B̌ X )   expand

         T     T   T                       T
min tr( X X - X B̌ B̌ X )                   B̌ B̌ = I

           T   T
min tr( - X B̌ B̌ X )                       remove constant

         T   T
max tr( X B̌ B̌ X )                         multiply -1

           T   T
max tr( B̌ B̌ X X )                         cyclic property of trace

           T     T
max tr( B̌ B̌ E D E )                       eigendecomposition

         T   T
max tr( E B̌ B̌ E D )                       cyclic property of trace

    T   T
令 E B̌ B̌ E = Ǐ 以得到最大值 Ď

令 B̌ = Ě 以得到最大值 Ď

   T
X X 较大的特徵值和特徵向量就是正解。

精髓

从不同角度看数据，时聚时散。不断滚动数据，重新找一个视角，让数据看起来最散开、最聚合。注意到，是聚散，不是宽窄。

延伸应用

whitening：实施 PCA，位移、旋转、缩放。将数据范围大致调整成单位圆，方便比对数据。观念类似于 normalization。

orientation：实施 PCA，位移、旋转、略过缩放。功效是重设座标轴，摆正数据。

alignment：两堆数据各自实施 PCA，特徵值由大到小排序，特徵向量一一对应，得知数据对应关係。

embedding：实施 Low-rank PCA，捨弃小的特徵值暨特徵向量，只做投影。功效是降维、压缩，而且是误差最小的方式。

应用：Alignment

Alignment（Registration）

“对齐”。找到一个函数，让一堆数据经过函数变换之后，尽量符合另一堆数据。已知原数据、新数据，求函数。

函数具有特殊限制，否则缺乏讨论意义。以下仅讨论相似形。

     rigid：刚体。函数仅包含位移、旋转。
similarity：相似形。函数仅包含位移、旋转、整体缩放（包含翻转）。
    affine：仿射。函数仅包含位移、旋转、缩放、歪斜。

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）

一、甲堆随机抓几点，乙堆随机抓几点，点数一样多，推定它们依序一一对应。
　　以这些点实施 PA，求得对齐函数。
二、计算“甲堆套用函数”与“乙堆”的误差大小。
三、重複一二，取误差最小者。

绝圣弃智，民利百倍。

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 = √DEᵀ。保留大的、捨弃小的特徵值和特徵向量，以符合新维度，同时也让误差最小。

YᵀY = EDEᵀ = E√D√DEᵀ = (E√D)(√DEᵀ) = (√DEᵀ)ᵀ(√DEᵀ)

还有另一种手法，YᵀY 实施 Cholesky 分解（对称矩阵的 LU 分解），得到答案 Y = Lᵀ。然而降维时误差大，没人这样做。

YᵀY = LLᵀ

最后顺带一提，此算法跟“ Mesh Smoothing ”似乎有关联。

演算法（Random Projection）

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

https://en.wikipedia.org/wiki/Johnson–Lindenstrauss_lemma

线性变换的本质是旋转、镜射、伸缩座标轴，数据的相对位置不变。随便找个线性变换做为嵌入函数，都能大致保留数据的相对位置，但是会大幅改变数据的相对距离。

而随机导致分散，分散导致正交，正交导致距离不变。

Johnson-Lindenstrauss Lemma。

绝圣弃智，民利百倍。

Manifold Embedding（Manifold Learning）

更换数据所在空间，从流形换成 N 维空间。

对付流形的手法主要是：一、取样。二、建立邻居关係（不见得是平面图、三角剖分、网格）。三、邻居关係尽量保持不变。

缺点是矩阵运算时间複杂度 O(N^3) 太高。

Isomap: All Pair Shortest Path + Metric Multidimensional Scaling
Locally Linear Embedding: Weighted Average of Neighbor
Self-organizing Map: 1-layer Neural Network
Auto-encoder: Neural Network

https://en.wikipedia.org/wiki/Nonlinear_dimensionality_reduction
http://www.cs.utah.edu/~piyush/teaching/spectral_dimred.pdf
http://www.cad.zju.edu.cn/reports/流形学习.pdf
http://web.stanford.edu/class/ee378b/lect-9.pdf
http://www.cs.cmu.edu/~liuy/distlearn.htm

Independent Component Analysis

Independent Component Analysis

Y = XA。已知数据 X，求新数据 Y、求转换矩阵 A。

YᵀY = I。令新数据 Y 是正规正交：相异数据的点积是 0，相同数据的点积是 1。

一笔数据视作一个向量；两个向量的点积，得到一个值；所有两两向量的点积，排列成矩阵。

推导过程

{ Y = XA
{ YᵀY = I
given X. find A and Y.
                       _____________________________________________ 
YᵀY             = I   | XᵀX is a square matrix.                     |
(XA)ᵀXA         = I   | XᵀX is a symmetric matrix.                  |
AᵀXᵀXA          = I   | XᵀX is a real matrix.                       |
Aᵀ(XᵀX)A        = I   | XᵀX is a positive semi-definite matrix.     |
Aᵀ(EDE⁻¹)A      = I   | thus XᵀX has real non-negative eigenvalues. |
Aᵀ(E√D)(√DE⁻¹)A = I   | thus XᵀX has orthonormal eigenbasis.        |
A = E √D⁻¹            | let XᵀX = EDE⁻¹ (E⁻¹ = Eᵀ)                  |
                       ‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾ 

PCA：Y = AX，YYᵀ = I，XXᵀ = EDEᵀ，A = √D⁻¹Eᵀ。

ICA：Y = XA，YᵀY = I，XᵀX = EDEᵀ，A = E√D⁻¹。

PCA 矩阵在左，维度独立。ICA 矩阵在右，数据独立。

如同 PCA，ICA 也是预先将 X 减去平均值。X 减去平均值之后，XᵀX 和 YᵀY 变成“数据的共变异矩阵”。

Low-rank Independent Component Analysis

重新表示版本：令新数据数量小于原数据数量，降低资讯量。

答案是 A = E√D⁻¹，保留大的、捨弃小的特徵值暨特徵向量。

分解版本：令原数据数量小于新数据数量，降低资讯量。

答案是 A = √DEᵀ，保留大的、捨弃小的特徵值暨特徵向量。

PCA 与 ICA 互为转置

原数据 X、新数据 Y、转换矩阵 A 通通转置之后实施 PCA，等同于直接实施 ICA。PCA 和 ICA 对调亦如是。

Y  = A X   , Y Yᵀ = I     PCA
--------------------------------------------------
Yᵀ = AᵀXᵀ  , YᵀYᵀᵀ= I     PCA: transpose X Y A
Yᵀ = (XA)ᵀ , YᵀY  = I     PCA: transpose X Y A
Y  = X A   , YᵀY  = I     ICA

写成数学公式就是 PCA(Xᵀ)ᵀ = ICA(X)、PCA(X) = ICA(Xᵀ)ᵀ、PCA(X)ᵀ = ICA(Xᵀ)，每个式子意义都相同。

当数据形成对称矩阵，那麽 ICA 与 PCA 完全相同。然而一般不会遇到这种事，比中彩券还难。

备注：直式改成横式

因为线性变换已经规定将矩阵乘在左边，所以矩阵乘在右边挺彆扭。于是数学家就想到，把左式右式一齐转置，使得矩阵乘在左边，方便讨论；但是缺点是数据和矩阵须排列成横条，更加彆扭。顺带一提，几乎所有 ICA 文献都採用此形式。

Y = XA。Yᵀ = (XA)ᵀ。Yᵀ = AᵀXᵀ。数据须改成横条。通常重新标记成 S = AX 或者进一步移项为 X = WS。

YᵀY = I。(YᵀY)ᵀ = Iᵀ。YᵀY = I。还是跟原本一样。通常重新标记成 SSᵀ = I。

另外，除了原本的 YᵀY = I，几乎所有 ICA 文献都额外规定原数据 X 也是正规正交 XᵀX = I，接著又推导出转换矩阵 A 也是正规正交 AᵀA = I。通常重新标记成 XXᵀ = I 以及 WWᵀ = I。

横式的坏处是不符合线性变换。矩阵以直条为主，矩阵存放数据是以直条排列。横式的好处是符合计算机资料结构。阵列以横条为主，阵列储存数据是以横条排列。

应用：Separation（Under Construction!）

Separation

“分离”。找到组成数据的关键因子。

例如麦克风录到一段演奏，尝试分离出每种乐器的声音。

例如相机拍到一个场景，尝试分离出光线的来源与强度。

演算法（Independent Component Analysis）

假定数据是关键因子的加权平均值。

演算法（Wavelet Analysis）

自订特殊数据，做为关键因子。

应用：Quantization（Under Construction!）

Quantization（Vector Quantization）

“量化”。找到一个函数，让数据经过函数变换之后，保留数值的重点，删除数值的细节。

简易的量化是四捨五入、无条件捨去、无条件进入。进阶的量化是区分数量级。经典的量化是 以图表相较并没有明显差异 。

如果预先知道数据大致分布，还可以自订量化结果。亦可推广到高维度。

量化只有一笔关键因子，分离则有许多笔关键因子。

演算法（Independent Component Analysis）

实施 ICA 分解版本，只做旋转、略过缩放、略过位移。数据分解成加权平均值的格式，取权重最大者做为量化结果。

缺点是量化种类无法超过数据维度。

演算法（Cluster Analysis）

实施分群，以群集中心做为量化结果；再实施分类，以分界线决定量化结果。

https://charlesmartin14.wordpress.com/2012/10/09/spectral-clustering/
https://en.wikipedia.org/wiki/Spectral_clustering

演算法（Location-Allocation Analysis）

分群时，直线距离改成了其他指标。

Principal Component Pursuit（Under Construction!）

Principal Component Pursuit

Y = X + A。已知新数据 Y，求数据 X、求转换矩阵 A。

low-rank X。令原数据 X 的 rank 尽量小。min nuclear norm。

sparse A。令转换矩阵 A 尽量稀疏。min L1 norm。

min ‖ X ‖⁎ + α ‖ A ‖₁     (α ≥ 0)

http://www.math.umn.edu/~lerman/Meetings/CPCP.pdf

物理意义

low-rank。常见的、平淡无奇的地方。数据是其他数据的线性组合。又或者数据重新组合，类似基因演算法。数据之间很相像。

sparse。罕见的、独树一格的地方。

应用：Recommendation（Under Construction!）

Recommendation

topic model = collaborative filtering
TF-IDF
word2vector
matrix completion

Sparse Coding（Under Construction!）

Sparsity

L⁰ norm regularization。

L⁰ norm 和 L¹ norm 效果相同。

abs、log(cosh(.))。

L¹ norm 与 L² norm 差异。

古代人的作法
假设 x 的出现机率呈 sigmoid function 或者 laplace distribution
因为独立  所以所有 x 的出现机率是连乘积
用 maximum likelihood + gradient descent 来解

现代人的做法
把 sigmoid 改成 unit step
机率连乘积  取 log 变连加  这些过程通通精简掉  直接变成 1-norm 最佳化
用 least squares (绝对值误差改成平方误差以利微分) + gradient descent 来解

感觉上
只要是 sigmoid / tanh
一定可以变成 unit step 之类的东西

如果我没猜错
1. 机率学所谓的的独立事件　是稀疏（1-norm 最佳化）退化到零维的结果
2. 不含机率的稀疏（1-norm 最佳化）
   是引入机率的稀疏 sigmoid/tanh/laplace + maximum likelihood 退化之后的结果

http://blog.sciencenet.cn/blog-261330-810156.html
http://www.zhihu.com/question/26602796/answer/33431062
http://blog.csdn.net/abcjennifer/article/details/7748833

Sparse Principal Component Analysis

PCA 可以改成最佳化问题：令新旧数据的距离平方总和尽量小，限制是新座标轴必须正规正交。

採用 regularization，藉由调整参数，以控制距离远近程度、正规正交程度。regularization 所求得的新座标轴通常不是正规正交。然而座标轴歪斜一些也无妨，说不定能让新旧数据距离平方和更小，找到更拟合数据的座标轴。

             T   2               T
min ‖ X - B̌ B̌ X ‖    subject to B̌ B̌ = I

             T   2            2
min ‖ X - B̌ B̌ X ‖  + α ( ‖ B̌ ‖ - I )     (α ≥ 0)

             T   2          2
min ‖ X - B̌ B̌ X ‖  + α ‖ B̌ ‖             (α ≥ 0)

乍看很理想，不过事情还没有结束。先前提到，正确答案的左方乘上正规正交矩阵，仍是正确答案。由于答案众多，于是再增加稀疏性限制：令新座标轴的 L¹ norm 总和尽量小。

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

新座标轴 B̌越靠近标准座标轴 I，L¹ norm 就越小。而 B̌是正规正交矩阵，几何意义是旋转暨镜射。也就是说，这是令旋转角度尽量小、尽量不旋转数据。

原始论文有提供专门的最佳化演算法，请读者迳行参阅。

似乎可以视作“点集合函数”最佳化，套用 Classification 演算法解决问题。

Reconstruction Independent Component Analysis

ICA 也可以改成最佳化问题，数据转置一下，仿效 PCA。由于答案很多，于是再增加稀疏性限制：令数据投影到座标轴上的 L¹ norm 总和尽量小。

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

数据越靠近新座标轴 B̌，L¹ norm 就越小。也就是说，这是令座标轴尽量贴近数据、尽量让数据外观像个十字。

Sparse Coding（Dictionary Learning）

ICA 分解版本，不限制正规正交。

演算法是 K-SVD，请读者自行查询。

http://www.cnblogs.com/salan668/p/3555871.html

Non-negative Matrix Factorization

进一步限制数据非负数、权重非负数，符合真实世界常见情况。

例如图片处理：数据是灰阶值，非负数；权重是图片合成比重，非负数。

演算法是 Projected Gradient Descent，请读者自行查询。

https://zhuanlan.zhihu.com/p/22043930