广西科学  2016, Vol. 23 Issue (5): 396-403   PDF    
二次半定规划一个原始对偶路径跟踪算法
黎健玲, 王培培     
广西大学数学与信息科学学院,广西 南宁 530004
摘要: 本文提出求解二次半定规划的一个基于H..K..M方向的原始对偶路径跟踪算法.文中首先导出确定H..K..M方向的线性方程组,并证明该搜索方向的存在唯一性;然后给出算法的具体步骤,并证明算法产生的迭代点列落在中心路径的某个邻域内.最后采用Matlab(R2011b)数学软件编程对算法进行数值试验.数值结果表明算法是有效的.
关键词: 二次半定规划原始对偶     算法     路径跟踪     中心路径    
A Primal-dual Path-following Algorithm for Quadratic Semi-definite Programming
LI Jianling , WANG Peipei     
College of Mathematics and Information Science,Guangxi University,Nanning,Guangxi,530004,China
Abstract: A primal-dual path-following algorithm based on H..K..M direction for quadratic semi-definite programming problems(QSDP) is proposed.Firstly,the system of linear equations yielding the H..K..M direction are derived,and the existence and uniqueness of the search direction are shown;Secondly,the algorithm is described in detail.We show that the iterates generated by the algorithm can fall into some neighborhood of the central path under some mild conditions.Finally,a preliminary numerical experiment is performed for the algorithm by using Matlab (R2011b) mathematical software,and the numerical results show that the proposed algorithm is effective.
Key words: quadratic semi-definite programming     primal-dual     algorithm     path-following     central path    
0 引言

研究意义】(线性)半定规划(简记SDP)既是定义在半正定矩阵锥上的规划,也是线性规划、凸二次规划、二阶锥规划等的一种推广,在控制论、组合优化等诸多领域有着广泛的应用.而二次半定规划是半定规划的拓展,它在最优控制、证券、金融风险分析等领域中都有着广泛的应用.【前人研究进展】原始对偶内点算法是求解半定规划的有效方法[1-10].由于确定搜索方向的方法不同,因此导出多种形式的原始对偶内点法.目前常用的搜索方向有以下3种:(i)AHO搜索方向[2];(ii)H..K..M搜索方向[3-5];(iii)Nesterov-Todd搜索方向(简称NT方向)[6-7].基于SDP的原始对偶内点算法,文献[11]提出一个求解二次半定规划的势减少(potential reduction)算法;文献[12]结合Dikin型方向和牛顿中心步提出一个求解二次半定规划的预估校正算法;文献[13]提出一个求解二次半定规划的投影收缩算法;文献[14]将基于NT方向的原始对偶内点算法推广到二次半定规划.【本研究切入点】本文考虑如下二次半定规划问题(简记为QSDP):

$\begin{array}{*{20}{l}} {minf\left( X \right) = \frac{1}{2}svec{{\left( X \right)}^T}{H^T}Hsvec\left( X \right)}\\ { - {a^T}Hsvec\left( X \right) + C \cdot X} \end{array}$ (0.1)
$\begin{array}{*{35}{l}} s.t.~{{A}_{i}}\cdot X={{b}_{i}},i=1,\cdots ,m; \\ {{X}^{T}}=X\succcurlyeq 0 \\ \end{array}$

其中HT=(svec(H1),…,svec(Hl))∈R $\frac{1}{2}$ n(n+1)×l,Hj(j=1,…,l),Ai(i=1,…,m),CX都是n×n阶实对称矩阵,aRl,bRm,算子svec(X)将n×n阶对称阵X的列堆叠(stack)成一个 $\frac{1}{2}$ n(n+1)维的列向量,它的逆算子记为smat.X$\succcurlyeq $($\succ $)0表示X是半正定(正定)矩阵.对任意ARp×q,B∈Rp×q,矩阵内积A· B=Tr(ATB)= $\sum\limits_{i=1}^{p}{{}}\sum\limits_{j=1}^{q}{{}}$ aijbij,其中Tr(·)表示矩阵的迹.若A,Bn阶对称阵,则有A·B=svec(A)Tsvec(B).分别记Sn,S+nS++nn阶实对称矩阵全体,n阶实对称半正定矩阵全体和n阶实对称正定矩阵的全体.因为HTH是对称半正定矩阵,所以(0.1)式的目标函数是凸的,而约束函数是非光滑但也是凸的,因此QSDP(0.1)是一个凸优化问题. 据此,本文基于文献[5]的思想,提出一个求解QSDP(0.1)的原始对偶路径跟踪算法.【拟解决的关键问题】首先,利用Wolfe对偶理论得到QSDP(0.1)的对偶规划,从而导出对偶间隙的表达式,并据此定义中心路径.然后结合最优性条件、矩阵对称化技术以及牛顿法导出确定H..K..M方向的线性方程组,并证明该搜索方向的存在唯一性.随后基于该搜索方向给出短步路径跟踪算法的具体步骤,并证明算法产生的迭代点列落在中心路径的某个邻域内.

1 二次半定规划问题的对偶规划

利用凸规划问题的Wolfe对偶理论,可得QSDP(0.1)的对偶规划(简记为DSDP)为

$\begin{array}{*{35}{l}} \begin{align} & max~D\left( X,y \right)=-\frac{1}{2}svec{{\left( X \right)}^{T}}{{H}^{T}}Hsvec\left( X \right)+ \\ & \sum\limits_{i=1}^{m}{{}}{{b}_{i}}~{{y}_{i}} \\ \end{align} \\ s.t.\sum\limits_{i=1}^{m}{{}}{{y}_{i}}~{{A}_{i}}+Z=C-\sum\limits_{j=1}^{l}{{}}{{a}_{j}}~{{H}_{j}}+\sum\limits_{j=1}^{l}{{}}{{H}_{j}}\left( {{H}_{j}}\cdot X \right), \\ X\succcurlyeq 0,Z\succcurlyeq 0. \\ \end{array}$ (1.1)

分别记问题(0.1),(1.1)的可行集为FP,FD,即

$\begin{array}{*{35}{l}} {{F}_{P}}=\left\{ X|{{A}_{i}}\cdot X\text{ }={{b}_{i}},X\succcurlyeq 0 \right\}, \\ \begin{align} & {{F}_{D}}=\left\{ \left( X,y,Z \right)|\sum\limits_{i=1}^{m}{{}}{{y}_{i}}~{{A}_{i}}+Z=C- \right. \\ & \left. \sum\limits_{j=1}^{l}{{}}{{a}_{j}}{{H}_{j}}+\sum\limits_{j=1}^{l}{{}}{{H}_{j}}\left( {{H}_{j}}\cdot X \right),X\succcurlyeq 0,Z\succcurlyeq 0 \right\} \\ \end{align} \\ \end{array}$

FPFD的严格内部为FP°和FD°.

易知DSDP(1.1)也是一个凸二次半定规划问题.对于任意的(X,y,Z)∈FD,若X∈FP,则对偶间隙为

$\begin{align} & f\left( X \right)-D\left( X,y \right)=svec{{\left( X \right)}^{T}}{{H}^{T}}Hsvec\left( X \right)-{{a}^{T}} \\ & Hsvec\left( X \right)+C\cdot X-\sum\limits_{i=1}^{m}{{}}1{{b}_{i}}~{{y}_{i}}, \\ & =svec{{\left( X \right)}^{T}}{{H}^{T}}Hsvec\left( X \right)-{{a}^{T}}Hsvec\left( X \right)-\sum\limits_{i=1}^{m}{{}}{{b}_{i}}~{{y}_{i}}~+ \\ & (\sum\limits_{i=1}^{m}{{}}{{y}_{i}}~{{A}_{i}}+Z+\sum\limits_{j=1}^{l}{{}}{{a}_{j}}{{H}_{j}}-\sum\limits_{j=1}^{l}{{}}{{H}_{j}}({{H}_{j}}\cdot \text{ }X))\cdot \text{ }X, \\ & =svec{{\left( X \right)}^{T}}~{{H}^{T}}Hsvec\left( X \right)-\sum\limits_{j=1}^{l}{{}}{{({{H}_{j}}\cdot \text{ }X)}^{2}}+Z\cdot \text{ }X= \\ & {{\left( Hsvec\left( X \right) \right)}^{T}}\left( Hsvec\left( X \right) \right)-\sum\limits_{j=1}^{l}{{}}{{({{H}_{j}}\cdot \text{ }X)}^{2}}+Z\cdot \text{ }X \\ & =Z\cdot X. \\ \end{align}$ (1.2)
2 H..K..M搜索方向

