浅谈张量分解（五）：稀疏张量的CP分解

在前面的文章中，我们已经讨论了稀疏张量的Tucker分解，并介绍了如何采用梯度下降训练出一个合理的分解结构，与Tucker分解略有不同，这篇文章将介绍在数学表达式上更为简洁的CP分解，同时讨论如何利用梯度下降训练出稀疏张量的CP分解结构。

需要注意的是，本文在撰写过程中尽量避开了繁琐的数学推导，阅读以下内容只需要具备以下几点即可：

• 首先，知道高阶张量实质上是多维数组，且仅仅是向量和矩阵的高阶泛化；
• 其次，对高等数学所学习的求偏导数较为熟悉；
• 最后，对矩阵计算较为熟悉且感兴趣。

1 张量的CP分解模型

一般而言，给定一个大小为n_1 \times n_2 \times n_3的张量\mathcal{X}，其CP分解可以写成如下形式，即

\mathcal{X} \approx \sum_{r=1} ^{R} {A(:,r) \otimes B(:,r) \otimes C(:,r)}

其中，矩阵A,B,C的大小分别为n_1 \times R, n_2 \times R, n_3 \times R，被称为因子矩阵（factor matrices），符号“\otimes ”表示张量积，对于向量而言，这里的张量积与外积是等价的。张量\mathcal{X}任意位置\left(i,j,k\right)上的元素估计值为

x_{ijk} \approx \sum _{r=1} ^{R} {a_{ir}b_{jr}c_{kr}}.

图1 CP分解示意图（图片来源：Scalable tensor factorizations for incomplete data）

若将该逼近问题转化为一个无约束优化问题，且目标函数为残差的平方和，即

f \left(A, B, C\right)=\frac{1}{2} \sum_{i=1}^{n_1} {\sum_{j=1}^{n_2} {\sum_{k=1}^{n_3} {\left(x_{ijk}-\sum_{r=1}^{R}{a_{ir}b_{jr}c_{kr}}\right)^2}}}.

知道了要优化的目标函数后，我们就能够通过梯度下降法训练出一个合理的分解结构。另外，对梯度下降法有所了解的读者都知道，梯度下降法的原理很简单，我们只需要掌握高等数学中的求导就可以了。

接下来，我们来看看如何对目标函数f \left(A,B,C\right)求偏导，并给出决策变量的梯度下降更新公式。

2 “元素级”的梯度下降更新公式

将张量\mathcal{X}被观测到的元素相应的索引集合记为\left(i,j,k\right) \in \Omega，前面的目标函数可以改写为

f \left(A,B,C\right)=\frac{1}{2} \sum_{\left(i,j,k\right) \in \Omega}{\left(x_{ijk}-\sum_{r=1}^{R}{a_{ir}b_{jr}c_{kr}}\right)^2}

对目标函数中的决策变量a_{ir},b_{jr},c_{kr}求偏导数，则我们能够很容易地得到

\frac{\partial f \left(A,B,C\right)}{\partial a_{ir}}=\sum_{j,k:(i,j,k) \in \Omega}{\left(x_{ijk}-\sum_{r=1}^{R}{a_{ir}b_{jr}c_{kr}}\right) \left(-b_{jr}c_{kr}\right)},

\frac{\partial f \left(A,B,C\right)}{\partial b_{jr}}=\sum_{i,k:(i,j,k) \in \Omega}{\left(x_{ijk}-\sum_{r=1}^{R}{a_{ir}b_{jr}c_{kr}}\right) \left(-a_{ir}c_{kr}\right)},

\frac{\partial f \left(A,B,C\right)}{\partial c_{kr}}=\sum_{i,j:(i,j,k) \in \Omega}{\left(x_{ijk}-\sum_{r=1}^{R}{a_{ir}b_{jr}c_{kr}}\right) \left(-a_{ir}b_{jr}\right)}.

其中，需要特殊指出的是，求和符号下面的j,k:(i,j,k) \in \Omega, i,k:(i,j,k) \in \Omega, i,j:(i,j,k) \in \Omega分别表示矩阵\mathcal{X}(i,:,:),\mathcal{X}(:,j,:),\mathcal{X}(:,:,k)被观测到的元素其索引构成的集合。进一步，决策变量a_{ir},b_{jr},c_{kr}的更新公式可以分别写成：

