（3）多线性代数基础

1.多线性代数中的张量单秩分解

1.1 多线性代数又称张量分解

思考：矩阵（二阶张量）存在本征值分解、奇异值分解等，对于更高阶的张量，是否存在类似的分解？
举例：考虑矩阵$M$（简要起见设$M$为实矩阵）的奇异值分解，记第$n$个左、右奇异向量分别为$u^{(n)}$和$v^{(n)}$，证明如下等式：
已知奇异值分解 $M=U\Lambda V^{†}$
——将$M$与第$n$个左奇异向量做矩阵乘（缩并），得到的刚好是第$n$个奇异值乘上第$n$个右奇异值向量：
$$\sum _a{u}_a^{(n)}{M}_{ab}={\Lambda }^{(n)}{v}_b^{(n)}$$
用图形表示为：

——将$M$与第$n$个右奇异向量做矩阵乘（缩并），得到的刚好是第$n$个奇异值乘上第$n$个左奇异值向量：

$$\sum _b{M}_{ab}{v}_b^{(n)}={\Lambda }^{(n)}{u}_a^{(n)}$$
用图形表示为：

其中，${\Lambda }^{(n)}$为第$n$个奇异值，且有${\Lambda }^{(n)}={\sum }_{ab}{u}_a^{(n)}{M}_{ab}{v}_b^{(n)}$
提示：利用变换矩阵的正交性

1.1.1 代码实例

import numpy as np

#验证第n组左右奇异向量与奇异值满足自洽方程

dim = 4
M = np.random.randn(dim,dim) #随机生成矩阵
U,S,V = np.linalg.svd(M) #进行svd分解，存储左右奇异矩阵和奇异值
n = 2 #检查第n个组左右奇异向量
u = U[:,n]
s = S[n]
v = V[n,:]
print(s*v) #第2个奇异值乘以第二个右奇异向量
print(u.dot(M))#第2个左奇异向量乘以矩阵

1.2 计算最大奇异值及奇异向量的迭代算法

类似于矩阵的幂级数算法
（a）随机初始化归一向量$u$和$v$；
（b）利用下面式子，计算$u$和$M$的收缩并归一化，更新$v$，归一化因子记为$\Lambda$；
$$\sum _a{u}_a{M}_{ab}=\Lambda v_b$$
（c）利用下面式子，计算$v$和$M$的收缩并归一化，更新$u$，归一化因子记为$\Lambda$；
$$\sum _b{M}_{ab}{v}_b=\Lambda u_a$$
（d）如果$u$和$v$（以及$\Lambda$）收敛，则返回$u、v和\Lambda$；否则，返回执行步骤（b）；
最后我们发现$u$和$v$刚好是$M$的最大的左右奇异向量，$\Lambda$刚好对于$M$的最大的奇异值
定义：
——仅最大奇异向量为上述迭代过程的稳定不动点
——当矩阵为是实对称时，上述算法可获得绝对值最大的本征值及对应的本征向量；

1.2.1 代码实例

——定义一个函数

#假设矩阵M为实矩阵
def svd0(mat,it_time = 100,tol = 1e-15):
	"""
	mat:输入矩阵（假设是实矩阵）
	it_time：最大迭代次数
	tol：错误阈值
	u：左奇异向量
	s：奇异值
	v：右奇异向量
	"""
	dim0,dim1 = mat.shape#dim0为行数，dim1为列数
	#随机初始化奇异向量
	u,v = np.random.randn(dim0,),np.random.randn(dim1,)
	#归一化初始向量
	u,v = u/np.linalg.norm(u),v/np.linalg.norm(v)
	s = 1
	
	for t in range(it_time):
		#更新v和s
		v1 = u.dot(mat)
		s1 = np.linalg.norm(v1)
		v1 /= s1 #归一化v1
		#更新u和s
		u1 = mat.dot(v1)
		s1 = np.linalg.norm(u1)
		u1 /= s1
		#计算收敛程度
		conv = np.linalg.norm(u - u1) / dim0 + np.linalg.norm(v - v1) / dim1
		u,s,v = u1,s1,v1
		#判断是否跳出循环
		if conv < tol:
		   break
	return u,s,v

验证上面算法的正确性
利用numpy库中的svd函数计算最大奇异值及对应的奇异向量

