张量CP分解原理及Python实现

张量分解

参考文献：Kolda TG, Bader BW. Tensor Decomposition and Application. SIAM Rev 2009;51:455–500. https://doi.org/10/dzcrv6.

张量可以视为多维数组，其“形状”取决于张量的阶（order）数。标量是第0阶张量，向量是第1阶张量，矩阵是第2阶张量，第3阶或阶数更高的张量被称为高阶张量（higher-order tensor），一般提到的张量都是特指高阶张量。

在矩阵中，我们需要用第$i$行第$j$列的形式定位一个元素，需要两个索引确定一个元素的位置。如果要确定高阶张量中某一元素的位置则需要更多的索引，所需索引数量与张量阶数相同。如下图所示的三阶张量，就可以用$(i,j,k)$的形式确定元素位置。

数学符号说明

1. 向量以加粗小写字母表示，如$\mathbf{a}$、$\mathbf{b}$
2. 矩阵以加粗大写字母表示，如$\mathbf{A}$、$\mathbf{B}$
3. 高阶张量以加粗花体字母显示， 如$\mathcal{X}$、$\mathcal{Y}$
4. 一个由$M$个矩阵（向量）组成的集合表示为${\left\{{\mathbf{A}}_{[m]}\in {\mathbb{R}}^{I_m\times N}\right\}}_m^M$（${\left\{{\mathbf{a}}_{[m]}\in {\mathbb{R}}^{I_m\times N}\right\}}_m^M$）
5. 定位张量中某一具体元素时以下标展示，如向量$\mathbf{a}$中第$i$个元素为$a_i$，矩阵$\mathbf{A}$的第$i$行第$j$列元素为$a_{ij}$,张量$\mathcal{X}$的某一元素可用$x_{ijk}$表示
6. 元素切片表示方法与numpy保持一致，用下标表示。矩阵$\mathbf{A}$切片如${\mathbf{a}}_{i:}$表示矩阵第$i$行，张量切片如下图所示。

6. ${\left|\left|\mathcal{X}\right|\right|}_p$表示张量的$p$范数，缺省状态下取$p=2$
   $${\left|\left|\mathcal{X}\right|\right|}_p={\left(\sum _i\sum _j\sum _k{\left|x_{ijk}\right|}^p\right)}^{\frac{1}{p}}$$

7. $⟨\mathcal{X}，\mathcal{Y}⟩$表示张量的内积
   $$⟨\mathcal{X}，\mathcal{Y}⟩=\sum _i\sum _j\sum _k{x}_{ijk}{y}_{ijk}$$

8. $\mathcal{X}\ast \mathcal{Y}$表示张量的哈达玛积，或element-wise product
   $$\mathcal{Z}=\mathcal{X}\ast \mathcal{Y}$$

   $$z_{ijk}=x_{ijk}{y}_{ijk}$$

9. $\mathbf{A}\otimes \mathbf{B}$表示Kronecker积，若$\mathbf{A}\in {\mathbb{R}}^{I\times J}，\mathbf{Y}\in {\mathbb{R}}^{M\times N}$，则 { A ⊗ B } ∈ R I M × J N \left\{\bold{A}\otimes\bold{B}\right\}\in\mathbb{R}^{IM\times JN} {A⊗B}∈RIM×JN，计算方式如下：
   $$\mathbf{A}\otimes \mathbf{B}=\left[\begin{array}{cccc}{a}_{11}\mathbf{B}& a_{12}\mathbf{B}& \cdots & a_{1m_2}\mathbf{B}\\ a_{21}\mathbf{B}& a_{22}\mathbf{B}& \cdots & a_{2m_2}\mathbf{B}\\ ⋮& ⋮& \ddots & ⋮\\ a_{m_{1}1}\mathbf{B}+& +a_{m_{1}2}\mathbf{B}& \cdots & a_{m_1{m}_2}\mathbf{B}\end{array}\right]$$

