【张量分解(二)】CP分解

本文是对论文Tensor Decompositions and Applications进行了翻译、整理、筛选和适当的补充，如何希望深入理解可以阅读原文。

相关文章：

【张量分解(一)】符号与基础知识
【张量分解(二)】CP分解

一、CP分解

1.1 定义

CP分解就是将一个张量分解成多个单秩张量的和。例如，给定一个三阶张量$\mathcal{X}\in\mathbb{R}^{I\times J\times K}$X∈RI×J×K，则CP分解可以写为
$\mathcal{X}\approx\sum_{r=1}^{R}\textbf{a}_r\circ\textbf{b}_r\circ\textbf{c}_r$X≈r=1∑R​ar​∘br​∘cr​
其中，$\circ$∘是指向量外积，$R$R是正整数且$\textbf{a}_r\in\mathbb{R}^I$ar​∈RI,$\textbf{b}_r\in\mathbb{R}^J$br​∈RJ,$\textbf{c}_r\in\mathbb{R}^K$cr​∈RK。下图展示了三阶张量的CP分解

将上面的CP分解展开，也可以写作
$x_{ijk}\approx\sum_{r=1}^R a_{ir}b_{jr}c_{kr},\quad i=1,2,\dots,I,j=1,2,\dots,J, k=1,2,\dots,K$xijk​≈r=1∑R​air​bjr​ckr​,i=1,2,…,I,j=1,2,…,J,k=1,2,…,K
此外，对于三阶张量来说，可以从通道切片(frontal slice)的角度表示CP分解
$\textbf{X}_k\approx\textbf{A}\textbf{D}^{(k)}\textbf{B}^T，\textbf{D}^{(k)}\equiv diag(\textbf{c}_{k:}),k=1,\dots,K$Xk​≈AD(k)BT，D(k)≡diag(ck:​),k=1,…,K
其中，$\textbf{X}_k$Xk​表示张量$\mathcal{X}$X的第k个通道切片。对于行切片和列切片也可以写出类似的公式。

1.2 张量矩阵化后的CP分解

在文章【张量分解(一)】符号与基础知识中介绍过张量的矩阵化。这里主要介绍将张量转换为矩阵后的CP分解。
首先，定义因子矩阵(factor matrices)为CP分解中组成单秩张量的同一维度的向量合并成的矩阵(这个表述有点绕)。具体来说，就是把所有的$\textbf{a}$a向量合并成一个矩阵$\textbf{A}=[\textbf{a}_1\quad\textbf{a}_2\quad\dots\quad\textbf{a}_R]$A=[a1​a2​…aR​]。同理，还可以合成因子矩阵$\textbf{B}$B和$\textbf{C}$C。那么矩阵化后的张量CP分解形式如下：
$\textbf{X}_{(1)}\approx\textbf{A}(\textbf{B}\odot\textbf{C})^T$X(1)​≈A(B⊙C)T
$\textbf{X}_{(2)}\approx\textbf{B}(\textbf{C}\odot\textbf{A})^T$X(2)​≈B(C⊙A)T
$\textbf{X}_{(3)}\approx\textbf{C}(\textbf{B}\odot\textbf{A})^T$X(3)​≈C(B⊙A)T
其中，$\odot$⊙表示Khatri-Rao积，$\textbf{X}_{(i)}$X(i)​表示张量$\mathcal{X}$X的模i矩阵化后的矩阵。

1.3 符号表示

为了更加简洁的表达，CP分解可以简写如下$\mathcal{X}\approx\lgroup\textbf{A},\textbf{B},\textbf{C}\rgroup$X≈⟮A,B,C⟯
实在是打不出空心方括号(摊手)，只能用$\lgroup\rgroup$⟮⟯代替了。


