二阶锥规划基于SOC函数的新光滑牛顿算法

A New Smoothing Newton Method for Se-cond-Order Cone Optimization Based on SOC Function

摘要: 本文在向量值Fischer-Burmeister (FB)函数和向量值Natural-Residual (NR)函数的基础上，提出一种求解二阶锥规划(SOCP)问题的光滑函数。用一个带扰动的牛顿方程组去获得搜索方向，在适当假设下，分析了算法的全局收敛和局部收敛速度，给出了数值实验结果。

Abstract: Based on the Fischer-Burmeister function and the Natural-Residual function, a new smoothing Newton method is proposed for solving the second-order cone programming. This algorithm adopts a Newton equation with disturbance to gain the search direction. Under suitable assumptions, we prove that the proposed method is globally and locally quadratically convergent. Finally, some numerical results are given.

1. 引言

本文主要考虑以下线性二阶锥规划问题：

$\begin{array}{l}\min {\displaystyle \underset{i}{\overset{r}{\sum }}{c}_i^{\text{T}}{x}_i}\\ \text{s}\text{.t}\text{.}{\displaystyle \underset{i}{\overset{r}{\sum }}{A}_i{x}_i=b}，\\ x_i\underset{\_}{\succ }{}_{K^{n_i}}0,i\in I=\left\{1,2,\cdots ,r\right\}\end{array}$ (1.1)

其中， $b\in R^m,c_i\in R^{n_i},A_i\in R^{m\times n_i{}_{,}}$ 是已知量， $x_i\underset{\_}{\succ }{}_{K^{n_i}}0\left(i\in I\right)$ 表示 $x_i\in K^{n_i}$ 是变量。 $K^{n_i}$ 是维数为 $n_i$ 的二阶锥，即

$K^{n_i}=\left\{x_i=\left(x_{i0};{\bar{x}}_i\right)\in R^{n_i}|x_{i0}\ge \Vert {\bar{x}}_i\Vert \right\},$

这里 $\Vert .\Vert$ 为向量的欧几里得范数。易知 $K^{n_i}$ 是自对偶的，即 $K^{n_i}={\left(K^{n_i}\right)}^{*}=\left\{s_i\in R_i^n|x_i^{\text{T}}{s}_i\ge 0,\forall x_i\in K^{n_i}\right\}$ 其对偶规划为：

$\begin{array}{l}\max b^{\text{T}}y\\ \text{s}\text{.t}\text{.}{A}_i^{\text{T}}y+s_i=c_i,i\in I,\\ s_i\underset{\_}{\succ }{}_{K^{n_i}}0,i\in I.\end{array}$ (1.2)

其中 $y\in R^m,s_i\in K^{n_i},I=\left\{1,2,\cdots ,r\right\}$ 。为便于问题的描述，使用如下记号：

$\begin{array}{l}n\triangleq {\displaystyle \underset{i=1}{\overset{r}{\sum }}{n}_i},K=K^{n_i}\times K^{n_i}\times \cdots \times K^{n_i},A=\left(A_1,A_2,\cdots ,A_r\right)\in R^{m\times n},\\ c=\left(c_1,c_2,\cdots ,c_r\right)\in R^n,x=\left(x_1;x_2;\cdots ;x_r\right)\in K,s=\left(s_1;s_2;\cdots ;s_r\right)\in K.\end{array}$

这里用 $R^n$ 表示n维实列向量，且将 $R^{n_1}\times R^{n_2}\times \cdots \times R^{n_r}$ 定义为 $R^{n_1+n_2+\cdots +n_r}$，因此 $\left(x_1;x_2;\cdots ;x_r\right)\in K$ 是 $R^{n_1+n_2+\cdots +n_r}$ 中的列向量。则上面(1.1)和(1.2)式可转化为如下简洁形式：

$\min \left\{c^{\text{T}}x|Ax=b,x\in K\right\}.$ (1.3)

其对偶问题：

$\max \left\{b^{\text{T}}y|A^{\text{T}}y+s=c,s\in K,y\in R^m\right\}.$ (1.4)

二阶锥规划是一类非光滑凸规划问题，它是在仿射空间与有限个二阶锥的笛卡尔积的交集上极小化或极大化一个线性函数。在现实生活中有着广泛的应用：许多工程，力学，投资组合，线线阵列等都可以转化为二阶锥问题求解。此外，许多数学问题也可转化为二阶锥求解，如线性规划，二次规划，凸二次约束二次规划，范数极小化，因而研究二阶锥规划具有重要的理论与现实意义。本文在向量值FB函数和向量值NR函数的基础上进行改进，得到了一个新的光滑函数，利用得到的光滑函数，把二阶锥互补函数转化为一个非线性方程组问题，并给出求解二阶锥问题的一个光滑化算法，在无严格互补条件的假设下，证明算法全局收敛和局部二次收敛。

2. 基础知识

二阶锥规划建立在与二阶锥相伴的代数基础上，为了更好地研究二阶锥规划，下面将简要介绍与二阶锥相伴的欧几里得若当代数的基本概念和与之相关的重要结论，二阶锥规划的对偶理论和最优性条件等。先介绍一些记号，其中U表示反射矩阵， $E_{n-1}$ 表示 $\left(n-1\right)$ 阶单位阵， $Arw\left(x\right)$ 表示一个箭形矩阵，本文用相应的大写字母表示箭形矩阵， $X=Arw\left(x\right)$， $S=Arw\left(s\right)$， $W=Arw\left(w\right)$ 。

$U\triangleq \left(\begin{array}{cc}1& 0^{\text{T}}\\ 0& -E_{n-1}\end{array}\right)\in R^{n\times n},Arw\left(x\right)\triangleq \left(\begin{array}{cc}{x}_1& {\bar{x}}_{}^{\text{T}}\\ \bar{x}& x_0{E}_{n-1}\end{array}\right),e\triangleq \left(1;0\right)\in R^n.$

