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

张量CP分解的详细数学推导以及Python实现！

文章目录

• ​​一. 张量的基本概念​​

• ​​1.1 张量定义以及张量空间​​
• ​​1.2 阶和纤维(fiber)及切片​​
• ​​1.3 内积和范数及秩一张量/可合张量​​
• ​​1.4 超对称和超对角及展开(unfolding/flattening)​​
• ​​1.5 n-mode(矩阵)乘积和n-mode(向量)乘积​​
• ​​1.6 矩阵的Kronecker乘积及Khatri-Rao乘积​​
• ​​1.7 矩阵的Hadamard乘积​​

• ​​二. 张量的CP分解​​

• ​​2.1. CP分解的表示​​
• ​​2.2. 交替最小二乘法求解​​
• ​​2.3. 交替最小二乘法Python代码​​
• ​​2.4. 梯度下降法求解​​
• ​​2.5. 梯度下降法Python代码​​

• ​​三. 拓展应用-数据修补​​

• ​​3.1. 矩阵分解(Matrix factorization)​​
• ​​3.2. 张量分解(Tensor factorization)​​

• ​​四. 参考文献​​

• 强烈推荐该PPT(特详细)：张量分解_张量CP分解_张量Tucker分解_详细介绍：​​javascript:void(0)​​
• 本文使用Python实现(欢迎各位朋友star,接下来还会继续更新项目!)，代码地址：​​https://github.com/zhangkaifang/cp_decomposition​​

一. 张量的基本概念

1.1 张量定义以及张量空间

• 张量(tensor)？简单地说，就是个多维数组。在本研究范围内，不考虑任何物理和工学领域内的张量定义，而仅仅考虑其数学领域。正式的说，应该叫张量域(tensor fields)。第一阶张量(first-order tensor)是个向量(vector)，第二阶张量(second-order tensor)是矩阵(matrix)，更多阶的张量我们称之为高阶张量(higher-order tensor)。

• 由若干个向量空间中的基底的外积张成的空间。

1.2 阶和纤维(fiber)及切片

• 阶​​(order/ways/modes/rank)​​；
• 张成所属张量空间的向量空间的个数；
• ；
• ；
• ；
• ；

• 纤维(fiber)表示从张量中抽取向量的操作， 规定其它维度只保持一个维度变化。

• 左前上角索引起始点。

1.3 内积和范数及秩一张量/可合张量

• 。
• 内积：
  
• (Frobenius)范数：
  
• 阶张量 是一个​​​秩一张量​​​，如果它能被写成 个向量的外积，即：
  

• 例子： 给定如下的两个张量

• 内积：
  
• 张量
  

1.4 超对称和超对角及展开(unfolding/flattening)

• 立方张量：各个​​mode​​的长度相等；
• 对称： 一个立方张量是对称的，如果其元素在下标的任意排列下是常数。如一个三阶立方张量是超对称的，如果
• 对角： 仅当 时，；

• 将 阶张量 沿​​​mode-n​​​展开成一个矩阵 。 (​​​mode-1​​​, 对应着张量的第一阶；​​mode-2,​​​ 对应着张量的第二阶；​​mode-3​​ 对应着张量的第三阶)。

• 例子： ​​如何简单地理解和实现「张量展开」？​​

1.5 n-mode(矩阵)乘积和n-mode(向量)乘积

• 张量与矩阵相乘(又称为模态积)相比矩阵与矩阵之间的相乘更为抽象，如何理解呢？
• 一个张量 和一个矩阵 ，则张量和矩阵的 模态积(-mode product)记为 ，其大小为。对于每个元素而言，有：
  其中 ，我们可以看出，模态积是张量、矩阵和模态（mode）的一种“组合”运算。另外
  
• 例子：上述给出张量与矩阵相乘的定义，为了方便理解，下面来看一个简单的示例，若给定张量 为： 其大小为 。另外给定矩阵 试想一下：张量 和矩阵 相乘会得到什么呢？，则张量 任意索引位置 上的值为： 这一运算规则也不难发现，张量 的大小为 ，以 位置为例：
  再以 位置为例：
  
• 也就是如下：
  
• （会得到大小为 的张量）或 （会得到大小为
• 需要注意的是， 有一个恒等的计算公式，即 ，由于 则 满足 ，即采用张量矩阵化的形式进行运算可以使问题变得更加简单，从这里也可以看出高阶张量进行矩阵化的优点。
• 性质：
• 例子：利用上一小节中的张量乘以矩阵
  最终组合成张量就是如下：
• 因此，乘法过程可用下图表示：先将张量矩阵化，再将张量和矩阵相乘。不同的mode-n矩阵化会使得相乘结果不同。

• 一个张量 和一个向量 的 ​​n-mode​​ 乘积，其元素定义为：
  
