数理解析研究所講究録
第 1773 巻 2012 年 186-197

186

0-1 整数変数を含む非凸 2 次最適化問題の非負半正定値緩和に対する面的縮小と効率的解法
中田 和秀 \d ag er

田中 未来 \dag er*

脇 隼人 \S概要

本論文では 0-1 整数変数を含む非凸 2 次最適化問題に対する非負半正定値緩和問題を効率よく解く方法を提案する．始めに，非負半正定値緩和問題が退化していることを用いて問題を縮小する方法を提案する．この縮小は面的縮小の部分的な適用とみなすことができるため，主双対パス追跡法の安定性が高まることが期待できる．続いて，縮小後の問題を主双対パス追跡法で解
く．探索方向の計算では，大規模な線形方程式系を Krylov 部分空間法の一種である PSQMR 法を用いて解くことで効率化を図る．さらに，PSQMR 法の収束性を高めるために 3 つの前処理を提案する．最後に，数値実験により提案手法の効率性を示す．

1

はじめに
本論文では次の 0-1 整数変数を含む非凸 2 次最適化問題について考える:

minimizesubject to$x^{T}\Gamma x+2\gamma^{T}x$$(i=1, \ldots, m)$ ,$\alpha_{i}^{T}x=\beta_{i}$$x\geq 0$,$(j\in B)$$x_{j}\in\{0,1\}$

ここで$x\in \mathbb{R}^{n}$$B\subset\{1, \ldots, n\}$は決定変数，$A=(\alpha_{1}, \ldots, \alpha_{m})^{T}\in \mathbb{R}^{m\cross n},$は与えられたデータである．なお，$S^{n}$(MBNQO).

$\beta=\phi_{1},$ $\ldots,\beta_{m})^{T}\in \mathbb{R}^{m},$$\gamma\in \mathbb{R}^{n},$$\Gamma\in S^{n}$,は $n$ 次対称行列の空間を表す．この問題は幅広

い応用範囲をもつ [2, 7]. しかし，この問題は離散構造と非凸性を有する問題であり，(MBNQO) の大域的最適解を求めることは NP 困難である．この問題に対して Burer [71 は，比較的ゆるい仮定の下で (MBNQO) の最適値が次の完全正値最適化問題の最適値と一致することを示している:

minimizesubject to$\Gamma\cdot X+2\gamma^{T}x$$\alpha_{i}^{T}x=\beta_{i}$$(i=1, \ldots, m)$ ,$\alpha_{i}^{T}X\alpha_{i}=\beta_{i}^{2}$$(i=1, \ldots, m)$ ,$x_{j}=X_{jj}$$(j\in B)$ ,$(\begin{array}{ll}1 x^{T}x X\end{array})\in(C^{1+n})^{*}$\dag er.

本研究は科研費 (22710136), 科研費 (22740056) の助成を受けたものである．
東京工業大学大学院社会理工学研究科経営工学専攻

1 東京工業大学大学院社会理工学研究科経営工学専攻\S 電気通信大学大学院情報理工学研究科情報通信工学専攻

(CPO)

187
ここで

$x\in \mathbb{R}^{n},$$X\in S^{n}$は決定変数，．は行列の内積，すなわち$\Gamma\cdot X=tr(\Gamma^{T}X)$を表し，$(C^{n})^{*}$

は

$n$

次

完全正値行列のなす錐 [1,4] を表す．(CPO) は閉凸錐上の線形最適化問題だが，最適解を求めることは依然として難しい．実際，与えられた行列が完全正値かどうかを判定することは NP 困難であることが Dickinson and Gijben [9] によって示されている．

Burer [7] は (CPO) の最適値に対する強い下界を求めるために非負半正定値緩和を提案している．
(CPO) に対する非負半正定値緩和問題は (CPO) における完全正値制約 $(1, x^{T};x, X)\in(C^{1+n})^{*}$ を非負半正定値制約 $(1, x^{T};x, X)\in S_{+}^{1+n}\cap S_{(+)}^{1+n}$ で置き換えることによって得られる．ここで鍔と
はそれぞれ

$n$

次半正定値対称行列のなす錐と$n$$S_{(+)}^{n}$次非負対称行列のなす錐であり，$(A, B;C, D)=(\begin{array}{ll}A BC D\end{array})$である．この非負半正定値緩和は単純な半正定値緩和よりも強い緩和となることが Ge and Ye[12]によって指摘されている．また，非負半正定値緩和問題は半正定値最適化問題に帰着することができるため，SeDuMi [20], SDPA [27], SDPT3 [24] のような半正定値最適化問題に対する既存のソルバを用いて解くことができる．

しかし，このアプローチには次の 2 つの問題点がある．1 つめは，非負半正定値緩和問題は実行可能内点解をもたないため，主双対内点法の数値的安定性が悪くなることである．2 つめは，中規模の非負半正定値緩和問題であっても，等価な半正定値最適化問題が非常に大規模な問題となることである．実際，非負半正定値緩和問題と等価な半正定値最適化問題は $(n+1)(n+2)/2+2m+|B|+1$ 個の線形等式制約をもっため，通常の主双対内点法の各反復では $(n+1)(n+2)/2+2m+|B|+1$ 次の線形方程式系を解くことになる．本論文ではこれら 2 つの問題点を解決する方法を提案する．
提案手法では，はじめに非負半正定値最適化問題が実行可能内点解をもたないことを用いて問題のサイズを縮小する．この縮小法は，面的縮小 [5,6,17,18,19,26] の部分的な適用とみなすことができる．本論文で提案する縮小法は不完全な面的縮小であるため，実行可能内点解を持つとは限らないが，主双対内点法の安定性が向上することが期待できる．

縮小後の非負半正定値緩和問題であっても等価な半正定値最適化問題は大規模になるため，通常の主双対内点法を用いて解く際には膨大な計算時間を要する．主双対内点法において特に計算時間のかかる計算は探索方向の計算である．実際，縮小後の非負半正定値緩和問題と等価な半正定値最適化問題は $(n+1)(n+2)/2+|B|+1$ 本の線形等式制約を持つため，各反復で $(n+1)(n+2)/2+|B|+1$次元の線形方程式系を解く必要がある．

さらに本論文では，縮小後の非負半正定値緩和問題を効率よく解く方法を提案する．提案手法では Toh [23] によって提案された凸 2 次半正定値最適化問題に対する主双対パス追跡法に基づき，縮小後の非負半正定値緩和問題を半正定値最適化問題に帰着させることなく主双対パス追跡法で解

く．各反復における探索方向の計算では，$(n-m+1)^{2}+|B|+1$ 次元の線形方程式系を$\infty 1ov$

部分空間法の一種である PSQMR 法 [11] を用いて解く．Krylov 部分空間法の収束性は悪条件な線形方程式系を解く場合に悪化することが知られている．今回の場合，主双対パス追跡法の反復点が最適解に近づくにつれて解くべき線形方程式系の係数行列は悪条件になる．そこで本研究では，PSQMR 法の収束性を高めるために Toh [23] によって提案された 3 つの前処理を今回の場合に拡張する．

188
本論文は次のような構成からなる．2 節では，非負半正定値緩和問題のサイズを縮小する方法を提案する．3 節では，縮小後の問題に対する主双対パス追跡法を提案する．4 節では，提案手法の効率性を示す数値実験の結果を示す．最後に 5 節で総括を行う．

に対し，Kronecker
以下，本論文で用いる記号の定義を行う．2 つの行列を積
に対し，
: $X\mapsto QXP^{T}$ と定義する．また，行列$P\in \mathbb{R}^{n\cross n},$$P\otimes Q:\mathbb{R}^{m\cross n}arrow \mathbb{R}^{m\cross n}$

diag$(X)\in \mathbb{R}^{n}$$Q\in \mathbb{R}^{m\cross m}$$X\in \mathbb{R}^{n\cross n}$$P\otimes Q$

で $X$ の対角成分からなるベクトルを表す．逆に，ベクトル$x\in \mathbb{R}^{n}$に対し，Diag $(x)\in S^{n}$である対角行列を表す．さらに，行列に対し，線形写像 Diag $(P)$ :を Diag $(P)$ : $X\mapsto P\circ X$ と定義する．ここで は Hadamard 積を表す．さらに，行列で対角成分が$P\in \mathbb{R}^{m\cross n}$$x$$X\in \mathbb{R}^{m\cross n}$$\circ$$p\in \mathbb{R}$に対し，$X^{(p)}\in \mathbb{R}^{m\cross n}$

は $(i, j)$ 成分が$X_{ij}^{p}$$\mathbb{R}^{m\cross n}$とである行列を表す．

2 面的縮小
この節では (CPO) に対する次の非負半正定値緩和問題について考える:

minimize$\gamma 0$

subject to$01$$\gamma^{T}\Gamma)\cdot(\begin{array}{ll}\xi x^{T}x X\end{array})$.