对于 $\forall x=\left(x_0;\bar{x}\right)\in R\times R^{n-1},s=\left(s_0;\bar{s}\right)\in R\times R^{n-1}$，定义如下若当积：

$x\circ s\triangleq \left(\begin{array}{c}{x}^{\text{T}}s\\ x_0\bar{s}+s_0\bar{x}\end{array}\right)=Arw\left(x\right)s=Arw\left(x\right)Arw\left(s\right)e=XSe.$

记 $x^2\triangleq x\circ x$， $x+y$ 表示通常的向量的向量加法。

定理2.1：(谱分解定理) (文 [1] )

对向量x定义与二阶锥有关的谱分解为： $x={\lambda }_1{c}_1+{\lambda }_2{c}_2$，其中x的谱值 ${\lambda }_i$ 和与之对应的谱向量 $c_i$ 如下：

$\begin{array}{l}{\lambda }_i=x_0+{\left(-1\right)}^{i+1}\Vert \bar{x}\Vert ,i=1,2.\\ c_i=\{\begin{array}{l}\frac{1}{2}\left(1;{\left(-1\right)}^{i+1}\frac{\bar{x}}{\Vert \bar{x}\Vert }\right),\bar{x}\ne 0;\\ \frac{1}{2}\left(1;{\left(-1\right)}^{i+1}\omega \right),\bar{x}=0;\end{array}i=1,2\end{array}$ 。 $\omega \in R^{n-1}$ 是满足 $\Vert \omega \Vert =1$ 的任意向量。

关于谱值，谱向量的相关性质归纳如下：

1) $x\in K\iff {\lambda }_2\ge {\lambda }_1\ge 0,x\in K^0\iff {\lambda }_2\ge {\lambda }_1>0$，

2) 对任意的 $x\in R^n$， $x^2={\lambda }_1^2{c}_1+{\lambda }_2^2{c}_2\in K$，

3) 对任意的 $x\in K$， $\sqrt{x}=\sqrt{{\lambda }_1}{c}_1+\sqrt{{\lambda }_2}{c}_2\in K$ 。

在线性规划中，互补性条件为：若 $x\ge 0,s\ge 0$，则 $x^{\text{T}}s=0$ 当且仅当 $x_i{s}_i=0,i=1,\cdots ,r$，二阶锥规划有着类似于线性规划的互补性条件。

引理2.1：(互补性条件)设 $x,s\in K$，即 $x_i,s_i\in K_i,i=1,\cdots ,r$，则 $x^{\text{T}}s=0$ 当且仅当 $x_i\circ s_i=0$ 。

在线性规划中，若原问题与对偶问题任意一个存在最优解，则另一个也存在最优解，且原问题与对偶问题最优值相等，但二阶锥规划不具有这种性质，若要二阶锥规划的强对偶理论成立，必须使得原问题与对偶问题同时存在严格可行解。因此二阶锥规划的最优性条件为：

定理2.2：(最优性条件) (文 [2] )如果(1.3)和(1.4)都有严格可行解，则 $\left(x,\left(y,s\right)\right)$ 是(1.3)和(1.4)的最优解对当且仅当

$\{\begin{array}{l}Ax=b,x\in K,\\ A^{\text{T}}y+s=c,s\in K,\\ x\circ s=0.\end{array}$ (2.1)

对二阶锥规划，其二阶锥互补函数定义如下：

定义2.1：(文 [3] )如果向量值函数 ${\varphi }^{\text{soc}}:R^n\times R^n\to R^n$ 满足

${\varphi }^{\text{soc}}\left(x,s\right)=0\iff x\circ s=0,x\in K,s\in K.$ (2.2)

则称 ${\varphi }^{\text{soc}}$ 为二阶锥互补函数。

结合问题(1.3)和(1.4)的最优性条件，定义函数 $H\left(x,y,s\right):R^{n+m+n}\to R^{m+n+n}$ 如下：

$H\left(x,y,s\right)=\left(\begin{array}{c}b-Ax\\ c-A^{\text{T}}y-s\\ {\varphi }^{\text{soc}}\left(x-s\right)\end{array}\right)，$ (2.3)

其中 ${\varphi }^{\text{soc}}\left(x-s\right)$ 是任意的二阶锥互补函数。易见 $H\left(x,y,s\right)=0$ 是与最优性条件等价的方程组，由此可知， $H\left(x,y,s\right)=0$ 的解 $\left(x^{*},y^{*},s^{*}\right)$ 满足最优性条件，是问题(1.3)和(1.4)的最优解对。由上易知，求解此方程组的关键在于构造合适的二阶锥互补函数。由于常见的二阶锥互补函数在(0,0) 点不可微(即非光滑)，例如向量值FB函数，向量值NR函数等，不能直接用牛顿法求解为了进一步研究，先给出光滑化的概念。

定义2.2：(文 [4] )对于不可微函数 $h:R^n\to R^m$，考虑带有参数 $\mu >0$ 的函数 $h_{\mu }:R^n\to R^m$，如果它具备如下性质：

1) 对于任意的 $\mu >0$， $h_{\mu }$ 是光滑的；

2) $\underset{\mu \to 0}{\lim }{h}_{\mu }\left(x\right)=h\left(x\right),\forall x\in R^n$，

则称 $h_{\mu }$ 是h的光滑函数。

3. 一个新的光滑函数及其性质

光滑函数在二阶锥规划光滑算法研究中起着重要作用(文 [5] [6] [7] [8] [9] )。由于向量值FB函数和向量值NR函数并不处处连续可微，大大影响了其实际应用。本文通过光滑化对称扰动得到一个新的向量值函数 $\varphi \left(\mu ,x,s\right):R_{+}\times R^n\times R^n\to R^n$ ：

