next up previous
Next: 4 制約つき最適化 Up: 計算天文学 II 第8回 最適化(1) Previous: 2 最適化手法の分類

Subsections


3 連続関数の最適化

ここで述べる手法は、とりえあえず定義域があまり変な形をしていなくて(と いうのをちゃんと定義することはできるが、やると話が長くなるので省略)、 関数はその定義域のなかで連続で極小値を1つだけ持つという場合のためのも のである。

この場合、目的関数$f$が2階微分を持てば、原理的には話は極めて簡単になっ て、$\nabla f=0$ という方程式をニュートン法で解けばいい。が、多くの場 合にこれはあんまりうまくいかない。というのは、変数の数 $n$ が多くなる と計算しないといけない2階微分の数は $n^2$ に比例して増えるわけで、計算 量が増えるだけでなく書かないといけないプログラムの量が増えるという問題 があるからである。

微分を数値的にやればいいのではとも考えられるが、これも丸め誤差等の影響 があって難しい場合が多い。

もっとも、微分を数値的にやるのではなく、数式的にやる、つまり、もとの関 数の数式から数式処理プログラムに生成させたり、あるいはプログラム自体を 「微分」する、つまり、目的関数を計算するプログラムから微分を計算するプ ログラムを自動生成するといった研究もかなり進んではいる。実際に問題を解 こうという時には、こういった手法が使えないかどうかも考えてはみるべきで あろう。

というわけで、以下はニュートン法とかではない方法。まず1変数、それから 多変数にいく。

3.1 1変数の場合

まあ、1変数ならニュートン法でいいのではないかとも思うわけだが、多変数 の場合のための準備ということで一応説明しておく。実際には、多変数の問題 を解く時に形式的には1変数の問題の繰り返しになって、そこでは簡単には2階 微分が計算できなくてニュートン法というわけにはいかないので、それ以外の 方法がいるわけである。

良く本に載っているのは、黄金分割法というものである。これは、以下のよう な方法である。区間$[a,b]$ のなかに関数 $f(x)$ の最小値があることはわかっ ているとする。

  1. $x_1$, $x_2$ を、$[a,b]$ をそれぞれ $\sqrt{5}-1:2$ およびその逆に 内分する点とする。
  2. それぞれの点で関数値 $f_1$, $f_2$ を計算する。
  3. $f_1>f_2$ なら最小値は $[x_1, b]$ にあるので、 $a$$x_1$ で置 き換え、1に戻る。細かいことをいえば $x_2$ が 次の $x_1$ になるので、使 い回せる。
  4. そうでなければ逆に $b$$x_2$ で置き換え、同様に 1に戻る。
  5. 上の全体を $\vert a-b\vert$ が十分小さくなるまで繰り返す。

なお、一般には別に黄金分割でなくても、区間内に適当に2点とってその大き い方を新しい区間の端にするというので構わない。黄金分割のミソは上の説明 の「細かいこと」、つまり、分割点の一方を使い回せるので反復一度について 関数計算が1度ですむということである。

一度ですむのはいいが、収束は遅い。一度反復した時に区間の幅が $0.62倍$ にしかならないからである。これはつまり一次収束で、反復毎に定数分の1に なるものである。

もうちょっと賢い方法としては、疑似ニュートン法的なものがある。要するに 3点あれば2次関数で近似できるので、その極値を求めようというものである。 最適化問題ではない非線形方程式では Secant 法と呼ばれているものの拡張で ある。 Secant 法程速くはないが、 一次よりも速い収束をすることが知られ ている。

3.2 多変数の場合

多変数の場合にも、黄金分割にあたるような直接探索法というのはあることは あるが、あんまり使えないので省く。大抵の本で最初に出ているのは最急降下 法というものである。とりあえず、関数$f$$N$次元ユークリッド空間全体で定 義されているとする。この方法の原理は、「ちょっと動いた時にもっとも関数 値が減る方向に動く」というのを繰り返すことである。もっとも減る方向は、


\begin{displaymath}
\Delta f = \nabla f\cdot \Delta \mbox{\boldmath$x$}
\end{displaymath} (1)