• 性质：
• 例子：利用上一小节中的张量乘以向量
  最终组合成张量就是如下：
• 很好的一篇文章：​​浅谈张量分解（二）：张量分解的数学基础​​

1.6 矩阵的Kronecker乘积及Khatri-Rao乘积

• 好文章推荐：​​代数基础 | Kronecker积​​
• 矩阵的​​Kronecker​​乘积
• 矩阵 ，则
  其中
• 性质：
  
• 性质： 矩阵的​​Kronecker​​​积还和张量和矩阵的​​n-mode​​​乘积有如下关系
• 例子：
  

矩阵的​​Khatri-Rao​​乘积

• (保证列相同)，则
  

• 其中： 都是列向量。 注意： 如果 都是向量，那么​​​Khatri-Rao​​​乘积和​​Kronecker​​​乘积相等，即：。
• 性质： 矩阵的​​Kronecker​​​积还和张量和矩阵的​​n-mode​​​乘积有如下关系
  
• 例子：
  
• 程序例子：

import scipy.linalg as lin
import numpy as np
# 自己实现
def kr_prod(a, b):
    return np.einsum('ir, jr -> ijr', a, b).reshape(a.shape[0] * b.shape[0], -1)

a = np.array([[1, 2],
              [4, 5]])
b = np.array([[10, 20],
              [30, 40],
              [50, 60]])
lin.kron(a, b)  
lin.khatri_rao(a, b)
# lin.kron(a, b)  
# array([[ 10,  20,  20,  40],
#        [ 30,  40,  60,  80],
#        [ 50,  60, 100, 120],
#        [ 40,  80,  50, 100],
#        [120, 160, 150, 200],
#        [200, 240, 250, 300]])
# # lin.khatri_rao(a, b) 
# array([[ 10,  40],
#        [ 30,  80],
#        [ 50, 120],
#        [ 40, 100],
#        [120, 200],
#        [200, 300]])

1.7 矩阵的Hadamard乘积

• 矩阵 ，则
• 性质：
  

二. 张量的CP分解

2.1. CP分解的表示

• CP分解的张量形式：
• 设 。将一个张量表示成有限个秩一张量之和，比如一个​​三阶张量​​​可以分解为：
  
• 表示向量的外积，其中向量 是因子矩阵 的第 列，并且向量 和 有相同的定义因子矩阵 和 中。 事实上，这些向量的外积是一个秩一张量，因此，我们可以用
• 元素方面，对于张量 中的任何第 元素，我们有：
  
• , , and 。

• CP组合示例：
• 给定矩阵 , 以及 。如果 ，我们有如下：
  
• 注意是近似：

• 

import numpy as np
def cp_combine(mat1, mat2, mat3):
    return np.einsum('ir, jr, tr -> ijt', mat1, mat2, mat3)

A = np.array([[1, 2], [3, 4]])
B = np.array([[1, 3], [2, 4], [5, 6]])
C = np.array([[1, 5], [2, 6], [3, 7], [4, 8]])
print(cp_combine(A, B, C))
print()
print('tensor size:')
print(cp_combine(A, B, C).shape)

[[[ 31  38  45  52]
  [ 42  52  62  72]
  [ 65  82  99 116]]

 [[ 63  78  93 108]
  [ 86 108 130 152]
  [135 174 213 252]]]

tensor size:
(2, 3, 4)

• CP分解的矩阵形式：
• 因子矩阵： 秩一张量中对应的向量组成的矩阵，如
  
• 利用因子矩阵，一个三阶张量的CP分解可以写成展开形式：
  
• 也就是得到如下近似：
  
• 与矩阵类似，​​张量也存在秩​​​的概念，使 时最小的 值，或者说能够完美拟合原始张量所需的最少秩1张量数量，即为张量
• 在实际进行​​CP分解​​​时，通常还会对 按列归一化，并将其范数作为权重存储为向量 ，则新的公式为：
  
• 对于更高阶的张量
  其中 ， diag表示对角矩阵。

2.2. 交替最小二乘法求解

• 交替最小二乘法(alternating least squares，ALS)
• 下面以三阶张量 为例介绍如何使用​​​ALS​​​计算​​CP分解​​​的因子矩阵 。首先我们的最终目标是让由 估计得到的张量张量 尽可能的接近原始张量 ，即：
  
• 首先固定 矩阵，来更新 矩阵，首先得到因子矩阵 ：
  
• 其中 表示张量 的模式 1 展开为矩阵，依此类推。 (????⊙????)???? 表示将 B 和 C 组合成单个矩阵的 Khatri-Rao 乘积。 一般来说，这是一个非凸问题； 然而，当我们当时优化一个矩阵时，这是一个凸问题。
  这样可以得到 ，这里 表示矩阵 的​​​广义逆(Moore–Penrose pseudoinverse)​​。
