非線形フィッティング Levenberg-Marquardt法
概要
Levenberg-Marquardt法とは関数の非線形フィッティングで用いられる手法の一つで関数の勾配を用い誤差を最小化する。
Gauss-Newton法にくらべ収束速度が遅いが、初期値に依らず安定したフィッティングができる。(ただし、あまりに最適値から離れているとうまくいかない。)
更新式
Levenberg-Marquardt法におけるフィッティングパラメータ\( \boldsymbol{r}\)の更新式は以下のようになる。
\( \quad \boldsymbol{r_{n+1}} = \boldsymbol{r_{n}} - (J^T J + \lambda_{n} I)^{-1} J^T \boldsymbol{e} \)
ここで\(J \)は\( \boldsymbol{r_{n}} \)におけるヤコビアン行列、\(I \)は単位行列、\(\boldsymbol{e} \)は誤差(近似値-真値)である。
\( \lambda_{n} \)は学習定数で\(0 \lt \lambda \leq 1 \)とし、勾配の大きさによって変える。
ソースコード
namespace DataFitting {
/// <summary>Levenberg-MarquardtMethod法</summary>
public class LevenbergMarquardtMethod : FittingMethod{
readonly FittingFunction func;
/// <summary>コンストラクタ</summary>
public LevenbergMarquardtMethod(FittingData[] data_list, FittingFunction func) : base(data_list, func.ParametersCount){
this.func = func;
}
/// <summary>フィッティング値</summary>
public override double FittingValue(double x, Vector parameters) {
return func.F(x, parameters);
}
/// <summary>フィッティング</summary>
public Vector ExecureFitting(Vector parameters, double lambda_init = 1, double lambda_decay = 0.9, int loop = 64) {
Vector errors, dparam;
Matrix jacobian, jacobian_transpose;
double lambda = lambda_init;
for(int j = 0; j < loop; j++) {
errors = Error(parameters);
jacobian = Jacobian(parameters);
jacobian_transpose = jacobian.Transpose;
dparam = (jacobian_transpose * jacobian + lambda * Matrix.Identity(ParametersCount)).Inverse * jacobian_transpose * errors;
if(!Vector.IsValid(dparam)) {
break;
}
parameters -= dparam;
lambda *= lambda_decay;
}
return parameters;
}
/// <summary>ヤコビアン行列</summary>
private Matrix Jacobian(Vector parameters) {
FittingData data;
Matrix jacobian = new Matrix(data_list.Length, func.ParametersCount);
for(int i = 0, j; i < data_list.Length; i++) {
data = data_list[i];
Vector df = func.DiffF(data.X, parameters);
for(j = 0; j < parameters.Dim; j++) {
jacobian[i, j] = df[j];
}
}
return jacobian;
}
}
}
単体テスト
namespace DataFitting.Tests {
[TestClass()]
public class LevenbergMarquardtTests {
[TestMethod()]
public void ExecureFittingTest() {
double[] x_list = { 1, 3, 4, 7 };
Vector p = new Vector(2, 3);
Func<double, Vector, double> fitting_func = (x, parameter) => {
double a = parameter[0], b = parameter[1];
return b / (1 + Math.Exp((-a) * x));
};
Func<double, Vector, Vector> fitting_diff_func = (x, parameter) => {
double a = parameter[0], b = parameter[1];
double v = Math.Exp(-a * x) + 1;
return new Vector((b * x * Math.Exp(-a * x)) / (v * v), 1 / v);
};
FittingData[] data_list = new FittingData[x_list.Length];
for(int i = 0; i < x_list.Length; i++) {
double x = x_list[i];
data_list[i].X = x;
data_list[i].Y = fitting_func(x, p);
}
LevenbergMarquardtMethod fitting = new LevenbergMarquardtMethod(data_list, new FittingFunction(2, fitting_func, fitting_diff_func));
Assert.AreEqual((fitting.ExecureFitting(new Vector(7, 2)) - p).Norm, 0, 1e-10);
}
}
}
関連項目
関数フィッティング手法基本クラス
ベクトルクラス
行列クラス
直線フィッティング
多項式フィッティング
非線形フィッティング Gauss-Newton法
重み付き直線フィッティング
重み付き多項式フィッティング
最小二乗法におけるロバスト推定