假设:

(A1) 设矩阵组{Ai,i=1,…,m}线性无关.

(A2) Slater条件成立,即存在X$\succ $0,Z$\succ $0和yRm使得XFP,(X,y,Z)∈FD.

由(1.2)式以及文献[13]的定理3.5知,在假设(A2)下,X∈FP是QSDP(0.1)的最优解当且仅当存在y∈Rm,ZSn,使得

$\sum\limits_{i=1}^{m}{{}}{{y}_{i}}{{A}_{i}}+Z=C-\sum\limits_{j=1}^{l}{{}}{{a}_{j}}~{{H}_{j}}+\sum\limits_{j=1}^{l}{{}}{{H}_{j}}({{H}_{j}}\cdot X)$ (2.1a)
$XZ=0$ (2.1b)

点对(X,Z)如果满足XFP,(X,y,Z)∈FD

$XZ=\sigma \mu I,$

其中参数μ>0,μ= $\frac{X\cdot Z}{n}$ ,σ是中心参数,In阶单位阵,则称点对(X,Z)在中心路径上,即中心路径上的点对(X,Z)满足:

${{A}_{i}}\cdot X={{b}_{i}},i=1,\cdots ,m,$ (2.2a)
$\sum\limits_{i=1}^{m}{{}}{{y}_{i}}{{A}_{i}}+Z=C-\sum\limits_{j=1}^{l}{{}}{{a}_{j}}{{H}_{j}}+\sum\limits_{j=1}^{l}{{}}{{H}_{j}}({{H}_{j}}\cdot \text{ }X)$ (2.2b)
$XZ=\sigma \mu I.$ (2.2c)

设当前迭代点为(X,y,Z),且XFP°,(y,Z)∈FD°,用一步Newton法求解非线性方程组(2.2)得到的系统如下:

${{A}_{i}}\cdot \Delta X=0,i=1,\cdots ,m$ (2.3a)
$\sum\limits_{i=1}^{m}{{}}\Delta {{y}_{i}}{{A}_{i}}+\Delta Z-\sum\limits_{j=1}^{l}{{}}{{H}_{j}}({{H}_{j}}\cdot \Delta X)=0,$ (2.3b)
$\Delta XZ+X\Delta Z=\sigma \mu I-XZ-\Delta X\Delta Z.$ (2.3c)

忽略(2.3c)式中的ΔXΔZ,得到

$\Delta XZ+X\Delta Z=\sigma \mu I-XZ$ (2.4)

注意到X,Z不可交换,即XZZX,而内点法的一个关键是产生的矩阵ΔXZ要满足对称性.文献[9]给出对称化矩阵M的一般通式:

${{H}_{P}}\left( M \right)=\frac{1}{2}(PM{{P}^{-1}}+{{(PM{{P}^{-1}})}^{T}}),$

其中P是可逆阵,有多种不同的取法[1-5].本文采用文献[1, 5]中的方法,取P=X $\frac{1}{2}$ .对称化方程(2.4),得到如下的线性方程组:

${{A}_{i}}\cdot \Delta X=0,i=1,\cdots ,m$ (2.5a)
$\sum\limits_{i=1}^{m}{{}}\Delta {{y}_{i}}~{{A}_{i}}+\Delta Z-\sum\limits_{j=1}^{l}{{}}{{H}_{j}}({{H}_{j}}\cdot \Delta X)=0,$ (2.5b)
$\begin{align} & {{X}^{-\frac{1}{2}}}\left( X\Delta Z+\Delta XZ \right){{X}^{\frac{1}{2}}}+{{X}^{\frac{1}{2}}}\left( \Delta ZX+Z\Delta X \right) \\ & {{X}^{-\frac{1}{2}}}~=2(\sigma \mu I-{{X}^{\frac{1}{2}}}Z{{X}^{\frac{1}{2}}}). \\ \end{align}$ (2.5c)

记W=σμI-X $\frac{1}{2}$ ZX $\frac{1}{2}$ ,求解线性方程组(2.5)得到的解(ΔXyZ)称为H..K..M搜索方向.记

${{G}^{T}}=(svec({{A}_{1}}),\cdots ,svec({{A}_{m}})),$ (2.6a)
${{r}_{c}}=svec(\sigma \mu I-{{X}^{\frac{1}{2}}}Z{{X}^{\frac{1}{2}}})=svecW,$ (2.6b)
$\begin{align} & E={{X}^{-\frac{1}{2}}}{{\otimes }_{s}}{{X}^{\frac{1}{2}}}Z=\left( I{{\otimes }_{s}}{{X}^{\frac{1}{2}}}Z{{X}^{\frac{1}{2}}} \right)\cdot \\ & ({{X}^{-\frac{1}{2}}}{{\otimes }_{s}}{{X}^{-\frac{1}{2}}}), \\ \end{align}$ (2.6c)
$F={{X}^{\frac{1}{2}}}{{\otimes }_{s}}{{X}^{\frac{1}{2}}}$ (2.6d)

其中X $\frac{1}{2}$ sX $\frac{1}{2}$ 表示两个矩阵的对称Kronecker积(文献[1]附录).使用算子svec,可把线性方程组(2.5)改写成以下形式:

$\left( {\begin{array}{*{20}{c}} 0&G&0\\ {{G^T}}&{ - {H^T}H}&I\\ 0&E&F \end{array}} \right)\left( \begin{array}{l} \Delta y\\ svec\left( {\Delta X} \right)\\ svec\left( {\Delta Z} \right) \end{array} \right) = \left( \begin{array}{l} 0\\ 0\\ {r_c} \end{array} \right),$ (2.7)

其中I $\frac{{n\left( {n + 1} \right)}}{2}$ 阶单位阵.根据假设(A1)知G是行满秩的.

引理2.1 假设(A1)和(A2)成立,X,ZS++n,则矩阵FHTH+E可逆.

证明 已知

$\begin{array}{l} E = {X^{ - \frac{1}{2}}}{ \otimes _s}{X^{\frac{1}{2}}}Z = \left( {I{ \otimes _s}{X^{\frac{1}{2}}}Z{X^{\frac{1}{2}}}} \right)\cdot\\ ({X^{ - \frac{1}{2}}}{ \otimes _s}{X^{ - \frac{1}{2}}}),\\ F = {X^{\frac{1}{2}}}{ \otimes _s}{X^{\frac{1}{2}}}, \end{array}$

X,Z是对称正定阵知X ${\frac{1}{2}}$ ,X- ${\frac{1}{2}}$ 也是对称正定阵,进而由对称Kronecker积的性质知E,F可逆.而

$\begin{array}{l} rank(F{H^T}H + E) = rank({F^{ - 1}}(F{H^T}H + E)) = \\ rank({H^T}H + {F^{ - 1}}E), \end{array}$ (2.8)

往证F-1E是对称正定矩阵.由对称Kronecker积的性质(文献[1]附录)知

