单侧障碍问题的交替方向乘子法.pdf

2024 年 4 月第 45 卷第 2 期湘南学院学报Journal of Xiangnan UniversityApr.,2024Vol.45 No.26单侧障碍问题的交替方向乘子法熊桂花，张守贵*（重庆师范大学 数学科学学院，重庆 401331）摘 要：基于最优罚参数的交替方向乘子算法可用于一类单侧障碍自由边界问题的数值求解。单侧障碍问题经过差分离散化后转变为有限维互补问题。表示区域位移的辅助变量和增广 Lagrangian 函数的引入将原问题等价转化为鞍点问题，然后交替方向乘子算法可用于求解。对此算法的辅助变量进行消除，可导出交替方向乘子算法的纯对偶算法，进而进行收敛分析和最优罚参数的讨论。数值算例结果验证了该算法的可靠性和有效性。关键词：单侧障碍问题;交替方向乘子算法;最优罚参数;增广 Lagrangian 函数中图分类号：O241.82 文献标志码：A DOI:10.3969/j.issn.1672-8173.2024.02.002收稿日期：2023-10-22基金项目：国家自然科学基金（11971085）；重庆市自然科学基金（cstc2020jcyj-msxmX0066）；重庆市研究生教育教学改革研究项目（yjg213071）；重庆市研究生科研创新项目（CYS22561）。作者简介：熊桂花（1997），女，重庆忠县人，硕士生，研究方向为微分方程的数值解法；通信作者：张守贵（1973），男，四川泸县人，教授，博士，研究方向为微分方程的数值解法。单侧障碍问题不同于普通的线性问题，它是用微分不等式和等式的组合形式表示的 1-2。对于这类定义在区域上的非线性问题，理论研究与数值计算都面临很大的挑战 3。交替方向乘子法作为一种重要的算子分裂算法，其收敛性和实际计算已得到深入研究，有很好的理论基础，且其应用十分广泛 4-9。交替方向乘子法求解变分不等式的主要想法是通过引入一个辅助变量，将问题的线性部分和非线性部分分开求解 10-13。交替方向乘子法计算简单，但该方法的收敛速度高度依赖于罚参数的取值，所以最优罚参数选择问题亟待解决。本文对交替方向乘子算法中的最优罚参数给出了近似计算方法。首先把单侧障碍问题用差分法离散化，得到有限维的交替方向乘子法。然后消除此算法的辅助变量，推导出交替方向乘子算法的纯对偶算法，对此算法进行收敛性分析，并得到用于数值计算的最优罚参数的近似计算公式，最后用数值模拟检验该最优罚参数对提升交替方向乘子算法性能的效果。1 单侧障碍问题和有限维交替方向乘子法将单侧障碍问题用差分法离散化为有限维形式，导出用矩阵和向量表示的问题，得到求解该问题的交替方向乘子法。1.1 单侧障碍问题的离散化形式二维平面单侧障碍问题可用式（1）表示为：()()()()()()()()0,()0u xf x xu xx xu xf xu xxxu xx=,（1）其中：是有界正方形区域，=为其边界，()f x 是 上的已知函数，()x是 上的已知障碍函数，()12,xx x=，222212/xx=+。将式（1）的问题利用等步长五点差分离散化，并代入已知函数()f x 和()x在网格结点上的函数值，假设区域内部网格结点数为 N，则得有限维单侧障碍问题，即700(,)0=AufuAuf u （2）其中：N 维列向量u 表示待求函数()u x 在网格结点上的差分解，向量 f 和 为已知函数()f x 和()x在网格结点上的值构成的 N 维列向量，N 阶稀疏矩阵 A 是五对角（分块三对角）且对称正定的。由文献 3 可知，有限维单侧障碍问题式（2）等价于极小值问题：求NuR，满足min(),s.t.Juu （3）其中，TT1()2J=uu Auf u。由于不等式（3）约束优化问题有且仅有一个解，因而式（2）的问题存在唯一解。1.2 有限维交替方向乘子法为了求得式（2）的数值解，引入一个辅助变量NpR，令=up。再引入集合,NC=pRp 和符号函数。则式（3）等价于，由此导出增广 Lagrangian 函数，即 （4）其中：罚参数0,NNNr vRqRR。鞍点问题可由式（4）表示为 （5）由文献 10 可知，式（5）的鞍点问题有以下结果。引理 1 10 若,u p 为在NNNRRR 上的一个鞍点，则u 是式（3）的解，并且=up。根据引理 1，要得到原式（3）的解，改为求解式（5）的鞍点,u p。由此出发，可用交替方向乘子算法求该鞍点。该算法只需求解一个线性方程组，得到(1)k+u，且(1)k+p和(1)k+均可显式求解 7，引入记号()max(,)+=0 xx，得算法 1。算法 1交替方向乘子算法：第一步，取初始值00,0,NNrpRR 令0k=；第二步，求解式（6），得(1)kN+uR；(1)(1)()()kkkkrr+=+uAufp （6）第三步，按式（7）显式计算(1)kN+pR；(1)(1)()()(1)1()kkkkkrr+=puu （7）第四步，按式（8）确定(1)kN+R；(1)()(1)(1)+()kkkkr+=pu （8）第五步，如果满足迭代终止条件，停止迭代并得到数值解(1)k+u，反之则令1kk=+，返回到第二步。2 改进交替方向乘子法消除交替方向乘子算法辅助变量，可推导出纯对偶的交替方向乘子算法，然后得到算法的收敛性分析结果。本文提出一种用于最优罚参数的数值近似的计算方法，得到基于最优罚参数的交替方向乘子法，对算法进行改进 1316。2.1 纯对偶交替方向乘子法考虑算法1，消去算法中的辅助变量。首先，令()(1)()(),kkkr+=u()(1)()min(,()kkkr=0u，从式（7）中得熊桂花，张守贵：单侧障碍问题的交替方向乘子法8湘南学院学报2024 年 4 月（第 45 卷）第 2 期(1)(1)()(1)()(1)()()kkkkkkkrrr+=puu （9）将式（9）代入式（8）得(1)()(1)()(1)kkkkk+=+=，从而有()()kk+=（10）其次，由式（8）有()()()()(1)()()()(1)()()(1)()=()kkkkkkkkkkkkrrrrrr=+=+ppuuupuu，可得()()()()(1)2kkkkkrr+=+pu （11）由于()()(1)()()()(1)()(1)()()()(1)()2,()02=2()2,()0kkkkkkkkkkkkkkrrrrrrr+=+uuuu，故由式（10）得()(1)()()()(1)()(1)(),()02+,()0kkkkkkkkkrrrrr+，令0k=；第二步，求解式（12）线性方程组，得(1)kN+uR；(1)()()()kkkrr+=+IA uf （12）第三步，按以下方式计算(1)(1),kkN+R；(1)()(1)(1)()(1)()()kkkkkkrr+=uu （13）第四步，如果满足迭代终止条件，停止迭代并得数值解(1)k+u，反之则令1kk=+，返回到第二步。从式（12）可得：(1)1()()()kkkrr+=+uAf （14）其中rr=+AIA。将式（14）代入到式（13）消去(1)k+u，令1()rr=+aAf，得到纯对偶的算法 3。算法 3纯对偶的交替方向乘子算法：第一步，取定0+，0，0r，0k=；第二步，按式（15）确定(1)(1),kkN+R；(1)1()1()(1)1()1()()()()()kkkrrkkkrrrrrrrr+=+=+AIAaAIAa （15）第三步，如果满足迭代终止条件，停止迭代得(1)(1),kk+，并由式（14）得数值解(2)k+u，反之则令1kk=+，返回到第二步。2.2 收敛性分析因为min(0,)max(0,)xx=，则式（15）可以改写成(1)1()1()(1)1()1()()()()()kkkrrkkkrrrrrrrr+=+=AIAaAIAa。因为+xx，NxR，所以得到(1)1()1()(1)1()1()()()()()kkkrrkkkrrrrrrrr+AIAaAIAa。将此两式相加得(1)(1)1()1()1()1()()()2()kkkkkkrrrrrrrrr+IAAIAAa （16）对 任 意 向 量12=yyy，1NyR，2NyR，则2NyR。定 义12=+yyy，则 由 式（16）得(1)()2()kkr+Ma，其中(1)(1)(1)kkk+=，1111 ()rrrrrrrr=IAAMIAA。9综合上述分析，矩阵 M 谱半径决定了交替方向乘子法的收敛性，且关于 M 的非零特征值，我们有以下的引理。引理 2 M 的非零特征值是12rrIA 的特征值，反之也成立。证明：令11,rr=MIA12rr=MA，则有1212 =MMMMM 。我们需证明：M 的非零特征值是1122rr=MMIA 的特征值，反之亦然。令为 M 的非零特征值，(，)为其对应的特征向量，则得 （17）由式（17）中两式相加得。由于非零，所以有 （18）将式（18）代入到式（17）中可得。由此可见，为12MM 的特征值。由以上证明过程易得1122rr=MMIA 的非零特征值，也是 M 的特征值。综上，引理得证。由引理2可知，通过考察矩阵12rrIA 的谱半径可以研究交替方向乘子法的收敛性，从而可得如下定理。定理 1 交替方向乘子法对任意正的罚参数r 都收敛。证明：考虑此序列(1)1()()1()(2)2kkkkrrrr+=IAA，两边同时乘以1A 得1(1)1()11()2kkkrr+=AAA A （19）因为11111()()rrr=+=+AAIIAA （20）将式（20）代入式（19）中得1(1)1()1111()2()kkkrr+=+AAAIAA （21）引入辅助变量(1)1(1)kk+=A，那么，式（21）变为(1)111()(2(),0kkrrk+=+IAIA。综合上述分析，序列()k 的收敛性决定了序列()k 和,kk+的收敛性，从而由式（14）决定了序列()ku的收敛性。为 了 证 明 序 列()k 的 收 敛 性，令1112()rr=+CIAIA，由 此 可 见，C 为 关 于1A形 如1()12(1)p trtrt=+的矩阵多项式函数，从而C 的特征值为1A 特征值的多项式函数1()()iip=CA，即1111()1()()12,1,2,1()1()iiiiirriNrr=+AACAA （22）由于1A 的特征值和特征向量 满足,=A （23）由于 A 是对称正定的，因此TT0=A，再由式（22）可得()1,1,2,iiN，令0k=；第二步，对1,2,iN=，把满足()()1()0NkkiijjijurA uf=，其中2212xx=+。由解析解可确定相应的 Dirichlet 边界条件。用最优罚参数计算方法，得到该问题在取步长 h=3/80 下的最优罚参数为113.069=r，用这个罚参数取值的交替方向乘子法进行计算，图 4 为数值结果，图 5 为自由边界（即接触区域0u=的边界）数值解的结果，显然，其与解析解的结果是相吻合的。图 4 算例 2 数值解图 5 算例 2 自由边界上的数值解图 6 为交替方向乘子算法（步长 h=3/40）选取不同罚参数的算例 2 迭代次数。从数值结果（表 2）同样发现：最优罚参数57.191=r的收敛速度较快，迭代次数较少；图 3 算例 1 迭代次数比较表 1 算例 1 两种算法性能比较h交替方向乘子法半光滑牛顿法迭代次数CPU 时间/s迭代次数CPU 时间/s1/10220.02480.0241/20320.03090.0321/40280.056100.0781/80410.410130.9801/160633.5371311.673图 6 算例 2 算法迭代次数比较表 2 算例 2 两种算法性能比较h交替方向乘子法半光滑牛顿法迭代次数CPU 时间/s迭代次数CPU 时间/s3/10240.02470.0333/20280.02980.0473/40330.079110.2283/80560.6061522.1333/160984.790151 195.969熊桂花，张守贵：单侧障碍问题的交替方向乘子法12湘南学院学报2024 年 4 月（第 45 卷）第 2 期交替方向乘子法比半光滑牛顿法需要的 CPU 时间更少。4 结论本文研究了求解单侧障碍问题的基于最优罚参数的交替方向乘子算法，并推导出了交替方向乘子算法的纯对偶算法，以此证明了交替方向乘子算法对于任意正的罚参数收敛。经过推导发现，最优罚参数取决于迭代矩阵的最大和最小特征值，从而得到最优罚参数的计算公式。数值实验表明，本文提出的最优罚参数是对交替方向乘子算法的有益改进。参考文献 1 韩渭敏,程晓良.变分不等式简介:基本理论、数值分析及应用 M.北京:高等教育出版社,2007.2 XUE L,CHENG X L.An algorithm for solving the obstacle problems J.Computers&Mathematics with Applications,2004,48(10/11):1651-1657.3 韩继业,修乃华,戚厚铎.非线性互补理论与算法 M.上海:上海科学技术出版社,2006.4 唐瑜,张守贵.可分离二次规划问题的自适应交替方向乘子法 J.重庆理工大学学报(自然科学),2022,36(5):284-289.5 赵建喜,易丹辉,刘晓倩.去除稀疏和结构化噪声的交替方向增强拉格朗日算法 J.数学的实践与认识,2017,47(16):164-170.6 GLOWINSKI R,LE T P.Augmented Lagrangian and Operator-Splitting Methods in Nonlinear Mechanics M.Philadelphia:Society for Industrial and Applied Mathematics,1989.7 GLOWINSKI R,MARINI L D,VIDRASCU M.Finite-element approximations and iterative solutions of a fourth-order elliptic variational inequality J.IMA Journal of Numerical Analysis,1984,4(2):127-167.8 LU X L,L X B.ADMM for image restoration based on nonlocal simultaneous sparse Bayesian coding J.Signal Processing:Image Communication,2019,70:157-173.9 苑雨薇.基于矩阵填充的抗癌药物敏感性预测研究 D.徐州:中国矿业大学,2021.10 张霖森,程兰,张守贵.曲率障碍下四阶变分不等式的交替方向乘子法 J.应用数学和力学,2023,44(5):595-604.11 郭楠馨.求解单侧障碍问题的自适应 Uzawa 块松弛算法 D.重庆:重庆师范大学,2020.12 KOKO J.Uzawa block relaxation method for the unilateral contact problem J.Journal of Computational and Applied Mathematics,2011,235(8):2343-2356.13 ESSOUFI E H,KOKO J,ZAFRAR A.Alternating direction method of multiplier for a unilateral contact problem in electro-elastostatics J.Computers&Mathematics with Applications,2017,73(8):1789-1802.14 CHORFI A,KOKO J.Alternating direction method of multiplier for the unilateral contact problem with an automatic penalty parameter selection J.Applied Mathematical Modelling,2020,78:706-723.15 DEMMEL J W.Applied Numerical Linear Algebra M.Philadelphia:Society for Industrial and Applied Mathematics,1997.16 史荣昌,魏丰.矩阵分析 M.第 3 版.北京:北京理工大学出版社,2010.17 STADLER G.Path-following and augmented Lagrangian methods for contact problems in linear elasticity J.Journal of Computational and Applied Mathematics,2007,203(2):533-547.