通常，假设矩阵$\textbf{A}$A,$\textbf{B}$B和$\textbf{C}$C的列向量是标准化后的向量，并且将提取出来的权重合并入向量$\mathrm{\lambda}\in\mathbb{R}^R$λ∈RR，因此CP分解还可以写成$\mathcal{X}\approx\sum_{r=1}^{R}\lambda_{r}\textbf{a}_r\circ\textbf{b}_r\circ\textbf{c}_r=\lgroup\mathrm{\lambda};\textbf{A},\textbf{B},\textbf{C}\rgroup$X≈r=1∑R​λr​ar​∘br​∘cr​=⟮λ;A,B,C⟯

1.4 高维扩展

先前主要介绍的是三阶张量的CP分解，主要是因为其具有广泛的适用性。对于N阶张量$\mathcal{X}\in\mathbb{R}^{I_1\times I_2\times \dots \times I_N}$X∈RI1​×I2​×⋯×IN​，其CP分解为
$\mathcal{X}\approx\sum_{r=1}^{R}\lambda_{r}\textbf{a}_r^{(1)}\circ\textbf{a}_r^{(2)}\circ\dots\circ\textbf{a}_r^{(N)}=\lgroup\mathrm{\lambda};\textbf{A}^{(1)},\textbf{A}^{(2)},\dots,\textbf{A}^{(N)}\rgroup$X≈r=1∑R​λr​ar(1)​∘ar(2)​∘⋯∘ar(N)​=⟮λ;A(1),A(2),…,A(N)⟯
其中，$\mathrm{\lambda}\in\mathbb{R}^R$λ∈RR且$\textbf{A}^{(n)}\in\mathbb{R}^{I_n\times R},n=1,2,\dots,N$A(n)∈RIn​×R,n=1,2,…,N

类似的，N阶张量$\mathcal{X}$X进行模n矩阵化后的CP分解为
$\textbf{X}_{(n)}\approx\textbf{A}^{(n)}\mathrm{\Lambda}(\textbf{A}^{(N)}\odot\dots\odot\textbf{A}^{n+1}\odot\textbf{A}^{n-1}\odot\dots\odot\textbf{A}^{(1)})^T$X(n)​≈A(n)Λ(A(N)⊙⋯⊙An+1⊙An−1⊙⋯⊙A(1))T
其中，对角矩阵$\mathrm{\Lambda}=diag(\mathrm{\lambda})$Λ=diag(λ)。

二、张量的秩(Tensor Rank)

2.1 张量秩的定义

用于生成张量$\mathcal{X}$X所需要的单秩张量的最小数量即为张量$\mathcal{X}$X的秩，用$rank{\mathcal{(X)}}$rank(X)表示。换个角度，张量的秩就是CP分解时单秩张量数量的最小值。

2.2 张量秩与矩阵秩

此外，张量的秩与矩阵秩的定义非常相似，但是二值的性质非常的不同。例如，实数张量的秩在实数域$\mathbb{R}$R和复数域$\mathbb{C}$C上可能会不同。另一个张量秩和矩阵秩的显著不同是，当前没有一个直接的方法来确定给定张量的秩。例如，Krushkal对特定的$9\times9\times9$9×9×9的张量进行分析，只能确定其秩在18到23之间。在实际应用中，张量的秩是通过CP分解来确定的。

2.3 张量的最大秩和典型秩

最大秩：一类张量能够达到的最大的秩称为张量的最大秩(maximum rank)。典型秩：一个从均匀连续分别中随机抽取元素所组成的张量中，出现概率大于0的任何秩。
具体来说，对于所有形状为$I\times J$I×J的矩阵，最大秩和典型秩均等于$min\{I,J\}$min{I,J}。但是对于张量来说，最大秩和典型秩可能不相同，而且典型秩可能不只一个。例如$2\times 2\times 2$2×2×2张量的典型秩为2或3，通过蒙特卡洛实验也可以发现秩为2的张量占79%，秩为3的张量占21%，秩为1的张量在理论上虽然可能，但是实际概率为0。
对于一般的三阶张量$\mathcal{X}\in\mathbb{R}^{I\times J\times K}$X∈RI×J×K，当前只知道其最大秩的一个弱上界
$rank(\mathcal{X})\leq min\{IJ,IK,JK\}$rank(X)≤min{IJ,IK,JK}
对于特定形状或类型的张量来说，有可能存在一些确定最大秩和典型秩的具体值或者范围的方法，可以参考原文Tensor Decompositions and Applications