$0^{T}O)\cdot(\begin{array}{ll}\xi x^{T}x X\end{array})=1$,$(\begin{array}{ll}\xi x^{T}x X\end{array})=2\beta_{i}$$(i=1, \ldots,m)$ ,$(\begin{array}{ll}0 0^{T}0 \alpha_{i}\alpha_{i}^{T}\end{array})\cdot(\begin{array}{ll}\xi x^{T}x X\end{array})=\beta_{i}^{2}$$(i=1, \ldots, m)$ ,$(_{-2e_{j}}1$$(j\in B)$ ,$(_{\alpha_{j}}^{0}$$\alpha_{O^{i)}}^{T}$(DNR)$4e_{j}e_{j}^{\dot{*})\circ}-2e^{T}(\begin{array}{ll}\xi x^{T}x X\end{array})=1$$Y=(\begin{array}{ll}\xi x^{T}x X\end{array})\in S_{+}^{1+n}\cap S_{(+)}^{1+n}$ここで$e_{j}\in \mathbb{R}^{n}$

は第

$j$

成分が 1 でそれ以外の成分が.

のベクトルである．この問題に対し，正定値か$0$つすべての成分が正であるような実行可能解を実行可能内点解という．

この問題に対し，$R=(-\beta, A)^{T}$ とおき，$S=RR^{T}$ とすると，(DNR) の任意の実行可能解$Y\cdot S=0$が成り立っ．これと

$S$

$Y$

に対し，

の半正定値性より， が正定値でないことがわかる．したがって，$Y$(DNR) は実行可能内点解をもたない．なお，(DNR) が実行可能内点解をもたないことを示す際には$Y$の半正定値性のみを用いている．そのため，(CPO) における$(C^{1+n})^{*}$