a_{ir} \Leftarrow a_{ir}- \alpha \sum_{j,k:(i,j,k) \in \Omega}{\left(x_{ijk}-\sum_{r=1}^{R}{a_{ir}b_{jr}c_{kr}}\right) \left(-b_{jr}c_{kr}\right)},

b_{jr} \Leftarrow b_{jr}- \alpha \sum_{i,k:(i,j,k) \in \Omega}{\left(x_{ijk}-\sum_{r=1}^{R}{a_{ir}b_{jr}c_{kr}}\right) \left(-a_{ir}c_{kr}\right)},

c_{kr} \Leftarrow c_{kr} - \alpha \sum_{i,j:(i,j,k) \in \Omega}{\left(x_{ijk}-\sum_{r=1}^{R}{a_{ir}b_{jr}c_{kr}}\right) \left(-a_{ir}b_{jr}\right)}.

在这些更新公式中，\alpha为梯度下降的学习率。到这里，我们会发现，上述过程与前面的浅谈张量分解（一）：如何简单地进行张量分解？一文中对稀疏张量进行Tucker分解的过程完全相同，只不过是目标函数和决策变量略有不同而已。那么问题来了，这篇文章的意义是什么呢？

接下来，我们来看看在运算程序时更为有效的一种梯度下降更新策略，即“矩阵化”的梯度下降更新公式。

3 "矩阵化"的梯度下降更新公式

上面给出的更新公式是以“元素级”的形式，熟悉MATLAB的读者可能也知道，这种形式的更新公式在运算相应的程序时效率是很低的，且不能发挥出MATLAB自身在矩阵计算方面的优势。不过，值得庆幸的是，这些更新公式都具有“矩阵化”的数学表达式。

定义一个与张量\mathcal{X}大小相同的张量\mathcal{W}，当(i,j,k) \in \Omega时，w_{ijk}=1，否则w_{ijk}=0；再令\mathcal{E}=\mathcal{X}-\sum_{r=1} ^{R} {A(:,r) \otimes B(:,r) \otimes C(:,r)}，则决策变量A,B,C的更新公式分别为

A \Leftarrow A + \alpha \left(\mathcal{W}\ast \mathcal{E}\right)_{(1)}\left(C \odot B\right),

B \Leftarrow B + \alpha \left(\mathcal{W}\ast \mathcal{E}\right)_{(2)}\left(C \odot A\right),

C \Leftarrow C + \alpha \left(\mathcal{W}\ast \mathcal{E}\right)_{(3)}\left(B \odot A\right).

其中，运算符“\ast ”表示点乘（element-wise multiplication），若大小均为n_1 \times n_2 \times n_3的张量\mathcal{B}和\mathcal{C}进行点乘，即\mathcal{A}=\mathcal{B} \ast \mathcal{C}，则a_{ijk}=b_{ijk} \cdot c_{ijk}；运算符“\odot ”表示Khatri-Rao积，其定义可参考前面的浅谈张量分解（二）：张量分解的数学基础一文。

---------------------一些简单的推导---------------------

以因子矩阵 A 的更新公式为例，令梯度H=-\left(\mathcal{W} \ast \mathcal{S} \right)_{(1)} \left(C \odot B\right)，由于

\left(\mathcal{W} \ast \mathcal{S} \right)_{(1)}=\left[ \begin{array}{ccccccc}g_{111} & \cdots & g_{1n_21} & \cdots & g_{11n_3} & \cdots & g_{1n_2n_3} \\\vdots & \ddots & \vdots & & \vdots & \ddots & \vdots \\g_{n_111} & \cdots & g_{n_1n_21} & \cdots & g_{n_11n_3} & \cdots & g_{n_1n_2n_3} \\\end{array} \right] , g_{ijk}=w_{ijk}e_{ijk},i=1,2,...,n_1,j=1,2,...,n_2,k=1,2,...,n_3.

