基于随机算法求解张量特征值

The Random Algorithm for Solving Eigenvalue of Tensors

摘要: 本文围绕张量H特征值问题展开研究。文中聚焦H特征对，将张量H特征值问题转化为非线性方程组问题，并利用三种Kaczmarz算法变体(NK、NRK、NURK)进行求解，这些算法的核心差异在于迭代过程中投影行的选择策略。为验证算法有效性，通过算例进行数值实验，结果显示三种算法均能有效求解张量最大H特征值，整体验证了算法对于求解张量特征值的可行性。

Abstract: This paper focuses on the H-eigenvalue problem of tensors. In this paper, the H-eigenvalue problem of tensors is transformed into a system of nonlinear equations by focusing on the H-eigenpair, and three nonlinear Kaczmarz algorithm variants (NK, NRK, NURK) are used to solve it. The core difference between these algorithms lies in the selection strategy of projection rows during the iteration process. To verify the effectiveness of the algorithms, numerical experiments are carried out with examples. The results show that all three algorithms can effectively solve the maximum H-eigenvalue of tensors, which generally verifies the feasibility of the proposed algorithms.

1. 引言

张量作为刻画多维数组与多维向量运算的基础数学工具，其特征值理论是张量分析的核心研究内容之一。张量特征值是通过特定运算规则从张量中求解得到的一类数值，其本质反映了张量在相应方向上所呈现的伸缩特性。在2005年之前，张量特征值领域的研究进展相对滞缓，尚未形成系统的理论框架。直至2005年，香港理工大学祁力群教授[1]与芝加哥大学林力行教授[2]分别独立完成了张量特征值及特征向量的严格定义，这一突破性成果为该领域的发展奠定了重要基础，不仅推动了张量特征值理论本身的快速完善，更促使其在超图谱理论[3] [4]、量子纠缠[5]、医疗影像等[6] [7]诸多跨学科领域获得了广泛应用，成为连接理论研究与实际问题的关键桥梁。

张量特征值问题是计算数学和张量分析领域的核心问题之一，在数据科学、信号处理、量子物理等领域有重要应用。由于其NP-hard的性质，开发高效且稳健的数值算法具有重要的理论和实践价值。本文将一个经典的迭代方法(Kaczmarz算法)应用于此问题，为该领域提供了一个新的求解思路。

对于$m$ 阶$n$ 维实张量$\mathcal{A}$ ，其元素用$a_{i_1}{{}_{i_2\cdots }}_{i_m}$ 表示，$n$ 维实向量$x$ 以及整数$0\le r\le m-1$ ，$\mathcal{A}$ 与$x$ 的$\left(m-r\right)$ 次积是一个$r$ 阶$n$ 维张量，记为$\mathcal{A}{x}^{m-r}$ ，其定义如下：

${\left(\mathcal{A}{x}^{m-r}\right)}_{i_1\cdots i_r}={\displaystyle \sum _{i_{r+1},\cdots ,i_m=1}^n{a}_{i_1\cdots i_m}{x}_{i_{r+1}}\cdots x_{i_m}},1\le i_1,\cdots ,i_r\le n$ (1)

特别地，$\mathcal{A}{x}^{m-1}$ 是一个$n$ 维向量，其分量为：

${\left(\mathcal{A}{x}^{m-1}\right)}_i={\displaystyle \sum _{i_2,\cdots ,i_m=1}^n{a}_{i{i}_2\cdots i_m}{x}_{i_2}\cdots x_{i_m}},1\le i\le n$ (2)

而$\mathcal{A}{x}^{m-2}$ 是一个$n\times n$ 的矩阵，其元素为：

${\left(\mathcal{A}{x}^{m-2}\right)}_{ij}={\displaystyle \sum _{i_3,\cdots ,i_m=1}^n{a}_{ij{i}_3\cdots i_m}{x}_{i_3}\cdots x_{i_m}},1\le i,j\le n$ (3)

当$\mathcal{A}$ 为半对称张量时，则$\mathcal{A}{x}^{m-1}$ 的梯度可以表示为：

${\nabla }_x\left(\mathcal{A}{x}^{m-1}\right)=\left(m-1\right)\mathcal{A}{x}^{m-2}$ (4)

