线性圆锥互补问题的非单调非精确光滑牛顿法

汪 洋, 张所滨, 迟晓妮, 李 坤

(1. 桂林电子科技大学 数学与计算科学学院 广西密码学与信息安全重点实验室, 广西 桂林 541004;2. 桂林电子科技大学 计算机与信息安全学院, 广西 桂林 541004; 3. 桂林电子科技大学 数学与计算科学学院 广西高校数据分析与计算重点实验室, 广西 桂林 541004; 4. 桂林电子科技大学 数学与计算科学学院 广西自动检测技术与仪器重点实验室, 广西 桂林 541004)

1 引言及预备知识

考虑线性圆锥互补问题(LCCCP):寻找(x,y)∈Rn×Rn,使得

(1)

(2)

其中n=n1+n2+…+nN.ni-维圆锥定义[1]为

cosθi‖xi‖≤xi0},

(3)

Kni={xi=(xi0,xi1)∈R×Rni-1:

‖xi1‖≤xi0}.

(4)

圆锥与二阶锥之间的代数关系[2]:对任意xi=(xi0,xi1)∈R×Rni-1和yi=(yi0,yi1)∈R×Rni-1有

⟺Hixi∈Kni,

(5)

其中

圆锥互补问题(CCCP)在金融、图像处理和工程等领域有广泛的应用,如多指机器人的最优抓取操作、四足机器人动力优化等实际问题[3-5].由于圆锥是非对称锥,一些经典的算法一般不能直接应用到CCCP,因此CCCP的研究有重要的理论意义和实际应用价值.近年来,国内外学者在圆锥的性质、谱分解[2]和凸二次圆锥优化的原-对偶内点算法[6]等方面做了一些研究.但目前关于圆锥互补问题的算法尚不多见.本文基于文献[7]给出了一个圆锥互补函数的光滑函数,运用一个新的非单调线搜索技术[8],给出求解LCCCP的非单调光滑非精确牛顿方法.该算法采用一个凸组合型非单调技巧,且在每步迭代只需求解一个线性方程组和进行一次线搜索.在适当假设下证明算法具有全局收敛性和局部二阶收敛速度.数值结果表明算法的有效性.

下面给出与二阶锥相伴的欧几里得约当代数[9].对任意x=(x0,x1)∈R×Rn-1和y=(y0,y1)∈R×Rn-1,定义约当积为

x∘y=(xTy,x0y1+y0x1),

其单位元e=(1,0,…,0)∈Rn.对任意x=(x0,x1)∈R×Rn-1,定义箭形矩阵

其中I是n-1阶单位矩阵.对任意x,y∈Rn有L(x)y=x∘y.

下给出向量的谱分解[9].对任意x=(x0,x1)∈R×Rn-1,其谱分解为

x=λ1u(1)+λ2u(2),

(6)

其特征值和特征向量分别为

λi=x0+(-1)i‖x1‖,

其中,ω∈Rn-1是满足‖ω‖=1的任意向量,u(1)和u(2)分别是属于特征值λ1和λ2的特征向量.

定义1.1[10]如果对任意(x,y)∈Rn×Rn有

〈x-y,F(x)-F(y)〉≥0,

则称函数F:Rn→Rn是单调的.

定义1.2设函数G:Rn→Rm在x∈Rn是局部Lipschitz连续的.若G在x处方向可导,且对任意V∈∂G(x+Δx)有

G(x+Δx)-G(x)-V(Δx)=o(‖Δx‖),

其中∂G表示G的广义雅可比矩阵[11],则称G是半光滑的.

若G在x处半光滑,且

G(x+Δx)-G(x)-V(Δx)=O(‖Δx‖1+p),

其中0

若函数G:Rn→Rm在任意x∈Rn处是半光滑的(强光滑的),则它是半光滑(强光滑)函数.

2 光滑函数及其性质

对任意(x,y)∈Rn×Rn和μ∈R++,Chen-Harker-Kanzow-Smale(CHKS)光滑函数[12]

是二阶锥互补函数

的一个光滑函数.

本文考虑函数φ:R++×Rn×Rn→Rn有

φ(μ,x,y)=

(7)

当μ=0时

可得(7)式定义的φ(μ,x,y)是圆锥互补函数

的一个光滑函数.

令z=(μ,x,y)∈R++×Rn×Rn,定义函数Φ(μ,x,y):R1+2n→R1+2n为

(8)

基于圆锥互补函数的一个光滑函数(7)式,在每步迭代中应用光滑牛顿方法求解方程组Φ(z)=0.令μ↓0时,可得LCCCP(1)的解.记z*:=(μ*,x*,y*).显然Φ(z*)=0⟺(x*,y*)是LCCCP(4)的解.

定理2.1设Φ(z)由(8)式定义,则下列结论成立.

