downhill simplex法による非線形最小二乗法について

概要

１．simplex法
　定量的なデータ解析において最小二乗法は欠くことのできない統計的手法である[1]が、その理論の解説は線形最小二乗法に基づくものが多く（例えば[2]）、現実に遭遇する非線形な系へ適用する際の障壁となっている。本稿では、非線形最小二乗法の中でも計算効率が優れて[3]実装が容易なdownhill simplex法[4]について、その計算原理と具体的な適用法について解説する。
　本来のdownhill simplex法（以下simplex法と呼ぶ）は、導関数を用いることなく多変数関数の最小値を求めるための手法で、しばしば線形計画法[5]のsimplex法と混同されるが、simplex（単体）を用いること以外に両者の共通点はない。n次元空間におけるsimplexとは、n+1個の頂点から構成される縮退していない凸な立体のことで、具体的には、2次元のsimplexは三角形、3次元のsimplexは四面体を指す。simplex法では、以下の手順[6]に従って、n次元の定義域でsimplexを変形しながら探索することにより、n変数関数の最小値を求めることができる（n=2の場合を図1に示す）。

1）関数の初期値を与える独立変数の値を座標として、頂点P0（位置ベクトル𝑝⃗0→）を1つ定める。
2）simplexの他のn個の頂点Pi（𝑖=1,⋯,𝑛）を、𝑝→𝜄=𝑝0→+𝜆𝑒→𝜄により定める（ただし𝑒1→,⋯, 𝑒𝑛→はn次元ベクトル空間の標準基底）。
3）simplexの各頂点Pi（𝑖=0,⋯,𝑛）において関数値を評価し、最も大きい値を与える頂点（最高点と呼ぶ）を求める。
4）simplexの最高点を、対面（n=2の場合は対辺）に関して等距離に反射し、反射点における関数値を評価する。
5）上記4)において、どの頂点よりも反射点の関数値の方が小さい場合は、さらに伸長する。
6）上記5)において伸長できない場合は、最高点を対面に近付けてsimplexを収縮する。
7）上記6)で収縮しても関数値が下がらない場合は、3)において最も小さい関数値を与える頂点（最低点と呼ぶ）に対して、残りの頂点を一様に近付ける（多重収縮）。
8）こうして得られた新たなsimplexに対して、上記3)以降を繰り返す。関数値の変化（またはsimplex内の関数値の分散）が閾値以下になったら、計算が収束したと判定する。

２．非線形最小二乗法への応用
　simplex法による非線形最小二乗法として、Gauss関数𝑦=𝑎𝑒−𝑏𝑥2に対するm個の実測データ(𝑥𝑖,𝑦𝑖)（ただし𝑖=1,⋯,𝑚）のフィッティングを取り上げる（図2左）。このときsimplex法で最小化する関数（目的関数と呼ぶ）は、独立変数がa,bの2変数関数として、
　𝑓(𝑎,𝑏)=𝑚∑𝑖=1(𝑦𝑖−𝑎𝑒−𝑏𝑥2)2/σ2𝑖
で与えられる[1]。ここで𝜎𝑖は𝑦𝑖の測定誤差を表す。
　なお𝑦=𝑎𝑒−𝑏𝑥2はlog𝑦=−𝑏𝑥2+log𝑎と変形できるから、横軸を𝑥2、縦軸をlog𝑦としてプロットする（Guinierプロットと呼ぶ）と、Gauss関数は直線で表される（図2右）。このためGuinierプロットは、フィッティング結果の診断や、a,bの初期値の推定によく用いられる。

３．アルゴリズムとしての実装
　まずsimplexの変形に関連した頂点の座標を、2次元配列として実装する。例えばn=2の場合、最初のsimplex（図1の網掛け）の3個の頂点は、simplex法の手順1)と2)より、(𝑎0,𝑏0),(𝑎0+𝜆,𝑏0),(𝑎0,𝑏0+𝜆)だから、それぞれ配列の第1列、第2列、第3列に置く。ここで𝑎0,𝑏0はGuinierプロット（図2右）等から推定したa,bの初期値で、𝜆は手順2)で定めたパラメータである。なおGauss関数のフィッティングではaとbに関する収束半径が異なるため、𝑒1→と𝑒2→で同一の𝜆を用いることが難しい。このため𝑎0,𝑏0の代わりに、実際にはlog𝑎0,log𝑏0を用いると、𝜆に対する影響を緩和することができる。
　次に、simplexの最高点の対辺（一般には対面）の重心（図1の点G）の座標を、上記配列の第4列に置く。ここで点Gの座標は、simplex全体の重心から最高点（図1の点H）の座標の寄与を引くことにより、簡単に求められる。この点Gに対して、点Hを反射した点（図1左上）の座標を配列の第5列に、反射したのちさらに伸長した点（図1左下）の座標を配列の第6列に置く。simplexの収縮（手順6）は、反射伸長（手順5）と排反なので、点Hを収縮した点（図1右上）の座標も配列の第6列に置くことができる。多重収縮（図1右下）では、元のsimplexの3頂点（配列の第1列、第2列、第3列）のうち最低点（図1の点L）以外の座標を、多重収縮後の座標に置き換えればよい。
　以上のように、simplex法の計算に必要な点の座標は、
　(𝑎0 𝑎0+𝜆　𝑎0　𝑎𝑔　𝑎′　𝑎′
　𝑏0　𝑏0　𝑏0+𝜆　𝑏𝑔　𝑏′　𝑏′′)
のように、2行6列（一般には𝑛行𝑛+4列）の2次元配列に全て保存することができる。ここで、(𝑎𝑔,𝑏𝑔)はsimplexの最高点の対辺（一般には対面）の重心の座標、(𝑎′,𝑏′)は最高点を反射した点の座標、(𝑎′′,𝑏′′)は最高点を反射伸長または収縮した点の座標である。なお第1列、第2列、第3列の成分は、simplexの変形を繰り返す過程（手順8）で逐次更新されていく点に注意されたい。さらに目的関数の値を評価する度に、最高点の座標を例えば第1列に、最低点の座標を第3列に保存しておくと、実際には都合がよい。
　なお上記で最高点Hを反射・反射伸長・収縮した点H’の座標は、分点の公式を用いて簡単かつ統一的に計算することができる。例えば図3において、反射のときは𝛼=1、反射伸長のときは𝛼>1、収縮のときは−1<𝛼<0とすればよい。多重収縮のときは、図3のHを最低点L以外のsimplexの各頂点に、GをLに置き換えたうえで、−1<𝛼<0を満たす値を適用すればよい。

