Shifted Symmetric Higher-Order Power Method for Computing Eigenvalues of Multiple Order Symmetric Tensors

本节首先回顾求解对称张量Z-特征值的带位移高阶幂法，随后给出求解不同阶对称张量组特征值的带位移高阶幂法.在此之前，先给出一个记号.设A∈R^[m，n]，x=(x₁，…，x_n)^T∈ℝⁿ，记

文献[16]提出了如下求解对称张量Z-特征值的带位移高阶幂法(SS-HOPM)：

算法1^[16]  求解对称张量Z-特征值的带位移高阶幂法(SS-HOPM).

输入：m阶n维对称张量A，位移α∈ℝ.

输出：A的Z-特征值λ_k+1.

初始化：任取非零向量x₀∈ℝⁿ，使得‖x₀‖₂=1.令λ₀=Ax₀^m.

for k=0，1，…，do

       x_k+1←Ax_k^m-1+αx_k

       x_k+1←x_k+1/‖x_k+1‖₂

       λ_k+1=Ax_k+1^m

end for

考虑(1)式，显然有

为了求λ，构造辅助函数

其中α∈ℝ，x∈Ω={x∈ℝⁿ：x^Tx=1}.由文献[16]的引理3.1知(2)式中f(x)的梯度为

由x∈Ω知

下面给出求解不同阶对称张量组特征值问题(1)的带位移高阶幂法：

算法2  求解不同阶对称张量组特征值的带位移高阶幂法.

输入：张量A₃，A₂，位移α，允许误差ε.

输出：{A_i}_i=2³的特征对(λ_k+1，x_k+1).

初始化：任取非零向量x₀∈ℝⁿ，使得‖x₀‖₂=1.令$f\left( {{\mathit{\boldsymbol{x}}_0}} \right) = \frac{1}{3}{\mathit{\boldsymbol{A}}_3}\mathit{\boldsymbol{x}}_0^3 + \frac{1}{2}{\mathit{\boldsymbol{A}}_2}\mathit{\boldsymbol{x}}_0^2 + \frac{1}{3}\alpha $.

for k=0，1，…，do

       x_k+1←A₃x_k²+A₂x_k+αx_k

       x_k+1←x_k+1/‖x_k+1‖₂

       f(x_k+1)=$\frac{1}{3}{\mathit{\boldsymbol{A}}_3}\mathit{\boldsymbol{x}}_{k + 1}^3 + \frac{1}{2}{\mathit{\boldsymbol{A}}_2}\mathit{\boldsymbol{x}}_{k + 1}^2 + \frac{1}{3}\alpha $

       if |f(x_k+1)-f(x_k)|＜ε，且‖x_k+1-x_k‖₂＜ε

             break.

       end if

end for

本节讨论算法2的收敛性.首先给出一些引理.

引理1^[16]  可微函数f(x)：Υ⊆ℝⁿℝ是凸函数当且仅当Υ是凸集，且对任意x，y∈Υ均有

引理2^[16]  二阶可微函数f(x)：Υ⊆ℝⁿ→ℝ是凸函数当且仅当Υ是凸集，且f(x)的Hessian矩阵在Υ上是半正定的，即对任意x∈Υ，$\nabla$²f(x)是半正定矩阵.

定理1  设A₃∈R^[3，n]是对称张量，A₂∈R^[2，n]是对称半正定矩阵，ρ(A₂)是A₂的谱半径，Φ(A₃)=$2\mathop {\max }\limits_{x \in \mathit{\Omega }} \rho \left( {{\mathit{\boldsymbol{A}}_3}\mathit{\boldsymbol{x}}} \right)$，若α＞Φ(A₃)+ρ(A₂)，则由算法2得到的序列$\left\{ {f\left( {{\mathit{\boldsymbol{x}}_k}} \right)} \right\}_{k = 0}^\infty $是单调不减序列且有上界，因此该序列有极限f(x_*)，即$\mathop {\lim }\limits_{k \to \infty } $f(x_k)=f(x_*)，x_*∈Ω.进一步，有

即(λ_*，x_*)是不同阶对称张量组{A₃，A₂}的特征对.

证  由文献[16]中引理3.3知，f(x)的Hessian矩阵为

其中I表示单位矩阵.对任意y∈Ω和$\mathit{\boldsymbol{\bar x}} = \frac{\mathit{\boldsymbol{x}}}{{{{\left\| \mathit{\boldsymbol{x}} \right\|}_2}}} \in \mathit{\Omega }$，由

知，f(x)的Hessian矩阵H(x)是半正定矩阵.因此，由引理2知f(x)是ℝⁿ上的凸函数.

由α＞Φ(A₃)+ρ(A₂)知，对任意x∈Ω，$\nabla$f(x)≠0.对x_k，x_k+1∈Ω且${\mathit{\boldsymbol{x}}_{k + 1}} = \frac{{\nabla f\left( {{\mathit{\boldsymbol{x}}_k}} \right)}}{{{{\left\| {f\left( {{\mathit{\boldsymbol{x}}_k}} \right)} \right\|}_2}}}$，由引理1得

故序列$\left\{ {f\left( {{\mathit{\boldsymbol{x}}_k}} \right)} \right\}_{k = 0}^\infty $是单调不减的.又由x_k∈Ω和