とテイラー展開の1次までとって、 $\vert\Delta \mbox{\boldmath$x$}\vert$が一定の時に$\vert\Delta f\vert$ を最大にすればよく、もちろん $\Delta x = \alpha \nabla f$、つまり一階 導関数自体が方向ということになる。 これで、あとは前節で述べた適当な方 法を使ってその方向での1次元最適化を適当にやり、あとはまたそこで新しく 勾配を求めて同じことを繰り返す。

最急降下法のよいところは、いつかは収束することであり、よくないところは 収束が必ずしも速くないことである。これは、2変数で目的関数が2次形式の場 合についていろいろ実験したり考えてみたりすればわかるが、結局直交する2 方向の繰り返しに陥ってしまうからである。この事情は多変数の場合でも同じ で、結局2方向になっていくということが証明されている。

というわけで、速く答を出そうとするなら、やはり疑似ニュートン法的なもの を考えたくなる。1変数の時のように簡単に2次形式の近似式が構成できるわけ ではないが、その辺をなんとかうまくやろうというわけでいろんな方法が提案 されている。

多分良くつかわれているのは DFP(Davidon-Fletcher-Powell)法や BFGS(Broyden-Fletcher-Goldfarb-Shanno)法で、どちらも考案者の名前である。 時々 Variable Metric method(可変計量法) と呼ばれることもある。

以下、基本的な考え方を説明する。

目的関数が、以下の2次形式

\begin{displaymath}
f(x) = \frac{1}{2}x^TQx
\end{displaymath} (2)

であるとする。 $x$$n$ 次元実ベクトル、 $Q$$n$次元正方行列で、 対称にとって良いので正定値である。

ニュートン法では、$f$ を2次形式で近似してそれの極値を求める。形式的に は、近似値 $x_0$ の回りでのテイラー展開

\begin{displaymath}
f(x-x_0) = f(x_0) +\nabla f(x_0)^T(x-x_0)+
\frac{1}{2}(x-x_0)^TH(x_0)(x-x_0)
\end{displaymath} (3)

の極値を求めるわけである。ここで $H$ はヘッシアンで、2階偏導関数を要素 とする正方行列である。つまり
\begin{displaymath}
h_{ij} = \frac{\partial^2f}{\partial x_i\partial x_j}
\end{displaymath} (4)

で、極値を与える点は、 $\Delta x = x-x_0$ として、
\begin{displaymath}
H(x_0)\Delta x = -\nabla f(x_0)
\end{displaymath} (5)

である。

目的関数が上の2次形式なら、これはもちろん $\Delta x = -Q^{-1}Qx_0 =
-x_0$で、正しい答が求まる。

そういうわけで、考え方としては、ヘッシアンなり $Q$ の近似値を作ってい こうというのが基本になる。

ここで、ちょっと定義を続ける。まず、共役という概念を定義する。2つのベ クトル $u, v$$Q$ について共役であるとは

\begin{displaymath}
u^TQv = 0
\end{displaymath} (6)

であることである。

さらに、互いに共役な$k$ 個のベクトル $p_0, ... p_{k-1}$を考え$(k<n)$、これら に対して、

\begin{displaymath}
H_k = \sum_{i=0}^{k-1}\frac{p_ip_i^T}{p_i^TQp_i}
\end{displaymath} (7)

を考えると、明らかに
\begin{displaymath}
H_kQp_i = p_i\quad(i=0,1, ..., k-1)
\end{displaymath} (8)

である。

これは、 $H_k$ $p_0, ..., p_{k-1}$ が張る部分空間では $Q^{-1}$ みた いなものであるということを意味している。したがって、$p_i$ を順番に発生 させる方法がなにかあれば、それを使って、

\begin{displaymath}
H_{k+1} = H_k + \frac{p_kp_k^T}{p_k^TQp_k}
\end{displaymath} (9)

$H_k$ を計算していけばよい。

