Apache commons math の Levenberg-Marquardt法で、パラメータに制約条件を付ける

１か月ぶりの更新。
かなりマニアックなネタ。

Apache commons mathは、数値計算用のJavaのライブラリで、Levenberg-Marqurdt法によるフィッティングや、機械学習など、いろんなメソッドが用意されている。Fijiにも外部ライブラリとして利用されている。
２次元のガウス分布でのフィッティング（例えば蛍光輝点のフィッティングなど）には、
optimizationライブラリの、LevenbergarquardtOptimizerを使う。
詳しいやり方は省くが、

1. モデル関数とそのモデル関数のjacobianを、それぞれMultivariateVectorFuncion, MultivariateMatrixFunctionとして作る
2. LevenbergMarquardtOptimizerのオブジェクトを作る
   

   LevenbergMarquardtOptimizer lmo = new LevenbergMarquardtOptimizer();

3. LeastSquaresProblemを作り、lmo.optimizeのパラメータとして渡す。返り値として、最適化されたパラメータを受け取る
   

   LeastSquaresOptimizer.Optimun lsoo = lmo.optimize(lsp); // lsp -> LeastSquaresProblemのオブジェクト

ここで、lspはLeastSquaresBuilderで作ることができる。

LeastSquaresBuilder lsb = new LeastSquaresBuilder();
lsb.model(modelfunction(), jacobianfunction()); // モデル関数とそのjacobian
lsb.target(data); // 入力データ
lsb.start(newStart); // パラメータの初期値
lsb.maxIterations(1000); // 最大itterationの回数

例えば有名なグラフ・統計解析ソフトのOriginでは、非線形フィッティングの際に、パラメータに制約条件を付けることができる。Apache commons mathでもできるが、詳しいやり方はちゃんと書いてない。一応、LeastSquaresのParameters Validationにそれっぽい記述がある。

In some cases, the model function requires parameters to lie within a specific domain. For example a parameter may be used in a square root and needs to be positive, or another parameter represents the sine of an angle and should be within -1 and +1, or several parameters may need to remain in the unit circle and the sum of their squares must be smaller than 1. The least square solvers available in Apache Commons Math currently don't allow to set up constraints on the parameters. This is a known missing feature. There are two ways to circumvent this.

Both ways are achieved by setting up a ParameterValidator instance. The input of the value and jacobian model functions will always be the output of the parameter validator if one exists.

どうやらこの、ParameterValidatorを使えばできそうな雰囲気だが、Javadocを見てもほとんど情報がない。ParameterValidatorはinterfaceで、validate(RealVector params)というメソッドを実装するインターフェースである。なんとなくこのparamsというのが怪しそうだが、RealVector型もよくわからん。

LeastSquaresBuilderで作ったlsbのメソッドで、parameterValidatorがあるのでなんとなく使ってみる。

lsb.parameterValidator(new ParameterValidator() { ... }

ここで、parametervValidatorの引数はParameterValidatorを実装するクラスだが、いちいちクラスを作るのは面倒臭いので内部クラスとして作る。

lsb.parameterValidator(new ParameterValidator() {
    @Override
    public RealVector validate(RealVector params) { ... }
});

Eclipseを使えば、インターフェースで実装するべきメソッドを自動で作ってくれる。ここで、先ほどの記述に戻ると、

The input of the value and jacobian model functions will always be the output of the parameter validator if one exists.

なんとなく訳すと、「値とjacobian モデル関数の入力は、Parameter Validatorが存在するときには常にその出力となります」というようなことが書いてあり、多分iterationのたんびにparameterValidatorが呼び出されることになりそう。あとは、パラメータにアクセスできればよい。 RealVector型であるparamsからパラメータにアクセスするには、

getEntry(int index)
setEntry(int index, double value)

メソッドを使えばよい。indexはパラメータを示すインデックスである。 例えば、

public RealVector validate(RealVector params) {
 System.out.println("parameter[0]: " + params.getEntry(0));
 return params;
}

としてフィッティングを行うと（ここではとあるモデル関数によるフィッティングを行っている）、

v0: 255.0

v0: 76.64496657297977

v0: 215.79832218070536

...

v0: 258.1698480422253

というように、iterationごとでv0がどのような値になったかを出力してくれるようになる。

従って、例えばv0は必ず0より大きな値をとる必要があるのであれば、

if (params.getEntry(0) < 0) params.setEntry(0, 0.00001);

とすれば、もしiterationの途中で0以下となったとしても、0.00001に強制的に引き戻すことができる（ここで0にすると、パラメータが０に張り付いてしまうことがあるので、微妙に0より大きな値にしている）。

これで、フィッティングのパラメータに制約を与えることができるようになった。

ちなみにこの方法をネット上で検索したが、英語のサイトでもまず見つからない。かなりレアな情報だと思う。