考虑二阶锥线性互补问题:求向量 $x\in R^{n}$, 使得

$ \begin{equation} x\in K, Ax-b\in K, x^{T}(Ax-b)=0, \end{equation} $

(1.1)

其中 $R^{n}$是 $n$维欧式空间, $A\in R^{n\times n}$是 $n\times n$阶实矩阵, $b\in R^{n}$是 $n$维实向量, $K$是二阶锥的笛卡尔乘积.也就是说

$ \begin{equation} K=K^{n_{1}}\times K^{n_{2}}\times\cdot\cdot\cdot\times K^{n_{m}}, \end{equation} $

(1.2)

其中 $m, n_{1}, \cdot\cdot\cdot, n_{m}\geq1, n_{1}+\cdot\cdot\cdot+n_{m}=n$, 且 $K^{n_{i}}\subset R^{n_{i}}$是 $n_{i}$维二阶锥.即

$ \begin{equation} K^{n_{i}}=\{(x_{1}, x_{2})\in R\times R^{n_{i}-1}\mid x_{1} \geq \|x_{2}\|\}, \end{equation} $

(1.3)

其中 $\|\cdot\|$表示欧几里得范数.为了方便用 $(x_{1}, x_{2})$表示 $(x_{1}, x_{2}^{T})^{T}$, 用 $u\succeq_{K^{n_i}} v$ (或 $v\preceq_{K^{n_i}} u$)表示 $u-v\in K^{n_i}$.显然 $K\subset R^{n}$是一个闭凸且自对偶锥.若 $n_{i}=1$, $K^{1}$为非负实数集 $R_{+}$.若 $n_{1}=\cdot\cdot\cdot=n_{m}=1$, 则 $K^{n}$退化为 $R^{n}_{+}$, 此时问题(1.1) 退化为一般的线性互补问题.

二阶锥互补问题是近十多年来研究的新问题, 它是二阶锥规划的推广, 在工程设计、金融、管理科学等领域都有广泛应用^{[1, 2, 3]}.在过去的十多年中, 对二阶锥规划和二阶锥线性互补问题已有广泛研究, 并提出了各种算法, 例如光滑牛顿算法^{[4, 5, 6, 7]}、半光滑牛顿算法^{[8, 9]}、内点算法^[10]、价值函数算法^[11]、矩阵分裂算法^[12]等, 其中有些算法需要 $O(n^{3})$浮点运算次数, 这在问题维数 $n$变得很大时可能不是有效的.

众所周知, 罚函数方法是求解约束优化问题的重要方法, 其主要思路是将约束优化问题转化为一系列无约束优化问题来求解.文献[13]提出了具有良好性质的 $l_{1}$罚函数.之后, 许多学者致力于精确罚函数算法的研究^{[14, 15, 16, 17]}, 相比经典的 $l_{1}$精确罚函数, 这些方法只需较弱的条件即可保证精确性. 2008年, Wang和Yang^[18]提出一种幂罚函数算法求解线性互补问题.之后, Huang和Wang^{[19, 20]}将幂罚函数算法推广到求解非线性互补问题和混合互补问题. 2015年, Hao和Wan^[21]等将幂罚函数算法推广到求解二阶锥线性互补问题.实际上, 幂罚函数算法就是一类低阶罚函数算法.文献[18, 21]中罚函数的幂参数仅是形如 $1/k(k\in Z^{+})$的特殊数.基于此, 本文将文献[21]中幂参数范围推广到(0, 1]中任意实数, 提出一般低阶罚函数算法.同时证明在矩阵 $A$正定时, 低阶罚函数方程组的解序列在罚参数趋于无穷时以指数速度收敛到二阶锥线性互补问题的解, 并将低阶罚函数算法的数值结果与著名的光滑Fischer-Burmeister (F-B)函数算法进行比较, 结果表明提出的算法是有效的.

本文共分5节, 第2节给出一些预备知识; 第3节探讨二阶锥线性互补问题的低阶罚函数算法, 并证明算法的收敛性; 第4节给出一些数值实验结果; 第5节进行总结.在本文中, 对任意 $p>1$, $\|\cdot\|_{p}$表示 $l_{p}$范数.特别地, 当 $p=2$时, 它代表欧几里得范数 $\|\cdot\|$.

本节给出单块二阶锥 $K^{n}$的一些基础知识.对任意的 $x=(x_{1}, x_{2})\in R\times R^{n-1}$, $y=(y_{1}, y_{2})\in R\times R^{n-1}$, 它们的约当积(Jordan product)^[1]定义为

$ \begin{equation} x\circ y= (x^{T}y, x_{1}y_{2}+y_{1}x_{2}). \end{equation} $

(2.1)

