降维系列之 LTSA 局部切空间排列

前面写的PCA，LE，LDA，LLE都是以前就比较熟悉的东西，从这篇开始写的都是之前不熟悉的甚至都不知道名字的算法，然而都还很经典。疫情期间在家里看看原文，学习学习，既是算法总结又是读论文笔记。这篇来写LTSA局部切空间排列。本篇符号尽量与原文保持一致，与前面几篇有所不同。

主要思路

LTSA（Local Tangent Space Alignment）的基本思路是用样本点的近邻区域的切空间来表示局部几何结构，然后对局部切空间进行重新排列得到非线性流形的用自然参数刻画的低维表示，是经典的流形学习与降维算法。

假设一个$d$维流形嵌于$m$维空间中（$d<m$），$m$维空间是一个包含噪声的高维空间。给定样本集合$X=[x_1,\cdots ,x_N],x_i\in {\mathbb{R}}^{m\times 1}$分布于这个含噪声的$m$维空间中，文章认为：
$$x_i=f({\tau }_i)+{ϵ}_i,i=1,\cdots ,N$$

${\tau }_i\in {\mathbb{R}}^{d\times 1}$是$x_i$的本征表示，$f$是一个映射函数，${ϵ}_i$表示噪声。

线性情况

先来看线性情况。线性时$f$是一个投影矩阵：
$$\begin{gathered} f({\tau }_i)=x_i=c+U{\tau }_i+{ϵ}_i \\ X=c{e}^T+UT+E \\ T=[{\tau }_1,\cdots ,{\tau }_N],E=[{ϵ}_1,\cdots ,{ϵ}_N] \end{gathered}$$

$c\in {\mathbb{R}}^{m\times 1}$表示常量参数，是个偏置项，$e$是全为1的列向量。线性降维与线性流形学习的目标就是最小化重建误差：
$$\min \Vert E\Vert =\underset{c,U,T}{\min }\Vert X-(c{e}^T+UT){\Vert }_F$$

这是个$F$范数的形式，希望噪声能够最小。文章说，基于观察认为：

• 如果能从$E$中移除$E$的行均值，就可以使其范数更小，这要求$c=Xe/N=\overline{x}$。其实就是对$X$做中心化的意思。
• 然后希望上式最小化就等价于希望$UT$是最接近于$X-\overline{x}{e}^T$的秩为$d$的矩阵，这可以用$X-\overline{x}{e}^T$的SVD分解来求解。这个确实如此，这是SVD分解的低秩近似性质。

SVD分解有的低秩近似性质：即给定秩为$r$的矩阵$A$并求出其SVD分解$A=U\Sigma V^T$，欲求一个秩为$k(k\le r)$的矩阵$\tilde{A}$，并使得$\tilde{A}与$A$最接近（差的F范数最小）：
$$\underset{\tilde{A}}{\min }\Vert A-\tilde{A}{\Vert }_F$$ 则$\tilde{A}$可以通过仅保留$A$的前$k$大奇异值得到。即$\tilde{A}=U_k{\Sigma }_k{V}_k^T$。

所以的话如果：
$$X-\overline{x}{e}^T=Q\Sigma V^T$$

则
$$UT=Q_d{\Sigma }_d{V}_d^T$$

$U$的最优解$U^{\ast }=Q_d$。因此：
$$\begin{gathered} f(\tau )=\overline{x}+U^{\ast }\tau =x \\ {U}^{\ast }T=X-\overline{x}{e}^T \\ T=(U^{\ast }{)}^T(X-\overline{x}{e}^T)={\Sigma }_d{V}_d^T \end{gathered}$$

因为$U^{\ast }$是正交矩阵，其逆矩阵就是其对称矩阵。这其实就是PCA，只是用SVD分解和F范数来解，思路也不一样，而且这告诉我们PCA降完维后的结果就是${\Sigma }_d{V}_d^T$.

继续读文章。$f$函数并不是唯一的，因为还可以被重新参数化。如果有$\tilde{\tau }$满足$\tau =P\tilde{\tau }$，则$U$^*应该变成$U^{\ast }P$。什么意思呢，意思就是：
$$f(\tau )=\overline{x}+U^{\ast }\tau =f(\tau )=\overline{x}+U^{\ast }P{P}^{-1}\tau$$