三、唯一性

高阶张量的一个有趣的特性是它的秩分解是唯一的，而通常矩阵分解不是。

3.1 矩阵分解的不唯一性

对于秩为$R$R的矩阵$\textbf{X}\in\mathbb{R}^{I\times J}$X∈RI×J，其秩分解可以写为
$\textbf{X}=\textbf{AB}^T=\sum_{r=1}^R\textbf{a}_r\circ\textbf{b}_r$X=ABT=r=1∑R​ar​∘br​
具体来说，对于矩阵$\textbf{X}$X的SVD分解为$\mathrm{U\Sigma V}^T$UΣVT，为了与上面的秩分解对于，令$\textbf{A}=\mathrm{U\Sigma}$A=UΣ且$\textbf{B}=\mathrm{V}$B=V。但是，如果令$\textbf{A}=\mathrm{U\Sigma W}$A=UΣW且$\textbf{B}=\mathrm{VW}$B=VW,其中$\mathrm{W}$W是$R\times R$R×R的正交矩阵($W^TW=E$WTW=E)，同样也满足矩阵秩分解的定义。
换句话说，我们可以轻易的构造两个完全不同的单秩矩阵集合，但是集合中的矩阵相加就等于原始矩阵。而SVD分解的唯一性仅仅是因为正交约束的加入。

3.2 张量分解的唯一性

通常，在十分微弱的约束条件下，张量的CP分解就是唯一的。对于秩为$R$R的三阶张量$\mathcal{X}\in\mathbb{R}^{I\times J\times K}$X∈RI×J×K，其CP分解为
$\mathcal{X}=\sum_{r=1}^{R}\textbf{a}_r\circ\textbf{b}_r\circ\textbf{c}_r=\lgroup\textbf{A},\textbf{B},\textbf{C}\rgroup$X=r=1∑R​ar​∘br​∘cr​=⟮A,B,C⟯
而唯一性就是指上面的分解中是唯一可能的单秩矩阵的组合。当然，这是排除了缩放和重新排列后的唯一性。例如这里使用置换矩阵对分解后的单秩矩阵的列进行重排列
$\mathcal{X}=\lgroup\textbf{A},\textbf{B},\textbf{C}\rgroup=\lgroup\textbf{A}\Pi,\textbf{B}\Pi,\textbf{C}\Pi\rgroup$X=⟮A,B,C⟯=⟮AΠ,BΠ,CΠ⟯
其中，$\Pi$Π是$R\times R$R×R的置换矩阵。同样，对于将CP分解中的向量进行缩放也不影响CP分解的结果，例如
$\mathcal{X}=\sum_{r=1}^R(\alpha_r\textbf{a}_r)\circ(\beta_r\textbf{b}_r)\circ(\gamma_r\textbf{c}_r)$X=r=1∑R​(αr​ar​)∘(βr​br​)∘(γr​cr​)
其中，$\alpha_r\beta_r\gamma_r=1,r=1,...,R$αr​βr​γr​=1,r=1,...,R

3.3 CP分解唯一性的充分条件

