半定规划(SDP)的算法包括内点算法和非内点算法,从初始点是否可行的角度,内点算法分为可行算法和不可行算法.它们之间存在这样一个矛盾: 可行内点算法理论复杂度小,但是实验效果不好; 不可行内点算法实验效果好,但是理论复杂度很大.目前,求解半定规划问题的可行内点算法已经很成熟,所以,不可行内点算法的研究越来越广泛,目的是在较好的实验效果基础上降低其理论复杂度.目前SDP不可行内点算法最好的迭代复杂度是O(n^1.5logL)^[1],其中n是SDP问题的维数,L=Tr(X⁰S⁰)/ε，ε是需要的精度,(X⁰,S⁰)是迭代起始点,最好的实验结果已在文献^[2]中给出.直接的算法已经很难取得更好的复杂度和实验结果,本文采用一个HMCP模型作为一个间接桥梁,提出一个新的不可行内点算法,并且采用新的宽邻域,比较起其他算法,此算法可以将复杂度降低到O(nlogL)， 并且实验效果比文献^[2]中的实验结果更好,另外一个优势是可以判别原问题是否是可行问题.所以,此算法在内点算法中是有效的.

Andersan和Ye^[3]研究线性规划的齐次算法,本文考虑半定规划的齐次算法.

考虑半定规划原问题和对偶问题的标准形式:

minTr(CX)

(P) s.t.Tr(A_iX)=b_i，i=1,…m，

X≥0，

maxb^Ty

$\left( D \right) s.t.\sum\limits_{i=1}^{m}{{{y}_{i}}}{{A}_{i}}+S=C$

S≥0.

其最优性条件是一个单调互补问题(MCP):

Tr(A_iX)=b_i，X≥0，i=1,…m，

(MCP) $\sum\limits_{i=1}^{m}{{{y}_{i}}}{{A}_{i}}+S=C$，S≥0，

XS=0.

以上问题等价于互补问题:

minTr(XS)

(MCP) s.t.Tr(A_iX)=b_i，X≥0，i=1,…m，

$\sum\limits_{i=1}^{m}{{{y}_{i}}}{{A}_{i}}+S=C$，S≥0.

令A=(vec(A₁),…,vec(A_m))^T，

h(y,vec(X))=Avec(X)－b，

g(y,vec(X))=vec(C)－A^Ty，

$f\left( y,vec\left( X \right) \right)=\left[ \begin{matrix} h\left( y,vec\left( X \right) \right) \\ g\left( y,vec\left( X \right) \right) \\ \end{matrix} \right].$

所以f(y,vec(X))是从R^m×Sⁿ₊→R^m×Rⁿ²的连续单调映射.

注: 1) vec(M)表示矩阵M的各列依次排列成的列向量; R^m表示m维欧式空间; Sⁿ₊表示n阶实对称半正定矩阵集; Sⁿ表示n阶实对称矩阵集.

2) 若∇f是f的雅克比矩阵,则∇f是半正定矩阵^[4].

那么，问题MCP等价于矩阵形式:

$\begin{align} & min vec{{\left( X \right)}^{T}}vec\left( S \right) \\ & s.t.\left[ \begin{matrix} 0 \\ vec\left( S \right) \\ \end{matrix} \right]=\left[ \begin{matrix} Avec\left( X \right)-b \\ vec\left( C \right)-{{A}^{T}}y \\ \end{matrix} \right]=f\left( y,vec\left( X \right) \right), \\ & S\ge 0,X\ge 0. \\ \end{align}$

本文将提出一个齐次模型,在新的宽邻域上来解决单调互补问题MCP. 据分析,比较起其他SDP不可行内点算法此算法具有以下优势:

1) 使迭代复杂度降低到O(nlogL).

2) 采用新的宽邻域,与其他SDP不可行算法相比,此算法具有较好的实验效果.

3) 判断原问题是否可行: 如果问题MCP有解,那么算法将会求得原问题最优解; 如果问题MCP无解,那么可以通过算法和齐次模型来证明.

本文中符号“⊗”表示2个矩阵的克罗内克积,即若A是一个m×n阶的矩阵,B是一个p×q阶的矩阵,则克罗内克积A⊗B是一个mp×nq的分块矩阵,表示为:

$A\otimes B=\left[ \begin{matrix} {{a}_{11}}B & \cdots & {{a}_{1n}}B \\ \vdots & \ddots & \vdots \\ {{a}_{m1}}B & \cdots & {{a}_{mn}}B \\ \end{matrix} \right].$

1 齐次MCP模型HMCP

提出问题MCP的齐次模型(HMCP):

minTr(XS)+τκ

(HMCP) s.t.∑mi=1y_iA_i+S=τC，S≥0，τ≥0，

Tr(A_iX)=τb_i，X≥0，i=1,…m，

－Tr(CX)+b^Ty－κ=0，κ≥0.

类似于问题MCP等价形式的分析,有

$\begin{align} & \tau f(y/\tau ,vec\left( X \right)/\tau )=\left[ \begin{matrix} \tau h(y/\tau ,vec\left( X \right)/\tau ) \\ \tau g(y/\tau ,vec\left( X \right)/\tau ) \\ \end{matrix} \right] \\ & =\left[ \begin{matrix} Avec\left( X \right)-\tau b \\ \tau vec\left( C \right)-{{A}^{T}}y \\ \end{matrix} \right]; \\ & -(y,vec\left( X \right))f(y/\tau ,~vec\left( X \right)/\tau ) \\ & =-{{y}^{T}}h(y/\tau ,vec\left( X \right)/\tau )\text{ }- \\ & vec{{\left( X \right)}^{T}}g(y/\tau ,vec\left( X \right)/\tau ) \\ & =-vec{{\left( C \right)}^{T}}vec\left( X \right)+{{b}^{T}}y, \\ \end{align}$

则问题HMCP等价于矩阵形式:

$\begin{align} & min~vec{{\left( X \right)}^{T}}vec\left( S \right)+\tau \kappa \\ & s.t.\left[ \begin{matrix} 0 \\ vec\left( S \right) \\ \kappa \\ \end{matrix} \right]=\left[ \begin{matrix} \tau f(y/\tau ,vec\left( X \right)/\tau ) \\ -(y,vec\left( X \right))f(y/\tau ,vec\left( X \right)/\tau ) \\ \end{matrix} \right] \\ & S\ge 0,\tau \ge 0,X\ge 0,\kappa \ge 0. \\ \end{align}$

定义

$\psi (y,vec\left( X \right),\tau ):=\left[ \begin{align} & \tau f(y/\tau ,vec\left( X \right)/\tau ) \\ & -(y,vec\left( X \right))f(y/\tau ,vec\left( X \right)/\tau ) \\ \end{align} \right]$

(1)

则(0,vec(S),κ)^T=ψ(y,vec(X),τ).那么可以得到

$\begin{align} & \Delta \psi = \\ & \left[ \begin{matrix} \Delta f & f-\frac{1}{\tau } \Delta f \\ {{f}^{T}}-\frac{1}{\tau } (y,vec\left( X \right)) \Delta f & -\frac{1}{\tau }(y/\tau ,vec\left( X \right)/\tau ) \Delta f \\ \end{matrix} \right] \\ \end{align}$

(2)

其中,f和ψ分别是f(y/τ,vec(X)/τ)和ψ(y,vec(X),τ)的简写形式.

引理1.1 定义ψ如(1)所示,则如下结论成立:

(a) (y,vec(X),τ)ψ(y,vec(X),τ)=0.

(b) Δψ(y,vec(X),τ)(y,vec(X),τ)^T=ψ(y,vec(X),τ).

(c) (y,vec(X),τ)Δψ(y,vec(X),τ)=－ψ(y,vec(X),τ)^T.

(d)Δf是定义在R^m×Rⁿ²上的半正定矩阵,Δψ是定义在R^m×Rⁿ²×R上的半正定矩阵.

证明 由ψ的定义和(2)易得(a),(b),(c).

因为f是从R^m×Sⁿ₊→R^m×Rⁿ²的连续单调映射,所以∇f是半正定矩阵^[4]. 将下式展开可得:

$\begin{align} & (dy,vec(dX),d\tau ) \nabla \psi (dy,vec(dX),d\tau ) \\ & {{(dy,vec(dX),d\tau )}^{T}} \\ & =[(dy,vec(dX))-d\tau (y,vec\left( X \right))/\tau ]\cdot \\ & \nabla f{{[(dy,vec(dX))-d\tau (y,vec\left( X \right))/\tau ]}^{T}} \\ & \ge 0 \\ \end{align}$

则∇ψ是半正定矩阵. (d)得证.□

定理1.1

(a) 若f是从R^m+n2到R^m+n2的连续单调映射,则ψ是从R^m+n2+1到R^m+n2+1的连续单调映射.

(b) 问题HMCP是渐进可行的,并且每一个渐进可行点是一个渐进互补解.

(c) 设问题HMCP的极大互补解是(vec(X^*),y^*,vec(S^*),τ^*,κ^*)，问题MCP有最优解当且仅当τ^*>0成立,这时,(vec(X^*)/τ^*,y^*/τ^*,vec(S^*)/τ^*)是问题MCP的一个互补解.

(d) 当κ^*>0时,问题MCP是强不可行的.

证明

(a) 任取(y₁,vec(X₁),τ₁)，(y₂,vec(X₂),τ₂)∈R^m+n2+1，因为 (y₁/τ₁,vec(X₁)/τ₁)，(y₂/τ₂,vec(X₂)/τ₂)∈R^m+n2， 由引理1.1得

$\begin{align} & ~({{y}_{1}}-{{y}_{2}},vec({{X}_{1}})-vec({{X}_{2}}),{{\tau }_{1}}-{{\tau }_{2}})\cdot \\ & (\psi ({{y}_{1}},vec({{X}_{1}}),{{\tau }_{1}})-\psi ({{y}_{2}},vec({{X}_{2}}),{{\tau }_{2}})) \\ & ={{\tau }_{1}}{{\tau }_{2}}{{(({{y}_{1}}/{{\tau }_{1}},vec({{X}_{1}})/{{\tau }_{1}})-({{y}_{2}}/{{\tau }_{2}},vec({{X}_{2}})/{{\tau }_{2}}))}^{T}}\cdot \\ & (f({{y}_{1}}/{{\tau }_{1}},vec({{X}_{1}})/{{\tau }_{1}})- \\ & f({{y}_{2}}/{{\tau }_{2}},vec({{X}_{2}})/{{\tau }_{2}})) \\ & \ge 0. \\ \end{align}$

则ψ也为单调映射.

(b) 取vec(X)^t=(1/2)^te₁(X≥0)，y^t=(1/2)^te₂，τ^t=(1/2)^t≥0，κ^t=(1/2)^t≥0，vec(S)^t=(1/2)^te₁(S≥0).

其中e₁为每个元素全为1的n²维列向量,e₂为每个元素全为1的m维列向量,当t→∞时,有

$\begin{align} & Avec{{\left( X \right)}^{t}}-{{\tau }^{t}}b={{\left( 1/2 \right)}^{t}}(A{{e}_{1}}-b)\to 0, \\ & vec{{\left( S \right)}^{t}}-{{\tau }^{t}}vec\left( C \right)+{{A}^{T}}{{y}^{t}}={{\left( 1/2 \right)}^{t}}\cdot \\ & [{{e}_{1}}-vec\left( C \right)+{{A}^{T}}{{e}_{2}}]\to 0, \\ & {{\kappa }^{t}}+vec{{\left( C \right)}^{T}}vec{{\left( X \right)}^{T}}-{{b}^{T}}{{y}^{t}}={{\left( 1/2 \right)}^{t}}\cdot \\ & [1+vec{{\left( C \right)}^{T}}{{e}_{1}}-{{b}^{T}}{{e}_{2}}]\to 0, \\ \end{align}$

则所取数列是渐进可行的,假设$(vec\left( {\tilde{X}} \right)\left( \tilde{X}\ge 0 \right),vec\left( {\tilde{S}} \right)\left( \tilde{S}\ge 0 \right),\tilde{y},\tilde{\tau }\ge 0,\tilde{\kappa }\ge 0)$是问题HMCP的任意渐进可行点,由引理1.1以及(0,vec(S),κ)^T=ψ(y,vec(X),τ)知,

$\begin{align} & vec{{\left( {\tilde{S}} \right)}^{T}}vec\tilde{X}+\tilde{\tau }\tilde{\kappa }=(\tilde{y},vec\left( {\tilde{X}} \right),\tilde{\tau }) \\ & {{(0,vec\left( {\tilde{S}} \right),\tilde{\kappa })}^{T}}~ \\ & =(\tilde{y},vec\left( {\tilde{X}} \right),\tilde{\tau })\psi (\tilde{y},vec\left( {\tilde{X}} \right),\tilde{\tau }) \\ & =0, \\ \end{align}$

即渐进可行点$(vec\left( {\tilde{X}} \right);\tilde{y};vec\left( {\tilde{S}} \right);\tilde{\tau };\tilde{\kappa })$是一个渐进互补解.

(c) “⇐”若τ^*>0，因为(vec(X^*),y^*,vec(S^*),τ^*,κ^*)是问题HMCP的极大互补解,则

$\begin{align} & vec({{S}^{*}})={{\tau }^{*}}g({{y}^{*}}/{{\tau }^{*}},vec({{X}^{*}})/{{\tau }^{*}}), \\ & vec{{({{X}^{*}})}^{T}}vec({{S}^{*}})=0, \\ & h({{y}^{*}}/{{\tau }^{*}},vec({{X}^{*}})/{{\tau }^{*}})=0, \\ \end{align}$

(3)

$所以vec({{S}^{*}})/{{\tau }^{*}}=g({{y}^{*}}/{{\tau }^{*}},vec({{X}^{*}})/{{\tau }^{*}}),$

(4)

$vec{{({{X}^{*}})}^{T}}vec({{S}^{*}})/{{({{\tau }^{*}})}^{2}}=0,$

(5)

由(3),(4),(5)得(vec(X^*)/τ^*,y^*/τ^*,vec(S^*)/τ^*)是问题MCP的一个最优互补解.

“⇒”假设问题MCP的一个解为$(vec\left( {\tilde{X}} \right),\tilde{y},vec\left( {\tilde{S}} \right))$，则可以找到$\tau =1,vec\left( X \right)=vec\left( {\tilde{X}} \right),vec\left( S \right)=vec\left( {\tilde{S}} \right),\kappa =0,y=\tilde{y}$是问题HMCP的解.

又τ^*是问题HMCP的一个极大互补解,且τ^*≥0， 则问题HMCP的每一个极大互补解中有τ^*>0 (否则,τ≡0， 矛盾).

(d) 因为(vec(X^*);y^*;vec(S^*);τ^*;κ^*)是问题HMCP的极大互补解,则

$Tr({{X}^{*}}{{S}^{*}})+{{\tau }^{*}}{{\kappa }^{*}}=0,$

(6)

$Avec({{X}^{*}})-{{\tau }^{*}}b=0,$

(7)

${{\tau }^{*}}vec\left( C \right)-{{A}^{T}}{{y}^{*}}=vec({{S}^{*}}),$

(8)

$-Tr(C{{X}^{*}})+{{b}^{T}}{{y}^{*}}={{\kappa }^{*}},$

(9)

因为κ^*>0，则τ^*=0.由式(7)—式(9)得

$Avec({{X}^{*}})=0,$

(10)

$-{{A}^{T}}{{y}^{*}}=vec({{S}^{*}}),$

(11)

$Tr(C{{X}^{*}})-{{b}^{T}}{{y}^{*}}<0,$

(12)

由式(12)得,Tr(CX^*)<0和b^Ty^*>0至少有一个成立.

当Tr(CX^*)<0成立时,由式(10)得,(P)有改良射线,则(D)强不可行,MCP强不可行.

当b^Ty^*>0成立时,由式(11)得,(D)有改良射线,则(P)强不可行,MCP强不可行.证毕.

2 齐次模型的中心路径

由定理1.1知,通过找问题HMCP的一个极大互补解来解问题MCP. 为了方便,取起始点X⁰=I,S⁰=I,τ⁰=1,κ⁰=1,y⁰=0，则X⁰S⁰=I,τ⁰κ⁰=1.

