对于非凸非光滑不可分离问题的并行式惯性近端非精确ADMM算法
Parallel Inertial Proximal Inexact ADMM Algorithm for Nonconvex Nonsmooth and Nonseparable Problems
摘要: 交替方向乘子法(ADMM)作为一种简单有效的求解可分离问题方法,已经被广泛应用于信号处理、并行计算、机器学习等许多领域。但是当目标函数较为复杂时,ADMM算法的收敛性就会难以保证。针对一类复杂的非凸非光滑不可分离的优化问题,提出了一种非精确的并行式惯性近端交替乘子法(PIPI-ADMM),在简单的参数假设下,算法产生的迭代序列收敛到增广拉格朗日函数的稳定点。另外,当增广拉格朗日函数满足Kurdyka-Lojasiewicz性质时,算法具有强收敛性。最后,数值实验结果表明,算法的迭代策略对求解具有显著的加速作用。
Abstract: The Alternating Direction Method of Multipliers (ADMM) has been widely applied in various fields such as signal processing, parallel computing, and machine learning due to its simplicity and effectiveness in solving separable problems. However, the convergence of the ADMM algorithm becomes challenging to guarantee when dealing with complex objective functions. For a class of complex, non-convex, non-smooth, and non-separable optimization problems, a novel Parallel Inexact Proximal Inertial Alternating Direction Method of Multipliers (PIPI-ADMM) is proposed. Under reasonable assumptions, the iterative sequence converges to a stationary point of the augmented Lagrange function. Moreover, with the help of Kurdyka-Lojasiewicz property, the proposed algorithm is strongly convergent. The numerical results show that the iterative strategy has significant acceleration.
文章引用:杨开元, 党亚峥, 杨灿. 对于非凸非光滑不可分离问题的并行式惯性近端非精确ADMM算法[J]. 理论数学, 2024, 14(12): 95-111. https://doi.org/10.12677/pm.2024.1412411

1. 引言

本文考虑如下问题:

min x,y,z f( x )+g( y )+h( x,y,z ) s.t.Ax+By+Cz=b, (1.1)

其中, x r ,y s ,z t 是函数变量, f: m + g: n + 是下半连续函数(可能非凸非光滑), h: r × s × t 是连续可微函数, A m×s ,B m×s ,C m×t 是给定的矩阵, b m 是给定的向量。问题(1.1)的增广拉格朗日函数(ALF)为

L β ( x,y,z,λ )=f( x )+g( y )+h( x,y,z ) λ, A 1 x 1 + A n x n +Byb + β 2 Ax+By+Czb 2 , (1.2)

其中, λ 是拉格朗日乘子, β>0 是罚参数。

问题(1.1)是一类经典的具有线性约束的优化问题,广泛应用于压缩感知、稀疏信号恢复、图像处理、机器学习、矩阵分解等多个领域。在实际应用中,许多优化问题涉及非凸非光滑函数。非凸函数通常具有多个最小值和鞍点,这使得求解过程变得十分困难;非光滑函数往往在一些点不可微,甚至无法求导。因此,研究非凸非光滑情形下优化问题的求解方法具有重要的理论和实践意义。

交替方向乘子法(ADMM)作为一种高效求解具有多变量的约束优化算法,因其灵活度高,简单,处理速度快而被广泛使用于图像处理、统计学习等领域[1] [2]。在问题(1.1)中,当目标函数为凸函数时,ADMM算法及其扩展研究已经较为完善[3]-[5]。近年来,人们更加关注目标函数非凸时ADMM算法的收敛性以及求解速度问题。一些关于两块非凸问题的ADMM算法的收敛性研究已经被证明[6]-[8],但随着变量的增加,对于三块和多块非凸问题,ADMM算法的收敛性证明非常具有挑战性。陈等[9]提出了一个反例说明了直接将ADMM算法从两块问题扩展到三块问题是不收敛的。因此许多学者通过一些更强的假设或者对算法进行适当的改进以提高算法的收敛性。邓等[10]发现,当矩阵 A i 是正交的并且列满秩时,通过扩展两块到多块的雅各比配置可以保证算法的收敛性,王等[11]通过在子问题中加入Bregman修正,来证明非凸情况下的收敛性。

与此同时,随着数据维度以及区块数量逐渐增加,基于ADMM的加速算法逐渐被深入研究。一些学者将惯性技术与ADMM算法结合起来对算法进行加速:惯性技术本质上是通过考虑前两次的迭代信息以计算下一次迭代点,这样就避免了相邻两次迭代之间剧烈的变化。Attouch [12]提出了惯性近端ADMM算法,其外推步结合了Nestrov的梯度加速算法。王等[13]使用了一种加权Frobenius范数的惯性形式,并使用了部分对称思想对算法进行加速。徐等[14]在ADMM算法中将Bregman距离和惯性技术结合起来解决非凸非光滑问题。

在求解过程中,当目标函数较为复杂时,ADMM算法可能无法推导出子问题的显示表达式或求解速度很慢。面对这种情况,一种思想是通过将一个困难的问题近似成为一个简单问题,通过牺牲求解的准确性以提高求解速度。目前,非精确的惯性ADMM算法[15]-[17]多用于求解可分离的凸优化问题,针对非凸同时具有不可分离结构的相关研究还很少,受以上文献启发,本文针对问题(1.1),提出了一种并行式惯性近端非精确ADMM (PIPI-ADMM)算法,该算法在交替方向法的基础上,通过对每个子问题的增广拉格朗日函数进行部分一阶近似,使得近似信息可以并行发送到每个子问题的求解过程中。同时结合惯性技术与近端算法,分别在 x,y,z 子问题中加入近端项,在 x,y 子问题中引入惯性项来加速求解过程。此外,在简单合理的的假设下,当正则化的增广拉格朗日函数满足Kurdyka-Lojasiewicz性质时,文章证明了算法针对非凸非光滑问题的强收敛性。最后,文章总结了算法在鲁棒性主成分分析(Robust PCA)问题中的表现,验证了该算法的有效性。

本文的剩余部分组织结构如下:第2节给出了一些基础的定义以及与收敛性证明相关的预备知识;第3节提出了PIPI-ADMM算法,并证明了算法的收敛性;第4节是数值实验部分,通过分析算法在实际问题中的性能验证算法的有效性;最后在第5节得出了结论。