另外如果我们限制$T$是标准正交的，即$T{T}^T=I$，则应取$T=V_d^T$，而且$f$应变为：
$$f(\tau )=\overline{x}+U^{\ast }{\Sigma }_d\tau$$

这两条说的都是PCA的一些变化了。大致如此。妙的地方就在于，如果这样处理PCA的话，不需要$f$函数的参与就能完成降维了，与我们之前写的那篇PCA不同。

非线性情况

非线性的情况更复杂。一般而言，全局的非线性结构来自于局部的线性分析和排列。然后引出LTSA。对于一个给定的样本点，LTSA使用其近邻区域来构建局部切空间来表征局部几何结构，局部切空间提供了非线性流形的局部几何结构的低维线性估计，通过局部切空间来保护近邻区域中的样本点的局部坐标。然后局部切坐标通过不同的局部仿射转换在低维空间重新排列，以获得更好的全局坐标系统。

仍然是假设$d$维流形通过未知函数$f$嵌于$m$维空间中，$d<m$。给定$N$个从无噪声模型中得来的$m$维样本$X=[x_1,\cdots ,x_N],x_i\in {\mathbb{R}}^{m\times 1}$，有
$$x_i=f({\tau }_i),i=1,\cdots ,N$$

${\tau }_i\in {\mathbb{R}}^{d\times 1}$是$x_i$降完维后的结果。非线性降维的目标就是从${\tau }_i$对应的$x_i$来重构${\tau }_i$而不显式地构建$f$函数。假设$f$足够光滑，在一个给定的$\tau$处做Taylor展开：
$$f(\overline{\tau })=f(\tau )+J_f(\tau )\cdot (\overline{\tau }-\tau )+O(\Vert \overline{\tau }-\tau {\Vert }^2)$$

这里$J_f(\tau )\in {\mathbb{R}}^{m\times d}$是$f$在$\tau$处的Jacobi矩阵：
$$f(\tau )=\left[\begin{array}{c}{f}_1(\tau )\\ ⋮\\ f_m(\tau )\end{array}\right],then{J}_f(\tau )=\left[\begin{array}{ccc}\partial f_1/\partial {\tau }_1& \cdots & \partial f_1/\partial {\tau }_d\\ ⋮& ⋮& ⋮\\ \partial f_m/\partial {\tau }_1& \cdots & \partial f_m/\partial {\tau }_d\end{array}\right]$$

$f$在$\tau$处的切空间${\tau }_{\tau }$是由$J_f(\tau )$的$d$个列向量做基底的生成空间，维度最高为$d$，${\tau }_{\tau }=span(J_f(\tau ))$。向量$\tau -\overline{\tau }$是$f(\tau )$在仿射空间$f(\tau )+{\tau }_{\tau }$的坐标（这一句我属实没懂）。因为不知道$f$，所以$J_f(\tau )$也没法求。如果$Q_{\tau }$是${\tau }_{\tau }$的一个标准正交基矩阵，如果基于$Q_{\tau }$可以知道$J_f(\tau )$，可以写作：
$$J_f(\tau )(\overline{\tau }-\tau )=Q_{\tau }{\theta }^{\ast }$$

然后
$${\theta }_{\tau }^{\ast }=Q_{\tau }^T{J}_f(\tau )(\overline{\tau }-\tau )\equiv P_{\tau }(\overline{\tau }-\tau )$$

从$\tau$到${\theta }_{\tau }^{\ast }$的映射表示局部仿射变换，这个仿射变换同样是未知的，因为$f$未知。然而向量${\theta }_{\tau }^{\ast }$有一个近似值${\theta }_{\tau }$，其可以正交地将$f(\overline{\tau })-f(\tau )$投影到${\tau }_{\tau }$：
$${\theta }_{\tau }\equiv Q_{\tau }^T(f(\overline{\tau })-f(\tau ))={\theta }_{\tau }^{\ast }+O(\Vert \overline{\tau }-\tau {\Vert }^2)$$

