# 写在前面

本文介绍一篇流形降维领域非常经典的文章[1]，利用SPD流形的几何关系分别构建监督和无监督两种范式的降维学习。

# 基础知识

# 仿射不变黎曼度量

SPD流形S++n\mathcal S_{++}^n研究最多的度量还是仿射不变黎曼度量。对于点P∈S++n\boldsymbol{P} \in \mathcal S_{++}^n和切向量v,w∈TPS++n\boldsymbol{v},\boldsymbol{w} \in T_{\boldsymbol{P}}{\mathcal S_{++}^n}，

⟨v,w⟩P≜⟨P−1/2vP−1/2,P−1/2wP−1/2⟩=tr(P−1vP−1w)\begin{aligned} \langle \boldsymbol{v}, \boldsymbol{w} \rangle_\boldsymbol{P} &\triangleq \langle \boldsymbol{P}^{-1/2}\boldsymbol{v}\boldsymbol{P}^{-1/2}, \boldsymbol{P}^{-1/2}\boldsymbol{w}\boldsymbol{P}^{-1/2} \rangle \\ &= \text{tr} \left( \boldsymbol{P}^{-1} \boldsymbol{v} \boldsymbol{P}^{-1} \boldsymbol{w}\right) \end{aligned}

该度量诱导出测地线距离δR(X,Y)=∥log⁡(X−1/2YX−1/2)∥F\delta_R(\boldsymbol{X},\boldsymbol{Y}) = \|\log(\boldsymbol{X}^{-1/2}\boldsymbol{Y}\boldsymbol{X}^{-1/2})\|_F具有仿射变换不变性，即δR(X,Y)=δR(AXAT,AYAT)\delta_R(\boldsymbol{X},\boldsymbol{Y}) = \delta_R(\boldsymbol{A}\boldsymbol{X}\boldsymbol{A}^T,\boldsymbol{A}\boldsymbol{Y}\boldsymbol{A}^T)。

# Bregman散度

记ζ:S++n→R\zeta : \mathcal{S}_{++}^{n} \rightarrow \mathbb{R}为严格凸且可微的函数，则Bregman矩阵散度为

dζ(X,Y)=ζ(X)−ζ(Y)−⟨∇ζ(Y),X−Y⟩,d_\zeta(\boldsymbol{X},\boldsymbol{Y}) = \zeta(\boldsymbol{X}) - \zeta(\boldsymbol{Y}) - \langle \nabla_{\zeta} (\boldsymbol{Y}) , \boldsymbol{X} - \boldsymbol{Y} \rangle \; ,

Bregman散度具有非对称、非负和正定性质。下面将Bregman散度对称化。

  • Stein散度

选择ζ(X)=−log⁡det⁡(X)\zeta(\boldsymbol{X}) = - \log \det(\boldsymbol{X})，通过Jensen-Shannon 对称化得到Stein散度

δS2(X,Y)≜12dζ(X,X+Y2)+12dζ(Y,X+Y2)=log⁡det⁡(X+Y2)−12log⁡det⁡(XY)\begin{aligned} \delta_S^2(\boldsymbol{X},\boldsymbol{Y}) &\triangleq \frac{1}{2} d_{\zeta}\left( \boldsymbol{X},\frac{\boldsymbol{X}+\boldsymbol{Y}}{2} \right) + \frac{1}{2} d_{\zeta}\left( \boldsymbol{Y},\frac{\boldsymbol{X}+\boldsymbol{Y}}{2} \right) \\ &= \log \det\bigg( \frac{\boldsymbol{X}+\boldsymbol{Y}}{2} \bigg) - \frac{1}{2} \log \det ( \boldsymbol{X}\boldsymbol{Y}) \end{aligned}

  • Jeffrey散度

选择ζ(X)=−ln⁡det⁡(X)\zeta(\boldsymbol{X}) = -\ln \det(\boldsymbol{X})，通过平均得到Jeffrey散度

δJ2(X,Y)≜12dζ(X,Y)+12dζ(Y,X)=12tr(X−1Y)−12log⁡det⁡(X−1Y)+12tr(Y−1X)−12log⁡det⁡(Y−1X)−n=12tr(X−1Y)+12tr(Y−1X)−n\begin{aligned} \delta_J^2(\boldsymbol{X},\boldsymbol{Y}) &\triangleq \frac{1}{2} d_\zeta(\boldsymbol{X},\boldsymbol{Y}) + \frac{1}{2}d_\zeta(\boldsymbol{Y},\boldsymbol{X}) \\ &= \frac{1}{2}\text{tr} (\boldsymbol{X}^{-1}\boldsymbol{Y}) - \frac{1}{2}\log \det(\boldsymbol{X}^{-1}\boldsymbol{Y}) \\ &+ \frac{1}{2}\text{tr} (\boldsymbol{Y}^{-1}\boldsymbol{X}) - \frac{1}{2}\log \det(\boldsymbol{Y}^{-1}\boldsymbol{X}) -n \\ &= \frac{1}{2}\text{tr} (\boldsymbol{X}^{-1}\boldsymbol{Y}) + \frac{1}{2}\text{tr} (\boldsymbol{Y}^{-1}\boldsymbol{X}) -n \end{aligned}