定义集合和差向量:

$\begin{align} & {{H}_{++}}={{S}^{n}}_{++}\times {{S}^{n}}_{++}\times {{R}_{++}}\times {{R}_{++}}\times {{R}^{m}}, \\ & {{R}_{D}}:=-\tau C+\underset{i=1}{\overset{m}{\mathop{\sum }}}\,{{A}_{i}}^{T}{{y}_{i}}+S, \\ & {{R}_{P}}:=\tau b-Avec\left( X \right), \\ & {{R}_{G}}:=vec{{\left( C \right)}^{T}}vec\left( X \right)-{{b}^{T}}y+\kappa . \\ \end{align}$

那么,可以定义齐次模型HMCP的不可行中心路径为:

$\begin{align} & C\left( \mu \right):=\left\{ \left( X,S,\tau ,\kappa ,y \right)\in {{H}_{++}}:\text{ }\left[ \begin{matrix} XS & 0 \\ 0 & \tau \kappa \\ \end{matrix} \right] \right.= \\ & \left. \theta I,\left[ \begin{matrix} {{R}_{P}} \\ {{R}_{D}} \\ {{R}_{G}} \\ \end{matrix} \right]=\theta \left[ \begin{matrix} R_{P}^{0} \\ R_{D}^{0} \\ R_{G}^{0} \\ \end{matrix} \right],\mu >0,0<\theta \le 1 \right\} \\ \end{align}$

(13)

表示为(X(θ),S(θ),τ(θ),κ(θ),y(θ)).

定理2.1

(a) 对任意的$0<\theta \le 1,\left[ \begin{matrix} {{R}_{P}} \\ vec({{R}_{D}}) \\ {{R}_{G}} \\ \end{matrix} \right]=\theta \left[ \begin{matrix} R_{P}^{0} \\ vec(R_{D}^{0}) \\ R_{G}^{0} \\ \end{matrix} \right]$存在一个严格正定点.

(b) 若起始点为(X⁰=I,S⁰=I,τ⁰=1,κ⁰=1,y⁰=0),则对任意的0<θ≤1， 中心路径(13)的解是存在唯一的.

(c) 中心路径是一条有界轨迹.

(d) θ→0时,存在唯一极限点(X(0),S(0),τ(0),κ(0),y(0))是问题HCMP的一个极大互补解.

证明

$\begin{align} & \left( a \right) R_{D}^{0}=-{{\tau }^{0}}C+\underset{i=1}{\overset{m}{\mathop{\sum }}}\,{{A}_{i}}^{T}{{y}^{0}}_{i}+{{S}^{0}}, \\ & R_{P}^{0}=\tau b-Avec({{X}^{0}}), \\ & R_{G}^{0}=vec{{\left( C \right)}^{T}}vec({{X}^{0}})-{{b}^{T}}{{y}^{0}}+{{\kappa }^{0}}. \\ \end{align}$

令${{J}_{++}}:=\left\{ \left[ \begin{matrix} 0 \\ vec\left( S \right) \\ \kappa \\ \end{matrix} \right]-\psi (y,vec\left( X \right),\tau ):X>0,S>0,\tau >0,\kappa >0 \right\}$，可证明,J₊₊是一个开凸集^[5-6],则(R⁰_P;vec(R⁰_D);R⁰_G)∈J₊₊， 又HMCP是渐进可行的,则0∈J₊₊，所以,对0<θ≤1，有 θ(R⁰_P;vec(R⁰_D);R⁰_G)∈J₊₊.

$\left( b \right) \left[ \begin{matrix} {{R}_{P}} \\ {{R}_{D}} \\ {{R}_{G}} \\ \end{matrix} \right]=\theta \left[ \begin{matrix} R_{P}^{0} \\ R_{D}^{0} \\ R_{G}^{0} \\ \end{matrix} \right]$展开后有如下形式:

$\begin{align} & \left[ \begin{matrix} 0 \\ vec\left( S \right) \\ \kappa \\ \end{matrix} \right]=\left[ \begin{matrix} 0 & A & -b \\ -{{A}^{T}} & 0 & vec\left( C \right) \\ {{b}^{T}} & -vec{{\left( C \right)}^{T}} & 0 \\ \end{matrix} \right]\cdot \\ & \left[ \begin{matrix} y \\ vec\left( X \right) \\ \tau \\ \end{matrix} \right]+\left[ \begin{matrix} \theta R_{P}^{0} \\ \theta vec(R_{D}^{0}) \\ \theta R_{G}^{0} \\ \end{matrix} \right]. \\ \end{align}$

则中心路径是线性互补问题的中心路径,文献^[7]已证明线性互补问题的中心路径是存在唯一的.

(c) 为了方便,我们令a=(y;vec(X);τ)，b=(0;vec(S);κ)，r=(R_p;vec(R_D);R_G). 则有r=b－ψ(a)=θr⁰，r⁰=b⁰－ψ(a⁰).

因为ψ单调,则(a⁰-a)^T(ψ(a⁰)－ψ(a))≥0，又a^Tψ(a)=0，(a⁰)^Tψ(a⁰)=0，则－(a⁰)^Tψ(a)-a^Tψ(a⁰)≥0. 所以有如下不等式:

$\begin{align} & {{a}^{T}}{{r}^{0}}={{a}^{T}}{{b}^{0}}-{{a}^{T}}\psi ({{a}^{0}}) \\ & ={{a}^{T}}{{b}^{0}}+{{b}^{T}}{{a}^{0}}-{{b}^{T}}{{a}^{0}}-{{a}^{T}}\psi ({{a}^{0}}) \\ & ={{a}^{T}}{{b}^{0}}+{{b}^{T}}{{a}^{0}}-{{({{a}^{0}})}^{T}}\cdot \\ & (\theta {{r}^{0}}+\psi \left( a \right))-{{a}^{T}}\psi ({{a}^{0}}) \\ & ={{a}^{T}}{{b}^{0}}+{{b}^{T}}{{a}^{0}}-\theta {{({{a}^{0}})}^{T}}{{r}^{0}}- \\ & {{({{a}^{0}})}^{T}}\psi \left( a \right)-{{a}^{T}}\psi ({{a}^{0}}) \\ & \ge {{a}^{T}}{{b}^{0}}+{{b}^{T}}{{a}^{0}}-\theta {{({{a}^{0}})}^{T}}{{r}^{0}} \\ & ={{a}^{T}}{{b}^{0}}+{{b}^{T}}{{a}^{0}}-\theta {{({{a}^{0}})}^{T}}({{b}^{0}}-\psi ({{a}^{0}})) \\ & ={{a}^{T}}{{b}^{0}}+{{b}^{T}}{{a}^{0}}-\theta {{({{a}^{0}})}^{T}}{{b}^{0}} \\ \end{align}$

又对任意的0<θ≤1， 有

$\begin{align} & \theta {{a}^{T}}{{r}^{0}}={{a}^{T}}\left( b-\psi \left( a \right) \right)={{a}^{T}}b \\ & =\theta \left( n+1 \right)=\theta {{({{a}^{0}})}^{T}}{{b}^{0}}. \\ \end{align}$

从以上两式可以得到,

${{a}^{T}}{{b}^{0}}+{{b}^{T}}{{a}^{0}}\le \left( 1+\theta \right){{({{a}^{0}})}^{T}}{{b}^{0}},$

即

$\begin{align} & {{(y;vec\left( X \right);\tau )}^{T}}(0;vec({{S}^{0}});{{\kappa }^{0}})+(0;vec\left( S \right);\text{ } \\ & \kappa {{)}^{T}}({{y}^{0}};vec({{X}^{0}});{{\tau }^{0}})\le \left( 1+\theta \right)({{y}^{0}};vec({{X}^{0}}); \\ & {{\tau }^{0}}{{)}^{T}}~(0;vec({{S}^{0}});{{\kappa }^{0}}) \\ \end{align}$

则(X,S,y,τ,κ)有界,即中心路径有界.

(d) 为方便,记$B=\left[ \begin{matrix} X & 0 \\ 0 & \tau \\ \end{matrix} \right],C=\left[ \begin{matrix} S & 0 \\ 0 & \kappa \\ \end{matrix} \right]$，则中心路径(13)中有BC=θI，根据文献^[8]可得中心路径存在唯一极限点(X(0),S(0),τ(0),κ(0),y(0))，并且是最优解,进而可以证明(X(0),S(0),τ(0),κ(0),y(0))是问题HCMP的一个极大互补解. □