４．誤差の評価
　まず2×𝑚型（一般には𝑛×𝑚型）のJacobi行列Jと、𝑚×𝑚型の重み行列𝑊を
　𝐽=(𝛛𝑦1/𝛛𝑎⋯𝛛𝑦𝑚/𝛛𝑎　𝛛𝑦1/𝛛𝑏⋯𝛛𝑦𝑚/𝛛𝑏),　𝑊=(1/𝜎2 1　𝑂　⋱　𝑂　1/𝛛2 𝑚)
として定義する。ここで𝑦𝑖=𝑎𝑒−𝑏𝑥2𝑖だから、
𝛛𝑦𝑖/𝛛𝑎=𝑒-𝑏𝑥2𝑖,　𝛛𝑦𝑖/𝛛𝑏=-𝑎𝑏𝑥2 𝑖𝑒-𝑏𝑥2 𝑖(-𝑥2 𝑖)
で与えられる。
　このとき、積𝑡𝐽𝑊𝐽（𝑡𝐽は𝐽の転置行列）は、2次（一般には𝑛次）の正方行列になり、最小二乗法によるa,bの推定値の誤差𝜎𝑎,𝜎𝑏は、
　σ2𝑎 =(𝑡𝐽𝑊𝐽)-1 11 𝑓(𝑎,𝑏)/𝑚-2,　σ2 𝑏 = (𝑡𝐽𝑊𝐽)-1 22 𝑓(𝑎,𝑏)/𝑚-2
として与えられる[1]。ここで、(𝑡𝐽𝑊𝐽)−1は𝑡𝐽𝑊𝐽の逆行列、(𝑡𝐽𝑊𝐽)−1は(𝑡𝐽𝑊𝐽)−1の(k,k)成分である（上式中の𝑚𝑚−2は一般には𝑚𝑚−𝑛𝑛）。
　与えられた正方行列の逆行列を求めるには、一般的な掃き出し法のアルゴリズム[7]を用いればよい。この際、掃き出しの要となる枢軸成分が0になる場合には、適当な行の交換が必要となる（例えば[8]）。
　以上の理論においては、𝑡𝐽𝑊𝐽が正則で、かつ(𝑡𝐽𝑊𝐽)−1の対角成分が全て非負であることが前提であるが、数値計算の実数の精度やm,nの多寡により、掃き出しの過程でrank(𝑡𝐽𝑊𝐽)<𝑛すなわち(𝑡𝐽𝑊𝐽)−1が存在しないこと（ランク落ち）が発生することがある。rank(𝑡𝐽𝑊𝐽)<𝑛では、n個のパラメータが独立でない可能性があるので、そもそもデータをGauss関数で一意的にフィッティングできるのか、またフィッティングするのに十分なデータ数が揃っているのか（nに対してmが十分に大きいか）再考する余地がある。

５．プログラムの作成と公開
　以上の背景を踏まえて、X線溶液散乱（SAXS）データから、分子の慣性半径と原点散乱強度（濃度と分子量に比例）を計算することができる[実例は9]。当研究室ではJavaScriptを用いて作成したWebアプリケーションを、ホームページで公開して国内外の研究者の利用に供している[10]。当アプリケーションは、2成分系のGauss関数𝑦=𝑎𝑒−𝑏𝑥2+𝑐𝑒−𝑑𝑥2（n=4の場合）にも対応している他、HTML5のCanvas APIを用いてプロットの度にブラウザ内で再描画するため、ユーザは条件を少しずつ変えながら、リアルタイムにフィッティング結果の正否を判断することができる（図4）。

この論文で使われている画像