10. 非線形最適化 (多変数)

10.1 概要

多変数の非線形関数 (一般関数) の最小点 (符号を逆にすれば最大点) を求める問題を考えます.
\[
min f(x_1, x_2, \dots, x_n)
\] \(x_1, x_2, \dots, x_n\) を要素とする列ベクトルを \(\boldsymbol{x}\) と表すことにすれば, 次のようにベクトル表現で書くことができます.
\[
min f(\boldsymbol{x})
\] 1 変数の場合と同様に, 対象とする領域内では 1 つの最小点しか持たないものとします. すなわち, 目的とする最小点の近傍だけに着目することにします. これを局所的最適化あるいは局所的最小点を求めるなどといいます.

10.2 多変数関数非線形最適化問題の解法

10.2.1 降下法

次のような反復を行って最小点を見つける方法を降下法といいます.
\[
\boldsymbol{x_{k+1}} = \boldsymbol{x_k} + \alpha_k\boldsymbol{d_k}
\] ここで, \(\boldsymbol{d_k}\) を方向ベクトル, \(\alpha_k\) をステップ長といいます.

この反復を図で表すと, まず現在の近似点から最小点があると思われる方向に向いた方向ベクトル \(\boldsymbol{d_k}\) を決め, 次に方向ベクトル上の最小点を求めて(1 変数の最適化を行って)ステップ長 \(\alpha_k\) を決めます. このようにして次々により最小点に近い近似点を求めていきます.

10.2.1.1 最急降下法

方向ベクトルとして最大勾配方向をとる方法を最急降下法といいます. 等高線が例えば円形に近ければすぐに最小点に近づいていきますが, 実際の問題では意外と収束は速くありません.

10.2.1.2 ニュートン法

関数 \(f\) が微分可能であるとき
\[
\nabla f(\boldsymbol{x}) =
\begin{pmatrix}
{\partial}f(\boldsymbol{x})/{\partial}x_1 \\
{\partial}f(\boldsymbol{x})/{\partial}x_2 \\
\vdots \\
{\partial}f(\boldsymbol{x})/{\partial}x_n \\
\end{pmatrix}
\] を \(f\) の勾配と呼びます.

さらに, \(f\) が 2 階微分可能であるとき
\[
\boldsymbol{H}(\boldsymbol{x}) =
\begin{pmatrix}
{\partial}^2f(\boldsymbol{x})/{{\partial}x_1{\partial}x_1} & {\partial}^2f(\boldsymbol{x})/{{\partial}x_2{\partial}x_1} & \cdots & {\partial}^2f(\boldsymbol{x})/{{\partial}x_n{\partial}x_1} \\
{\partial}^2f(\boldsymbol{x})/{{\partial}x_1{\partial}x_2} & {\partial}^2f(\boldsymbol{x})/{{\partial}x_2{\partial}x_2} & \cdots & {\partial}^2f(\boldsymbol{x})/{{\partial}x_n{\partial}x_2} \\
\vdots & \vdots & \ddots & \vdots \\
{\partial}^2f(\boldsymbol{x})/{{\partial}x_1{\partial}x_n} & {\partial}^2f(\boldsymbol{x})/{{\partial}x_2{\partial}x_n} & \cdots & {\partial}^2f(\boldsymbol{x})/{{\partial}x_n{\partial}x_n} \\
\end{pmatrix}
\] を \(f\) のヘッセ行列といいます. \(f\) の 1 階と 2 階の導関数が連続であれば \(\boldsymbol{H}(\boldsymbol{x})\) は正定値対称行列になります.

局所的最小点を求めるために方程式 \(\nabla f(\boldsymbol{x}) = 0\) を多変数非線形方程式の要領でニュートン法で解くことにすると, 反復式は次のようになります.
\[
\begin{align}
& \boldsymbol{d_k} = -\boldsymbol{H}(\boldsymbol{x_k})^{-1}\nabla f(\boldsymbol{x_k}) \\
& \boldsymbol{x_{k+1}} = \boldsymbol{x_k} + \boldsymbol{d_k} \\
\end{align}
\] \(\boldsymbol{d_k}\) は降下法の方向ベクトルになり, ニュートン方向と呼ばれることがあります. ステップ長 \(\alpha_k\) はニュートン法の式では 1 になりますが, 実際の問題ではニュートン方向上で 1 変数の最適化を行って定める必要があります.