Stein散度和Jeffrey散度同样具有仿射变换不变性

δS2(X,Y)=δS2(AXAT,AYAT).δJ2(X,Y)=δJ2(AXAT,AYAT).\begin{aligned} \delta_S^2(\boldsymbol{X},\boldsymbol{Y}) & = \delta_S^2(\boldsymbol{A}\boldsymbol{X}\boldsymbol{A}^T,\boldsymbol{A}\boldsymbol{Y}\boldsymbol{A}^T). \\ \delta_J^2(\boldsymbol{X},\boldsymbol{Y}) & = \delta_J^2(\boldsymbol{A}\boldsymbol{X}\boldsymbol{A}^T,\boldsymbol{A}\boldsymbol{Y}\boldsymbol{A}^T). \end{aligned}

此外，Stein散度、Jeffrey散度和AIRM度量得到的测地线距离具有倍数关系。

# 优化框架

一般都是提出优化问题后再设计求解算法，不过本文集成了监督学习和无监督学习，但都是使用梯度下降法，所以把优化前置了，后面在模型里直接放结果。该框架可表示为如下优化框架：

min⁡W∈Rn×mf(W)s.t.WTW=Im\min_{\boldsymbol{W}\in \mathbb{R}^{n \times m}}~f(\boldsymbol{W}) \quad {\rm s.t.}\quad\boldsymbol{W}^T\boldsymbol{W} = \mathbf{I}_m

在仿射变换不变性质下，f(W)f(\boldsymbol{W})与由W\boldsymbol{W}张成的基选择无关，即对正交矩阵R∈O(m)\boldsymbol{R} \in \mathcal{O}(m)，有f(W)=f(WR)f(\boldsymbol{W}) = f(\boldsymbol{W}\boldsymbol{R})，则该问题落在了Grassmann流形G(m,n)\mathcal G(m,n)上，下面使用Grassmann优化技术中的黎曼共轭梯度法完成优化。

# 黎曼共轭梯度法

记∇W(f)\nabla_\boldsymbol{W}(f)为f(W){f}(\boldsymbol{W})关于W\boldsymbol{W}的雅可比矩阵，G(m,n)\mathcal G(m,n)上的梯度表示为

grad⁡f(W)=(In−WWT)∇W(f)\operatorname{grad}f(\boldsymbol{W}) = (\mathbf{I}_n - \boldsymbol{W}\boldsymbol{W}^T)\nabla_\boldsymbol{W}(f)

黎曼共轭梯度法包含以下三个步骤：

  1. 计算当前的梯度grad⁡f(W)\operatorname{grad}f(\boldsymbol{W})
  2. 通过平行传输先前的搜索方向并将其与grad⁡f(W)\operatorname{grad}f(\boldsymbol{W})组合来确定搜索方向HH。
  3. 在方向HH上沿WW的测地线进行直线搜索。

在停止前反复迭代这三步可得到局部极小点。

# 平行传输

上面有一步τ(H,W,V)\tau(\boldsymbol{H},\boldsymbol{W},\boldsymbol{V})表示从W\boldsymbol{W}到V\boldsymbol{V}的切向量的的并Δ\Delta行传输。

与欧式空间不同，在流形上不能通过简单的平移将切向量Δ\Delta从一点传递到另一点。在流形M\mathcal M上将Δ\Delta从W\boldsymbol{W}运输到V\boldsymbol{V}需要减去法向分量Δ⊥\Delta_{\perp}到的矢量为切向量。这种切向量的转移称为平行传输(parallel transporting)。

# 流形S++n\mathcal S_{++}^n降维

通过双线性映射实现高维流形S++n\mathcal S_{++}^n到低维流形S++m\mathcal S_{++}^m的嵌入，即对X∈S++n\boldsymbol{X}\in\mathcal S_{++}^n

fW(X)=WTXW∈S++mf_\boldsymbol{W}(\boldsymbol{X}) = \boldsymbol{W}^T\boldsymbol{X}\boldsymbol{W}\in \mathcal S_{++}^m

# 监督降维