2. 预备知识

本节给出文中所用记号的意义以及相关概念与性质。

n 表示 n 维欧氏空间, · 表示欧式范数, x,y = x T y= i=1 n x i y i ,对于矩阵 A ,其像空间被定义为 ImA:=( Ax:x n ) 表示在像空间 ImA 上的欧式投影。文章使用 δ min( A T A) 表示矩阵 A T A 的最小特征值。设集合 S n ,点 x n ,将点 x 到集合 S 的距离记作 d( x,S )= inf yS yx . S= 时,定义 d( x,S )=+ 。对于函数 f ,其定义域记作 domf ,即 domf:={ x n :f( x )<+ } 。函数 f 称为正常函数当且仅当 domf f>

定义2.1 [18]若函数 f x 0 处满足 f( x 0 ) liminf x x 0 f( x ) ,则称函数 f x 0 处下半连续。若函数 f domf 内每点处均下半连续,则称 f 为下半连续函数。

定义2.2 [18]设函数 f: n { + } 正常下半连续,则

i) 函数 f xdomf 处的Fréchet次微分记为

^ f( x )={ x * n : liminf yx,yx f( y )f( x ) x * ,yx yx 0 },

xdomf 时,记 ^ f( x )=

ii) 函数在 x domf 处的极限次微分记为

f( x )={ x * n : x k x,f( x k )f( x ), x ^ k ^ f( x k ), x ^ k x * }.

[18]:设 f: n { + } g: n { + } 为正常下半连续函数,则下列结论成立

1) x n ,有 f ^ ( x )f( x ) ,其中 f ^ ( x ) 是闭凸的,而 f( x ) 是凸的。

2) 若 x k * f( x k ) lim k+ ( x k , x k * )=( x, x * ) ,则 x * f( x ) ,这表明 f( x ) 是闭的。

3) 设 x n f 的极小值点,则 0f( x ) 。若 0f( x ) ,称 x f 的稳定点, f 的稳定点集记作 critf.

4) 若 g: n 连续可微,则 ( f+g )=f( x )+g( x ),xdomf.

定义2.3 ( x * , y * , λ * ) 为问题(1.1)的增广拉格朗日函数的稳定点,即 0 L β ( ω * ) ,当且仅当

{ x h( x , y , z * )+ A T λ f( x ), y h( x , y , z * )+ B T λ =g( y ), z h( x , y , z * )+ C T λ * =0, A x +B y +C z * b=0.

定义2.4 [19] (Kurdyka-Lojasiewicz性质)设函数 f: n { + } 为正常下半连续函数,若存在 η>0, x * 的邻域 U 及连续凹函数 φ:[ 0,η ) + 满足如下条件:

1) φ( 0 )=0

2) φ ( 0,η ) 连续可微的且在0处连续;

3) t( 0,η ) φ ( t )>0

4) xU[ f( x * )<f( x )<f( x * )+η ] ,有如下Kurdyka-Lojasiewicz不等式成立:

φ ( f( x )f( x * ) )d( 0,f( x ) )1 ,

则称函数 f x domf 处具有Kurdyka-Lojasiewicz (KL)性质。

Kurdyka-Lojasiewicz (KL)性质主要适用于非凸非光滑优化问题。此类问题由于目标函数结构复杂,通常缺乏传统凸优化中的全局收敛性保证。KL性质通过对目标函数的局部增长特性进行限制,使得即便在非凸优化中也能获得一定的收敛性保证。对于特定结构的函数(例如半代数函数和半解析函数),KL性质在其定义域上表现出良好的适用性。在实际应用中,很多函数都满足KL性质。例如 L p 范数,Relu激活函数等都满足KL性质,这使得KL性质在非凸非光滑问题的收敛性证明中具有重要作用。

引理2.1 [20] A m×n 为非零矩阵, δ min( A T A) 为该矩阵得最小特征值,则

(2.1)

其中,表示在像空间上的欧氏投影。

引理2.2 若函数连续可微,且其梯度-Lipschitz连续的,则

(2.2)

证明

如果取,由不等式,此时有

引理2.3 [21] (一致KL性质)设是紧集,函数是正常下半连续函数,若函数上是常数且在上得每个点都满足KL性质,那么存在,使得对于任意以及

3. 算法与其收敛性分析

本节将给出针对问题(1.1)的并行式惯性近端非精确ADMM算法(PIPI-ADMM),并证明其收敛性。

3.1. 算法及假设

算法通过使用一阶近似的方法,非精确求解每一个子问题的值,进而达到对算法的加速效果。对于每个子问题,定义其增广拉格朗日函数在第k次迭代处的一种一阶近似为:

由此,算法的迭代框架如下:

(3.1)

(3.2)

(3.3)

算法1:并行式惯性近端非精确ADMM算法

a) 对于给定的使,其中,

b) 依次迭代

(3.4)

c) 终止准则:计算

c 2 =A x k+1 +B y k+1 +C z k+1 b ε 2 .

在下文中,为简便证明过程,记 { ω k } 为算法产生的迭代序列,且

u=( x,y,z ), u k =( x k , y k , z k ),ω=( x,y,z,λ ), ω k =( x k , y k , z k , λ k ), ω * =( x * , y * , z * , λ * ), ω ^ =( x,y,z,λ, x , y , z ), ω ^ k+1 =( x k+1 , y k+1 , z k+1 , λ k+1 , x k , y k , z k ) Δ x k+1 = x k+1 x k 2 ,Δ y k+1 = y k+1 y k 2 ,Δ z k+1 = z k+1 z k 2 ,Δ u k+1 = u k+1 u k 2 (3.5)

注:在算法设计中,我们采用了非精确求解思想和惯性思想来加速算法的迭代求解。非精确思想如(3.1)~(3.3)所示,在求解子问题过程中,不追求精确地求解每一个子问题的解,而是通过求解耦合项的一阶近似,来获得一个近似解,以此来加快求解速度。这种方法通常用于大型或复杂的优化问题中,能够在计算量较小的情况下实现快速求解。惯性近端思想被应用在了(3.4)的子问题中。其中,惯性近端项为:

τ x k 2 x x k 2 + ρ x k ( x, x k x k1 )  τ y k 2 y y k 2 + ρ y k ( y, y k y k1 ).