假设$Q_{\tau }^T$在每一个$\tau$处都是已知的，忽略二阶项，则全局坐标$\tau$应满足：
$$\int d\tau {\int }_{\Omega (\tau )}\Vert P_{\tau }(\overline{\tau }-\tau )-{\theta }_{\tau }\Vert d\overline{\tau }\approx 0$$

$\Omega (\tau )$是$\tau$的近邻区域的定义。因此，一种自然的优化目标就是寻找合适的$\tau$和$P_{\tau }$以最小化以下误差函数：
$$\int d\tau {\int }_{\Omega (\tau )}\Vert P_{\tau }(\overline{\tau }-\tau )-{\theta }_{\tau }\Vert d\overline{\tau }$$

这表示的是降维问题的非线性排列方法。

线性排列方法可以有如下思路。如果$J_f(\tau )$是列满秩的，矩阵$P_{\tau }$就应该是非奇异的并且：
$$\overline{\tau }-\tau \approx P_{\tau }^{-1}{\theta }_{\tau }\equiv L_{\tau }{\theta }_{\tau }$$

这时候应该寻找全局坐标$\tau$和局部仿射变换$L_{\tau }$来最小化下式：
$$\int d\tau {\int }_{\Omega (\tau )}\Vert (\overline{\tau }-\tau )-{\theta }_{\tau }\Vert d\overline{\tau }$$

如果不是列满秩，那就太复杂了，不讨论。

局部切空间排列LTSA

给定具有潜在非线性流形结构的包含噪声的样本集$X=[x_1,\cdots ,x_N]$：
$$x_i=f({\tau }_i)+{ϵ}_i,i=1,\cdots ,N$$

令$X_i=[x_{i1},\cdots ,x_{ik}]$是$x_i$用欧氏距离度量下的$k$近邻（包含$x_i$自身）。为$X_i$中的样本计算最佳的$d$维近似仿射空间：
$$\underset{x,\Theta ,Q}{\min }=\sum _{j=1}^k\Vert x_{ij}-(x+Q{\theta }_j){\Vert }_2^2=\underset{x,\Theta ,Q}{\min }\Vert X_i-(x{e}^T+Q\Theta ){\Vert }_F^2$$

这里$Q$是$d$列的标准正交矩阵，$\Theta =[{\theta }_1,\cdots ,{\theta }_k]$。注意，这里我写$F$范数平方的地方原文写的是二范数的平方，我觉得不太对劲，所以写成了$F$范数。
这个问题在限行情况那部分已经解过了，最优的$x$值应取${\overline{x}}_i$，最优的$Q$值$Q_i$应该取$X_i(I-e{e}^T/k)$（其实就是对$X_i$做了中心化的结果）前$d$大的奇异值对应的左奇异向量。然后$\Theta$应该由${\Theta }_i$组成，且：
$${\Theta }_i=Q_i^T{X}_i(I-\frac{1}{k}e{e}^T)=[{\theta }_1^{(i)},\cdots ,{\theta }_k^{(i)}],{\theta }_j^{(i)}=Q_i^T(x_{ij}-{\overline{x}}_i)$$

意思就是在$x_i$的包含自身的$k$个近邻组成的$X_i$范围内做了个PCA，只是不降维，维度仍为$m$，做完PCA后的$X_i$为${\Theta }_i$，单个样本为${\theta }_j^{(i)}$.

然后有：
$$\begin{gathered} x_{ij}={\overline{x}}_i+Q_i{\theta }_j^{(i)}+{\xi }_j^{(i)} \\ {\xi }_j^{(i)}=x_{ij}-{\overline{x}}_i-Q_i{Q}_i^T(x_{ij}-{\overline{x}}_i)=(I-Q_i{Q}_i^T)(x_{ij}-\overline{x}) \end{gathered}$$

这里的${\xi }_j^{(i)}$就是重构误差。${\theta }_j^{(i)}$可以表征局部结构，现在基于${\theta }_j^{(i)}$来得到${\tau }_i,i=1,\cdots ,N$. LTSA希望${\tau }_{ij}$满足如下形式：
$${\tau }_{ij}={\overline{\tau }}_i+L_i{\theta }_j^{(i)}+{ϵ}_j^{(i)},i=1,\cdots ,N$$