10. $\mathbf{A}\odot \mathbf{B}$表示Khatri-Rao积，若$\mathbf{A}\in {\mathbb{R}}^{I\times J}，\mathbf{Y}\in {\mathbb{R}}^{K\times J}$，则 { A ⊙ B } ∈ R I K × J \left\{\bold{A}\odot\bold{B}\right\}\in\mathbb{R}^{IK\times J} {A⊙B}∈RIK×J，计算方式如下：
    $$\mathbf{A}\odot \mathbf{B}=\left({\mathbf{a}}_{:1}\otimes {\mathbf{b}}_{:1},{\mathbf{a}}_{:2}\otimes {\mathbf{b}}_{:2},\dots ,{\mathbf{a}}_{:\mathbf{n}}\otimes {\mathbf{b}}_{:\mathbf{n}}\right)$$

    $${\mathbf{A}}_{[1]}\odot {\mathbf{A}}_{[2]}\odot \cdots \odot {\mathbf{A}}_{[\mathbf{M}]}\dot{=}\underset{m=1}{\overset{M}{⨀}}{\mathbf{A}}_{[\mathbf{m}]}$$

11. $\mathbf{a}\circ \mathbf{b}$表示向量的外积
    $$\mathbf{a}\circ \mathbf{b}=\mathbf{a}{\mathbf{b}}^T$$

12. $\mathcal{X}{\bar{\times }}_n\mathbf{Y}$，表示张量$\mathcal{X}$与矩阵$\mathbf{Y}$的$n$模态积，若$\mathcal{X}\in {\mathbb{R}}^{I_1\times I_2\times \cdots \times I_N}，\mathbf{Y}\in {\mathbb{R}}^{m\times I_n}$，则 { X × ˉ n Y } ∈ R I 1 × ⋯ × I n − 1 × m × I n + 1 × ⋯ × I N \left\{\boldsymbol{\mathcal{X}}\bar{\times}_n\boldsymbol{Y}\right\}\in\mathbb{R}^{I_1\times\dots\times I_{n-1}\times m \times I_{n+1}\times\dots\times I_N} {X×ˉ​n​Y}∈RI1​×⋯×In−1​×m×In+1​×⋯×IN​
    $$\mathcal{X}{\bar{\times }}_1{\mathbf{Y}}_{[1]}{\bar{\times }}_2{\mathbf{Y}}_{[2]}{\bar{\times }}_3\cdots {\bar{\times }}_M{\mathbf{Y}}_{[\mathbf{M}]}=\mathcal{X}\prod _{m=1}^M{\bar{\times }}_m{Y}_{[m]}$$

13. ${\mathbf{X}}^{†}$表示矩阵$\mathbf{X}$的广义逆（Moore–Penrose pseudoinverse）

14. ${\mathcal{X}}_{(n)}$表示张量以第$n$模态（$n-mode$）展开得到的矩阵，若$\mathcal{X}\in {\mathbb{R}}^{4\times 3\times 2}$，则
    $$\begin{gathered} {\mathcal{X}}_{\left(1\right)}=[\mathcal{X}\left(:,:,1\right),\mathcal{X}\left(:,:,2\right)] \\ {\mathcal{X}}_{\left(2\right)}=[\mathcal{X}{\left(:,:,1\right)}^T,\mathcal{X}{\left(:,:,2\right)}^T] \\ {\mathcal{X}}_{\left(3\right)}=[\mathcal{X}{\left(:,1,:\right)}^T,\mathcal{X}{\left(:,2,:\right)}^T,\mathcal{X}{\left(:,3,:\right)}^T] \end{gathered}$$

CP分解

基本概念

首先我们将可以由多个向量外积得到的张量成为秩1（rank-one）张量，如$\mathcal{Y}=\mathbf{a}\circ \mathbf{b}\circ \mathbf{c}$就是一个秩1矩阵，CP分解就是将一个张量分解为多个形状相同的秩1张量的和，或者说用多个秩1张量去近似表示原始张量，如下图所示：

