二、张量 CP 分解

1. 张量的几个概念

---

1. 张量的内积

          相同大小的两个张量 $\mathcal{X},\mathcal{Y}\in {\mathbb{R}}^{I_1\times I_2\times \cdots \times I_N}$ ，其内积为对应位置的元素相乘后，将所有位置的乘积累加，可以用公式表示为：
   $$⟨\mathcal{X},\mathcal{Y}⟩=\sum _{i_1=1}^{I_1}\sum _{i_2=1}^{I_2}\cdots \sum _{i_N=1}^{I_N}{x}_{i_1{i}_2\dots i_N}{y}_{i_1{i}_2\dots i_N}$$

2. 张量的范数

          张量的范数可以类比矢量的模，它表示为张量自身的内积再开平方，对张量 $\mathcal{X}\in {\mathbb{R}}^{I_1\times I_2\times \cdots \times I_N}$ 其范数为：
   $$||\mathcal{X}||=\sqrt{\sum _{i_1=1}^{I_1}\sum _{i_2=1}^{I_2}\cdots \sum _{i_N=1}^{I_N}{x}_{i_1{i}_2\dots i_N}^2}$$

3. 矢量的外积

          两个矢量的外积是一个矩阵，三个矢量的外积是一个三阶矩阵。计算三个矢量的外积，我们可以用前两个矢量外积的矩阵，分别乘以第三个矢量的每个元素，得到第若干个矩阵按照正面切片拼接为新的张量。

          对矢量 $\vec{a}=(1,2{)}^T，\vec{b}=(3,4{)}^T，\vec{c}=(5,6,7{)}^T$ 其外积记为
   $$\mathcal{X}=\vec{a}\circ \vec{b}\circ \vec{c}$$

   先来看 $\vec{a}\circ \vec{b}$，可以得到
   $$\vec{a}\circ \vec{b}=\left[\begin{array}{c}1\\ 2\end{array}\right]\left[\begin{array}{cc}3& 4\end{array}\right]=\left[\begin{array}{cc}3& 4\\ 6& 8\end{array}\right]$$

   我们用 $\vec{c}$ 的第一个元素和上面的矩阵相乘，得到第一个切片，即
   $${\mathcal{X}}_{::1}=5\left[\begin{array}{cc}3& 4\\ 6& 8\end{array}\right]=\left[\begin{array}{cc}15& 20\\ 30& 40\end{array}\right]$$

   $${\mathcal{X}}_{::2}=6\left[\begin{array}{cc}3& 4\\ 6& 8\end{array}\right]=\left[\begin{array}{cc}18& 24\\ 36& 48\end{array}\right]$$

   $${\mathcal{X}}_{::3}=7\left[\begin{array}{cc}3& 4\\ 6& 8\end{array}\right]=\left[\begin{array}{cc}21& 28\\ 42& 56\end{array}\right]$$

4. 秩一张量

          对于一个 $N$ 阶张量 $\mathcal{X}\in {\mathbb{R}}^{I_1\times I_2\times \cdots \times I_N}$ ，如果可以被写成 $N$ 个向量的外积，则这个张量的秩为 1 ，这个张量称为秩一张量，记为：
   $$\mathcal{X}=a^{(1)}\circ a^{(2)}\circ \cdots \circ a^{(N)}$$

          每个张量的元素都可以写成这些向量对应位置元素之积：
   $$x_{i_1{i}_2\dots i_N}=a_{i_1}^{(1)}\circ a_{i_2}^{(2)}\circ \cdots \circ a_{i_N}^{(N)}$$

5. 张量的秩

          张量 $\mathcal{X}$ 的秩定义为用秩一张量之和来精确表示 $\mathcal{X}$ 所需要的秩一张量的最少个数，记为 $rank(\mathcal{X})$ 。

