光学薄膜設計ソフトの設計 その7「gslの最適化関数とのインタフェイス」

さっさと寝よ、とか言ってながらついつい。
実装の詳細に突っ込んでしまう前に、破綻が無いかいくつか確認。まず最適化をgslの関数でやることに問題ないかどうか。

まずちょっと整理。gslの多次元関数の最適化として35 Multidimensional Minimizationがある。
これはいくつかのアルゴリズムを同じインターフェイスで呼べるようになっていてマニュアルからコピーしてくるとver1.10では

1. Fletcher-Reeves conjugate gradient algorithm
2. Polak-Ribiere conjugate gradient algorithm
3. Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm
4. steepest descent algorithm
5. Simplex algorithm of Nelder and Mead

の五つがあるらしい。どう違うのかよくわからないけど、Simplex algorithm以外は最適化する関数の値と、その微分値が必要となる。光学薄膜の計算はそれほど素性が悪くないけど、微係数を求めるのは簡単ではない。なので、Simplex algorithmを使うことにする。

使うにあたってアルゴリズムのインスタンスを作ってやる必要がある。

gsl_multimin_fdfminimizer * gsl_multimin_fdfminimizer_alloc(
		const gsl_multimin_fdfminimizer_type * T, size_t n)

あるいは

gsl_multimin_fminimizer * gsl_multimin_fminimizer_alloc (
		const gsl_multimin_fminimizer_type * T, size_t n)

を呼んでやる。この二つの違いは最適化が微係数を必要とするときはひとつ目、微係数を使わないときはふたつ目を呼ぶ。Simplex algorithmの場合はふたつ目を使い、それ以外はひとつ目を使うことになる。

最適化するメリット関数（この関数の返す値をなるべく小さくする）は

double (* f) (const gsl_vector * x, void * params)

みたいな形でなければならない。つまり、最適化の変数をベクトルxで受け、そのほかのパラメータをparamsで受けて、値を返す。paramsの中身をgslは関知しないので計算に必要なパラメータを詰め込んでわたすとよい。最適化アルゴリズムによってはメリット関数の微係数を必要とするものがある（というかほとんどそう）ので微係数を計算する関数dfも定義しなければいけないけど、それは

void (* df) (const gsl_vector * x, void * params, gsl_vector * g)

の形でなければならない。xとparamsは同じ意味で、それぞれの変数で偏微分した値（grad）をベクトルgに入れて返す。メリット関数の形によっては微係数を計算するときにメリット関数の値も計算できてしまう場合があるので

void (* fdf) (const gsl_vector * x, void * params, double * f,gsl_vector * g)

と言う呼び方も用意されている。

このメリット関数をアルゴリズムのインスタンスに教えてやるために

int gsl_multimin_fdfminimizer_set (gsl_multimin_fdfminimizer * s,
	gsl_multimin_function_fdf * fdf, const gsl_vector * x, double step_size, double tol)

を呼ぶ。sは最適化のインスタンス、fdfはメリット関数をさす構造体、xは初期ベクトルである。 fdfは５個のメンバを持つ構造体で、メリット関数の値を返す関数や微係数へのポインタなんかを与えてやる。微係数を使わない場合、gsl_multimin_fminimizer_setを呼ぶ。これはgsl_multimin_function_fdfのかわりにgsl_multimin_functionと言う構造体を渡すようになっている。これは微係数無しのメンバがみっつの構造体。

実際の最適化は

int gsl_multimin_fdfminimizer_iterate (gsl_multimin_fdfminimizer * s)

を呼ぶ（微係数無しはgsl_multimin_fminimizer_iterate）ことで最適化のひとステップ分進めることができる。この関数の返す値は最適化の状況を表す。外でこの値を毎回チェックしながらループを回す。定数はgsl_errno.hに定義されている。帰る可能性のあるのを列挙すると

GSL_SUCCESS	収束した
GSL_CONTINUE	収束しなかったのでもっと回してくれ
GSL_ENOPROG	改善しないのでやめた方がいい
GSL_FAILURE	何らかのエラーが起きた
GSL_ETOL	十分小さくならなかった
GSL_ETOLG	微係数が十分小さくならなかった

などがありそう。全部に出会ったことはないけど。