ここでの実際的な問題は、互いに共役な $p_k$ を作るよい方法があるかどうかよく わからないことである。一つの考え方は、反復による修正量 $s_k =
x_{k+1}-x_k$ を使うことである。 $s_k$ は共役ではないが、それでもなんと か

\begin{displaymath}
H_{k+1}Qs_k = s_k
\end{displaymath} (10)

が成り立つように、式(9)を適当に変更する。 変更のために右辺に $C_k$ という項を追加することにすると、 これの満たす べき式は
\begin{displaymath}
(H_k+C_k)y_k = 0
\end{displaymath} (11)

ここで $y_k$
\begin{displaymath}
y_k = Qs_k = \nabla f(x_{k+1}) - \nabla f(x_{k})
\end{displaymath} (12)

である。

このように $C_k$ をとる一つの方法(DFP法)は、

\begin{displaymath}
C_k = - \frac{H_ky_ky_k^TH_k}{y_k^TH_ky_k}
\end{displaymath} (13)

とするものである。

まあ、自分でプログラムを書くならこれをつかうので OK であろう。ライブラ リとかだと BFGS のほうが少しよいようである。

3.3 CG法

さて、さっき共役な $p_k$ を作る方法はよくわからないと書いたが、これは 原理的には知られている。但し、いつでも上手くいくわけではないのでそれを 使わない方法も研究されているわけである。

共役な $p_k$ を直接作る方法が共役勾配法、すなわち CG(Conjugate gradient)法である。これは形式的には簡単である。いま、 $g_k = \nabla f(x_k)$ と書くことにして、

$\displaystyle p_0$ $\textstyle =$ $\displaystyle -g_0$ (14)
$\displaystyle p_k$ $\textstyle =$ $\displaystyle -g_k + \beta_kp_{k-1}\quad (k=1,2, ...)$ (15)

というベクトル列を考える。ここで $\beta_k$
\begin{displaymath}
\beta_k = \frac{g_k^TQp_{k-1}}{p_{k-1}^TQp_{k-1}}
\end{displaymath} (16)

である。

$f$ が2次関数の時には、これで求まる $p_k$ は互いに共役であることを証明 できる。さらに、ちょっと変形すると、

\begin{displaymath}
\beta_k = \frac{\vert g_k\vert^2}{\vert g_{k-1}\vert^2}
\end{displaymath} (17)

となって、これは簡単に計算できる。

なお、疑似ニュートン法、CG法のどちらの場合でも、方向は決まるが1次元問 題を繰り返し毎に解く必要はある。これは黄金分割なり疑似ニュートン法なり を使うことになる。多次元の反復では勾配を求めないといけないので少なくと も$n$個の関数を計算することになるが、1次元問題では反復毎に1度 $f$ を計 算するだけなのでここではそれほど効率に気を使う必要はないことに注意。

3.4 CG 法の応用:連立1次方程式

CG 法は現在最適化よりも大規模な線型方程式を解くのに広く使われている。 今、

\begin{displaymath}
Ax = b
\end{displaymath} (18)

という線型連立方程式を考える。で、A は対称行列であるとする。この時、上 の方程式は、以下の2次形式
\begin{displaymath}
f(x) = \frac{1}{2}x^tAx - bx
\end{displaymath} (19)

を最小化するものと考えることができる。ここで CG 法を使うというのが基本 的な考え方である。CG法なので1次元方向の最小化がいるが、これは線型問題 なので答がわかっている。

何故こんなややこしいことをするのかと思うかもしれない。理由はちゃんとあっ て、一つはこの方法は反復法であるにもかかわらず原理的には有限回で収束す ること、つまり、ガウスザイデルや SOR とは違って、(計算精度が無限に高い なら)収束が保証されていることである。もう一つはいろいろ工夫することで 収束を非常に速くできることが多いということである。

収束を速くするための工夫は色々な「前処理」と言われるもので、それだけを扱った本 がいっぱいあるのでここではこれ以上は扱わない。


next up previous
Next: 4 制約つき最適化 Up: 計算天文学 II 第8回 最適化(1) Previous: 2 最適化手法の分類
Jun Makino
平成16年12月6日