6. 对称和对角

   立方张量：各 mode 的长度相等的张量。

   对称张量：如果一个立方张量的各元素再下标的任意排列下是常数，称该张量是对称的。

   超对称张量：当立方张量中的任何一个元素的索引被置换后元素值不变时，称这个张量为超对称的。如，对于一个三阶张量，如果各元素满足以下等式，则被称之为超对称。
   $$x_{ijk}=x_{ikj}=x_{jik}=x_{jki}=x_{kij}=x_{kji}i,j,k=1,\dots ,I$$

   对角张量：如果一个张量 $\mathcal{X}\in {\mathbb{R}}^{I_1\times I_2\times \cdots \times I_N}$ 的任何元素只有在 $i_1=i_i=\cdots =i_N$ 的时候不为零，被称为对角张量。

2. 矩阵的分解

---

矩阵分解是指将一个矩阵分解成两个或者多个矩阵的乘积。这里我们先介绍几种矩阵分解的方式。

2.1 非负矩阵分解

定义：对于矩阵 $X\in R^{D\times N}$ ，非负矩阵分解的目的是找到两个非负 $U^{D\times R},V^{R\times N}$ ，使得它们的乘积近似等于 $X$ 即：
$$X\approx U{V}^T$$

可以看出，非负矩阵的目的就是找到 $U,V$ 使得 $U{V}^T$ 最大可能的和 $X$ 相等。这里我们定义一种方式来衡量 $U{V}^T$ 与 $X$ 接近的程度，称为代价函数。如果两个非负矩阵 $A$ 和 $B$ ，我们定义它们之间的代价函数为
$$\Vert A-B{\Vert }^2=\sum _{ij}(A_{ij}-B_{ij}{)}^2$$

那么非负矩阵的分解的目的就可以演变为求
m i n { ∥ X − U V T ∥ } min\{\|X-UV^T\|\} min{∥X−UVT∥}

即使以上等式最接近 0 的 $U,V$的最优解。

上述关于 $U,V$ 的代价函数是非凸的，基本的做法是交替优化 U 和 V 从而得到 一个局部最优解。

2.2 本征值分解（EVD）

       给定一个 $m\times m$ 的矩阵 $A$ 设 $m$ 维归一向量 $v\text{v}v$ 与标量 $\lambda$ ，当其满足 $Av\text{v}v=\lambda v\text{v}v$ 时，称 $v\text{v}v$ 与 $\lambda$ 分别为 $A$ 的本征值和本征向量，我们可以通过线性代数中求特征值和特征向量的方式求出所有的本征值和本征向量。

       本征值分解就是对于一个 $m\times m$ 的矩阵 $A$ （即 $A=A^T$） ，可以将其分解为如下的形式
$$A=Q\Lambda Q^T$$

其中 $Q$ 是 $A$ 的本征向量组成的正交矩阵，$\Lambda$ 是本征值为对角线元素构成的对角矩阵，也称为本征谱。

       在特征值分解中有一个 最大本征值问题 ：假设矩阵本征值为实数，求解给定矩阵的最大本征值及其对应的本征向量。该问题可以转化为如下问题，给定矩阵 $A$ ，求解归一化向量 $v\text{v}v$ ，使得函数
$$f=|v^{T}Av|$$

的值极大化。

从线性代数可以证明，该问题的解为 $A$ 的绝对值最大的本征向量，而且对应的 $f$ 的值为该本征向量对应的本征值，也是最大的一个本征值。

最大本征值问题有一种求解方法为幂级数法，通过证明我们可以得到如下等式：
$$\underset{K\to \infty }{\lim }{M}^K={\Gamma }_0^K{u}_{(0)}{u}_{(0)}^T$$

其中，${\Gamma }_0$ 与 $u_{(0)}$ 是绝对值最大的本征值与对应的本征向量。

用 python 求解本征值和本征向量

我们可以用 numpy 库中的函数来求解本征值和本征向量，代码如下：

import numpy as np
dim = 4          #设置矩阵的维数
M = np.random.randn(dim, dim)      #随机化一个 dim × dim 大小的矩阵
M = M + M.T  # 对称化以保证本征分解存在; .T 是对矩阵进行转置
print(M)    

