非線形最適化問題に対する主双対内点法の収束性について

矢部 博（東京理科大学理学部応用数学科）

1. はじめに

本講演では，非線形最適化問題

minimize f(x), x ∈ R^{n},
subject to g(x) = 0, x_{i} ≥ 0, i ∈ I_{P} ⊂ {1, 2, ..., n}

を解くための主双対内点法を考える．ただし， f : R^{n} → R ， g : R^{n} → R^{m}， I_{P} は添字集合で p = |I_{P} | とする．ここで，行ベクトルが e_{i}^{t}, i ∈ I_{P} から成る p × n 行列を Eとすれば，非負条件は x′ ≡ Ex ≥ 0 と書ける． Lagrange 関数を

L(w) = f(x) − y^{t}g(x) − z^{t}Ex

と定義したとき，上記の最適化問題に対する Karush-Kuhn-Tucker(KKT) 条件は

r_{0}(w) = ( \begin{matrix}∇_{x}L(w) \\ g(x) \\ X′Ze \end{matrix} ) = 0, x′ ≥ 0, z ≥ 0, (1)

で表される．ただし， w = (x, y, z)^{t} とし, y ∈ R^{m} と z ∈ R^{p} はそれぞれ等式制約条件と非負条件に対応する Lagrange 乗数であり，また， A(x) をその行ベクトルが ∇g_{i}(x)^{t}からなる m×n 行列， X′ = diag (x′_{1}, ⋯, x′_{p}) , Z = diag (z_{1}, ⋯, z_{p}) , e = (1, ⋯, 1)^{t} ∈ R^{p}とする．

主双対内点法の基本的な考えは，バリヤパラメータ μ を導入して相補条件 (XZe = 0) を XZe = μe に置き換えて得られるバリヤ KKT 条件にニュートン法を適用するものである [10] ．その際，直線探索法の枠組みで大域的収束性を実現するために， Yamashita[12] が提案した l_{1} 型バリヤペナルティ関数が主変数 x に対するメリット関数としてよく使われているが， (y, z) 空間では別々にステップ幅を選ぶ必要がある．そこで本講演では，主双対空間 ((x, z) 空間) におけるメリット関数を構築することによって，ステップ幅の選び方を統一することを考える．提案するメリット関数は l_{1}型バリヤペナルティ関数を利用したものである．そして提案する主双対内点法の大域的収束性を示す．最後に，局所的な収束速度についても触れる

なお，直線探索法を用いた場合の大域的収束性については [1],[2],[12] 等を，信頼領域法を用いた場合の大域的収束性については [3],[4],[6],[14] 等を参照されたい． l_{2}型ペナルティ関数に基づいたメリット関数を用いた場合のアプローチが [8],[15] 等で研究されている．さらに，局所的な収束性（収束速度など）に関する研究は [5],[7],[9],[11], [13] 等を参照されたい．

以下では，添字 k はアルゴリズムの外部反復や内部反復の反復回数を表し， ∥ · ∥ は l_{2} ベクトルノルムおよびそれから導かれる作用素ノルムを表す．また S = {x ∈ R^{n}| x′ > 0}, R_{+}^{p} = {z ∈ R^{p}| z > 0} とおく.

2. バリヤ KKT 条件と主双対メリット関数

F_{B}(x; μ) = f(x) − μ ∑_{i=1}^{p} log x′_{i}, x ∈ S

で定義されるバリヤ関数を考える．ただし， μ > 0 はバリヤパラメータである．このとき，もとの最適化問題の代わりに次の問題

minimize F_{B}(x; μ), x ∈ S
subject to g(x) = 0

(2)

を考えると，この問題に対する最適性の１次の必要条件は

∇f(x) − A(x)^{T} y − μE^{t}(X′)^{−1}e = 0, g(x) = 0, (x′ > 0)

で与えられる．ここで z = μ(X′)^{−1}e ∈ R^{p} とおけば， ∇f(x) − A(x)^{T} y − z = ∇_{x}L(w), XZe = μe となるので，上記の最適性条件は