3 具有O(nlogL)复杂性的齐次不可行算法

下面用齐次不可行内点算法计算在中心路径附近产生的迭代点.

3.1 计算搜索方向

对单调互补问题MCP采用牛顿法，可得以下线性方程组:

$\begin{align} & \left[ \begin{matrix} 0 & vec\left( C \right) & -{{A}^{T}} \\ -vec{{\left( C \right)}^{T}} & 0 & {{b}^{T}} \\ A & -b & 0 \\ \end{matrix} \right]~\left[ \begin{matrix} vec(dX) \\ d\tau \\ dy \\ \end{matrix} \right]- \\ & \left[ \begin{matrix} vec(dS) \\ d\kappa \\ 0 \\ \end{matrix} \right]=\eta \left[ \begin{matrix} vec(R_{D}^{k}) \\ R_{G}^{k} \\ R_{P}^{k} \\ \end{matrix} \right] \\ \end{align}$

(14a)

和

${{X}^{k}}dS+dX{{S}^{k}}=\gamma {{\mu }^{k}}I-{{X}^{k}}{{S}^{k}}$

(14b)

${{\kappa }^{k}}d\tau +{{\tau }^{k}}d\kappa =\gamma {{\mu }^{k}}-{{\tau }^{k}}{{\kappa }^{k}}$

(14c)

其中，${{\mu }^{k}}=(Tr({{X}^{k}}{{S}^{k}})+{{\tau }^{k}}{{\kappa }^{k}})/\left( n+1 \right),\eta \in \left( 0,1 \right),\gamma \in \left( 0,1 \right).$

从式(14)中能得到对称的dS，但是不能保证得到对称的dX，为了得到对称的dX， 使用Zhang^[9]提出的方法,将(14b)替换为

${{H}_{P}}({{X}^{k}}dS+dX{{S}^{k}})={{H}_{P}}(\gamma {{\mu }^{k}}I-{{X}^{k}}{{S}^{k}}).$

(14b1)

其中,H(·)是一个对称线性变换,定义如下:

${{H}_{P}}\left( M \right)=12[PM{{P}^{-1}}+{{P}^{-T}}{{M}^{T}}{{P}^{T}}],M\in {{R}^{n\times n}}.$

P是给定的非奇异矩阵,称为尺度矩阵,本文取P∈P(X,S)，P(X,S):={P∈Sⁿ₊₊:PXSP^－1∈Sⁿ}.

在文献^[9]中Zhang证明了若P是一个非奇异矩阵,则

${{H}_{P}}\left( M \right)=\gamma {{\mu }^{k}}I\Leftrightarrow M=\gamma {{\mu }^{k}}I,$

则(14b1)等价于

${{H}_{P}}({{X}^{k}}dS+dX{{S}^{k}})=\gamma {{\mu }^{k}}I-{{H}_{P}}({{X}^{k}}{{S}^{k}}).$

(14b2)

文献^[10]已证明,由(14a),(14b2),(14c)得到的方向dX,dS为对称方向,且是存在唯一的.

定义(X^k0,y^k,S^k0,τ^k>0,κ^k>0)是当前第k个迭代点,并且第k+1个迭代点为

(X^k+1,y^k+1,S^k+1,τ^k+1,κ^k+1)=(X^k,y^k,S^k,τ^k,κ^k)+α(dX,dy,dS,dτ,dκ).

其中,α是搜索步长,(dX,dτ,dy,dS,dκ)为搜索方向.

引用文献^[15]给出的矩阵以及向量正部和负部的概念,本文采用如下系统计算搜索方向

$\begin{align} & (dX,d\tau ,dy,dS,d\kappa ): \\ & \left[ \begin{matrix} 0 & vec\left( C \right) & -{{A}^{T}} \\ -vec{{\left( C \right)}^{T}} & 0 & {{b}^{T}} \\ A & -b & 0 \\ \end{matrix} \right]~\left[ \begin{matrix} vec(dX) \\ d\tau \\ dy \\ \end{matrix} \right]- \\ & \left[ \begin{matrix} vec(dS) \\ d\kappa \\ 0 \\ \end{matrix} \right]=\eta \left[ \begin{matrix} vec(R_{D}^{k}) \\ R_{G}^{k} \\ R_{P}^{k} \\ \end{matrix} \right], \\ \end{align}$

(15a)

$\begin{align} & {{H}_{P}}({{X}^{k}}dS+dX{{S}^{k}})={{[\gamma {{\mu }^{k}}I-{{H}_{P}}({{X}^{k}}{{S}^{k}})]}^{-}}+ \\ & \sqrt{n+1}{{[\gamma {{\mu }^{k}}I-{{H}_{P}}({{X}^{k}}{{S}^{k}})]}^{+}}, \\ \end{align}$

(15b)

$\begin{align} & {{\kappa }^{k}}d\tau +{{\tau }^{k}}d\kappa ={{[\gamma {{\mu }^{k}}-{{\tau }^{k}}{{\kappa }^{k}}]}^{-}}+\sqrt{n+1}\cdot \\ & {{[\gamma {{\mu }^{k}}-{{\tau }^{k}}{{\kappa }^{k}}]}^{+}}. \\ \end{align}$

(15c)

其中,μ^k=(Tr(X^kS^k)+τ^kκ^k)/(n+1)，η∈(0,1)，γ∈(0,1).

下面用克罗内克积的形式等价表示(15b):

$\begin{align} & \left( 15b \right)\Leftrightarrow [P({{X}^{k}}dS+dX{{S}^{k}}){{P}^{-1}}+ \\ & {{P}^{-T}}({{S}^{k}}dX+dS{{X}^{k}}){{P}^{T}}]/2 \\ & ={{[\gamma {{\mu }^{k}}I-{{H}_{P}}({{X}^{k}}{{S}^{k}})]}^{-}}+ \\ & \sqrt{n+1}{{[\gamma {{\mu }^{k}}I-{{H}_{P}}({{X}^{k}}{{S}^{k}})]}^{+}} \\ & \Leftrightarrow \text{ }{{E}^{k}}vec(dX)+{{F}^{k}}vec(dS)=vec(R_{C}^{k}), \\ \end{align}$

(15b1)

其中,

$\begin{align} & E:=({{P}^{-T}}S\otimes P+P\otimes {{P}^{-T}}S)/2, \\ & F:=(PX\otimes {{P}^{-T}}+{{P}^{-T}}\otimes PX)/2, \\ & R_{C}^{k}:={{[\gamma {{\mu }^{k}}I-{{H}_{P}}({{X}^{k}}{{S}^{k}})]}^{-}}+ \\ & \sqrt{n+1}{{[\gamma {{\mu }^{k}}I-{{H}_{P}}({{X}^{k}}{{S}^{k}})]}^{+}}. \\ \end{align}$

因为本文取P∈P(X,S)，文献^[17]已证明E^k，F^k为半正定矩阵(不一定对称).与系统(14)相同,可以证明系统(15)的解也是存在唯一的.

下面求解由系统(15a),(15b1),(15c)得到的迭代搜索方向(dX,dy,dS,dτ,dκ).

由(15b1)得,

$vec(dS)={{({{F}^{k}})}^{-1}}[vec({{R}^{k}}_{C})-{{E}^{k}}vec(dX)],$

(16)

由(15c)得,

$d\kappa ={{({{\tau }^{k}})}^{-1}}({{r}_{C}}^{k}-{{\kappa }^{k}}d\tau ),$

(17)

其中,${{r}_{C}}^{k}={{[\gamma {{\mu }^{k}}-{{\tau }^{k}}{{\kappa }^{k}}]}^{-}}+\sqrt{n+1}{{[\gamma {{\mu }^{k}}-{{\tau }^{k}}{{\kappa }^{k}}]}^{+}}.$

将式(16)、式(17)代入式(15a)减少牛顿方程,得牛顿搜索方向(dX,dy,dS,dτ,dκ)为如下系统的解:

$\begin{align} & \left[ \begin{matrix} {{({{F}^{k}})}^{-1}}{{E}^{k}} & vec\left( C \right) & -{{A}^{T}} \\ -vec{{\left( C \right)}^{T}} & {{({{\tau }^{k}})}^{-1}}{{\kappa }^{k}} & {{b}^{T}} \\ A & -b & 0 \\ \end{matrix} \right]\left[ \begin{matrix} vec(dX) \\ d\tau \\ dy \\ \end{matrix} \right]= \\ & \left[ \begin{matrix} \eta vec({{R}^{k}}_{D})+{{({{F}^{k}})}^{-1}}vec(R_{C}^{k}) \\ \eta R_{G}^{k}+{{({{\tau }^{k}})}^{-1}}{{r}_{C}}^{k} \\ \eta R_{P}^{k} \\ \end{matrix} \right]. \\ \end{align}$

(18)

令$K:=\left[ \begin{matrix} {{({{F}^{k}})}^{-1}}{{E}^{k}} & -{{A}^{T}} \\ A & 0 \\ \end{matrix} \right]$，其中A是m×n²阶矩阵,E^k和F^k均是n²阶方阵,则K是m+n²阶方阵.

引理3.1 若K有如上定义,则

${{K}^{-1}}=\left[ \begin{matrix} {{Q}^{-1}}-{{Q}^{-1}}{{A}^{T}}{{(A{{Q}^{-1}}{{A}^{T}})}^{-1}}A{{Q}^{-1}} & {{Q}^{-1}}{{A}^{T}}{{(A{{Q}^{-1}}{{A}^{T}})}^{-1}} \\ {{(A{{Q}^{-1}}{{A}^{T}})}^{-1}}A{{Q}^{-1}} & -{{(A{{Q}^{-1}}{{A}^{T}})}^{-1}} \\ \end{matrix} \right],$

其中,Q=(F^k)/^－1E^k.

证明 易验证K^－1K=KK^－1=I_m+n².□

由上述可知,牛顿搜索方向(18)等价于

$\begin{align} & \left[ \begin{matrix} K & (vec\left( C \right);-b) \\ {{(-vec\left( C \right);b)}^{T}} & {{({{\tau }^{k}})}^{-1}}{{\kappa }^{k}} \\ \end{matrix} \right]\cdot \\ & \left[ \begin{matrix} (vec(dX);dy) \\ d\tau \\ \end{matrix} \right]=\left[ \begin{matrix} {\tilde{r}} \\ \eta (R_{G}^{k})+{{({{\tau }^{k}})}^{-1}}{{r}_{C}}^{k} \\ \end{matrix} \right], \\ \end{align}$

(19)

其中,$\tilde{r}:=\left[ \begin{matrix} \eta vec(R_{D}^{k})+{{({{F}^{k}})}^{-1}}vec({{R}^{k}}_{C}) \\ \eta R_{P}^{k} \\ \end{matrix} \right].$

定义p和q是如下线性系统的解:

$Kp=(vec\left( C \right);-b);\text{ }Kq=\tilde{r}$

由式(19)知,

(第二行)

$\begin{align} & {{({{\tau }^{k}})}^{-1}}{{\kappa }^{k}}d\tau +{{(-vec\left( C \right);b)}^{T}}(vec(dX);dy) \\ & =\eta (R_{G}^{k})+{{({{\tau }^{k}})}^{-1}}{{r}_{C}}^{k}, \\ \end{align}$

(20)

(第一行)

$K(vec(dX);dy)+(vec\left( C \right);-b)d\tau =\tilde{r}$

即 K(vec(dX);dy)+Kpdτ=Kq.

又由引理3.1知,K是非奇异矩阵,则

$(vec(dX);dy)=q-pd\tau ,$

(21)

将其代入(20)得

(τ^k)^－1κ^kdτ+(－vec(C);b)^T(q-pdτ)

=η(R^k_G)+(τ^k)^－1r_C^k,

即

$\begin{align} & [{{({{\tau }^{k}})}^{-1}}{{\kappa }^{k}}-(-vec\left( C \right),b)p]d\tau \\ & =\eta R_{G}^{k}+{{({{\tau }^{k}})}^{-1}}{{r}_{C}}^{k}-(-vec\left( C \right),b)q, \\ \end{align}$

(22)

此式即为搜索方向系统(15)的等价形式.

因为搜索方向系统的解是存在唯一的,所以dτ可以从式(22)中获得,进而可以从式(21)中得到vec(dX),dy，从式(16),式(17)中得到vec(dS),dκ.

通过以上分析,计算搜索方向的主要工作在于计算向量p和q.

定义$\bar{K}:=\left[ \begin{matrix} {{({{F}^{k}})}^{-1}}{{E}^{k}} & {{A}^{T}} \\ A & 0 \\ \end{matrix} \right]$为对称矩阵,则p和q是如下线性系统的解

$K\left( p;q \right)=\bar{K}\left( p;-q \right)=(vec\left( C \right);-b;\tilde{r}).$

K是对称矩阵,据文献有Bunch-Parlett分解K=LDL^T，其中L是下三角矩阵,D是一个非奇异的分块对角矩阵,那么,可以求出p和q，进而求得方向(dX,dy,dS,dτ,dκ).

3.2 不可行性和对偶间隙的下降关系

本文将限制关系γ=τ₁，τ₁≤1/4，β≤1/2成立.由文献^[15]引理3.5知,若(X,y,S,τ,κ)∈N(τ₁,β)，τ₁≤1/4，β≤1/2， 有如下不等式成立:

$\begin{align} & \left( \tau -1 \right){{\mu }^{k}}\left( n+1 \right)\le Tr(R_{C}^{k})+r_{C}^{k}\le \\ & \left( \tau +\beta \tau /\sqrt{n+1}-1 \right){{\mu }^{k}}\left( n+1 \right)<0. \\ \end{align}$

所以本文定义

$\eta =\frac{-Tr({{R}^{k}}_{C})+r_{C}^{k}}{\left( n+1 \right){{\mu }^{k}}}.$

(23)

其中关于N(τ₁,β)的介绍在下节给出说明.

引理3.2 算法更新后,迭代点的不可行性有如下结论：

(a) R^k+1_P=(1－αη)R^k_P.

(b) R^k+1_D=(1－αη)R^k_D.

(c) R^k+1_G=(1－αη)R^k_G.

证明 (a) 根据R_P,R_D,R_G的定义及(15a)可得,

R^k+1_P=(τ^k+αdτ)b-A(vec(X^k)+αvec(dX))

=τ^kb-Avec(X^k)+α[bdτ-Avec(dX)]

=R^k_P-αηR^k_P=(1-αη)R^k_P.

同理,可得(b),(c).

引理3.3 若τ₁≤1/4，β≤1/2， 则由系统(15)定义的方向满足

vec(dX)^Tvec(dS)+dτdκ=0.

证明 由系统(15)得

$\begin{align} & vec{{(dX)}^{T}}vec(dS)+d\tau d\kappa \\ & =vec{{(dX)}^{T}}[-{{A}^{T}}dy+vec\left( C \right)d\tau -\eta vec(R_{D}^{k})]+ \\ & d\tau [-vec{{\left( C \right)}^{T}}vec(dX)+{{b}^{T}}dy-\eta R_{G}^{k}] \\ & =[-vec{{(dX)}^{T}}{{A}^{T}}+d\tau {{b}^{T}}]dy- \\ & vec{{(dX)}^{T}}vec(R_{D}^{k})\eta -d\tau R_{G}^{k}\eta ~ \\ & =-[{{(R_{P}^{k})}^{T}}dy+vec{{(dX)}^{T}}vec(R_{D}^{k})+d\tau R_{G}^{k}]\eta . \\ \end{align}$

(24)

由R_P,R_D,R_G的定义知

$\begin{align} & {{({{R}_{P}})}^{T}}y+vec{{\left( X \right)}^{T}}vec({{R}_{D}})+\tau {{R}_{G}} \\ & =vec{{\left( X \right)}^{T}}vec\left( S \right)+\tau \kappa \\ \end{align}$

对(y^k,X^k,S^k,τ^k,κ^k)和(y^k+dy,X^k+dX,S^k+dS,τ^k+dτ,κ^k+κ)，有