lm, u = np.linalg.eig(M)      #该函数计算矩阵 M 的本征值和本征向量，并将本征值组成列表返回给 lm ，将每个本征向量作为 u 的一列，组成矩阵
print('\n Eigenvalues:')
print(lm)
print('\n Eigenvectors:')
print(u)

输出结果如下

2.3 奇异值分解（SVD）

       特征值分解对矩阵有着较高的要求，它需要被分解的矩阵 $A$ 为实对称矩阵，但是现实中，我们所遇到的问题一般不是实对称矩阵，奇异值分解就是将一个一般性的 $m\times n$ 的矩阵 $A$ ，分解为
$$A=U\Sigma V^T$$

的形式，其中 $U,V$ 均为单位正交矩阵，即有 $U{U}^T=I$ 和 $V{V}^T=I$ ，$U$ 称为 左奇异矩阵，$V$ 称为 右奇异矩阵，$\Sigma$ 是将 $A$ 的特征值降序排列在主对角线上的对角矩阵，对角线上的值称为 奇异值 。其中矩阵 $U\in R^{m\times m},\Sigma \in R^{m\times n},V\in R^{n\times n}$ 。

我们如何对奇异值和奇异矩阵求解呢？通过上面的式子，我们可以得到如下的等式：
$$A{A}^T=U\Sigma V^{T}V{\Sigma }^T{U}^T=U\Sigma {\Sigma }^T{U}^T$$

$$A^{T}A=V{\Sigma }^T{U}^{T}U\Sigma V^T=V{\Sigma }^T\Sigma V^T$$

因为 $A{A}^T$ 和 $A{A}^T$ 是对称矩阵，则通过对上面两个等式求它们的特征值和特征向量，就可以求出 $U,V,\Sigma$ 。

奇异值分解的性质：

1. 奇异值分解可以降维
          如果 $A$ 表示 $n$ 个 $m$ 维向量，可以通过奇异值分解表示成 $m+n$ 个 $r$ 维向量。若A的秩远远小于 $m$ 和 $n$ ，则通过奇异值分解可以降低 $A$ 的维数。可以计算出，当 $r<\frac{mn}{m+n+1}$ 时，可以达到降维的目的，同时可以降低计算机对存储器的要求。
2. 奇异值对矩阵的扰动不敏感
          特征值对矩阵的扰动是敏感的。在数学上可以证明，奇异值的变化不会超过相应矩阵的变化，即对任何的相同阶数的实矩阵 $A、B$ 按从大到小排列的奇异值 ${\alpha \text{α}\alpha }_i$ 和 ${\omega \text{ω}\omega }_i$ 有
   $$\sum |{\alpha \text{α}\alpha }_i-{\omega \text{ω}\omega }_i|\le \Vert A-B{\Vert }_2$$
3. 奇异值具有比例不变性。
          即 $\alpha A$ 的奇异值是 $A$ 的奇异值的 $|\alpha |$ 倍。
4. 奇异值具有旋转不变性。
          若 $P$ 是正交矩阵，$PA$ 的奇异值与 $A$ 的奇异值相同。

矩阵的低秩近似问题

首先我们来定义矩阵的秩，在奇异值分解中，非零奇异值的个数称为该矩阵的秩。那么我们便可以引出矩阵的低秩近似问题，该问题可以描述为：对于给定的 $m\times n$ 的矩阵 $A$ ，设其秩为 $R$ ，求解秩为 $R^{{}^{\prime }}$ 的矩阵 $M^{{}^{\prime }}$ ，其中 $R>R^{{}^{\prime }}>0$ ，使两矩阵之间的范数最小，即求尽可能小的 $\epsilon$
$$\epsilon =\Vert M-M^{{}^{\prime }}\Vert =\sqrt{\sum _{ij}(M_{ij}-M_{ij}^{{}^{\prime }}{)}^2}\sim \Vert {\Sigma }_{R^{{}^{\prime }}:R-1}\Vert$$

该 $\epsilon$ 又称为裁剪误差，上式后面的式子为 $\Sigma$ 的后面几个奇异值的范数，我们可以直接求出裁剪误差。