设$\mathcal{X}\in {\mathbb{R}}^{I\times J\times K}$，$\mathbf{A}=\left({\mathbf{a}}_1,{\mathbf{a}}_2,\dots ,{\mathbf{a}}_R\right)\in {\mathbb{R}}^I$，$\mathbf{B}=\left({\mathbf{b}}_1,{\mathbf{b}}_2,\dots ,{\mathbf{b}}_R\right)\in {\mathbb{R}}^J$，$\mathbf{C}=\left({\mathbf{c}}_1,{\mathbf{c}}_2,\dots ,{\mathbf{c}}_R\right)\in {\mathbb{R}}^K$，则CP分解可以公式化为：
$$\mathcal{X}\approx \sum _{r=1}^R{\mathbf{a}}_r\circ {\mathbf{b}}_r\circ {\mathbf{c}}_r$$
CP分解公式还可以进一步矩阵化为：
$${\mathcal{X}}_{(1)}\approx \mathbf{A}{\left(\mathbf{C}\odot \mathbf{B}\right)}^T$$

$${\mathcal{X}}_{(2)}\approx \mathbf{B}{\left(\mathbf{C}\odot \mathbf{A}\right)}^T$$

$${\mathcal{X}}_{(3)}\approx \mathbf{C}{\left(\mathbf{B}\odot \mathbf{A}\right)}^T$$

与矩阵类似，张量也存在秩的概念，使$\mathcal{X}={\sum }_{r=1}^R{\mathbf{a}}_r\circ {\mathbf{b}}_r\circ {\mathbf{c}}_r$时最小的$R$值，或者说能够完美拟合原始张量所需的最少秩1张量数量，即为张量$\mathcal{X}$的秩。

在实际进行CP分解时，通常还会对$\mathbf{A}、\mathbf{B}、\mathbf{C}$按列归一化，并将其范数作为权重存储为向量$\mathbf{\lambda }\in {\mathbb{R}}^R$，则新的公式为：
$$\mathcal{X}\approx \sum _{r=1}^R{\mathbf{\lambda }}_r{\mathbf{a}}_r\circ {\mathbf{b}}_r\circ {\mathbf{c}}_r$$
对于更高阶的张量$\mathcal{X}\in {\mathbb{R}}^{I_1\times I_2\times \cdots \times I_N}$，通用表达式为：
$${\mathcal{X}}_{(n)}\approx {\mathbf{A}}_{[n]}\mathbf{\Lambda }{\left({\mathbf{A}}_{[N]}\odot \cdots \odot {\mathbf{A}}_{[n+1]}\odot {\mathbf{A}}_{[n-1]}\odot \cdots \odot {\mathbf{A}}_{[1]}\right)}^T$$
其中$\mathbf{\Lambda }=diag(\mathbf{\lambda })$，$diag$为对角矩阵化。

ALS-CP

ALS（alternating least squares、交替最小二乘法）

公式推导

下面以三阶张量$\mathcal{X}$为例介绍如何使用ALS计算CP分解的因子矩阵$\mathbf{A}、\mathbf{B}、\mathbf{C}$。首先我们的最终目标是让由

$\mathbf{A}、\mathbf{B}、\mathbf{C}$估计得到的张量$\hat{\mathcal{X}}$尽可能的接近原始张量$\mathcal{X}$，即：
$$min\left|\left|\hat{\mathcal{X}}-\mathcal{X}\right|\right|\text{ with }\hat{\mathcal{X}}=\sum _{r=1}^R{\mathbf{\lambda }}_r{\mathbf{a}}_r\circ {\mathbf{b}}_r\circ {\mathbf{c}}_r$$
为此，我们利用公式(12)(16)可以得到如下估计因子矩阵$\hat{\mathbf{A}}=\mathbf{A}\mathbf{\Lambda }$的方法：
$$min\left|\left|{\hat{\mathcal{X}}}_{(1)}-{\mathcal{X}}_{(1)}\right|\right|=min\left|\left|\hat{\mathbf{A}}{\left(\mathbf{C}\odot \mathbf{B}\right)}^T-{\mathcal{X}}_{(1)}\right|\right|$$