$\varphi \left(\mu ,x,s\right)=x+s-\sqrt{\mu {\left(x-s\right)}^2+\left(1-\mu \right)\left(x^2+s^2\right)+2{\mu }^{2}e}$ (3.1)

为了减少变量的个数，使证明简单，令 $z=\left(\mu ,x,y\right)$，定义函数 $H\left(z\right)$ 如下：

$H\left(z\right)=\left(\begin{array}{c}{\text{e}}^{\mu }-1\\ b-Ax\\ \varphi \left(\mu ,x,c-A^{\text{T}}y\right)\end{array}\right)，$ (3.2)

显然，若 $\left(0,x^{*},y^{*}\right)$ 是 $H\left(z\right)=0$ 的解，则 $\left(x^{*},y^{*},c-A^{\text{T}}{y}^{*}\right)$ 是原问题(1.3)及其对偶问题(1.4)的最优解。因此，可考虑用牛顿法求解 $H\left(z\right)=0$，但对于所构造的 $H\left(z\right)$，其雅可比矩阵必须是非奇异的。在讨论 $H\left(z\right)$ 的性质之前先介绍相关结论。

引理3.1：(文 [3] )对于任意的 $a^1,\cdots ,a^p\in R^n$，令 $\chi \left(a^1,\cdots ,a^p\right)=\sqrt{{\displaystyle \underset{i=1}{\overset{p}{\sum }}{\left(a^i\right)}^2}}$ 。则 $\chi$ 是全局利普希茨连续的；若 $v_0\ne \Vert \bar{v}\Vert$，其中 $v=\left(v_0;\bar{v}\right)\in R\times R^{n-1}$ 且 $v={\displaystyle \underset{i=1}{\overset{p}{\sum }}{\left(a^i\right)}^2}$，则 $\chi$ 在 $\left(a^1,\cdots ,a^p\right)$ 的任意邻域是连续可微的。

定理3.1：设函数 $H\left(z\right)$ 由(3.2)定义，令 $\omega =\sqrt{\mu {\left(x-s\right)}^2+\left(1-\mu \right)\left(x^2+s^2\right)+2{\mu }^{2}e},\mu \in (0,1)$，则对任意的 $\varphi \left(\mu ,x,s\right)\in R_{+}\times R^n\times R^n$，有以下结论：

1) $\varphi \left(\mu ,x,s\right)$ 是全局利普希茨连续的处处光滑函数，且在任意 $\left(\mu ,x,s\right)$ 处连续可微，有

${{\varphi }^{\prime }}_{\mu }\left(\mu ,x,s\right)=-W^{-1}\left(x\circ s+2\mu e\right),$ (3.3)

${{\varphi }^{\prime }}_x\left(\mu ,x,s\right)=E_n-W^{-1}\left(X-\mu S\right),$ (3.4)

${{\varphi }^{\prime }}_s\left(\mu ,x,s\right)=E_n-W^{-1}\left(S-\mu X\right).$ (3.5)

2) $\varphi \left(\mu ,x,s\right)$ 在 $\mu =0$ 处是光滑的。

证明：1) 由引理3.1可知 $\varphi \left(\mu ,x,s\right)$ 是全局利普希茨连续且在任意 $\left(\mu ,x,s\right)$ 处连续可微。由 ${{\omega }^{\prime }}_{\mu }\left(\mu ,x,s\right)=W^{-1}\left(x\circ s\text{+}2\mu e\right)$，知

${{\varphi }^{\prime }}_{\mu }\left(\mu ,x,s\right)=-W^{-1}\left(x\circ s\text{+}2\mu e\right),$ (3.6)

同理可得：

$\begin{array}{l}{{\omega }^{\prime }}_x\left(\mu ,x,s\right)=W^{-1}\left(X-\mu S\right),\\ {{\omega }^{\prime }}_s=W^{-1}\left(S-\mu X\right).\end{array}$

由此可得：

${{\varphi }^{\prime }}_x\left(\mu ,x,s\right)=E_n-W^{-1}\left(X-\mu S\right),$ (3.7)

${{\varphi }^{\prime }}_s\left(\mu ,x,s\right)=E_n-W^{-1}\left(S-\mu X\right).$ (3.8)

2) 对于任意的 $x=\left(x_0;\bar{x}\right),s=\left(s_0;\bar{s}\right)\in R\times R^{n-1}$，由谱分解定理有

${\omega }^2={\lambda }_1\left(\mu \right)u_1\left(\mu \right)+{\lambda }_2\left(\mu \right)u_2\left(\mu \right),$

因此

$\varphi \left(\mu ,x,s\right)=x+s-\left(\sqrt{{\lambda }_1\left(\mu \right)}{u}_1\left(\mu \right)+\sqrt{{\lambda }_2\left(\mu \right)}{u}_2\left(\mu \right)\right),$

其中

${\lambda }_i\left(\mu \right)=\mu {\left(x-s\right)}^2+\left(1-\mu \right)\left(x^2+s^2\right)+2{\mu }^2+2{\left(-1\right)}^{i+1}\Vert v\left(\mu \right)\Vert ,i=1,2.$

