Apache Commons Math3のLevenbergMarquardtOptimizerの使い方 (1変数関数)


##はじめに
[12 Optimization](http://commons.apache.org/proper/commons-math/userguide/optimization.html)の項目に書いてある内容は古いクラスを使ったものであり、将来的にDeprecatedされてしまうので、新しいクラスを用いたやり方を理解しておく必要がある。にもかかわらず、新しいクラスを使った具体的な例が見当たらず、私のような素人にはお手上げ状態である。

そこで、LM法による非線形最小二乗を行ってモデル関数のパラメータ推定を行うために、`org.apache.commons.math3.fitting.leastsquares.LevenbergMarquardtOptimizer`の働きを理解する必要がある。Javadocを眺めながらなんとか理解をしていくことを目的とする。

## LevenbergMarquardtOptimizerを使う
`org.apache.commons.math3.fitting.leastsquares.LevenbergMarquardtOptimizer`のJavadocにいろいろ書いてあるが適当に読み飛ばし、メソッドのoptimizeをみてみる。

```java
public LeastSquaresOptimizer.Optimum optimize(LeastSquaresProblem problem)
	Solve the non-linear least squares problem.
	Specified by:
		optimize in interface LeastSquaresOptimizer
	Parameters:
		problem - the problem definition, including model function and convergence criteria.
	Returns:
		The optimum.
```
結局のところ、optimizeを実行することにより非線形最小二乗問題を解くことができる。しかし、このoptimizeの形が[12 Optimization](http://commons.apache.org/proper/commons-math/userguide/optimization.html)に載っている形とだいぶ違う。

```java
PointVectorValuePair optimum = optimizer.optimize(100,
                                                   problem,
                                                   problem.calculateTarget(),
                                                   weights,
                                                   initialSolution);
```
返り値も違うし、なにより引数が大幅に違う。唯一似ている点としては、problemという名前の変数を引数として使うところだけ。javadocをたぐってoptimizeの定義を見てみると（上位クラスからの継承だった）、

```java
protected PointVectorValuePair optimize(int maxEval,
                            FUNC f,
                            OptimizationData... optData)
                                 throws TooManyEvaluationsException,
                                        DimensionMismatchException
Deprecated. 
Optimize an objective function.
Parameters:
	maxEval - Allowed number of evaluations of the objective function.
	f - Objective function.
	optData - Optimization data. The following data will be looked for:
	Target
	Weight
	InitialGuess
```
となっている。変数名problemとなっていた引数はFUNC型の関数fで、モデル関数（目的関数）を指している。Apache Commons Math3の例題プログラムでは、ここでQuadraticProblemというクラスを作ってその中でモデル関数とかヤコビアンとかいろいろ書いている。このクラスは **Deprecated**の指定があるので将来的には使えなくなる。

最新のクラスの方のoptimizerに戻ると、唯一の引数であるproblemは`LeastSquaresProblem`型らしい。この`LeastSquareProblem`型のjavadocを見てみると、

```java
The data necessary to define a non-linear least squares problem.
Includes the observed values, computed model function, and convergence/divergence criteria. Weights are implicit in LeastSquaresProblem.Evaluation.getResiduals() and LeastSquaresProblem.Evaluation.getJacobian().

Instances are typically either created progressively using a builder or created at once using a factory.
```
とあるので、古い方のクラスのoptimzeメソッドの引数の内容が、この`LeastSquareProblem`に含まれているようだ。そして、インスタンスを作るにはbuilderかfactoryを使え、とかいてある。builderの方をみてみると、`LeastSquaresBuilder`が用意されている。メソッドを眺めてみると、例えば

```java
model(MultivariateJacobianFunction newModel)
model(MultivariateVectorFunction value, MultivariateMatrixFunction jacobian)
	newModel - the model function value and Jacobian
	value - the model function value
	jacobian - the Jacobian of value
```
このあたりでモデル関数を定義するらしい。また、

```java
target(RealVector newTarget)
target(double[] newTarget)
	newTarget - the observed data.
```
これでターゲットデータ（生データ）を放り込むようだ。さらに、

```java
start(RealVector newStart)
start(double[] newStart)
	newStart - the initial guess.
``` 

これで初期値を放り込める。

EvaluationとIterationの回数は次のように設定する。

```java
maxEvaluations(int newMaxEvaluations)
	newMaxEvaluations - the maximum number of evaluations permitted.
maxIterations(int newMaxIterations)
	newMaxIterations - the maximum number of iterations permitted.
```
は古いライブラリでは個別の引数として渡してたこれらのオブジェクトを、新しいライブラリでは全てLeastSquareProblemの中に内包するようだ。

### model関数
順番に上から見ていく。modelは二種類あるが、そのうち、

```java
model(MultivariateVectorFunction value, MultivariateMatrixFunction jacobian)
```

のほうを眺めてみる。`MultivariateMatrixFunction jacobian`の方は、まんまQuadraticProblemにあるものと同じである。したがって、

```java
MultivariateMatrixFunction jacobian = new MultivariateMatrixFunction() {
	@Override
	public double[][] value(double[] point)
			 throws IllegalArgumentException {
				// TODO Auto-generated method stub
            return jacobian(point);
			}

	private double[][] jacobian(double[] variables) {
		double[][] jacobian = new double[x.size()][3];
		for (int i = 0; i < jacobian.length; ++i) {
		    jacobian[i][0] = x.get(i) * x.get(i);
		    jacobian[i][1] = x.get(i);
		    jacobian[i][2] = 1.0;
		    }
		return jacobian;
	}
			
};
```
と書くことができる。もう一方の`MultivariateVectorFunction`の方も、メソッドの`double[] value(double[] point)`が同じなので、

```java
MultivariateVectorFunction value = new MultivariateVectorFunction() {
	@Override
	public double[] value(double[] variables)
			 throws IllegalArgumentException {
		 double[] values = new double[x.size()];
		 for (int i = 0; i < values.length; ++i) {
		     values[i] = (variables[0] * x.get(i) + variables[1]) * x.get(i) + variables[2];
		 }
		 return values;
	}			
};
```
と書くことができる。

### LeastSquaresProblemをLeastSquaresBuilderで作る
LeastSquaresProblemはLeastSquaresBuilderで作ることができる。

```java
LeastSquaresBuilder lsb = new LeastSquaresBuilder();
```

さっき作ったモデル関数の値とヤコビアンは

```java
lsb.model(value, jacobian);
```
で設定できる。以下同様にして、

```java
double[] newTarget = calculateTarget();
//ターゲットデータを設定
lsb.target(newTarget);
double[] newStart = {1,1,1};
//初期値を設定
lsb.start(newStart);
//評価回数を設定
lsb.maxEvaluations(100);
//Iterationの上限を設定
lsb.maxIterations(1000);
```
モデル関数の評価回数の上限と、Iterationの上限はデフォルトでは **"0"**なので、設定をしておかないとあっさり上限を超えてエラーが出る。
`calculateTarget()`もQuadraticProblemにあるものと全く同じだ。

```java
public double[] calculateTarget() {
        double[] target = new double[y.size()];
        for (int i = 0; i < y.size(); i++) {
            target[i] = y.get(i).doubleValue();
        }
        return target;
}
```

### LevenbergMarquardtOptimizerの実装
LeastSquareProblemが”ほぼ”できたので、LevenbergMaqruardtOptimizerを作って最適化を行う。

```java
LevenbergMarquardtOptimizer lmo = new LevenbergMarquardtOptimizer();
try{
	LeastSquaresOptimizer.Optimum lsoo = lmo.optimize(lsb.build());
	final double[] optimalValues = lsoo.getPoint().toArray();			
	System.out.println("A: " + optimalValues[0]);
	System.out.println("B: " + optimalValues[1]);
	System.out.println("C: " + optimalValues[2]);
	System.out.println("Iteration number: "+lsoo.getIterations());
	System.out.println("Evaluation number: "+lsoo.getEvaluations());
} catch (Exception e) {
	System.out.println(e.toString());
}
```
1行目で普通にコンストラタで作成する。
3行目がoptimizerを実行しているところである。`LeastSquresOptimier.Optimum`型としてlsooを宣言し、optimizeメソッドを実行している。EvaluationやIterationの上限を超えるとエラーを投げるので、try~catchでくくっている。引数の`lsb.build()`は、先ほどのbuilderに設定した項目をビルドメソッドを使うことで、`LeastSquresProblem`の実体化をしている。
4行目で最適化後のパラメータの受け取りをしている。`getPoint()`の返り値は`RealVector`型なので、RealVectorのtoArrayメソッドで配列に渡している。
5行目から結果の表示だ。実行すると

```java
A: 7.99883217237272
B: 10.001841530162494
C: 16.324008168386484
Iteration number: 2
Evaluation number:2
```
となり、QuadraticProblemの例と同じ結果を得ることができた!!

### サンプルコード

作ったプログラムはGithubにアップした。
[QuadraticProblemExample2](https://github.com/arayoshipta/QuadraticProblemExample2)

## まとめ
Apache Commons Math3 (API3.4)を使ってのLevenberg Marquardtによる最適化を行うことができた。作業としては、

- モデル関数、ヤコビアンを作る
- LeastSquaresBuilderでLeastSquresProblemを作る
- LebenbergMarquardtOptimizerで最適化する

という感じ。ビルダーの代わりにファクトリーでもLeastSquresProblemを作ることができるようなので、好みで選べばいいのかな。

SimpleCurveFitterと異なり、iterationとevaluationの回数も得ることができるので、1変数関数の場合もこっちを使って記述したほうが良い気がする。

## References

- [12 Optimization](http://commons.apache.org/proper/commons-math/userguide/optimization.html)
- [反復と関数計算回数](http://jp.mathworks.com/help/optim/ug/iterations-and-function-counts.html)