$$\hat{\mathbf{A}}={\mathcal{X}}_{(1)}{\left[{\left(\mathbf{C}\odot \mathbf{B}\right)}^T\right]}^{†}$$

根据Khatri-Rao的性质${\left(\mathbf{A}\odot \mathbf{B}\right)}^T={\mathbf{A}}^T\mathbf{A}\ast {\mathbf{B}}^T\mathbf{B}$、${\left(\mathbf{A}\odot \mathbf{B}\right)}^{†}={\left({\mathbf{A}}^T\mathbf{A}\ast {\mathbf{B}}^T\mathbf{B}\right)}^{†}{\left(\mathbf{A}\odot \mathbf{B}\right)}^T$，公式(19)可转化为：
$$\hat{\mathbf{A}}={\mathcal{X}}_{(1)}\left(\mathbf{C}\odot \mathbf{B}\right){\left({\mathbf{C}}^T\mathbf{C}\ast {\mathbf{B}}^T\mathbf{B}\right)}^{†}$$
在得到$\hat{\mathbf{A}}$后，我们又可以通过$\hat{\mathbf{A}}$和$\mathbf{C}$得到$\hat{\mathbf{B}}$，再通过$\hat{\mathbf{A}}$和$\hat{\mathbf{B}}$计算$\hat{\mathbf{C}}$，从而完成一轮$\mathbf{A}、\mathbf{B}、\mathbf{C}、\mathbf{\Lambda }$的更新，重复此过程即可使$\hat{\mathcal{X}}$不断接近原始张量$\mathcal{X}$，最终完成CP分解。
$$\hat{\mathbf{B}}={\mathcal{X}}_{(2)}\left(\mathbf{C}\odot \mathbf{A}\right){\left({\mathbf{C}}^T\mathbf{C}\ast {\mathbf{A}}^T\mathbf{A}\right)}^{†}$$

$$\hat{\mathbf{C}}={\mathcal{X}}_{(2)}\left(\mathbf{B}\odot \mathbf{A}\right){\left({\mathbf{B}}^T\mathbf{B}\ast {\mathbf{A}}^T\mathbf{A}\right)}^{†}$$

对高阶张量$\mathcal{X}\in {\mathbb{R}}^{I_1\times I_2\times \cdots \times I_N}$以及参数矩阵集合${\left\{{\mathbf{A}}_{[n]}\in {\mathbb{R}}^{I_n\times R}\right\}}_n^N$的通用表达式为:
$$\mathbf{V}=\left({\mathbf{A}}_{[N]}^T{\mathbf{A}}_{[N]}\ast \cdots {\mathbf{A}}_{[n+1]}^T{\mathbf{A}}_{[n+1]}\ast {\mathbf{A}}_{[n-1]}^T{\mathbf{A}}_{[n-1]}\ast \cdots \ast {\mathbf{A}}_{[1]}^T{\mathbf{A}}_{[1]}\right)$$

$${\hat{\mathbf{A}}}_{[n]}={\mathcal{X}}_{(n)}\left({\mathbf{A}}_{[N]}\odot \cdots \odot {\mathbf{A}}_{[n+1]}\odot {\mathbf{A}}_{[n-1]}\odot \cdots \odot {\mathbf{A}}_{[1]}\right){\mathbf{V}}^{†}$$

程序实现

程序参照《Tensor Decompositions and Applications》实现，流程如下：

• Step1: 初始化参数矩阵集合${\left\{{\mathbf{A}}_{[n]}\in {\mathbb{R}}^{I_n\times R}\right\}}_n^N=\text{A}$以及权重向量$\mathbf{\lambda }=\text{lbd}$:

# 使用tensorly提供的函数初始化
lbd, A = initialize_cp(tensor, rank, init='svd', svd='numpy_svd',
                       random_state=0,
                       normalize_factors=True)
# 或自己进行随机初始化
A = []
lbd = tl.ones(rank)
for n in range(N):
    A.append(tl.tensor(np.random.random((tensor.shape[n], rank))))

• Step2: 通过公式(23)计算$\mathbf{V}$：

# 使用None
V = None
for i in range(N):
	if i != n:
		if V is None:
			V = np.matmul(A[i].T, A[i])
		else:
			V = np.matmul(A[i].T, A[i]) * V

# 或将V初始化为RxR的1矩阵
V = np.ones((R, R))
for i in range(N):
	if i != n:
		V = np.matmul(A[i].T, A[i]) * V

• Step3: 通过公式(24)计算${\hat{\mathbf{A}}}_{[n]}=\text{A[n]}$：

T = khatri_rao(A, skip_matrix=n)
A[n] = np.matmul(np.matmul(tl.unfold(tensor, mode=n), T), np.linalg.pinv(V))

$$T=\left({\mathbf{A}}_{[N]}\odot \cdots \odot {\mathbf{A}}_{[n+1]}\odot {\mathbf{A}}_{[n-1]}\odot \cdots \odot {\mathbf{A}}_{[1]}\right)$$

$${\mathcal{X}}_{(n)}=\text{tl.unfold(tensor, mode=n)}$$

$${\mathbf{V}}^{†}=\text{np.linalg.pinv(V)}$$

• Step4: 对${\hat{\mathbf{A}}}_{[n]}$的每一列做归一化得到权重向量$\mathbf{\lambda }$以及${\mathbf{A}}_{[n]}$：

for r in range(R):
    lbd[r] = tl.norm(A[n][:, r])
A[n] = A[n] / tl.reshape(lbd, (1, -1))

• Step5: 结束迭代的条件包括损失值足够小或不再变小，因子矩阵的变化很小，目标值接近于零，或者超过预定义的最大迭代次数。下面只实现了第一种作为示例：

tensor_pred = tl.fold(np.matmul(np.matmul(A[0], np.diag(lbd)), 
                                khatri_rao(A, skip_matrix=0).T),
                      mode=0,
                      shape=tensor.shape)
if tl.norm(tensor - tensor_pred) <= 1e-7:
    return A, lbd

$$loss=\left|\left|\hat{\mathcal{X}}-\mathcal{X}\right|\right|=\text{tl.norm(tensor - tensor\_pred)}$$

$$\hat{\mathcal{X}}=\text{tensor\_pred},\text{ 计算方式参考公式(16)}$$

$$\mathbf{\Lambda }=\text{np.diag(lbd)}$$

完整程序如下：

import numpy as np
import tensorly as tl
from tensorly.decomposition._cp import initialize_cp

from tensorly.tenalg import khatri_rao
from tqdm import tqdm

def cp_als(tensor: np.ndarray, R=1, max_iter=100):
    N = tl.ndim(tensor)
    # Step 1
    lbd, A = initialize_cp(tensor, R, init='svd', svd='numpy_svd',
                           random_state=0,
                           normalize_factors=True)
    # A = []
    # for n in range(N):
    #     np.random.seed(N)
    #     A.append(tl.tensor(np.random.random((tensor.shape[n], rank))))
    # lbd = tl.ones(rank)

    for epoch in tqdm(range(max_iter)):
        for n in range(N):
            # Step 2
            V = np.ones((R, R))
            for i in range(N):
                if i != n:
                    V = np.matmul(A[i].T, A[i]) * V
            # V = None
            # for i in range(N):
            #     if i != n:
            #         if V is None:
            #             V = np.matmul(A[i].T, A[i])
            #         else:
            #             V = np.matmul(A[i].T, A[i]) * V


            # Step 3
            T = khatri_rao(A, skip_matrix=n)
            A[n] = np.matmul(np.matmul(tl.unfold(tensor, mode=n), T), np.linalg.pinv(V))

            # Step 4
            for r in range(R):
                lbd[r] = tl.norm(A[n][:, r])
            A[n] = A[n] / tl.reshape(lbd, (1, -1))
		# Step 5
        tensor_pred = tl.fold(np.matmul(np.matmul(A[0], np.diag(lbd)),
                                        khatri_rao(A, skip_matrix=0).T),
                              mode=0,
                              shape=tensor.shape)
        if tl.norm(tensor - tensor_pred) <= 1e-7:
            return A, lbd, epoch

    return A, lbd, max_iter