M = np.random.randn(3,5)
U,S,V = np.linalg.svd(M) #自带的svd中都是非升序排列的
u0,s0,v0 = U[:,0],S[0],V[0,:]
#利用上述算法计算最大奇异值及对应的奇异向量
u1，s1,v1 = svd0(M)
print(u0)
print(u1)
print(v0)
print(v1)
print(s0)
print(s1)

1.3 自然地将上述自洽方程组推广到高阶张量的单秩分解（rank-1分解；以三阶张量T为例）

假设：$$T=\Lambda u⨂v⨂w$$
满足如下自洽方程组：
$$\sum _{ab}{T}_{abc}{u}_a{v}_b=\Lambda w_c$$$$\sum _{ac}{T}_{abc}{u}_a{w}_c=\Lambda v_b$$$$\sum _{bc}{T}_{abc}{v}_b{w}_c=\Lambda u_a$$
注：第一个式子是$u_a,v_b$与$T_{abc}$进行收缩，以下两个同理
其中$u,v,w$为归一化向量，方程组成立时有$\Lambda ={\sum }_{abc}{T}_{abc}{u}_a{v}_b{w}_c$被称为rank-1分解系数
类似于奇异值分解的迭代算法，上述方程组可迭代求解。

1.4 Rank-1分解

Rank-1分解对应于如下优化问题：
给定张量$T$，求解张量$\tilde{T}=\Lambda {\prod }_{⨂n}{v}^{(n)}$，使得$T$和$\tilde{T}$之间的范数$f=|T-\tilde{T}|$极小。（注：张量的范数这里定义为张量元平方和之后开方）
可以证明：Rank-1分解得到的向量与系数构成的$\tilde{T}$,使范数$f$极小
思考：当$T$为二阶张量时，Rank-1分解与矩阵奇异值分解对应的最优问题之间的关系。

1.4.1 代码实例

import numpy as np
import copy

def rank1decomp(x,it_time = 100,tol = 1e-15)
	"""
	x:待分解的张量
	it_time:最大迭代步数
	tol:迭代终止的阈值
	vs:存储rank-1分解各个向量的list
	k:rank_1系数
	"""
	ndim = x.ndim #读取张量x的阶数
	dims = x.shape #读取张量x各个指标的维数

	#初始化vs中的各个向量并归一化
	vs = list() #vs用以存储rank-1分解得到的各个向量
	for n in range(ndim):
        _v = np.random.randn(dims[n])
        vs.append(_v / np.linalg.norm(_v))
    k = 1
	for t in range(it_time):
		vs0 = copy.deepcopy(vs) #暂存各个向量以计算收敛情况
		for _ in range(ndim):
			#收缩前（ndim-1）个向量，更新最后一个向量
			x1 = copy.deepcopy(x)
			for n in range(ndim - 1)：
				x1 = np.tensordot(x1,vs[n],[[0],[0]])
			#归一化得到的向量，并更新常数k
			k = np.linalg.norm(x1)
			x1 /= k
			#将最后一个向量放置到第0个位置
			vs.pop()
			vs.insert(0,x1)
			#将张量最后一个指标放置到第0位置
			x = x.transpose([ndim - 1] + list(range(ndim - 1))) #transpose（）函数的作用是调换数组的行列值的索引值，类似于求矩阵的转置
		#计算收敛情况
		conv = np.linalg.norm(np.hstack(vs0) - np.hstack(vs))#np.hstack（）在水平方向上平铺
		if conv < tol:
			break
	return vs,k
		

测试rank-1分解程序

tensor = np.random.randn(2,2,2,3,4)
vecs,coeff = rank1decomp(tensor)
print('The vectors by rank-1 decomposition = ')
for x in vecs:
	print(x)
print('\nThe rank-1 coefficient = ' + str(coeff))

注：以上代码是对于任意阶张量而言的

2.高阶奇异值分解

在rank-1分解中，我们将最大奇异值及奇异向量的方法推广到了高阶张量，下面我们考虑将完整的奇异值分解进行推广
——高阶奇异值分解（Higher-order singular value decomposition,简称HOSVD）,又称Tucker分解
定义如下（以三阶实张量T为例）：
$$T_{abc}=\sum _{ijk}{G}_{ijk}{U}_{ai}{V}_{bj}{W}_{ck}（从右向左叫Tucker积）$$
将一个高阶张量分解为一个同阶的张量和三个矩阵的缩并
Tucker分解的图形表示：

$$\Downarrow$$

