本文中, $R _+^n$和 $R_{++}^n$分别表示 $n$维列向量空间 $R^n$中的非负向量和正向量.对任意的 $x, y\in R^n$, $(x, y)$表示列向量 $(x^{\scriptsize\rm{T}}, y^{\scriptsize\rm{T}})^{\scriptsize\rm{T}}, $ $I$表示单位矩阵, $||\cdot||$表示向量 $x$的欧几里得范数, 即 $||x||=\sqrt{x^{\scriptsize\rm{T}}x}$.

二阶锥规划(SOCP)是在有限个二阶锥的笛卡尔乘积的仿射子空间之交上极大化或极小化一个线性函数.考虑标准形式的二阶锥规划问题

$ {{{\rm{(P)}}}}\min\big\{c^{\scriptsize\rm{T}}x:Ax=b, x\in {K}\big\} $

(1)

及其对偶问题

$ {{{\rm{(D)}}}}\max\big\{b^{\scriptsize\rm{T}}y:A^{\scriptsize\rm{T}}y+s=c, s\in {K}\big\}, $

(2)

其中 $A \in {R}^{m\times n}, c\in {R}^n$和 $b\in {R} ^m$是已知的量, ${K}\subset R^n$是有限个二阶锥的笛卡尔乘积, 即 ${K} = {K}^{n_1} \times {K}^{n_2}\times\cdots\times {K}^{n_r}, $这里 $\sum\limits_{i=1}^rn_i=n, $而 $n_i$维的二阶锥 ${K}^{n_i}$定义为

$ {K}^{n_i}:=\{(x_1, \bar{x})\in R\times R^{n_i-1}: x_{1}\geq\|\bar{x}\|\}. $

为分析简单, 本文假设 ${K}={K}^n$.所得结论可以很容易推广到一般情形.

二阶锥规划作为数学规划领域的一个重要分支, 有着非常重要而又广泛的应用背景和实际意义, 其研究问题涉及控制、金融、组合优化、工程技术、神经网络等诸多领域^[1].许多数学规划问题, 如凸二次规划、二次约束的凸二次规划、矩阵分式优化、范数极小化问题、鲁棒最小二乘问题等均可以转化为二阶锥规划求解^[1].因此, 近年来二阶锥规划成为数学规划领域一个值得关注的方向.

由于光滑牛顿法(非内点连续法)不但在理论上具有好的收敛性质, 而且在具体实施中具有很好的实际计算效果, 因而近年来得到了迅速的发展, 人们陆续提出了许多有效的光滑算法^[2-4].本文基于一个新的最小值函数的光滑函数, 给出一个求解二阶锥规划的光滑牛顿算法.算法可以起始于任意点, 在每步迭代中只需解一个线性方程组并进行一次线性搜索.在不需要满足严格互补假设条件下, 证明了所给算法是全局收敛和局部二阶收敛的.数值试验表明算法是有效的.

对任意的 $ x=(x_1, \bar{x})\in {R}\times {R}^{n-1}, s=(s_1, \bar{s})\in {R}\times {R}^{n-1}, $与二阶锥 $K$相伴的欧几里得约当代数定义为 $x\circ s=(x^{\scriptsize\rm{T}}s, x_{1}\bar{s}+s_{1}\bar{x}).$易知 $({R}^n, \circ)$是一个欧几里得约当代数, 其单位元为 ${{\bf{e}}}:=(1, 0, \cdots, 0)^{\scriptsize\rm{T}}\in R^{n}$.

对任意的 $x=(x_1, \bar{x})\in {R}\times {R}^{n-1}, $定义对称矩阵

$ L_x:= \left[\begin{array}{cc} x_1& \bar{x}^{\scriptsize\rm{T}}\\ \bar{x}& x_1I \end{array}\right], $

其中 $I$为 $(n-1)\times (n-1)$维的单位矩阵, 则有 $L_xs=x\circ s, \forall x, s\in R^n.$此外, $L_x$是半正定矩阵(正定矩阵)当且仅当 $x\in K(x\in K^0)$, 这里 $K^0$为 $K$的内部, 即

