0 引 言

定义1(单调)^[1] 如果映射H:D⊆Rⁿ在D₀⊆D上满足

则称H在D₀上单调.

定义2(Slater约束规格)^[2] 对于凸优化问题

若存在x∈relint D，这里relint D是D的相对内部，使得

(1)

则称优化问题满足Slater约束规格.

考虑变分不等式问题VI(Ω,F)，即求x^*∈Ω，使

(2)

Ω={x∈Rⁿ₊|g(x)≤0}，这里F：Ω→Rⁿ连续单调，g(x)=(g₁(x),g₂(x),…,g₃(x))T，且g_i:Rⁿ₊→R是可微凸函数.

易证x^*是变分不等式问题VI(Ω,F)的解的充分必要条件是x^*是问题的最优解.假设这个优化问题满足Slater条件，由约束最优性条件参考文献^[2]，则必存在乘子向量λ^*，μ^*使得

(3)

在(3)中消去μ^*，问题归结为求互补问题

(4)

令ω=(x,λ)T,，问题(4)与变分不等式(5)等价.

(5)

本文将设计算法求解变分不等式(5)，文献^[3]在G(ω)单调并且零点存在时，通过添加一个临近正则项ω－ω^k，将求变分不等式(5)转化为求解方程组

(6)

其中ξ^k是一个误差序列，满足是适当参数序列，η_k∈(0,L)这里L>0，它的精确形式是ξ^k=0,k=1,2,….可以证明{ω^k}最终收敛于G(ω)的一个零点，这个算法被称作临近点算法.

当G(ω)在Rⁿ₊×R^m₊上零点不存在时，需要对临近点算法做一些改进，通常的做法是使用一个非线性项φ(ω,ω^k)替换ω－ω^k，文献^[4]中作者采用了正则项ω－ω^k+μ(ω－Ω²_kω^－1),Ωk= diag(ω^k₁,ω^k₂,…,ω^k_n+m),ω^－1=(ω^－1₁,ω^－1₂,ω^－1_n+m)T，它是函数f的梯度

这种方法被称为“对数二次临近点” (logrithmic-quadratic proximal,LQP)方法^[5].若ω^k+1是新产生的点，则必为内点.算法满足时才能收敛于变分不等式(5)的一个解，由于ω^k+1的值在ξ_k给定后才能计算出，因此这个条件在实际应用中难以控制.为克服这个困难，本文在“算法”部分给出上述算法的一种改进算法，并在“数值实验”部分给出验证结果.

1 算 法

文献^[6~8]中假设向量值函数G(ω)各分量为单调函数，据此给出基于LQP方法的交替方向法.本文只证明了G(ω)在Rⁿ₊×R^m₊中单调，虽不满足各分量均为单调函数，但可以参考它们对参数的处理方法，算法中用对角矩阵R替换η_k，由于R的对角元素可以不同，比用单一η_k更加灵活.另外我们在算法中使用了预测——校正技术^{[9, 10]}，证明了算法每一步都是可行的，并且只要初始点ω⁰>0，算法都将收敛于变分不等式(5)的解，下面是具体的算法.

1) 初始化：给定ε>0,μ∈(0,+∞),ω⁰=(x⁰,λ⁰)T>0，令k：=0；

2) 预估：求解下列方程组得到解

(7)

R,S为适当选取的正定对角矩阵，且满足

3) 收敛性检验：如果停止，否则进行下一步；

4) 校正：取θ>0，求解下列方程组得到解ω^k+1=((x^k+1)T,(λ^k+1)T)T

(8)

令k：=k+1，返回第一步.

引理1 假设R=diag(r₁,r₂,…,r_n)为正定对角矩阵，向量函数q(u)在Rⁿ₊上关于u单调，对于任意的μ>0,uk=(u^k₁,u^k₂,…,u^k_n)T>0，如果Uk=diag(u^k₁,u^k₂,…,u^k_n)，那么方程组

有且仅有惟一正解u.对于任意的v≥0

证 第一个结论在文献^[4]的推论2中有详细论证，这里略过证明.后面部分证明如下

引理2 G(ω)在Rⁿ×R^m₊上关于ω单调.

证 对于ω₁,ω₂∈Rⁿ×R^m₊

故G(ω)在Rⁿ×R^m₊上关于ω单调.

定理1 若ω^k>0，则预估方程组(7)存在惟一正解，对于这个，校正方程组(8)也存在惟一正解ω^k+1.

证 由引理 2知G(ω)在Rⁿ₊×R^m₊上单调，应用引理 1可知，方程组(7)存在惟一正解.同理，方程组(8)是(7)的特殊情况，所以也存在惟一正解.