其中，变换矩阵满足正交性$U{U}^T=V{V}^T=W{W}^T=I$，张量$G$被称为核张量（core tensor）,满足如下性质：
①定义$G$的键约化矩阵（bond reduced matrix）,以指标$i$为例，其键约化矩阵定义为：
$$J_{i{i}^{{}^{\prime }}}=\sum _{jk}{G}_{ijk}{G}_{i^{{}^{\prime }}jk}$$
（只露出一个指标，其它指标进行缩并）
核张量的任意$G$的任意键约化矩阵为非负定对角阵，且元素按非升序排列（$J_{00}\ge J_{11}\ge \dots \dots \ge 0$）
②高阶奇异值分解的HOSVD算法（以三阶实张量$T$为例）
（a）计算各个指标的键约化矩阵；
$$I_{a{a}^{{}^{\prime }}}=\sum _{bc}{T}_{abc}{T}_{a^{{}^{\prime }}bc}$$$$J_{b{b}^{{}^{\prime }}}=\sum _{ac}{T}_{abc}{T}_{a{b}^{{}^{\prime }}c}$$$$K_{c{c}^{{}^{\prime }}}=\sum _{ab}{T}_{abc}{T}_{ab{c}^{{}^{\prime }}}$$
（b）计算每个键约化矩阵的本征值分解：
键约化矩阵一定是实对称的
$$I=U\Omega U^T$$$$J=V\Pi V^T$$$$K=W{\rm Y}{W}^T$$
（c）计算核张量：$G_{ijk}={\sum }_{abc}{T}_{abc}{U}_{ai}{V}_{bj}{W}_{ck}$
最终得到高阶奇异值分解：$T_{abc}={\sum }_{ijk}{G}_{ijk}{U}_{ai}{V}_{bj}{W}_{ck}$
注：
各个本征谱中，非零值得个数被称为张量的Tucker秩
张量的Tucker低秩近似可通过HOSVD,HOOL等算法实现

2.1 代码实例

import numpy as np
import torch as tc
import copy

#求解hosvd的算法
def hosvd(x):
    """
    :param x: 待分解的张量
    :return G: 核张量
    :return U: 变换矩阵
    :return lm: 各个键约化矩阵的本征谱
    """
    ndim = x.ndim
    U = list()  # 用于存取各个变换矩阵
    lm = list()  # 用于存取各个键约化矩阵的本征谱
    x = tc.from_numpy(x)#将x 转换为torch类型
    for n in range(ndim):
        index = list(range(ndim))
        index.pop(n)
        _mat = tc.tensordot(x, x, [index, index])
        _lm, _U = tc.symeig(_mat, eigenvectors=True)#求对称矩阵的特征值和特征向量 如果eigenvectors为False，则仅计算特征值。如果为True，则同时计算特征值和特征向量。
        lm.append(_lm.numpy())#_lm.numpy()将_lm又转换为numpy()类型
        U.append(_U)
    # 计算核张量
    G = tucker_product(x, U)
    U1 = [u.numpy() for u in U]#转化类型
    return G, U1, lm


#从右到左的tucker积
def tucker_product(x, U, dim=1):
    """
    :param x: 张量
    :param U: 变换矩阵
    :param dim: 收缩各个矩阵的第几个指标
    :return G: 返回Tucker乘积的结果
    """
    ndim = x.ndim
    if type(x) is not tc.Tensor:
        x = tc.from_numpy(x)
        
    U1 = list()
    for n in range(len(U)):
        if type(U[n]) is not tc.Tensor:
            U1.append(tc.from_numpy(U[n]))
        else:
            U1.append(U[n])
    
    ind_x = ''
    for n in range(ndim):
        ind_x += chr(97 + n)
    ind_x1 = ''
    for n in range(ndim):
        ind_x1 += chr(97 + ndim + n)
    contract_eq = copy.deepcopy(ind_x)
    for n in range(ndim):
        if dim == 0:
            contract_eq += ',' + ind_x[n] + ind_x1[n]
        else:
            contract_eq += ',' + ind_x1[n] + ind_x[n]
    contract_eq += '->' + ind_x1
    # print(x.shape, U[0].shape, U[1].shape, U[2].shape)
    # print(type(contract_eq), contract_eq)
    G = tc.einsum(contract_eq, [x] + U1)#爱因斯坦求和公式
    G = G.numpy()
    return G

检查Tucker分解

tensor = np.random.randn(3, 4, 2)
Core, V, LM = hosvd(tensor)

print('检查Tucker分解等式是否成立：')
tensor1 = tucker_product(Core, V, dim=0)
error = np.linalg.norm(tensor - tensor1)
print('Tucker分解误差 = ' + str(error))

