目录

  • 张量分解
    • 数学符号说明
    • CP分解
      • 基本概念
      • ALS-CP
        • 公式推导
        • 程序实现
        • 附加说明
      • BP-CP
        • 公式推导
        • 程序实现

张量分解

参考文献：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 i i行第 j j j列的形式定位一个元素，需要两个索引确定一个元素的位置。如果要确定高阶张量中某一元素的位置则需要更多的索引，所需索引数量与张量阶数相同。如下图所示的三阶张量，就可以用 ( i , j , k ) (i,j,k) (i,j,k)的形式确定元素位置。

数学符号说明

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

  1. ∣ ∣ X ∣ ∣ p \left|\left| \boldsymbol{\mathcal{X}}\right|\right|_p ∣∣X∣∣p​表示张量的 p p p范数，缺省状态下取 p = 2 p=2 p=2
    ∣ ∣ X ∣ ∣ p = ( ∑ i ∑ j ∑ k ∣ x i j k ∣ p ) 1 p \left|\left| \boldsymbol{\mathcal{X}}\right|\right|_p= \left(\sum_i\sum_j\sum_k{\left| x_{ijk}\right|^p}\right)^{\frac{1}{p}} ∣∣X∣∣p​=(i∑​j∑​k∑​∣xijk​∣p)p1​

  2. < X ， Y > \left<\boldsymbol{\mathcal{X}}，\boldsymbol{\mathcal{Y}}\right> ⟨X，Y⟩表示张量的内积
    < X ， Y > = ∑ i ∑ j ∑ k x i j k y i j k \left<\boldsymbol{\mathcal{X}}，\boldsymbol{\mathcal{Y}}\right>= \sum_i\sum_j\sum_k{ x_{ijk}y_{ijk}} ⟨X，Y⟩=i∑​j∑​k∑​xijk​yijk​

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

    z i j k = x i j k y i j k z_{ijk}=x_{ijk}y_{ijk} zijk​=xijk​yijk​

  4. A ⊗ B \bold{A}\otimes\bold{B} A⊗B表示Kronecker积，若 A ∈ R I × J ， Y ∈ R M × N \bold{A}\in\mathbb{R}^{I\times J}，\boldsymbol{Y}\in\mathbb{R}^{M\times N} A∈RI×J，Y∈RM×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，计算方式如下：
    A ⊗ B = [ a 11 B a 12 B ⋯ a 1 m 2 B a 21 B a 22 B ⋯ a 2 m 2 B ⋮ ⋮ ⋱ ⋮ a m 1 1 B + + a m 1 2 B ⋯ a m 1 m 2 B ] \bold{A}\otimes \bold{B}= \left[\begin{array}{cccc} a_{11}\bold{B}&a_{12}\bold{B}&\cdots&a_{1m_2}\bold{B}\\ a_{21}\bold{B}&a_{22}\bold{B}&\cdots&a_{2m_2}\bold{B}\\ \vdots&\vdots&\ddots&\vdots\\ a_{m_11}\bold{B}+&+a_{m_12}\bold{B}&\cdots&a_{m_1m_2}\bold{B} \end{array}\right] A⊗B=⎣⎢⎢⎢⎡​a11​Ba21​B⋮am1​1​B+​a12​Ba22​B⋮+am1​2​B​⋯⋯⋱⋯​a1m2​​Ba2m2​​B⋮am1​m2​​B​⎦⎥⎥⎥⎤​

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

    A [ 1 ] ⊙ A [ 2 ] ⊙ ⋯ ⊙ A [ M ] = ˙ ⨀ m = 1 M A [ m ] \boldsymbol{A_{[1]}}\odot\boldsymbol{A_{[2]}}\odot\cdots\odot\boldsymbol{A_{[M]}} \dot{=}\bigodot^M_{m=1}\boldsymbol{A_{[m]}} A[1]​⊙A[2]​⊙⋯⊙A[M]​=˙m=1⨀M​A[m]​

  6. a ∘ b \bold{a}\circ\bold{b} a∘b表示向量的外积
    a ∘ b = a b T \bold{a}\circ\bold{b}=\bold{a}\bold{b}^T a∘b=abT

  7. X × ˉ n Y \boldsymbol{\mathcal{X}}\bar{\times}_n\boldsymbol{Y} X×ˉ​n​Y，表示张量 X \boldsymbol{\mathcal{X}} X与矩阵 Y \boldsymbol{Y} Y的 n n n模态积，若 X ∈ R I 1 × I 2 × ⋯ × I N ， Y ∈ R m × I n \boldsymbol{\mathcal{X}}\in\mathbb{R}^{I_1\times I_2\times\dots\times I_N}，\boldsymbol{Y}\in\mathbb{R}^{m\times I_n} X∈RI1​×I2​×⋯×IN​，Y∈Rm×In​，则 { 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​
    X × ˉ 1 Y [ 1 ] × ˉ 2 Y [ 2 ] × ˉ 3 ⋯ × ˉ M Y [ M ] = X ∏ m = 1 M × ˉ m Y [ m ] \boldsymbol{\mathcal{X}}\bar{\times}_1\boldsymbol{Y_{[1]}}\bar{\times}_2\boldsymbol{Y_{[2]}}\bar{\times}_3\cdots\bar{\times}_M\boldsymbol{Y_{[M]}}= \boldsymbol{\mathcal{X}}\prod_{m=1}^M\bar{\times}_mY_{[m]} X×ˉ​1​Y[1]​×ˉ​2​Y[2]​×ˉ​3​⋯×ˉ​M​Y[M]​=Xm=1∏M​×ˉ​m​Y[m]​

  8. X † \boldsymbol{X}^{\dagger} X†表示矩阵 X \boldsymbol{X} X的广义逆（Moore–Penrose pseudoinverse）

  9. X ( n ) \boldsymbol{\mathcal{X}}_{(n)} X(n)​表示张量以第 n n n模态（ n − m o d e n-mode n−mode）展开得到的矩阵，若 X ∈ R 4 × 3 × 2 \boldsymbol{\mathcal{X}}\in\mathbb{R}^{4\times 3\times 2} X∈R4×3×2，则
    X ( 1 ) = [ X ( : , : , 1 ) , X ( : , : , 2 ) ] X ( 2 ) = [ X ( : , : , 1 ) T , X ( : , : , 2 ) T ] X ( 3 ) = [ X ( : , 1 , : ) T , X ( : , 2 , : ) T , X ( : , 3 , : ) T ] {\mathcal{X}}_{\left(1\right)}=\Big[{\mathcal{X}}\left(:,:,1\right),{\mathcal{X}}\left(:,:,2\right)\Big]\\ {\mathcal{X}}_{\left(2\right)}=\Big[{\mathcal{X}}\left(:,:,1\right)^T,{\mathcal{X}}\left(:,:,2\right)^T\Big]\\ {\mathcal{X}}_{\left(3\right)}=\Big[{\mathcal{X}}\left(:,1,:\right)^T,{\mathcal{X}}\left(:,2,:\right)^T, {\mathcal{X}}\left(:,3,:\right)^T\Big] X(1)​=[X(:,:,1),X(:,:,2)]X(2)​=[X(:,:,1)T,X(:,:,2)T]X(3)​=[X(:,1,:)T,X(:,2,:)T,X(:,3,:)T]

CP分解

基本概念

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

设 X ∈ R I × J × K \boldsymbol{\mathcal{X}}\in\mathbb{R}^{I\times J\times K} X∈RI×J×K， A = ( a 1 , a 2 , … , a R ) ∈ R I \bold{A}=\left(\bold{a}_1,\bold{a}_2,\dots,\bold{a}_R\right)\in\mathbb{R}^I A=(a1​,a2​,…,aR​)∈RI， B = ( b 1 , b 2 , … , b R ) ∈ R J \bold{B}=\left(\bold{b}_1,\bold{b}_2,\dots,\bold{b}_R\right)\in\mathbb{R}^J B=(b1​,b2​,…,bR​)∈RJ， C = ( c 1 , c 2 , … , c R ) ∈ R K \bold{C}=\left(\bold{c}_1,\bold{c}_2,\dots,\bold{c}_R\right)\in\mathbb{R}^K C=(c1​,c2​,…,cR​)∈RK，则CP分解可以公式化为：
X ≈ ∑ r = 1 R a r ∘ b r ∘ c r \boldsymbol{\mathcal{X}}\approx \sum_{r=1}^R \bold{a}_r\circ\bold{b}_r\circ\bold{c}_r X≈r=1∑R​ar​∘br​∘cr​
CP分解公式还可以进一步矩阵化为：
X ( 1 ) ≈ A ( C ⊙ B ) T \boldsymbol{\mathcal{X}}_{(1)}\approx \bold{A}\left(\bold{C}\odot\bold{B}\right)^T X(1)​≈A(C⊙B)T

X ( 2 ) ≈ B ( C ⊙ A ) T \boldsymbol{\mathcal{X}}_{(2)}\approx \bold{B}\left(\bold{C}\odot\bold{A}\right)^T X(2)​≈B(C⊙A)T

X ( 3 ) ≈ C ( B ⊙ A ) T \boldsymbol{\mathcal{X}}_{(3)}\approx \bold{C}\left(\bold{B}\odot\bold{A}\right)^T X(3)​≈C(B⊙A)T

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

在实际进行CP分解时，通常还会对 A 、 B 、 C \bold{A}、\bold{B}、\bold{C} A、B、C按列归一化，并将其范数作为权重存储为向量 λ ∈ R R \boldsymbol{\lambda}\in\mathbb{R}^R λ∈RR，则新的公式为：
X ≈ ∑ r = 1 R λ r a r ∘ b r ∘ c r \boldsymbol{\mathcal{X}}\approx \sum_{r=1}^R \boldsymbol{\lambda}_r\bold{a}_r\circ\bold{b}_r\circ\bold{c}_r X≈r=1∑R​λr​ar​∘br​∘cr​
对于更高阶的张量 X ∈ R I 1 × I 2 × ⋯ × I N \boldsymbol{\mathcal{X}}\in\mathbb{R}^{I_1\times I_2\times\cdots\times I_N} X∈RI1​×I2​×⋯×IN​，通用表达式为：
X ( n ) ≈ A [ n ] Λ ( A [ N ] ⊙ ⋯ ⊙ A [ n + 1 ] ⊙ A [ n − 1 ] ⊙ ⋯ ⊙ A [ 1 ] ) T \boldsymbol{\mathcal{X}}_{(n)}\approx \bold{A}_{[n]}\bold{\Lambda}\left( \bold{A}_{[N]}\odot\cdots\odot\bold{A}_{[n+1]}\odot\bold{A}_{[n-1]}\odot\cdots\odot\bold{A}_{[1]} \right)^T X(n)​≈A[n]​Λ(A[N]​⊙⋯⊙A[n+1]​⊙A[n−1]​⊙⋯⊙A[1]​)T
其中 Λ = d i a g ( λ ) \bold{\Lambda}=diag(\boldsymbol{\lambda}) Λ=diag(λ)， d i a g diag diag为对角矩阵化。

ALS-CP

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

公式推导

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

A 、 B 、 C \bold{A}、\bold{B}、\bold{C} A、B、C估计得到的张量 X ^ \hat{\boldsymbol{\mathcal{X}}} X^尽可能的接近原始张量 X \boldsymbol{\mathcal{X}} X，即：
m i n ∣ ∣ X ^ − X ∣ ∣ with  X ^ = ∑ r = 1 R λ r a r ∘ b r ∘ c r min\left|\left| \hat{\boldsymbol{\mathcal{X}}}-\boldsymbol{\mathcal{X}} \right|\right| \text{ with } \hat{\boldsymbol{\mathcal{X}}}= \sum_{r=1}^R \boldsymbol{\lambda}_r\bold{a}_r\circ\bold{b}_r\circ\bold{c}_r min∣∣∣​∣∣∣​X^−X∣∣∣​∣∣∣​ with X^=r=1∑R​λr​ar​∘br​∘cr​
为此，我们利用公式(12)(16)可以得到如下估计因子矩阵 A ^ = A Λ \hat{\bold{A}}=\bold{A}\bold{\Lambda} A^=AΛ的方法：
m i n ∣ ∣ X ^ ( 1 ) − X ( 1 ) ∣ ∣ = m i n ∣ ∣ A ^ ( C ⊙ B ) T − X ( 1 ) ∣ ∣ min\left|\left| \hat{\boldsymbol{\mathcal{X}}}_{(1)}-\boldsymbol{\mathcal{X}}_{(1)} \right|\right|= min\left|\left| \hat{\bold{A}}\left(\bold{C}\odot\bold{B}\right)^T -\boldsymbol{\mathcal{X}_{(1)}} \right|\right| min∣∣∣​∣∣∣​X^(1)​−X(1)​∣∣∣​∣∣∣​=min∣∣∣​∣∣∣​A^(C⊙B)T−X(1)​∣∣∣​∣∣∣​

A ^ = X ( 1 ) [ ( C ⊙ B ) T ] † \hat{\bold{A}}=\boldsymbol{\mathcal{X}_{(1)}}\left[\left(\bold{C}\odot\bold{B}\right)^T\right]^{\dagger} A^=X(1)​[(C⊙B)T]†

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

C ^ = X ( 2 ) ( B ⊙ A ) ( B T B ∗ A T A ) † \hat{\bold{C}}=\boldsymbol{\mathcal{X}_{(2)}} \left(\bold{B}\odot\bold{A}\right) \left(\bold{B}^T\bold{B}\ast\bold{A}^T\bold{A}\right)^{\dagger} C^=X(2)​(B⊙A)(BTB∗ATA)†

对高阶张量 X ∈ R I 1 × I 2 × ⋯ × I N \boldsymbol{\mathcal{X}}\in\mathbb{R}^{I_1\times I_2\times\cdots\times I_N} X∈RI1​×I2​×⋯×IN​以及参数矩阵集合 { A [ n ] ∈ R I n × R } n N \left\{\bold{A}_{[n]}\in\mathbb{R}^{I_n\times R}\right\}_n^N {A[n]​∈RIn​×R}nN​的通用表达式为:
V = ( A [ N ] T A [ N ] ∗ ⋯ A [ n + 1 ] T A [ n + 1 ] ∗ A [ n − 1 ] T A [ n − 1 ] ∗ ⋯ ∗ A [ 1 ] T A [ 1 ] ) \bold{V}= \left( \bold{A}_{[N]}^T\bold{A}_{[N]}\ast\cdots\bold{A}_{[n+1]}^T\bold{A}_{[n+1]}\ast\bold{A}_{[n-1]}^T\bold{A}_{[n-1]}\ast\cdots\ast\bold{A}_{[1]}^T\bold{A}_{[1]} \right) V=(A[N]T​A[N]​∗⋯A[n+1]T​A[n+1]​∗A[n−1]T​A[n−1]​∗⋯∗A[1]T​A[1]​)

A ^ [ n ] = X ( n ) ( A [ N ] ⊙ ⋯ ⊙ A [ n + 1 ] ⊙ A [ n − 1 ] ⊙ ⋯ ⊙ A [ 1 ] ) V † \hat{\bold{A}}_{[n]}=\boldsymbol{\mathcal{X}}_{(n)} \left( \bold{A}_{[N]}\odot\cdots\odot\bold{A}_{[n+1]}\odot\bold{A}_{[n-1]} \odot\cdots\odot\bold{A}_{[1]} \right)\bold{V}^{\dagger} A^[n]​=X(n)​(A[N]​⊙⋯⊙A[n+1]​⊙A[n−1]​⊙⋯⊙A[1]​)V†

程序实现

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

  • Step1: 初始化参数矩阵集合 { A [ n ] ∈ R I n × R } n N = A \left\{\bold{A}_{[n]}\in\mathbb{R}^{I_n\times R}\right\}_n^N=\text{A} {A[n]​∈RIn​×R}nN​=A以及权重向量 λ = lbd \boldsymbol{\lambda}=\text{lbd} λ=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)计算 V \bold{V} 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)计算 A ^ [ n ] = A[n] \hat{\bold{A}}_{[n]}=\text{A[n]} A^[n]​=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 = ( A [ N ] ⊙ ⋯ ⊙ A [ n + 1 ] ⊙ A [ n − 1 ] ⊙ ⋯ ⊙ A [ 1 ] ) T=\left( \bold{A}_{[N]}\odot\cdots\odot\bold{A}_{[n+1]}\odot\bold{A}_{[n-1]} \odot\cdots\odot\bold{A}_{[1]} \right) T=(A[N]​⊙⋯⊙A[n+1]​⊙A[n−1]​⊙⋯⊙A[1]​)

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