本文用 $x^{2}$表示 $x\circ x$, 用 $x+y$表示普通向量加法, 这样 $“\circ "$, $“+ ”$以及 $e=(1, 0, \cdot\cdot\cdot, 0)\in R\times R^{n-1}$构成了 $K^{n}$的约当代数.易知约当积满足交换律, 但一般不满足结合律, 即 $(x\circ y)\circ z$与 $x\circ(y\circ\ z)$并不总相等, 这是研究和分析二阶锥优化算法的主要困难之一.但约当积在幂意义下结合律成立, 即 $(x\circ x)\circ x=x\circ(x\circ x)$对任意 $x\in R^{n}$成立.因此对任意正整数 $p$, $x^{p}$有意义, 且对任意正整数 $m$和 $n$, 有 $x^{m+n}=x^{m}\circ x^{n}$.关于与二阶锥相关的欧几里得约当代数的其他概念, 如迹、行列式、可逆向量、平方根、绝对值向量等, 详细内容可参考文献[1, 4, 22, 23].

与矩阵的谱分解类似, 以 $K^{n}$为基础, $R^{n}$中的任意向量都可以在二阶锥意义下进行谱分解^{[4, 22]}.对任意 $x=(x_{1}, x_{2})\in R\times R^{n-1}$, $x$的谱分解为

$ \begin{eqnarray} &&x=\lambda_{1}u^{(1)}+\lambda_{2}u^{(2)}, \\ \end{eqnarray} $

(2.2)

$ \begin{eqnarray} &&\lambda_{i}=x_{1}+(-1)^{i}\|x_{2}\|, u^{(i)}= \begin{cases} \frac{1}{2}\left(1, (-1)^{i}\frac{x_{2}}{\|x_{2}\|}\right), &\mbox{如果}\;x_{2}\neq 0, \\ \frac{1}{2}\left(1, (-1)^{i}\omega\right), &\mbox{如果}\;x_{2}=0, \end{cases}i=1, 2, \end{eqnarray} $

(2.3)

其中 $\lambda_{1}, \lambda_{2}$和 $u^{(1)}, u^{(2)}$分别表示 $x$的特征值和相应的特征向量, $\omega$是 $R^{n-1}$中满足 $\|\omega\|=1$的任意向量.易知 $e=u^{(1)}+u^{(2)}$, $\lambda_{1}\leq\lambda_{2}$, 且 $x_{2}\neq0$时分解式(2.2) 唯一.

下面给出谱分解的一些基本性质.

命题2.1^{[4, 22]}  对任意的 $x=(x_{1}, x_{2})\in R\times R^{n-1}$, $\lambda_{1}, \lambda_{2}$和 $u^{(1)}, u^{(2)}$是按(2.2) 和(2.3) 式给出的 $x$的特征值和相应的特征向量, 则

(1) $u^{(1)}\circ u^{(2)}=0$;

(2) $\|u^{(1)}\|=\|u^{(2)}\|=\frac{1}{\sqrt{2}}$;

(3) $u^{(i)}\circ u^{(i)}=u^{(i)}, i=1, 2$;

(4) $\lambda_{1}$是非负的(正的)当且仅当 $x\in K^{n}(x\in \text{int}K^{n})$, 其中 $\text{int}K^{n}$表示 $K^{n}$内部;

(5) $\text{tr}(x)=\lambda_{1}+\lambda_{2}=2x_{1}, \text{det}(x)=\lambda_{1}\lambda_{2}= x_{1}^{2}-\|x_{2}\|^{2}, 2\|x\|=\lambda_{1}^{2}+\lambda_{2}^{2}$.

由命题2.1, 对任意的 $x, y\in R^{n}$且 $x\in K^{n}, y\in K^{n}$, $\langle x, y\rangle=0$当且仅当 $x\circ y=0$.因此互补条件可以由约当积等价描述.根据 $x$的谱分解(2.2) 和(2.3), 一个标量函数 $\widehat{f}:R\rightarrow R$可以被扩展为与 $K^{n}(n\geq1)$相关的二阶锥向量值函数^{[4, 24]}

$ \begin{equation} f(x)=\widehat{f}(\lambda_{1})u^{(1)}+\widehat{f}(\lambda_{2})u^{(2)}, \end{equation} $

(2.4)

其中 $\lambda_{1}, \lambda_{2}$是 $x$的特征值, $u^{(1)}, u^{(2)}$是相应的特征向量.

式(2.4) 给出了定义二阶锥函数的一种形式, 但其含义还需进一步分析.对任意 $x\in R^{n}$, $x$到 $K^{n}$的投影定义为在 $K^{n}$上离 $x$最近的点, 记作 $[x]_{+}$, 即 $[x]_{+}\in K^{n}$且

$ \begin{equation} \parallel x-[x]_{+}\parallel=\text{min}\{\parallel x-y\parallel\;\mid y\in K^{n}\}. \end{equation} $

(2.5)

显然, 当 $n=1$时, 上述投影退化为

$ [t]_{+}=\text{max}\{0, t\} (t\in R). $

下述命题指出, $|x|$和 $[x]_{+}$具有(2.4) 式的形式(参见文献[4]命题3.3).