$u_i\left(\mu \right)=\{\begin{array}{l}\frac{1}{2}\left(1;{\left(-1\right)}^{i+1}\frac{v\left(\mu \right)}{\Vert v\left(\mu \right)\Vert }\right),v\left(\mu \right)\ne 0;\\ \frac{1}{2}\left(1;{\left(-1\right)}^{i+1}w\right),v\left(\mu \right)=0.\end{array}i=1,2.$

其中 $v\left(\mu \right)=\mu \left(x_0-s_0\right)\left(\bar{x}-\bar{s}\right)+\left(1-\mu \right)\left(x_0\bar{x}+s_0\bar{s}\right)$， $\omega \in R^{n-1}$ 是满足 $\Vert \omega \Vert =1$ 的任意向量。同理有 $\varphi \left(0,x,s\right)=x+s-\left(\sqrt{{\lambda }_1}{u}_1+\sqrt{{\lambda }_2}{u}_2\right)$，其中 ${\lambda }_i={\Vert x\Vert }^2+{\Vert s\Vert }^2+2{\left(-1\right)}^{i+1}\Vert v\Vert ,i=1,2$ 。

$u_i=\{\begin{array}{l}\frac{1}{2}\left(1;{\left(-1\right)}^{i+1}\frac{v}{\Vert v\Vert }\right),v\ne 0;\\ \frac{1}{2}\left(1;{\left(-1\right)}^{i+1}w\right),v=0.\end{array}i=1,2.$

其中 $v=x_0{s}_0+\bar{x}\bar{s}$， $\omega \in R^{n-1}$ 是满足 $\Vert \omega \Vert =1$ 的任意向量。

不失一般性，取 $\omega \in R^{n-1}$ 和 $v\left(\mu \right)$ 中的一致，显然 $\underset{\mu \to 0}{\lim }v\left(\mu \right)=v$，于是 $\underset{\mu \to 0}{\lim }{\lambda }_i\left(\mu \right)={\lambda }_i,\underset{\mu \to 0}{\lim }{u}_i\left(\mu \right)=u_i$，因此 $\underset{\mu \to 0}{\lim }\varphi \left(\mu ,x,s\right)=\varphi \left(0,x,s\right)$，结合 $\varphi \left(\mu ,x,s\right)$ 在任意 $\left(\mu ,x,s\right)$ 处连续可微，知 $\varphi \left(\mu ,x,s\right)$ 是光滑的，因此由光滑函数的定义知 $\varphi \left(\mu ,x,s\right)$ 是 $\varphi \left(0,x,s\right)$ 的一个光滑函数。

4. 算法的描述

记 $z\left(\mu ,x,y\right):R_{+}\times R^n\times R^n,\theta \left(z\right)={\Vert H\left(z\right)\Vert }^2,\eta =\sqrt{\theta \left(z\right)}+1$，令 $\gamma \in \left(0,1\right)$，当 $H\left(z\right)\ne 0$ 时，定义函数 $\beta ,\bar{\beta }:R^{n+m+1}\to R_{+}:\beta \left(z\right)={\text{e}}^{\mu }\gamma \min \left\{1,\theta \left(z\right)\right\},\bar{\beta }\left(z\right)=\gamma \min \left\{1,\theta \left(z\right)\right\}$ 。在矩阵A行满秩的假设下，可证明算法是适定的。

算法4.1：(二阶锥规划的一个光滑型牛顿法)

步骤0：给出常数 $\delta ,\sigma ,{\mu }_0\in \left(0,1\right),\gamma \in \left(0,1\right)$ 满足 $\gamma {\mu }_0<\frac{1}{2},\gamma \eta <\frac{1}{2}$，任取 $\left(x^0,y^0\right)\in R^n\times R^m$， $z^0=\left({\mu }_0,x^0,y^0\right),\bar{z}=\left({\mu }_0,0,0\right)$，令 $k:=0$ 。

步骤1：若 $H\left(z^k\right)=0$，停止，其中 $\varphi \left(z^k\right)$ 及 $H\left(z^k\right)$ 分别由(2.5)和(2.6)式所定义，否则令 ${\beta }_k:=\beta \left(z^k\right)$ 。

步骤2：(确定搜索方向)解线性方程组 $H\left(z^k\right)+H^{\prime }\left(z^k\right)\Delta z^k={\beta }_k\bar{z}$ 得到搜索方向 $\Delta z^k$ 。

步骤3：(确定步长)记 ${\alpha }_k$ 是满足下式的最小非负整数 $\alpha$ ：

$\theta \left(z^k+{\delta }^{\alpha }\Delta z^k\right)\le \left[1-\sigma \left(1-2\gamma {\mu }_0\eta \right){\delta }^{\alpha }\right]\theta \left(z^k\right),$

令步长 $t_k={\delta }^{{\alpha }_k}$ 。

步骤4：令 $z^{k+1}=z^k+t_k\Delta z^k,k:=k+1$，返回步骤1。

在证明算法的适定性之前，先介绍相关性质。

引理4.1：(文 [10] ) 对任意 $\mu >0$，

$-\mu \le \frac{1-{\text{e}}^{\mu }}{{\text{e}}^{\mu }}\le -\mu {\text{e}}^{-\mu }.$ (4.1)

引理4.2：(文 [11] )对于任意的 $x,y\in R$ 和 $\omega \succ {}_{K^n}0$，如果 ${\omega }^2\succ {}_{K^n}x{}^2+y^2$，则

$\left[Arw\left(\omega \right)-Arw\left(x\right)\right]\left[Arw\left(\omega \right)-Arw\left(y\right)\right]\succ 0$， $Arw\left(\omega \right)-Arw\left(x\right)\succ 0$， $Arw\left(\omega \right)-Arw\left(y\right)\succ 0$，且当 $\succ$ 变为 $\underset{\_}{\succ }$ 时，上述结论仍然成立。

定理4.1：令 $z=\left(\mu ,x,y\right)\in R_{+}\times R^n\times R^m$， $H\left(z\right):R^{1+n+m}\to R^{1+m+n}$ 由(3.2)定义，则有以下结论成立：

1) 在任意点 $z=\left(\mu ,x,y\right)$ 处， $H\left(z\right)$ 全局利普希茨连续，并且连续可微，其雅克比矩阵表示为：

$H^{\prime }\left(z\right)=\left(\begin{array}{ccc}{\text{e}}^{\mu }& 0& 0\\ 0& -A& 0\\ {{\varphi }^{\prime }}_{\mu }\left(z\right)& {{\varphi }^{\prime }}_x\left(z\right)& -{{\varphi }^{\prime }}_s\left(z\right)A^{\text{T}}\end{array}\right),$ (4.2)

2) 对任意的 $\mu >0$， $H^{\prime }\left(z\right)$ 非奇异。

