张量分解算法

基础知识

张量

在介绍张量分解算法（Tensor Decomposition Algorithms）之前，我们先介绍一下什么是张量（tensor）。张量是一个多维的数据存储形式（多维数组），数据的的维度被称为张量的阶。它可以看成是向量和矩阵在多维空间中的推广，更正式地说，一个N维或N阶张量是N个向量空间的张量积的一个元素，每个向量空间都有自己的坐标系。其中零阶张量是标量（scalar），一阶张量是一个向量（vector），二阶的张量是矩阵（matrix），三阶及三阶以上的就是我们通常所说的张量，也叫高阶张量。下图是一个三阶张量$\mathcal{X}\in {\mathbb{R}}^{I\times J\times K}$，有三个指标。
应该注意的是，这里的张量不是物理和工程中的张量(如应力张量)，后者在数学中通常被称为张量场。

秩一张量

若一个$N$阶张量$\mathcal{X}\in {\mathbb{R}}^{I_1\times I_2\times \cdots \times I_N}$能够被表示为$N$个向量的外积，则称$\mathcal{X}$为秩一张量，即$$\mathcal{X}={\mathbf{a}}^{(1)}⨂{\mathbf{a}}^{(2)}\cdots ⨂{\mathbf{a}}^{(\mathbf{N})}$$这意味着张量的每个元素都是对应的向量元素的乘积：$$x_{i_1{i}_2\cdots i_N}=a_{i_1}^{(1)}{a}_{i_2}^{(2)}\cdots a_{i_N}^{(N)}$$其中$1\le i_n\le I_n$ 。下图是一个三阶的秩一张量$\mathcal{X}=\mathbf{a}⨂\mathbf{b}⨂\mathbf{c}$
张量$\mathcal{X}$的元素$(i,j,k)$由$x_{ijk}=a_i{b}_j{c}_k$给出。

张量分解的意义

张量分解可以看作是矩阵奇异值分解(SVD)和主成分分析(PCA)，即矩阵分解的高阶泛化。它能够解决维度灾难问题，而传统的方法（例如ICA、PCA、SVD和NMF）对于维数比较高的数据，一般将数据展成二维的数据形式（矩阵），从而进行降维处理、缺失数据填补（或者说成“稀疏数据填补”）和隐性关系挖掘。但是这种方式会使数据的结构信息丢失，而采用张量对数据进行存储，能够保留数据的结构信息，因此近些年在图像处理以及CV等领域得到了一些广泛的应用。

张量分解中常见的两种分解是CP分解(Canonical Polyadic Decomposition，CPD)和Tucker分解(Tucker Decomposition)，下面将重点介绍这两种分解。

张量分解算法

CP分解

介绍

1927年Hitchcock基于秩一张量的定义，首次提出了将张量拆成有限个秩为1的张量相加的形式。后来在1970年Carrol和Chang将这个概念以CANDECOMP (canonical decomposition，标准分解)的形式在引入到心理测验学当中，与此同时 Harshman也提出了PARAFAC(parallel factors，平行因素) 的形式，该方法才得以流行，所以今天我们称这个将张量分解为有限个秩一张量和的形式的分解方法为CP分解，即是CANDECOMP/PARAFAC分解的缩写。

简单来说，CP分解就是将一个张量分解成多个秩为一的张量之和的形式。例如，一个三维张量$\mathcal{X}\in {\mathbb{R}}^{I\times J\times K}$分解成$R$个秩为一的张量和的形式为$$\mathcal{X}\approx \sum _{r=1}^R{a}_r⨂b_r⨂\cdots c_r$$其中$R$是正整数，$a_r\in {\mathbb{R}}^I,b_r\in {\mathbb{R}}^J,c_r\in {\mathbb{R}}^K(r=1,2,\dots ,R)$,$\mathcal{X}$中的元素为$$x_{ijk}\approx \sum _{r=1}^R{a}_{ir}{b}_{jr}{c}_{kr},i=1,2,\dots I,j=1,2,\dots J,k=1,2,\dots K.$$它的CP分解示意图如下图所示