该低秩问题的最优解为：
$$M^{{}^{\prime }}=U[:,0:R^{{}^{\prime }}-1]\Sigma [0:R^{{}^{\prime }}-1,0:R^{{}^{\prime }}-1]V[:,0:R^{{}^{\prime }}-1]$$

即取 $U$ 的前 $R^{{}^{\prime }}$ 个列向量组成的矩阵，$\Sigma$ 的前 $R^{{}^{\prime }}$ 个特征值组成的对角阵，$V$ 的前 $R^{{}^{\prime }}$ 个列向量组成的矩阵的转置，它们的乘积得到的矩阵就是 $M^{{}^{\prime }}$ 。即下图所示：

用 python 实现奇异值分解

奇异值分解我们可以直接使用 np.linalg.svd 函数来实现，具体代码如下：

import numpy as np
M=M = np.random.randn(5, 5)
U, S, V = np.linalg.svd(M)  #返回三个矩阵使得 M=USV
print(U)
print(S)
print(V)

运行结果为：

通过奇异值的低秩近似可以用来压缩图片，一张图片可以视为一个矩阵 $M$，我们对该矩阵进行奇异值分解，用上面的低秩近似的方法求秩为 $R^{{}^{\prime }}$ （$R^{{}^{\prime }}$ 的大小可以自由选择，$R^{{}^{\prime }}$ 越大，图片越清晰，但是压缩率越小）的矩阵，也就是对求得的 $U,\Sigma ,V$ 取前 $R^{{}^{\prime }}$ 个向量和值的乘积构成新的图片，也就是压缩后的图片。下面的代码展示了这个应用：

import numpy as np
import cv2
import matplotlib.pyplot as plt
import scipy.sparse.linalg as la

img = cv2.imread('./Imgs/example2.jpg')  # 读取RGB图片
img = np.sum(img, axis=2) / 3      #将RGB图片的三张分量R,G,B矩阵相加后每个元素除以三，得到灰度图

def img_compress(img, k):        
    u, lm, v = la.svds(img, k=k)   #该函数可以直接求出取k个特征值和特征向量时的奇异值分解
    img1 = u.dot(np.diag(lm)).dot(v)   #将奇异值分解得到的矩阵相乘得到压缩后的图片
    return img1

img1 = img_compress(img, 20)     #取前 20 个特征值进行压缩后的图片
img2 = img_compress(img, 200)    #取前 200 个特征值进行压缩后的图片

plt.subplot(1,3,1)               #用 1 × 3 的排列方式在第一个位置显示图片
plt.imshow(img)
plt.subplot(1,3,2)               #用 1 × 3 的排列方式在第二个位置显示图片
plt.imshow(img1)
plt.subplot(1,3,3)               #用 1 × 3 的排列方式在第三个位置显示图片
plt.imshow(img2)

plt.show()                       #显示整张图

得到的结果图片如下，可以发现，随着选取特征值个数的增多，图片更加清晰，但是压缩率会更小。

3. CP 分解

---

       CP 分解是一种对张量进行拆分的方法，其核心思想是用有限个秩一张量的和来 近似 地表示该张量。

       如上图所示，如果要把一个三阶张量 $\mathcal{X}\in {\mathbb{R}}^{I\times J\times K}$ 进行分解，我们期望结果如下：
$$\begin{array}{rl}\mathcal{X}\approx & a_1\circ b_1\circ c_1+a_2\circ b_2\circ c_2+\cdots +a_R\circ b_R\circ c_R\\ =& \sum _{r=1}^R{a}_r\circ b_r\circ c_r\end{array}$$

如果我们记
$$\mathrm{A}=\left[\begin{array}{cccc}{a}_1& a_2& \dots & a_R\end{array}\right]$$

$$\mathrm{B}=\left[\begin{array}{cccc}{b}_1& b_2& \dots & b_R\end{array}\right]$$

$$\mathrm{C}=\left[\begin{array}{cccc}{c}_1& c_2& \dots & c_R\end{array}\right]$$