V † = np.linalg.pinv(V) \bold{V}^{\dagger}=\text{np.linalg.pinv(V)} V†=np.linalg.pinv(V)

  • Step4: 对 A ^ [ n ] \hat{\bold{A}}_{[n]} A^[n]​的每一列做归一化得到权重向量 λ \boldsymbol{\lambda} λ以及 A [ n ] \bold{A}_{[n]} 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

l o s s = ∣ ∣ X ^ − X ∣ ∣ = tl.norm(tensor - tensor_pred) loss=\left|\left| \hat{\boldsymbol{\mathcal{X}}}-\boldsymbol{\mathcal{X}} \right|\right| =\text{tl.norm(tensor - tensor\_pred)} loss=∣∣∣​∣∣∣​X^−X∣∣∣​∣∣∣​=tl.norm(tensor - tensor_pred)

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

Λ = np.diag(lbd) \bold{\Lambda}=\text{np.diag(lbd)} Λ=np.diag(lbd)

完整程序如下：

import numpy as np
import tensorly as tl
from tensorly.decomposition._cp import initialize_cpfrom tensorly.tenalg import khatri_rao
from tqdm import tqdmdef cp_als(tensor: np.ndarray, R=1, max_iter=100):N = tl.ndim(tensor)# Step 1lbd, 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 2V = 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 3T = khatri_rao(A, skip_matrix=n)A[n] = np.matmul(np.matmul(tl.unfold(tensor, mode=n), T), np.linalg.pinv(V))# Step 4for r in range(R):lbd[r] = tl.norm(A[n][:, r])A[n] = A[n] / tl.reshape(lbd, (1, -1))# Step 5tensor_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, epochreturn A, lbd, max_iterif __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. 为了初始化 A [ n ] \bold{A}_{[n]} A[n]​，需要先对 X ( n ) \boldsymbol{\mathcal{X}}_{(n)} X(n)​做SVD得到左奇异矩阵 U ∈ R M × N \bold{U}\in\mathbb{R}^{M\times N} U∈RM×N:

X ( n ) = U Σ V \boldsymbol{\mathcal{X}}_{(n)}=\bold{U}\bold{\Sigma}\bold{V} X(n)​=UΣV

  1. 若选取的 R ≤ N \bold{R}\leq N R≤N，取左奇异矩阵 U \bold{U} U的前 R \bold{R} R列作为 A [ n ] \bold{A}_{[n]} A[n]​即可：

A [ 1 ] = [ u : 1 , u : 2 , … , u : R ] \bold{A}_{[1]}=\left[\bold{u}_{:1},\bold{u}_{:2},\dots,\bold{u}_{:R}\right] A[1]​=[u:1​,u:2​,…,u:R​]

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

BP-CP

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

公式推导

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

Θ = X ( n ) − A [ n ] ( A [ N ] ⊙ ⋯ ⊙ A [ n + 1 ] ⊙ A [ n − 1 ] ⊙ ⋯ ⊙ A [ 1 ] ) T \bold{\Theta}=\boldsymbol{\mathcal{X}}_{(n)}-\bold{A}_{[n]}\left(\bold{A}_{[N]}\odot\cdots\odot\bold{A}_{[n+1]}\odot\bold{A}_{[n-1]} \odot\cdots\odot\bold{A}_{[1]}\right)^T Θ=X(n)​−A[n]​(A[N]​⊙⋯⊙A[n+1]​⊙A[n−1]​⊙⋯⊙A[1]​)T