• 根据​​Khatri-Rao的性质​​​，这里就是参考公式 。最终可以得到：
  
• 类似地可以更新
  
• 对高阶张量 以及参数矩阵集合
  
• 这里参考一些这篇博客：​​张量CP分解原理及Python实现！​​

2.3. 交替最小二乘法Python代码

#!/usr/bin/env python
# -*- coding: UTF-8 -*-
'''=====================================
@Author :Kaifang Zhang
@Time   :2021/7/5 1:31
@Contact: [e-mail]
========================================'''
import os
import numpy as np
import scipy.linalg as lin
import tensorly as tl

np.random.seed(2021)


def col_normalize(X):
    return X / np.sqrt(np.sum(X ** 2, axis=0))


def tensor2matrix(X, mode):
    """
    X shape: (40, 80, 120)
    Return mu-mode matricization from a given tensor
    """
    num_dim = len(X.shape)
    n = X.shape[num_dim - mode]
    X = np.moveaxis(X, num_dim - mode, -1)
    return np.reshape(X, (-1, n)).T


def matrix2tensor(X1, out_shape):
    """
    Input: 1-mode matricization
    Output: tensor with size like (n3, n2, n1)
    """
    return np.reshape(X1.T, out_shape)


def ALS_solver(X, r, nmax=1000, err_tol=1e-4):
    """
    Parameters
    X : tensor like B1
    r : tensor rank
    nmax : maximum number of iterations
    err_tol : tolerance for relative residual error, optional
    Returns：approximated tensor with same shape as X
    """
    ############################# 第一种 ###############################
    # n1, n2, n3 = X.shape
    # B = np.random.random((n2, r))
    # C = np.random.random((n3, r))
    # X1 = tl.unfold(X, mode=0)  # (4, 96)这里调用tensorly包里面的，结果是一样的。
    # X2 = tl.unfold(X, mode=1)  # (8, 48)这里调用tensorly包里面的，结果是一样的。
    # X3 = tl.unfold(X, mode=2)  # (12, 32)这里调用tensorly包里面的，结果是一样的。

    # B = np.random.normal(0, 1, (n2, r))
    # C = np.random.normal(0, 1, (n3, r)) # 这种和下面都是可以的。
    # X1 = tensor2matrix(X, 3)
    # X2 = tensor2matrix(X, 2)
    # X3 = tensor2matrix(X, 1)
    ############################# 第二种 ###############################
    n3, n2, n1 = X.shape
    B = np.random.normal(0, 1, (n2, r))
    C = np.random.normal(0, 1, (n3, r)) # 这种和下面都是可以的。
    X1 = tensor2matrix(X, 1)
    X2 = tensor2matrix(X, 2)
    X3 = tensor2matrix(X, 3)

    X_norm = lin.norm(X1, 'fro')
    err = np.inf
    B = col_normalize(B)
    i = 0
    while (err >= err_tol) and i < nmax:
        C = col_normalize(C)
        tem1 = lin.khatri_rao(C, B)
        A, res, rnk, s = lin.lstsq(tem1, X1.T)
        A = A.T

        A = col_normalize(A)
        tem2 = lin.khatri_rao(C, A)
        B, res, rnk, s = lin.lstsq(tem2, X2.T)
        B = B.T

        B = col_normalize(B)
        tem3 = lin.khatri_rao(B, A)
        C, res, rnk, s = lin.lstsq(tem3, X3.T)
        C = C.T

        X_hat1 = A.dot(lin.khatri_rao(C, B).T)
        err = lin.norm(X_hat1 - X1, 'fro') / X_norm
        i += 1
        print('Relative error at iteration ', i, ': ', err)
    X_hat = matrix2tensor(X_hat1, X.shape)
    print('Finished!')
    return A, B, C, X_hat


if __name__ == "__main__":
    tensor = np.random.rand(4, 8, 12)
    A, B, C, X_hat = ALS_solver(tensor, r=4)

• 解释​​np.moveaxis​​函数：

2.4. 梯度下降法求解

• 反向传播(Back Propagation，BP)，利用梯度下降法求解CP参数 。
• 估计得到的张量 尽可能的接近原始张量 ，因此将损失/目标函数设置为：
  
• 乘以 是为了方便求导,为了简化表达式，设 ，则公式 中的 ，则对矩阵 求偏导数可以得到：
  
• 同理可以算出 ，然后用梯度下降法更新参数即可，下式中 表示学习率：
  
• 综上所述，高阶张量的通用公式如下，损失函数为：
  
• 同样令 偏导为:
  
• 参数更新公式为：
  
• 采用梯度下降法的思想也可以参考本篇文章：​​浅谈张量分解（五）：稀疏张量的CP分解​​