命题2.2  对任意 $x=(x_{1}, x_{2})\in R\times R^{n-1}$, $\lambda_{1}, \lambda_{2}$和 $u^{(1)}, u^{(2)}$是按式(2.2) 和(2.3) 给出的 $x$的特征值和相应的特征向量, 则

(1) $|x|=(x^{2})^{\frac{1}{2}}=|\lambda_{1}|u^{(1)}+|\lambda_{2}|u^{(2)}$,

(2) $[x]_{+}=(x+|x|)/2=[\lambda_{1}]_{+}u^{(1)}+[\lambda_{2}]_{+}u^{(2)}$.

类似于在 $K^{n}$上投影的概念, 定义^[6]

$ \begin{equation} [x]_{-}=[\lambda_{1}]_{-}u^{(1)}+[\lambda_{2}]_{-}u^{(2)}, \end{equation} $

(2.6)

显然 $[x]_{-}$是 $x$在 $-K^{n}$上投影的相反向量, 且 $x=[x]_{+}-[x]_{-}, [x]_{+}, [x]_{-}\in K^{n}$, $[x]_{+}\circ[x]_{-}=0$.

本节提出二阶锥线性互补问题(1.1) 的低阶罚函数算法, 并给出相应的收敛性分析.不失一般性, 假设向量由一个二阶锥块组成, 即 $K=K^{n}$, 因为收敛性分析很容易推广到二阶锥的一般形式(1.2).

下面考虑低阶罚函数方程组:求 $x_{\eta}\in R^{n}$, 使得

$ \begin{equation} Ax_{\eta}-\eta[-x_{\eta}]_{+}^{r}=b, \end{equation} $

(3.1)

其中 $0 < r\leq1$是幂参数, $\eta\geq1$是罚参数.当 $-x_{\eta}$的谱分解为 $-x_{\eta}=\gamma_{1}v^{(1)}+\gamma_{2}v^{(2)}$时, 定义

$ \begin{equation} [-x_{\eta}]_{+}^{r}=[\gamma_{1}]_{+}^{r}v^{(1)}+ [\gamma_{2}]_{+}^{r}v^{(2)}. \end{equation} $

(3.2)

当 $x\succeq_{K}0$不成立时, 方程组(3.1) 中的罚项 $\eta[-x_{\eta}]_{+}^{r}$惩罚 $x_{\eta}$的“负”项.由于 $\eta[-x_{\eta}]_{+}^{r}\in K$, 故 $Ax_{\eta}-b\succeq_{K}0$总成立.当 $\eta\rightarrow+\infty$时, 我们希望低阶罚函数方程组(3.1) 的解序列收敛于二阶锥线性互补问题(1.1) 的解.因此二阶锥线性互补问题(1.1) 被转化为低阶罚函数方程组序列.基于此, 给出如下低阶罚函数算法.

算法3.1  步骤1  取任意向量 $b\in R^{n}$, 矩阵 $A\in R^{n\times n}$, 令 $\widetilde{x}=0\in R^{n}$.

步骤2  若 $b\preceq_{K}0$, 停止, 转步骤7;否则, 转步骤3.

步骤3  若 $A^{-1}b\succeq_{K}0$, 令 $\widetilde{x}=A^{-1}b$, 停止, 转步骤7;否则, 转步骤4.

步骤4  设置幂参数 $0 < r\leq1$, 罚参数 $\eta\geq1$, 误差限为eps, 倍数参数 $c>1$, 选择一个初始点 $(x_{01}, x_{02})\in R\times R^{n-1}$, 且当 $x_{01} < 0$时, $x_{02}\neq0$.

步骤5  对罚参数 $\eta$和初始点 $x_{0}$, 解非线性方程组

$ \begin{equation} Ax-\eta[-x]_{+}^{r}=b, \end{equation} $

(3.3)

假设 $x_{\eta}$是(3.3) 式的解, 令 $\text{Tol}=\mid x_{\eta}^{T}(Ax_{\eta}-b)\mid$.

步骤6  如果 $\text{Tol}\leq\text{eps}$, 令 $\widetilde{x}=x_{\eta}$, 停止, 转步骤7;否则, 令 $x_{0}=x_{\eta}, \eta=c\eta$, 转步骤5.

步骤7  向量 $\widetilde{x}$是二阶锥线性互补问题(1.1) 的近似最优解.

下面讨论算法3.1的收敛性分析, 为此给出如下假设.

假设3.1  矩阵 $A$正定, 但不一定对称, 即存在正常数 $\omega$, 使得

$ \begin{equation} y^{T}Ay\geq\omega\|y\|^{2}, \forall y\in R^{n}. \end{equation} $

(3.4)

众所周知, 二阶锥线性互补问题(1.1) 与下述变分不等式等价:求 $x^*\succeq_{K}0$, 使得

$ \begin{equation} (y-x^*)^{T}Ax^*\geq(y-x^*)^{T}b, \;\forall y\succeq_{K}0. \end{equation} $

(3.5)

由假设3.1及文献[3]中定理2.3.3知变分不等式(3.5) 有唯一解, 因而二阶锥线性互补问题(1.1) 亦有唯一解.同时, 由于方程组可看作无约束变分不等式, 故(3.1) 式也有唯一解.