if __name__ == '__main__':
    np.random.seed(10086)
    inpt = tl.tensor(np.random.random((3, 3, 3)), dtype=np.float32)
    A, lbd, epoch = cp_als(inpt, R=5, max_iter=1000)
    tensor_pred = tl.fold(np.matmul(np.matmul(A[0], np.diag(lbd)),
                                    khatri_rao(A, skip_matrix=0).T),
                          mode=0,
                          shape=inpt.shape)

    print(tl.norm(inpt - tensor_pred), epoch)

附加说明

initialize_cp中使用SVD初始化的方式：

1. 为了初始化${\mathbf{A}}_{[n]}$，需要先对${\mathcal{X}}_{(n)}$做SVD得到左奇异矩阵$\mathbf{U}\in {\mathbb{R}}^{M\times N}$:

$${\mathcal{X}}_{(n)}=\mathbf{U}\mathbf{\Sigma }\mathbf{V}$$

2. 若选取的$\mathbf{R}\le N$，取左奇异矩阵$\mathbf{U}$的前$\mathbf{R}$列作为${\mathbf{A}}_{[n]}$即可：

$${\mathbf{A}}_{[1]}=\left[{\mathbf{u}}_{:1},{\mathbf{u}}_{:2},\dots ,{\mathbf{u}}_{:R}\right]$$

3. 若选取的$\mathbf{R}>N$，则还需要生成随机矩阵$\mathbf{F}\in {\mathbb{R}}^{M\times \left(R-N\right)}$拼接到矩阵$\mathbf{U}$的右侧作为${\mathbf{A}}_{[n]}$

BP-CP

BP（Back Propagation，反向传播），利用梯度下降法求解CP参数${\left\{{\mathbf{A}}_{[n]}\in {\mathbb{R}}^{I_n\times R}\right\}}_{n=1}^N$

公式推导