证明：由定理3.1知(1)成立。下证(2)成立，对任意给定的 $\mu >0$，要证 $H^{\prime }\left(z\right)$ 非奇异，只需证明方程组 $H^{\prime }\left(z\right)\Delta z=0$ 只有零解，即 $\Delta z=\left(\Delta \mu ,\Delta x,\Delta y\right)=0$ 。将(4.2)式代入 $H^{\prime }\left(z\right)\Delta z=0$ 中可知

$\Delta \mu =0,A\Delta x=0,{{\varphi }^{\prime }}_x\left(z\right)\Delta x-{{\varphi }^{\prime }}_s\left(z\right)A^{\text{T}}\Delta y=0$ (4.3)

将(3.7)，(3.8)式中的 ${{\varphi }^{\prime }}_x\left(z\right)$， ${{\varphi }^{\prime }}_x\left(z\right)$ 代入(4.3)式并将等式左右的两端同时左乘以W，得：

$\left[W-\left(X-\mu S\right)\Delta x\right]-\left[W-\left(S-\mu X\right)\right]A^{\text{T}}\Delta y\text{=}0.$ (4.4)

因 $w^2-\mu {\left(x-s\right)}^2-\left(1-\mu \right)\left(x^2+s^2\right)=2{\mu }^{2}e\in K^0$，由引理4.2可知：矩阵 $W-\left(S-\mu X\right)$ 正定，从而 $W-\left(S-\mu X\right)$ 可逆。将方程(4.4)左乘 $\Delta x^{\text{T}}{\left[W-\left(S-\mu X\right)\right]}^{-1}$，结合 $A\Delta x=0$，得： $\Delta x^{\text{T}}{\left[W-\left(S-\mu X\right)\right]}^{-1}\left[W-\left(X-\mu S\right)\right]\Delta x-\Delta x^{\text{T}}{\left[W-\left(S-\mu X\right)\right]}^{-1}\left[W-\left(S-\mu X\right)\right]A^{\text{T}}\Delta y=0.$ 令 $\Delta \tilde{x}{\left[W-\left(S-\mu X\right)\right]}^{-1}\Delta x，$ 则

$\Delta {\tilde{x}}^{\text{T}}\left[W-\left(X-\mu S\right)\right]\left[W-\left(S-\mu X\right)\right]\Delta \tilde{x}=0.$ (4.5)

由引理4.2可知矩阵 $\left[W-\left(X-\mu S\right)\right]\left[W-\left(S-\mu X\right)\right]$ 正定。因此 $\Delta \tilde{x}=0$，进而 $\Delta x=0$ 。结合矩阵A行满秩的假设和(4.4)式，知 $\Delta y=0$ 。这就说明 $H^{\prime }\left(z\right)\Delta z=0$ 只有零解，故 $H^{\prime }\left(z\right)$ 非奇异。

定理4.2：设矩阵A行满秩，如果 ${\mu }_k>0$，则算法4.1是适定的。

证明：因为矩阵A行满秩且 ${\mu }_k>0$，由定理4.1知 $H^{\prime }\left(z\right)$ 非奇异，所以算法4.1的步骤2是适定的。令 $\Delta z^k=\left(\Delta {\mu }_k,\Delta x^k,\Delta y^k\right)\in R\times R^n\times R^m$ 是步骤2中方程组的解，则对任意的 $t\in \left(0,1\right]$，有

${\mu }_{k+1}={\mu }_k+t\Delta {\mu }_k={\mu }_k+t\left(\frac{1-{\text{e}}^{{\mu }_k}}{{\text{e}}^{{\mu }_k}}+\frac{{\beta }_k{\mu }_0}{{\text{e}}^{{\mu }_k}}\right)\ge \left(1-t\right){\mu }_k+t{\mu }_0\frac{{\beta }_k}{e^{{\mu }_k}}>0.$ (4.6)

由 $\beta \left(z\right)$ 的定义有 $\beta \left(z\right)\le {\text{e}}^{\mu }\gamma \sqrt{\theta \left(z\right)}$，

$\begin{array}{c}{\text{e}}^{{\mu }_k+t\Delta {\mu }_k}-1={\text{e}}^{{\mu }_k}{\text{e}}^{t\Delta {\mu }_k}-1\\ ={\text{e}}^{{\mu }_k}\left(1+t\Delta {\mu }_k+O\left(t^2\right)\right)-1\\ ={\text{e}}^{{\mu }_k}-1+t\left(1-{\text{e}}^{{\mu }_k}\right)+t{\beta }_k{\mu }_0+O\left(t^2\right)\\ =\left(1-t\right)\left({\text{e}}^{{\mu }_k}-1\right)+t{\beta }_k{\mu }_0+O(t2)\end{array}$

从而

$\begin{array}{c}{\left({\text{e}}^{{\mu }_k+t\Delta {\mu }_k}-1\right)}^2={\left(1-t\right)}^2{\left({\text{e}}^{{\mu }_k}-1\right)}^2+2t\left(1-t\right){\beta }_k\left({\text{e}}^{{\mu }_k}-1\right){\mu }_0+t^2{\beta }_k{\mu }_0+O\left(t^2\right)\\ \le \left(1-2t\right)){\left({\text{e}}^{{\mu }_k}-1\right)}^2+2t{\text{e}}^{{\mu }_k}\gamma \sqrt{\theta \left(z^k\right)}\left({\text{e}}^{{\mu }_k}-1\right){\mu }_0+O\left(t^2\right)\\ \le \left(1-2t\right){\left({\text{e}}^{{\mu }_k}-1\right)}^2+2t\gamma {\mu }_0\eta \theta \left(z^k\right).\end{array}$

在上述最后一个不等式利用 ${\text{e}}^{{\mu }_k}-1\le \sqrt{\theta \left(z^k\right)},{\text{e}}^{{\mu }_k}\le \eta$ 可得。