$ K^0:=\{(x_1, \bar{x})\in R\times{R} ^{n-1}: x_{1}>\|\bar{x}\|\}. $

对任意的 $x=(x_1, \bar{x})\in {R}\times {R}^{n-1}, $定义它的与二阶锥 $K$关联的特征值分解为 $x=\lambda_{1}u_1+\lambda_{2}u_2, $其中 $x$的特征值 $\lambda_i$及其对应的特征向量 $\mu_i$为

$ \lambda_i:=x_1+(-1)^i||\bar{x}||, u_{i}=\left\{ \begin{array}{c} \frac{1}{2}\bigg(1, (-1)^i\frac{\bar{x}}{||\bar{x}||}\bigg), \bar{x}\neq0, \\ \frac{1}{2}\bigg(1, (-1)^i\omega\bigg), \ \ \bar{x}=0, \\ \end{array} \right.i=1, 2, $

这里 $\omega\in R^{n-1}$是满足 $||\omega||=1$的任意向量.

利用 $x$的特征值分解可以定义:

$ \begin{eqnarray*}&&\sqrt{x}:=\sqrt{\lambda_{1}}u_1+\sqrt{\lambda_{2}}u_2, \forall x\in K;\\ &&x^{-1}:=\lambda_{1}^{-1}u_1+\lambda_{2}^{-1}u_2, \forall x\in K^0;\\ &&x^2:=\lambda_{1}^2u_1+\lambda_{2}^2u_2, \forall x\in R^n.\end{eqnarray*} $

容易证明 $x^2=x\circ x, (\sqrt{x})^2=x$, 并且对任意的 $x\in R^n, $有 $x^2\in K$.

本文作如下假设:二阶锥规划问题(1) 和(2) 存在严格可行解.

在此假设条件下, 由二阶锥规划的强对偶定理^[1]可知, 解(1) 和(2) 等价于解最优性条件

$ Ax=b, A^{\scriptsize\rm{T}}y+s=c, x\circ s=0, x, s\in K. $

(3)

文献[5]已经证明向量最小值函数 $\phi_{\scriptsize{{min}}}(x, s): R^n\times R^n \rightarrow R^n $:

$ \phi_{\scriptsize{{min}}}(x, s)=x+s-\sqrt{(x-s)^2} $

满足重要性质

$ \phi_{\scriptsize{{min}}}(x, s)=0 \Longleftrightarrow x\in K, s\in K, x\circ s=0. $

(4)

然而, $\phi_{\scriptsize{{min}}}$是非光滑的, 因为它在 $(0, 0)\in R^n\times R^{n}$点处是不可微的, 这大大影响了其实际应用.

本文通过光滑 $\phi_{\scriptsize{{min}}}(x, s), $得到一个新的向量值函数 $\phi(\mu, x, s): {R} _{+}\times {R} ^n\times {R} ^n\rightarrow {R} ^n$:

$ \phi(\mu, x, s)=(e^{\mu}+\mu)(x+s)-\sqrt{(e^{\mu}-\mu)^2(x-s)^2+4\mu^2{{\bf{e}}}}. $

(5)

定理3.1  设 $\phi(\mu, x, s)$由式(5) 定义, 则下列结论成立.

(ⅰ)函数 $\phi(\mu, x, s)$在 $R^{1+2n}$上处处强半光滑.此外, $\phi(\mu, x, s)$在任意的 $(\mu, x, s)\in R_{++}\times R^n\times R^n$处是连续可微的, 其雅可比矩阵为

$ \phi'(\mu, x, s)=\left( \begin{array}{c}(e^{\mu}+1)(x+s)-L_{Q}^{-1}[(e^{\mu}-\mu)(e^{\mu}-1)q^2+4\mu {\bf{e}}]\\ (e^{\mu}+\mu)I-(e^{\mu}-\mu)^2L_{Q}^{-1}L_{q}\\ (e^{\mu}+\mu)I+(e^{\mu}-\mu)^2L_{Q}^{-1}L_{q} \end{array} \right), $

(6)

其中 $q:=x-s, Q:=\sqrt{(e^{\mu}-\mu)^2q^2+4\mu^2{\bf{e}}}.$

(ⅱ)对任意的 $(x, s)\in R^n\times R^n$, 有

$ \lim\limits_{\mu\downarrow0}\phi(\mu, x, s)=\phi_{\min}(x, s). $

因此, $\phi(\mu, x, s)$是 $\phi_{\min}(x, s)$的光滑函数.

证  由文献[6]中的定理3.2易证 $\phi(\mu, x, s)$是 $R^{1+2n}$上的强半光滑函数, 且在任意的 $(\mu, x, s)\in R_{++}\times R^n\times R^n$处是连续可微的.下证(6) 式成立.对任意的 $(\mu, x, s)\in R_{++}\times R^n\times R^n$有

$ Q^2=(e^{\mu}-\mu)^2q^2+4\mu^2{\bf{e}}. $

对上式左右两端同时求导得

$ Q'(\mu, x, s)=\left( \begin{array}{c}L_{Q}^{-1}[(e^{\mu}-\mu)(e^{\mu}-1)q^2+4\mu {\bf{e}}]\\ (e^{\mu}-\mu)^2L_{Q}^{-1}L_{q}\\ -(e^{\mu}-\mu)^2L_{Q}^{-1}L_{q} \end{array} \right), $

从而由复合函数的求导法则可得(6) 式.结论(ⅱ)的证明类似于文献[4]中的定理2.4, 在此省略.证毕.

记 $z:=(\mu, x, y)\in R\times R^n\times R^m$.利用光滑函数(5), 定义函数 $H(z):R^{1+n+m}\rightarrow R^{1+m+n}:$

$ H(z)=\left( \begin{array}{c} \mu\\ b-Ax\\ \phi(\mu, x, c-A^{\scriptsize\rm{T}}y) \end{array} \right). $

(7)

由式(3) 和(4) 以及定理3.1中的(ⅱ)可知 $z^*=(\mu^*, x^*, y^*)$是 $H(z)=0$的解当且仅当 $( x^*, y^*, c-A^{\scriptsize\rm{T}}y^*)$满足最优性条件(3), 即 $( x^*, y^*, c-A^{\scriptsize\rm{T}}y^*)$是(1) 和(2) 的最优解.因此我们可以利用牛顿法求解方程组 $H(z) =0$.

令 $\gamma\in(0, 1)$, 对任意的 $z=(\mu, x, y)\in {R}\times {R}^n\times {R}^m$, 定义函数

$ \rho(z):=\gamma\min\{1, \|H(z)\|^2\}. $

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

步骤0  取常数 $\delta, \sigma\in (0, 1)$, $\mu_0\in {R}_{++}$.记 $\bar{z}:=(\mu_0, 0, 0)\in{R}_{++}\times {R}^n\times {R}^m$.选取任意的初始点 $(x_0, y_0)\in{R}^n\times {R}^m$.令 $z_0:=(\mu_0, x_0, y_0)$.选取参数 $\gamma\in (0, 1)$满足 $\mu_0\gamma<1$.令 $k:=0.$

步骤1  如果 $||H(z_k)||=0$, 则停止迭代.

步骤2  通过求解方程组

$ H(z_k)+H'(z_k)\Delta z_k=\rho(z_k) \bar{z} $

(8)

得到搜索方向 $\Delta z_k:=(\Delta \mu_k, \Delta x_k, \Delta y_k)\in {R}\times {R}^n\times {R}^m$.

步骤3  设 $l_k$是满足

$ ||H(z_k+\delta^l\Delta z_k)||\leq[1-\sigma(1-\mu_0\gamma)\delta^l]||H(z_k)|| $

(9)

最小的非负整数.令 $\alpha_k:=\delta^{l_k}.$

步骤4  令 $z_{k+1}:=z_k+\alpha_k\Delta z_k$和 $k:=k+1.$转步骤1.

定理4.1  设 $H(z)$由式(7) 定义, 则

(ⅰ) $H(z)$是 $R^{1+n+m}$上的强半光滑函数, 且在任意点 $z=(\mu, x, y)\in{R}_{++}\times {R}^n\times {R}^m$处是连续可微的, 其雅可比矩阵为

$ H'(z)=\left( \begin{array}{c} &1\ \ \ \ \ \ &0\ \ \ \ \ \ &0\\ &0\ \ \ \ \ \ &-A\ \ \ \ \ \ &0\\ &P(z)\ \ \ \ \ \ &M(z)\ \ \ \ \ \ &-N(z)A^{\scriptsize{{\rm{T}}}}\\ \end{array} \right), $

其中

$ \begin{eqnarray*}&&P(z)=(e^{\mu}+1)(x+c-A^{\scriptsize{{\rm{T}}}}y)-L_{\bar{Q}}^{-1}[(e^{\mu}-\mu)(e^{\mu}-1)\bar{q}^2+4\mu {\bf{e}}], \\ &&M(z)=(e^{\mu}+\mu)I-(e^{\mu}-\mu)^2L_{\bar{Q}}^{-1}L_{\bar{q}}, \\ &&N(z)=(e^{\mu}+\mu)I+(e^{\mu}-\mu)^2L_{\bar{Q}}^{-1}L_{\bar{q}}, \\ &&\bar{q}:=x-(c-A^{\scriptsize{{\rm{T}}}}y), \\ &&\bar{Q}:=\sqrt{(e^{\mu}-\mu)^2\bar{q}^2+4\mu^2{\bf{e}}}.\end{eqnarray*} $

(ⅱ)如果矩阵 $A$行满秩, 则对任意的 $\mu> 0$有 $H'(z)$可逆.

证  由定理3.1易知(ⅰ)成立.结论(ⅱ)的证明类似于文献[4]中的定理3.1, 在此省略.证毕.

定理4.2  设矩阵 $A$行满秩, 如果 $\mu_k>0$, 则算法4.1有好的定义.

证  因为矩阵 $A$行满秩且 $\mu_k>0$, 故由定理4.1可知 $H'(z_k)$可逆, 所以在第 $k$步迭代步骤2有好的定义.令 $\Delta z_k=(\Delta \mu_k, \Delta x_k, \Delta y_k)\in {R}\times {R}^n\times {R}^m$是方程组(8) 的解, 则对任意的 $\alpha\in (0, 1)$, 有

$ \mu_k+\alpha\Delta\mu_k=\mu_k+\alpha(-\mu_k+\mu_0\rho(z_k)) =(1-\alpha)\mu_k+\alpha\mu_0\rho(z_k)>0. $

(10)

由(10) 式及定理4.1可知 $H(\cdot)$在 $z_k$附近连续可微.对任意的 $\alpha\in (0, 1)$, 定义

$ g(\alpha):=H(z_k+\alpha \Delta z_k)-H(z_k)-\alpha H'(z_k) \Delta z_k, $

(11)

则 $||g(\alpha)||=o(\alpha)$.由函数 $\rho(\cdot)$的定义知如果 $\|H(z^k)\|\geq1$, 则 $\rho(z_k)=\gamma\leq\gamma||H(z_k)||$, 否则 $\rho(z_k)=\gamma||H(z_k)||^2\leq\gamma||H(z_k)||$.这表明对任意的 $k\geq0$, 有 $\rho(z_k)\leq\gamma||H(z_k)||$.由式(8) 和(11) 可知对任意的 $\alpha\in (0, 1)$有

$ \begin{eqnarray*} ||H(z_k+\alpha \Delta z_k)|| &=&||H(z_k)+\alpha H'(z_k) \Delta z_k+g(\alpha)||\\ &\leq&(1-\alpha)||H(z_k)||+\alpha\mu_0\rho(z_k)+||g(\alpha)||\\ &\leq&(1-\alpha)||H(z_k)||+\alpha\mu_0\gamma||H(z_k)||+o(\alpha)\\ &=&[1-(1-\mu_0\gamma)\alpha]||H(z_k)||+o(\alpha). \end{eqnarray*} $

因为 $\mu_0\gamma < 1$, 所以存在一个常数 $\bar{\alpha}\in (0, 1)$使得对任意的 $\alpha\in(0, \bar{\alpha}]$, 有

$ ||H(z_k+\alpha \Delta z_k)||\leq[1-\sigma(1-\mu_0\gamma)\alpha]||H(z_k)||, $

即步骤3在第 $k$步迭代有好的定义.综上所述可知, 算法4.1有好的定义.证毕.

为证明算法4.1的全局收敛性, 首先给出如下引理.

引理5.1  设矩阵 $A$行满秩且 $\{z_k:=(\mu_k, x_k, y_k)\}$是算法4.1生成的无穷迭代点列, 则对任意的 $k\geq0$有 $\mu_k>0$且 $z_k\in \Omega$, 其中

$ \Omega:=\Big\{z=(\mu, x, y)\in {R}_{++}\times {R}^n \times {R}^m| \mu\geq \mu_0\rho(z)\Big\}. $

证  假设 $\mu_k>0$, 则由(10) 式得

$ \mu_{k+1} =(1-\alpha_k)\mu_k+\alpha_k\mu_0\rho(z_k)>0. $

(12)

因此, 由(12) 式及 $\mu_0>0$知对任意的 $k\geq0$有 $\mu_k>0$.

下面用数学归纳法证明对任意的 $k\geq0$有 $z_k\in \Omega$.因为 $\rho(z_0)\leq \gamma < 1, $故 $\mu_0\geq\mu_0\rho(z_0)$, 从而可知 $z_0\in \Omega$.假设 $z_k\in\Omega$, 则有 $\mu_k\geq\mu_0\rho(z_k).$由算法4.1的步骤3和步骤4可知 $||H(z_{k+1})||\leq||H(z_k)||$, 从而可得

$ \begin{eqnarray*} \mu_{k+1}&=&(1-\alpha_k)\mu_k+\alpha\mu_0\rho(z_k)\\ &\geq&(1-\alpha_k)\mu_0\rho(z_k)+\alpha\mu_0\rho(z_k)\\ &=&\mu_0\gamma\min\{1, ||H(z_k)||^2\}\\ &\geq&\mu_0\gamma\min\{1, ||H(z_{k+1})||^2\}. \end{eqnarray*} $

综上所述可知任意的 $k\geq0$有 $z_k\in \Omega$.证毕.

定理5.1  假设矩阵 $A$行满秩且 $\{z_k\}$是由算法4.1产生的无穷点列, 则 $\{z_k\}$的任意聚点 $z^*=(\mu^*, x^*, y^*)$都是 $H(z)=0$的解.

证  不失一般性, 假设当 $k\rightarrow+\infty$时 $z_k$收敛于 $z^*$.由算法4.1的步骤3和步骤4可知

$ ||H(z_{k+1})||\leq||H(z_k)||. $

又因为 $z_k\in\Omega$, 由引理5.1可知

$ \mu_{k+1} =(1-\alpha_k)\mu_k+\alpha_k\mu_0\rho(z_k)\leq\mu_k, $

故 $\{||H(z_k)||\}$和 $\{\mu_k\}$都是单调递减序列且有下界0.因此当 $k\rightarrow+\infty$时, $||H(z_k)||\rightarrow||H(z^*)||$且 $\mu_k\rightarrow\mu^*$.若 $||H(z^*)||=0, $定理得证.假设 $||H(z^*)||>0, $则由 $\rho(\cdot)$定义可知 $\rho(z_k)$收敛于 $\rho^*:=\gamma\min\{1, ||H(z^*)||^2\}>0.$因为 $z_k\in\Omega$, 故 $\mu^*\geq\mu_0\rho^*>0$, 从而由定理4.1知 $H'(z^*)$存在且可逆, 因此存在 $z^*$的一个闭的邻域 $N(z^*)$使得对任意的 $z\in N(z^*)$, 有 $\mu\in R_{++}$且 $H'(z)$可逆.对 $z\in N(z^*)$, 令 $\Delta z=(\Delta\mu, \Delta x, \Delta y)$是方程组

$ H(z)+H'(z)\Delta z=\rho(z)\bar{z} $

的唯一解.定义

$ g_z(\alpha):=H(z+\alpha \Delta z)-H(z)-\alpha H'(z) \Delta z. $

类似于定理4.2的证明, 存在 $\bar{\alpha}\in (0, 1)$使得对任意的 $\alpha\in(0, \bar{\alpha}]$及 $z\in N(z^*)$, 有

$ ||H(z+\alpha \Delta z)||\leq[1-\sigma(1-\mu_0\gamma)\alpha]||H(z)|| . $

因此对于充分大的 $k$, 存在一个非负整数 $\nu$使得 $\delta^{\nu}\in (0, \bar{\alpha}]$且

$ ||H(z_k+\delta^{\nu}\Delta z_k)||\leq[1-\sigma(1-\mu_0\gamma)\delta^{\nu}]||H(z_k)||. $

对于充分大的 $k, $因为 $\alpha_k=\delta^{l_k}\geq \delta^{\nu}$, 故由算法4.1的步骤3和步骤4可知

$ ||H(z_{k+1})||\leq[1-\sigma(1-\mu_0\gamma)\delta^{l_k}]||H(z_k)||\leq[1-\sigma(1-\mu_0\gamma)\delta^{\nu}]||H(z_k)||, $

即 $||H(z_{k+1})||\leq C||H(z_k)||$, 其中 $C=1-\sigma(1-\mu_0\gamma)\delta^{\nu}<1$为常数.因此, 当 $k\rightarrow\infty$时, $\|H(z_k)\|\rightarrow0$.这与 $\|H(z_k)\|\rightarrow\|H(z^*)\|>0$矛盾.因此, $\|H(z^*)\|=0$.证毕.

为了给出算法4.1的局部二阶收敛性, 需要下面的假设.

假设5.1  假设 $z^*$满足非奇异假设, 即所有的 $V\in \partial H(z^*)$都是非奇异的, 其中 $\partial H(z)$表示函数 $H(z)$在 $z$点处的Clark意义下的广义雅克比矩阵^[7].

定理5.2  如果矩阵 $A$行满秩并且假设5.1成立, $z^*$是算法4.1产生的迭代点列 $\{z_k\}$的任意聚点, 则 $\{z_k\}$二阶收敛到 $z^*$, 即 $||z_{k+1}-z^*||=O(||z_k-z^*||^2), $并且 $\mu_{k+1}=O(\mu_k^2).$

证  证明过程类似于文献[4]中的定理4.3, 在此省略.

为检验算法4.1的实算效果, 我们用Matlab7.1编程, 在Intel(R) Pentium(R)4 CPU, 2.00 GHz, 512MB内存, Windows XP操作系统的电脑上做数值试验.测试问题是随机生成的规模 $n(=2m)$从100到500不等且 $r=1$的二阶锥规划问题.具体步骤为:首先生成随机的行满秩矩阵 $A\in R^{m\times n}$和随机向量 $x, s\in K^0$, 然后令 $b:=Ax, c:=A^{\scriptsize{{\rm{T}}}} y+s$, 则得到的二阶锥规划的原问题和对偶问题都存在最优解且最优值相等.在测试过程中, 参数取值为 $\mu_0=0.01, \sigma=0.25, \delta=0.75, \gamma=0.95.$停止准则为 $||H(z_k)||\leq 10^{-7}.$计算结果列于表 1, 其中 ${{ones}}(m, 1):=(1, \cdots, 1)^{\scriptsize\rm{T}}\in R^m$, ${{\bf{e}}}_m:=(1, 0, \cdots, 0)^{\scriptsize\rm{T}}\in R^m$, IT和CPU分别表示算法4.1所需的迭代次数和CPU时间(单位:秒), FV和MU分别表示算法4.1终止时的 $||H(z_k)||$和 $\mu_k$的值.

由表 1可以看出, 算法4.1对初始点的选取没有任何限制, 并且能够求解大规模的二阶锥规划问题, 它只需很少的迭代次数和CPU时间就可以得到满足终止条件的解.此外, 我们还发现算法4.1求解问题所需的迭代次数几乎不受问题规模的影响, 但所需的CPU时间会随问题规模的变大而增加.