$\begin{array}{l} {F^{ - 1}}E = {\left( {{X^{\frac{1}{2}}}{ \otimes _s}{X^{12}}} \right)^{ - 1}}\left( {{X^{ - \frac{1}{2}}}{ \otimes _s}{X^{\frac{1}{2}}}Z} \right) = \\ ({X^{ - \frac{1}{2}}}{ \otimes _s}{X^{ - \frac{1}{2}}})({X^{ - \frac{1}{2}}}{ \otimes _s}{X^{\frac{1}{2}}}Z) = \frac{1}{2}[({X^{ - 1}}{ \otimes _s}Z) + \\ (Z{ \otimes _s}{X^{ - 1}})]{\rm{ }} = {X^{ - 1}}{ \otimes _s}Z, \end{array}$ (2.9)

X-1和Z是对称正定阵,故X-1sZ,即F-1E是对称正定阵.又因为HTH是对称半正定阵,因此HTH+F-1E是对称正定阵,结合(2.8)式即知矩阵FHTH+E是可逆的.

引理2.2 假设(A1)和(A2)成立,X,ZS++n,则矩阵(FHTH+E)-1FS++n.

证明 因为

$\begin{array}{l} {\left( {F{H^T}H + E} \right)^{ - 1}}F = {\left( {{F^{ - 1}}\left( {F{H^T}H + E} \right)} \right)^{ - 1}} = \\ {({H^T}H + {F^{ - 1}}E)^{ - 1}}, \end{array}$

由引理2.1的证明过程知矩阵HTH+F-1E是对称正定阵,所以(FHTH+E)-1F也是对称正定阵,即(FHTH+E)-1FS++n.

引理2.3 假设(A1)和(A2)成立,XFP°,(X,y,Z)∈FD°,则线性方程组(2.7)有唯一解.

证明首先证明线性方程组(2.7)的系数矩阵可逆,即证明方程组(2.7)对应的齐次线性方程组

$\left( {\begin{array}{*{20}{c}} 0&G&0\\ {{G^T}}&{ - {H^T}H}&I\\ 0&E&F \end{array}} \right)\left( \begin{array}{l} \Delta y\\ svec\left( {\Delta X} \right)\\ svec\left( {\Delta Z} \right) \end{array} \right) = \left( \begin{array}{l} 0\\ 0\\ 0 \end{array} \right),$ (2.10)

只有平凡解.记

$\left( {\begin{array}{*{20}{c}} {G{{(F{H^T}H + E)}^{ - 1}}F{G^T}}&0&0\\ {{G^T}}&{ - {F^{ - 1}}E - {H^T}H}&0\\ 0&E&F \end{array}} \right),$

利用块高斯消去法对(2.10)式进行化简,得到

$B\left( \begin{array}{l} \Delta y\\ svec\left( {\Delta X} \right)\\ svec\left( {\Delta Z} \right) \end{array} \right) = \left( \begin{array}{l} 0\\ 0\\ 0 \end{array} \right),$ (2.11)

于是有

$(G{(F{H^T}H + E)^{ - 1}}F{G^T})\Delta y = 0.$ (2.12)

由引理2.2知(FHTH+E)-1F是对称正定阵,又因为G为行满秩,所以G(FHTH+E)-1FGT也是对称正定阵,于是由(2.12)式知Δy=0.

将Δy=0代入(2.11)式,即得

$\begin{array}{*{20}{l}} {svec\left( {\Delta X} \right) = {{\left( {F{H^T}H + E} \right)}^{ - 1}}F{G^T}\Delta y = 0,}\\ {svec\left( {\Delta Z} \right) = - {F^{ - 1}}Esvec\left( {\Delta X} \right) = 0,} \end{array}$

所以(2.10)式只有平凡解,从而(2.7)式的解存在并且唯一.

引理2.4 假设(A1)和(A2)成立,设X∈FP°,(X,y,Z)∈FD°,(ΔXyZ)是线性方程组(2.7)对于某个给定的矩阵W(∈Rn×n)的解,则以下结论成立:

① ΔZ·ΔX≥0;

X·ΔZ+Z·ΔX=Tr(W);

③ 如果W=σμI-X ${\frac{1}{2}}$ ZX ${\frac{1}{2}}$ ,σR,μ= $\frac{{X \cdot Z}}{n}$ ,则

(X+αΔX)·(Z+αΔZ)=(1-α+α σ)(X·Z)+α2X·ΔZ).

证明 ① 根据(2.5b)式和(2.5a)式,并结合矩阵内积的性质可得

$\begin{array}{l} \Delta Z\cdot\Delta X = \sum\limits_{j = 1}^l {} {H_j}({H_j}\cdot\Delta X) - \sum\limits_{i = 1}^m {} \Delta {y_i}{A_i})\cdot\\ \Delta X = svec{\left( {\Delta X} \right)^T}{H^T}Hsvec\left( {\Delta X} \right) - \sum\limits_{i = 1}^m {} \Delta {y_i}({A_i}\cdot\\ \Delta X){\rm{ }} = svec{\left( {\Delta X} \right)^T}{H^T}Hsvec\left( {\Delta X} \right) \ge 0, \end{array}$

最后一个不等式成立是因为HTH是对称半正定矩阵.

② 由(2.5c)式知X- ${\frac{1}{2}}$ (XΔZXZ)X ${\frac{1}{2}}$ +X ${\frac{1}{2}}$ ZX+ZΔX)X- ${\frac{1}{2}}$ =2W,则由矩阵迹的性质有

$\begin{array}{l} 2Tr\left( W \right){\rm{ }} = Tr\left( {{X^{ - \frac{1}{2}}}\left( {X\Delta Z + \Delta XZ} \right){X^{\frac{1}{2}}}} \right) + \\ Tr\left( {{X^{\frac{1}{2}}}\left( {\Delta ZX + Z\Delta X} \right){X^{ - \frac{1}{2}}}} \right) = Tr\left( {X\Delta Z + \Delta XZ} \right) + \\ Tr\left( {\Delta ZX + Z\Delta X} \right) = 2Tr\left( {X\Delta Z + \Delta XZ} \right) = \\ 2\left( {X\cdot\Delta Z + Z\cdot\Delta X} \right), \end{array}$

X·ΔZ+Z·ΔX=Tr(W).

③ 根据矩阵内积的线性性质以及矩阵迹的性质,再结合结论(2)得到

$\begin{array}{*{20}{l}} {\left( {X + \alpha \Delta X} \right) \cdot \left( {Z + \alpha \Delta Z} \right) = X \cdot Z + \alpha \left( {X \cdot \Delta Z + Z \cdot \Delta X} \right) + }\\ {{\alpha ^2}\left( {\Delta X \cdot \Delta Z} \right) = X \cdot Z + \alpha Tr(\sigma \mu I - {X^{12}}Z{X^{12}}) + }\\ {{\alpha ^2}\left( {\Delta X \cdot \Delta Z} \right) = X \cdot Z + \alpha Tr\left( {\sigma \mu I} \right) - }\\ {\alpha Tr({X^{\frac{1}{2}}}Z{X^{\frac{1}{2}}}) + {\alpha ^2}\left( {\Delta X \cdot \Delta Z} \right),} \end{array}$

又因为Tr(σμI)=σμTr(I)=σμn,而μ= $\frac{{X \cdot Z}}{n}$ ,所以有

$\begin{array}{*{20}{l}} {\left( {X + \alpha \Delta X} \right) \cdot \left( {Z + \alpha \Delta Z} \right) = X \cdot Z + \alpha \sigma \mu n - }\\ {\alpha Tr\left( {XZ} \right) + {\alpha ^2}\left( {\Delta X \cdot \Delta Z} \right) = X \cdot Z + \alpha \sigma \mu n - }\\ {\alpha \left( {X \cdot Z} \right) + {\alpha ^2}\left( {\Delta X \cdot \Delta Z} \right) = \left( {1 - \alpha + \alpha \sigma } \right)\left( {X \cdot Z} \right) + }\\ {{\alpha ^2}\left( {\Delta X \cdot \Delta Z} \right),} \end{array}$