知f(x)是有界函数.因此序列$\left\{ {f\left( {{\mathit{\boldsymbol{x}}_k}} \right)} \right\}_{k = 0}^\infty $有极限f(x_*)，即有

因为x_*是f(x)的极值点，所以$\nabla$f(x_*)=0.因此x_*是不同阶对称张量组{A₃，A₂}的特征向量，λ_*=A₃x_*³+A₂x_*²是相应的特征值，故(λ_*，x_*)是一个特征对.

在下一节的数值试验中，可以选取$\alpha = 2\sum\limits_{i, j, k \in [n]} {\left| {{a_{ijk}}} \right|} + \rho \left( {{\mathit{\boldsymbol{A}}_2}} \right)$.原因是∀y∈Ω，由三角不等式有

所以

这保证了α大于Φ(A₃)+ρ(A₂)，满足定理1中的收敛条件.事实上，选取适当的α可以使算法收敛得更快，详见下一节的数值例子.

本节给出3个数值算例验证算法2的有效性.共作100次数值试验，每次试验随机选取不同的初始向量x₀，且元素服从[-1, 1]的均匀分布，迭代终止的条件为ε=10^-8或达到最大迭代步数4 000. Occurrences表示在进行的100次试验中求出特征值的次数，Iter表示平均迭代次数.

例1  设A₃=(a_ijk)∈R^[3,3]是对称张量，A₂=(b_ij)∈R^[2,3]是对称正定矩阵，且

求解不同阶对称张量组{A₃，A₂}的实特征对等价于求解非线性方程组

用Mathematic软件中的NSolve命令求解方程组(3)，得到14个不同的实特征对，见表 1.

下面验证算法2的有效性.数值结果见表 2.

由表 2知，在进行的100次数值试验中，不带位移(α=0)的算法2计算出了表 1中的两个特征值.而算法2在分别取α=1和$\alpha = 2\sum\limits_{i, j, k \in [3]} {\left| {{a_{ijk}}} \right|} + \rho \left( {{\mathit{\boldsymbol{A}}_2}} \right)$=9.357 7时，100次试验都能有效地找到14个特征值中的4个.

例2  设A₃=(a_ijk)∈R^{[3, 3]}是对称张量，A₂=(b_ij)∈R^{[2, 3]}是对称矩阵，且

用NSolve命令求解方程组(3)，得到不同阶对称张量组{A₃，A₂}的8个不同的实特征对，见表 3.

下面验证算法2的有效性.数值结果见表 4.

由表 4知，在进行的100次数值试验中，不带位移(α=0)的算法2仅有65次收敛到特征值λ=2.750 1，其余35次在迭代到4 000步时计算不出特征值.而算法2在分别取α=3和$\alpha = 2\sum\limits_{i, j, k \in [3]} {\left| {{a_{ijk}}} \right|} + \rho \left( {{\mathit{\boldsymbol{A}}_2}} \right)$=28.735 6时，100次试验都能有效地找到8个特征值中最大的2个.

图匹配中的不同阶对称张量组大多是大型稀疏张量.下面给出高维稀疏不同阶对称张量组的例子用以检验算法2的有效性.

例3  设A₃=(a_ijk)∈R^{[3, 100]}和A₂=(b_ij)∈R^{[2, 100]}是随机生成的稀疏对称张量，非零元的个数都是150个，其中每个元素服从[0, 1]的均匀分布.由算法2得到的数值结果见表 5.

由表 5可见，在进行的100次数值试验中，不带位移(α=0)的算法2在迭代到4 000步时计算不出特征值.在算法2中分别取α=2和$\alpha = 2\sum\limits_{i, j, k \in [100]} {\left| {{a_{ijk}}} \right|} + \rho \left( {{\mathit{\boldsymbol{A}}_2}} \right)$=28.803 1，每次试验均有效地求出特征值λ=1.157 9或λ=1.157 5，且α=2时的迭代次数比α=28.803 1时的迭代次数少.例3表明，算法2在高维情形下效果也较好.

以上数值例子表明，带位移的算法2对任意对称张量A₃和A₂构成的不同阶对称张量组{A₃，A₂}都是有效的.由文献[16]的引理4.1知$2\sum\limits_{i, j, k \in [n]} {\left| {{a_{ijk}}} \right|} > \mathit{\Phi }\left( {{\mathit{\boldsymbol{A}}_3}} \right)$.取$\alpha = 2\sum\limits_{i, j, k \in [n]} {\left| {{a_{ijk}}} \right|} + \rho \left( {{\mathit{\boldsymbol{A}}_2}} \right)$，算法2收敛，此时的α远大于Φ(A₃)+ρ(A₂)，因此取适当的α，算法2计算得更快.而且在100次试验中总能计算出不同阶对称张量组{A₃，A₂}的最大特征值.

本文给出了求解不同阶对称张量组{A_i}_i=2³特征对的带位移高阶幂法，即算法2.该算法通过构造辅助函数f(x)将问题(1)中求特征向量问题转化为求辅助函数的极值点问题.然后，将此辅助函数f(x)的梯度$\nabla$f(x)作为迭代格式，使得算法2收敛.与不带位移(α=0)的算法相比，带位移的算法2能更有效地求出不同阶对称张量组{A_i}_i=2³的实特征对.但是α的选取对迭代次数影响很大，如何选取最佳的α将是我们下一步的研究工作.