2.5. 梯度下降法Python代码

• 实现方式1： 逐个元素进行更新

# !/usr/bin/env python
# -*- encoding: utf-8 -*-
"""
@author : Kaifang Zhang
@time   : 2021/6/28
@Contact: [e-mail]
"""

import os
import numpy as np
import scipy.linalg as lin
from scipy.linalg import toeplitz

np.random.seed(2021)


def col_normalize(X):
    return X / np.sqrt(np.sum(X ** 2, axis=0))


def tensor2matrix(X, mode):
    """
    Return mu-mode matricization from a given tensor
    """
    num_dim = len(X.shape)
    n = X.shape[num_dim - mode]
    X = np.moveaxis(X, num_dim - mode, -1)
    return np.reshape(X, (-1, n)).T


def matrix2tensor(X1, out_shape):
    """
    Input: 1-mode matricization
    Output: tensor with size like (n3, n2, n1)
    """
    return np.reshape(X1.T, out_shape)


def als_solver_vector(X, r, nmax=1000, err_tol=1e-3):
    n3, n2, n1 = X.shape
    A = np.random.normal(0, 1, (n1, r))  # shape: (22, 10)
    B = np.random.normal(0, 1, (n2, r))  # shape: ( 7, 10)
    C = np.random.normal(0, 1, (n3, r))  # shape: ( 4, 10)

    X1 = tensor2matrix(X, 1)  # shape: (22, 28)
    X2 = tensor2matrix(X, 2)  # shape: (7, 88)
    X3 = tensor2matrix(X, 3)  # shape: (4, 154)

    X_norm = lin.norm(X1, 'fro')  # Frobenius范数
    err = np.inf  # 正无穷大的浮点表示

    # 构造两个系数。
    row1, col1 = np.zeros(B.shape[0]), np.zeros(B.shape[0])  # 生成两个toeplitz矩阵
    row1[0], row1[1], col1[0] = 1, -1, 1
    toe_D = toeplitz(col1, row1)  # shape:(7, 7)
    row2, col2 = np.zeros(C.shape[0]), np.zeros(C.shape[0])
    row2[0], row2[2], col2[0] = 1, -1, 1
    toe_T = toeplitz(col2, row2)  # shape:(4, 4)

    B = col_normalize(B)
    count, alpha, lambda_1, penalty_2, penalty_3 = 0, 2e-4, 4e-3, 4e-3, 1e-3  # 迭代次数以及惩罚系数
    while count < nmax:
        ######################################## optimize A: (100, 3)
        # C = col_normalize(C)
        D = lin.khatri_rao(C, B).T  # shape: (10, 28)

        # H, Q = toe_D.dot(B), toe_T.dot(C)  # H: (7, 10), Q: (4, 10)
        # M, S = lin.khatri_rao(C, H), lin.khatri_rao(B, Q)  # M: (28, 10); S: (28, 10)

        for u in range(X1.shape[0]):
            for i in range(X1.shape[1]):
                if X1[u][i] > 0:
                    eui = np.dot(A[u, :], D[:, i]) - X1[u][i]
                    # 带入公式，按照梯度下降算法更新当前的Au
                    for k in range(A.shape[1]):
                        A[u][k] = A[u][k] - 2 * alpha * (eui * D[k][i] + lambda_1 * A[u][k])

        D = lin.khatri_rao(C, A).T  # shape: ()
        for u in range(X2.shape[0]):
            for i in range(X2.shape[1]):
                if X2[u][i] > 0:
                    eui = np.dot(B[u, :], D[:, i]) - X2[u][i]
                    # 带入公式，按照梯度下降算法更新当前的Au
                    for k in range(B.shape[1]):
                        B[u][k] = B[u][k] - 2 * alpha * (eui * D[k][i] + lambda_1 * B[u][k])

        D = lin.khatri_rao(B, A).T  # shape: ()
        for u in range(X3.shape[0]):
            for i in range(X3.shape[1]):
                if X3[u][i] > 0:
                    eui = np.dot(C[u, :], D[:, i]) - X3[u][i]
                    # 带入公式，按照梯度下降算法更新当前的Au
                    for k in range(C.shape[1]):
                        C[u][k] = C[u][k] - 2 * alpha * (eui * D[k][i] + lambda_1 * C[u][k])

        count = count + 1
        print("当前迭代次数：", count)
        # Reconstruct tensor
        X_hat1 = A.dot(lin.khatri_rao(C, B).T)
        X_hat = matrix2tensor(X_hat1, X.shape)
    return X_hat