参数更新公式为：
A [ n ] ← A [ n ] − l r ∗ ∂ L ∂ A [ n ] \bold{A}_{[n]}\leftarrow \bold{A}_{[n]}-lr\ast\frac{\partial\bold{L}}{\partial\bold{A}_{[n]}} A[n]​←A[n]​−lr∗∂A[n]​∂L​

程序实现

  • Step1: 初始化参数矩阵集合 { A [ n ] ∈ R I n × R } n N = A \left\{\bold{A}_{[n]}\in\mathbb{R}^{I_n\times R}\right\}_n^N=\text{A} {A[n]​∈RIn​×R}nN​=A、梯度矩阵（Jacobian矩阵）集合 { L A [ n ] ∈ R I n × R } n N = grad_A \left\{\frac{\bold{L}}{\bold{A}_{[n]}}\in\mathbb{R}^{I_n\times R}\right\}_n^N=\text{grad\_A} {A[n]​L​∈RIn​×R}nN​=grad_A以及权重向量 λ = lbd \boldsymbol{\lambda}=\text{lbd} λ=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 \Theta=\text{theta} Θ=theta，并将梯度矩阵（Jacobian矩阵）集合 { L A [ n ] ∈ R I n × R } n N = grad_A \left\{\frac{\bold{L}}{\bold{A}_{[n]}}\in\mathbb{R}^{I_n\times R}\right\}_n^N=\text{grad\_A} {A[n]​L​∈RIn​×R}nN​=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: 计算每个参数矩阵 A [ n ] \bold{A}_{[n]} A[n]​对应的梯度矩阵 L A [ n ] \frac{\bold{L}}{\bold{A}_{[n]}} A[n]​L​并存储，对应公式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的损失值并输出，方便观察收敛情况，这里对 L \bold{L} 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_cpfrom tensorly.tenalg import khatri_rao