首先需要将标签信息编码至流形点的几何关系中，下面使用类内距离和类间距离定义一个仿射关系，利用双线性映射同时最小化类间距离且最大化类内距离，完成判别分析。

  • 类内相似度gw:S++n×S++n→R+g_w: \mathcal S_{++}^n \times \mathcal S_{++}^n \to \mathbb{R}_+

    gw(Xi,Xj)={1,ifXi∈Nw(Xj)orXj∈Nw(Xi)0,otherwiseg_w(\boldsymbol{X}_i,\boldsymbol{X}_j) = \left\{ \begin{matrix} 1, & \text{if} \; \boldsymbol{X}_i \in N_w(\boldsymbol{X}_j ) \; \text{or} \; \boldsymbol{X}_j \in N_w(\boldsymbol{X}_i ) \\ 0, & \text{otherwise} \end{matrix} \right.

    Nw(Xi)N_w(\boldsymbol{X}_i)是Xi\boldsymbol{X}_i与共享yiy_i相同标签的νw\nu_w个最近邻居的集合。

  • 类间相似度gb:S++n×S++n→R+g_b: \mathcal S_{++}^n \times \mathcal S_{++}^n \to \mathbb{R}_+

    gb(Xi,Xj)={1,ifXi∈Nb(Xj)orXj∈Nb(Xi)0,otherwiseg_b(\boldsymbol{X}_i,\boldsymbol{X}_j) = \left\{ \begin{matrix} 1, & \text{if} \; \boldsymbol{X}_i \in N_b(\boldsymbol{X}_j ) \; \text{or} \; \boldsymbol{X}_j \in N_b(\boldsymbol{X}_i ) \\ 0, & \text{otherwise} \end{matrix} \right.

    Nb(Xi)N_b(\boldsymbol{X}_i)是νb\nu_b个包含Xi\boldsymbol{X}_i的最近邻居，并带有不同的标签。

  • 仿射函数a(⋅,⋅)a(\cdot,\cdot)类似于最大间隔准则(Maximum Margin Criterion)

    a(Xi,Xj)=gw(Xi,Xj)−gb(Xi,Xj)a(\boldsymbol{X}_i,\boldsymbol{X}_j) = g_w(\boldsymbol{X}_i,\boldsymbol{X}_j) - g_b(\boldsymbol{X}_i,\boldsymbol{X}_j)

  • 目标函数

    L(W)=∑i,j=1j≠ipa(Xi,Xj)δ2(WTXiW,WTXjW)L(\boldsymbol{W}) = \sum_{\substack{i,j=1 \\ j \neq i}}^p a(\boldsymbol{X}_i,\boldsymbol{X}_j) \delta^2 \left(\boldsymbol{W}^T\boldsymbol{X}_i\boldsymbol{W},\boldsymbol{W}^T\boldsymbol{X}_j\boldsymbol{W}\right)

    其中δ\delta可以选择δR\delta_R、δS\delta_S或δJ\delta_J。

  • 梯度计算

算法流程

# 无监督降维

无监督场景无法利用标签信息对样本的进行区分。下面使用最大化方差的准则来完成流形无监督降维，即

记所有样本点{Xi}i=1p\{\boldsymbol{X}_i\}_{i = 1}^{p}关于度量δ\delta的几何均值为M\boldsymbol{M}，通常选择Frechet 形式来计算

M∗≜arg⁡min⁡M∈S++n∑i=1pδ2(Xi,M)\boldsymbol{M}^\ast \triangleq \arg\min_{\boldsymbol{M} \in \mathcal S_{++}^n}\sum_{i=1}^p \delta^2 (\boldsymbol{X}_i, \boldsymbol{M})

  • AIRM度量，利用指数算子和对数算子计算黎曼(Karcher)均值
  • Stein散度，通过迭代凸凹步骤(CCCP)计算均值
  • Jeffrey散度，均值具有显式表达

最大化方差意味着所有投影点与均值投影点的距离最大化，即

L(W)=∑i=1pδ2(WTXiW,WTMW)L(\boldsymbol{W}) = \sum_{i = 1}^p \delta^2\left( \boldsymbol{W}^T\boldsymbol{X}_i\boldsymbol{W}, \boldsymbol{W}^T\boldsymbol{M}\boldsymbol{W} \right)

注意：WTMW\boldsymbol{W}^T \boldsymbol{M} \boldsymbol{W}并不是投影后的几何均值点，即

WT(arg⁡min⁡M∈S++n∑i=1pδ2(Xi,M))W≠arg⁡min⁡F∈S++m∑i=1pδ2(WTXiW,F)\boldsymbol{W}^T\left(\arg\min_{\boldsymbol{M} \in \mathcal S_{++}^n}\sum_{i=1}^p \delta^2 (\boldsymbol{X}_i, \boldsymbol{M}) \right)\boldsymbol{W}\neq \arg\min_{\boldsymbol{F} \in \mathcal S_{++}^m}\sum_{i=1}^p \delta^2(\boldsymbol{W}^T\boldsymbol{X}_i\boldsymbol{W},\boldsymbol{F})

# Log-Euclidean度量与监督降维

除了AIRM和两个Bregman散度，Log-Euclidean度量也是流形S++n\mathcal S_{++}^n上常用的度量

δlE(X,Y)=∥log⁡(X)−log⁡(Y)∥F\delta_{lE}(\boldsymbol{X},\boldsymbol{Y}) = \|\log(\boldsymbol{X}) - \log(\boldsymbol{Y})\|_F

监督降维的度量δ\delta改为Log-Euclidean度量，得到如下优化目标

L(W)=∑i,j=1pa(Xi,Xj)∥log⁡(WTXiW)−log⁡(WTXjW)∥F2L(\boldsymbol{W}) = \sum_{i,j=1}^p a(\boldsymbol{X}_i,\boldsymbol{X}_j) \Big \| \log(\boldsymbol{W}^T\boldsymbol{X}_i \boldsymbol{W}) -\log(\boldsymbol{W}^T \boldsymbol{X}_j \boldsymbol{W}) \Big \|_F^2

通过将log⁡(WTXW)\log(\boldsymbol{W}^T \boldsymbol{X} \boldsymbol{W})近似为WTlog⁡(X)W\boldsymbol{W}^T \log(\boldsymbol{X}) \boldsymbol{W}，目标函数可改写为

L(W)=∑i,j=1pa(Xi,Xj)∥WTlog⁡(Xi)W−WTlog⁡(Xj)W∥F2=tr(WTF(W)W)\begin{aligned} L(\boldsymbol{W}) &= \sum_{i,j=1}^p a(\boldsymbol{X}_i,\boldsymbol{X}_j) \Big \| \boldsymbol{W}^T \log(\boldsymbol{X}_i)\boldsymbol{W} - \boldsymbol{W}^T \log(\boldsymbol{X}_j)\boldsymbol{W} \Big \|_F^2\\ &= \text{tr}{ \Big( \boldsymbol{W}^T \boldsymbol{F}(\boldsymbol{W}) \boldsymbol{W} \Big)} \end{aligned}

其中

F(W)=∑i,j=1pa(Xi,Xj)(log⁡(Xi)−log⁡(Xj))WWT×(log⁡(Xi)−log⁡(Xj))\boldsymbol{F}(\boldsymbol{W}) = \sum_{i,j=1}^p a(\boldsymbol{X}_i,\boldsymbol{X}_j) \Big( \log(\boldsymbol{X}_i) - \log(\boldsymbol{X}_j) \Big) \boldsymbol{W}\boldsymbol{W}^T \times \Big(\log(\boldsymbol{X}_i) - \log(\boldsymbol{X}_j) \Big)

则优化问题可改写为特征分解问题

min⁡Wtr(WTF(W)W) s.t. WTW=Im.\min_{\boldsymbol{W}} \text{tr}{\Big( \boldsymbol{W}^T \boldsymbol{F}(\boldsymbol{W}) \boldsymbol{W} \Big)} ~\text{s.t.}~ \boldsymbol{W}^T\boldsymbol{W} = \mathbf{I}_m\;.

采样二阶段迭代方式更新W\boldsymbol{W}和F(W)\boldsymbol{F}(\boldsymbol{W})，算法步骤如下：

  • Log-Euclidean度量的监督模型求解算法对比发现，特征分解法比黎曼共轭梯度法快10倍。

  • 把Bregman散度的函数选为ζ(X)=∥X∥F2\zeta(\boldsymbol{X}) = \|\boldsymbol{X}\|_F^2，即使用∥⋅∥F2\|\cdot\|_F^2作为相似性度量，可以得到类似的特征分解问题。

# 小结

后面还有一些东西就不细说了。作者M. Harandi[2]、M. Salzmann[3]和R. Hartley[4]做了很多机器学习与几何的结合。

# References

  1. M. Harandi, M. Salzmann and R. Hartley, Dimensionality Reduction on SPD Manifolds: The Emergence of Geometry-Aware Methods, in IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 40, no. 1, pp. 48-62, 1 Jan. 2018, doi: 10.1109/TPAMI.2017.2655048. ↩︎

  2. Mehrtash Harandi homepage: https://sites.google.com/site/mehrtashharandi/ (opens new window) ↩︎

  3. Mathieu Salzmann homepage: https://people.epfl.ch/mathieu.salzmann (opens new window) ↩︎

  4. Richard Hartley google scholar: https://scholar.google.com/citations?hl=zh-CN&user=cHia5p0AAAAJ&view_op=list_works&sortby=pubdate (opens new window) ↩︎

基于Frechet均值的Grassmann判别分析