print('\n检查各个变换矩阵为正交阵：')
for n, v in enumerate(V):
    print('\n对于第' + str(n) +'个变换矩阵，VV.T = ')
    print(v.dot(v.T))

3.拓展：张量的CP分解与CP秩

回顾：单秩张量（rank-1分解）满足的形式为$\tilde{T}=\Lambda {\prod }_{⨂n}{v}^{(n)}$，易得，单秩张量的Tucker秩为（1，1，…，1）
——将单秩张量的形式稍作扩展，定义CP积，设N阶张量$T$为$R$个单秩张量的求和：
$$T=\sum _{r=0}^{R-1}{\tilde{T}}^{(r)}=\sum _{r=0}^{R-1}{\Lambda }^{(r)}\prod _{⨂n=0}^{N-1}{v}^{(r,n)}$$
注：从右向左是CP积
——将$\Lambda$看成向量，$v^{(r,n)}$看成矩阵，则上式的CP形式也可以写成如下的形式：
$$T_{a_1{a}_2...}=\sum _{r=0}^{R-1}{\Lambda }_r\prod _{⨂n=0}^{N-1}{V}_{r{a}_n}^{(n)}$$
其中,$\Lambda$为R 维向量， $V^{(n)}$为对应于第n个指标的矩阵，由R个列向量组成
——进一步改写CP积的形式
$$T_{a_1{a}_2...}=\sum _{r=0}^{R-1}{\Lambda }_r\prod _{⨂n=0}^{N-1}{V}_{r{a}_n}^{(n)}=\sum _{r=0}^{R-1}\sum _{r_1{r}_2...}{\delta }_{r,r_1{r}_2...}{\Lambda }_r\prod _{⨂n=0}^{N-1}{V}_{r{a}_n}^{(n)}$$
其中，$\delta$为（N+1）阶推广的超对角张量（N为张量$T$的阶数），第一个指标的维数为R，其余指标的维数与$T$指标的维数相等
——推广的超对角张量$\delta$满足：对于任意$r$（逗号之前的指标），张量${\delta }_{r,r_1{r}_2...}$满足：
δ r , r 1 r 2 . . . = { 0   其 它 情 况 下 1   当 r 1 = r 2 = . . . = r N ， ∀ r \delta_{r,r_1r_2...} = \{^{1\ 当r_1 = r_2 = ...=r_N}_{0\ 其它情况下}，\forall r δr,r1​r2​...​={0 其它情况下1 当r1​=r2​=...=rN​​，∀r
在这个情况下三阶CP积的图形表示为：

$$\Downarrow$$

——定义CP分解：
给定张量$T$,求解$\Lambda$和$V$，满足$T_{a_1{a}_2...}={\sum }_{r=0}^{R-1}{\Lambda }_r{\prod }_{⨂n=0}^{N-1}{V}_{r{a}_n}^{(n)}$
——定义CP秩：
使得CP分解严格成立的最小的R取值
——定义CP低秩近似：
给定R，使得CP分解式中左右之间的范数极小
CP分解中还存在尚未解决的关键科学问题：
①计算任意给定张量的CP秩是一个NP难问题
②目前存在多种CP低秩近似的算法，但是都无法保证获得的近似是最优的
③当取R = 1时，CP分解退化成rank-1分解；但CP分解中贡献最大的项（如$\Lambda$绝对值最大的元素对应的rank-1张量）不一定对应于rank-1分解所得的张量

4.重要内容总结

1. 张量的基本定义：是用指标标记的一堆数
2. 不同的标记方式由reshape、transpose等命令相互转换；
3. 张量的图形表示：张量本身由节点（圆圈或方块等）表示，指标由连接节点的边
   表示；不同张量的共有指标由连接对应节点的边表示，默认进行求和；
4. 本征值分解的形式与性质；
5. 最大本征问题对应的最优化问题；
6. 最大本征问题的幂级数求解法；
7. 奇异值分解的形式与性质，与矩阵秩的定义；
8. 基于奇异值分解的矩阵低秩近似问题，与裁剪误差；
9. 最大奇异值/向量的迭代算法；
10. 张量的单秩分解与最优单秩近似问题；
11. 高阶奇异值分解与Tucker秩；
12. 高阶奇异值分解的HOSVD算法；
13. 张量的Tucker低秩近似问题（HOOI）；
14. CP分解与CP秩