下述命题3.1给出低阶罚函数方程组(3.1) 解的有界性, 命题3.2给出 $\|[-x_{\eta}]_{+}\|$的一个上界, 其证明只需将文献[21]命题3.1和3.2中的幂参数 $1/k(k\in Z^{+})$换为式(3.1) 中的实数 $r\in(0, 1]$即可, 这里不再赘述.

命题3.1  对任意的 $\eta\geq1, 0 < r\leq1$, 低阶罚函数方程组(3.1) 的解有界, 即存在正常数 $W$, 独立于 $x_{\eta}, \eta$和 $r$, 使得 $\|x_{\eta}\|\leq W$.

命题3.2  对任意 $\eta\geq1, 0 < r\leq1$, 存在正常数 $M$, 独立于 $x_{\eta}$和 $\eta$, 使得

$ \begin{equation} \|[-x_{\eta}]_{+}\|\leq\frac{M}{\eta^{1/r}}. \end{equation} $

(3.6)

利用命题3.1和命题3.2, 可以得到如下重要结果, 它给出了算法3.1求解二阶锥线性互补问题(1.1) 的收敛性分析.

定理3.1  对任意 $\eta\geq1, 0 < r\leq1$, 令 $x_{\eta}$是低阶罚方程组(3.1) 的解, 且 $x^{*}$是二阶锥线性互补问题(1.1) 的解.则存在独立于 $x^{*}, x_{\eta}$和 $\eta$的正常数 $M$, 使得

$ \begin{equation} \|x^{*}-x_{\eta}\|\leq\frac{M}{\eta^{1/r}}. \end{equation} $

(3.7)

证  利用 $-x_{\eta}=[-x_{\eta}]_{+}-[-x_{\eta}]_{-}$, 向量 $x^{*}-x_{\eta}$可以分解为

$ \begin{equation} x^{*}-x_{\eta}=x^{*}-[-x_{\eta}]_{-}+[-x_{\eta}]_{+}=r_{\eta}+[-x_{\eta}]_{+}, \end{equation} $

(3.8)

其中 $[-x_{\eta}]_{-}$如(2.6) 式的形式所定义, 且

$ \begin{equation} r_{\eta}=x^{*}-[-x_{\eta}]_{-}. \end{equation} $

(3.9)

下面考虑 $r_{\eta}$的估计值.如果 $r_{\eta}=0$, 由(3.6) 和(3.8) 式可知不等式(3.7) 显然成立, 故下面只讨论 $r_{\eta}\neq0$的情形.显然 $[-x_{\eta}]_{-}\succeq_{K}0$, 由二阶锥线性互补问题与变分不等式(3.5) 的等价性及(3.9) 式, 在(3.5) 式中用 $[-x_{\eta}]_{-}$代替 $y$可得

$ \begin{equation} -r_{\eta}^{T}Ax^{*}\geq-r_{\eta}^{T}b. \end{equation} $

(3.10)

在(3.1) 式两端同时左乘 $r_{\eta}$, 可得

$ \begin{equation} r_{\eta}^{T}Ax_{\eta}-\eta r_{\eta}^{T}[-x_{\eta}]_{+}^{r}=r_{\eta}^{T}b. \end{equation} $

(3.11)

将(3.10) 和(3.11) 式相加得

$ \begin{equation} r_{\eta}^{T}A(x_{\eta}-x^{*})-\eta r_{\eta}^{T}[-x_{\eta}]_{+}^{r}\geq0. \end{equation} $

(3.12)

由(3.2) 式知 $[-x_{\eta}]_{+}^{r}\succeq_{K}0$, 又由于 $x^{*}\succeq_{K}0$, 根据(3.9) 式可得

$ \begin{equation} \begin{array}{rcl} r_{\eta}^{T}[-x_{\eta}]_{+}^{r} & = & (x^{*}-[-x_{\eta}]_{-})^{T}[-x_{\eta}]_{+}^{r}\\ & = & (x^{*})^{T}[-x_{\eta}]_{+}^{r}-[-x_{\eta}]_{-}^{T}[-x_{\eta}]_{+}^{r}\\ & = & (x^{*})^{T}[-x_{\eta}]_{+}^{r}\geq0. \end{array} \end{equation} $

(3.13)

将不等式(3.13) 应用于(3.12) 式得 $r_{\eta}^{T}A(x_{\eta}-x^{*})\geq0$, 即

$ \begin{equation} r_{\eta}^{T}A(x^{*}-x_\eta)\leq0. \end{equation} $

(3.14)

由(3.8) 式知 $r_{\eta}=x^{*}-x_{\eta}-[-x_{\eta}]_{+}$, 代入(3.14) 式得

$ \begin{equation} (x^*-x_{\eta})^{T}A(x^*-x_{\eta})\leq [-x_{\eta}]_+^{T}A(x^*-x_{\eta}). \end{equation} $