惯性技术通过在第 k+1 次迭代过程中考虑前两次迭代的信息 x k , y k , x k1 , y k1 ,使得算法可以加速向最优点逼近。特别是在面对非凸目标函数时,惯性技术可以更快速的跳过鞍点或局部最小值,在不增加计算量的前提下减少整体迭代次数,提升算法的性能。

下面,我们证明算法1的收敛性。首先,算法1的收敛性分析需要满足如下假设:

假设3.1

1) 函数 f: m + g: n + 是正常下半连续函数;

2) 函数 h: r × s × t 是连续可微函数且其梯度 h 是模 l h -Lipschiz连续,即

h( a )h( b ) ab ;

3) 矩阵 A,B,C 满足 C0 ,且 Im( C ){ b }Im( A )Im( B )

4) 算法参数选取满足: β> l h + ( τ z k ) 2 / δ min( C T C) min{ τ x , τ y , τ z } l h +2max{ ρ x , ρ y }+4

注:注意到,我们的算法的参数选择仅仅需要满足假设3.1(4)即可。该参数假设被使用在后文引理3.2的证明中,用以保证增广拉格朗日函数序列的单调性。在算法设计中, l h , ρ x 0 , ρ y 0 ,均为已知或已设置的常数,因此在参数选择中,我们只需令 τ x 0 , τ y 0 , τ z 0 β 足够大,即可保证算法是收敛的。

3.2. 收敛性证明

引理3.1 在算法1中,对于任意的 k ,下列不等式成立:

λ k+1 λ k 2 3 l h δ min( C T C) Δ x k+1 + 3 l h δ min( C T C) Δ y k+1 + 3( l h + ( τ z k ) 2 ) δ min( C T C) Δ z k+1 + 3 ( τ z k ) 2 δ min( C T C) Δ z k . (3.6)

证明 根据(3.4)中 z 子问题的最优条件得

0= z h( x k , y k , z k ) C T λ k +β C T ( A x k+1 +B y k+1 +C z k+1 b )+ τ z k ( z k+1 z k ). (3.7)

将(3.4)中第四式代入(3.7)中,得到

C T λ k+1 C T λ k 2 = z h( x k , y k , z k ) z h( x k1 , y k1 , z k1 )+ τ z k Δ z k+1 τ z k1 Δ y k 2 = z h( x k , y k , z k ) z h( x k1 , y k1 , z k1 ), τ z k1 Δ z k 2 + τ z k Δ z k+1 2 + τ z k1 Δ z k 2 2 τ z k Δ z k+1 , τ z k1 Δ z k 2 z h( x k , y k , z k ) z h( x k1 , y k1 , z k1 ), τ z k1 Δ z k

+2 z h( x k , y k , z k ) z h( x k1 , y k1 , z k1 ), τ z k Δ z k+1 3 l h 2 Δ x k 2 +3 l h 2 Δ y k 2 +3 l h 2 Δ z k 2 +3 τ z k 2 Δ z k+1 2 +3 τ z k 2 Δ z k 2

根据引理2.1可得

λ k+1 λ k 2 3 l h δ min( C T C) Δ x k+1 + 3 l h δ min( C T C) Δ y k+1 + 3( l h + ( τ z k ) 2 ) δ min( C T C) Δ z k+1 + 3 ( τ z k ) 2 δ min( C T C) Δ z k .

由此,引理得证。 □

为了证明收敛性,确定增广拉格朗日函数得单调性是非常必要的。但是由于算法添加了惯性项,直接证明增广拉格朗日函数单调非常困难。因此通过构建正则化的增广拉格朗日函数 L ^ β ( ω ^ k+1 ) ,借助 L ^ β ( ω ^ k+1 ) 的单调性来证明算法的收敛性。正则化的增广拉格朗日函数定义如下:

L ^ β ( ω ^ k+1 )= L ^ β ( u k+1 , λ k+1 , u k )= L β ( u k+1 , λ k+1 )+α( u k+1 u k 2 ). (3.8)

其中, ω ^ 的定义见式(3.5)。

引理3.2 对于算法1,选择足够大的 β, τ x 0 , τ y 0 , τ z 0 , 使得

β> l h + ( τ z k ) 2 / δ min( C T C) , min{ τ x , τ y , τ z } l h +2max{ ρ x , ρ y }+4, (3.9)

则对于 k ,有

L ^ β ( u k+1 , λ k+1 , u k )= L β ( u k+1 , λ k+1 )+ α 1 Δ u k+1 L β ( u k , λ k )+ α 2 Δ u k L β ( u k , λ k )+ α 1 Δ u k = L ^ β ( u k , λ k , u k1 ), (3.10)

其中, α 1 =min{ τ x k l h ρ x k 2 1, τ y k l h ρ y k 2 1, τ z k l h 2 1 }, α 2 =max{ ρ x k 2 , ρ y k 2 ,1 }.

证明 根据(3.4)第一式得:

f( x k+1 )+g( y k )+ x h( x k , y k , z k ), x k+1 x k λ k ,A x k+1 +B y k +C z k b + β 2 A x k+1 +B y k +C z k b 2 + τ x k 2 Δ x k+1 + ρ x k x k+1 , x k x k1 f( x k )+g( y k ) λ k ,A x k +B y k +C z k b + β 2 A x k +B y k +C z k b 2 + ρ x k x k , x k x k1 ,

根据(3.4)第二式得:

f( x k+1 )+g( y k+1 )+ y h( x k , y k , z k ), y k+1 y k λ k ,A x k+1 +B y k+1 +C z k b + β 2 A x k+1 +B y k+1 +C z k b 2 + τ y k 2 Δ y k+1 + ρ y k y k+1 , y k y k1 f( x k+1 )+g( y k ) λ k ,A x k+1 +B y k +C z k b + β 2 A x k+1 +B y k +C z k b 2 + ρ y k y k , y k y k1 ,

根据(3.4)第三式得:

f( x k+1 )+g( y k+1 )+ z h( x k , y k , z k ), z k+1 z k λ k ,A x k+1 +B y k+1 +C z k+1 b + β 2 A x k+1 +B y k+1 +C z k+1 b 2 + τ z k 2 Δ z k+1 f( x k+1 )+g( y k+1 ) λ k ,A x k+1 +B y k+1 +C z k b + β 2 A x k+1 +B y k+1 +C z k b 2 ,

将上述三个不等式相加,得到