同时，我们也可以将张量$\mathcal{X}$矩阵化（利用Khatri–Rao product进行转换，符号$⨀$,其实就是列元素匹配），并引入了权重向量$\lambda \in {\mathbb{R}}^R$将A，B，C的列标准化长度为1，此时的CP分解模型可以简单地表示为$$\mathcal{X}\approx \sum _{r=1}^R{\lambda }_r{a}_r⨂b_r⨂c_r=\left[\left[\lambda ;A,B,C\right]\right]$$
对于一般的$N$阶张量$\mathcal{X}\in {\mathbb{R}}^{I_1\times I_2\times \cdots \times I_N}$分解成$R$个秩为一的张量和的形式为$$\mathcal{X}\approx \left[\left[\lambda ;A^{(1)},A^{(2)},\dots ,A^{(N)}\right]\right]\equiv \sum _{r=1}^R{\lambda }_r{a}_r^{(1)}⨂a_r^{(2)}⨂\cdots a_r^{(N)}$$通常情况下，$a_r^{(n)}$是一个单位向量且$A^{(n)}\in {\mathbb{R}}^{I_n\times R}(n=1,2,\dots ,N)$，$A^{(n)}=\left[a_1^n{a}_2^n\cdots a_R^n\right]，\Lambda =diag（\lambda ）$,则上式可以写成$$\mathcal{X}=\Lambda {\times }_1{A}^{(1)}{\times }_2{A}^{(2)}\cdots {\times }_N{A}^{(N)}$$矩阵表达形式为$$X_{(n)}\approx A^{(n)}\Lambda (A^{(N)}⨀\cdots A^{(n+1)}⨀A^{(n-1)}\cdots ⨀A^1{)}^T$$

CP秩与CP分解求解

张量$\mathcal{X}$的秩是指，用秩一张量之和来精确表示所需要的秩一张量的最少个数R的最小值为张量的秩，记作$Rank(\mathcal{X})=R$，这样的CP分解也称为张量的秩分解。 R a n k c p ( X ) = m i n { R ∣ X = ∑ r = 1 R X ( r ) , X ( r ) 是 秩 张 量 } Rank_{cp}(\mathcal{X})=min\{R\mid\mathcal{X}=\sum \limits_{r=1}^{R}\mathcal{X}^{(r)},\mathcal{X}^{(r)}是秩张量\} Rankcp​(X)=min{R∣X=r=1∑R​X(r),X(r)是秩张量}
显然，这样定义的张量秩和矩阵的秩的定义十分相似，换句话说这种定义可以退化为矩阵的秩（即对于矩阵而言，秩的定义与普通的矩阵秩的定义等价）。但是在代数性质上，张量的秩与矩阵秩相比则表现出很多不一样的特性。张量秩的个数求解是一个NP问题，研究者只能通过一些数值方法近似的找到数个该张量的近似分解，从而估计该张量的秩。那么我们怎么去优化所有参数，使得分解后跟分解前尽可能的相似呢？

通常我们通过迭代的方法对$R$从1开始遍历直到找到一个合适的解。当数据无噪声时，重构误差为0所对应的解即为CP分解的解，当数据有噪声的情况下，可以通过CORCONDIA算法估计$R$。当分解的秩一张量的个数确定下来之后，可以通过交替最小二乘方法（the Alternating Least Squares,ALS）对CP分解进行进行求解。

对于$N$阶张量$\mathcal{X}\in {\mathbb{R}}^{I_1\times I_2\times \cdots \times I_N}$完整的交替最小二乘法优化CP分解参数的过程如下图
其中，它假设分解出秩一张量的个数$R$是已知的。
在一些应用中，为了使得CP分解更加的鲁棒和精确，可以在分解出的因子上加上一些先验知识即约束。比如说平滑约束(smooth)、正交约束、非负约束(nonegative) 、稀疏约束(sparsity)等。

Tucker 分解

介绍

Tucker分解最早由Tucker在1963年提出的，Levin和Tucker在随后的文章中对其进行了改进。
张量的Tucker分解是指，将一个张量表示成一个核张量（Core Tensor）沿着每一个模式乘上一个矩阵的形式，下面以一个三阶$\mathcal{X}\in {\mathbb{R}}^{I\times J\times K}$示意图为例