其中${\overline{\tau }}_i$是$k$个${\tau }_{ij}$的均值，$L_i$是未知的仿射变换，起到一个排列的作用。上式的矩阵形式表示为：
$$T_i=\frac{1}{k}{T}_{i}e{e}^T+L_i{\Theta }_i+E_i$$

$T_i=[{\tau }_{i1},\cdots ,{\tau }_{ik}],E_i=[{ϵ}_1^{(i)},\cdots ,{ϵ}_k^{(i)}]$。重构残差$E_i$为：
$$E_i=T_i(I-\frac{1}{k}e{e}^T)-L_i{\Theta }_i$$

为了在低维空间中保护尽可能多的局部几何结构，LTSA希望降维后得到的样本表示${\tau }_i$以及局部仿射变换$L_i$，可以最小化重构残差${ϵ}_j^{(i)}$：
$$\sum _i\Vert E_i{\Vert }^2\equiv \sum _i\Vert T_i(I-\frac{1}{k}e{e}^T)-L_i{\Theta }_i{\Vert }^2=\min \text{\hspace{1em}(1)}$$

显然，能够最小化误差$\Vert E_i{\Vert }_F$的排列矩阵$L_i$的最优解应为：
$$L_i=T_i(I-\frac{1}{k}e{e}^T)=T_i{\Theta }_i^{+}，and{E}_i=T_i(I-\frac{1}{k}e{e}^T)(I-{\Theta }_i^{+}{\Theta }_i)$$

这里的${\Theta }_i^{+}$是${\Theta }_i$的Moor-Penrose广义逆矩阵。$T=[{\tau }_1,\cdots ,{\tau }_N]$，令$S_i$为满足$T{S}_i=T_i$的0-1选择矩阵（比方说要筛选第1，3个样本，$S_i$的第1行的第1个数就为1，第2列的第3个数就为1，其余位置全为0），注意前面定义了$T_i=[{\tau }_{i1},\cdots ,{\tau }_{ik}],E_i=[{ϵ}_1^{(i)},\cdots ,{ϵ}_k^{(i)}]$. 现在需要寻找$T$使得总体所有样本的重构误差最小化：
$$\underset{T}{\min }\Vert E_i{\Vert }_F^2=\underset{T}{\min }\Vert TSW{\Vert }_F^2$$

其中$S=[S_1,\cdots ,S_N]$，$W=diag(W_1,\cdots ,W_N)$且：
$$W_i=(I-\frac{1}{k}e{e}^T)(I-{\Theta }_i^{+}{\Theta }_i)$$

为了使得$T$有唯一解，LTSA限制$T{T}^T=I_d$，$I_d$是$d$维的单位矩阵。

求解

现在基本搞明白了，给定$X$，LTSA先给每个样本找欧氏距离度量下的$k$个近邻（包含自身），为每个样本$x_i$构成一个包含自身的近邻区域$X_i$，然后在$X_i$这个区域内做PCA但却不降维，然后$X_i$通过PCA变成了${\Theta }_i$，其中的$x_{ij}$变成了${\theta }_j^{(i)}$。然后又认为降维的结果$T_i$与$X_i$之间有线性关系，希望二者之间的误差最小，仿射关系$L_i$被$T_i$表示了，残差$E_i$也用$T_i$表示了，所以变成了一个非线性的方法。局部PCA就是所谓的局部切空间，后面这个非线性降维就是所谓的排列，最终变成如下形式的问题：

$$\begin{gathered} \underset{T}{\min }\Vert TSW{\Vert }_F^2 \\ s.t.T{T}^T=I \end{gathered}$$

求解思路也很清晰，F范数的平方化为向量二范数的平方之和即可。这里重新定义$T_i$为$T$的第$i$行，与上一部分有所不同，则：
$$\begin{gathered} \Vert TSW{\Vert }_F^2 \\ =\sum _{i=1}^d\Vert T_{i}SW{\Vert }_2^2 \\ =\sum _{i=1}^d(T_{i}SW)(T_{i}SW{)}^T \\ =\sum _{i=1}^d{T}_i(SW{W}^{T}S)T_i^T \\ =tr(TSW{W}^T{S}^T{T}^T) \\ s.t.T{T}^T=I \end{gathered}$$