对于CP分解$\mathcal{X}=\lgroup\textbf{A},\textbf{B},\textbf{C}\rgroup$X=⟮A,B,C⟯，令$k_A$kA​、$k_B$kB​、$k_C$kC​分别表示矩阵$\textbf{A}$A、$\textbf{B}$B、$\textbf{C}$C的秩，那么CP分解唯一的充分条件是$k_A+k_B+k_C\ge2R+2$kA​+kB​+kC​≥2R+2
将上面的条件扩展至N维，对于张量$\mathcal{X}=\sum_{r=1}^R\textbf{a}_{r}^{(1)}\circ \textbf{a}_{r}^{(2)}\circ\dots\circ\textbf{a}_{r}^{(N)}=\lgroup\textbf{A}^{(1)},\textbf{A}^{(2)},\dots,\textbf{A}^{(N)}\rgroup$X=r=1∑R​ar(1)​∘ar(2)​∘⋯∘ar(N)​=⟮A(1),A(2),…,A(N)⟯，其CP分解唯一性的充分条件为
$\sum_{n=1}^N k_{\textbf{A}^{(n)}}\ge2R+(N-1)$n=1∑N​kA(n)​≥2R+(N−1)

3.4 CP分解唯一性的必要条件

上面的充分条件在$R=2$R=2或$R=3$R=3的条件下，也是CP分解唯一性的必要条件，但是当$R>3$R>3则不成立。更加广泛的CP分解唯一性的必要条件为
$min\{rank(\textbf{A}\odot\textbf{B}),rank(\textbf{A}\odot\textbf{C}),rank(\textbf{B}\odot\textbf{C})\}=R$min{rank(A⊙B),rank(A⊙C),rank(B⊙C)}=R
推广的N维情况下，则
$min_{n=1,\dots,N}rank(\textbf{A}^{(1)}\odot\dots\odot\textbf{A}^{(n-1)}\odot\textbf{A}^{(n+1)}\odot\dots\odot\textbf{A}^{(N)})=R$minn=1,…,N​rank(A(1)⊙⋯⊙A(n−1)⊙A(n+1)⊙⋯⊙A(N))=R
但是，由于性质
$rank(\textbf{A}\odot\textbf{B})\le rank(\textbf{A}\otimes\textbf{B})\le rank(\textbf{A})\cdot rank(\textbf{B})$rank(A⊙B)≤rank(A⊗B)≤rank(A)⋅rank(B)
因此，N维下的必要条件可以扩展为
$min_{n=1,\dots,N}\Big(rank(\textbf{A}^{(1)})\cdot\dots\cdot rank(\textbf{A}^{(n-1)})\cdot rank(\textbf{A}^{(n+1)})\cdot\dots\cdot rank(\textbf{A}^{(N)})\Big)\ge R$minn=1,…,N​(rank(A(1))⋅⋯⋅rank(A(n−1))⋅rank(A(n+1))⋅⋯⋅rank(A(N)))≥R

3.5 CP分解唯一性的判断标准

对于秩为$R$R的三阶张量$\mathcal{X}\in\mathbb{R}^{I\times J\times K}$X∈RI×J×K，当满足条件
$R\le K并且R(R-1)\le I(I-1)J(J-1)/2$R≤K并且R(R−1)≤I(I−1)J(J−1)/2
则，其CP分解是唯一的。
类似的，对于秩为R的四阶张量$\mathcal{X}\in\mathbb{R}^{I\times J\times K\times L}$X∈RI×J×K×L，其CP分解唯一的条件是
$R\le L并且R(R-1)\le IJK(3IJK-IJ-IK-JK-I-J-K+3)/4$R≤L并且R(R−1)≤IJK(3IJK−IJ−IK−JK−I−J−K+3)/4

四、低秩近似与边界秩(border rank)

4.1 矩阵的低秩近似

给定一个秩为$R$R的矩阵$\textbf{A}$A，那么该矩阵的SVD分解可以写作：
$\textbf{A}=\sum_{r=1}^R\sigma_r\textbf{u}_r\circ\textbf{v}_r，其中\sigma_1\ge\sigma_2\ge\dots\ge\sigma_R$A=r=1∑R​σr​ur​∘vr​，其中σ1​≥σ2​≥⋯≥σR​
那么该矩阵的秩k近似，可以直接使用SVD分解中前k个部分，即
$\textbf{B}=\sum_{r=1}^k\sigma_r\textbf{u}_r\circ\textbf{v}_r$B=r=1∑k​σr​ur​∘vr​