ニュートン法は, ヘッセ行列を計算するのが面倒なこと, 反復の途中でヘッセ行列が正定値でなくなることがあるなどに注意が必要です.

10.2.1.3 準ニュートン法

ニュートン法の問題点を解決するためにヘッセ行列 \(\boldsymbol{H}(\boldsymbol{x})\) を適当な正定値行列 \(\boldsymbol{H}\) で近似する方法です.

代表的なものとして次の BFGS (Broyden, Fletcher, Goldfarb, Shanno) 公式があります.
\[
\boldsymbol{H_{k+1}} = \boldsymbol{H_k} – \boldsymbol{H_k}\boldsymbol{s_k}(\boldsymbol{H_k}\boldsymbol{s_k})^T/(\boldsymbol{s_k}^T\boldsymbol{H_k}\boldsymbol{s_k}) + \boldsymbol{y_k}\boldsymbol{y_k}^T/(\boldsymbol{s_k}\boldsymbol{y_k})
\] ここで, \(\boldsymbol{s_k} = \boldsymbol{x_{k+1}} – \boldsymbol{x_k}\), \(\boldsymbol{y_k} = \nabla f(\boldsymbol{x_{k+1}}) – \nabla f(\boldsymbol{x_k})\) を表します.

\(\boldsymbol{H_0}\) としては適当な正定値対称行列, 例えば単位行列をとり, 以下この漸化式を適用します.

10.2.2 信頼領域法

\(\boldsymbol{x_k}\) を中心とする半径 \(\delta_k\) の円を信頼領域と呼ぶことにします.

\(f\) を \(\boldsymbol{x_k}\) のまわりでテイラー展開し, 2 次までの項を採用して近似関数 \(\boldsymbol{q_k}(\boldsymbol{s})\) を作ります.
\[
\boldsymbol{q_k}(\boldsymbol{s}) = f(\boldsymbol{x_k}) + \boldsymbol{s}^T \nabla f(\boldsymbol{x_k}) + (1/2)\boldsymbol{s}^T\boldsymbol{H_k}\boldsymbol{s}
\] ここで, \(\boldsymbol{H_k}\) は \(\boldsymbol{x_k}\) における \(f\) のヘッセ行列または準ニュートン法のような近似行列とします.

半径 \(\delta_k\) の信頼領域はこの 2 次近似がおおむね妥当と思われる範囲を表すものです.

信頼領域法では, 領域内で近似関数 \(\boldsymbol{q_k}(\boldsymbol{s})\) の最小点 \(\boldsymbol{s_k}\) を求め, それを元の問題の次の点 \(\boldsymbol{x_{k+1}}\) として反復します. すなわち, 傾斜法における方向ベクトルとステップ長を一度に決めてしまいます.

\(\boldsymbol{H_k}\) が正定値対称で \(||\boldsymbol{H_k}^{-1}\nabla f(\boldsymbol{x_k})|| \leq \delta_k\) であれば, 信頼領域内で \(\boldsymbol{s_k} = \boldsymbol{H_k}^{-1}\nabla f(\boldsymbol{x_k})\) となります. そうでなければ \(\boldsymbol{s_k}\) は信頼領域外にある可能性があるので, 信頼領域の円周上で最小となる点を探して \(\boldsymbol{s_k}\) とします.

実際のプログラムでは収束状況に応じて信頼領域の半径 \(\delta_k\) を拡大・縮小する制御も行います.

10.3 XLPack を使った多変数関数非線形最適化問題の解き方

非線形最適化問題には, 準ニュートン法の VBA サブルーチン Optif0 を使うことができます. これは微分 (勾配とヘッセ行列) をそれぞれ有限差分近似とBFGS公式で計算しユーザーによる微分計算を不要にしたサブルーチンです. Optif0 は XLPack ソルバーから使うこともできます.