然后用经典的Lagrangian乘子法：
$$L(T,\lambda )=tr(TSW{W}^T{S}^T{T}^T)+tr(\lambda (I-T{T}^T))$$

如同之前几篇博客，这里的$\lambda$为对角矩阵。然后求导并令导数为0：
$$\begin{gathered} \frac{\partial }{\partial X}L(T,\lambda )=2T(SW{W}^T{S}^T)-2\lambda T=0 \\ T(SW{W}^T{S}^T)=\lambda T \\ (SW{W}^T{S}^T)T^T=T^T\lambda \\ \frac{\partial }{\partial \lambda }L(T,\lambda )=I-T{T}^T=0 \\ T{T}^T=I \end{gathered}$$

然后$(SW{W}^T{S}^T)T^T=T^T\lambda$就是告诉我们，$T^T$的每一列，即$T$的每一行$T_i$都是$SW{W}^T{S}^T$的特征向量，对应的特征值分布在$\lambda$的对角线的对应位置上。因为是个必要条件，再代回去：
$$tr(TSW{W}^T{S}^T{T}^T)=tr(T{T}^T\lambda )=tr(\lambda )$$

这告诉我们，原始最小化就等价于$tr(\lambda )$最小化，因此要选$SW{W}^T{S}^T$最小的前$d$个特征值对应的特征向量来构成$T$.

这里还有一个问题，就是全1列向量$e$也是$SW{W}^T{S}^T$的特征向量。首先$S$的每一列只有一个数字是1，其余都为0，因此$S^{T}e=e$.
$$SW{W}^T{S}^{T}e=SW{W}^{T}e$$

然后根据前面定义的：$W=diag(W_1,\cdots ,W_N),W_i=(I-\frac{1}{k}e{e}^T)(I-{\Theta }_i^{+}{\Theta }_i)$：
$$\begin{gathered} W_i^{T}e=(I-{\Theta }_i^{+}\Theta {)}^T(I-\frac{1}{k}e{e}^T)e \\ =(I-{\Theta }_i^{+}\Theta {)}^T(e-\frac{1}{k}ek) \\ =0 \end{gathered}$$

因此
$$SW{W}^T{S}^{T}e=SW{W}^{T}e=0\cdot e$$

即无论对什么样的$X$，$SW{W}^T{S}^T$都有一个特征值为0，对应的特征向量为$e$，这显然不是我们想要的，因此要选$SW{W}^T{S}^T$前$d$小的非0特征值对应的特征向量构成$T$。

到这里其实LTSA的定义和求解就结束了，但是论文还没有结束，还讨论了很多误差分析之类的东西，我目前并不感兴趣，也就不看了写了。

关于其他

根据前面讨论非线性情况时候得到的结论：
$$\begin{gathered} f(\overline{\tau })=f(\tau )+J_f(\tau )\cdot (\overline{\tau }-\tau )+O(\Vert \overline{\tau }-\tau {\Vert }^2) \\ f(\overline{\tau })-f(\tau )=J_f(\tau )\cdot (\overline{\tau }-\tau )+O(\Vert \overline{\tau }-\tau {\Vert }^2) \end{gathered}$$

可知，如果去掉二阶项，则有：
$$X_i(I-e{e}^T/k)\approx J_f^i{T}_i(I-e{e}^T/k)$$

0-1选择矩阵$S_i,S$的定义和前面一样，我们希望寻找合适的$J_f^i$使得二者更接近：
$$\underset{J,T}{\min }\sum _{i=1}^N\Vert (X-J_f^{(i)})S_i(I-e{e}^T/k){\Vert }_F^2=\underset{J,T}{\min }E(J,T)$$

其中$J=[J_f^{(1)},\cdots ,J_f^{(N)}]$. 这个问题可以用交替最小二乘法来解：固定$J$，通过调节$T$使得$E$最小；然后固定$T$，通过调节$J$使$E$最小。$T$的初值可以取LTSA解出来的$T$。这是另一篇文章的内容。

另一方面，公式$(1)$约束$Y{Y}^T$是为了让该问题“well-posed”，“适定”。因为如果取$T_i,L_i$都为0也有更小的解，但是这不是我们想要的。而且$Y{Y}^T$这个约束是众多可以避免这个问题的约束之一。