f( x k+1 )+g( y k+1 ) λ k ,A x k+1 +B y k+1 +C z k+1 b + β 2 A x k+1 +B y k+1 +C z k+1 b 2 y h( x k , y k , z k ), y k+1 y k + x h( x k , y k , z k ), x k+1 x k + z h( x k , y k , z k ), z k+1 z k + τ x k 2 Δ x k+1 + τ y k 2 Δ y k+1 + τ x k 2 Δ x k+1 f( x k )+g( y k ) λ k ,A x k +B y k +C z k b + β 2 A x k +B y k +C z k b 2 + ρ x k x k x k+1 , x k x k1 + ρ y k y k y k+1 , y k y k1 .

根据三角不等式 a,b 1 2 a 2 + 1 2 b 2 得到

f( x k+1 )+g( y k+1 ) λ k ,A x k+1 +B y k+1 +C z k+1 b + β 2 A x k+1 +B y k+1 +C z k+1 b 2 y h( x k , y k , z k ), y k+1 y k + x h( x k , y k , z k ), x k+1 x k + z h( x k , y k , z k ), z k+1 z k + τ x k 2 Δ x k+1 + τ y k 2 Δ y k+1 + τ z k 2 Δ z k+1 f( x k )+g( y k ) λ k ,A x k +B y k +C z k b + β 2 A x k +B y k +C z k b 2 + ρ x k 2 Δ x k+1 + ρ x k 2 Δ x k + ρ y k 2 Δ y k+1 + ρ y k 2 Δ y k .

通过变换得到

f( x k )+g( y k )+h( x k , y k , z k ) λ k ,A x k +B y k +C z k b + β 2 A x k +B y k +C z k b 2 + ρ x k 2 Δ x k+1 + ρ x k 2 Δ x k + ρ y k 2 Δ y k+1 + ρ y k 2 Δ y k . f( x k+1 )+g( y k+1 )+h( x k+1 , y k+1 , z k+1 ) λ k ,A x k+1 +B y k+1 +C z k+1 b + β 2 A x k+1 +B y k+1 +C z k+1 b 2 ( h( x k+1 , y k+1 , z k+1 )h( x k , y k , z k ) ( x h( x k , y k , z k ) y h( x k , y k , z k ) z h( x k , y k , z k ) ) T ,( x k+1 x k y k+1 y k z k+1 z k ) )+ τ x k 2 Δ x k+1 + τ y k 2 Δ y k+1 + τ z k 2 Δ z k+1

f( x k+1 )+g( y k+1 )+h( x k+1 , y k+1 , z k+1 ) λ k ,A x k+1 +B y k+1 +C z k+1 b + β 2 A x k+1 +B y k+1 +C z k+1 b 2 l h 2 ( Δ x k+1 + Δ y k+1 + Δ z k+1 )+ τ x k 2 Δ x k+1 + τ y k 2 Δ y k+1 + τ z k 2 Δ z k+1 ,

其中,最后一个不等式是由引理2.2得到的。由此得到

L β ( x k+1 , y k+1 , z k+1 , λ k )+ τ x k l h ρ x k 2 Δ x k+1 + τ y k l h ρ y k 2 Δ y k+1 + τ z k l h 2 Δ z k+1 L β ( x k , y k , z k , λ k )+ ρ x k 2 Δ x k + ρ y k 2 Δ y k . (3.11)

由(1.2)式,

L β ( x k+1 , y k+1 , z k+1 , λ k+1 ) L β ( x k+1 , y k+1 , z k+1 , λ k ) = λ k λ k+1 ,A x k+1 +B y k+1 +C z k+1 = 1 β λ k λ k+1 2 l h β δ min( C T C) Δ x k+1 + l h β δ min( C T C) Δ y k+1 + l h + ( τ z k ) 2 β δ min( C T C) Δ z k+1 + ( τ z k ) 2 β δ min( C T C) Δ z k . (3.12)

将(3.11)式与(3.12)式相加,得到

L β ( x k+1 , y k+1 , z k+1 , λ k+1 )+( τ x k l h ρ x k 2 l h β δ min( C T C) ) Δ x k+1 +( τ y k l h ρ y k 2 l h β δ min( C T C) ) Δ y k+1 +( τ z k l h 2 l h + ( τ z k ) 2 β δ min( C T C) ) Δ z k+1 L β ( x k , y k , z k , λ k )+ ρ x k 2 Δ x k + ρ y k 2 Δ y k + ( τ z k ) 2 β δ min( C T C) Δ z k .

由于 3( l h 2 + ( τ z k ) 2 ) β δ min( C T C) < 3( l h 2 + ( τ z 0 ) 2 ) β δ min( C T C) <1 ,则

L β ( x k+1 , y k+1 , z k+1 , λ k+1 )+( τ x k l h ρ x k 2 1 ) Δ x k+1 +( τ y k l h ρ y k 2 1 ) Δ y k+1 +( τ z k l h 2 1 ) Δ z k+1 L β ( x k , y k , z k , λ k )+ ρ x k 2 Δ x k + ρ y k 2 Δ y k + Δ z k .

由于 min{ τ x 0 , τ y 0 , τ z 0 } l h +2max{ ρ x 0 , ρ y 0 }+4, 由算法参数定义, min{ τ x k , τ y k , τ z k } l h +2max{ ρ x k , ρ y k }+4, α 1 =min{ τ x 0 l h ρ x 0 2 1, τ y 0 l h ρ y 0 2 1, τ z 0 l h 2 1 }, α 2 =max{ ρ x 0 2 +1, ρ y 0 2 +1 } ,根据(3.9)式有

L β ( u k+1 , λ k+1 )+ α 2 Δ u k+1 2 +ν Δ u k+1 2 L β ( u k , λ k )+ α 2 Δ u k 2 .

由于 α 2 α 1 υ= α 1 α 2 ,则 ν>0 。在(3.8)式中,令 α= α 2 , 则(3.10)式成立,引理得证。 □

引理3.3 考虑算法1生成的序列 ω k ,存在常数 C>0 使得

d( 0, L β ( u k+1 , λ k+1 ) )C( Δ x k+1 + Δ y k+1 + Δ z k+1 + Δ x k + Δ y k + Δ z k ) (3.13)

证明 根据增广拉格朗日函数的定义,得到