在本文我们考虑的是H特征对。

若复数$\lambda \in ℂ$ 与非零向量$x\in {ℂ}^n$ 满足下述齐次多项式方程组

${\left(\mathcal{A}{x}^{m-1}\right)}_i=\lambda x_i^{m-1},\forall i=1,\cdots ,n$ (5)

则称$\lambda$ 为$\mathcal{A}$ 的H特征值，则称该解向量$x$ 为$\mathcal{A}$ 对应于特征值$\lambda$ 的特征向量。如果$\left(x,\lambda \right)$ 都是实的，则称为$\mathcal{A}$ 的H特征对。此外，若记$x^{\left[m-1\right]}=\left(x_1^{m-1},\cdots ,x_n^{m-1}\right)\in {ℂ}^n$ ，则(5)可简化表示为

$\mathcal{A}{x}^{m-1}=\lambda x^{\left[m-1\right]}$ (6)

2. 基于Kaczmarz算法求解张量H特征值问题

自Qi [1]与Lim 开创性地提出高阶张量特征值及特征向量的定义以来，非负不可约张量的最大特征值问题便迅速成为学术界的研究热点。Lim 与Chang等[8]学者通过严格定义不可约张量，成功将经典的Perron-Frobenius定理推广至非负不可约张量的最大H特征值研究领域，为该方向的理论发展奠定了重要基础。Ng等人[9]则创新性地设计了求解非负不可约张量最大H特征值的NQZ方法，为数值计算提供了有效工具。

在此基础上，Ni与Qi [10]敏锐地揭示了张量特征值问题与齐次多项式映射之间的对应关系——每个张量特征值问题均可等价转化为一个齐次多项式映射问题，这一突破性发现首次建立了张量特征值问题与非线性方程组之间的桥梁。基于此，他们将牛顿法巧妙应用于非负齐次多项式映射的最大H特征值求解，并严格证明了该方法具有二阶收敛性，为高效求解此类问题提供了理论保障与算法支撑。

在本文中将张量H特征值问题转化为非线性方程组问题，并利用随机Kaczmarz方法进行求解，数值结果表明算法可行。将(6)转化为非线性方程组

$F\left(x\right)=\mathcal{A}{x}^{m-1}-\lambda x^{\left[m-1\right]}$ (7)

函数$F\left(x\right)$ 对$x$ 是连续可微的。

若$\mathcal{A}$ 是半对称张量，则$F\left(x\right)$ 的雅可比矩阵为

$F^{\prime }\left(x\right)=\left(m-1\right)\mathcal{A}{x}^{m-2}-\lambda \left(m-1\right)D$ (8)

其中$D=diag\left(x^{\left[m-2\right]}\right)$ 。

对(7)采用非线性Kaczmarz算法进行求解。

经典Kaczmarz算法遵循矩阵行序逐次迭代，在方程组满足相容性条件时展现出优良的收敛性能。基于此，本文采用三种非线性Kaczmarz算法变体进行求解，其核心差异在于迭代过程中投影行的选择策略：

1. 非线性Kaczmarz算法(NK)：以方程组自身的行顺序依次确定投影行的选择；

2. 非线性随机Kaczmarz算法(NRK)：以和残差分量绝对值的平方呈正相关的概率确定投影行的选择；

3. 非线性均匀随机Kaczmarz算法(NURK)：以均匀分布概率随机确定投影行的选择。

为了将Kaczmarz算法应用到张量特征值问题上，本文对其进行一定的修改和调整。首先在每次迭代中，通过Kaczmarz方法更新特征向量，并根据当前特征向量更新特征值。

具体算法步骤如下：

1. 根据不同策略选择投影行$i$

1) 策略1：按照顺序依次选择

$i=k\left(\mathrm{mod}m\right)+1$

2) 策略2：按照和残差分量绝对值的平方呈正相关的概率选择

$i=\frac{F_i\left(x\right)}{{\Vert F\left(x\right)\Vert }^2}$

3) 策略3：按照均匀分布概率选择

$i=unidrnd\left(n,1,1\right)$

2. 计算该方程的梯度：

$g_i=\nabla F_i\left(x_k\right)$

3. 按照公式更新特征向量：