这样我们就证明了算法的每一步都可以惟一确定，下面我们来证明算法的收敛性.

定理2 给定ω^k>0，是预估方程组(7)的解，ω^k+1是校正方程组(8)的解.那么对于变分不等式(5)的解ω^*有

这里.

证 对于方程组(7)令，应用引理 1得到

(9)

对于方程组(8)令，可得

(10)

由不等式组(10)可知

(11)

因为，由引理2得

(12)

又因为

(13)

由不等式(9)和等式(13)得到

(14)

这里〈x,y〉=x^TMy,x,y∈Rⁿ.结合不等式(12)和(14)由不等式(11)得：

令，则Ψ(θ)=－a_kθ²+2b_kθ+c_k是一个二次函数且在

(15)

时取得最大值，下面证明Ψ(θ^*)具有有限下界

定理3 给定ω^k>0，由预估方程组(7)产生，则

证 应用Cauchy-Schwarz不等式得

并且

(16)

所以

(17)

由不等式(16)和(17)得

(18)

定理4 给定ω^k>0，ω^*是变分不等式(5)的解，ω^k是预估方程组(7)的解，ω^k+1是校正方程组(8)的解，则均有界且

这里

证 由定理2和定理 3显然有

(19)

因为ρ>0，因此得到

(20)

所以序列{ω^k}有界.由不等式(19)得

因此

又因为{ω^k}有界，所以{也有界.

定理5 由算法产生的序列{ω^k}收敛于变分不等式问题(5)的一个解.

证 由定理1可知，，令应用引理1得

整理上式，并应用Cauchy-Schwarz不等式得

(21)

由定理 4知，是有限维空间中的有界序列，所以存在收敛子列，设

是它的一个收敛子列，那么ω^∞≥0，由不等式(21)可知

因此ω^∞是变分不等式(5)的一个解.又因为

所以对于ε>0，分别存在N₁,N₂∈N+使得当k>N₁和k_i>N₂时分别有

对上述ε，存在N=N₁+N₂，对于任意的k>N，令k>k_i>max{N₁,N₂}并结合不等式(20)得

故{ω^k}收敛于变分不等式问题(5)的解ω^∞.

2 数值结果

我们接下来通过数值实验来验证算法确实是可行的. 例1 Fx=2Ax－2b，其中A∈Rⁿ×n为一对称正定矩阵，b∈Rⁿ为一向量，考虑变分不等式问题.求x^*∈Ω使得

为了测试上述结论的正确性，取μ=s=0.25，θ=0.8，λ⁰=1.A=(aij)_n×n和R=(rij)_n×n为对角矩阵，b和x⁰为n维列向量，它们的元素值在{1,2}中随机取得，精度ε=1.0×10^－8.

预估步可以使用matlab中的solve和fsolve函数求解方程组(7)，fsolve使用数值方法求解方程组，每次只能求出一组与初始值有关的解.当这组解不满足迭代条件时，可以使用solve函数求方程组(7)的解析解，solve函数使用符号运算，可以求出方程组的所有解，但是求解速度比fsolve慢.在校正步中，由于已知，相应方程组可以转化为二次方程组，直接求出解的显式表示，

这里，表 1是我们的测试结果,通过分析表 1我们发现实验结果与理论是相符的，第一种算法多了求数值解的步骤，在问题不是很大时，可加快求解速度.但是当问题规模比较大时，求出的解往往是不满足迭代条件的，从而增加了求解时间.图 1是二维情形的散点图，图 2是图 1在收敛点附近的情况.这里A=diag(2,2)，R=diag(2,1)，b=(2,2)T，x⁰=(2,2)T.算法最终收敛于点(0.707 106 77，0.707 1067 7)，可以认为与理论上的点(是相符的.算例中当A为任意正定对称矩阵时，对角化后类似处理即可.文献^[11]中提到了解非线性方程组(7)的一种数值解法，由于已经证明，可以令x_i=l_iλ,l_i>0,i=1,2,…,n.方程组(7)中可以消去λ，得到关于l_i,i=1,2,…,n的n元方程组.依此方法继续，最终得到一元非线性方程，可运用二分法求相应方程，得到解后做相应替换，就可得到所需的非线性方程组(7)的解.在方程组(7)的解析解不存在的情况下可以考虑使用此方法，以上算例说明我们的算法在实际应用中是可行的.本文没有考虑R在迭代过程中变化的情况，也没有与其他算法做比较，这是我们接下来要做的工作.