与ALS-CP相同，先以三阶张量为例，我们的目标是让由$\mathbf{A}、\mathbf{B}、\mathbf{C}$估计得到的张量$\hat{\mathcal{X}}$尽可能的接近原始张量$\mathcal{X}$，因此将损失/目标函数设置为：
$$\begin{array}{rl}\mathbf{L}& =\frac{1}{2}{\left[{\mathcal{X}}_{(1)}-\mathbf{A}(\mathbf{C}\odot \mathbf{B}{)}^T\right]}^2\\ & =\frac{1}{2}{\left[{\mathcal{X}}_{(2)}-\mathbf{B}(\mathbf{C}\odot \mathbf{A}{)}^T\right]}^2\\ & =\frac{1}{2}{\left[{\mathcal{X}}_{(3)}-\mathbf{C}(\mathbf{B}\odot \mathbf{A}{)}^T\right]}^2\end{array}$$
乘以$\frac{1}{2}$是为了方便求导,为了简化表达式，设$\mathbf{\Theta }={\mathcal{X}}_{(1)}-\mathbf{A}(\mathbf{C}\odot \mathbf{B}{)}^T$，则公式(34)可以化作$loss=\frac{1}{2}{\left(\mathbf{\Theta }\right)}^2$，对$\mathbf{A}$求偏导的结果为：
$$\begin{array}{rl}\frac{\partial \mathbf{L}}{\partial \mathbf{A}}& =2\times \frac{1}{2}\mathbf{\Theta }\times \frac{\partial \mathbf{\Theta }}{\partial \mathbf{A}}\\ & =\mathbf{\Theta }\frac{\partial {\mathcal{X}}_{(1)}-\mathbf{A}(\mathbf{C}\odot \mathbf{B}{)}^T}{\partial \mathbf{A}}\\ & =-\mathbf{\Theta }\left(\mathbf{C}\odot \mathbf{B}\right)\end{array}$$
同理可以算出$\frac{\partial \mathbf{L}}{\partial \mathbf{B}}$与$\frac{\partial \mathbf{L}}{\partial \mathbf{C}}$，然后用梯度下降法更新参数即可，下式中$lr$表示学习率：
$$\begin{array}{r}\mathbf{A}\leftarrow \mathbf{A}-lr\ast \frac{\partial \mathbf{L}}{\partial \mathbf{A}}\\ \mathbf{B}\leftarrow \mathbf{B}-lr\ast \frac{\partial \mathbf{L}}{\partial \mathbf{B}}\\ \mathbf{C}\leftarrow \mathbf{C}-lr\ast \frac{\partial \mathbf{L}}{\partial \mathbf{C}}\end{array}$$
综上所述，高阶张量的通用公式如下，损失函数为：
$$\begin{array}{rl}\mathbf{L}& =\frac{1}{2}{\left[{\mathcal{X}}_{(n)}-{\mathbf{A}}_{[n]}{\left({\mathbf{A}}_{[N]}\odot \cdots \odot {\mathbf{A}}_{[n+1]}\odot {\mathbf{A}}_{[n-1]}\odot \cdots \odot {\mathbf{A}}_{[1]}\right)}^T\right]}^2\end{array}$$
偏导为：
$$\begin{array}{rl}\frac{\partial \mathbf{L}}{\partial {\mathbf{A}}_{[n]}}& =2\times \frac{1}{2}\mathbf{\Theta }\times \frac{\partial \mathbf{\Theta }}{\partial {\mathbf{A}}_{[n]}}\\ & =\mathbf{\Theta }\frac{\partial {\mathcal{X}}_{(n)}-{\mathbf{A}}_{[n]}{\left({\mathbf{A}}_{[N]}\odot \cdots \odot {\mathbf{A}}_{[n+1]}\odot {\mathbf{A}}_{[n-1]}\odot \cdots \odot {\mathbf{A}}_{[1]}\right)}^T}{\partial {\mathbf{A}}_{[n]}}\\ & =-\mathbf{\Theta }\left({\mathbf{A}}_{[N]}\odot \cdots \odot {\mathbf{A}}_{[n+1]}\odot {\mathbf{A}}_{[n-1]}\odot \cdots \odot {\mathbf{A}}_{[1]}\right)\end{array}$$

$$\mathbf{\Theta }={\mathcal{X}}_{(n)}-{\mathbf{A}}_{[n]}{\left({\mathbf{A}}_{[N]}\odot \cdots \odot {\mathbf{A}}_{[n+1]}\odot {\mathbf{A}}_{[n-1]}\odot \cdots \odot {\mathbf{A}}_{[1]}\right)}^T$$

参数更新公式为：
$${\mathbf{A}}_{[n]}\leftarrow {\mathbf{A}}_{[n]}-lr\ast \frac{\partial \mathbf{L}}{\partial {\mathbf{A}}_{[n]}}$$

程序实现

• Step1: 初始化参数矩阵集合${\left\{{\mathbf{A}}_{[n]}\in {\mathbb{R}}^{I_n\times R}\right\}}_n^N=\text{A}$、梯度矩阵（Jacobian矩阵）集合${\left\{\frac{\mathbf{L}}{{\mathbf{A}}_{[n]}}\in {\mathbb{R}}^{I_n\times R}\right\}}_n^N=\text{grad\_A}$以及权重向量$\mathbf{\lambda }=\text{lbd}$:

# 使用tensorly提供的函数初始化
lbd, A = initialize_cp(tensor, rank, init='random', svd='numpy_svd',
                       random_state=0,
                       normalize_factors=True)
# 或自己进行随机初始化
A = []
lbd = tl.ones(rank)
for n in range(N):
    A.append(tl.tensor(np.random.random((tensor.shape[n], rank))))

• Step2: 在每个epoch开始时计算$\Theta =\text{theta}$，并将梯度矩阵（Jacobian矩阵）集合${\left\{\frac{\mathbf{L}}{{\mathbf{A}}_{[n]}}\in {\mathbb{R}}^{I_n\times R}\right\}}_n^N=\text{grad\_A}$归零:

tensor_pred = tl.fold(np.matmul(np.matmul(A[0], np.diag(lbd)), 
                                khatri_rao(A, skip_matrix=0).T),
                      mode=0,
                      shape=tensor.shape)
theta = (tensor - tensor_pred)
grad_A = []

• Step3: 计算每个参数矩阵${\mathbf{A}}_{[n]}$对应的梯度矩阵$\frac{\mathbf{L}}{{\mathbf{A}}_{[n]}}$并存储，对应公式47：

for n in range(N):
    grad_A.append(np.zeros_like(A[n]))
    grad_A[n] = np.matmul(tl.unfold(theta, n), khatri_rao(A, skip_matrix=n))

• Step4: 更新参数矩阵，对应公式49，减法变成加法是因为由公式47求得的梯度矩阵带有负号：

for n in range(N):
    A[n] = A[n] + lr * grad_A[n]

• Step5: 最后计算以下每个epoch的损失值并输出，方便观察收敛情况，这里对$\mathbf{L}$内元素求和输出：

loss = np.sum(0.5 * np.square(tl.unfold(theta, 0)))
print("epoch {}: loss={}".format(epoch, loss))

完整程序如下：

import numpy as np
import tensorly as tl
from tensorly.decomposition._cp import initialize_cp

from tensorly.tenalg import khatri_rao
from tqdm import tqdm


def cp_bp(tensor: np.ndarray, R=1, lr=1e-2, max_iter=100):
    N = tl.ndim(tensor)
    # Step 1
    lbd, A = initialize_cp(tensor, R, init='random', svd='numpy_svd',
                         random_state=0,
                         normalize_factors=True)
    
    for epoch in range(max_iter):
        # Step 2
        tensor_pred = tl.fold(np.matmul(np.matmul(A[0], np.diag(lbd)), 
        								khatri_rao(A, skip_matrix=0).T),
                              mode=0,
                              shape=tensor.shape)
        theta = (tensor - tensor_pred)
        grad_A = []
        
        # Step 3
        for n in range(N):
            grad_A.append(np.zeros_like(A[n]))
            grad_A[n] = np.matmul(tl.unfold(theta, n), khatri_rao(A, skip_matrix=n))
            
        # Step 4
        for n in range(N):
            A[n] = A[n] + lr * grad_A[n]
            
        # Step 5
        loss = np.sum(0.5 * np.square(tl.unfold(theta, 0)))
        print("epoch {}: loss={}".format(epoch, loss))
        
    return A, lbd

if __name__ == '__main__':
    np.random.seed(10086)
    inpt = tl.tensor(np.random.random((3, 3, 3)), dtype=np.float32)
    A, lbd = cp_bp(inpt, R=5, lr=1e-2, max_iter=100)
    tensor_pred_cp = tl.fold(np.matmul(np.matmul(A[0], np.diag(lbd)),
                                    khatri_rao(A, skip_matrix=0).T),
                          mode=0,
                          shape=inpt.shape)

    print("tensor_pred_cp: ", tl.norm(inpt - tensor_pred_cp), epoch)