from tqdm import tqdmdef cp_bp(tensor: np.ndarray, R=1, lr=1e-2, max_iter=100):N = tl.ndim(tensor)# Step 1lbd, A = initialize_cp(tensor, R, init='random', svd='numpy_svd',random_state=0,normalize_factors=True)for epoch in range(max_iter):# Step 2tensor_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 3for 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 4for n in range(N):A[n] = A[n] + lr * grad_A[n]# Step 5loss = np.sum(0.5 * np.square(tl.unfold(theta, 0)))print("epoch {}: loss={}".format(epoch, loss))return A, lbdif __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)

张量CP分解原理及Python实现相关推荐

  1. 『矩阵论笔记』张量CP分解的详细数学推导以及Python实现

    张量CP分解的详细数学推导以及Python实现! 文章目录 一. 张量的基本概念 1.1 张量定义以及张量空间 1.2 阶和纤维(fiber)及切片 1.3 内积和范数及秩一张量/可合张量 1.4 超 ...

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

    本文是对论文Tensor Decompositions and Applications进行了翻译.整理.筛选和适当的补充,如何希望深入理解可以阅读原文. 相关文章: [张量分解(一)]符号与基础知识 ...

  3. 浅析张量分解（Tensor Decomposition）

    一般一维数组,我们称之为向量(vector),二维数组,我们称之为矩阵(matrix);三维数组以及多位数组,我们称之为张量(tensor). 在介绍张量分解前,我们先看看矩阵分解相关知识概念. 一. ...

  4. 从15000个Python开源项目中精选的Top30，Github平均star为3707，赶紧收藏！

    翻译 | AI科技大本营(ID:rgznai100) 参与 | SuiSui 继推出2017年机器学习开源项目Top 30榜单后,Mybridge AI又推出了一个Python开源项目Top 30榜单 ...

  5. 收藏！15000个Python开源项目中精选Top30！

    来源:授权自AI科技大本营(ID:rgznai100) 本文长度为1700字,建议阅读6分钟 本文基于项目质量.用户参与度以及其他因素为你列出Python开源项目Top 30,建议收藏. 继推出201 ...

  6. 从Python调用外部命令

    您如何在Python脚本中调用外部命令(就像我在Unix Shell或Windows命令提示符下键入的一样)? #1楼 os.system不允许您存储结果,因此,如果您要将结果存储在某个列表中或sub ...

  7. mac如何导入python第三方库_Mac系统中python idle导入第三方模块成功，ecplise导入python第三方模块失败解决方法...

    遇到一个比较纠结了4个月的问题,一直没有在意,今天实在忍受不了,尝试各种解决办法,终于把这个烦人的问题完美解决,不敢独享,写出来和各位大神共享. 问题:在mac OSx操作系统下,安装了python第 ...

  8. 计算成本缩减100倍！港中文提出语义分割新方法：张量低秩重建｜ECCV2020

    原文链接:https://bbs.cvmart.net/articles/3099 专注计算机视觉前沿资讯和技术干货 微信公众号:极市平台 官网:https://www.cvmart.net/ --- ...

  9. 怎么查看python是多少位_python+位数

    广告关闭 腾讯云11.11云上盛惠 ,精选热门产品助力上云,云服务器首年88元起,买的越多返的越多,最高返5000元! 程序中的所有数在计算机内存中都是以二进制的形式储存的. 位运算就是直接对整数在内 ...