{ x L β ( u k+1 , λ k+1 )=f( x k+1 )+ x h( u k+1 ) A T ( λ k+1 β r k+1 ), y L β ( u k+1 , λ k+1 )=g( x k+1 )+ y h( u k+1 ) B T ( λ k+1 β r k+1 ), z L β ( u k+1 , λ k+1 )= z h( u k+1 ) C T ( λ k+1 β r k+1 ), λ L β ( u k+1 , λ k+1 )= r k+1 .

根据算法迭代的最优条件(3.4),推出

{ 0 x f( x k+1 )+ x h( x k , y k , z k ) A T λ k β A T ( A x 1 k+1 +B y k +C z k b )+ τ x k ( x k+1 x k )+ ρ x k ( x k1 x k ), 0 y g( y k+1 )+ y h( x k , y k , z k ) B T λ k β B T ( A x 1 k+1 +B y k+1 +C z k b )+ τ y k ( y k+1 y k )+ ρ y k ( y k1 y k ), 0= z h( x k , y k , z k ) C T λ k β C T ( A x 1 k+1 +B y k+1 +C z k+1 b )+ τ z k ( z k+1 z k ), λ k+1 = λ k β( A x k+1 +B y k+1 +C z k+1 b ) (3.14)

将上述两个式子结合在一起,得到 ( ζ x k+1 , ζ y k+1 , ζ z k+1 , ζ λ k+1 ) L β ( u k+1 , λ k+1 ) ,其中,

{ ζ x k+1 = x g( u k+1 ) x g( u ) A T ( λ k λ k+1 )+β A T B( y k+1 y k )+β A T C( z k+1 z k ) + τ x k ( x k+1 x k )+ ρ x k ( x k1 x k ) ζ y k+1 = y g( u k+1 ) y g( u ) B T ( λ k λ k+1 )+β B T C( z k+1 z k )+ τ y k ( y k+1 y k )+ ρ y k ( y k1 y k ) ζ z k+1 = z g( u k+1 ) z g( u ) C T ( λ k λ k+1 ) ζ λ k+1 = 1 β ( λ k λ k+1 )

由于 h( ) Lipschitz连续,根据假设3.1,结合引理3.1,存在 C>0 使得 d( 0, L β ( u k+1 , λ k+1 ) ) C( Δ x k+1 + Δ y k+1 + Δ z k+1 + Δ x k + Δ y k + Δ z k )

引理3.4 根据假设3.1,如果 L β 是强制的,对于算法1生成的序列 { ω k } ,满足:

1) 序列 { ω k } 是有界的。

2) k=0 + ω k+1 ω k 2 <+

证明1) 根据引理3.2和算法的定义,得到

L β ( ω k )( u k+1 , λ k , u k )( u 0 , λ 0 , u 1 )= L β ( ω 0 ).

根据 L β 的强制性以及 inf ω L β ( ω )> ,序列 { ω k } 是有界的。

2) 因为 { ω k } 是有界的,则 { ω ^ k } 也是有界的。因此其至少有一个聚点。假设 ω ^ * { ω ^ k } 的一个聚点,则必然存在一个子列 { ω ^ k j } 使得 lim j ω ^ k j = ω ^ * . 由于 f,g 是下半连续的, h 是连续函数,则

lim j L ^ β ( ω ^ k j ) L ^ β ( ω ^ * ).

因此, { L ^ β ( ω ^ k j ) } 是有下界的。根据引理3.2, { L ^ β ( ω ^ k ) } 是单调递减的,因此函数列 { L ^ β ( ω ^ k ) } 收敛且 L ^ β ( ω ^ * ) L ^ β ( ω ^ k ) 。同时有

ν Δ u k+1 2 L ^ β ( u k+1 , λ k+1 , u k ) L ^ β ( u k , λ k , u k1 ),

将上述不等式从 k=0,,M 相加,得

k=0 M ν Δ u k+1 2 L β ( u 0 , λ 0 ) L ˜ β ( u M+1 , λ M+1 , u M ) L β ( u 0 , λ 0 )inf L ˜ β ( ω ¯ )

M ,有 k=0 M ν Δ u k+1 2 <+ ,结合引理3.1,即可得 k=0 ω k+1 ω k 2 <+ . 引理得证。 □

定理3.1 对于算法1生成的序列 { ω k } ,令 S, S ^ 分别表示序列 { ω k } { ω ^ k } 聚点的集合。如果假设3.1成立,则下述结论成立

1) S, S ^ 是非空紧集且 lim k d( ω k ,S )= lim k d( ω ^ k , S ^ )=0.

2) 如果 u * = u *' 并且 ω * S ,则

3) Scrit L β .

4)序列 { L ^ β ( ω ^ k ) } 收敛,且

证明:1) 根据 S, S ^ 的定义,(1)显然成立。

2) 根据引理3.4和 { ω k } { ω ^ k } 的定义,显然(2)成立。

3) 对于 ω * S ,可以找到 { ω k } 的一个子列 { ω ^ k j } ,这个子列收敛到 ω * 。根据引理3.4,得到这时有

A x * +B y * +C z =b.

因此 ( x * , y * , z * ) 是问题(1.1)的可行点。由于 x k+1 x 子问题的最小值,则

f( x k j +1 )+g( y k j )+ x h( x k j , y k j , z k j ), x k j +1 x k j λ k j ,A x k j +1 +B y k j +C z k j b + β 2 A x k j +1 +B y k j +C z k j b 2 + τ x k j 2 Δ x k j +1 + ρ x k j x k j +1 , x k j x k j 1 f( x * )+g( y k j )+ x h( x k j , y k j , z k j ), x * x k j λ k j ,A x * +B y k j +C z k j b + β 2 A x * +B y k j +C z k j b 2 + τ x k j 2 x * x k j 2 + ρ x k j x * , x k j x k j 1 ,

考虑到得到

limsup j f( x k j )f( x ).

根据的下半连续性,有 f( x )lim inf j f( x k j ). 因此

lim j f( x k j )=f( x ). (3.15)

同理,有

lim j g( y k j )=g( y ). (3.16)

由于集合 f,g 是闭的以及 h 的Lipschitz连续性,在式(3.14)取 k= k j ,得到