令 $\phi \left(z^k\right)=\left(\begin{array}{l}b-A{x}^k\\ \varphi \left(z^k\right)\end{array}\right),$ 并记 $\Phi \left(z^k\right)={\Vert \varphi \left(z^k\right)\Vert }^2,$ 可得 ${\left({\Phi }^{\prime }\left(z^k\right)\right)}^{\text{T}}\Delta z^k=-2{\Vert \varphi \left(z^k\right)\Vert }^2=-2\Phi \left(z^k\right)$，

从而

$\begin{array}{c}\Phi \left(z^k+t\Delta z^k\right)=\Phi \left(z^k\right)+t{\left({\Phi }^{\prime }\left(z^k\right)\right)}^{\text{T}}\Delta z^k+o\left(t\right)\\ =\Phi \left(z^k\right)-2t\Phi \left(z^k\right)+o\left(t\right)\\ =\left(1-2t\right)\Phi \left(z^k\right)+o\left(t\right).\end{array}$

因此

$\begin{array}{c}\theta \left(z^k+t\Delta z^k\right)={\Vert H\left(z^k+t\Delta z^k\right)\Vert }^2\\ \le \left(1-2t\right){\left({\text{e}}^{{\mu }_k}-1\right)}^2+2t\gamma {\mu }_0\eta \theta \left(z^k\right)+\left(1-2t\right)\Phi \left(z^k\right)+o\left(t\right)\\ \le \left(1-2t\right)\theta \left(z^k\right)+2t\gamma {\mu }_0\eta \theta \left(z^k\right)+o\left(t\right)\\ \le \left[1-\left(1-2\gamma {\mu }_0\eta \right)t\right]\theta \left(z^k\right)+o\left(t\right).\end{array}$

因为 $\gamma {\mu }_0\eta <1$，所以存在一个常数 $\bar{t}\in \left(0,1\right]$，使得对任意的 $t\in \left(0,\bar{t}\right]$，有

$\theta \left(z^k+t\Delta z^k\right)\le \left[1-\sigma \left(1-2\gamma {\mu }_0\eta \right)t\right]\theta \left(z^k\right),$ (4.7)

即步骤3可在有限步终止，故步骤3是适定的。综上可知，算法4.1是适定的。

5. 收敛性分析

定理5.1：1) 设矩阵A行满秩且 $\left\{z^k:=\left({\mu }_k,x^k,y^k\right)\right\}$ 是算法4.1生成的无穷迭代点列，则对 $\forall k\ge 0$，若 ${\mu }_k>0,z^k\in \Lambda ,\theta \left(z^k\right)\le \theta \left(z^0\right)$，则 $z^{k+1}\in \Lambda$，其中

$\Lambda :=\left\{z:=\left(\mu ,x,y\right)\in R_{+}\times R^n\times R^m|\mu \ge \beta \left(z\right){\mu }_0\right\},$ (5.1)

2) (全局收敛性) $\left\{z^k\right\}$ 的任意聚点 $z^{*}=\left({\mu }^{*},x^{*},y^{*}\right)$ 都是 $H\left(z^k\right)=0$ 的解。

证明：1) 根据 $\beta \left(z\right)$ 的定义分两种情况讨论：

i) 若 $\theta \left(z^k\right)\le 1$，则 ${\beta }_k=\gamma {\text{e}}^{{\mu }_k}\Vert \theta \left(z^k\right)\Vert$ 。对任意的 $t\in \left(0,1\right]$，从(4.7)式知

$\theta \left(z^k+t\Delta z^k\right)\le \left[1-\sigma \left(1-2\gamma {\mu }_0\eta \right)t\right]\theta \left(z^k\right)$，故 $\bar{\beta }\left(z^k+t\Delta z_k\right)=\gamma \theta \left(z^k+t\Delta z_k\right).$ 。由于 $z^k\in \Lambda$，可得

$\begin{array}{l}{\mu }_k+t\Delta {\mu }_k-\bar{\beta }\left(z^k+t\Delta z^k\right){\mu }_0\\ \ge \left(1-t\right){\mu }_k+t\gamma {\mu }_0\theta \left(z^k\right)-\bar{\beta }\left(z^k+t\Delta z^k\right){\mu }_0\\ \ge \left(1-t\right){\bar{\beta }}_k{\mu }_0+t\gamma {\mu }_0\theta \left(z^k\right)-\gamma \theta \left(z^k+t\Delta z^k\right){\mu }_0\\ \ge \left(1-t\right)\gamma \theta \left(z^k\right){\mu }_0+t\gamma {\mu }_0\theta \left(z^k\right)-\gamma \left[1-\sigma \left(1-2\gamma {\mu }_0\eta \right)t\right]\theta \left(z^k\right){\mu }_0\\ =2\gamma \sigma \left(1-\gamma {\mu }_0\eta \right)t\theta \left(z^k\right){\mu }_0\ge 0.\end{array}$

ii) 若 $\theta \left(z^k\right)>1$，则 ${\beta }_k=\gamma {\text{e}}^{{\mu }_k},{\bar{\beta }}_k:=\bar{\beta }\left(z^k\right)=\gamma$ 。对任意的 $z\in R^{2n+1}$，有 $\bar{\beta }\left(z\right)\le \gamma$，从而对任意的 $t\in \left(0,1\right]$，有

$\begin{array}{l}{\mu }_k+t\Delta {\mu }_k-\bar{\beta }\left(z^k+t\Delta z^k\right){\mu }_0\\ ={\mu }_k+t\left(\frac{1-{\text{e}}^{{\mu }_k}}{{\text{e}}^{{\mu }_k}}+\frac{{\beta }_k{\mu }_0}{{\text{e}}^{{\mu }_k}}\right)-\bar{\beta }{\left(z^k+t\Delta z\right)}^k{\mu }_0\\ \ge \left(1-t\right){\mu }_k+t\gamma {\mu }_0-\gamma {\mu }_0\ge \left(1-t\right){\bar{\beta }}_k{\mu }_0+t\gamma {\mu }_0-\gamma {\mu }_0\\ =\left(1-t\right)\gamma {\mu }_0+t\gamma {\mu }_0-\gamma {\mu }_0=0.\end{array}$