最新文章

  1. 全面解读：腾讯 CDB 内核特性与优化实践
  2. android yuv加水印_Android Camera添加预览水印
  3. 防火墙工作原理—Vecloud微云
  4. 微软中国发布“IE8浏览器性能解密”，为金山网盾辟谣
  5. 让PIP源使用国内镜像，提升下载速度和安装成功率。
  6. another rejection from Cambridge MPhil in Management
  7. 计算机更改家庭组密码,Win10系统怎么修改家庭组密码 win10修改家庭组密码的方法...
  8. RocketMQ的Producer详解之分布式事务消息（原理分析）
  9. SpringMVC 控制器默认支持GET和POST两种方式
  10. Delphi 的一些函数（Windows相关）
  11. wpf之界面控件MaterialDesignInXAML
  12. python定义常量、装饰器_Python 装饰器
  13. 利润表模板excel_年薪60w财务总监：工作八年，这10个Excel必备财务系统，效率翻倍...
  14. 20sccm_SCCM 2016 使用PXE 部署操作系统（一）
  15. 图片怎样编辑文字？分享三个图片编辑修改文字的方法
  16. HBuilderX、微信开发者工具、VScode之间运行微信公众号
  17. C# 以MP3的格式将录制的音频数据写入文件流
  18. 自己总结关于浏览器证书安全的二点小技巧
  19. [2019长沙长郡中学集训]加法
  20. 组建计算机网络目的三个,计算机网络概述（其一）

热门文章

  1. 微信公众号{errcode:40164,errmsg:invalid ip 113.70.100.150, not in whitelist hint: [i50PvA0089sha6]}
  2. Python装逼神器，5 行 Python 代码 实现一键批量扣图
  3. 在Scrum开发模式下，为Sprint起名字的艺术
  4. SpringBoot系列教程JPA之指定id保存
  5. git 本地仓库关联到远程仓库
  6. 新型冠状病毒个人防护报告：如何过个安全健康年？
  7. 算法交易系统的数据库之路
  8. matlab fis编辑器在哪,基本FIS编辑器
  9. Python基础编程：利用列表实现简单的先进后出、先进先出
  10. Latex中ACM-Reference-Format顺序与论文引用顺序不一致solution