if __name__ == '__main__':
    # 1、加载数据集合
    tensor = np.array([[[1.1, 1.4, 1.6, 1.8, 0, 0, 4.3, 2.5],
                        [0.8, 3.3, 0, 0, 0, 0, 6.3, 5.9],
                        [2.1, 3.4, 0, 0, 1.1, 1.2, 4.3, 3.5],
                        [0, 0, 1.8, 2.3, 0, 0, 0.3, 0.6]],

                       [[11.1, 11.4, 0, 0, 1.7, 2.8, 9.3, 4.5],
                        [1.8, 4.3, 0, 0, 0.5, 3.5, 1.3, 9.9],
                        [2.1, 3.4, 0, 0, 1.1, 1.2, 2.3, 4.5],
                        [0, 0, 1.8, 2.3, 0, 0, 2.3, 7.6]],

                       [[11.1, 11.4, 0, 0, 1.7, 2.8, 9.3, 4.5],
                        [1.8, 4.3, 0, 0, 0.5, 3.5, 1.3, 9.9],
                        [2.1, 3.4, 0, 0, 1.1, 1.2, 2.3, 4.5],
                        [0, 0, 1.8, 2.3, 0, 0, 2.3, 7.6]]])
    # result = als_solver(tensor, 10)
    result = als_solver_vector(tensor, 10)
    tensor2matrix = np.reshape(tensor, (-1, tensor.shape[2]))
    result2matrix = np.reshape(result, (-1, result.shape[2]))
    print(result.shape)

• 实现方式2： 矩阵形式进行更新

import time
import numpy as np
from scipy.linalg import toeplitz, norm
import matplotlib.pyplot as plt


def kr_prod(a, b):
    """Return khatri-rao product of matrix a and b"""
    return np.einsum('ir, jr -> ijr', a, b).reshape(a.shape[0] * b.shape[0], -1)


def ten2mat(tensor, mode):
    """Return mu-mode matricization from a given tensor"""
    return np.reshape(np.moveaxis(tensor, mode, 0), (tensor.shape[mode], -1), order='F')


def mat2ten(mat1, mat2, mat3):
    """Return tensor based on matrix"""
    return np.einsum('ir, jr, tr -> ijt', mat1, mat2, mat3)


def alter_optimization(X, r, nmax=5000000, lambda_1=1e-4, lambda_d=0, lambda_t=0, alpha=1e-4):
    n1, n2, n3 = X.shape
    A = 0.1 * np.random.rand(n1, r)  # shape: (4, 10)
    B = 0.1 * np.random.rand(n2, r)  # shape: (7, 10)
    C = 0.1 * np.random.rand(n3, r)  # shape: (22, 10)
    ## 构造非0元素对应的张量
    pos = np.where(X != 0)  # where内只有一个参数时，那个参数表示条件，当条件成立时，where返回的是每个符合condition条件元素的坐标,返回的是以元组的形式
    bin_X = np.zeros((n1, n2, n3))
    bin_X[pos] = 1
    for iters in range(nmax):
        ######################################## optimize A
        var1 = kr_prod(C, B)
        xi = X - mat2ten(A, B, C)
        grad = lambda_1 * A - np.dot(ten2mat(bin_X * xi, 0), var1)
        A = A - alpha * grad / np.linalg.norm(grad)
        ######################################## optimize B
        var1 = kr_prod(C, A)
        xi = X - mat2ten(A, B, C)
        grad = lambda_1 * B - np.dot(ten2mat(bin_X * xi, 1), var1)
        B = B - alpha * grad / np.linalg.norm(grad)
        ######################################## optimize C
        var1 = kr_prod(B, A)
        xi = X - mat2ten(A, B, C)
        grad = lambda_1 * C - np.dot(ten2mat(bin_X * xi, 2), var1)
        C = C - alpha * grad / np.linalg.norm(grad)
        ######################################## Reconstruct tensor
        X_hat = mat2ten(A, B, C)
        loss = np.sum(np.square(X[pos] - X_hat[pos])) / X[pos].shape[0]
        if (iters + 1) % 200 == 0:
            print('迭代次数:', iters, '代价函数:', loss)
        # if (iters + 1) % 5000 == 0:
        #     alpha = alpha * 0.5
    return X_hat

三. 拓展应用-数据修补

3.1. 矩阵分解(Matrix factorization)

• 在​​MF​​​中，观测到的 维时间序列数据被组织在矩阵 中，矩阵 由维度特性矩阵 的与时间特性矩阵 的组合进行​​​低秩逼近​​​，从而修补缺失数据：
• 使用​​最小二乘法、梯度下降​​​等方法求解下述最小化问题，从而对矩阵 与矩阵 进行逼近：
  
• ，， 是原矩阵 中非零元所处位置的集合； 为残差矩阵​​​F范数​​​的平方，用来描述 矩阵 与原矩阵 的差异； 与 分别是 与 的正则项，用来防止过拟合。
  

• 以及 的过程。更新 需要最小化：
  