{ 0 x f( x * )+ x h( x * , y * , z * ) A T λ * , 0 y g( y * )+ y h( x * , y * , z * ) B T λ * , 0= z h( x * , y * , z * ) C T λ * , A x * +B y * +C z * b=0

根据定义2.1, ω * L β 的一个稳定点,即 ω * Scrit L β .

4) 由于 ω ^ k S ^ ,存在 { ω ^ k j } 使得根据式(3.15)和(3.16),以及 h 的连续性,得到

(3.17)

考虑到, { L ^ β ( ω ^ k ) } 是单调的,所以 { L ^ β ( ω ^ k ) } 收敛,因此

定理3.2 对于算法1生成的序列 { ω k } ,如果假设3.1成立且 lim j L ^ β ( u k j , λ k j , u k j 1 ) ( u * , λ * , u * ) 具有Kurdyka-Lojasiewicz性质,这时,有

(3.18)

并且序列 { ω k } 收敛到 L β ( ) 的稳定点。

证明 以下从两个方面来证明该定理

1) 基于式(3.17),如果存在 k ˜ >0 使得 L ^ β ( ω ^ k ˜ )= L ^ β ( ω ^ * ) ,对(3.10)式移项得

α Δ u k+1 L β ( ω ^ k ˜ ) L β ( ω ^ k ˜ +1 ) L β ( ω ^ k ˜ ) L β ( ω ^ ).

由此可得,对任意 k> k ˜ ,都有 u k+1 = u k 。进一步,根据式(3.4),可得对于任意 k> k ˜ +1, 都有 ω k+1 = ω k ,因此可以得到(3.18)。

2) 由于 { L ^ β ( ω ^ k ) } 单调递减,因此有 L ^ β ( ω ^ k+1 )> L ^ β ( ω ^ * ) 。从对任意常数 η ,存在常数 k 1 >0 ,当 k> k 1 时, L ^ β ( ω ^ k )< L ^ β ( ω ^ )+η. 根据定理3.1(2),对于任意 ε>0, 存在 k 2 >0 使得对任意 k> k 2 d( ω ^ k , S ^ )ε. 因此,对于任意 ε,η>0 ,当 k>k'=max{ k 1 , k 2 } ,得到

因为 S ^ 是非空紧集并且 { L ^ β ( ω ^ k ) } S ^ 。根据引理2.3,得到

φ ( L β ( ω ^ k ) L β ( ω ^ * ) )d( 0, L β ( ω ^ k ) )1,k>k'.

由于 φ ( L β ( ω ^ k ) L β ( ω ^ * ) )>0 ,则

1 φ ( L β ( ω ^ k ) L β ( ω ^ * ) ) <d( 0, L β ( ω ^ k ) ) (3.19)

由函数 φ 是凹函数,得到

φ( L β ( ω ^ k ) L β ( ω ^ * ) )φ( L β ( ω ^ k+1 ) L β ( ω ^ * ) ) φ ( L β ( ω ^ k ) L β ( ω ^ * ) )( L β ( ω ^ k+1 ) L β ( ω ^ k ) )

t k = Δ x k+1 + Δ y k+1 + Δ z k+1 根据式(3.13)和式(3.19)得到

L β ( ω ^ k+1 ) L β ( ω ^ k ) φ( L β ( ω ^ k ) L β ( ω ^ * ) )φ( L β ( ω ^ k+1 ) L β ( ω ^ * ) ) φ'( L β ( ω ^ k ) L β ( ω ^ * ) ) ( φ( L β ( ω ^ k ) L β ( ω ^ * ) )φ( L β ( ω ^ k+1 ) L β ( ω ^ * ) ) )d( 0, L β ( ω ^ k ) ) ( φ( L β ( ω ^ k ) L β ( ω ^ * ) )φ( L β ( ω ^ k+1 ) L β ( ω ^ * ) ) )C( t k+1 + t k ).

由(3.10)和上述不等式,推出

v Δ u k+1 2 φ( L β ( ω ^ k ) L β ( ω ^ * ) )φ( L β ( ω ^ k+1 ) L β ( ω ^ * ) )C( t k+1 + t k )

根据不等式 a+b+c 3 a 2 +3 b 2 +3 c 2 得到

t k+1 3 Δ u k+1 2 3C ν φ( L β ( ω ^ k ) L β ( ω ^ * ) )φ( L β ( ω ^ k+1 ) L β ( ω ^ * ) )( t k+1 + t k )

根据不等式 ab a+ 1 4 b 得到

t k+1 3C v φ( L β ( ω ^ k ) L β ( ω ^ ) )φ( L β ( ω ^ k+1 ) L β ( ω ^ ) )+ 1 4 ( t k+1 + t k )

将上述不等式由 k=k'+2, , M 相加,得到

k=k'+2 M t k+1 C ν φ( L β ( ω ^ k'+2 ) L β ( ω ^ * ) )φ( L β ( ω ^ M ) L β ( ω ^ * ) )+ 1 4 k=k'+2 M ( t k+1 + t k )

M+ ,得到

k=k'+2 + t k+1 2( c ν φ( L β ( ω ^ k +2 ) L β ( ω ^ * ) )+ t k +2 ).

由于 ν,C,>0 且为常数,得到

k=0 x k+1 x k <+, k=0 y k+1 y k <+, k=0 z k+1 z k <+.

结合引理3.1,得到

则(3.18)式得证,因此序列 { ω k } 是收敛的,根据定理3.1, { ω k } 收敛到 L β ( ) 的稳定点 ω * 。 □

4. 算法的应用

鲁棒主成分分析(Robust PCA)

Robust PCA问题[13]-[22]可以概括为

min X +α Y 1 + ω 2 ZO F 2 ,s.t.X+Y=Z, (4.1)

其中, O p×n 是观测矩阵, X,Y,Z p×n 是决策变量,定义 X = i=1 min(p,n) | σ i ( X ) | 1 2 为核范数, 1 范数为 Y 1 = i=1 n i=1 p | Y ij | , ω 是罚参数, α X Y 1 之间的权重参数,用 λ 表示拉格朗日乘子,则问题(4.1)的增广拉格朗日函数被定义为:

L β ( X,Y,Z,λ )= X +α Y 1 + ω 2 ZO F 2 λ,X+YZ + β 2 X+YZ F 2 , (4.2)

问题(4.1)可以被算法1通过如下迭代方式求解:

{ X k+1 argmin X β x k ( x, y k , z k , λ k )+ τ x k 2 X X k 1 τ x k ρ X k ( X k X k1 ) F 2 , Y k+1 argmin Y β y k ( x k+1 ,y, z k , λ k )+ τ y k 2 Y Y k 1 τ x k ρ Y k ( Y k Y k1 ) F 2 , Z k+1 = argmin Z β z k ( x k+1 , y k+1 ,z, λ k )+ τ z k 2 Z Z k F 2 , λ k+1 = λ k β( X k+1 + Y k+1 Z k+1 ),

其中, β ( ) L β ( ) 的一阶近似,其定义见公式(3.1)~(3.3)。上述方程组的显式表达式为

{ X k+1 =( β( Z k Y k )+ λ k + τ X k X k + ρ X k ( X k X k1 ) β+ τ X k , 1 β+ τ X k ), Y k+1 =S( β( Z k X k+1 )+ λ k + τ Y k Y k + ρ Y k ( Y k Y k1 ) β+ τ Y k , α β+ τ Y k ), Z k+1 = τ Z k Z k +β( X k+1 + Y k+1 )+ωO λ k β+ω+1 , λ k+1 = λ k β( X k+1 + Y k+1 Z k+1 ),

其中, ( ,μ ) 是奇异值阈值算子[23] S( ,μ ) 是软阈值算子[24]

在下述结果中, r.( X ),spr.( Y ) 分别表示矩阵 X 的秩和矩阵 Y 的稀疏率,实验选取了三组不同维度的观测矩阵 O 和两组不同的 ( r.( X ),spr.( Y ) ) 来进行矩阵分解。实验中,取 α= 0.5 p ,ω=1.0, 初始矩阵 X 0 , X 1 , Y 0 , Y 1 , Z 0 , Z 1 均被设置为零矩阵。算法1的参数设置为 β=200, ρ x 0 = ρ y 0 = ρ z 0 =1 τ x 0 = τ y 0 = τ z 0 =6. 实验的终止准则为

RelChg:= ( X k+1 , Y k+1 , Z k+1 ) ( X k , Y k , Z k ) F ( X k , Y k , Z k ) F +1 10 6 .

实验分别使用PIPI-ADMM和线性化ADMM算法(LADMM) [25]求解问题(4.1)。表1给出了两个算法求解该问题的数值结果,其中 Iters,Times( s ),( r.( X ^ ),spr.( Y ^ ) ) 分别表示迭代次数,求解时间以及分解后矩阵的秩和稀疏率。另外,使用

RelErr:= ( X ^ , Y ^ , Z ^ )( X * , Y * , Z * ) F ( X * , Y * , Z * ) F +1 ,

来衡量矩阵的数值解与原始矩阵的误差,其中, X * , Y * , Z * 表示原始矩阵, X ^ , Y ^ , Z ^ 表示问题(4.1)的数值解。从表1可以看出,PIPI-ADMM的迭代步数和计算时间明显小于LADMM算法,为更加直观展示结果,图1绘制了 p=100,n=200 式目标函数及其相关误差的变化趋势。上述实验结果验证了PIPI-ADMM算法的有效性,并且可以看出,PIPI-ADMM算法的数值效果明显优于LADMM算法。

Table.1. Numerical results of Robust PCA problem

1. Robust PCA问题的数值结果

p

n

( r.( X * ),spr.( Y * ) )

PIPI-ADMM

LADMM

Iters

Times (s)

( r.( X ^ ),spr.( Y ^ ) )

Iters

Times (s)

( r.( X ^ ),spr.( Y ^ ) )

100

200

(2, 0.05)

561

5.359

(2, 0.05)

624

5.901

(2, 0.05)

(2, 0.1)

858

7.933

(2, 0.1)

905

9.315

(2, 0.1)

(10, 0.05)

654

5.657

(10, 0.05)

686

5.849

(10, 0.05)

(10, 0.1)

923

11.246

(10, 0.1)

1271

21.045

(10, 0.1)

500

1000

(2, 0.05)

1743

72.432

(2, 0.05)

1865

75.847

(2, 0.06)

(2, 0.1)

2763

153.543

(2, 0.1)

3140

172.624

(2, 0.1)

(10, 0.05)

1822

74.469

(10, 0.05)

1970

87.570

(10, 0.05)

(10, 0.1)

3109

198.734

(10, 0.1)

3932

225.623

(10, 0.12)

Figure 1. The trends of RelChg,RelErr,r.( X k ), and spr.( Y k ) vary with the changes in number of Iters when p=100,n=200 and r.( X )=10,spr.( Y )=0.1

1. p=100,n=200,r.( X )=10,spr.( Y )=0.1 时两种算法的 RelChg,RelErr,r.( X k ) spr.( Y k ) 随迭代次数的变化情况

此外,我们还对PIPI-ADMM算法中惯性因子 ρ x 0 = ρ y 0 的大小对算法的加速效果进行了探究。图2绘制了不同的惯性因子对算法收敛速度的影响。可以看出,当惯性因子越大时,达到终止准则所需要的迭代步数越少。因此,惯性技术对于ADMM算法的加速是非常有效的,惯性因子越大,收敛速度越快,算法性能越强。

Figure 2. The trends of RelChg and RelErr vary with the changes in number of Iters when p=100 , and n=200 . r.( X * )=10,spr.( Y * )=0.1

2. p=100,n=200,r.( X )=10,spr.( Y )=0.1 时不同的惯性因子的 RelChg,RelErr ,随迭代次数的变化

5. 结论

本文针对具有不可分离结构的非凸非光滑问题,提出了并行式惯性近端非精确ADMM算法(PIPI-ADMM)。此算法结合了惯性技术与近端思想,通过非精确求解每个子问题的增广拉格朗日函数,以产生新的迭代点。同时,文章利用Kurdyka-Lojasiewicz性质证明了算法在非凸情形下的强收敛性。最后,数值实验结果说明,所提出的PIPI-ADMM算法具有更高的效率和更快的收敛速度。

NOTES

*通讯作者。

参考文献

[1] Liu, J., Ma, R., Zeng, X., Liu, W., Wang, M. and Chen, H. (2021) An Efficient Non-Convex Total Variation Approach for Image Deblurring and Denoising. Applied Mathematics and Computation, 397, Article 125977.
https://doi.org/10.1016/j.amc.2021.125977
[2] Zhou, S. and Li, G.Y. (2023) Federated Learning via Inexact ADMM. IEEE Transactions on Pattern Analysis and Machine Intelligence, 45, 9699-9708.
https://doi.org/10.1109/TPAMI.2023.3243080
[3] Han, D. (2022) A Survey on Some Recent Developments of Alternating Direction Method of Multipliers. Journal of the Operations Research Society of China, 10, 1-52.
https://doi.org/10.1007/s40305-021-00368-3
[4] Bai, J., Li, J., Xu, F. and Zhang, H. (2017) Generalized Symmetric ADMM for Separable Convex Optimization. Computational Optimization and Applications, 70, 129-170.
https://doi.org/10.1007/s10589-017-9971-0
[5] 薛中会, 殷倩雯, 党亚峥. 求解可分离凸优化问题的惯性近似松弛交替方向乘子法[J]. 上海理工大学学报, 2022, 44(2): 204-212.
[6] Guo, K., Han, D.R. and Wu, T.T. (2016) Convergence of Alternating Direction Method for Minimizing Sum of Two Nonconvex Functions with Linear Constraints. International Journal of Computer Mathematics, 94, 1653-1669.
https://doi.org/10.1080/00207160.2016.1227432
[7] Wu, Z., Li, M., Wang, D.Z.W. and Han, D. (2017) A Symmetric Alternating Direction Method of Multipliers for Separable Nonconvex Minimization Problems. Asia-Pacific Journal of Operational Research, 34, Article 1750030.
https://doi.org/10.1142/s0217595917500300
[8] Chao, M.T., Zhang, Y. and Jian, J.B. (2020) An Inertial Proximal Alternating Direction Method of Multipliers for Nonconvex Optimization. International Journal of Computer Mathematics, 98, 1199-1217.
https://doi.org/10.1080/00207160.2020.1812585
[9] Chen, C., He, B., Ye, Y. and Yuan, X. (2014) The Direct Extension of ADMM for Multi-Block Convex Minimization Problems Is Not Necessarily Convergent. Mathematical Programming, 155, 57-79.
https://doi.org/10.1007/s10107-014-0826-5
[10] Deng, W., Lai, M., Peng, Z. and Yin, W. (2016) Parallel Multi-Block ADMM with O(1/K) Convergence. Journal of Scientific Computing, 71, 712-736.
https://doi.org/10.1007/s10915-016-0318-2
[11] Wang, F., Cao, W. and Xu, Z. (2018) Convergence of Multi-Block Bregman ADMM for Nonconvex Composite Problems. Science China Information Sciences, 61, Article No. 122101.
https://doi.org/10.1007/s11432-017-9367-6
[12] Attouch, H. (2020) Fast Inertial Proximal ADMM Algorithms for Convex Structured Optimization with Linear Constraint. Minimax Theory and Its Applications, 6, 1-24.
[13] Wang, X., Shao, H., Liu, P. and Wu, T. (2023) An Inertial Proximal Partially Symmetric ADMM-Based Algorithm for Linearly Constrained Multi-Block Nonconvex Optimization Problems with Applications. Journal of Computational and Applied Mathematics, 420, Article 114821.
https://doi.org/10.1016/j.cam.2022.114821
[14] Xu, J. and Chao, M. (2021) An Inertial Bregman Generalized Alternating Direction Method of Multipliers for Nonconvex Optimization. Journal of Applied Mathematics and Computing, 68, 1-27.
https://doi.org/10.1007/s12190-021-01590-1
[15] Hager, W.W. and Zhang, H. (2020) Convergence Rates for an Inexact ADMM Applied to Separable Convex Optimization. Computational Optimization and Applications, 77, 729-754.
https://doi.org/10.1007/s10589-020-00221-y
[16] Liu, Q., Shen, X. and Gu, Y. (2019) Linearized ADMM for Nonconvex Nonsmooth Optimization with Convergence Analysis. IEEE Access, 7, 76131-76144.
https://doi.org/10.1109/access.2019.2914461
[17] Huang, Z., Hu, R., Guo, Y., Chan-Tin, E. and Gong, Y. (2020) DP-ADMM: ADMM-Based Distributed Learning with Differential Privacy. IEEE Transactions on Information Forensics and Security, 15, 1002-1012.
https://doi.org/10.1109/tifs.2019.2931068
[18] Rockafellar, R.T. and Wets, R.J.B. (2009) Variational Analysis. Springer Science & Business Media.
[19] Attouch, H., Bolte, J., Redont, P. and Soubeyran, A. (2010) Proximal Alternating Minimization and Projection Methods for Nonconvex Problems: An Approach Based on the Kurdyka-Łojasiewicz Inequality. Mathematics of Operations Research, 35, 438-457.
https://doi.org/10.1287/moor.1100.0449
[20] Goncalves, M.L.N., Melo, J.G. and Monteiro, R.D.C. (2017) Convergence Rate Bounds for a Proximal ADMM with Over-Relaxation Stepsize Parameter for Solving Nonconvex Linearly Constrained Problems. arXiv:1702.01850.
https://doi.org/10.48550/arXiv.1702.01850
[21] Bolte, J., Sabach, S. and Teboulle, M. (2013) Proximal Alternating Linearized Minimization for Nonconvex and Nonsmooth Problems. Mathematical Programming, 146, 459-494.
https://doi.org/10.1007/s10107-013-0701-9
[22] Xu, H., Caramanis, C. and Mannor, S. (2013) Outlier-Robust PCA: The High-Dimensional Case. IEEE Transactions on Information Theory, 59, 546-572.
https://doi.org/10.1109/tit.2012.2212415
[23] Cai, J., Candès, E.J. and Shen, Z. (2010) A Singular Value Thresholding Algorithm for Matrix Completion. SIAM Journal on Optimization, 20, 1956-1982.
https://doi.org/10.1137/080738970
[24] Bonesky, T. and Maass, P. (2009) Iterated Soft Shrinkage with Adaptive Operator Evaluations. Journal of Inverse and Ill-Posed Problems, 17, 337-358.
https://doi.org/10.1515/jiip.2009.023
[25] Hien, L.T.K., Phan, D.N. and Gillis, N. (2022) Inertial Alternating Direction Method of Multipliers for Non-Convex Non-Smooth Optimization. Computational Optimization and Applications, 83, 247-285.
https://doi.org/10.1007/s10589-022-00394-8