$\begin{align} & {{({{y}^{k}})}^{T}}{{R}^{k}}_{P}+vec{{({{X}^{k}})}^{T}}vec(R_{D}^{k})+{{\tau }^{k}}R_{G}^{k} \\ & =vec{{({{X}^{k}})}^{T}}vec({{S}^{k}})+{{\tau }^{k}}{{\kappa }^{k}} \\ \end{align}$

(25)

和

$\begin{align} & {{({{y}^{k}}+dy)}^{T}}R_{P}^{k}+{{(vec({{X}^{k}})+vec(dX))}^{T}}\cdot \\ & vec(R_{D}^{k})+({{\tau }^{k}}+d\tau )R_{G}^{k} \\ & ={{(vec({{X}^{k}})+vec(dX))}^{T}}(vec({{S}^{k}})+ \\ & vec(dS))+({{\tau }^{k}}+d\tau )({{\kappa }^{k}}+d\kappa ). \\ \end{align}$

(26)

由式(24)—式(26)得如下等式

$\begin{align} & [vec{{({{X}^{k}})}^{T}}vec({{S}^{k}})+{{\tau }^{k}}{{\kappa }^{k}}]\left( 1-\eta \right)+ \\ & [{{(R_{P}^{k})}^{T}}dy+vec{{(dX)}^{T}}vec(R_{D}^{k})+d\tau R_{G}^{k}]\cdot \\ & \left( 1-\eta \right) \\ & =[vec{{({{X}^{k}})}^{T}}vec({{S}^{k}})+{{\tau }^{k}}{{\kappa }^{k}}]+ \\ & [vec{{(dX)}^{T}}vec(dS)+d\tau d\kappa ]+ \\ & vec{{({{X}^{k}})}^{T}}vec(dS)+vec{{(dX)}^{T}}vec({{S}^{k}})+ \\ & {{\tau }^{k}}d\kappa +d\tau {{\kappa }^{k}}. \\ \end{align}$

将式(24)代入上式,得

$\begin{align} & {{(R_{P}^{k})}^{T}}dy+vec{{(dX)}^{T}}vec(R_{D}^{k})+d\tau R_{G}^{k} \\ & =[vec{{({{X}^{k}})}^{T}}vec({{S}^{k}})+{{\tau }^{k}}{{\kappa }^{k}}]\eta + \\ & [vec{{({{X}^{k}})}^{T}}vec(dS)+vec{{(dX)}^{T}}vec({{S}^{k}})+ \\ & {{\tau }^{k}}d\kappa +d\tau {{\kappa }^{k}}]. \\ \end{align}$

(27)

由式(15b),式(15c)得

$\begin{align} & vec{{({{X}^{k}})}^{T}}vec(dS)+vec{{(dX)}^{T}}vec({{S}^{k}})+{{\tau }^{k}}d\kappa + \\ & d\tau {{\kappa }^{k}}=Tr(R_{C}^{k})+r_{C}^{k}. \\ \end{align}$

(28)

由式(24),式(27),式(28)得

$\matrix{ \matrix{ vec{\left( {dX} \right)^T}vec\left( {dS} \right) + d\tau d\kappa \hfill \cr = - [\left( {vec{{\left( {{X^k}} \right)}^T}vec\left( {{S^k}} \right) + {\tau ^k}{\kappa ^k}} \right)\eta + \hfill \cr} \hfill \cr \matrix{ Tr(R_C^k) + r_C^k\eta ] \hfill \cr = - [\left( {n + 1} \right){\mu ^k}\eta + Tr(R_C^k) + r_C^k] \hfill \cr \eta = 0 \hfill \cr} \hfill \cr } $

□

引理3.4 若τ₁≤1/4，β≤1/2， 则算法更新后对偶间隙满足如下等式

μ^k+1=(1－αη)μ^k.

证明

$\begin{align} & {{\mu }^{k+1}}=[vec{{({{X}^{k+1}})}^{T}}vec({{S}^{k+1}})+{{\tau }^{k+1}}{{\kappa }^{k+1}}]/\left( n+1 \right) \\ & =[{{(vec({{X}^{k}})+\alpha vec(dX))}^{T}}(vec({{S}^{k}})+ \\ & \alpha vec(dS))+({{\tau }^{k}}+\alpha d\tau )({{\kappa }^{k}}+\alpha d\kappa )]/\left( n+1 \right) \\ & =\{[vec{{({{X}^{k}})}^{T}}vec({{S}^{k}})+{{\tau }^{k}}{{\kappa }^{k}}]+ \\ & \alpha [vec{{({{X}^{k}})}^{T}}vec(dS)+ \\ & vec{{(dX)}^{T}}vec({{S}^{k}})+{{\tau }^{k}}d\kappa +d\tau {{\kappa }^{k}}]+ \\ & {{\alpha }^{2}}[vec{{(dX)}^{T}}vec(dS)+d\tau d\kappa ]\}/\left( n+1 \right) \\ & ={{\mu }^{k}}+\alpha Tr({{R}^{k}}_{C}+{{r}^{k}}_{C})/\left( n+1 \right) \\ & =\left( 1-\alpha \eta \right){{\mu }^{k}} \\ \end{align}$

其中,第4个等式用到引理3.3.□

引理3.2和引理3.4说明当$\eta =\frac{-Tr(R_{C}^{k}+r_{C}^{k})}{\left( n+1 \right){{\mu }^{k}}},{{\tau }_{1}}\le 1/4,\beta \le 1/2$时,不可行性和最优性以相同的速度降低,下降系数均为(1－αη)， 这是定义中心路径(13)的原因,以及下节将证明的迭代点始终保持在一个宽邻域内,都是本文提出新算法的关键点.

3.3 齐次不可行内点算法

上一节已限制尺度矩阵:

P(X,S):={P∈Sⁿ₊₊:PXSP^－1∈Sⁿ}.

本文取P=W^1/2， 得到的方向称为NT方向. 其中矩阵$W={{S}^{\frac{1}{2}}}{{({{S}^{\frac{1}{2}}}X{{S}^{\frac{1}{2}}})}^{-\frac{1}{2}}}{{S}^{\frac{1}{2}}}$或者$W={{X}^{-\frac{1}{2}}}{{({{X}^{\frac{1}{2}}}S{{X}^{\frac{1}{2}}})}^{\frac{1}{2}}}{{X}^{-\frac{1}{2}}}$.可以得到E^k，F^k的值^[14],容易验证P∈P(X,S)和PXP=P^－1SP^－1.

定义F₊₊=Sⁿ₊₊×R^m×Sⁿ₊₊， SDP中常用的一个宽邻域是－∞邻域,定义为:

N^－_∞(1－τ):={(X,y,S)∈F₊₊:λ_min(XS)≥τμ}.

其中,0<τ<1，μ^k=Tr(X^kS^k)/(n+1).

后来,Ai和Zhang^[11]提出求解线性互补问题的新的宽邻域,Li和Terlaky^[12]将其推广到半定规划,Liu等^[1]又将其推广到对称锥规划,刘新泽^[13]将这种新的宽邻域推广到关于1范数的宽邻域. 本文采用Li和Terlaky^[12]提出的新的宽邻域:

N(τ₁,τ₂)={(X,y,S)∈F₊₊:

‖[τ₁μI－X¹²SX¹²]⁺‖_F≤(τ₁－τ₂)μ}.

其中,0<τ₂<τ₁<1.

文献^[16]已证明,下面关系成立:

$\begin{align} & N_{\infty }^{-}(1-{{\tau }_{1}})N({{\tau }_{1}},{{\tau }_{2}})N_{\infty }^{-}(1-{{\tau }_{2}}), \\ & \forall 0<{{\tau }_{2}}<{{\tau }_{1}}<1. \\ \end{align}$

为方便,记β:=(τ₁－τ₂)/τ₁，则β∈(0,1)且

$\begin{align} & N\left( {{\tau }_{1}},\beta \right)=\{\left( X,y,S \right)\in {{F}_{++}}: \\ & \|{{\left[ {{\tau }_{1}}\mu I-{{X}^{\frac{1}{2}}}S{{X}^{\frac{1}{2}}} \right]}^{+}}{{\|}_{F}}\le \beta {{\tau }_{1}}\mu \} \\ \end{align}$

