应用改进的Levenberg-Marquardt方法求解一类多线性系统

王 丽，陈 震，刘奇龙

（贵州师范大学数学科学学院，贵州贵阳550025）

1 预备知识

设A是一个m阶n维的张量，表示为

其中，［n］＝｛1，2，…，n｝.称张量A是对称的，如果ai1i2…im＝aπ（i1i2…im），∀π∈Πm，其中Πm是指标向量［1，2，…，m］的所有排列的全体.称张量A是半对称的，如果

其中Πm-1是指标向量［2，3，…，m］的所有排列的全体.

本文考虑如下多线性系统的数值求解问题

其中，x＝［x1，x2，…，xn］T∈Rn是未知向量，Axm-1是列向量，它的第i个元素定义为

多线性系统（1）来源于数据挖掘、微分方程数值解和张量补等实际问题［1-3］，近年来受到广泛关注.2015年，Li等［1］提出了求解数据挖掘背景下产生的系数为稀疏非负张量的多线性系统的一种迭代方法，并证明了算法的线性收敛性；2016年，Ding等［2］研究了系数张量为M-张量的多线性系统，在张量A是非奇异的M-张量，b是正向量的情形下证明了方程有唯一的正解，并提出了求解张量方程的几种迭代算法，推广了经典的Jacobi方法与Gauss-Seidel方法；Han［4］提出了确定这个唯一正解的同伦方法；He等［5］提出了二次收敛的Newton型算法；2017年，Li等［6］将一些经典的分裂方法推广到求解对称张量方程，在适当的条件下证明了算法的全局收敛性和局部r-线性收敛性；对于系数为强M-张量的张量方程，Liu等［7］提出了一些基于张量分裂的算法，并将该算法应用到高阶马尔可夫链模型；Lü 等［8］将经典的 Levenberg-Marquardt（LM）方法应用到系数为半对称张量的多线性系统，并证明了在局部误差界下的全局收敛性与局部二次收敛性；Li等［9］提出了求解多线性系统的混合交替投影算法，并证明了在适当条件下算法的局部线性收敛性.

综上，现有方法通常是针对系数为对称或半对称的张量或具有一定结构的非奇异M-张量的多线性系统进行讨论.本文则是考虑一般情形下的多线性系统（1）.首先讨论了应用改进的LM方法求解多线性系统的数值算法，然后证明了该方法在局部误差界条件下的全局收敛性和局部二次收敛性，最后通过数值算例验证了该方法的有效性.

2 迭代算法

3 收敛性分析

4 数值算例

选取几个例子验证算法2.2的有效性.所有程序在配置为intel（R）Core（TM）i5-6200U CPU @2.30 GHz的笔记本电脑环境下使用Matlab 2015b编写，涉及到张量计算的部分使用了工具箱Tensor toolbox 2.5［16］.迭代过程中的容许误差取为，其中算法2.2相关参数的选取为：α0＝0.5，αmin＝1.0×10-5，a1＝4，a2＝1/4，h0＝0.02，h1＝0.3，h2＝0.6，最大迭代次数N＝20.

例4.1采用文献［8］中的例5.1，考虑多线性系统（1）.首先生成一个m阶n维随机非负张量A，其元素在（0，1）上均匀分布.为了得到向量b，先选取x*＝2*ones（n，1）为准确解.在实验中选取初值x0＝3*ones（n，1），通过算法2.2 计算多线性系统的解，并与文献［8］的算法相比较.数值结果如表1所示，其中，k表示迭代次数，t表示运算时间，∈res表示相对误差‖x-x*‖/‖x*‖.因为这里的张量A没有对称性，所以不能直接使用文献［8］的LM方法.为了与文献［8］中的LM方法比较，在计算Jacobi矩阵时不能利用文献［8］中的（3.4）式，而是利用本文的（8）式.由表1可以看出：算法2.2与LM方法相比较，前者所需的迭代次数更少，时间更短.

表1 选取不同阶数不同维数的多线性系统的数值结果Tab.1 The numerical results of multilinear systems with different orders and dimensions

图1是选取m＝4，n＝30时，算法2.2的收敛性演示.显然，在局部误差界条件下，算法2.2产生的点列是二次收敛的.

下面利用数值实验说明该方法可应用于如下广义的多线性系统

其中，Ap是一个p阶n维的张量，p＝2，3，…，m.

例4.2随机产生一个三阶张量A1和二阶矩阵A2，其元素都在（0，1）上均匀分布，选取

图1 算法2.2的收敛性演示Fig.1 Convergence demonstration of Algorithm 2.2

为准确解，从而确定右端向量b.利用算法2.2求解上述的多线性系统，结果如表2所示.表2表明算法2.2同样可以有效求解广义的多线性系统.

表2 选取不同维数的广义多线性系统的数值结果Tab.2 The numerical results of generalized multilinear systems with different dimensions

例4.3参考文献［2］中的例4.3，考虑引力作用下的质点运动方程

Dirichlet边界条件为

其中重力常数G≈6.67×10-11Nm2/kg2和地球质量M≈5.98×1024kg.

离散化后可以得到如下的多项式方程组

上述方程组可以改写成多线性系统Ax3＝b，其中系数A是一个四阶张量，其元素为

右端向量b中的元素为

另一方面，质点运动轨迹可以用抛物线近似地表示

其中，c0和c1表示地球半径，等于6.37×106km，重力加速度g≈9.8 m/s2，α和β是由边界条件确定的常数.

图2是利用算法2.2求解上述多线性系统的数值结果与抛物线的对比，显然，数值结果完全符合实际情况.

图2 在引力作用下质点运动的轨迹Fig.2 Trajectory of particle motion under the action of gravity