r(w; μ) ≡ ( \begin{matrix}∇_{x}L(w) \\ g(x) \\ XZe − μe \end{matrix} ) = ( \begin{matrix}0 \\ 0 \\ 0 \end{matrix} ) , (x > 0, z > 0) (3)

で表される．これをバリヤ KKT 条件という．これはちょうど KKT 条件 (1) で，相補条件 XZe = 0 を XZe = μe で置き換えたものである．

他方，問題 (2) に対して次の l_{1}型バリヤペナルティ関数

F_{0}(x; μ) = f(x) − μ ∑_{i=1}^{p} log x′_{i} + ρ ∑_{i=1}^{m} |g_{i}(x)|

が提案されている [12] ．ただし， ρ > 0 はペナルティパラメータである．この l_{1}型バリヤペナルティ関数を無制約最小化する問題の最適性の１次の必要条件は，ペナルティパラメータρが十分に大きい場合にはバリヤ KKT 条件 (3) で表されることが知られている．

バリヤ KKT 条件にニュートン法を適用すれば，ニュートン方向Δw = (Δx, Δy, Δz)^{t}は次の連立１次方程式（ニュートン方程式）の解として定義される．

J(w)Δw = −r(w, μ), J(w) = ( \begin{matrix}G & −A(x)^{t} & −E^{t} \\ A(x) & 0 & 0 \\ ZE & 0 & X′ \end{matrix} ) ,

ただし，行列 G はヘッセ行列 ∇_{x}^{2}L(w) もしくはその準ニュートン近似である．バリヤペナルティ関数 F_{0}(x + s; μ) の１次近似を

F_{l}(x; s) ≡ f(x) + ∇f(x)^{T} s − μ ∑_{i=1}^{p} (log x′_{i} + \frac{(Es)_{i}}{x′_{i}} ) + ρ ∑_{i=1}^{m} |g_{i}(x) + ∇g_{i}(x)^{T} s|

と定義したとき，任意方向 s に対する関数 F_{0}とその１次近似との差ΔF_{l}は

ΔF_{l}(x; s) ≡ F_{l}(x; s) − F_{0}(x; μ)
= ∇f(x)^{T} s − μe^{T} (X′)^{−1}Es + ρ ∑_{i=1}^{m} |g_{i}(x) + ∇g_{i}(x)^{T} s| − ρ ∑_{i=1}^{m} |g_{i}(x)|

で与えられる．これは方向微係数に相当するものである．

次に直線探索法で用いるメリット関数について考える．我々は，次の l_{1} 型主双対メリット関数を提案する．ただし ρ′ > 0 である．

F (x, z; μ) = F_{0}(x; μ) + ρ′ϕ(x, z), ϕ(x, z) = log \frac{(x′)^{t}z/p + ∥X′z − μe∥^{2}}{( ∏_{i=1}^{p} x′_{i}z_{i})^{1/p}}

3. アルゴリズムと大域的収束性バリヤ KKT 条件を用いた主双対内点法のアルゴリズムは以下の通りである．

[アルゴリズム IP] （外部反復） Step 0. (初期設定) パラメータ ε > 0, M_{c} > 0 および，正の数列 {μ_{k}} , μ_{k} ↓ 0 を与える． k = 0 とおく． Step 1. もし ∥r_{0}(w_{k})∥ ≤ ε ならば停止する． Step 2. 近似バリヤ KKT 条件 ∥r(w_{k+1}, μ_{k})∥ ≤ M_{c}μ_{k} を満たす点 w_{k+1}を求める． Step 3. k := k + 1 とおいて Step 1 へ行く. □

このアルゴリズムの Step 2 (近似バリヤ KKT 点の計算) の部分は，次のアルゴリズムとして記述出来る．