(3.15)

由假设3.1, 式(3.15), Cauchy-Schwarz不等式及范数相容性知, 存在常数 $b_{0}>0$, 使得

$ \begin{equation} \begin{array}{rcl} b_{0}\|x^*-x_{\eta}\|^2&\leq & (x^*-x_{\eta})^{T}A(x^*-x_{\eta})\\ &\leq & [-x_{\eta}]_+^{T}A(x^*-x_{\eta})\\ &\leq & \|[-x_{\eta}]_+\|\|A(x^*-x_{\eta})\|\\ &\leq & \|[-x_{\eta}]_+\|\|A\|\|x^*-x_{\eta}\|, % %b_{0}\|r_{\eta}\|^{2} & \leq & r_{\eta}^{T}Ar_{\eta} \leq r_{\eta}^{T}A[-x_{\eta}]_{+}\\ % & \leq & \|r_{\eta}\|\cdot\|A[-x_{\eta}]_{+}\|\\ % & \leq & \|r_{\eta}\|\cdot\|A\|\cdot\|[-x_{\eta}]_{+}\|, \end{array} \end{equation} $

(3.16)

其中 $\|A\|$表示矩阵 $A$的2-范数.用 $\|x^*-x_{\eta}\|$同除(3.16) 式两端得

$ \begin{equation} \|x^*-x_{\eta}\|\leq\frac{\|A\|}{b_{0}}\|[-x_{\eta}]_{+}\|. %\triangleq M^{'}\|[-x_{\eta}]_{+}\|. \end{equation} $

(3.17)

最后由(3.6) 及(3.17) 式知结论(3.7) 成立.证毕.

根据定理3.1, 在假设3.1下, 当 $\eta\rightarrow+\infty$时, 低阶罚函数方程组(3.1) 的解序列 $\{x_{\eta}\}$以指数速度收敛到二阶锥线性互补问题(1.1) 的解.但若假设3.1不成立, 又会有怎样的收敛性结果呢? 下述定理3.2证明了此种情形下低阶罚函数方程组(3.1) 解序列的任意极限点都是二阶锥线性互补问题(1.1) 的解.

定理3.2  对任意 $\eta_{i}$, 假设低阶罚函数方程组(3.1) 有解 $x_{\eta_{i}}$, 对任意幂参数 $0 < r\leq1$, 令 $\{\eta_{i}\}$是一列满足 $\eta_{i}\geq1$的正数序列, 且 $\lim\limits_{i\rightarrow{\infty}}\eta_{i}=+\infty$, 则序列 $\{x_{\eta_{i}}\}$的任意极限点都是二阶锥线性互补问题(1.1) 的解.

证   对任意 $\eta_{i}\geq1$, 因为 $x_{\eta_{i}}$是低阶罚函数方程组(3.1) 关于 $\eta_{i}$的解, 故

$ \begin{equation} Ax_{\eta_{i}}-b-\eta_{i}[-x_{\eta_{i}}]_{+}^{r}=0. \end{equation} $

(3.18)

设 $\overline{x}$是序列 $\{x_{\eta_{i}}\}$的任意极限点, 为了方便, 不妨假设当 $i\rightarrow\infty$时 $x_{\eta_{i}}\rightarrow\overline{x}$, 否则用 $\{x_{\eta_{i}}\}$的一个收敛到 $\overline{x}$子序列代替 $\{x_{\eta_{i}}\}$即可.因而存在正数 $P$, 使得

$ \begin{equation} \parallel x_{\eta_{i}}\parallel\leq P, i=1, 2, \cdot\cdot\cdot. \end{equation} $

(3.19)

注意到

$ \begin{equation} \lim\limits_{i\rightarrow{\infty}}x_{\eta_{i}}=\overline{x}, \end{equation} $

(3.20)

从式(3.18) 得 $Ax_{\eta_{i}}-b=\eta_{i}[-x_{\eta_{i}}]_{+}^{r}\in K$, 在此式中令 $i\rightarrow\infty$, 由 $K$的闭性知

$ \begin{equation} A\overline{x}-b=\lim\limits_{i\rightarrow{\infty}}Ax_{\eta_{i}}-b\in K. \end{equation} $

(3.21)

下面证明

$ \begin{equation} \overline{x}\in K. \end{equation} $

(3.22)

反设 $\overline{x}\notin K$, 此时令 $\varepsilon:=\parallel[-\overline{x}]_{+}^{r}\parallel$, 则 $\varepsilon>0$.故由(3.20) 式知当 $i$充分大时, $x_{\eta_{i}}\notin K$且 $\parallel[x_{\eta_{i}}]_{+}^{r}\parallel\geq\frac{1}{2}\varepsilon>0$.同时, 由(3.18) 式得当 $i\rightarrow\infty$时,

$ \parallel Ax_{\eta_{i}}-b\parallel=\parallel\eta_{i}[-x_{\eta_{i}}]_{+}^{r}\parallel \geq\frac{1}{2}\eta_{i}\varepsilon\rightarrow\infty, $

这与式(3.19) 矛盾, 故式(3.22) 成立.

下面再证明

$ \begin{equation} \overline{x}^{T}(A\overline{x}-b)=0. \end{equation} $

(3.23)

因为 $\overline{x}\in K$, 所以有以下三种情形.

情形1  如果 $\overline{x}=0$, 显然 $\overline{x}^{T}(A\overline{x}-b)=0$.

情形2  如果 $\overline{x}\in \text{int}K$, 则由(3.20) 式知当 $i$充分大时, $[-x_{\eta_{i}}]_{+}^{r}=0$.由(3.18) 式知

$ Ax_{\eta_{i}}-b=\eta_{i}[-x_{\eta_{i}}]_{+}^{r}, $

在此式中令 $i\rightarrow\infty$, 则 $A\overline{x}-b=0$, 所以 $\overline{x}^{T}(A\overline{x}-b)=0$.

情形3  如果 $\overline{x}\in \text{bd}K\backslash\{0\}$, 其中 $\text{bd}K$表示 $K$的边界, 此时令 $\overline{x}=(z_{1}, z_{2})\in R\times R^{n-1}$, 则 $z_{1}=\|z_{2}\|>0$, 根据向量的谱分解(2.2)-(2.3), $\overline{x}$可以分解为 $\overline{x}=\mu_{1}v^{(1)}+\mu_{2}v^{(2)}$, 这里 $\mu_{1}, \mu_{2}$和 $v^{(1)}, v^{(2)}$分别表示 $\overline{x}$的特征值和相应的特征向量.因此

$ \mu_{1}=z_{1}-\|z_{2}\|=0, \mu_{2}=z_{1}+\|z_{2}\|>0, v^{(1)}= (\frac{1}{2}, \frac{-z_{2}}{2\| z_{2}\|}), v^{(2)}=(\frac{1}{2}, \frac{z_{2}}{2\|z_{2}\|}), $

即 $\overline{x}=\mu_{2}v^{(2)}$.假设 $x_{\eta_{i}}$的谱分解为

$ x_{\eta_{i}}=\mu_{1}(i)v^{(1)}(i)+\mu_{2}(i)v^{(2)}(i), $

由于对任意 $y=(y_{1}, y_{2})\in R\times R^{n-1}$, 当 $\|y_{2}\|>0$时,

$ f_{1, 2}(y)=y_{1}\pm\| y_{2}\|, f_{3, 4}(y)=(\frac{1}{2}, \frac{\pm y_{2}}{2\|y_{2}\|}) $

都连续, 故由(3.20) 式知当 $i\rightarrow\infty$时,

$ \mu_{1}(i)\rightarrow \mu_{1}=0, \mu_{2}(i)\rightarrow \mu_{2}>0, v^{(1)}(i)\rightarrow v^{(1)}, v^{(2)}(i)\rightarrow v^{(2)}. $

特别地, 当 $i$充分大时有 $\mu_{2}(i)\geq\frac{1}{2}\mu_{2}>0$, 此时由(3.18) 式得

$ \begin{equation} \begin{array}{rcl} Ax_{\eta_{i}}-b&=&\eta_{i}[-x_{\eta_{i}}]_{+}^{r}\\ & = & \eta_{i}([-\mu_{1}(i)]_{+}^{r}v^{(1)}(i)+[-\mu_{2}(i)]_{+}^{r}v^{(2)}(i))\\ & = & \eta_{i}[-\mu_{1}(i)]_{+}^{r}v^{(1)}(i). \end{array} \end{equation} $

(3.24)

在(3.24) 式两边令 $i\rightarrow\infty$, 并记 $\lim\limits_{i\rightarrow{\infty}}(\eta_{i}[-\mu_{1}(i)]_{+}^{r})=c\geq0$, 则 $A\overline{x}-b=cv^{(1)}$.因此

$ \overline{x}^{T}(A\overline{x}-b)=\mu_{2}c(v^{(2)})^{T}v^{(1)}=0. $

由以上3种情形知式(3.23) 成立.

最后, 由(3.21)-(3.23) 式可得 $\overline{x}$是二阶锥线性互补问题(1.1) 的解.由极限点 $\overline{x}$的任意性, 可知结论成立.证毕.

本节用一些数值算例来检验算法3.1的有效性, 所有例子均采用离散牛顿法求解其对应于(3.1) 式的非线性方程组, 所有的数值实验运行环境为Intel(R) Core(TM) i5-3570K CPU 3.40GHz, 4G of RAM, and windows7 with MATLAB R2009a.

例4.1  考虑在 $K^{2}$上的二阶锥线性互补问题(1.1), 其数据如下

$ A= \left( \begin{array}{ccc} 1 & 1\\ 0 & 2\\ \end{array} \right), b= \left( \begin{array}{ccc} 0 \\ 4 \\ \end{array} \right). $

例4.1的精确解为 $x^{*}=(1, 1)^{T}$.下面考虑用算法3.1求解, 首先取幂参数 $r=1$, 初始点 $x_{0}=(-1, 1)^{T}$, 取罚参数 $\eta=10\times2^{i}, i=2, \cdot\cdot\cdot, 7$, 其数值实验结果如表 4.1所示, 其中NS代表数值解, Err代表数值解NS与 $x^{*}$的距离.事实上, 本例对应的低阶罚函数方程组(3.1) 可以利用分析方法求解, 且 $x_{\eta}=((\eta-2)/(1+\eta), (\eta+2)/(1+\eta))$, 容易得到

$ \parallel x^{*}-x_{\eta}\parallel=\sqrt{10}/(1+\eta)\leq\sqrt{10}/\eta. $

因此当 $\eta\rightarrow+\infty$时, $x_{\eta}$以指数速度 $O(\eta^{-1})$收敛到 $x^{*}$, 这与式(3.7) 的结论一致.由表 4.1可以看出数值实验结果和分析结果是一致的.

其次令 $\eta=10\times2^{i}, i=1, \cdot\cdot\cdot, 5, x_{0}=(-1, 1)^{T}$, 且 $r=\frac{3}{5}, \frac{2}{5}, \frac{\sqrt{2}}{5}$, 其数值解的误差结果如表 4.2所示.从表 4.2可以看出对于 $0<r\leq1$内一个固定的值, 当 $\eta$增大时, 误差 $\|x^{*}-x_{\eta}\|$趋于0.对固定的 $\eta>1$, 随着 $r$的减小, 误差 $\|x^{*}-x_{\eta}\|$逐渐趋于0.这些与式(3.7) 的结论一致.

下面考虑Chen和Ma^{[25, 26]}已经测试过的两个例子.在测试中令 $r=\frac{\sqrt{3}}{4}, c=10$且初始 $\eta=1000$.例4.2和例4.3的终止准则分别设为 $\text{eps}=10^{-8}$和 $\text{eps}=10^{-7}$.

例4.2  考虑在 $K^{5}$上的二阶锥线性互补问题(1.1), 其数据如下

$ A= \left( \begin{array}{ccccc} 15 & -5 & -1 & 4 & -5\\ 0 & 5 & 0 & 0 & 1\\ -1 & -3 & 8 & 2 & -3\\ 2 & -4 & 2 & 9 & -4\\ 0 & -5 & 0 & 0 & 10\\ \end{array} \right), b= \left( \begin{array}{c} 0 \\ 0 \\ 0 \\ 0 \\ 1 \\ \end{array} \right). $

本例的解 $x^{*}\approx(0.049185,-0.0030997,0.0096024,0.0031883,0.048033)^{T}$^[25], 且矩阵 $A$满足假设3.1.选取不同的初始点, 测试结果如表 4.3所示.

例4.3  考虑在 $K^{3}$上的二阶锥线性互补问题(1.1), 其数据如下

$ A= \left( \begin{array}{ccc} 21 & -9 & 18\\ -9 & 4 & -7\\ 18 & -7 & 19\\ \end{array} \right), b= \left( \begin{array}{c} -3 \\ -7 \\ -1 \\ \end{array} \right). $

本例的解 $x^{*}\approx(0.183606,-0.154346,-0.099440)^{T}$^[25], $A$是对称半正定矩阵.选取不同初始点, 测试结果如表 4.4所示.

从表 4.3和4.4可以看出, 例4.2和例4.3将文献[21]中形如 $1/k(k\in Z^{+})$的幂参数扩展为 $(0, 1]$中的一般实数.而且取幂参数 $r=\frac{\sqrt{3}}{4}$时, 其数值精确度比文献[21]中 $r=\frac{1}{2}$时的情形有明显改进.

众所周知, 当 $K=K^{n}$时, F-B函数定义为 $\phi_{FB}(x, y)=x+y-\sqrt{x^{2}+y^{2}}$, 满足

$ \phi_{FB}(x, y)=0\Leftrightarrow x\in K, y\in K, \langle x, y\rangle=0. $

二阶锥互补问题的光滑F-B函数 $\phi(\mu, x, y):R\times R^{n}\times R^{n}\rightarrow R^{n}$为^[4]

$ \begin{equation} \phi(\mu, x, y)=x+y-\sqrt{x^{2}+y^{2}+2\mu^{2}e}. \end{equation} $

(4.1)

上述平方和平方根都是指约当积意义下的平方和平方根.将 $y=Ax-b$代入(4.1) 式, 则二阶锥线性互补问题转化为光滑F-B方程, 即求 $x\in R^{n}$, 使得

$ \begin{equation} x+(Ax-b)-\sqrt{x^{2}+(Ax-b)^{2}+2\mu^{2}e}=0, \end{equation} $

(4.2)

其中 $\mu>0$是光滑参数.

众所周知, 当 $\mu\downarrow0$时, 方程组(4.2) 的解序列趋向于二阶锥线性互补问题(1.1) 的解^[4].下面给出相应的光滑F-B算法.

光滑F-B算法

步骤1  取任意向量 $b\in R^{n}$, 矩阵 $A\in R^{n\times n}$, 令 $\widetilde{x}=0$且 $\widetilde{x}\in R^{n}$.

步骤2  如果 $b\preceq_{K}0$, 停止, 转步骤7;否则, 转步骤3.

步骤3  如果 $A^{-1}b\succeq_{K}0$, 令 $\widetilde{x}=A^{-1}b$, 停止, 转步骤7;否则, 转步骤4.

步骤4  取光滑参数 $\mu>0$, 误差限为eps, 倍数参数 $0<d<1$, 初始点 $x_{0}$.

步骤5  对光滑参数 $\mu$和初始点 $x_{0}$, 解非线性方程组(4.2), 假设 $x_{\mu}$是(4.2) 式的解, 令 $\text{Tol}=|x_{\mu}^{T}(Ax_{\mu}-b)|$.

步骤6  如果 $\text{Tol}\leq \text{eps}$, 令 $\widetilde{x}=x_{\mu}$, 停止, 转步骤7;否则, 令 $x_{0}=x_{\mu}, \mu=d\mu$, 转步骤5.

步骤7  向量 $\widetilde{x}$是二阶锥线性互补问题(1.1) 的近似最优解.

例4.4用随机产生的二阶锥线性互补问题来比较算法3.1和光滑F-B算法的数值效果, 为了方便, 例4.4中随机产生的矩阵取对称正定矩阵.

例4.4  考虑 $K$为100个块的二阶锥线性互补问题, 即 $K=K^{n_{1}}\times K^{n_{2}}\times\cdot\cdot\cdot\times K^{n_{100}}$, 且 $A=\text{diag}(A_{1}, \cdots, A_{100}), b=(b_{1}, \cdots, b_{100})$.每一个块矩阵 $A_{i}(i=1, \cdot\cdot\cdot, 100)$是随机产生的对称正定矩阵.每一个向量块 $b_{i}(i=1, \cdot\cdot\cdot, 100)$满足 $q_{i}^{T}(Aq_{i}-b_{i})=0$, 其中 $q_{i}\in K^{n_{i}}$是随机产生的.

本例中 $A=\text{diag}(A_{1}, \cdots, A_{100})$是块对角矩阵, $q=(q_{1}, \cdots, q_{100})$是构造的精确最优解.首先固定 $r=\frac{\sqrt{2}}{5}$, 取 $n_{i}=2, 3, 4, 5(i=1, \cdots, 100)$, 在算法3.1中令 $c=10$, 初始 $\eta=10^{3}$, 而光滑F-B算法中令 $d=0.1$, 初始 $\mu=0.001$, 其所能达到的数值精度如表 4.5所示, 其中 $m$-Val代表 $\text{max}_{i}\{|\widetilde{x}_{i}^{T}(A\widetilde{x}_{i}-b_{i})|\}$, $a$-Val代表 $(\sum\limits_{i=1}^{100}|\widetilde{x}_{i}^{T}(A\widetilde{x}_{i}-b_{i})|)/100$, $m$-Err代表 $\text{max}_{i}\{\|\widetilde{x}_{i}-q_{i}\|\}$, $a$-Err代表 $(\sum\limits_{i=1}^{100}\|\widetilde{x}_{i}-q_{i}\|)/100$, $\widetilde{x}=(\widetilde{x}_{1}, \cdots, \widetilde{x}_{100})$表示数值最优解.

其次固定 $n_{i}=8$, 且 $r=\frac{\sqrt{3}}{2}, \frac{\sqrt{2}}{3}, \frac{3}{10}$, 在算法3.1中令 $c=10$, 初始 $\eta=10^{3}$, 而光滑F-B方法中令 $d=0.1$, 初始 $\mu=0.001$.取共同的终止准则为 $\text{eps}=10^{-6}$, 则数值结果如表 4.6所示, 其中m-Iter和a-Iter分别代表最大迭代次数和平均迭代次数.

由表 4.5可知两种算法对例4.4都是可行的, 但算法3.1的精度明显优于光滑F-B算法; 而表 4.6说明, 算法3.1与光滑F-B算法的结果相当, 但算法3.1所需时间明显少于光滑F-B算法.由表 4.6还可看出, 算法3.1对不同参数 $r$所得结果也不同, 但一般来说, 参数较小时数值结果更优, 这与结论(3.7) 一致.

本文将文献[21]中二阶锥线性互补问题的幂罚函数算法推广到一般的低阶罚函数算法.当矩阵 $A$正定时, 低阶罚函数方程组的解序列以指数速度收敛到二阶锥线性互补问题的解.当矩阵 $A$非正定时, 低阶罚函数方程组解序列的任意极限点是二阶锥线性互补问题的解, 但收敛速度未知.数值实验结果表明, 在矩阵正定条件下, 低阶罚函数算法的数值精度和计算时间优于光滑F-B算法, 而且在低阶罚函数算法中, 选取较小的幂参数 $r$, 数值结果更优.