(i)Φ(z)在任意z=(μ,x,y)∈R++×Rn×Rn处连续可微,其雅可比矩阵为

Φ′(z)=

(9)

其中

H-L-1(w)L(w1)H,

(10)

H-1+L-1(w)L(w1)H-1,

(11)

w1:=Hx-H-1y,

(ii) 设M是半正定矩阵,则对任意z=(μ,x,y)∈R++×Rn×Rn,Φ′(z)可逆.

证明类似于文献[13]的定理2.4的证明,知(i)成立.

(ii) 令Δz:=(Δμ,Δx,Δy)∈R×Rn×Rn是Φ′(z)零空间中任一向量,只需证Δz=0即可.由(9)式知

Δμ=0,

(12)

-MΔx+Δy=0,

(13)

C(μ,x,y)Δμ+D(μ,x,y)Δx+

E(μ,x,y)Δy=0.

(14)

由(12)和(14)式得

D(μ,x,y)Δx+E(μ,x,y)Δy=0.

(15)

在(15)式的左右两边同时左乘L(w)得

L(w)D(μ,x,s)Δx+

L(w)E(μ,x,y)Δy=0.

(16)

由D(μ,x,y)和E(μ,x,y)的定义及H是对角阵有

L(w)D(μ,x,y)=L(w)H-L(w1)H=

L(w)E(μ,x,y)=L(w)H-1+L(w1)H-1=

又因为

≻Kn0.

由文献[14]的引理3.1得L(w)D(μ,x,y)和L(w)E(μ,x,y)都是正定矩阵必可逆,且

[L(w)D(μ,x,y)][L(w)E(μ,x,y)]

的对称部分正定.在(16)式两边同时左乘ΔxT[L(w)E(μ,x,y)]-1有

ΔxT[L(w)E(μ,x,y)]-1[L(w)D(μ,x,y)]×

Δx+ΔxTΔy=0.

由(13)式和M半正定得ΔxTΔy≥0,从而有

ΔxT[L(w)E(μ,x,y)]-1×

[L(w)D(μ,x,y)]Δx≤0.

(17)