[アルゴリズム LS] （内部反復）バリヤパラメータμ > 0 と正定数 M_{c}が与えられる． Step 0. (初期設定) w_{0} ∈ S × R^{m} × R_{+}^{p}, ρ > 0, γ ∈ (0, 1), β ∈ (0, 1), ε_{0} ∈ (0, 1) を与える. k = 0 とおく. Step 1. もし ∥r(w_{k}, μ)∥ ≤ M_{c}μ ならばアルゴリズム IP(外部反復) に戻る． Step 2. ニュートン方程式 J(w_{k})Δw_{k} = −r(w_{k}, μ) の解Δw_{k}を探索方向に選ぶ．ただし

J(w_{k}) = ( \begin{matrix}G_{k} & −A(x_{k})^{t} & −E^{t} \\ A(x_{k}) & 0 & 0 \\ Z_{k}E & 0 & X′_{k} \end{matrix} ) ,

G_{k}はヘッセ行列 ∇_{x}^{2}L(w_{k}) もしくはその近似行列である． Step 3. ᾱ_{k} = min {γα_{kmax}, 1},

α_{kmax} = min \{min_{i} \{−\frac{(x′_{k})_{i}}{(Δx′_{k})_{i}} | (Δx′_{k})_{i} < 0\} , min_{i} \{−\frac{(z_{k})_{i}}{(Δz_{k})_{i}} | (Δz_{k})_{i} < 0\}\}

として，次のアルミホ基準

F (x_{k} + ᾱ_{k}β^{l_{k}}Δx_{k}, z_{k} + ᾱ_{k}β^{l_{k}}Δz_{k}; μ) − F (x_{k}, z_{k}; μ) ≤ ε_{0}ᾱ_{k}β^{l_{k}}ΔF (x_{k}, z_{k}; Δx_{k}, Δz_{k})

を満足する最小の非負整数 l_{k}を求めてα_{k} = ᾱ_{k}β^{l_{k}}とおく．ただし，

ΔF (x, z; Δx, Δz) = ΔF_{l}(x; Δx) + ρ′∇ϕ(x, z)^{t} ( \begin{matrix}Δx \\ Δz \end{matrix} ) .

Step 4. w_{k+1} = w_{k} + α_{k}Δw_{k}, k := k + 1 とおいて Step 1 へ行く. □

大域的収束性示すために次の条件を仮定する． (G1) f と g は２回連続的微分可能である. (G2) 初期点における F (x, z; μ) の準位集合 {(x, z) ∈ S× R_{+}^{p}|F (x, z; μ) ≤ F (x_{0}, z_{0}; μ)} はコンパクトである. (G3) (G2) で定義された準位集合において rank A(x) = m が成り立つ． (G4) 任意の kに対して ρ > ∥y_{k} + Δy_{k}∥_{∞}が成り立つ． (G5) 行列 G_{k}は一様に正定値かつ有界である． □

このとき，次の定理によって探索方向が F (x, z; μ) の降下方向になることが評価される．

[figure]

[定理 1] 探索方向 (Δx_{k}, Δz_{k}) は次式を満たす．

ΔF (x_{k}, z_{k}; Δx_{k}, Δz_{k}) ≤ −Δx_{k}^{t}(G_{k} + E^{t}(X′_{k})^{−1}Z_{k}E)Δx_{k} − (ρ − ∥y_{k} + Δy_{k}∥_{∞})∥g(x_{k})∥_{1}

−ρ′ \frac{∥X′_{k}Z_{k}e − μe∥^{2}}{(x′_{k})^{t}z_{k}/p + ∥X′_{k}z_{k} − μe∥^{2}}

[figure]

この定理より，固定された μ > 0 に対するアルゴリズム LS の大域的収束性が示せる．

[figure]

{[定理 2] アルゴリズム LS によって生成される点列 \{w_{k}\} には少なくとも一つ集積点が存在する．さら \atop に，任意の集積点はバリヤ KKT 条件を満足する．}

[figure]

したがってμ_{k} → 0 であることに注意すれば，定理 2 よりアルゴリズム IP の大域的収束性が証明出来る．

4. 局所的超１次収束性局所的収束性を調べるために，次のアルゴリズムを考える．（以下では簡単のために p = n とする）