故矩阵 \left(\mathcal{W} \ast \mathcal{S} \right)_{(1)} 的第 i 行：

\left[ \begin{array}{ccccccc}w_{i11}e_{i11} & \cdots & w_{in_21}e_{in_21} & \cdots & w_{i1n_3}e_{i1n_3} & \cdots & w_{in_2n_3}e_{in_2n_3} \\\end{array} \right].

由于\left(C \odot B\right)=\left[ \begin{array}{ccc} c_{11}b_{11} & \cdots & c_{1R}b_{1R} \\ \vdots & \ddots & \vdots \\c_{11}b_{n_21} & \cdots & c_{1R}b_{n_2R} \\ \vdots & & \vdots \\ c_{n_31}b_{11} & \cdots & c_{n_3R}b_{1R} \\ \vdots & \ddots & \vdots \\c_{n_31}b_{n_21} & \cdots & c_{n_3R}b_{n_2R} \\\end{array} \right] ，故该矩阵的第 r 列为：

\left[ \begin{array}{ccccccc}c_{1r}b_{1r} & \cdots & c_{1r}b_{n_2r} & \cdots & c_{n_3r}b_{1r} & \cdots & c_{n_3r}b_{n_2r} \\\end{array} \right]^T .

综上，我们可以计算出矩阵 H 的第 i 行第 r 列元素为：

H(i,r)=-\sum_{(i,j,k) \in \Omega}{w_{ijk}e_{ijk} c_{kr}b_{jr}}=-\sum_{j,k:(i,j,k) \in \Omega}{e_{ijk} b_{jr}c_{kr}} .

从而，我们可以发现，上述“矩阵级”的更新公式表达式与“元素级”等价。

----------------------------------------------------------------

有了这些“矩阵化”的更新公式后，我们还需要考虑一个问题，即对于阶数大于3的稀疏张量如何进行CP分解呢？是否存在类似的梯度下降更新公式呢？

4 延伸：更高阶稀疏张量的CP分解

给定一个大小为n_1 \times n_2 \times ... \times n_d的张量\mathcal{X}，并将其CP分解写成如下形式：

\mathcal{X} \approx \sum_{r=1} ^{R} {A^{(1)}(:,r) \otimes A^{(2)}(:,r) \otimes ... \otimes A^{(d)}(:,r)}

其中，A^{(1)},A^{(2)},...,A^{(d)}分别是大小为n_1 \times R,n_2 \times R,..., n_d \times R的因子矩阵，采用梯度下降，令\mathcal{E}=\mathcal{X}-\sum_{r=1} ^{R} {A^{(1)}(:,r) \otimes A^{(2)}(:,r) \otimes ... \otimes A^{(d)}(:,r)}，则决策变量A^{(1)},A^{(2)},...,A^{(d)}的更新公式可以分别写成

A^{(1)} \Leftarrow A^{(1)}+\alpha \left(\mathcal{W} \ast \mathcal{E} \right)_{(1)} \left(A^{(d)} \odot A^{(d-1)} \odot ... \odot A^{(3)}\odot A^{(2)}\right),

A^{(2)} \Leftarrow A^{(2)}+\alpha \left(\mathcal{W} \ast \mathcal{E} \right)_{(2)} \left(A^{(d)} \odot A^{(d-1)} \odot ... \odot A^{(3)} \odot A^{(1)} \right),

... ...

A^{(d)} \Leftarrow A^{(d)}+\alpha \left(\mathcal{W} \ast \mathcal{E} \right)_{(d)} \left(A^{(d-1)} \odot A^{(d-2)} \odot ... \odot A^{(2)} \odot A^{(1)} \right).

到这里，我们已经介绍完了稀疏张量的CP分解，有兴趣的读者此时就可以设计实验了，不过，需要额外说明的是，CP分解中的R一般表示张量的秩（tensor rank），由于对于张量的秩的求解往往是一个NP-hard问题，所以读者在设计实验时要预先给定R的值。

5 参考文献

本文主要参考了Scalable tensor factorizations for incomplete data一文，文中介绍了如何利用CP分解来解决张量填补问题（tensor completion problem）。