普通はGSL_CONTINUEが帰ったらもう一度回して、GSL_SUCCESSが帰ったらループを止めて結果を収集する。それ以外だとエラーなり何なりでうまく行かなかったということになるので、諦めるなり、ステップサイズやトレランスを再調整してやり直すなりする。

さて、ここからが一番の問題。

普通、最適化ではいわゆるconstraintを導入する必要がある。constraintとは最適化で動かす変数に値の範囲の制約を与えること。例えば、光学薄膜の最適化の場合、少なくとも各層の膜厚は負になってはならないし、誘電体の屈折率は1より小さくなってはならない。また何十μmもあるあまりに厚い層や、高すぎる屈折率（誘電体で3を超えるなど）も解として望ましくない。また、場合によってはさらに狭い範囲に制限したいこともあり得る。例えば膜応力のバランスのために膜厚を50nm±10nmに入れたい、なんてことはよくある。

一般的には「制約付き最適化」などと呼ばれるが、普通の最適化ルーチンにはそういった機能は備わっていない。gslも同じ。

一般的にconstraintを導入する方法としてまず、「ラグランジュの未定乗数法」と言う方法がある。詳細はWikipediaを見てもらうとして制約条件を含めた形式に変形する方法。制約条件がg(x)=0の形に書けるとき、使える。今回のような「値がある範囲内に入る」というようなconstraintの場合、ちょっと工夫がいる。それと陽に微分がわかっている場合は適用が簡単だがそうでない場合、制約の無い最適化ルーチンに組み込むのは難しい。

もう一つの方法として、まず範囲の外に出てしまった変数はループに戻る前にその範囲の端の値に戻してループに入るようにする。そうすればその変数に関してだけそこで足踏みすることになる。これはステップを大きく変えないような最適化の場合、うまく行く。しかし普通の最適化アルゴリズムはループのたびに外で強制的に値を変更することを考えて作られていないので実装に依存する。また、gslは可能だけど、一般のライブラリではループの途中で変数ベクトルに手を加えることがそもそもできないものが多い。ということでこれが一番簡単だけど自分で最適化アルゴリズムを書くのでなければやめた方がいい。

また他には、変数が必要な値の範囲を出た場合、メリット関数の値として大きな数を返す、というのがある。これを「ペナルティ関数」と呼ぶらしい。 これも一見簡単だけど、スクエアな箱形（範囲を出ると単に大きな値を足す）ではメリット関数の素性が悪くなり、範囲の端にはりついたとき振動したりして収束性に悪影響を与える可能性がある。値の範囲の外でも連続で、できるなら微分可能な形が良い。つまり、バスタブ型が望ましいが、そのときどのようなカーブにすればいいかと言う問題が残る。範囲の外をあまり急にすると不連続にしたときと変わらなくなるし、また緩やかだと範囲の外でたまたまメリット関数が最小値となってしまう場合もある。その頃合いをどうするか。

なるべく緩やかなバスタブにしておいて結果として範囲の外に出てしまった場合、バスタブの角度を急にする、というのもある。しかしその場合、メリット関数が最適化の履歴に依存するということになる。普通の最適化のアルゴリズムはそこまで考えて作られていない。うまく行くかも知れないし、ダメかも知れない。最適化という作業そのものが強い非線形性を持っているのでその辺の見通しは悪い。

もし、最適化アルゴリズムがメリット関数の履歴依存を許すなら、このペナルティ関数を範囲の内側にも設定できる。外側でメリット関数は発散（と言うのは大袈裟だけど、極端に大きな値）して、内側でその壁に近づけばペナルティ関数の値が大きくなるようにしておくと最適化アルゴリズムは壁に近づく解を出さなくなる。ループを回しながらペナルティ関数の値をだんだん小さくして行けば、もし範囲の端でメリット関数が最小になる場合でも安定に収束する場合がある。

外側でバスタブになるのを「外点ペナルティ関数」、あとの内側に設定するのを「内点ペナルティ関数」などと呼ぶらしい。どちらにしても使いこなすのは難しい。この他にいろいろな制約を導入する方法が研究されている。

さいわい、光学薄膜の最適化はメリット関数の素性は素直（滑らかで微係数の大きな変化は少ない）なので上手いやり方が考えられる。メリット関数の構成をどうするかということにもつながるのでこの先は次回に。

gslの説明で終わってしまったな。