又[L(w)D(μ,x,y)][L(w)(E(μ,x,y)]对称部分正定,从而

ΔxT[L(w)E(μ,x,y)]-1[L(w)D(μ,x,y)]Δx=

(18)

其中

故由(17)和(18)式得Δx=0.由于L(w)E(μ,x,y)可逆,故由(16)式得Δy=0.证毕.

3 LCCCP的非单调非精确光滑牛顿法

定义函数f:R++×Rn×Rn→R++为

f(z):=‖Φ(z)‖2=μ2+

‖y-(Mx+q)‖2+‖φ(μ,x,y)‖2.

(19)

算法3.1LCCCP的非单调非精确光滑牛顿法.

步1若‖Φ(zk)‖=0,停止,否则,令

ρ(zk):=γmin{1,f(z0),…,f(zk)}.

(20)

步2解方程组

(21)

其中,hk=(0,ϑk)∈R×R2n,‖hk‖≥εμ0min{1,f(z0),…,f(zk)}.求得搜索方向Δzk:=(Δμk,Δxk,Δyk)∈R×Rn×Rn.

步3确定步长λk=δlk.这里lk是满足下列不等式的最小非负整数

f(zk+δlΔzk)≤

Ck-2σ(1-γμ0-εμ0)μlCk.

(22)

步4令zk+1:=zk+λkΔzk,且

Ck+1=f(zk+1)+ηk(Ck-f(zk+1)).

(23)

置k:=k+1.转步1.

注3.1(i) 算法3.1运用了非单调线搜索技术[8];

(ii)Ck是Ck-1和f(zk)的凸组合;

(iii) 当ηk=0时为单调线搜索;

(iv) 参数ηk的选取对线搜索的非单调程度有影响.

定义集合

Γ={z=(μ,x,y)∈R++×Rn×Rn:

μ≥ρ(z)μ0},

其中ρ(z)由(20)式定义.

引理3.1设M是半正定矩阵且{zk=(μk,xk,yk)}是算法3.1生成的迭代序列,则对任意k≥0有

(i) {ρ(zk)}和{μk}都是单调递减序列;

(ii) 对任意k≥0有μk>0和zk∈Γ;

(iii) {Ck}是单调递减的,且对任意k≥0有

f(zk+1)≤Ck+1≤Ck.

证明由文献[15]的引理4.5易证(i)和(ii).下证(iii).对任意k≥0,由(22)式有

f(zk+1)-Ck≤Ck-2σ(1-γμ0-εμ0)=

δlkCk-Ck=

-2σ(1-γμ0-εμ0)δlkCk≤0.

(24)

由(23)式有

Ck+1-Ck=f(zk+1)+ηk(Ck-f(zk+1))-Ck=

(1-ηk)f(zk+1)+(ηk-1)Ck=

(1-ηk)(f(zk+1)-Ck)≤0.

(25)

又由(23)和(24)式得

f(zk+1)-Ck+1=ηk(f(zk+1)-Ck)≤0,

所以{Ck}单调递减,且f(zk+1)≤Ck+1≤Ck成立.

定理3.2设M是半正定矩阵,则算法3.1具有适定性.

证明由定理2.1知对任意μ>0有Φ′(zk)可逆,从而步2是适定的.下证步3是适定的.根据ρ(zk)的定义(20)式,当f(zk)<1时有

当f(zk)≥1时,有

因此对任意k≥0有

(26)

同理可得

‖Φ(zk)‖·‖hk‖≤εμ0f(zk).

(27)

对任意λ∈(0,1],定义

rk(λ):=f(zk+λΔzk)-

f(zk)-λf′(zk)Δzk.

(28)

由于f(·)在zk∈R++×Rn×Rn处连续可微,故

|rk(λ)|=o(λ).

(29)

由(19)、(21)、(26)～(29)式和引理3.1得

f(zk+λΔzk)=f(zk)+λf′(zk)Δzk+rk(λ)=

f(zk)+2λΦT(zk)Φ′(zk)Δzk+o(λ)=

2λΦT(zk)hk+o(λ)≤

f(zk)+2λγμ0f(zk)-2λf(zk)+

2λ‖Φ(zk)‖·‖hk‖+o(λ)≤

f(zk)-2λ(1-γμ0-εμ0)f(zk)+o(λ)≤

Ck-2λ(1-γμ0-εμ0)Ck+o(λ).

f(zk+λΔzk)≤Ck-2σ(1-γμ0-εμ0)λCk.

所以步3在第k次迭代是适定的,从而算法3.1适定.

4 收敛性分析

下面分析算法3.1的全局收敛性和局部二阶收敛性.

定理4.1设M是半正定矩阵,{zk=(μk,xk,yk)}是算法3.1生成的序列,则{zk}的任一聚点z*:=(μ*,x*,y*)都是Φ(z)=0的解,从而(x*,y*)是LCCCP(1)的解.

证明由引理3.1知{ρ(zk)}单调递减必收敛,因此存在常数ρ(z*)≥0,使得

设ρ(z*)>0,由引理3.1(ii)有

0<μ0ρ(z*)≤μ*.

故由引理3.1得

(30)

不失一般性,令

(31)

f(zk+λΔzk)≤Ck-2σ(1-γμ0-εμ0)λCk.

(32)

f(zk+1)=f(zk+λkΔzk)≤

Ck-2σ(1-γμ0-εμ0)λkCk≤

(33)

由(25)式得

Ck+1-Ck=(1-ηk)(f(zk+1)-Ck)≤0.

又由引理3.1知{Ck}单调递减且有下界,且1-ηk≥1-ηmax>0,故当k→∞时

从而

对(33)式两边取极限得

f(z*)≤C*-

(34)

定理4.2(局部二阶收敛) 设M是半正定矩阵,z*是算法3.1产生序列{zk}的任意聚点,且任意V∈∂Φ(z*)非奇异,则{zk}局部二阶收敛到{z*},即

‖zk+1-z*‖=O(‖zk-z*‖2),

证明类似于文献[16]中定理8的证明,略去.

5 数值算例

运用Matlab(R2012a)编程,在Intel(R) Core(TM) i5 6200U CPU @ 2.30GHz 2.40 GHz 8 GB内存,Windows 10 操作系统的计算机上做数值算例,测试算法3.1的数值效果.

算法3.1中参数取值为:μ0=0.1,δ=0.75,ηk=0.7,σ=0.225,γ=0.2,ε=0.1.当‖Φ(z)‖≤10-6时,算法3.1停止.Iter指平均迭代次数,ACPU指平均CPU时间.LCCCP(4)中矩阵M是通过M=NTN所得到的,N是方阵,N和q中元素是随机产生的位于区间[0,1]的数.令x0=e∈Rn,y0=0∈Rn为初始点.针对每种情形均随机产生10个问题进行求解,数值结果见表1～4.

表 1 运用算法3.1求解不同规模和旋转角的CCCP的数值结果

表 2 当时,算法3.1单调线搜索与非单调线搜索性能比较

表2～4给出运用算法3.1求解不同旋转角度和不同规模的LCCCP时单调线搜索和非单调线搜索所需的CPU时间和迭代次数.数值结果表明非单调线搜索要比单调线搜索所需的CPU时间和迭代次数要少.从表1～4的数值结果看出,运用算法3.1求解LCCCP所需的CPU时间和迭代次数较少,且相对稳定,从而表明算法3.1的有效性.

表 3 当时,算法3.1单调线搜索与非单调线搜索性能比较

表 4 当时,算法3.1单调线搜索与非单调线搜索性能比较