从而 $z^{k+1}\in \Lambda$ 。

2) 不失一般性，设 $z^{*}$ 是点列 $\left\{z^k\right\}$ 的任意聚点，要证 $H\left(z^{*}\right)=0$，用反证法。假设 $\Vert H\left(z^{*}\right)\Vert >0$，由算法4.1的步骤3知 $\left\{\Vert H\left(z^k\right)\Vert \right\}$ 单调下降且有界，结合 $H\left(z^k\right)$ 的连续性，有 $\underset{k\to +\infty }{\lim }\Vert H\left(z^k\right)\Vert =\Vert H\left(z^{*}\right)\Vert$ 。由 $\beta (\cdot )$ 的定义知 $\beta \left(z^k\right)$ 单调下降且趋于 $\beta \left(z^{*}\right):{\beta }_{*}=\beta \left(z^{*}\right)={\text{e}}^{{\mu }_{*}}\gamma \min \left\{1,\theta \left(z^{*}\right)\right\}>0$ 。另外，由(5.1)和(4.6)式得到：

$\begin{array}{c}0<{\mu }_{k+1}={\mu }_k+t\Delta {\mu }_k\\ ={\mu }_k+t\left(\frac{1-{\text{e}}^{{\mu }_k}}{{\text{e}}^{{\mu }_k}}+\frac{{\beta }_k{\mu }_0}{{\text{e}}^{{\mu }_k}}\right)\\ \le {\mu }_k+t\left(-{\mu }_k{\text{e}}^{{\mu }_k}+\frac{{\text{e}}^{{\mu }_k}{\bar{\beta }}_k{\mu }_0}{{\text{e}}^{{\mu }_k}}\right)\\ \le {\mu }_k+t\left(-{\mu }_k{e}^{{\mu }_k}+{\mu }_k\right)\le {\mu }_k.\end{array}$

即 ${\mu }_k$ 是单调下降的，则 $\underset{k\to +\infty }{\lim }{\mu }_k={\mu }_{*}$ 。对 ${\mu }_k\ge {\beta }_k{\mu }_0$ 两边取极限，有 ${\mu }_{*}\ge {\beta }_{*}{\mu }_0>0$ 。由定理4.1知 $H^{\prime }\left(z^{*}\right)$ 存在且非奇异。由引理4.1知：存在 $z^{*}$ 的一个闭邻域 $N\left(z^{*}\right)$ 和正数 $\bar{t}\in \left(0,1\right]$，使得对于任意的 $z=\left(\mu ,x,y\right)\in N\left(z^{*}\right)$ 和所有的 $t\in \left[0,\bar{t}\right]$，都有 $\mu >0$， $H^{\prime }\left(z\right)$ 非奇异且

$\theta \left(z^k+t\Delta z^k\right)\le \left[1-\sigma \left(1-2\gamma {\mu }_0\eta \right)t\right]\theta \left(z^k\right).$

可以找到一个非负整数 $\tau$，使得 ${\delta }^{\tau }\in \left[0,\bar{t}\right]$，当k充分大时，有

$\theta \left(z^{k+1}\right)\le \left[1-\sigma \left(1-2\gamma {\mu }_0\eta \right)t\right]\theta \left(z^k\right).$

由步骤3知，对于充分大的k，步长 $t_k={\delta }^{{\alpha }_k}\ge {\delta }^{\tau }$，因此有

$\theta \left(z^{k+1}\right)\le \left[1-\sigma \left(1-2\gamma {\mu }_0\eta \right)t_k\right]\theta \left(z^k\right)\le \left[1-\sigma \left(1-2\gamma {\mu }_0\eta \right){\delta }^{\tau }\right]\theta \left(z^k\right).$

这与序列 $\left\{\theta \left(z^k\right)\right\}$ 的极限 $\theta \left(z^{*}\right)>0$ 矛盾，故 $\theta \left(z^{*}\right)=0$，即 $H\left(z^{*}\right)=0$，因此算法4.1是全局收敛的。

接下来证明算法的一个局部收敛性，在证明之前先介绍相关引理。

引理5.1 (文 [12] )：如果 $V\in \partial F\left(z^{*}\right)$ 都是非奇异的，那么存在 $z^{*}$ 的一个邻域 $N\left(z^{*}\right)$ 和一个常数 $C>0$，使得对于任意的 $z^k\in N{z}^{*}$ 和任意的 $V\in \partial F\left(z^k\right)$ 。V非奇异且 $\Vert V^{-1}\Vert \le C$ 。

定理5.2 (局部二次收敛)：设矩阵A行满秩且 $z^{*}$ 是算法4.1产生的迭代点列 $\left\{z^k\right\}$ 的任意聚点，则当k充分大时， $\left\{z^k\right\}$ 二阶收敛到 $z^{*}$，即 $\Vert z^{k+1}-z^{*}\Vert =O\left({\Vert z^k-z^{*}\Vert }^2\right)$，且 ${\mu }_{k+1}=O\left({\mu }_k^2\right)$ 。

证明：由定理5.1的(2)知 $H\left(z^{*}\right)=0$，又由定理4.1知 $H^{\prime }\left(z^{*}\right)$ 存在且非奇异。故由引理5.1知，存在一个常数 $C>0$，对于充分接近 $z^{*}$ 的所有点 $z^k$，有 $\Vert {\left(H^{\prime }\left(z^k\right)\right)}^{-1}\Vert \le C$ 。由定理3.1和定理4.1知 $H\left(z\right)$ 在 $z^k$ 处是全局利普希茨连续的强半光滑函数，因为 $H\left(x+\Delta x\right)-H\left(x\right)-V\left(\Delta x\right)=O\left({\Vert \Delta x\Vert }^{1+p}\right)$，当 $z^k$ 充分接近 $z^{*}$ 时，由半光滑函数的定义知，当 $z^k$ 充分接近点 $z^{*}$ 时，有