[アルゴリズム IPlocal] Step 0. 初期点 w_{0} ∈ R_{+}^{n} × R^{m} × R_{+}^{n}とパラメータμ_{0} > 0, 0 < γ_{0} < 1, ε > 0 を与える． k = 0 とおく． Step 1. もし ∥r_{0}(w_{k})∥ ≤ εならば停止する． Step 2. パラメータμ_{k} > 0, 0 < γ_{k} < 1 を与える． Step 3. ニュートン方程式

J(w_{k})Δw_{k} = −r(w_{k}; μ_{k})

を解いて探索方向Δw_{k}を求める． Step 4. γ_{k}に対してステップ幅

α_{k} = min \{ 1, γ_{k} min_{i} \{−\frac{(x_{k})_{i}}{(Δx_{k})_{i}} | (Δx_{k})_{i} < 0\} ,

γ_{k} min_{i} \{−\frac{(z_{k})_{i}}{(Δz_{k})_{i}} | (Δz_{k})_{i} < 0\}\}

を計算する． Step 5. w_{k+1} = w_{k} + α_{k}Δw_{k} とおく． Step 6. k := k + 1 とおいてステップ 1 へ行く．

ここでは，中心パスの近傍をたどりながら KKT 点に近づいて行った場合の点列の挙動について考える．以下では，反復回数 kは十分に大きいとし，パラメータμ_{k}の値は十分に 0 に近いとし，さらに次の条件を仮定する．

仮定 CP

(CP1) 点列 {w_{k}} は w^{∗}に収束する．

(CP2) 関数 f, g_{i}, i = 1, ..., m, のヘッセ行列が x^{∗}においてリプシッツ連続である． (CP3) w^{∗}において１次独立制約想定，２次の十分条件，狭義相補条件が成り立つ．

(CP4) min(1, τ_{2}) > τ_{1}, 0 < κ < 1 となる正定数τ_{1}, τ_{2}, κに対して，パラメータμ_{k}とγ_{k} が次式を満足する．

μ_{k} = ξ_{k}∥r_{0}(w_{k})∥^{1+τ_{1}}, 1 − γ_{k} = κξ_{k}∥r_{0}(w_{k})∥^{τ_{2}}
ただし， ξ_{k}は定数 M′に対して\frac{1}{M′} ≤ ξ_{k} ≤ M′を満たす正数である．

(CP5) 0 < M_{c} < \sqrt{n}.

(CP6) すべての kに対して，行列 G_{k}が

∥G_{k} − ∇_{x}^{2}L(w^{∗})∥ < δ,
∥(G_{k} − ∇_{x}^{2}L(w_{k}))Δx_{k}∥ = O(∥Δw_{k}∥^{1+τ_{3}})

を満足する．ただし， δ > 0 は十分に小さな数であり， τ_{3}はτ_{3} > τ_{1}を満たす正定数である．

以上の仮定のもとで次の定理を得る．

[figure]

[定理 3]([14]) 仮定 CP が成り立つとする．このとき，アルゴリズム IPlocal で生成される点列 {w_{k}} に対して次のことが成り立つ．

(1) もし点ŵ ∈ R_{+}^{n} × R^{m} × R_{+}^{n} が ∥r(ŵ; μ_{k})∥ ≤ M_{c}μ_{k}を満たすならば，適当な正定数ν_{1}, ν_{2}がとれて

ν_{1}∥r_{0}(w_{k})∥^{1+τ_{1}} ≤ ∥r_{0}(ŵ)∥ ≤ ν_{2}∥r_{0}(w_{k})∥^{1+τ_{1}}

が成り立つ．

(2) もし ∥r(w_{k}; μ_{k−1})∥ ≤ M_{c}μ_{k−1}ならば，ステップ幅がα_{k} = 1 となり，次式が成り立つ．

∥r(w_{k} + Δw_{k}; μ_{k})∥ ≤ M_{c}μ_{k} [figure]

この定理によれば，もし w_{k}がμ_{k−1}に対して近似バリヤ KKT 条件を満たすならば，アルゴリズム IPlocal