• 表示两个向量的内积。使用最小二乘法求解上述规划问题，推导 迭代式 ，首先将上式对 求导：
  
• 需要注意的是： 表示 。令该导数为零，即
  
• 又由：
  
• 变为：
  
• 的更新迭代式为：
  
• 可以参考这个Github：​​Transportation data online Prediction And Imputation(TransPAI).
  ​
• 下面代码部分的一些测试，为了方便理解：

# !/usr/bin/env python
# -*- encoding: utf-8 -*-
"""=====================================
@author : kaifang zhang
@time   : 2021/7/16 3:14 下午
@contact: [e-mail]
====================================="""
import numpy as np
from numpy.linalg import inv as inv


def kr_prod(a, b):
    return np.einsum('ir, jr -> ijr', a, b).reshape(a.shape[0] * b.shape[0], -1)


def CP_ALS(sparse_mat, m, n, k=5, lam=0.0001, nmax=10000):
    W = np.random.rand(m, k)
    X = np.random.rand(n, k)
    binary_mat = np.zeros((m, n))
    position = np.where((sparse_mat != 0))
    binary_mat[position] = 1
    for iter in range(nmax):
        var1 = X.T
        var2 = kr_prod(var1, var1)
        var3 = np.matmul(var2, binary_mat.T)
        var4 = np.matmul(var1, sparse_mat.T)
        for i in range(m):
            W[i, :] = np.matmul(inv((var3[:, i].reshape([k, k])) + lam * np.eye(k)), var4[:, i])

        var1 = W.T
        var2 = kr_prod(var1, var1)
        var3 = np.matmul(var2, binary_mat)
        var4 = np.matmul(var1, sparse_mat)
        for i in range(n):
            X[i, :] = np.matmul(inv((var3[:, i].reshape([k, k])) + lam * np.eye(k)), var4[:, i])
    return np.dot(W, X.T)


if __name__ == '__main__':
    """ R(n, T)
    评分矩阵R=PQ
    P: 初始化用户特征矩阵n * k
    Q：初始化物品特征矩阵k * T
    """
    R = np.array([[4, 0, 2, 0, 1],
                  [0, 0, 2, 3, 1],
                  [4, 1, 2, 0, 1],
                  [4, 1, 2, 5, 1],
                  [3, 0, 5, 0, 2],
                  [1, 0, 3, 0, 4]])
    mat = CP_ALS(R, R.shape[0], R.shape[1])
    print(mat)

/Users/devinzhang/miniconda3/envs/env1/bin/python /Users/devinzhang/codes/datasets/test_one.py
[[3.99995849 0.99998186 2.0000103  2.50675217 0.99998922]
 [2.49918409 0.65339661 1.99999078 2.99999346 1.00000338]
 [3.99996808 0.99998414 2.00001203 2.50675927 0.99998894]
 [4.         1.         2.00000498 4.99996453 0.99999818]
 [3.00001217 0.88680175 4.99993989 2.2969244  2.00003505]
 [0.9999976  0.32825754 3.00000097 0.88609748 3.99992235]]

Process finished with exit code 0

• 对比下面的梯度下降法实现：

#!/usr/bin/env python
# -*- coding: UTF-8 -*-
'''=====================================
@Author :Kaifang Zhang
@Time   :2021/7/5 1:31
@Contact: [e-mail]
========================================'''
import numpy as np


def LFM_grad_desc(R, K, max_iter, alpha=1e-4, lamda=1e-4):
    # 基本维度参数定义
    M = len(R)
    N = len(R[0])

    # P、Q初始值，随机生成
    P = np.random.rand(M, K)
    Q = np.random.rand(N, K)
    Q = Q.T
    # 开始迭代
    for step in range(max_iter):
        # 对所有的用户u、物品i做遍历，对应的特征向量Pu，Qi梯度下降
        for u in range(M):
            for i in range(N):
                # 对于每一个大于0的评分，求出预测的评分误差
                if R[u][i] > 0:
                    eui = np.dot(P[u, :], Q[:, i]) - R[u][i]

                    # 带入公式，按照梯度下降算法更新当前的Pu与Qi
                    for k in range(K):
                        P[u][k] = P[u][k] - alpha * (2 * eui * Q[k][i] + 2 * lamda * P[u][k])
                        Q[k][i] = Q[k][i] - alpha * (2 * eui * P[u][k] + 2 * lamda * Q[k][i])

        # u、i遍历完成，所有的特征向量更新完成，可以得到P、Q，可以计算预测评分矩阵
        predR = np.dot(P, Q)

        # 计算当前损失函数
        cost = 0
        for u in range(M):
            for i in range(N):
                if R[u][i] > 0:
                    cost += (np.dot(P[u, :], Q[:, i]) - R[u][i]) ** 2
                    # 加上正则化项
                    for k in range(K):
                        cost += lamda * (P[u][k] ** 2 + Q[k][i] ** 2)
        if step % 100 == 0:
            print("迭代次数：", step, "损失函数：", cost)
        if cost < 0.001:
            break

    return P, Q.T, cost