$x_{k+1}=x_k+\frac{r_i}{{\Vert g_i\Vert }^2}{g}_i$

4. 利用特征向量计算特征值的上下界，取上下界的平均值作为新的特征值${\lambda }_{k+1}$ 。

3. 数值实验

在本节中，我们将通过数值算例来验证算法1在求解张量最大H特征值时的有效性。针对所有测试问题，利用MATLAB内置的rand函数生成初始向量$x_0$ 。迭代停止条件设定为：当第$k$ 次迭代的残差误差满足$Res\le {10}^{-10}$ 时，终止迭代过程；若未满足该条件，但迭代次数已达到200次，也会停止迭代，并将此情形判定为未能成功求解该问题。

为评估算法的性能表现，对于每个测试问题，我们在100个不同的初始点上对所有方法进行了测试。测试所采用的评估指标包括平均迭代次数(记为“IT”)，平均迭代误差(记为“Res”)，以秒为单位的平均计算时间(记为“CPU”)以及出现最大H特征值的次数与实验总次数的比值(记为“Occur”)。

为对比NURK、NRK、NK 三种算法的性能，表1，表2，表3依次给出了各算法在三个不同张量上的数值计算结果。

算例1：张量$\mathcal{A}$ 是一个4阶2维的不可约非负对称张量，其元素为：

$\begin{array}{l}{a}_{1111}=a_{2222}=\frac{4}{\sqrt{3}},a_{1112}=a_{1121}=a_{1211}=a_{2111}=1,\\ a_{1222}=a_{2122}=a_{2212}=a_{2221}=1\end{array}$

其余元素$a_{i_1{i}_2{i}_3{i}_4}$ 均为0。

Table 1. Numerical results of Example 1

表1. 算例1数值结果

Method	IT	Res	CPU	Occur	λ
NURK	6.86	4.56E−12	0.0043	93.00%	6.3094
NRK	6.11	6.28E−12	0.0044	92.00%	6.3094
NK	6.79	7.89E−12	0.0048	96.00%	6.3094

算例2：张量$\mathcal{A}$ 是一个4阶2维的不可约非负对称张量，其元素为：

$a_{1112}=30,a_{1212}=1,a_{1222}=1,a_{2111}=6,a_{2112}=13,a_{2122}=37,$

其余元素$a_{i_1{i}_2{i}_3{i}_4}$ 均为0。

Table 2. Numerical results of Example 2

表2. 算例2数值结果

Method	IT	Res	CPU	Occur	λ
NURK	14.14	4.11E−11	0.0085	91.00%	41.0049
NRK	13.10	3.72E−11	0.0090	86.00%	41.0049
NK	14.13	4.34E−11	0.0090	92.00%	41.0049

算例3：张量$\mathcal{A}$ 是一个4阶2维张量，其元素为：

$a_{1111}=1,a_{1122}=a_{1221}=a_{1212}=a_{2121}=a_{2211}=a_{2112}=\frac{1}{3},a_{2222}=1$

其余元素$a_{i_1{i}_2{i}_3{i}_4}$ 均为0。

Table 3. Numerical results of Example 3

表3. 算例3数值结果

Method	IT	Res	CPU	Occur	λ
NURK	6.70	8.33E−12	0.0042	100%	2
NRK	5.98	6.31E−12	0.0042	100%	2
NK	6.60	6.10E−12	0.0047	100%	2

4. 总结

数值实验结果表明，三种算法都能够求解出张量最大H特征值，展现出各自的优势。其中，NURK算法在平均迭代时间(CPU)方面较为突出；NRK算法在平均迭代次数(IT)方面表现较为突出，能够相对高效地收敛到结果；NK算法在最大H特征值出现的占比(Ocurr)上具备一定优势，反映出其在特定场景下对目标特征值的捕捉能力。这不仅验证了所提算法在求解张量H特征值问题上的可行性，也为该领域的数值计算方法提供了新的思路和参考。

然而，本研究亦存在一定局限性。例如，实验主要围绕4阶2维张量展开，对于更高阶数、更大维度的张量，算法的性能及适用性有待进一步验证。同时，在算法效率提升方面仍有优化空间，例如如何进一步降低计算复杂度，提升算法在大规模张量计算中的效率。