求められる解は局所的最小点であり, 一般に大域的最小点を求めることはできません. 局所的最小点はいくつも存在する可能性があり, 求めたい最小点に近い初期値を与えないと別の最小点に収束することがあります.

例題

次の関数の最小点を求める.
\[
f(x_1, x_2) = 100(x_2 – x_1^2)^2 + (1 – x_1)^2
\] 関数の形は次のようになります (縦軸は対数値で表しています).

この関数は, ローゼンブロックの関数と呼ばれ, (-1.2, 1) を出発点として, 最小点 (1, 1) を求めるテスト問題としてよく使用されます (なるべく近い初期値を与えるという原則には沿いませんがテスト用としてあえて遠くの初期値を与えています). 深い曲がった谷の中に最小点があり, 最小点に到達するのが難しいとされています.

10.3.1 VBA プログラムを使用した解き方 (1)

Optif0 を使ったプログラム例を示します.

Sub F(N As Long, X() As Double, Fval As Double)
    Fval = 100 * (X(1) - X(0) ^ 2) ^ 2 + (1 - X(0)) ^ 2
End Sub

Sub Start()
    Const NMax = 10, Ndata = 5
    Dim N As Long, X(NMax) As Double, Xpls(NMax) As Double, Fpls As Double
    Dim Info As Long, I As Long
    N = 2
    For I = 0 To Ndata - 1
        '--- Input data
        X(0) = Cells(6 + I, 1): X(1) = Cells(6 + I, 2)
        '--- Find min point of equation
        Call Optif0(N, X(), AddressOf F, Xpls(), Fpls, Info)
        '--- Output result
        Cells(6 + I, 3) = Xpls(0): Cells(6 + I, 4) = Xpls(1): Cells(6 + I, 5) = Info
    Next
End Sub

目的関数を外部サブルーチンとして用意し, 初期値を指定して Optif0 を呼び出します. Optif0 は, 初期値として与えられた近似解を出発点に反復計算を行い, その近傍の最小点を求めます. ここでは, 5 つの異なる初期値を与えたときにどのような解が得られるかを見るプログラムとしています.

初期値としては, (-1.2, 1) の他に (-1, 1), (-1, -1), (0, 1), (0, 0) を試してみます.

このプログラムを実行すると, 次の結果が得られました.

(-1, 1) のように初期値によっては収束しない場合がありました.

10.3.2 VBA プログラムを使用した解き方 (2)

リバースコミュニケーション版 (RCI) の Optif0_r を使ったプログラム例を示します.

Sub Start()
    Const NMax = 10, Ndata = 5
    Dim N As Long, X(NMax) As Double, Xpls(NMax) As Double, Fpls As Double
    Dim Info As Long, I As Long
    Dim XX(NMax - 1) As Double, YY As Double, IRev As Long
    N = 2
    For I = 0 To Ndata - 1
        '--- Input data
        X(0) = Cells(6 + I, 1): X(1) = Cells(6 + I, 2)
        '--- Find min point of equation
        IRev = 0
        Do
            Call Optif0_r(N, X(), Xpls(), Fpls, Info, XX(), YY, IRev)
            If IRev = 1 Then YY = 100 * (XX(1) - XX(0) ^ 2) ^ 2 + (1 - XX(0)) ^ 2
        Loop While IRev <> 0
        '--- Output result
        Cells(6 + I, 3) = Xpls(0): Cells(6 + I, 4) = Xpls(1): Cells(6 + I, 5) = Info
    Next
End Sub

目的関数を外部関数として与えるのではなく, IRev = 1 のときに XX() の値を使って関数値を計算し YY に入れて再度 Optif0_r を呼び出します. RCI の詳細については こちら を参照してください.

このプログラムを実行すると上と同じ結果が得られます.

10.3.3 ソルバーを使用した解き方

XLPack ソルバーアドインの「多変数非線形最適化」を使って解くこともできます. B9セルに数式 (=100*(C8-B8^2)^2+(1-B8)^2) が入力されています.

ソルバーについては こちら も参照ください.