半定规划(Semidefinite Programming, SDP)是线性规划(Linear Programming, LP)的一种拓广.与线性规划的约束条件不同, 半定规划是在满足约束条件“对称矩阵的仿射组合半正定”的条件下, 求目标函数极大(极小)化的问题.在半定规划模型下, 线性规划、凸二次规划(Convex Quadratic Programming, CQP)、二阶锥规划(Second-order Cone Programming, SOCP)等规划问题均可看作其特殊形式[1].于此同时, 半定规划在特征值优化、控制论、组合优化以及工业工程设计等许多领域都具有广泛的应用.目前, 半定规划的研究主要集中在求解方法上, 线性规划内点法的巨大成功, 启发人们将线性规划的内点法推广到半定规划上.在理论上, 内点算法已经非常成熟, 被证明是求解半定规划问题的有效方法之一[2], 但是内点法的运行通常需要一对严格的原始一对偶初始可行解, 而选择合适的严格原始一对偶初始可行解非常困难, 更为重要的是内点算法的每步迭代涉及到大量矩阵运算, 这在很大程度上限制了其在大规模问题中的应用.
变分不等式问题(Variational Inequalities, VIP)不仅是数学规划领域的一个重要的组成部分, 而且它还与其它数学规划问题有着紧密的联系.线性规划、非线性规划及互补问题都是变分不等式的特殊情形[3].对于变分不等式问题的求解算法研究在近十几年中取得了很多优秀的研究成果. Pang和Chan在文献[4]中全面地介绍了变分不等式问题的解法及其收敛性, 主要包括线性化求解方法, 非线性化求解方法(如非线性Jacobi方法)等. 20世纪90年代, 何炳生教授提出了投影收缩算法[5], 韩乔明基于投影收缩算法的思想提出了解半定规划的二次摄动方法[6], 关秀翠证明了二次半定规划的最优性条件与线性对称变分不等式的等价性, 并利用投影收缩算法求解[7], 徐凤敏等转化半定规划的最优性条件, 提出了一种新的投影算法[8].本文把半定规划问题转化为等价的变分不等式, 在变分不等式满足单调性和Lipschitz连续的前提下, 提出了一种邻近外梯度算法求解半定规划问题, 通过数值实验验证本文的邻近外梯度算法对于半定规划问题的求解是可行的.
考虑如式(2.1) 所示的标准形式的半定规划问题和如式(2.2) 所示的其对偶问题.
其中矩阵 $\boldsymbol{C}, {\boldsymbol{A}_i} \in {R^{n \times n}}$, $\boldsymbol{X} \in {{\bf{S}}^{n \times n}}$表示对称矩阵, $\boldsymbol{X}\underline \succ 0$表示 $\boldsymbol{X}$为半正定矩阵, “ $\cdot$”表示Frobenius内积, $\boldsymbol{C}\boldsymbol{X} = {\rm Tr}\left( {\boldsymbol{C}\boldsymbol{X}} \right) = \sum\limits_{i, j = 1}^{ n} {{\boldsymbol{C}_{ij}}{\boldsymbol{X}_{ij}}} $表示矩阵 $\boldsymbol{C}\boldsymbol{X}$的迹,
为给定的 $m$维向量.
半定规划原问题(2.1) 和其对偶问题(2.2) 的最优性充分条件如下:
转化最优性条件(2.3) 为变分不等式, 利用邻近外梯度算法对SDP问题进行求解.
定义2.1[9] 设 $\boldsymbol{C}$为有限维欧氏空间中的非空闭凸子集, 映射 $F:\boldsymbol{C} \to {R^n}$, 若对任意的 $x \in \boldsymbol{C}$, 寻求 ${x^*} \in \boldsymbol{C}$使得下式恒成立
称式(2.5) 为变分不等式问题, 记为VIP $(\Omega, F)$.
Noor在文献[10]中把上述变分不等式的经典定义推广到实Hilbert空间, 给出如下定义.
定义2.2[10] 设H为实Hilbert空间, 其内积记为 $\left\langle { \cdot, \cdot } \right\rangle $, 范数记为 $\left\| { \cdot } \right\|$, $\Omega\subset H$为一非空闭凸集, 映射 $F:H \to H$, 对任意的 $v\in\Omega$, 求 $u\in\Omega$, 使下式恒成立
定义2.3 设H为实Hilbert空间, $F:H \to H$为一映射,
(1) 若对任意的 $u, v \in \Omega $, 都有 $\left\langle {F(u) - F(v), v - u} \right\rangle \ge 0$成立, 则称 $F$在 $\Omega$上是单调的.
(2) 若对任意的 $u, v \in \Omega $, 存在 $\tau>0$, 使得
成立, 则称 $F$在 $\Omega$上强单调, 其中 $\tau$被称为映射 $F$的强单调常数.
(3) 若对任意的 $u, v \in \Omega $, 存在 $L>0$, 使得
则称 $F$在 $\Omega$上Lipschitz连续, 其中 $L$为常数, 称其为映射 $F$的Lipschitz常数.
定义2.4 设 $\Omega$为实Hilbert空间 $H$上一非空闭凸集, 记为 $\Omega\subset H$, 给定其中一点 $u\subset H$, 如果存在 $z\subset\Omega$, 使得
则称 $z$是 $u$在 $\Omega$上的投影, 记为 $z = {P_\Omega }[u]$.
考虑半定规划原问题(2.l)的对偶问题(2.2), 去掉其松弛矩阵 $\boldsymbol{S}$可写为
从而(2.3) 式的最优性条件可表示为
定义有限维欧式空间 ${R^{n \times n}} \times {R^n}$上的内积和范数, 若令 $u = {(X, y)^{\mathop{\rm T}\nolimits} }$, 则对任意的 ${u_1} = {({X_1}, {y_1})^{\mathop{\rm T}\nolimits} }$, ${u_2} = {({X_2}, {y_2})^{\mathop{\rm T}\nolimits} }$, 定义内积为 $\left\langle {{u_1}, {u_2}} \right\rangle = {X_1} \cdot {X_2} + y_1^{\mathop{\rm T}\nolimits} {y_2}$, 则由此内积诱导出的范数为 $\left\| u \right\| = \sqrt {X \cdot X + {y^T}y} $, 用vec $(X)$表示矩阵 $X$的列依次按列的顺序尾首连接成的 $n^2$维向量.
由矩阵论[11]理论易知, 任一个对称矩阵 $\boldsymbol{A}$可以写成 $\boldsymbol{\boldsymbol{A}} = {\boldsymbol{W}^{\mathop{\rm T}\nolimits} }\boldsymbol{D}\boldsymbol{W}$的形式, 其中 $\boldsymbol{W}$为正交矩阵, $\boldsymbol{D}$为对角阵.定义 ${d_i} = {\boldsymbol{D}_{ii}}$, $d_i^ + = \max \left\{ {{d_i}, 0} \right\}$, $d_i^ - = \min \left\{ {{d_i}, 0} \right\}$, 而且 $d_i^ +$, $d_i^ - $分别为对角阵 ${\boldsymbol{D}^ + }$, ${\boldsymbol{D}^ - }$的对角元素, 则
并且 $P_1[\boldsymbol{A}]=\boldsymbol{W}^T\boldsymbol{D}^+\boldsymbol{W} $表示对称矩阵 $\boldsymbol{A}$在 $S^n_+$上的正交投影, $P_2[\boldsymbol{A}]=\boldsymbol{W}^T\boldsymbol{D}^-\boldsymbol{W}$表示对称矩阵 $\boldsymbol{A}$在 $PS^n_+$上的正交投影[12].在凸紧集上的投影具有一个重要的性质[13]:对于任意的 $v \in {R^{n \times n}} \times {R^m}$和 $w\in \Omega$, 有
下面给出的引理和定理建立了互补松弛条件和投影方程之间的等价关系.
引理2.1 $\boldsymbol{A}\in \boldsymbol{S}^n, \boldsymbol{B}\in \boldsymbol{S}_+^n$, 如果 $\langle \boldsymbol{A}, \boldsymbol{B}\rangle=0$当且仅当 $\boldsymbol{A}\boldsymbol{B}=0$.
定理2.1[14] $\boldsymbol{X}^0\in \boldsymbol{S}_+^n$, $\boldsymbol{S}^0\in \boldsymbol{S}_+^n $, 若 $\boldsymbol{X}^0\boldsymbol{S}^0=0$当且仅当
证 必要性:令 $\boldsymbol{X}=\boldsymbol{X}^0-\boldsymbol{S}^0 $, 对于 $\forall \boldsymbol{X}^T=\boldsymbol{X} $, 则在如下的不等式
得到
因为 $P_{\boldsymbol{S}_+^n}[\boldsymbol{X}^0-\boldsymbol{S}^0]\in \boldsymbol{S}^n_+$, 则满足 $\boldsymbol{X}^0\boldsymbol{S}^0=0$.
由引理2.1可知
又由 $\langle P_ {\boldsymbol{S}_+^n}[\boldsymbol{X}^0-\boldsymbol{S}^0], \boldsymbol{S}^0\rangle \geq0$得到
通过不等式(2.9) 和(2.10), 得到 $E(\boldsymbol{X}^0)=0$.
充分性:由 $E(\boldsymbol{X}^0)=0$, 得到 ${\boldsymbol{X}^0} = {P_{\boldsymbol{S}_ + ^n}}\left[{{\boldsymbol{X}^0}-{\boldsymbol{S}^0}} \right] \in \boldsymbol{S}_ + ^n$.由
成立, 从而得到
记
$\forall \boldsymbol{X}\in \boldsymbol{S}^n_+ $均成立, 则在这个不等式中若令 $\boldsymbol{X}^0=P_{\boldsymbol{S}_+^n}[\boldsymbol{X}^0-\boldsymbol{S}^0]$, 得到 $\forall \boldsymbol{X}\in \boldsymbol{S}_+^n$,
在不等式(2.13) 中, 用 $\boldsymbol{X}=0\in \boldsymbol{S}_+^n$和 $\boldsymbol{X}=2\boldsymbol{X}^0 \in \boldsymbol{S}_+^n$替换 $\boldsymbol{X}$得到
于是由引理2.1得到 $\boldsymbol{X}^0\boldsymbol{S}^0=0$.
定理2.2 设 $u=(\boldsymbol{X}, y)^T\in \Omega$, 那么最优性条件(2.7) 与如下的投影方程等价
其中 $ Mu = \left( {\begin{array}{*{20}{c}} { - \sum\limits_{i = 1}^m {{\boldsymbol{y}_i}{\boldsymbol{A}_i}} } \\ {{\boldsymbol{A}_1} \cdot \boldsymbol{X}} \\ \vdots \\ {{\boldsymbol{A}_m} \cdot \boldsymbol{X}} \\ \end{array}} \right)$, $ \boldsymbol{q} = \left( {\begin{array}{*{20}{c}} \boldsymbol{C} \\ { - \boldsymbol{b}} \\ \end{array}} \right)$, ${\Omega ^ + } = \boldsymbol{S}_ + ^n \times {R^m}$.
证 最优性条件(2.3) 等价于
于是SDP的最优性条件等价于式(2.15) 的变分不等式问题.所以求解半定规划最优解的问题, 可以通过求解变分不等式(2.15) 实现, 能够解决该变分不等式的方法, 都可能同时获得SDP的最优解.而在本文中提出邻近外梯度算法解决该问题.
有许多求解变分不等式的算法和研究成果, 如投影法[15]、Wiener-Hopf方程法、逐点逼近法等[16-18].将半定规划问题转化为变分不等式的等价形式后, 许多学者利用投影一收缩算法对其进行求解, 但是在求解时, 每次迭代都需要计算两次投影, 计算效果并不是很理想: $VI(\Omega, F)$的外梯度法最早是由Korpelevich提出, 何炳生教授为了提高投影算法的求解效率做了许多研究工作[19], 通过每次迭代中增加一个投影来克服一般投影算法限制太强的缺点.并且, 何炳生教授在2004年提出了一种近似邻近外梯度类型算法[20].本文针对迭代步骤中不规则闭凸区域上投影计算难的问题, 结合近似邻近点算法, 外梯度算法以及次梯度半空间投影算法, 提出邻近外梯度算法求解半定规划问题, 在每次迭代步中定义误差限 $\varepsilon$, 保证算法具有很好的收敛性.
通过最优性条件(2.15), SDP问题可等价转化为一个定义在有限维欧氏空间中, 并在闭凸集 $\Omega=\boldsymbol{S}^n_+\times R^m$上的变分不等式问题
其中
由 $F(u)$的定义可知 $F(u)$单调且Lipschitz连续[21, 22].本文结合外梯度算法, 构造包含原投影区域的半空间, 将投影建立在半空间上, 简化投影的求解过程, 并对新的邻近点序列作相应限制.下面给出邻近外梯度算法具体的实现过程.
步骤1 给定 $\beta_0=1$, $\varepsilon\ge 0$, $0<u<v<1$, $0<\rho<1$, $0<\gamma<2$, $u^0 \in \Omega$, $k=0$.
步骤2 计算预测点
停, 否则转步骤3.
步骤3 令 ${r_k} = {\beta _k}\left\| {F\left( {{u^k}-{{\bar u}^{k + 1}}} \right)} \right\|{\rm{/}}\left\| {e\left( {{u^k}, {\beta _k}} \right)} \right\|$, 若 $r_k\leq v$, 计算
步骤4 构造半空间 $H_k$,
计算在半空间上 $H_k$投影
若 $\frac{{{\alpha _k}F({y^k}) - {\varepsilon ^k}}}{{{\beta _k}F({x^k})}} > \frac{{{\tau _k}}}{{{\beta _k}L}}$, 修正
计算下一个迭代点
步骤5 若 $e(F({u^k}), {\beta _k}) = 0$, 计算
否则 $k=k+1$, 转步骤1.
教学测试稳定性的统计研究中有一类重要问题即教育测试问题(ETP).为了验证算法的可行性, 在Matlab2012中执行该算法, 在Intel(R), Pentium(R), Core(TM) i7-3520M CPU, 2.9GHz, 4.00GB内存, Windows 7操作系统上进行数值实验.在实施算法时, 为了与外梯度算法和半定规划软件包SDP packuser内点算法的计算结果进行比较[23-24], 实验参数的设置相同.教育测试问题模型如下.
给定一个对称正定矩阵 $\boldsymbol{C}$, 求
其中 ${\rm Diag}(\boldsymbol{y})$表示由向量 $\boldsymbol{y}$为对角元素得到的对角矩阵, $\boldsymbol{e}$表示元素全为1的列向量. Iter表示求解问题过程中算法得到最终结果所需要的迭代次数, CPU-time表示完成计算所需要的时间, 单位 $s$, 问题规模从 $n=4$到 $n=180$, $\varepsilon=10^{-4}$为停止准则.下表给出了求解SDP问题的迭代次数和CPU计算时间的结果.
表 1和2分别给出了求解半定规划问题一次计算的迭代次数和CPU计算时间, 与解半定规划的内点法相比, 本文算法编程容易, 单步计算量小, 能充分利用原始数据的稀疏性, 占用内存小.可以看出, 在精度要求相同的情况下, 本文提出的邻近外梯度算法在求解半定规划问题时运行速度比内点算法快得多, 迭代次数明显减少, 算法运行时间短.与外梯度算法相比, 本文算法的迭代次数和运行时间更优.从数值实验结果可以看出, 提出的算法对于半定规划问题的求解是可行的.
本文通过把求解半定规划问题有效的转化为求解变分不等式问题, 把变分不等式的数值解法从欧式空间推广到有限维空间, 给出了一个求解半定规划问题的邻近外梯度算法.并在算法中考虑不规则闭凸区域上投影计算难的问题来构造半空间上的投影, 在算法中引入特殊的迭代步长因子, 增加了迭代点的灵活性, 通过数值实验说明该算法是求解半定规划问题的一种可行方法.半定规划与许多优化问题有着紧密的联系, 将半定规划和互补问题结合是今后的研究重点.