if __name__ == '__main__':
    '''
    @输入参数
    R:M*N的评分矩阵
    K:隐特征向量维度
    max_iter:最大迭代次数
    alpha:步长
    lamda:正则化系数
    @输出
    分解之后的P、Q
    P:初始化用户特征矩阵M*k
    Q：初始化物品特征矩阵N*K
    '''
    # 评分矩阵R
    R = np.array([[4, 0, 2, 0, 1],
                  [0, 0, 2, 3, 1],
                  [4, 1, 2, 0, 1],
                  [4, 1, 2, 5, 1],
                  [3, 0, 5, 0, 2],
                  [1, 0, 3, 0, 4]])

    # 给定超参数
    K = 5
    max_iter = 100000
    alpha = 1e-4
    lamda = 1e-3
    P, Q, cost = LFM_grad_desc(R, K, max_iter, alpha, lamda)
    predR = P.dot(Q.T)
    # 预测矩阵
    print(predR)

迭代次数： 99900 损失函数： 0.17000001031232792
[[3.99772084 1.04193669 1.99994564 3.82581908 0.99952361]
 [2.88106432 0.93827679 1.99879708 2.99945399 1.00023044]
 [3.99709708 0.99985647 1.99995041 3.57072602 0.99946976]
 [3.99849106 0.99963024 2.00082615 4.99726539 1.00047217]
 [3.00003145 2.04377479 4.99574099 4.57912227 2.0004716 ]
 [1.00055115 1.93433    2.99934762 3.40790853 3.99565708]]

Process finished with exit code 0

3.2. 张量分解(Tensor factorization)

• 事实上，​​CP分解​​​公式是​​Tucker分解​​的特例。在数学上，对于给定三阶张量 的任何 项，​​​CP分解​​​的形式可以写成：
  
• 其中核心张量 的超对角线项为1。换句话说, 对于任意的 可以知 ，否则 。
• 将CP分解作为一个机器学习问题，我们可以通过最小化因子矩阵上的损失函数来执行学习任务，即:
  