我们称 $A,B,C$ 为因子矩阵，则张量可以表示为
$$\mathcal{X}\approx [[\mathrm{A},\mathrm{B},\mathrm{C}]]\equiv \sum _{r=1}^R{\mathrm{a}}_r\circ {\mathrm{b}}_r\circ {\mathrm{c}}_r$$

这个公示就是张量的 CP 分解。

这里有个性质：$Rank(A)+Rank(B)+Rank(C)\ge 2R+2$ 。

       通常为了计算便利，我们假设 $A,B$ 和 $C$ 的列向量都是归一化的，我们引入一个权重向量 $\lambda \in {\mathbb{R}}^R$，来表示每个秩一矩阵所占权重，则分解可以记为如下形式：
$$\mathcal{X}\approx [[\lambda ;\mathrm{A},\mathrm{B},\mathrm{C}]]\equiv \sum _{r=1}^R{\lambda }_r{\mathrm{a}}_r\circ {\mathrm{b}}_r\circ {\mathrm{c}}_r$$

       三阶张量是应用当中最为广泛也往往是足够满足我们需求的张量。但是对于一个N阶的张量 $\mathcal{X}\in {\mathbb{R}}^{I_1\times I_2\times \cdots \times I_N}$ ，其 CP 分解可以被写为：
$$\mathcal{X}\approx [[\lambda ;{\mathrm{A}}^{(1)},{\mathrm{A}}^{(2)},\dots ,{\mathrm{A}}^{(N)}]]\equiv \sum _{r=1}^R{\lambda }_r{\mathrm{a}}_r^{(1)}\circ {\mathrm{a}}_r^{(1)}\circ \cdots \circ {\mathrm{a}}_r^{(N)}$$

其中 $\lambda \in {\mathbb{R}}^R$ 且 ${\mathrm{A}}^{(n)}\in {\mathbb{R}}^{I_n\times R}$ 。

分解的切片表示

       利用因子矩阵，一个三阶张量的 CP 分解可以被等价写作以下矩阵形式，其左侧都是张量的对应的 mode 的矩阵化：
$$\begin{gathered} {\mathrm{X}}_{(1)}\approx \mathrm{A}(\mathrm{C}\odot \mathrm{B}{)}^{\mathsf{T}} \\ {\mathrm{X}}_{(2)}\approx \mathrm{B}(\mathrm{C}\odot \mathrm{A}{)}^{\mathsf{T}} \\ {\mathrm{X}}_{(3)}\approx \mathrm{C}(\mathrm{B}\odot \mathrm{A}{)}^{\mathsf{T}} \end{gathered}$$

以上的模型还可以用张量的正面切片方式表示为：
$${\mathcal{X}}_{(k)}\approx \mathrm{A}{\mathrm{D}}^{(k)}{\mathrm{B}}^{\mathsf{T}}$$

其中 ${\mathrm{D}}^{(k)}\equiv \text{diag}(c_{k:}),k=1,...,K$ ，整个表示可以用下图表示：

对于高维的张量，其 mode - n 矩阵化可以写为：
$${\mathrm{X}}_{(n)}\approx {\mathrm{A}}^{(n)}\Lambda ({\mathrm{A}}^{(N)}\circ \cdots \circ {\mathrm{A}}^{(n+1)}\circ {\mathrm{A}}^{(n-1)}\cdot \cdots \cdot {\mathrm{A}}^{(1)}{)}^{\mathsf{T}}$$

其中 $\Lambda =diag(\lambda )$ 。

张量的秩分解

       上面我们已经介绍了张量的秩，其值为为表示张量所需秩一张量的最小数目，精确确定表示张量的每个秩一张量就是 秩分解 。CP 分解是对张量的近似表示，因此可以将秩分解视为 CP 分解的特例。

       虽然张量秩的定义和矩阵类似，但他们的性质之间存在很多不同。其中一个不同便是在 $\mathbb{R}$ 和 $\mathbb{C}$ 之下, 张量可以存在不同的秩。如，对于正面切片如下的张量：
$${\mathrm{X}}_1=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]，{\mathrm{X}}_2=\left[\begin{array}{cc}0& 1\\ -1& 0\end{array}\right]$$