基于以上分析,给出算法的框架:

算法3.1

输入精度ε>0，以及参数τ₁≤1/4，β≤1/2.

初始点满足(X⁰,y⁰,S⁰,τ⁰,κ⁰)∈N(τ₁,β)，并假设k:=0，μ⁰=(Tr(X⁰S⁰)+τ⁰κ⁰)/(n+1).

步骤1 如果μ^k≤εμ⁰， 则停止.

步骤2 计算NT尺度矩阵P=W^1/2，由(23)计算参数η.

步骤3 由系统(15)计算方向(dX^k,dy^k,dS^k,dτ^k,dκ^k).

步骤4 计算使(X(α^k),y(α^k),S(α^k),τ(α^k),κ(α^k))∈N(τ₁,β)

成立的步长${{\alpha }^{k}}=\frac{2\beta {{\tau }_{1}}^{2}}{(1+\beta {{\tau }_{1}})\sqrt{n+1}}.$.

步骤5 令

$({{X}^{k+1}},{{y}^{k+1}},{{S}^{k+1}},{{\tau }^{k+1}},{{\kappa }^{k+1}})=(X({{\alpha }^{k}}),y({{\alpha }^{k}}),S({{\alpha }^{k}}),\tau ({{\alpha }^{k}}),\kappa ({{\alpha }^{k}})),$${{\mu }^{k+1}}=(Tr({{X}^{k+1}}{{S}^{k+1}})+{{\tau }^{k+1}}{{\kappa }^{k+1}})/\left( n+1 \right),k:=k+1,$ 并转步骤1.

3.4 算法的复杂性分析

本节在给出一些引理之后,将证明算法在宽邻域内是收敛的,且迭代复杂度为O(nlogL).

将问题元素尺度化如下:

$\begin{align} & \left( \hat{X},\hat{y},\hat{S} \right)=\left( PXP,y,{{P}^{\text{-}1}}S{{P}^{\text{-}1}} \right), \\ & (\hat{C},{{{\hat{A}}}_{i}},{{{\hat{b}}}_{i}})=({{P}^{\text{-}1}}C{{P}^{\text{-}1}},{{P}^{\text{-}1}}Ai{{P}^{\text{-}1}},{{b}_{i}})\text{,} \\ & (d\hat{X},d\hat{y},d{{{\hat{S}}}^{)=(}}PdXP,dy,{{P}^{\text{-}1}}dS{{P}^{\text{-}1}}). \\ \end{align}$

可以得到

$\begin{align} & P(X,S)=\left\{ P\in S_{++}^{n}:{{(PXS{{P}^{\text{-}1}})}^{T}}=PXS{{P}^{\text{-}1}} \right\} \\ & =\left\{ P\in S_{++}^{n}:\hat{X}\hat{S}=\hat{S}\hat{X} \right\} \\ \end{align}$

为方便,记$X=\left[ \begin{matrix} X & 0 \\ 0 & \tau \\ \end{matrix} \right]\in S_{++}^{n+1},S=\left[ \begin{matrix} S & 0 \\ 0 & \kappa \\ \end{matrix} \right]\in S_{++}^{n+1},$那么

$(15b)\text{和}(15c)\Leftrightarrow {{{\hat{X}}}^{k}}d\hat{S}+d\hat{X}{{{\hat{S}}}^{k}}=R_{C}^{k}\text{,}$

(29)

$\Leftrightarrow {{H}_{P}}({{X}^{k}}dS+dX{{S}^{k}})=R_{C}^{k},$

(30)

其中,

$\begin{align} & {{{\hat{R}}}_{C}}={{\left[ \gamma \mu I-\hat{X}\hat{S} \right]}^{-}}+\sqrt{n+1}{{\left[ \gamma \mu I-\hat{X}\hat{S} \right]}^{+}}, \\ & \mu =Tr\left( XS \right)/\left( n+1 \right). \\ \end{align}$

由式(29)知,方向系统(15b),(15c)可以转化为

${{{\hat{E}}}^{k}}vec(d\hat{X})+{{{\hat{F}}}^{k}}vec(d\hat{S})=vec(\hat{R}_{C}^{k}).$

(31)

其中,

${{{\hat{E}}}^{=}}(\hat{S}\otimes I+I\otimes \hat{S})/2\text{,}{{{\hat{F}}}^{=}}(\hat{X}\otimes I+I\otimes {{{\hat{X}}}^{)}}/2.$

(32)

注意到,对于NT尺度可得$\hat{X}=\hat{S}\text{,}\hat{E}=\hat{F}$，并且矩阵$\hat{X}\hat{S},\hat{S}\hat{X},XS,SX,{{X}^{1/2}}S{{X}^{1/2}},{{S}^{1/2}}X{{S}^{1/2}}$都相似.

引理3.5^[16] 设U,V∈Sⁿ， 则下面不等式成立:

‖(U+V)⁺‖_F≤‖U⁺+V⁺‖_F≤

‖U⁺‖_F+‖V⁺‖_F.

引理3.6^[16] 设(X,y,S)∈N(τ₁,β)，矩阵${{{\hat{E}}}^{k}},{{{\hat{F}}}^{k}}$由式(32)定义,并且τ₁≤1/4，β≤1/2， 则以下不等式成立

$\begin{align} & \left\| {{{\hat{E}}}^{-1}} \right.vec[\gamma \mu I\text{-}\hat{X}\hat{S}){{\left. {{]}^{-}} \right\|}^{2}}\le (n+1)\mu \text{,} \\ & \left\| {{{\hat{E}}}^{-1}} \right.vec[\gamma \mu I\text{-}\hat{X}\hat{S}){{\left. {{]}^{+}} \right\|}^{2}}\le \beta {{\tau }_{1}}\mu . \\ \end{align}$

引理3.7^[16] 设u,v,r∈Rⁿ与A,B∈R^n×n满足Au+Bv=r.若BA^T∈Sⁿ₊₊， 则

$\begin{align} & \|{{(B{{A}^{T}})}^{-1/2}}Au{{\|}^{2}}+\|{{(B{{A}^{T}})}^{-1/2}}Bv{{\|}^{2}}+ \\ & 2{{u}^{T}}v=\|{{(B{{A}^{T}})}^{-1/2}}r{{\|}^{2}}. \\ \end{align}$

引理3.8 假设(X,y,S)∈N(τ₁,β)，矩阵${{{\hat{E}}}^{k}},{{{\hat{F}}}^{k}}$由式(32)定义,$\eta =-\frac{Tr(R_{C}^{k}+r_{C}^{k})}{\left( n+1 \right){{\mu }^{k}}},{{\tau }_{1}}\le 1/4,\beta \le \frac{1}{2}$，若$\alpha =\frac{2\beta {{\tau }_{1}}^{2}}{(1+\beta {{\tau }_{1}})\sqrt{n+1}}$则以下不等式成立

$\frac{\alpha }{\sqrt{n+1}}{{\left\| {{\left[ d\hat{X}d\hat{S} \right]}^{-}} \right\|}_{F}}\le \beta {{\tau }_{1}}\mu \left( \alpha \right).$

证明

$\begin{align} & {{\left\| {{\left[ d\hat{X}d\hat{S} \right]}^{-}} \right\|}_{F}}\le {{\left\| d\hat{X}d\hat{S} \right\|}_{F}}\le {{\left\| d\hat{X} \right\|}_{F}}\le {{\left\| d\hat{S} \right\|}_{F}} \\ & =\left\| vec(d\hat{X}) \right\|\left\| vec(d\hat{S}) \right\| \\ & \le \frac{1}{2}\left( \left\| vec(d\hat{X}) \right\|+\left\| vec(d\hat{S}) \right\| \right), \\ \end{align}$

利用引理3.6,引理3.7,引理3.3正交性以及式(31)可得

$\begin{align} & {{\left\| vec(d\hat{X}) \right\|}^{2}}+{{\left\| vec(d\hat{S}) \right\|}^{2}} \\ & \le {{\left\| {{{\hat{E}}}^{-1}}vec{{\left[ \gamma \mu I-\hat{X}\hat{S} \right]}^{-}} \right\|}^{2}}+ \\ & n+1{{\left\| {{{\hat{E}}}^{-1}}vec{{\left[ \gamma \mu I-\hat{X}\hat{S} \right]}^{+}} \right\|}^{2}} \\ & \le \left( 1+\beta {{\tau }_{1}} \right)\left( n+1 \right)\mu , \\ \end{align}$