即结论(3)成立.

3 H..K..M搜索方向的计算

H..K..M搜索方向(ΔXyZ)可以通过求解线性方程组(2.7)得到,而方程组(2.7)包含m+n(n+1)个线性方程,通过块高斯消去法将其化简为如下的方程

$\begin{array}{l} (G{(F{H^T}H + E)^{ - 1}}F{G^T})\Delta y = - G(F{H^T}H + \\ E{)^{ - 1}}{r_c}, \end{array}$ (3.1a)
${G^T}\Delta y - ({F^{ - 1}}E + {H^T}H)svec\left( {\Delta X} \right) = - {F^{ - 1}}{r_c},$ (3.1b)
$Esvec\left( {\Delta X} \right) + Fsvec\left( {\Delta Z} \right) = {r_c},$ (3.1c)

M=G(FHTH+E)-1FGT,h=-G(FHTH+E)-1rc,则(3.1a)式可改写为

$M\Delta y = h.$ (3.2)

可证M是对称正定阵.事实上,因为

$\begin{array}{l} M = G{({F^{ - 1}}(F{H^T}H + E))^{ - 1}}{G^{T}} = G({H^T}H\\ + {F^{ - 1}}E{)^{ - 1}}{G^T} = G{({H^T}H + ({X^{ - 1}}{ \otimes _s}Z))^{ - 1}}{G^T}, \end{array}$

上述最后一个等式成立是由(2.9)式得到.由X-1,Z是对称正定阵知X-1sZ也是对称正定阵,而HTH是对称半正定阵,因此(HTH+(X-1sZ))-1是对称正定阵,故M是对称正定阵.

于是可利用Cholesky分解求解线性方程组(3.2)得到Δy,进而由(3.1b)式和(3.1c)式可求出

$\Delta X = smat({(F{H^T}H + E)^{ - 1}}(F{G^T}\Delta y + {r_c})),$ (3.3a)
$\Delta Z = smat({H^T}Hsvec\left( {\Delta X} \right) - {G^T}\Delta y).$ (3.3b)
4 短步原始对偶路径跟踪算法

这一节我们将给出基于H..K..M方向的短步原始对偶路径跟踪算法.

算法的具体步骤如下:

算法A

步骤0(初始步) 选取X0FP°,(X0,y0,Z0)∈FD°,参数η>1,令k=0,μ0= $\frac{{{X^0} \cdot {Z^0}}}{n}$ .

μk>2-ημ0时执行以下步骤:

步骤1 记X=Xk,(y,Z)=(yk,Zk),μ=μk.

步骤2 选择合适参数σ=σk∈[0,1],令W=σμI-X ${\frac{1}{2}}$ ZX ${\frac{1}{2}}$ .

步骤3 求解线性方程组(2.7)得到H..K..M搜索方向(ΔXyZ).

步骤4 选择α≥0,使得

$\begin{array}{*{20}{l}} {\hat X = X + \alpha \Delta X \in S_{ + + }^n,\hat y = y + \alpha \Delta y \in {R^m},}\\ {\hat Z = Z + \alpha \Delta Z \in S_{ + + }^n.} \end{array}$

步骤5 令Xk+1= ${\hat X}$ ,(yk+1,Zk+1)=( $\hat y,\hat Z$ ),μ=μk+1= $\frac{{{X^{k + 1}} \cdot {Z^{k + 1}}}}{n}$ .令k=k+1.

注:算法A是一个短步路径跟踪算法[5, 8],对所有的k≥0,取αk=1,σk=1-δ/ $\sqrt n $ ,其中δ>0是一个常量.

算法产生的迭代点列将落在中心路径的如下邻域内:

$\begin{array}{l} {N_F}\left( \gamma \right){\rm{ }} \equiv \{ \left( {X,y,Z} \right)|X \in F_P^ \circ ,\left( {X,y,Z} \right) \in \\ F_D^ \circ ,{\left\| {{X^{\frac{1}{2}}}Z{X^{\frac{1}{2}}} - \mu I} \right\|_F} \le \gamma \mu \} = \{ \left( {X,y,Z} \right)|X \in \\ F_P^ \circ ,\left( {X,y,Z} \right) \in F_D^ \circ ,{\rm{ }}{\left( {\sum\limits_{i = 1}^n {} {{\left( {{\lambda _i}\left( {XZ} \right) - \mu } \right)}^2}} \right)^{1/2}} \le \gamma \mu , \end{array}$ (4.1)

其中μ= $\frac{{X \cdot Z}}{n}$ ,γ∈(0,1),‖·‖F表示矩阵的Frobenius范数.

引理4.1 假设(A1)和(A2)成立,设X∈FP°,(X,y,Z)∈FD°,(ΔXZy)是线性方程组(2.7)的解,其中W=σμI-X ${\frac{1}{2}}$ ZX ${\frac{1}{2}}$ ,对任意的α≥0,令

$\left( {X\left( \alpha \right),Z\left( \alpha \right),y\left( \alpha \right)} \right) \equiv \left( {X,Z,y} \right) + \alpha \left( {\Delta X,\Delta Z,\Delta y} \right),$ (4.2a)
$\mu \left( \alpha \right) \equiv \left( {X\left( \alpha \right)\cdot{\rm{ }}Z\left( \alpha \right)} \right)/n,$ (4.2b)
$Q\left( \alpha \right) \equiv {X^{ - \frac{1}{2}}}\left( {X\left( \alpha \right)Z\left( \alpha \right) - \mu \left( \alpha \right)I} \right){X^{\frac{1}{2}}}$ (4.2c)

则有

$\begin{array}{l} Q\left( \alpha \right) + Q{\left( \alpha \right)^T} = 2\left( {1 - \alpha } \right)\left( {{X^{\frac{1}{2}}}Z{X^{\frac{1}{2}}} - \mu I} \right) - \\ \frac{{2{\alpha ^2}}}{n}\left( {\Delta X\cdot\Delta Z} \right)I + {\alpha ^2}({X^{ - \frac{1}{2}}}\Delta X\Delta Z{X^{\frac{1}{2}}}\\ + {X^{\frac{1}{2}}}\Delta Z\Delta X{X^{ - \frac{1}{2}}}). \end{array}$ (4.3)

其中In阶单位阵.

证明 假设α≥0是给定的.根据引理2.4(3)有

$\begin{array}{*{20}{l}} {X\left( \alpha \right) \cdot Z\left( \alpha \right) = \left( {X + \alpha \Delta X} \right) \cdot \left( {Z + \alpha \Delta Z} \right) = \left( {1 - \alpha + \alpha \sigma } \right)}\\ {\left( {X \cdot Z} \right) + {\alpha ^2}\left( {\Delta X \cdot \Delta Z} \right),} \end{array}$

所以由(4.2b)即知

$\mu \left( \alpha \right) = \left( {1 - \alpha + \alpha {\rm{ }}\sigma } \right)\mu + \frac{{{\alpha ^2}}}{n}\left( {\Delta X\cdot\Delta Z} \right),$

因此,

$\begin{array}{l} X\left( \alpha \right)Z\left( \alpha \right) - \mu \left( \alpha \right)I = \left( {X + \alpha \Delta X} \right)\left( {Z + \alpha \Delta Z} \right) - \\ \left( {\left( {1 - \alpha + \alpha {\rm{ }}\sigma } \right)\mu {\rm{ }} + \frac{{{\alpha ^2}}}{n}\left( {\Delta X\cdot\Delta Z} \right)} \right)I{\rm{ }} = XZ + \alpha X\Delta Z + \alpha \Delta XZ + \\ {\alpha ^2}\Delta X\Delta Z - \left( {1 - \alpha } \right)\mu I{\rm{ }} - \alpha {\rm{ }}\sigma \mu I - \frac{{{\alpha ^2}}}{n}\left( {\Delta X\cdot\Delta Z} \right)I = \\ XZ - \alpha XZ - \left( {1 - \alpha } \right)\mu I - \alpha {\rm{ }}\sigma \mu I + \alpha XZ + \\ \alpha \left( {X\Delta Z + \Delta XZ} \right) + {\alpha ^2}\Delta X\Delta Z - \frac{{{\alpha ^2}}}{n}\left( {\Delta X\cdot\Delta Z} \right)I = \\ \left( {1 - \alpha } \right)\left( {XZ - \mu I} \right) + \alpha \left( {XZ - \sigma \mu I} \right) + \alpha (X\Delta Z + \Delta XZ)\\ + {\alpha ^2}\Delta X\Delta Z - \frac{{{\alpha ^2}}}{n}\left( {\Delta X\cdot\Delta Z} \right)I \end{array}$ (4.4)

类似可以得到

$\begin{array}{l} Z\left( \alpha \right)X\left( \alpha \right) - \mu \left( \alpha \right)I = \left( {1 - \alpha } \right)\left( {ZX - \mu I} \right) + \\ \alpha \left( {ZX - \sigma \mu I} \right) + \alpha \left( {\Delta ZX + Z\Delta X} \right) + {\alpha ^2}\Delta Z\Delta X - \\ \frac{{{\alpha ^2}}}{n}\left( {\Delta X\cdot\Delta Z} \right)I. \end{array}$ (4.5)

于是,根据(4.2c)式,(4.4)式,(4.5)式得到

$\begin{array}{l} Q\left( \alpha \right) + Q{\left( \alpha \right)^T} = {X^{ - \frac{1}{2}}}\left( {X\left( \alpha \right)Z\left( \alpha \right) - \mu \left( \alpha \right)I} \right){X^{\frac{1}{2}}} + \\ {X^{\frac{1}{2}}}\left( {Z\left( \alpha \right)X\left( \alpha \right) - \mu \left( \alpha \right)I} \right){X^{ - \frac{1}{2}}} = 2\left( {1 - \alpha } \right)\left( {{X^{\frac{1}{2}}}Z{X^{\frac{1}{2}}} - \mu I} \right)\\ + 2\alpha \left( {{X^{\frac{1}{2}}}Z{X^{\frac{1}{2}}} - \sigma \mu I} \right){\rm{ }} + \alpha ({X^{ - \frac{1}{2}}}(X\Delta Z + \Delta XZ){X^{\frac{1}{2}}}\\ + {X^{\frac{1}{2}}}\left( {\Delta ZX + Z\Delta {\rm{ }}X} \right){X^{ - \frac{1}{2}}}) + \\ {\alpha ^2}({X^{ - \frac{1}{2}}}\Delta X\Delta Z{X^{\frac{1}{2}}} + {X^{\frac{1}{2}}}\Delta Z\Delta X{X^{ - \frac{1}{2}}}) - \\ {X^{ - \frac{1}{2}}}(\frac{{{\alpha ^2}}}{n}\cdot\left( {\Delta X\cdot\Delta Z} \right)I){X^{\frac{1}{2}}} - {X^{\frac{1}{2}}}(\frac{{{\alpha ^2}}}{n}\left( {\Delta X\cdot\Delta Z} \right)I){X^{ - \frac{1}{2}}} = \\ 2\left( {1 - \alpha } \right)({X^{\frac{1}{2}}}Z{X^{\frac{1}{2}}} - \mu I) + \\ {\alpha ^2}({X^{ - \frac{1}{2}}}\Delta X\Delta Z{X^{\frac{1}{2}}} + {X^{\frac{1}{2}}}\Delta Z\Delta X{X^{ - \frac{1}{2}}}) - \frac{{2{\alpha ^2}}}{n}\left( {\Delta {\rm{ }}X\cdot\Delta {\rm{ }}Z} \right)I, \end{array}$

最后一个等式成立是根据(2.5c)式.

接下来的引理给出scaling方向X- ${\frac{1}{2}}$ Δ XX- ${\frac{1}{2}}$ ,X ${\frac{1}{2}}$ Δ ZX ${\frac{1}{2}}$ 的界,对任意(X,y,Z)∈FP°×FD°都成立.证明过程需要使用如下结论,对任意的A,BRn×n,都有(详见文献[15]5.6节习题20).

$\begin{array}{l} {\left\| {A + B} \right\|_F} \le {\left\| A \right\|_F} + {\left\| B \right\|_F},{\left\| {AB} \right\|_F} \le \\ {\left\| A \right\|_F}{\left\| B \right\|_F},{\rm{ }}{\left\| {AB} \right\|_F} \le {\left\| A \right\|_F}\left\| B \right\| \end{array}$ (4.6)

引理4.2 假设(A1)和(A2)成立,设X∈FP°,(X,y,Z)∈FD°,使得‖X ${\frac{1}{2}}$ ZX ${\frac{1}{2}}$ -υI‖≤υγ,对任意的γ∈[0,1),υ>0成立.(ΔXyZ)是线性方程组(2.7)对于某个确定矩阵W的解.记δx=υX- ${\frac{1}{2}}$ ΔXX- ${\frac{1}{2}}$ F,δz=υX ${\frac{1}{2}}$ ΔZX ${\frac{1}{2}}$ F,则下列结论成立:

${\delta _x}{\delta _z} \le \frac{1}{2}(\delta _x^2 + \delta _z^2) \le {\rm{ }}\frac{{\left\| W \right\|_F^2}}{{2{{\left( {1 - \gamma } \right)}^2}}}.$ (4.7)

证明 根据(2.5c)式得

$\begin{array}{l} W = \frac{1}{2}({X^{ - \frac{1}{2}}}\left( {X\Delta Z + \Delta XZ} \right){X^{\frac{1}{2}}} + {X^{\frac{1}{2}}}(\Delta ZX + Z\Delta X)\\ {X^{ - 12}}){\rm{ }} = \frac{1}{2}{X^{\frac{1}{2}}}\Delta Z{X^{\frac{1}{2}}} + \frac{1}{2}{X^{ - \frac{1}{2}}}\Delta XZ{X^{\frac{1}{2}}} + \\ 12{X^{\frac{1}{2}}}\Delta Z{X^{\frac{1}{2}}} + \frac{1}{2}{X^{\frac{1}{2}}}Z\Delta X{X^{ - \frac{1}{2}}} = {X^{\frac{1}{2}}}\Delta Z{X^{\frac{1}{2}}} + \\ \upsilon {X^{ - \frac{1}{2}}}\Delta X{X^{ - \frac{1}{2}}} - \frac{1}{2}\upsilon {X^{ - \frac{1}{2}}}\Delta X{X^{ - \frac{1}{2}}} + \\ \frac{1}{2}{X^{ - \frac{1}{2}}}\Delta X{X^{ - \frac{1}{2}}}{X^{\frac{1}{2}}}Z{X^{\frac{1}{2}}} + \frac{1}{2}{X^{\frac{1}{2}}}Z{X^{\frac{1}{2}}}{X^{ - \frac{1}{2}}}\Delta X{X^{ - \frac{1}{2}}} - \\ \frac{1}{2}\upsilon {X^{ - \frac{1}{2}}}\Delta X{X^{ - \frac{1}{2}}} = {X^{\frac{1}{2}}}\Delta Z{X^{\frac{1}{2}}} + \upsilon {X^{ - \frac{1}{2}}}\Delta X{X^{ - \frac{1}{2}}} + \\ \frac{1}{2}{X^{ - \frac{1}{2}}}\Delta X{X^{ - \frac{1}{2}}}({X^{\frac{1}{2}}}Z{X^{\frac{1}{2}}} - \upsilon I) + \frac{1}{2}({X^{\frac{1}{2}}}Z{X^{\frac{1}{2}}} - \\ \upsilon I){X^{ - \frac{1}{2}}}\Delta X{X^{ - \frac{1}{2}}}, \end{array}$

由引理2.4(1)知ΔZ·ΔX≥0,再根据(4.6)式以及引理4.2的已知条件得到

$\begin{array}{l} {\left\| W \right\|_F} \ge {\left\| {{X^{\frac{1}{2}}}\Delta Z{X^{\frac{1}{2}}} + \upsilon {X^{ - \frac{1}{2}}}\Delta X{X^{ - \frac{1}{2}}}} \right\|_F}\\ - \frac{1}{2}\left\| {{X^{ - \frac{1}{2}}}\Delta X{X^{ - \frac{1}{2}}}\left( {{X^{\frac{1}{2}}}Z{X^{\frac{1}{2}}} - \upsilon I} \right) + } \right.\\ {\left. {({X^{\frac{1}{2}}}Z{X^{\frac{1}{2}}} - \upsilon I){X^{ - \frac{1}{2}}}\Delta X{X^{ - \frac{1}{2}}}} \right\|_F} \ge \left\| {{X^{\frac{1}{2}}}\Delta Z{X^{\frac{1}{2}}} + } \right.\\ {\left. {\upsilon {X^{ - \frac{1}{2}}}\Delta X{X^{ - \frac{1}{2}}}} \right\|_F} - \left\| {{X^{\frac{1}{2}}}Z{X^{\frac{1}{2}}} - {\rm{ }}\upsilon I} \right\|\cdot\\ {\left\| {{X^{ - \frac{1}{2}}}\Delta X{X^{ - \frac{1}{2}}}} \right\|_F} \ge ({\left\| {{X^{\frac{1}{2}}}\Delta Z{X^{\frac{1}{2}}}} \right\|^2}_F + \\ {\upsilon ^2}{\left\| {{X^{ - \frac{1}{2}}}\Delta X{X^{ - \frac{1}{2}}}} \right\|^2}_F{)^{\frac{1}{2}}} - \upsilon \gamma \cdot\frac{{{\delta _x}}}{\upsilon } = \sqrt {\delta _x^2 + \delta _z^2} - \\ \gamma \cdot{\delta _x} \ge \left( {1 - \gamma } \right)\sqrt {\delta _x^2 + \delta _z^2} , \end{array}$

所以

$\frac{{\left\| W \right\|_F^2}}{{{{\left( {1 - \gamma } \right)}^2}}} \ge \delta _x^2 + \delta _z^2.$

注意到δx2+δz2≥2δxδz,于是(4.7)式成立.

以下定理说明本文的短步路径跟踪算法产生的迭代点列将落在中心路径邻域NF(γ)内.

定理4.1 假设(A1)和(A2)成立,设常数γ∈(0,1),δ∈[0, $\sqrt n $ ),并且满足

$\frac{{{\gamma ^2} + {\delta ^2}}}{{2{{\left( {1 - \gamma } \right)}^2}\left( {1 - \delta /\sqrt n } \right)}} \le \gamma ,\gamma \le \frac{1}{2}.$ (4.8)

又设(X,y,Z)∈NF(γ),(ΔXyZ)是线性方程组(2.7)当W=σμI-X ${\frac{1}{2}}$ ZX ${\frac{1}{2}}$ 时的解,σ=1-δ/ $\sqrt n $ .则

( $\hat X,\hat y,\hat Z$ )=(X+ΔX,yy,ZZ)∈NF(γ).

证明 因为

$\begin{array}{*{20}{l}} {({X^{\frac{1}{2}}}Z{X^{\frac{1}{2}}} - \mu I) \cdot I = Tr(({X^{\frac{1}{2}}}Z{X^{\frac{1}{2}}} - \mu I)I) = }\\ {Tr({X^{\frac{1}{2}}}Z{X^{\frac{1}{2}}} - \mu I) = Tr({X^{\frac{1}{2}}}Z{X^{\frac{1}{2}}}) - }\\ {n\mu = Tr\left( {XZ} \right) - n\mu = X \cdot Z - n\mu = 0,} \end{array}$

并结合σ=1-δ/ $\sqrt n $ ,和NF(γ)的定义,则有

$\begin{array}{l} \left\| W \right\|_F^2 = \left\| {\sigma \mu I - {X^{\frac{1}{2}}}Z{X^{\frac{1}{2}}}} \right\|_F^2 = \left\| {\left( {\sigma - 1} \right)\mu I + } \right.\\ \left. {\mu I - {X^{\frac{1}{2}}}Z{X^{\frac{1}{2}}}} \right\|_F^2 = \left\| {\left( {\sigma - 1} \right)\mu I} \right\|_F^2 + \left\| {\mu I} \right.\\ \left. { - {X^{\frac{1}{2}}}Z{X^{\frac{1}{2}}}} \right\|_F^2 \le {\left( {1 - \sigma } \right)^2}{\mu ^2}\left\| I \right\|_F^2 + {\gamma ^2}{\mu ^2} = \\ ({\delta ^2} + {\gamma ^2}){\mu ^2}. \end{array}$ (4.9)

因为‖X ${\frac{1}{2}}$ ZX ${\frac{1}{2}}$ -μIFγμ,根据引理4.2取υ=μ,W=σμI-X ${\frac{1}{2}}$ ZX ${\frac{1}{2}}$ ,则有

$\begin{array}{*{20}{l}} {\delta _x^2 \le \delta _x^2 + \delta _z^2 \le \frac{{\left\| {} \right\|W_F^2}}{{{{\left( {1 - \gamma } \right)}^2}}},}\\ {{\delta _x}{\delta _z} \le \frac{1}{2}\left( {\delta _x^2 + \delta _z^2} \right) \le \frac{{\left\| W \right\|_F^2}}{{2{{\left( {1 - \gamma } \right)}^2}}},} \end{array}$

从而

$\begin{array}{l} {\left\| {{X^{ - \frac{1}{2}}}\Delta X{X^{ - \frac{1}{2}}}} \right\|_F} = \frac{{{\delta _x}}}{\upsilon } \le \frac{{{{\left\| W \right\|}_F}}}{{\left( {1 - \gamma } \right)\upsilon }} = \\ \frac{{{{\left\| {\sigma \mu I - {X^{\frac{1}{2}}}Z{X^{\frac{1}{2}}}} \right\|}_F}}}{{\left( {1 - \gamma } \right)\mu }}, \end{array}$ (4.10)
$\begin{array}{l} {\left\| {{X^{ - \frac{1}{2}}}\Delta X{X^{ - \frac{1}{2}}}} \right\|_F}{\left\| {{X^{\frac{1}{2}}}\Delta Z{X^{\frac{1}{2}}}} \right\|_F} = \frac{{{\delta _x}\cdot{\delta _z}}}{\upsilon } \le \\ \left\| W \right\|_F^22\upsilon {\left( {1 - \gamma } \right)^2} = \frac{{\left\| {\sigma \mu {\rm{ }}I - {X^{\frac{1}{2}}}Z{X^{\frac{1}{2}}}} \right\|_F^2}}{{2\mu {{\left( {1 - \gamma } \right)}^2}}}. \end{array}$ (4.11)

根据(4.2b)式,(4.2c)式和(4.3)式得

$\begin{array}{*{20}{l}} {\hat Q = Q\left( 1 \right) = {X^{ - \frac{1}{2}}}\left( { - I} \right){X^{\frac{1}{2}}},\hat \mu = \frac{{\hat X\cdot\hat Z}}{n},}\\ \begin{array}{l} Q\left( 1 \right) + Q{\left( 1 \right)^T} = ({X^{ - \frac{1}{2}}}\Delta X\Delta Z{X^{\frac{1}{2}}} + \\ {X^{\frac{1}{2}}}\Delta Z\Delta X{X^{ - \frac{1}{2}}}) - \frac{2}{n}\left( {\Delta X\cdot\Delta Z} \right)I. \end{array} \end{array}$

再根据(4.6)式,有

$\begin{array}{l} \frac{1}{2}{\left\| {\hat Q + {{\hat Q}^T}} \right\|_F} = \frac{1}{2}\left\| {{X^{ - \frac{1}{2}}}\Delta X\Delta Z{X^{\frac{1}{2}}} - } \right.\\ \frac{1}{n}\left( {\Delta X\cdot\Delta Z} \right)I + {X^{\frac{1}{2}}}\Delta {\rm{ }}Z\Delta X{X^{ - \frac{1}{2}}} - \frac{1}{n}\\ {\left. {\left( {\Delta X\cdot\Delta Z} \right)I} \right\|_F} \le {\left\| {{X^{ - \frac{1}{2}}}\Delta X\Delta Z{X^{\frac{1}{2}}} - 1n\left( {\Delta X\cdot\Delta Z} \right)I} \right\|_F}, \end{array}$ (4.12)

由矩阵Frobenious范数的定义知

$\begin{array}{l} {\left\| {{X^{ - \frac{1}{2}}}\Delta X\Delta Z{X^{\frac{1}{2}}} - 1n\left( {\Delta X\cdot\Delta Z} \right)I} \right\|^2}_F = \\ Tr({X^{\frac{1}{2}}}\Delta Z\Delta X{X^{ - \frac{1}{2}}}{X^{ - \frac{1}{2}}}\Delta X\Delta Z{X^{\frac{1}{2}}}) - \frac{2}{n}(\Delta X\cdot\\ \Delta Z)Tr({X^{\frac{1}{2}}}\Delta Z\Delta X{X^{ - \frac{1}{2}}}) + \frac{1}{{{n^2}}}{\left( {\Delta X\cdot\Delta Z} \right)^2}Tr\left( I \right) = \\ Tr({X^{\frac{1}{2}}}\Delta Z\Delta X{X^{ - \frac{1}{2}}}{X^{ - \frac{1}{2}}}\Delta X\Delta Z{X^{\frac{1}{2}}}) - \frac{2}{n}(\Delta X\cdot\\ \Delta Z)\left( {\Delta X\cdot\Delta Z} \right) + \frac{1}{{{n^2}}}{\left( {\Delta X\cdot\Delta Z} \right)^2}n = \\ Tr({X^{\frac{1}{2}}}\Delta Z\Delta X{X^{ - \frac{1}{2}}}{X^{ - \frac{1}{2}}}\Delta X\Delta Z{X^{\frac{1}{2}}}) - \\ 1n{\left( {\Delta X\cdot\Delta Z} \right)^2} \le \\ Tr({X^{\frac{1}{2}}}\Delta Z\Delta X{X^{ - \frac{1}{2}}}{X^{ - \frac{1}{2}}}\Delta X\Delta Z{X^{\frac{1}{2}}}) = \\ ({X^{ - \frac{1}{2}}}\Delta X\Delta Z{X^{\frac{1}{2}}})\cdot({X^{ - \frac{1}{2}}}\Delta X\Delta Z{X^{\frac{1}{2}}}) = \\ \left\| {{X^{ - \frac{1}{2}}}\Delta X\Delta Z{X^{\frac{1}{2}}}} \right\|_F^2, \end{array}$

$\begin{array}{l} {\left\| {{X^{ - \frac{1}{2}}}\Delta X\Delta Z{X^{\frac{1}{2}}} - \frac{1}{n}\left( {\Delta X\cdot\Delta Z} \right)I} \right\|_F} \le \\ {\left\| {{X^{ - \frac{1}{2}}}\Delta X\Delta Z{X^{\frac{1}{2}}}} \right\|_F}. \end{array}$

将上式代入(4.12)式,再结合(4.11)式,(4.9)式和(4.8)式,得

$\begin{array}{l} \frac{1}{2}{\left\| {\hat Q + {{\hat Q}^T}} \right\|_F} \le {\left\| {{X^{ - \frac{1}{2}}}\Delta X\Delta Z{X^{\frac{1}{2}}}} \right\|_F} = \\ {\left\| {{X^{ - \frac{1}{2}}}\Delta X{X^{ - \frac{1}{2}}}{X^{\frac{1}{2}}}\Delta Z{X^{\frac{1}{2}}}} \right\|_F} \le \\ {\left\| {{X^{ - \frac{1}{2}}}\Delta X{X^{ - \frac{1}{2}}}} \right\|_F}{\left\| {{X^{\frac{1}{2}}}\Delta Z{X^{\frac{1}{2}}}} \right\|_F} \le \\ \frac{{\left\| {\sigma \mu I - {X^{\frac{1}{2}}}Z{X^{\frac{1}{2}}}} \right\|_F^2}}{{2\mu {{\left( {1 - \gamma } \right)}^2}}} \le \frac{{({\delta ^2} + {\gamma ^2})\mu }}{{2{{\left( {1 - \gamma } \right)}^2}}} \le \gamma (1 - \\ \delta /\sqrt n )\mu \le \gamma \hat \mu , \end{array}$ (4.13)

其中 $\hat \mu = \frac{{\hat X \cdot \hat Z}}{n} \ge \left( {1 - \delta /\sqrt n } \right)\frac{{X \cdot Z}}{n} = \left( {1 - \delta /\sqrt n } \right)\mu $ .

根据(4.10)式,(4.9)式和定理4.1的条件,得

$\begin{array}{l} {\left\| {{X^{ - \frac{1}{2}}}\Delta X{X^{ - \frac{1}{2}}}} \right\|_F} \le \frac{{{{\left\| {\sigma \mu I - {X^{\frac{1}{2}}}Z{X^{\frac{1}{2}}}} \right\|}_F}}}{{\left( {1 - \gamma } \right)\mu }} \le \\ \frac{{{{({\delta ^2} + {\gamma ^2})}^{\frac{1}{2}}}\mu }}{{\left( {1 - \gamma } \right)\mu }} \le {\left( {2\gamma \left( {1 - \delta /\sqrt n } \right)} \right)^{\frac{1}{2}}} < 1, \end{array}$

从而知I+X- ${\frac{1}{2}}$ ΔXX- ${\frac{1}{2}}$ 是正定矩阵.注意到X ${\frac{1}{2}}$ =(X ${\frac{1}{2}}$ )T,于是=XX=X ${\frac{1}{2}}$ (I+X- ${\frac{1}{2}}$ ΔXX- ${\frac{1}{2}}$ )X ${\frac{1}{2}}$ 也是正定矩阵.而

$\begin{array}{l} \hat Q = {X^{ - 12}}\left( {\hat X\hat Z - \hat \mu I} \right){X^{\frac{1}{2}}} = ({X^{ - \frac{1}{2}}}{{\hat X}^{\frac{1}{2}}})({{\hat X}^{\frac{1}{2}}}\hat Z{{\hat X}^{\frac{1}{2}}} - \\ \hat \mu I){({X^{ - \frac{1}{2}}}{{\hat X}^{\frac{1}{2}}})^{ - 1}}, \end{array}$

$D = {X^{ - \frac{1}{2}}}{{\hat X}^{\frac{1}{2}}},G = {{\hat X}^{\frac{1}{2}}}\hat Z{{\hat X}^{\frac{1}{2}}} - \hat \mu I$ ,则 ${\hat Q}$ =DGD-1.由文献[5]的引理3.3以及(4.13)式可得

${\left\| {{{\hat X}^{\frac{1}{2}}}\hat Z{{\hat X}^{\frac{1}{2}}} - \hat \mu I} \right\|_F} \le \frac{1}{2} + {\left\| {\hat Q + {{\hat Q}^T}} \right\|_F} \le \gamma \hat \mu ,$ (4.14)

于是λmin( ${{{\hat X}^{\frac{1}{2}}}\hat Z{{\hat X}^{\frac{1}{2}}}}$ )≥(1-γ) ${\hat \mu }$ >0,所以 ${{{\hat X}^{\frac{1}{2}}}\hat Z{{\hat X}^{\frac{1}{2}}}}$ 是正定矩阵.注意到 ${{\hat X}^{\frac{1}{2}}} = {\left( {{{\hat X}^{\frac{1}{2}}}} \right)^T}$ ,从而知 ${\hat Z}$ 也是正定矩阵.根据(2.3a)式,(2.3b)式以及(X,y,Z)∈NF(γ)可得到( $\hat X,\hat y,\hat Z$ )∈FD°且∈FP°,再结合(4.14)式,即知( $\hat X,\hat y,\hat Z$ )∈NF(γ).

5 数值试验

采用Matlab(R2011b)数学软件编程对算法A进行数值试验,程序运行环境为Windows7(64bite),Intel(R)Core(TM)i3-2330M CPU @ 2.20GHZ,RAM:4G.

考虑如下两类测试问题:(1)随机二次半定规划问题;矩阵Ai,Hj,C都是随机生成的对阵矩阵,初始点(X,y,Z)的选取需在FP,FD内.记该类测试问题为RP.(2)最近相关矩阵Nearest correlation matrix(简记为NCM)问题[16-17].

$\begin{array}{l} min\frac{1}{2}{\left\| {X - G} \right\|^2}\\ \begin{array}{*{20}{l}} {s.t.{\rm{ }}diag\left( X \right) = e,}\\ {X \in S_ + ^n,} \end{array} \end{array}$

其中G是对称正定阵,对角线上的元素为1,对角线外的元素在-1到1之间.e是所有元素都为1的n维向量.

NCM问题具有广泛的应用,例如在市场营销以及经济学方面.但是遗憾的是,由于缺乏数据等信息,不能得到一个完整的矩阵G,即矩阵G的一些元素是未知的.于是通过求解NCM问题得到一个有效的、且与G最近的相关矩阵.在数值测试中,取G是随机矩阵且G:=(1-α)B+α E[18],其中α∈(0,1),E是元素属于[-1,1]的随机对称阵,B是具有指定特征值的测试矩阵,在NCM问题的数值测试中取n=10,20,30.

利用对称Kronecker积的矩阵形式

$A{ \otimes _s}B = \frac{1}{2}Q\left( {A \otimes B + B \otimes A} \right){Q^T}.$

将矩阵的对称Kronecker积转化为矩阵的Kronecker积计算,其中Q的取法详见文献[1]附录.

在数值测试中,取算法A中的参数δ=0.3,γ=0.3,终止准则是Zk· Xk≤10-6.

测试的数值结果见表 1(所有结果都是每个维数测试10次,然后取平均值得到的),其中表中符号含义如下: prob,测试问题的编号;n,矩阵X的维数;m,不等式约束的个数; Itr,算法迭代次数;Nf,目标函数值的计算次数;time(s),算法迭代所需CPU时间(s).

表 1 数值结果 Table 1 Numerical results

上述数值结果表明本文提出的基于H..K..M方向的路径跟踪算法是有效的,能在合理的时间内求解中等规模的NCM问题.为进一步提高本文算法的数值效果,我们将在后续的研究中深入探讨步长α的更有效的确定方法,即线搜索技术.

参考文献
[1]
TODD M J, TOH K C, TVTVNCV R H. On the nesterov-todd direction in SDP[J]. SIAM Journal on Optimization, 1998, 8(3): 769-796. DOI:10.1137/S105262349630060X
[2]
ALIZADEH F, HAEBERLY J P, OVERTON M L. Primal-dual interior-point methods for semidefinite programming:Convergence rates,stability and numerical results[J]. SIAM Journal on Optimization, 1998, 8: 746-768. DOI:10.1137/S1052623496304700
[3]
KOJIMA M, SHINDOH S, HARA S. Interior-point me-thods for the monotone semidefinite linear complementarity problem in symmetric matrices[J]. SIAM Journal on Optimization, 1997, 7(1): 86-125. DOI:10.1137/S1052623494269035
[4]
HELMBERG C, RENDL F, VANDERBEI R J, et al. An interior-point methods for semidefinite programming[J]. SIAM Journal on Optimization, 1996, 6(2): 342-361. DOI:10.1137/0806020
[5]
MONTEIRO R D C. Primal-dual path-following algorithms for semidefinite programming[J]. SIAM Journal on Optimization, 1997, 7(3): 663-678. DOI:10.1137/S1052623495293056
[6]
NESTEROV Y E, TODD M J. Primal-dual interiorpoin-t methods for self-scaled cones[J]. SIAM Journal on Optimization, 1998, 8(2): 324-364. DOI:10.1137/S1052623495290209
[7]
NESTEROV Y E, TODD M J. Self-scaled barriers and interior-point methods for convex programming[J]. Mathematics of Operations Research, 1997, 22(1): 1-42. DOI:10.1287/moor.22.1.1
[8]
STURM J F, ZHANG S Z. Symmetric primal-dual path-following algorithms for semidefinite programming[J]. Applied Numerical Mathematics, 1999, 29(3): 301-315. DOI:10.1016/S0168-9274(98)00099-3
[9]
ZHANG Y. On extending some primal-dual interior-point algorithms from linear programming to semidefinite programming[J]. SIAM Journal on Optimization, 1998, 8(2): 365-386. DOI:10.1137/S1052623495296115
[10]
NESTEROV Y, NEMIROVSKII A S. Interior-Point Polynomial Methods in Convex Programming[M]. Philadelphia PA: SIAM, 1994.
[11]
NIE J W, YUAN Y X. A potential reduction algorithm for an extended SDP problem[J]. Science in China Series A Mathematics, 2000, 43(1): 35-46. DOI:10.1007/BF02903846
[12]
NIE J W, YUAN Y X. A Predictor-corrector algorith-m for QSDP combining dikin-type and newton centering steps[J]. Annals of Operations Research, 2001, 103.
[13]
关秀翠, 刁在筠. 二次半定规划问题及其投影收缩算法[J]. 全国高等学校计算数学学报, 2002, 24(2): 97-108.
GUAN X C, SIAO Z Y. The quadratic semi-definite programming problem and its projection and contraction algorithm[J]. Numerical Mathematics a Journal of Chinese Universities, 2002, 24(2): 97-108.
[14]
徐凤敏, 徐成贤. 求解二次半定规划的原对偶内点算法[J]. 工程数学学报, 2006, 23(4): 590-598.
XU F M, XU C X. Primal-dual algorithm for quadratic semi-definite programming[J]. Chinese Journal of Engineering Mathematics, 2006, 23(4): 590-598.
[15]
HORN R A, JOHNSON C R. Matrix Analysis[M]. New York: Cambridge University Press, 1985.
[16]
HIGHAM N J. Computing the nearest correlation matrixa problem from finance[J]. IMA Journal of Numerical Analysis, 2002, 22(3): 329-343. DOI:10.1093/imanum/22.3.329
[17]
TOH K C. An inexact primal-dual path following algorithm for convex quadratic SDP[J]. Mathematical Programming, 2008, 112(1): 221-254.
[18]
ZHAO X Y.A Semismooth Newton-CG Augmented Lagrangian Method for Large Scale Linear and Convex Quadratic SDPs[D].Singapore:National University of Singapore,2009.