を$S_{+}^{1+n}$に含まれる錐で置き換えた問題も実行可能内点解をもたないことを示すことができる．また，Burer[71 は (CPO) のサイ$m=1$ の場合に実行可能内点解が出現する例を示している．しズを縮小する方法を提案しており，$m\geq 2$ の場合はこの方法を適用しても実行可能内点解が出現しないことが Tanaka et al. [21,かし，Appendix$A$] により示されている．

Kobayashi et al. [13] が指摘しているように，実行可能内点解をもたない錐最適化問題を主双対内点法で解く場合，数値的に不安定になる．そこで，この$S$を用いて (DNR) のサイズを縮小する．この縮小法は Waki- and Muramatsu [26, Section 4] で提案された面的縮小に基づく．錐最適化問題に対す

る面的縮小は有限回の反復で元の問題と等価で実行可能内点解をもっ問題を出力する．一般的な面

189
的縮小についての詳細は [5,6,17,18,19,26] を参照されたい．本論文で提案する縮小法はこの面的縮小の 1 反復に相当するため，提案する縮小法を適用して得られる非負半正定値緩和問題は実行可能内点解をもっとは限らないが，主双対内点法の数値的安定性を高めることができると期待できる．
そのために $V=(L, R)$ を正則にする $L\in \mathbb{R}^{(1+n)\cross(1+n-m)}$ をとる．この $L$ のとり方には自由度があるが，以下では Tanaka et al. [21, Section 2.1] で提案されている を用いた を用いて (DNR) のサ
イズを縮小することを考える．いま，$V^{T}YV=(L^{T}YL, L^{T}YR;R^{T}YL, R^{T}YR)$ に対し，(DNR) の任意$R^{T}YR=O_{m\cross m}$ が成り立つ．このこととの実行可能解 に対して $Y\cdot S=0$ が成り立っことから，$V^{T}YV$の半正定値性より，$L^{T}YR=O_{(1+n-m)\cross m}$ が成り立っ．このことを用いると，(DNR) と等価な次
$L$

$V$

$Y$

の問題を得る:mini 而 ze$V^{-1}(\begin{array}{ll}0 \gamma^{T}\gamma \Gamma\end{array})V^{-T}\cdot V^{-1}YV$

subject to$V^{-1}(\begin{array}{ll}1 0^{T}0 O\end{array})V^{-T}\cdot V^{T}YV=1$,$V^{-1}(\begin{array}{ll}0 \alpha_{i}^{T}\alpha_{i} O\end{array})V^{-T}\cdot V^{T}YV=2\beta_{i}$$(i=1, \ldots,m)$ ,$V^{-1}(\begin{array}{ll}0 0^{T}0 \alpha_{i}\alpha_{i}^{T}\end{array})V^{-T}\cdot V^{T}YV=\beta_{i}^{2}$$(i=1, \ldots, m)$ ,$V^{-1}$$4eje_{j}-2e^{T}\dotplus)V^{-T}\cdot V^{T}YV=1$$(_{-21}1_{e}$$V^{T}YV=(\begin{array}{ll}X O_{(1+n-m)\cross m}O_{m\cross(1+n-m)} O_{m\cross m}\end{array})$$X\in S_{+}^{1+n-m}$

さらに，$U=(I_{1+n-m}, O_{(1+n-m)\cross m})V^{-1}$,$Y\in S_{(+)}^{1+n}$(DNR2)$(j\in B)$ ,,.

, $C=U(O, \gamma^{T};\gamma,\Gamma)U^{T},$ $A_{j}=U(1, -2e_{j}^{T};-2e_{j}, 4e_{j}e_{j}^{T})U^{T}(j\in B)$

とおくことにより，(DNR2) のサイズを次のように縮小することができる:
mlnlmlze$C\cdot X$

subject to$X_{11}=1$ ,$(j\in B)$ ,$A_{j}\cdot X=1$

(FRDNR)$Y=U^{T}XU$ ,$X\in S_{+}^{1+n-m}$,$Y\in S_{(+)}^{1+n}$.

ここで (FR-DNR) のすべての係数行列のサイズが $(1 +n-m)\cross(1+n-m)$と縮小されていることに注意する．また，この縮小により (DNR2) の 2 番目と 3 番目の線形等式制約に対応する線形等式制約が冗長になるため，(FRDNR) ではそれらを削除している．詳細については Tanaka et al. [21,Section 2.1] を参照されたい．

3 主双対パス追跡法
この節では，前節で述べた面的縮小を適用した非負半正定値最適化問題の主問題:

minimizesubject to

Ce

$X$

$AX=b$ ,$Y-U^{T}XU=O$ ,$X\in S_{+}^{n-r}$,$Y\in S_{(+)}^{n}$

(P)

190
とその双対問題:

maximizesubject to$b^{T}y$$\ovalbox{\tt\small REJECT}^{T}y+S+UTU^{T}=C$$S\in S_{+}^{n-r}$

,

(D)

,

$T\in S_{(+)}^{n}$を同時に解く主双対パス追跡法を提案する．ここで，fl : $X\mapsto(A_{1}\cdot X, \ldots, A_{m}\cdot X)^{T}$ は与えられた線形写像であり， は図の随伴写像，すなわち詔 :$T$$\ovalbox{\tt\small REJECT}^{T}$$y \mapsto\sum_{i=1}^{m}A$助である．

まず，(P) と (D) に対する中心パスを次のように定義する:$\{(X, Y,y,S, T):\ovalbox{\tt\small REJECT} X=b,Y-U^{T}XU=O\ovalbox{\tt\small REJECT}^{T}y+S+UTU^{T}=C\exists\mu>0’$$XS=\mu wIX\in S_{++}^{n-r}S\in S_{++}^{n-r},$$Y\circ T=\mu ET\in S_{(++)}^{n}Y\in S_{(++)}^{n},$’$\}\cdot$,
ここで$S_{++}^{n}$と$S_{(++)}^{n}$はそれぞれ$S_{+}^{n}$と$w>0$ は重み付けパラメータ (定数)の内部を表す．また，$S_{(+)}^{n}$である．中心パスは滑らかな曲線となっており，$\mu\downarrow 0$のとき最適解の 1 つに収束する．そのため，この曲線を数値的に追跡することにより最適解を得ることができる．
続いて，(P) と (D) に対する主双対パス追跡法の概要を述べる．このアルゴリズムでは次のように中心パスを数値的に追跡する:(Step ) 初期点 $(X, Y,y,S, T)\in S$ 紅 r$0$$\cross$

Sn

$(++$

Rm Sn X Sn$\cross$$\cross$$)\cross$

$(++$

$+$

$)$

とパラメータ $\sigma\in(0,1),$ $\gamma\in(0,1)$を定める．

(Step 1) 現在の反復点 $(X, Y,y,S, T)$ が最適解に十分近ければ終了する．

(Step 2) スケーリングされた Newton 方程式を近似的に解き，最適解に近い中心パス上の点に近づく探索方向$(\Delta X, \Delta Y, \Delta y, \Delta S, \Delta T)$を求める．
$:= \min\{\max\{\alpha$:$:= \min\{\max\{\alpha:S+\alpha\Delta S\in S_{+}^{n-r}\},$ $\max\{\alpha$:(Step 3) 次の反復点が錐に入るような最大のステップサイズ，すなわち$X+\alpha\Delta X\in S_{+}^{n-r}\},$ $\max\{\alpha:Y+\alpha\Delta Y\in S_{(+)}^{n}\}\},$$T+\alpha\Delta T\in S_{t+)}^{n}\}\}$を求める．