这个张量在 $\mathbb{R}$ 下的秩为 3 ，但是在 $\mathbb{C}$ 下的秩为 2 。其 $\mathbb{R}$ 下的秩分解 $\mathcal{X}=[[\mathrm{A},\mathrm{B},\mathrm{C}]]$ 为：
$$\mathrm{A}=\left[\begin{array}{ccc}1& 0& 1\\ 0& 1& -1\end{array}\right],\mathrm{B}=\left[\begin{array}{ccc}1& 0& 1\\ 0& 1& 1\end{array}\right],\mathrm{C}=\left[\begin{array}{ccc}1& 1& 0\\ -1& 1& 1\end{array}\right]$$

在 $\mathbb{C}$ 之下的秩分解为：
$$\mathrm{A}=\frac{1}{\sqrt{2}}\left[\begin{array}{cc}1& 1\\ -i& i\end{array}\right],\mathrm{B}=\frac{1}{\sqrt{2}}\left[\begin{array}{cc}1& 1\\ i& -i\end{array}\right],\mathrm{C}=\frac{1}{\sqrt{2}}\left[\begin{array}{cc}1& 1\\ i& -i\end{array}\right]$$

注意：张量的秩没有一个简单直接的求法，直接求秩本身已被证明是 NP - hard 问题。

CP 分解的计算

       张量的秩没有直接的求法，但是如果想要进行 CP 分解，必须要知道张量的秩，那么我们怎么计算 CP 分解呢？从理想上说，如果数据是没有噪音的而且在确定了张量的秩后，其 CP 分解的计算方法是已知的，那么我们可以假设张量的秩依次为 $R=1,2,3,\cdots$ ，用该秩去计算 CP 分解，直到得到分解的拟合精度达到 100%。

       那么如果我们假设了张量的秩为 R ，那么计算方法是什么样的呢？现在在给定秩之后计算 CP 分解的方法有很多，其中的一种方法 ALS(alternating least squares)(交替最小方差法)，是一类比较有效的算法。

       在数学最优化问题中，拉格朗日乘数法是用来寻找多元函数的极值的方法。这种方法将一个有 $n$ 个变量与 $k$ 个约束条件的最优化问题转换为一个有 $n+k$ 个变量的方程组的极值问题，其变量不受任何约束。这种方法引入了一种新的标量未知数，即拉格朗日乘数：即形成方程的线性组合里每个矢量的系数。CP 分解的计算可以转换为 ALS 的一个子问题。

       我们用非负矩阵里面的代价函数来衡量 CP 分解得到结果的相似度，以三阶张量为例，我们令
$$\hat{\mathcal{X}}=[[\lambda ;\mathrm{A},\mathrm{B},\mathrm{C}]]\equiv \sum _{r=1}^R{\lambda }_r{\mathrm{a}}_r\circ {\mathrm{b}}_r\circ {\mathrm{c}}_r$$

那么 CP 分解的计算就转变为求
$$\min \Vert \mathcal{X}-\hat{\mathcal{X}}\Vert$$

求使上式尽可能小的 $\hat{\mathcal{X}}$ 就能得到最终的结果。

上面的问题可以作为 ALS 的一个子问题，因为分解可以表示为
$${\mathrm{X}}_{(1)}\approx \mathrm{A}(\mathrm{C}\odot \mathrm{B}{)}^{\mathsf{T}}$$

则问题可以转变为固定 $B,C$ ，求解
$$\min \Vert {\mathrm{X}}_{(1)}-A\text{diag}(\lambda )(C\odot B{)}^T\Vert$$

可以得到
$$A\text{diag}(\lambda )={\mathrm{X}}_{(1)}[(\mathrm{C}\odot \mathrm{B}{)}^{\mathsf{T}}{]}^H={\mathrm{X}}_{(1)}(\mathrm{C}\odot \mathrm{B})({\mathrm{C}}^{\mathsf{T}}\mathrm{C}\ast {\mathrm{B}}^{\mathsf{T}}\mathrm{B}{)}^H$$

再通过归一化分别求出 $A$ 和 $\lambda$ 。