上图中$A\in {\mathbb{R}}^{I\times P}$，$B\in {\mathbb{R}}^{J\times Q}$,$C\in {\mathbb{R}}^{K\times R}$是三个因子矩阵（通常是正交的），$\mathcal{G}\in {\mathbb{R}}^{P\times Q\times R}$是核张量。具体公式如下$$\mathcal{X}=\mathcal{G}{\times }_{1}A{\times }_{2}B{\times }_{3}C=\sum _{p=1}^P\sum _{q=1}^Q\sum _{r=1}^R{g}_{pqr}{a}_p⨂b_q⨂c_r=\left[\left[\mathcal{G};A,B,C\right]\right]$$
在元素级别上，可表示为$$x_{ijk}\approx \sum _{p=1}^P\sum _{q=1}^Q\sum _{r=1}^R{g}_{pqr}{a}_{ip}{b}_{jq}{c}_{kr},i=1,2,\dots I,j=1,2,\dots J,k=1,2,\dots K.$$
这里$P、Q、R$是分量的数量(即因子矩阵$A、B$和$C$中的列向量。如果$P、Q、R$小于$I、J、K，$则核心张量$\mathcal{G}$可以看作是$\mathcal{X}$的压缩版本，在某些情况下，分解后的张量的存储显著小于原始张量。
虽然Tucker模型是在三维数组下引入的，但它可以推广到$N$阶张量：$$\mathcal{X}=\mathcal{G}{\times }_1{A}^{(1)}{\times }_2{A}^{(2)}\cdots {\times }_{(N)}{A}^{(N)}=\left[\left[\mathcal{G};A^{(1)},A^{(2)},\cdots ,A^{(N)}\right]\right]$$
其矩阵形式为$$X_{(n)}=A^{(n)}{G}_{(n)}(A^{(N)}\otimes \cdots \otimes A^{(n+1)}\otimes A^{(n-1)}\cdots \otimes A^{(1)}{)}^T$$
同时由图可以看出Tucker分解实质上是一种高阶的主成分分析（PCA，每个mode上的因子矩阵称为张量在每个mode上的基矩阵或者是主成分），因此Tucker 分解又称为高阶PCA，高阶SVD等。实际上，当核张量是对角化的且$P=Q=R$时，CP分解便是Tucker分解的一种特殊形式且CP分解。

Tucker秩与Tucker分解求解

对于$N$阶张量$\mathcal{X}\in {\mathbb{R}}^{I_1\times I_2\times \cdots \times I_N}$，定义张量的$n$-秩为其按照第$n$个模式矩阵化展开后得到的矩阵$X_{(n)}$的秩，表示为$Ran{k}_{(n)}(X)$。若令$R_n=Ran{k}_n(X)，n=1,2,\dots ,N$，那么我们称张量$\mathcal{X}$为 R a n k − { R 1 , R 2 , … , R N } Rank-\{R_1,R_2,\ldots,R_N\} Rank−{R1​,R2​,…,RN​}张量，同时张量$\mathcal{X}$的Tucker秩可定义为 R a n k T u c k e r ( X ) = m a x { R n ， n = 1 , 2 , … , N } Rank_{Tucker}(\mathcal{X})=max\{R_n，n=1,2,\ldots,N \} RankTucker​(X)=max{Rn​，n=1,2,…,N}
对于固定的$n$-秩，Tucker分解的唯一性不能保证，通常情况下会加上一些约束，如分解得到的因子单位正交约束等。对于Tucker分解优化的方法有HOSVD（High Order SVD）和HOOI（High-order orthogonal iteration）求解算法,如下所示
其中HOSVD算法通过对张量的每一个模式上做SVD分解，对各个模式上的因子矩阵进行求解，最后计算张量在各个模式上的投影之后的张量作为核张量。但是HOSVD 不能保证得到一个较好的近似。而HOSVD算法是HOOI算法的一个很好的初始化，HOOI算法是将张量分解看作是一个优化的过程，不断迭代得到分解结果。

BTD 分解

引入

当$\mathbb{R}$和$\mathbb{C}$的区别不重要时，我们用$\mathbb{K}$来表示$\mathbb{R}$或$\mathbb{C}$。现在定义一个张量$\mathcal{T}\in {\mathbb{K}}^{I\times J\times K}$,在CP分解中，张量$\mathcal{T}$分解为$R$个秩一张量的形式，公式为$$\mathcal{T}=\sum _{r=1}^R{a}_r⨂b_r⨂\cdots c_r=\sum _{r=1}^R(a_r{b}_r^T)⨂c_r$$可视化如下图所示
从公式可以看出，$a_r{b}_r^T$是一个秩-1的矩阵（低秩矩阵）。而低秩往往是一些局部的特征，但是在一些应用中，我们想要得到多尺度的特征，这个时候CP分解就无法直观地做到了。如果$a_r{b}_r^T$不是一个秩-1️的矩阵，而是一个秩为$K$的矩阵，那么显然$a_r$和$b_r$就不是向量了，而是一个秩为$L_r$的矩阵，因此我们就得到了BTD模型的一种形式即$$\mathcal{T}=\sum _{r=1}^K{A}_r{B}_r^T⨂c_r$$然后将这个模型进行推广得出BTD的一般形式。

BTD分解与求解

首先，我们给出一些定义：
1.张量$\mathcal{T}\in {\mathbb{K}}^{I_1\times I_2\times I_3}$的$mode$-$n$向量是由张量$\mathcal{T}$通过改变指标$i_n$并保持其他指标不变而得到的$I_n$维向量。$mode$-$n$向量泛化行向量和列向量。
2.张量$\mathcal{T}$的$mode$-$n$秩是它的$mode$-$n$向量张成子空间的维数，高阶张量的$mode$-$n$秩是矩阵列(行)秩的明显推广。
3.如果一个三阶张量的$mode$-1秩，$mode$-2秩，$mode$-3秩分别等于$L,M和N$，那么它就是$rank-(L,M,N）$。一个$rank$-$(1,1,1）$张量简称$rank$-1。
4.如果一个三阶张量$\mathcal{T}$等于三个向量的外积，那么它的秩为1，秩(与 $mode$-$n$秩相反)现在定义如下：
5.张量$\mathcal{T}$的秩是在线性组合中生成$\mathcal{T}$的$rank$-1张量的最小数目。

2008年Lieven De Lathauwer 等人提出了一种Block term decomposition（BTD）张量分解方法，它将一个$N$阶张量分解为$R$个成员张量的形式。定义一个张量$\mathcal{T}\in {\mathbb{K}}^{I\times J\times K}$，$\mathcal{T}$在秩($L,M,N$)下分解的可视化如下图所示，每一个成员张量是一个$Rank-(L,M,N）$的Tucker分解

其数学表达式为$$\mathcal{T}=\sum _{r=1}^R{\mathcal{D}}_r{\times }_1{A}_r{\times }_2{B}_r{\times }_3{C}_r$$其中，${\mathcal{D}}_r\in {\mathbb{K}}^{K\times M\times N}$是满秩$Rank-(L_r,M_r,N_r)$张量，而$A_r\in {\mathbb{K}}^{I\times L}（I\ge L）$，$B_r\in {\mathbb{K}}^{J\times M}（J\ge M）$，$C_r\in {\mathbb{K}}^{K\times N}（K\ge N）$都是全列秩，且$1\le r\le R$。

从公式和图解我们可以看出，BTD分解可以看作是我们之前介绍的CP分解和Tucker分解的组合形式。当$R=1$的时候，成员张量只有一个，此时这个分解就是$Rank-(L,M,N)$的Tucker分解。而当每一个成员张量是一个秩一分解的时候，该分解退化为一个CP分解 (CP分解是将一个张量$\mathcal{T}$分解为$R$个秩为1的张量形式),从而表现出BTD分解具有很强的泛化能力。

BTD分解现在的求解算法有很多，比如说非线性最小二乘法,交替最小二乘算法（ALS），张量对角算法等。$\mathcal{T}\in {\mathbb{K}}^{I\times J\times K}$采用ALS算法的求解步骤如下图所示
其中，算法中的QR分解是为了防止算法出现数值计算错误。

特殊的BTD分解

在此之前，我已经介绍了BTD分解的两种特殊形式：Tucker和CP分解，下面介绍两种特殊的张量分解。

$rank$-$(L_r,L_r,1)$分解

定义一个张量$\mathcal{T}\in {\mathbb{K}}^{I\times J\times K}$，$\mathcal{T}$在$rank$-$(L_r,L_r,1)$块和下分解为$$\mathcal{T}=\sum _{r=1}^R{\mathbf{E}}_{\mathbf{r}}\circ {\mathbf{c}}_{\mathbf{r}}$$其中$I\times J$矩阵${\mathbf{E}}_{\mathbf{r}}$是$rank$-$L_r$，且$1\le r\le R$。如果我们将${\mathbf{E}}_{\mathbf{r}}$因式分解为${\mathbf{A}}_{\mathbf{r}}\cdot {\mathbf{B}}_{\mathbf{r}}^{\mathbf{T}}$，其中矩阵${\mathbf{A}}_{\mathbf{r}}\in {\mathbb{K}}^{\mathbf{I}\times {\mathbf{L}}_{\mathbf{r}}}$和矩阵${\mathbf{B}}_{\mathbf{r}}\in {\mathbb{K}}^{\mathbf{J}\times {\mathbf{L}}_{\mathbf{r}}}$是$rank$-$L_r$，因此，可以将上式写成$$\mathcal{T}=\sum _{r=1}^R({\mathbf{A}}_{\mathbf{r}}\cdot {\mathbf{B}}_{\mathbf{r}}^{\mathbf{T}})\circ {\mathbf{c}}_{\mathbf{r}}$$同时，如果对角矩阵${\mathbf{D}}_{\mathbf{r}}$能被解释为一个$(L_r,L_r,1)$张量，那么还可以写成$$\mathcal{T}=\sum _{r=1}^R{\mathbf{D}}_{\mathbf{r}}{\cdot }_1{\mathbf{A}}_{\mathbf{r}}{\cdot }_2{\mathbf{B}}_{\mathbf{r}}{\cdot }_3{\mathbf{c}}_{\mathbf{r}}$$其中$c_r$是单位范数，张量分解如下图所示求解算法与BTD算法类似，相应的是没有最后求解core tensor部分。

$rank$-$(L,L,1)$分解

当$rank$-$(L_r,L_r,1)$的每一个$L_r$都取值为$L$时i，即可得到$rank$-$(L,L,1)$分解即：$$\mathcal{T}=\sum _{r=1}^R{\mathbf{E}}_{\mathbf{r}}\circ {\mathbf{c}}_{\mathbf{r}}$$其中$I\times J$矩阵${\mathbf{E}}_{\mathbf{r}}$是$rank$-$L$，图解为：