4.2 张量的低秩近似

上面对于矩阵的结果并不适用于张量。给定一个秩为$R$R的三阶张量，其CP分解为
$\mathcal{X}=\sum_{r=1}^R\lambda_r\textbf{a}_r\circ\textbf{b}_r\circ\textbf{c}_r$X=r=1∑R​λr​ar​∘br​∘cr​
按上面矩阵的低秩近似来看，三阶张量的秩k近似也应该是其中k个部分的和，但实际情况并非如此。
Kolda提供过一个例子，对于一个三阶张量的单秩近似并不是秩2近似的组成部分(在矩阵的低秩分解中一定成立)。因此会得出一个推论，一个张量的最优秩k近似中的k个组成部分并不是按顺序求得的，而是需要同时被发现的。
总的来说，这个问题比较复杂，有时一个张量的最优秩k近似不一定存在。如果一个张量可以通过低秩的因式分解任意逼近，那么该张量就是一个退化张量。
举一个具体的例子来说，给定一个秩为3的具体三阶张量$\mathcal{X}\in\mathbb{R}^{I\times J\times K}$X∈RI×J×K为
$\mathcal{X}=\textbf{a}_1\circ\textbf{b}_1\circ\textbf{c}_2+\textbf{a}_1\circ\textbf{b}_2\circ\textbf{c}_1+\textbf{a}_2\circ\textbf{b}_1\circ\textbf{c}_1$X=a1​∘b1​∘c2​+a1​∘b2​∘c1​+a2​∘b1​∘c1​
其中,$\textbf{A}\in\mathbb{R}^{I\times 2},\textbf{B}\in\mathbb{R}^{J\times 2},\textbf{C}\in\mathbb{R}^{K\times 2}$A∈RI×2,B∈RJ×2,C∈RK×2是由于上式中对应的向量组成的，且这三个矩阵的列向量线性无关。
上面描述的张量可以使用下面的下面的秩2张量进行任意的近似
$\mathcal{Y}=n\Big(\textbf{a}_1+\frac{1}{n}\textbf{a}_2\Big)\circ\Big(\textbf{b}_1+\frac{1}{n}\textbf{b}_2\Big)\circ\Big(\textbf{c}_1+\frac{1}{n}\textbf{c}_2\Big)-n\textbf{a}_1\circ\textbf{b}_1\circ\textbf{c}_1$Y=n(a1​+n1​a2​)∘(b1​+n1​b2​)∘(c1​+n1​c2​)−na1​∘b1​∘c1​
原始的秩3张量$\mathcal{X}$X和近似的秩2张量$\mathcal{Y}$Y之间的误差为
$\Vert\mathcal{X}-\mathcal{Y}\Vert=\frac{1}{n}\Big\Vert\textbf{a}_2\circ\textbf{b}_2\circ\textbf{c}_1+\textbf{a}_2\circ\textbf{b}_1\circ\textbf{c}_2+\textbf{a}_1\circ\textbf{b}_2\circ\textbf{c}_2+\frac{1}{n}\textbf{a}_2\circ\textbf{b}_2\circ\textbf{c}_2\Big\Vert$∥X−Y∥=n1​∥∥∥​a2​∘b2​∘c1​+a2​∘b1​∘c2​+a1​∘b2​∘c2​+n1​a2​∘b2​∘c2​∥∥∥​
当然，这个误差可以任意的小。

4.3 边界秩(border rank)

在不存在最优低秩近似的情况下，可以考虑边界秩。其定义为，能够以任意非零误差充分近似给定张量的最小单秩张量的数量。形式化的定义为
$\widetilde{rank}(\mathcal{X})=min\{r|对于任意\epsilon>0,均存在一个张量\mathcal{E}满足\Vert\mathcal{E}\Vert<\epsilon且rank(\mathcal{X}+\mathcal{E})=r\}$rank(X)=min{r∣对于任意ϵ>0,均存在一个张量E满足∥E∥<ϵ且rank(X+E)=r}
显然，$\widetilde{rank}(\mathcal{X})\le rank(\mathcal{X})$rank(X)≤rank(X)