$\Vert H\left(z^k\right)\Vert =\Vert H\left(z^k\right)-H\left(z^{*}\right)\Vert =O\left(\Vert z^k-z^{*}\Vert \right)$

和

$H\left(z^{*}\right)-H\left(z^k\right)-H^{\prime }\left(z^k\right)\left(z^{*}-z^k\right)=O\left({\Vert z^k-z^{*}\Vert }^2\right).$

又 ${\beta }_k={\text{e}}^{{\mu }_k}\gamma \min \left\{1,{\Vert H\left(z^k\right)\Vert }^2\right\}\le {\Vert H\left(z^k\right)\Vert }^2$，即 ${\beta }_k=O\left({\Vert H\left(z^k\right)\Vert }^2\right)=O\left({\Vert z^k-z^{*}\Vert }^2\right)$ 。因此，

$\begin{array}{c}\Vert z^k+\Delta z^k-z^{*}\Vert =\Vert z^k+{\left(H^{\prime }\left(z^k\right)\right)}^{-1}\left({\beta }_k\bar{z}-H\left(z^k\right)\right)-z^{*}\Vert \\ =\Vert {\left(H^{\prime }\left(z^k\right)\right)}^{-1}\left[H^{\prime }\left(z^k\right)\left(z^k-z^{*}\right)+{\beta }_k\bar{z}-H\left(z^k\right)\right]\Vert \\ \le \Vert {\left(H^{\prime }\left(z^k\right)\right)}^{-1}\Vert \Vert H\left(z^{*}\right)-H\left(z^k\right)-H^{\prime }\left(z^k\right)\left(z^{*}-z^k\right)+{\beta }_k\bar{z}\Vert \\ \le C\left[\Vert H\left(z^{*}\right)-H\left(z^k\right)-H^{\prime }\left(z^k\right)\left(z^{*}-z^k\right)\Vert +{\beta }_k{\mu }_0\right]\\ =O\left({\Vert z^k-z^{*}\Vert }^2\right).\end{array}$ (5.2)

进一步有

$\Vert H\left(z^k+\Delta z^k\right)\Vert =O\left(\Vert z^k+\Delta z^k-z^{*}\Vert \right)=O\left({\Vert z^k-z^{*}\Vert }^2\right)=O\left({\Vert H\left(z^k\right)\Vert }^2\right)$ (5.3)

即当 $z^k$ 充分接近 $z^{*}$ 时，有 $z^{k+1}=z^k+\Delta z^k$ 。由(4.5)知

$\Vert z^{k+1}-z^{*}\Vert =O\left({\Vert z^k-z^{*}\Vert }^2\right)$ (5.4)

即算法4.1是二次收敛的，其次，当k充分大时， $H\left(z^k\right)$ 单调减少趋向于 $H\left(z^{*}\right)$，即 $H\left(z^k\right)\to 0$，所以 ${\beta }_k={\text{e}}^{{\mu }_k}\gamma {\Vert H\left(z^k\right)\Vert }^2$，又因为 $z^{k+1}=z^k+\Delta z^k$，所以当k充分大时有

${\mu }_{k+1}={\mu }_k+\Delta {\mu }_k=\frac{{\beta }_k{\mu }_0}{e^{{\mu }_k}}=\gamma {\Vert H\left(z^k\right)\Vert }^2{\mu }_0.$ (5.5)

结合(5.3)式知，当 $z^k$ 充分接近 $z^{*}$ 时，有

$\frac{{\mu }_{k+1}}{{\mu }_k}=\frac{{\Vert H\left(z^k\right)\Vert }^2}{{\Vert H\left(z^{k-1}\right)\Vert }^2}=\frac{O\left({\Vert H\left(z^{k-1}\right)\Vert }^4\right)}{{\Vert H\left(z^{k-1}\right)\Vert }^2}=O\left({\Vert H\left(z^{k-1}\right)\Vert }^2\right)=O\left({\mu }_k\right),$

即 ${\mu }_{k+1}=O\left({\mu }_k^2\right)$ 。

6. 数值试验

为了检验算法4.1的有效性，本文用Matlab 2014a编程，在Inter (R) Pentium (R) 4CPU，3.40 GHz，2 GB内存，Windows 7操作系统的电脑上做数值试验。测试问题是随机生成的规模n = 2 m从20到800维不等且r = 1的二阶锥规划问题。具体步骤为：首先生成随机的行满秩矩阵 $A\in R^{m\times n}$ 和随机向量 $x,s\in K_0$，然后令 $b=Ax,c=A^{\text{T}}y+s,$，则得到的二阶锥规划的原问题和对偶问题都存在最优解且最优值相等。在测试问题过程中，初始点选取 $x^0=e\in R^n,y^0=0\in R^m$ 。参数取值为： ${\mu }_0=0.01,\delta =0.65,\sigma =0.35,\gamma =0.90$ 。终止准则为 $\Vert H\left(z^k\right)\Vert \le {10}^{-6}$ 。计算结果列于表1，其中IT和CPU(s)表示每个测试问题进行10次所求得最优解时所用的平均迭代次数和平均时间。

Table 1. The numerical experiment results of the algorithm

表1. 算法的数值实验结果

由表1可以看出，算法是有效的，并且可以求解大规模的二阶锥规划问题，它只需较少的CPU时间就可以得到满足终止条件的解。并且本文中的算法4.1求解问题所用时间大部分比文献 [13] 所用时间少。

基金项目

国家自然科学基金(11561005)。