(Step 4) ステップサイズを$\alpha_{p}\Delta X,$$\alpha_{d}$$\alpha_{p}$$\alpha_{p}$$:= \min\{1,\gamma\alpha_{p}\},$$\alpha_{d}$$:= \min\{1, \gamma\alpha_{d}\}$$Y+\alpha_{p}\Delta Y,y+\alpha_{d}\Delta y,S+\alpha_{d}\Delta S,$ $T+\alpha_{d}\Delta T)$$(X, Y,y,S, T)$ $:=(X+$とし，と反復点を更新し，(Step 1) に戻る．

一般的な主双対パス追跡法では，(Step2) においてスケーリングされた Newton 方程式を縮小した線形方程式系を直接法を用いて正確に解く．しかし今回の場合，計算量の観点から直接法を適用することは困難である．そのため，(Step2) ではスケーリングされた Newton 方程式を縮小した拡大方程式系と呼ばれる線形方程式系を解く．この線形方程式系の次元は依然として大きいため，直接法では現実的な計算量で解くことができない．そのため，反復法を用いて計算量を削減することを試みる．この方法は，Toh [231 によって提案された凸 2 次半正定値最適化問題に対する主双対パス
追跡法に基づく．

本論文では，前処理の構築しやすさの観点から探索方向として NT (Nesterov-Todd) 方向 [22] を考える．現在の反復点 $(X, Y,y,S, T)$ に対し，探索方向$(\Delta X, \Delta Y,\Delta y, \Delta S, \Delta T)$の計算は次の拡大方程式系の求解に帰着できる:$(\begin{array}{ll}-\prime H \ovalbox{\tt\small REJECT}^{T}\ovalbox{\tt\small REJECT} O\end{array})(\begin{array}{l}\Delta X\Delta y\end{array})=(\begin{array}{l}R_{a}r_{p}\end{array})$

.

(AS)

191
ここで，
${}^{t}H=W^{-1}\otimes W^{-1}+(U\otimes U)Diag(Y^{t-1)}\circ T)(U\otimes U)^{T}$,$R_{a}=R_{d}-R_{psd}-U(Y^{(-1)}\circ R_{nn})U^{T}+U(Y^{(-1)}\circ T\circ R_{p})U^{T}$である．

$H^{-1}$の計算に膨大な計算量を要することから，(AS) から$\Delta X$を消去することは難しい．そのため，(AS) を解くことにより探索方向を計算する．(P) と等価な半正定値最適化問題を一般的な$n(n+1)/2+m$ 次元の Schur 方程式を解く．一方，提案手法で解く
主双対内点法を用いて解く場合，(AS) の次元は $(n-r)^{2}+m$ と小さくなっていることに注意する．
$O([(n-r)^{2}+m]^{3})$ の計算時間と $O([(n-r)^{2}+m]^{2})$ の記憶領域を(AS) を直接法で解くためには，要する．そのため，中規模以上の問題を解く際に直接法を用いることは現実的ではない．したがって，(AS) を反復法を用いて解くことを考える．Toh [23] は PSQMR 法 [11] を用いると，(AS) と構造の近い線形方程式系を効率良く解くことができることを示している．そのため，以下では PSQMR 法を用いることとする．

PSQMR 法は Freund and Nachtigal [11] [こよって提案された Krylov 部分空間法の 1 つである．各反復で計算する係数行列とベクトルの積は次のように効率よく計算できる:.

$(\begin{array}{ll}-\prime H \ovalbox{\tt\small REJECT}^{T}\ovalbox{\tt\small REJECT} t O\end{array})(\begin{array}{l}Qq\end{array})=(\begin{array}{ll}-W^{-1}QW^{-1}-U[Y^{(- 1)}\circ T\circ(U^{T}QU)]U^{T}+\ovalbox{\tt\small REJECT}^{T}q\ovalbox{\tt\small REJECT} Q \end{array})$

この右辺は $O(n^{2}(n-r)+m(n-r)^{2})$ の計算時間と $O(n^{2}+m)$ の記憶領域で計算することができる．

したがって，PSQMR 法の反復回数を抑えることができれば，(AS) を少ない記憶領域で高速に解くことができる．

反復法を用いる欠点は (AS) の残差が O にならないことにある．提案手法では，このことが主双対

パス追跡法の収束性に支障をきたさないように Kojima et al. [14] が提案している補正を施した．詳
細については Tanaka et al. [21, Section 3. 1] を参照されたい．

次に，PSQMR 法の収束性を高めるために (AS) に対する 3 つの前処理を提案する．拡大方程式系の係数行列は，主双対パス追跡法の反復点が最適解に近づくにつれて悪条件になる．一般に，Krylov
部分空間法の収束性は係数行列が悪条件であるときに悪化するため，拡大方程式系の係数行列に対して適切な前処理を行うことが不可欠である．Toh [231 は (AS) と構造の近い線形方程式系に対す
る 3 つの効果的な前処理を提案している．以下で提案する前処理は，これらの前処理の$rV^{*},$$\Psi^{*},$$\Sigma$拡張である．

以下で提案する前処理はすべて$(-\prime K,\ovalbox{\tt\small REJECT}^{T};\ovalbox{\tt\small REJECT}, O)$

易な賀の近似である．

$(R;r)\in S^{n-r}\cross \mathbb{R}^{m}$という形をもつ．ここで穴はに対し，$(Q;q)=(-(K,ffl^{T};\ovalbox{\tt\small REJECT}, O)^{-1}(R;r)$$K^{-1}$の計算が容は簡単に計算できる．以下では， の構成方法の概要について述べる．前処理についての詳細は Tanaka et al. [21,Section 3.21 を参照されたい．
$’\kappa$$\blacksquare$

Kro ecke 積近似
$n$

$r$

によって提案された$K_{1},$$\ldots,$$K_{q}\in R^{m\cross m}$

まず，Kronecker 積近似に基づく方法にっいて述べる．この方法は Toh [23]$\gamma$の拡張である．K-ronecker 積近似とは，与えられた行列に対し，次のように$q$$G_{1},$$\ldots,G_{q}\in \mathbb{R}^{n\cross n}$,個の K-ronecker 積の和を 1 個の Kronecker 積で近似する

192
ことである:$\min_{v_{\text{し}}\in R^{n\cross n},V_{R}\in R^{m\cross m}}\Vert\sum_{p=1}^{q}G_{p}\otimes K_{p}-V_{L}\otimes V_{R}\Vert_{F}^{2}$.

Langville and Stewart [15] と Toh et al. [251 はこの問題の最適解が効率よく計算できることを示している．この近似を用いるために，賀を Kronecker 積の和で表す．そのために$Y^{(-1)}\circ T$

を次のように固有値分解する:$Y^{(-1)}\circ T=\sum_{i=1}^{n}\lambda_{i}q_{i}q_{i}^{T}$.

これらの固有値と固有ベクトルを用いて，次のように Kronecker 積近似を施す:$cH=W^{-1} \otimes W^{-1}+\sum_{i=1}^{n}\sqrt{\sigma_{i}\lambda_{i}}U$

ここで $\sigma=$ $(sgn(\lambda_{1}),$$\ldots$, sgn$(\lambda_{n}))^{T}$Diag$(q_{i})U^{T}\otimes\sigma_{i}\sqrt{\sigma_{i}\lambda_{i}}U$である．なお，$V$

ブロック Kronecke 積近似$(q_{i})U^{T}\simeq V\otimes V=’\kappa$$:=V_{L}^{*}=V_{R}^{*}\in S_{++}^{n-r}$積近似が存在することが Tanaka et al. [21, Appendix$\blacksquare$

Diag

$B$

.

を満たす最適な Kronecker] によって示されている．

次に，ブロック Kronecker 積近似に基づく方法について述べる．こ$r$の方法は Toh [23] によって提案されたの固有値の大きさの拡張である．この前処理ではに基づいて係数行列の (1, 1) ブロックを分割する．Toh [23] が述べているように， を最大相補解に$W^{-1}$$\Psi$$l$$W^{-1}$はおける $X$ の階数としたとき，のオーダーをもつ 1 個の小さな固有値と，$\Theta(1/\sqrt{\mu})$$\Theta(\sqrt{\mu})$

オーダーをもつ $n-r-l$ 個の大きな固有値をもつ．これに基づき$W^{-1}$$W^{-1}=PDP^{T}=P_{1}D_{1}P_{1}^{T}+P_{2}D_{2}P_{2}^{T}$

ここで $D_{1}=$ Diag $(d_{1})\in S_{++}^{l}$ と$D_{2}=$

Diag$(d_{2})\in S_{++}^{n-r-l}$固有値を対角成分にもっ対角行列であり，$P_{1}\in \mathbb{R}^{(n-r)\cross l}$に対応する$W^{-1}$を次のように分割する:

.

はそれぞれ

$W^{-1}$の小さな固有値と大きなと $P_{2}\in \mathbb{R}^{(n-r)\cross(n-r-l)}$ はそれぞれの固有ベクトルからなる行列である．すなわち $D=$ Diag$(d_{1};d_{2}),$$D_{1}$と$P=(P_{1}, P_{2})$ある．この表記を用いると，賀を次のように表すことができる:.

$\prime H=(P\otimes P)[D\otimes D+\sum_{i=1}^{n}\sqrt{\sigma_{i}\lambda_{i}}P^{T}UDiag(q_{i})U^{T}P\otimes\sigma_{i}\sqrt{\sigma_{i}t_{i}}P^{T}UDiag(q_{i})U^{T}P](P\otimes P)^{T}$

ここで，次のようにブロック対角近似を施す:$P^{T}U$Diag$(q_{i})U^{T}P\simeq(^{P_{\iota_{O}^{UDiag(q_{j})U^{T}P_{1}}}^{T}}$$P_{2}^{T}UDiag(q_{i})U^{T}P_{2}O)$$(i=1, \ldots, n)$

さらに，次のように Kronecker 積近似を施す:${}^{t}\overline{H}_{ij}^{*} \simeq D_{j}\otimes D_{i}+\sum_{i=1}^{n}Z_{kj}\otimes\sigma_{i}Z_{ki}$$(i,j=1,2)$ .

ここで，$Z_{ki}=\sqrt{\sigma_{k}\lambda_{k}}P_{i}^{T}U$

Diag$(q_{k})U^{T}P_{i}$の$(k=1, \ldots, n;i=1,2)$.

$D_{2}$で193である．これらを用いて次のように賀の第 2 因子を近似する:$\overline{\mu}*:(\begin{array}{ll}M_{11} M_{21}^{T}M_{21} M_{22}\end{array})\mapsto$$( \frac{}{H}2^{*}1(M_{21})$$[ \overline{H}_{1}^{*}(M_{21})]^{T}\frac{2}{H}2^{*}2(M_{22}))$.

ここで $M_{11}\in S^{l},$ $M_{21}\in \mathbb{R}^{(n-r-l)\cross l},$ $M_{22}\in S^{(n-r-l)}$ である．最終的に次の近似を得る$’\kappa=(P\otimes P)\text{賀^{}*}(P\otimes P)^{T}-$$\blacksquare$

対角スケーリング

によって提案された$W^{-1}=PDP^{T}$:.

最後に，一種の対角スケーリングによる方法を述べる．この方法は Toh [23]$\Phi_{\pm}$の拡張とみなすことができる．この前処理では$W^{-1}$の固有値分解を用いて賀を次のように表現する:

賀 $=(P\otimes P)[D\otimes D+(P^{T}U\otimes P^{T}U)$ Diag $(Y^{(-1)}\circ T)(P^{T}U\otimes P^{T}U)^{T}](P\otimes P)^{T}$ .

この第 2 因子を次のように対角写像で近似する:$D\otimes D+(P^{T}U\otimes P^{T}U)Diag(Y^{(-1)}\circ T)(P^{T}U\otimes P^{T}U)^{T}\simeq$

Diag [diag $(D)$ diag$(D)^{T}+\Delta$].

以上より次の近似を得る:$’\kappa=(P\otimes P)$

Diag diag $(D)diag(D)^{T}+\Delta](P\otimes P)^{T}$ .
$[$対角写像近似の簡単な計算方法については Tanaka et al. [21, Section 3.2] を参照されたい．

4 数値実験
この節では，提案手法の効率性を示すために行った数値実験の結果について述べる．そのために，
2 次多次元ナップサック問題とポートフォリオ最適化問題の非負半正定値緩和問題を以下に挙げる
3 つの解法を用いて解き，比較を行った:(Sl) (DNR) と等価な半正定値最適化問題を既存のソルバで解く．

(S2) (FRDNR) と等価な半正定値最適化問題を既存のソルバで解く．
(S3) (FRDNR) を提案手法で解く．

以下の実験は，OS が Red Hat Enterprise Linux 55, CPU が Intel Xeon 266 GHz, メモリが 24 GBである計算機上で行った．また，半正定値最適化問題に対するソルバは SDPA(version 73.1) を用い，すべて規定値のパラメータを用いた．なお，計算時間が 20,000 秒に達した時点で計算を打ち切った．
提案手法は MATLAB (version 7.11) 上で実装した．予備実験の結果，$\Psi^{*}$と$\Sigma$は現在の反復点が最適解から遠いときはあまり効果的でないことがわかったため，以下の実験では次の規則に従って前処理を行った:(Pl) すべての反復で$\gamma*$を用いる．

(P2) err6 $>1$ のとき$\gamma*$

を，err6$\leq 1$のとき

(P3) err6 $>1$ のとき$\gamma*$

を，err6$\leq 1$のとき$\Psi^{*}$$\Sigma$

を用いる．

を用いる．

194
提案手法の実装に関する詳細は Tanaka et al. [21, Section 3.31 を参照されたい．
まず，2 次多次元ナップサック問題:

maximizesubject to$\frac{1}{2}x^{T}Qx+c^{T}x$$(i=1, \ldots, m)$ ,$a_{i}^{T}x\leq b_{i}$$x\geq 0$(QMKP),$(j=1, \ldots,n)$$x_{j}\in\{0,1\}$の非負半正定値緩和問題を解いた結果を Table 1 に示す．errl, err3, err6 はそれぞれ主問題の残差，双対問題の残差，双対ギャップを DIMACS errors [16] で測ったものである．また，time は計算
時間 (秒), solver は使った解法，precond は解法として (S3) を用いた場合に用いた前処理を表す．
(S2) を用いた場合，$n=100,$ $m=20$ の場合を除いて時間切れとなった．これは，面的縮小により係数行列の疎性が失われたことによると考えられる．$n=100,$ $m=20$ の場合，(S2) を用いて得た解の双対ギャップが最も小さいが，非常に長い計算時間がかかっている．したがって，SDPA は (QMKP) の(FRDNR) と等価な半正定値最適化問題を高速に解くことはできないといえる．一方，(Sl) と (S3)を用いた場合は多くの問題で解を得ることができた． が小さいとき，(Sl) は (S3) よりも高速だっ$m$

た．しかし， が大きくなるにつれて (Sl) の計算時間が増加し，特に $n=150,$ $m=90$ のときは時間$m$切れとなった．一方，提案手法を用いた場合の計算時間は$m$に大きく依存することがなかった．また，解の精度については (S3) を用いた場合のほうが優れていた．
表1(QMKP) に対する非負半正定値緩和問題を 3 つの解法を用いて解いた結果$\overline{\frac{nmsolver/precond(S1)/-(S2)/-(S3)/(P1)(S3)/(P2)(S3)/(P3)}{1\mathfrak{d}02Qerr11.9e-079.9e-044.0e-114.4e-114.1e-11}}$

err3
err6

lOQ

5. $2e-04$

time $[\sec]$

197.76errl

2. $4e-07$
3. $7e-96$
1.5e-03$6\Theta$

err3
err6

time $[\sec]$
$150$

30

err3
$[\sec]err6$

1. $8e-15$

2. $0e-05$
606.88403.16
1. Se-101. $6e-$ lQ
1. $0e-12$
9. $9e-13$
9.le-05444.73
S. $7e-11$
5. $9e-11$
7. $4e-15$
7. $1e-15$$4.4e-\mathfrak{d}5$658.53$7\mathfrak{d}4.39$$5.\mathfrak{d}e-11$–4.8e-lS

2. $8e-05$

4. $1e-Q3$$1.7e-\mathfrak{d}5$2,756.81$1,593.8\mathfrak{d}$$2.3\mathfrak{d}8.51$—-9. $3e-10$
3 . le-14

1. $Ne-N9$
3 . Oe-14

1. $2e-09$
3 . Qe-14$4,119.3541e-05$$2,532.2263e-\mathfrak{d}5$

1,912.14

errl

time

–

$5.2e-\mathfrak{d}4$time $[\sec]$$9Q$$>2\emptyset,$

2. $7e-Q7$
3. $5e-96$

errl

1. $4e-15$$4.6e-\mathfrak{d}5$$4.4e-\emptyset 5$$\mathfrak{d}\mathfrak{d}\mathfrak{d}\overline{.0\emptyset}$$989.\Theta 7$

err3
err6

150

1. $3e-12$
2. $1e-15$
8. $3e-06$
4.5e-05
5,767.85
443.39–1. $9e-10$–9. $9e-13$$2.3e-\mathfrak{d}6$$>2\mathfrak{d},0\mathfrak{d}0^{-}-$$>20,00\mathfrak{d}\overline{.00}$$>20,0O^{-}Q.QQ$

OQ$3,2\mathfrak{d}7.9326e-05$

次に，ポートフォリオ最適化問題:

minimizesubject to$x^{T}\Sigma x-2(\mu^{T}x)^{2}$$e^{T}x=1$ ,$a_{i}^{T}x\leq b_{i}$$(i=1, \ldots, m)$ ,(Port)$x\geq 0$の非負半正定値緩和問題を解いた結果を Table 2 に示す．
$m=0$ の場合，すべての解法で同程度の精度の解を得ることができた．計算時間を比較すると，(S3) は (Sl), (S2) に比べて速いといえる．
$n=100,$ $m=40$ や $n=100,$ $m=80$ の場合，(Sl) と (S3) では解を得ることができたが (S2) は時間195
切れとなった．(S3) を用いて得た解は (Sl) を用いて得た解よりも同程度かやや高精度であり，計算時間を比較しても (S3) は (Sl) よりも速い．
$n=150,$ $m=60$ や $n=150,$ $m=120$ の場合は (S3) のみが制限時間内に解を求めることができた．(S3) によって得られた解は十分に高精度である．また，(Port) の (FRDNR) に対する前処理を比較すると，計算時間の短い (P3) が最も優れた前処理だといえる．
表2(Port) に対する非負半正定値緩和問題を 3 つの解法を用いて解いた結果$\overline{\frac{nmsolver/precond(S1)/-(S2)/--(S3)/(P1)(S3)/(P2)(S3)/(P3)}{1\emptyset\emptyset\emptyset err14.2e-121.1e-1\emptyset 1.1e-164.4e-161.7e-16}}$

err3
err6

time $[\sec]$$1\otimes Q$

errl

$4Q$

err3
err6

2.4e-l2
8.le-07$78.\emptyset 7$

2. $6e-1Q$
l.Ne-lN

1. $1e-12$
3. $9e-07$
112.852.8e-l5
9. $8e-07$
78.19

2.8e-l5
9. $8e-07$

62.36

52.23

–2.le-l5–1.6e-l4

1. $8e-15$

1.6e-l4

1. $7e-15$
1.6e-l4$9.8e-\emptyset 7$$4.1e-\emptyset 6$$9.4e-\emptyset 7$$3.6e-\emptyset 6$

2. $8e-15$$9.8e-\emptyset 7$$\frac{time[\sec]774.44>2\emptyset,\emptyset\emptyset Q.QQ175.95172.3711626}{1\emptyset\emptyset 8\emptyset err13.7e-\emptyset 9--2.6e-162.4e-162.3e-16}$err3$2.\emptyset e-1\emptyset$–l.Ne-14

1.

$\emptyset e-14$$1.\mathfrak{G}e-14$$\frac{time[\sec]2,247.93>2\emptyset,\emptyset\emptyset^{-}\otimes\overline{.\emptyset\emptyset}324.7\emptyset 323.\mathfrak{G}92\emptyset 621err69.2e-\emptyset 49.3e-\emptyset 78.3e-\otimes 78.5e-\mathfrak{d}7}{15\emptyset\emptyset err12.1e-123.9e-112.\emptyset e-163.4e-164.2e-17}$

err3
err6

$4.\emptyset e-12$$2.\emptyset e-12$

7. $5e-07$

3.8e-l5

3. $8e-15$

3.8e-l5$5.9e-\emptyset 7$$6.6e-\emptyset 7$$6.6e-\emptyset 7$$6.6e-\emptyset 7$$\frac{time[\sec]7\emptyset 4.51887.36325.28311.8315851}{1506Q1}$
err3$3.3e\prime 14$

3.3e-l4

3. $3e-14$$\frac{time[\sec]>2\emptyset,\emptyset \mathfrak{G}^{--}\emptyset\overline{.\emptyset\emptyset}>2\emptyset,\emptyset\emptyset\emptyset\overline{.\emptyset\emptyset}96\emptyset.21err65.9e-\emptyset 71.2e-\emptyset 65.4e-\emptyset 7}{15\emptyset 120err1----2.8e-162.6e-162.8e-16}$

978.65

err3

time

5

$[\sec]err6$$>2\emptyset,\mathfrak{G}\emptyset Q.QQ-$$593.9\emptyset$

2. $9e-14$

2. $9e-14$

2. $9e-14$$2,436.417.6e-\emptyset 7$$2,634.788.5e-\emptyset 7$$1,528.557.2e-07$$>2\emptyset,\emptyset\otimes^{-}\emptyset\overline{.\emptyset\emptyset}$

おわりに
本論文では 0-1 整数変数を含む非凸 2 次最適化問題 (MBNQO) に対する非負半正定値緩和問

題 (DNR) を解く際に現れる次の 2 つの問題点を解決する方法を提案した．1 つめの問題点は，(DNR) は実行可能内点解をもたないため，主双対内点法を用いて精度良く解くことが難しいことである．2 つめの問題点は，大規模な (DNR) を半正定値最適化問題に帰着させた場合，問題のサイズが非常に大きくなるため，既存のソフトウェアで解くことが難しいことである．
本研究では 1 つめの問題点を克服するために問題の退化を用いて (DNR) のサイズを縮小した．
この縮小は面的縮小 [26] の 1 反復とみなすことができるため，主双対内点法の安定性を高めることができる．面的縮小の適用後も等価な半正定値最適化問題のサイズが大きいという問題点は残る
上，面的縮小により問題の疎性が悪化する場合がある．本研究ではこの問題点を解決するために縮小後の問題に対する主双対パス追跡法を提案した．この手法は Toh [23] によって提案された凸 2 次
半正定値最適化問題に対する主双対パス追跡法に基づく．この手法では探索方向を計算するために大規模な線形方程式系を PSQMR 法 [11] を用いて解いた．また，PSQMR 法の収束性を高めるため，
Toh [23] によって提案された前処理を拡張し，3 つの前処理を提案した．数値実験の結果，提案手法

196
を用いることで非負半正定値緩和問題を速く正確に解くことができることがわかった．
以下に今後の課題を 2 っ挙げる．面的縮小に関しては，完全正値最適化問題に対して面的縮小を直接適用することにより，より強い下界を得ることができる可能性がある．そのためには，完全正値行列のなす錐の幾何学的構造 [8] などを利用する必要がある．主双対パス追跡法に関しては，係数行列の低ランク性を用いることが考えられる．(FRDNR) の係数行列図が低ランク構造をもっことから，SDPT3[24] で実装されているような計算方法を用いることにより，前処理の計算を高速化することができる．

参考文献[1] A. Berman and N. Shaked-Monderer, Completely Positive Matrices, World Scientific, Singapore,

2003.
[2] I. M. Bomze, On standard quadratic optimization problems, Joumal of Global optimization, 13(1998), 369-387.
[3] I. M. Bomze, Copositive optimization–recent developments and applications, Technical Repon(2011).

[4] I. M. Bomze, W. Schachinger and G. Uchida, Think co(mpletely )positive!

Matrix properties,

examples and a clustered bibliography on copositive optimization, Technical Report (2011).
[5] J. M. Borwein and H. Wolkowicz, Facial reduction for a cone-convex programming problem,

Journal of the Australian Mathematical Society, 30, 369-380, 1981.
[6] J. M. Borwein and H. Wolkowicz, Regularizing the abstract convex program, Joumal of Mathe-

matical Analysis and Applications, 83, 495-530, 1981.
[71 S. Burer, On the copositive representation of binary and continuous nonconvex quadratic pro-grams, Mathematical Programming, 120 (2009), 479-495.
[8] P. J. C. Dickinson, Geometry of the copositive and completely positive cones, Joumal of Mathematical Analysis and Applications, 380 (2011), 377-395.
[9] P. J. C. Dickinson and L. Gijben, On the computational complexity of membership problems forthe completely positive cone and its dual, Technical Report, (2011).

a survey. In M. Diehl, F. Glineur, E. Jarlebring and
Copositive programming
W. Michiels (eds.), Recent Advances in optimization and its Applications in Engineering (pp.
3-20), Springer, New York, 2010.
[11] R. W. Freund and N. M. Nachtigal, A new Krylov-subspace method for symmetric indefinite
linear systems, Proceedings of the 14th IMACS World Congress on Computational and Applied[10] M. D\"ur,–Mathematics, 1994, 1253-1256.
[12] D. Ge and Y. Ye, On doubly positive semidefinite programming relaxations, Technical Report(2010).

[13] K. Kobayashi, K. Nakata and M. Kojima, A conversion of an SDP having free variables into the197
standard form SDP, Computational optimization and Applications, 36 (2007), 289-307.
[14] M. Kojima, M. Shida and S. Sindoh, Search directions in the SDP and monotone SDLCP: gener-

alization and inexact computation, Mathematical Progmmming, 85 (1999), 51-80.
[15] A. N. Langville and W. J. Stewart, A Kronecker product approximate preconditioner for SANs,
Numerical Linear Algebra with Applications, 11 (2004), 723-752.
[16] H. D. Mittelmann, An independent benchmarking of SDP and SOCP solvers, Mathematical Pro-

gramming, 95 (2003), 407-430.
[17] G. Pataki, A simple derivation of a facial reduction algorithm and extended dual systems, Tech-nical Repon (2000).
[18] M. V. Ramana, An exact duality theory for semidefinite programming and its complexity impli-

cations, Mathematical Progmmming, 77, 129-162, 1997.
[19] M. V. Ramana, L. Tun\cael and H. Wolkowicz, Strong duality for semidefinite programming, SIAM

Joumal on optimization, 7, 3, 641-662, 1997.
[20] J. F. Sturm, Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones,

optimization Methods and Software, $11$ & $12$ (1999), 625-653.
[21] M. Tanaka, K. Nakata and H. Waki, Application of a facial reduction algorithm and an inexact

primal-dual interior-point method for doubly nonnegative relaxation for mixed binary nonconvex

quadratic optimization problems, Technical Report (2011).
[22] M. J. Todd, K. C. Toh and R. H. Tut\"untu, On the Nesterov-Todd direction in semidefinite pro-gramming, SIAM Joumal on optimization, 8 (1998), 769-796.
[23] K. C. Toh, An inexact primal-dual path following algorithm for convex quadratic SDP, Mathe-

matical Progmmming, 112 (2008), 221-254.
[24] K. C. Toh, M. J. Todd and R. H. T\"ut\"uncu, On the implementation and usage of SDPT3–a MAT-

LAB software package for semidefinite-quadratic-linear programming, version 4.0, Technical
Report (2010).
[25] K. C. Toh, R. H. Tutunc\"u and M. J. Todd, Inexact primal-dual path-following algorithms for a

special class of convex quadratic SDP and related problems, Pacific Joumal of optimization, 3(2007), 135-164.
[26] H. Waki and M. Muramatsu, Facial reduction algorithm for conic optimization problems, Tech-nical Report (2010).
[27] M. Yamashita, K. Fujisawa, K. Nakata, M. Nakata, M. Fukuda, K. Kobayashi and K. Goto,

A high-performance software package for semidefinite programs: SDPA 7, Technical Report(2010).