五、计算CP分解

本小节介绍怎么计算一个张量的CP分解。
在前面的小节中提到过，没有一个有限的算法可以确定张量的秩。而CP分解则是将待分解张量分解成$R$R个单秩张量，其中$R$R就是待分解张量的秩。因此，计算CP分解的第一个问题就是如何确定张量的秩。
大多数的CP求解思路是尝试不同的$R$R值来拟合待分解张量，直至找到一个最佳的分解。对于无噪声的数据，那么可以对$R$R的值从1,2,…这样逐步尝试，从而得到一个最优的CP分解。但是，前面介绍了张量的低秩近似，一个张量可以被一个更低秩的张量任意逼近，这在实际中有一些问题。

5.1 计算三阶张量的CP分解

假设CP分解中的$R$R取值已经确定，那么这里介绍一种求解CP分解的ALS(交替最小二乘法)算法。
令$\mathcal{X}\in\mathbb{R}^{I\times J\times K}$X∈RI×J×K，该算法的目标是计算一个包含$R$R个单秩张量的CP分解，使其尽量近似$\mathcal{X}$X，即
$min_{\hat{\mathcal{X}}}\Vert\mathcal{X}-\hat{\mathcal{X}}\Vert,其中\hat{\mathcal{X}}=\sum_{r=1}^{R}\lambda_r\textbf{a}_r\circ\textbf{b}_r\circ\textbf{c}_r=\lgroup\mathrm{\lambda};\textbf{A},\textbf{B},\textbf{C}\rgroup$minX^​∥X−X^∥,其中X^=r=1∑R​λr​ar​∘br​∘cr​=⟮λ;A,B,C⟯
交替最小二乘法(ALS)就是固定B和C，求解A；再固定A和C，求解B；再固定A和B，求解C。重复上面的过程，直至满足收敛条件。这就是ALS的思路。
固定两个张量来求解另外一个张量，这就变成了线性最小二乘的问题。例如，B和C固定，那么依照1.2节中张量矩阵化后的CP分解，那么就能把上面的最小化问题重写为
$min_{\hat{\textbf{A}}}\Vert\textbf{X}_{(1)}-\hat{\textbf{A}}(\textbf{C}\odot\textbf{B})^T\Vert_F,其中\hat{\textbf{A}}=\textbf{A}\cdot diag(\lambda)$minA^​∥X(1)​−A^(C⊙B)T∥F​,其中A^=A⋅diag(λ)
上面的最小化问题的最优解为
$\hat{\textbf{A}}=\textbf{X}_{(1)}[(\textbf{C}\odot\textbf{B})^T]^{-1}$A^=X(1)​[(C⊙B)T]−1
其中，-1是指张量的伪逆，而Khatri-Rao积的伪逆可以进行变换，因此上面的最优解还可以写作
$\hat{\textbf{A}}=\textbf{X}_{(1)}(\textbf{C}\odot\textbf{B})(\textbf{C}^T\textbf{C}*\textbf{B}^T\textbf{B})^{-1}$A^=X(1)​(C⊙B)(CTC∗BTB)−1
这个版本写法的最优解有一个优势，仅需要求解一个$R\times R$R×R矩阵的伪逆，而不是$JK\times R$JK×R矩阵的伪逆。最后，对矩阵$\hat{\textbf{A}}$A^的列进行标准化后就得到了矩阵$\textbf{A}$A。

5.2 ALS算法在高维张量上的应用

给定N阶张量$\mathcal{X}\in\mathbb{R}^{I_1\times I_2\times\dots\times I_N}$X∈RI1​×I2​×⋯×IN​，使用CP分解将其分解为$R$R个单秩矩阵的ALS算法。