のステップ 4 でα_{k} = 1 が採用されて，次の点 w_{k} + Δw_{k}が新しいパラメータμ_{k}に対する近似バリヤ KKT

条件を満足する．したがって，

ν_{1}∥r_{0}(w_{k})∥^{1+τ_{1}} ≤ ∥r_{0}(w_{k+1})∥ ≤ ν_{2}∥r_{0}(w_{k})∥^{1+τ_{1}}

が成り立つので，点列 {w_{k}} の超１次収束性が保証される．

なお， (x, z) 空間における主双対内点法の局所的超１次収束性についての議論は Yabe and Yamashita[11]

を参照されたい．

参考文献

[1] I.Akrotirianakis and B.Rustem, A globally convergent interior point algorithm for general nonlinear programming problems, Technical Report 97-14, Department of Computing, Imperial College of Science, Technology and Medicine, March 1999.

[2] M.Argaez and R.A.Tapia, On the global convergence of a modified augmented Lagrangian linesearch interior point Newton method for nonlinear programming, Technical Report TR95-38, Department of Computational and Applied Mathematics, Rice University, October 1995 (revised February 1997).

[3] R.H.Byrd, J.C.Gilbert and J.Nocedal, A trust region method based on interior point techniques for nonlinear programming, Technical Report OTC 96/02, Optimization Technology Center, Argonne National Laboratory, June, 1996.

[4] R.H.Byrd, M.E.Hribar and J.Nocedal, An interior point algorithm for large scale nonlinear programming, Technical Report OTC 97/05, Optimization Technology Center, Argonne National Laboratory, August, 1997.

[5] R.H.Byrd, G.Liu and J.Nocedal, On the local behaviour of an interior point method for nonlinear programming, in Numerical analysis 1997, D.F.Griffiths, D.J.Higham and G.A Watson eds., Longman (1998), pp.37-56.

[6] J.E.Dennis, Jr., M.Heinkenschloss and L.N.Vicente, Trust-region interior-point SQP algorithms for a class of nonlinear programming problems, TR94-45, Dept. of Computational and Applied Mathematics, Rice University, Houston, Texas, USA, December 1994 (revised November 1995).

[7] A.S.El-Bakry, R.A.Tapia, T.Tsuchiya and Y.Zhang, On the formulation and theory of the Newton interior-point method for nonlinear programming, Journal of Optimization Theory and Applications, 89 (1996), pp.507-541.

[8] A.Forsgren and P.E.Gill, Primal-dual interior methods for nonconvex nonlinear programming, SIAM J. on Optimization, 8 (1998) pp.1132-1152.

[9] H.J.Martinez, Z.Parada and R.A.Tapia, On the characterization of Q-superlinear convergence of quasi-Newton interior-point methods for nonlinear programming, Bol. Soc. Mat. Mexicana, 1 (1995), pp.137-148.

[10] 矢部博，八巻直一「非線形計画法」，朝倉書店, 1999.

[11] H.Yabe and H.Yamashita, Q-superlinear convergence of primal-dual interior point quasi-Newton methods for constrained optimization, Journal of the Operations Research Society of Japan, 40 (1997), pp.415-436.

[12] H.Yamashita, A globally convergent primal-dual interior point method for constrained optimization, Optimization Methods and Software, 10 (1998), pp.443-469.

[13] H.Yamashita and H.Yabe, Superlinear and quadratic convergence of some primal-dual interior point methods for constrained optimization, Mathematical Programming, 75 (1996), pp.377-397.

[14] H.Yamashita, H.Yabe and T.Tanabe, A globally and superlinearly convergent primal-dual interior point trust region method for large scale constrained optimization, Technical Report, Mathematical Systems, Inc., Tokyo, Japan, July 1997 (revised July 2000).

[15] H.Yamashita and H.Yabe, An interior point method with a primal-dual l_{2} barrier penalty function for nonlinear optimization Technical Report, Mathematical Systems, Inc., Tokyo, Japan, March 1999.