• 其中 表示一组 3 元组索引。在这个优化问题中，三个因子矩阵之间的乘法(作为参数)使这个问题变得困难。因此，我们将 ALS算法应用于CP分解。
• 特别地，固定矩阵 ，针对因子矩阵 中的每一行 ，对应的优化问题为(这里可以参考公式：
  
• 需要注意的是：， 以及 分别表示矩阵 , 以及 位置上 所有非零元素的位置索引构成的集合 。
• 具体推导为： 令损失函数对 求偏导数，令偏导数为0。最终可以得到这个优化的最小二乘是：
  
• 同理，可以得到 and 的最小二乘为：
  
• 定理 2：假设矩阵 为 。 如果张量 有这样的CP模型：
  
• 可以得到：
  
• 补充： ​​np.einsum('ir, jr -> ijr', a, b)​​

In [1]: import numpy as np
In [2]: a = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
In [3]: b = np.array([[1, 2, 3], [1, 2, 3], [1, 2, 3]])
In [4]: np.einsum('ir, jr -> ijr', a, b)
Out[4]: 
array([[[ 1,  4,  9],
        [ 1,  4,  9],
        [ 1,  4,  9]],

       [[ 4, 10, 18],
        [ 4, 10, 18],
        [ 4, 10, 18]],

       [[ 7, 16, 27],
        [ 7, 16, 27],
        [ 7, 16, 27]]])

• 补充：向量的外积。这里参考的是作者：​​矩阵乘法核心思想(5)：内积与外积​​
• 矩阵乘法也可看成是向量外积之和。外积对于数据科学非常重要，因为它能帮助我们发现矩阵中最重要的东西。
• ，如：
  
• 的列空间是 维的，向量方向都与 相同； 行空间也是 维，向量方向都与 相同。所有非零 都是秩 矩阵，它们是构建任何矩阵完美的基础砖块。上图体现了：矩阵乘法可看成为 个秩

In [1]: import numpy as np
In [3]: a = np.array([[2, 1, 1], [2, 2, 1], [1, 2, 3]])
In [4]: b = np.array([[3, 4, 5], [1, 2, 1], [1, 2, 2]])
In [5]: np.einsum('ir, rj -> irj', a, b)
Out[5]: 
array([[[ 6,  8, 10],
        [ 1,  2,  1],
        [ 1,  2,  2]],
       [[ 6,  8, 10],
        [ 2,  4,  2],
        [ 1,  2,  2]],
       [[ 3,  4,  5],
        [ 2,  4,  2],
        [ 3,  6,  6]]])

# !/usr/bin/env python
# -*- encoding: utf-8 -*-
"""=====================================
@author : kaifang zhang
@time   : 2021/7/13 3:50 下午
@contact: [e-mail]
====================================="""
import numpy as np
from scipy.linalg import khatri_rao


def ten2mat(tensor, mode):
    """Return mu-mode matricization from a given tensor"""
    return np.reshape(np.moveaxis(tensor, mode, 0), (tensor.shape[mode], -1), order='F')


def mat2ten(mat1, mat2, mat3):
    """Return tensor based on matrix"""
    return np.einsum('ir, jr, tr -> ijt', mat1, mat2, mat3)


def alter_optimization(X, r, nmax=50000):
    n1, n2, n3 = X.shape
    A = np.random.normal(0, 1, (n1, r))  # shape: (4, 10)
    B = np.random.normal(0, 1, (n2, r))  # shape: (7, 10)
    C = np.random.normal(0, 1, (n3, r))  # shape: (22, 10)

    pos = np.where(X != 0)  # where内只有一个参数时，那个参数表示条件，当条件成立时，where返回的是每个符合condition条件元素的坐标,返回的是以元组的形式
    bin_X = np.zeros((n1, n2, n3))
    bin_X[pos] = 1
    X_hat = np.zeros((n1, n2, n3))

    for iters in range(nmax):
        ######################################## optimize A
        var1 = khatri_rao(C, B).T
        var2 = khatri_rao(var1, var1)
        var3 = np.matmul(var2, ten2mat(bin_X, 0).T).reshape([r, r, n1])
        var4 = np.matmul(var1, ten2mat(X, 0).T)
        for i in range(n1):
            var_Lambda = var3[:, :, i]
            inv_var_Lambda = np.linalg.inv((var_Lambda + var_Lambda.T) / 2)
            A[i, :] = np.matmul(inv_var_Lambda, var4[:, i])
        ######################################## optimize B
        var1 = khatri_rao(C, A).T
        var2 = khatri_rao(var1, var1)
        var3 = np.matmul(var2, ten2mat(bin_X, 1).T).reshape([r, r, n2])
        var4 = np.matmul(var1, ten2mat(X, 1).T)
        for j in range(n2):
            var_Lambda = var3[:, :, j]
            inv_var_Lambda = np.linalg.inv((var_Lambda + var_Lambda.T) / 2)
            B[j, :] = np.matmul(inv_var_Lambda, var4[:, j])
        ######################################## optimize C
        var1 = khatri_rao(B, A).T
        var2 = khatri_rao(var1, var1)
        var3 = np.matmul(var2, ten2mat(bin_X, 2).T).reshape([r, r, n3])
        var4 = np.matmul(var1, ten2mat(X, 2).T)
        for t in range(n3):
            var_Lambda = var3[:, :, t]
            inv_var_Lambda = np.linalg.inv((var_Lambda + var_Lambda.T) / 2)
            C[t, :] = np.matmul(inv_var_Lambda, var4[:, t])
        ######################################## Reconstruct tensor
        X_hat = mat2ten(A, B, C)
        loss = np.sum(np.square(X[pos] - X_hat[pos])) / X[pos].shape[0]

        if (iters + 1) % 100 == 0:
            print('迭代次数:', iters, '代价函数:', loss)

    return X_hat


if __name__ == '__main__':
    tensor = np.load('./datasets/tensor.npy')
    X_hat = alter_optimization(tensor, 10)
    print(tensor.shape)

四. 参考文献

文章参考以下文献，这里表示感谢！

• 浅谈张量分解（一）：如何简单地进行张量分解？
• Low-Rank Tensor Networks for Dimensionality Reduction and Large-Scale Optimization Problems: Perspectives and Challenges PART 1
• ​​隐语义模型（Latent Factor Model, LFM）原理以及代码实现​​
• Tensor Learning (张量学习)
• 重要参考：tensor-decomposition-in-python
• 重要参考：https://github.com/YuxuanLongBeyond/Low_rank_tensor_approx
• GitHub项目tensor-learning新增内容 (1)：基于交替最小二乘法的张量分解技术
• NumPy中的维度(dimension)、轴(axis)、秩(rank)的含义
• ​​np.dot()、np.multiply()、np.matmul()方法以及*和@运算符的用法总结​​
• ​​张量CP分解原理及Python实现​​
• ​​Python中的张量分解！​​
• 重要参考： GitHub项目tensor-learning新增内容 (1)：基于交替最小二乘法的张量分解技术
• 矩阵乘法核心思想(5)：内积与外积
• ​​向量、矩阵、张量基础知识​​