$\begin{align} & \text{又} Tr(\hat{X}(\alpha )\hat{S}(\alpha )) \\ & =Tr(\hat{X}\hat{S})+\alpha [Tr({{\tau }_{1}}\mu I\text{-}\hat{X}\hat{S})+ \\ & (\sqrt{n+1}\text{-}1)Tr{{({{\tau }_{1}}\mu I\text{-}\hat{X}\hat{S})}^{+}}]. \\ & \mu (\alpha )=\mu +\alpha [({{\tau }_{1}}-1)\mu + \\ & \frac{\sqrt{n+1}-1}{n+1}Tr{{({{\tau }_{1}}\mu I\text{-}\hat{X}\hat{S})}^{+}}] \\ & \ge \mu +\alpha ({{\tau }_{1}}-1)\mu \\ & \ge \mu +{{\tau }_{1}}\mu \text{-}\mu ={{\tau }_{1}}\mu . \\ \end{align}$

所以,$\begin{align} & \frac{\alpha }{\sqrt{n+1}}{{\left\| {{\left[ d\hat{X}d\hat{S} \right]}^{-}} \right\|}_{F}}-\beta {{\tau }_{1}}\mu \left( \alpha \right)\le \\ & \frac{\alpha }{2\sqrt{n+1}}\left( 1+\beta {{\tau }_{1}} \right)\left( n+1 \right)\mu -\beta {{\tau }_{1}}^{2}\mu \\ & =\left[ \frac{2\beta {{\tau }_{1}}^{2}}{\left( 1+\beta {{\tau }_{1}} \right)\sqrt{n+1}}\cdot \right.\frac{1}{2\sqrt{n+1}}\cdot \\ & \left. \left| \left( 1+\beta {{\tau }_{1}} \right)\left( n+1 \right)\mu -\beta {{\tau }_{1}}^{2} \right. \right]\mu \\ & =0. \\ \end{align}$ □

引理3.9 假设(X,y,S)∈N(τ₁,β)，矩阵${{{\hat{E}}}^{k}},{{{\hat{F}}}^{k}}$由式(32)定义,$\eta =-\frac{Tr(R_{C}^{k})+r_{C}^{k}}{\left( n+1 \right){{\mu }^{k}}},{{\tau }_{1}}\le 1/4,\beta \le \frac{1}{2}$，若$\alpha =\frac{2\beta {{\tau }_{1}}^{2}}{\left( 1+\beta {{\tau }_{1}} \right)\sqrt{n+1}}$ 则

(X(α),y(α),S(α))∈N(τ₁,β).

证明 因为$\alpha =\frac{2\beta {{\tau }_{1}}^{2}}{\left( 1+\beta {{\tau }_{1}} \right)\sqrt{n+1}}<\frac{1}{\sqrt{n+1}}$，由文献^[15]引理4.9得,

$\begin{align} & {{\left\| {{\left[ {{\tau }_{1}}\mu (\alpha )I\text{-}(\hat{X}\hat{S}+\alpha {{{\hat{R}}}_{C}}) \right]}^{+}} \right\|}_{F}}\le \\ & (1\text{-}\alpha \sqrt{n+1})\beta {{\tau }_{1}}\mu (\alpha ). \\ \end{align}$

所以,

$\begin{align} & {{\left\| {{\left[ {{\tau }_{1}}\mu (\alpha )I\text{-X}{{(\alpha )}^{1/2}}\text{S}(\alpha )\text{X}{{(\alpha )}^{1/2}} \right]}^{+}} \right\|}_{F}} \\ & ={{\left\| {{\left[ {{\tau }_{1}}\mu (\alpha )I\text{-}\hat{X}(\alpha )\hat{S}(\alpha ) \right]}^{+}} \right\|}_{F}} \\ & ={{\left\| {{\left[ {{\tau }_{1}}\mu (\alpha )I\text{-}(\hat{X}+\alpha d\hat{X})\hat{S}(\hat{S}+\alpha d\hat{S}) \right]}^{+}} \right\|}_{F}} \\ & ={{\left\| {{\left[ {{\tau }_{1}}\mu (\alpha )I\text{-}(\hat{X}\hat{S}+\alpha {{{\hat{R}}}_{C}}) \right]}^{+}} \right\|}_{F}}+ \\ & {{\alpha }^{2}}{{\left\| {{\left[ d\hat{X}d\hat{S} \right]}^{-}} \right\|}_{F}} \\ & \le \left( 1-\alpha \sqrt{n+1} \right)\beta {{\tau }_{1}}\mu \left( \alpha \right)+\alpha \sqrt{n+1}\beta {{\tau }_{1}}\mu \left( \alpha \right) \\ & =\beta {{\tau }_{1}}\mu \left( \alpha \right). \\ \end{align}$

其中不等式利用了引理3.8.

引理3.10^[15] 若τ₁≤1/4，β≤1/2， 则

μ(α)≤(1－ξα)μ，

其中ξ=1－τ₁－βτ₁.

定理3.1 当使用NT搜索方向时,算法4.1的迭代复杂度是$O\left( \sqrt{n}\log L \right)$.

证明 由引理3.10得μ^k≤(1－ξα)μ^k－1≤(1－ξα)^k－1μ⁰，即Tr(X^kS^k)≤(1－ξα)^kTr(X⁰S⁰).

所以,要使Tr(X^kS^k)≤ε成立,即算法终止,只需保证

(1－ξα)^kTr(X⁰S⁰)≤ε

即可,即

$k\log \left( 1-\xi \alpha \right)\le -\log \frac{Tr{{X}^{0}}{{S}^{0}}}{\varepsilon }.$

又 log(1－ξα)≤－ξα，

所以,

$\begin{align} & k\ge \frac{1}{\xi \alpha }\frac{logTr({{X}^{0}}{{S}^{0}})}{\varepsilon }=\frac{(1+\beta {{\tau }_{1}})\sqrt{n+1}}{2\beta {{\tau }_{1}}^{2}\xi }\cdot \\ & \frac{logTr({{X}^{0}}{{S}^{0}})}{\varepsilon } \\ \end{align}$

即算法4.1至多经过$\left\lceil \frac{(1+\beta {{\tau }_{1}})\sqrt{n+1}}{2\beta {{\tau }_{1}}^{2}\xi }\frac{logTr({{X}^{0}}{{S}^{0}})}{\varepsilon } \right\rceil $步迭代后停止. □

4 数值实验

本节将比较算法3.1和Rangarajan^[2]中算法3.2. 所有测试是在Windows XP系统下,使用MATLAB R2010(b),测试一些随机产生的半定规划问题,包括随机半定规划(表 1),最大割问题(表 2),教育测试问题(表 3),模极小化问题(表 4),关于4个问题的具体理论详见文献^[17].

算法3.1中定义:

normp=‖τb-AX‖_F，

normd=‖τC-A^Ty-S‖_F.

算法3.2中定义:

normp=‖b-AX‖_F，

normd=‖C-A^Ty-S‖_F.

对算法选取最优参数为τ₁=0.05，β=0.01. 使用NT尺度化方向,对同一个m和n， 运行10次取平均值结果列在下表.从表 1—表 4可以看出算法3.1的平均迭代次数比算法3.2减少了38.35%.虽然这2个程序是粗糙的,但整体上不影响2个算法的比较.

5 结论

通过以上讨论,算法迭代后得到τ和κ的值,结合定理1.1可得知可行点的存在性,即原问题是可行问题或者不可行问题.若原问题是可行问题,结合定理1.1可以得到原问题的最优解.除此结论还可得到如下结论: 采用NT尺度化矩阵和新的宽邻域得到的复杂度比一般算法的复杂度小,同时也具有更好的实验结果,从而证明本文研究的算法是有效的.

本文提出的是齐次不可行内点算法,此外,也可以研究齐次可行内点算法,具体过程与本文十分相似,进而可以求解辅助问题HMCP,计算迭代复杂度,测试实验结果,与本文所提出的算法进行比较,具体可以参考文献^[16].