非線形フィッティング Gauss-Newton法

最適化
Gauss-Newton法

概要

Gauss-Newton法とは関数の非線形フィッティングで用いられる手法の一つで関数の勾配を用い誤差を最小化する。
初期値によって誤差の局所最小値にすら収束しないことがある。

更新式

Gauss-Newton法におけるフィッティングパラメータ\( \boldsymbol{r}\)の更新式は以下のようになる。
\( \quad \boldsymbol{r_{n+1}} = \boldsymbol{r_{n}} - \lambda J^{-1} \boldsymbol{e} \)

ここで\(J \)は\( \boldsymbol{r_{n}} \)におけるヤコビアン行列、\(\boldsymbol{e} \)は誤差(近似値-真値)である。
\( \lambda \)は学習定数で\(0 \lt \lambda \leq 1 \)とする。

ソースコード


namespace DataFitting {

    /// <summary>Gauss-Newton法</summary>
    public class GaussNewtonMethod : FittingMethod{
        readonly FittingFunction func;

        /// <summary>コンストラクタ</summary>
        public GaussNewtonMethod(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 = 0.75, int loop = 32) {
            Vector errors, dparam;
            Matrix jacobian;

            for(int j = 0; j < loop; j++) {
                errors = Error(parameters);
                jacobian = Jacobian(parameters);
                dparam = jacobian.Inverse * errors;

                if(!Vector.IsValid(dparam)) {
                    break;
                }

                parameters -= dparam * lambda;
            }

            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 GaussNewtonMethodTests {
        [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);
            }

            GaussNewtonMethod fitting = new GaussNewtonMethod(data_list, new FittingFunction(2, fitting_func, fitting_diff_func));

            Assert.AreEqual((fitting.ExecureFitting(new Vector(3, 4)) - p).Norm, 0, 1e-10);
        }
    }
}