posted with カエレバ

金森敬文,鈴木大慈,竹内一郎,佐藤一誠講談社 2016-12-07

Amazonで探す

楽天市場で探す

Yahooショッピングで探す

はじめに

スタンフォード大学には

機械学習を学ぶ上での第一歩として、

Introduction to Matrix Methods (EE103)という授業があります。

今回の記事では、

この授業の教科書である

Introduction to Applied Linear Algebraを

読んだ際の技術メモです。

この教科書は下記のリンクのページから

pdfをダウンロードすることができます。

Introduction to Applied Linear Algebra – Vectors, Matrices, and Least Squares

本記事では、

上記の教科書の非線形最小二乗法の部分のみのメモです。

他の部分に関しては、下記の記事を参照下さい。

myenigma.hatenablog.com

非線形最小二乗法の解法2: レーベンバーグ・マーカート法

レーベンバーグ・マーカート法の概要

レーベンバーグ・マーカート法は、

こちらの記事で紹介した

myenigma.hatenablog.com

ガウス・ニュートン法の問題点を解決するための

アルゴリズムです。

Levenberg–Marquardt algorithm - Wikipedia

レーベンバーグ・マーカート法は

画像処理における非線形最適化などで広く利用されています。

A Levenberg–Marquardt Scheme for Nonlinear Image Registration | SpringerLink

ガウス・ニュートン法の一番の問題点として、

近似点から遠い場所に補正した場合、残差が大きくなる方法に補正することがある。（本当は小さくしたい)

というのがあります。

前述の通り、ガウス・ニュートン法は、一つ前の状態量を使って、

その周辺を線形化しているため、その状態量から離れた地点では、

線形近似誤差が大きくなってしまうためです。

そこでレーベンバーグ・マーカート法では、

ヤコビ行列の二乗を最小化するだけではなく、

近似点からの距離も評価関数に組み込み、

複数のコスト関数の最小二乗法を解くことにより、

元の近似点から離れすぎないようなxの更新を実施します。

f:id:meison_amsl:20181015111933p:plain

下記がレーベンバーグ・マーカート法の更新式です。

f:id:meison_amsl:20181015111258p:plain

この式において、

ラムダは元の近似点から離れすぎないようにする

重みパラメータです。

このパラメータが大きすぎると、

最適値の近くで、最適値へ向かう速度が遅くなり、

小さすぎると、

先ほどのガウスニュートン法の問題点が発生してしまいます。

また、上記の更新式で単位行列の部分を、

ヤコビ行列の二乗の対角行列にする方法もよく使われます。

Levenberg–Marquardt algorithm - Wikipedia

そこで、レーベンバーグ・マーカート法では、

ラムダの値を、

最適化中に自動チューニングするようになっています。

一般的には、最適化のイタレーション中に、

コスト関数が改善した(小さくなった)場合、

ラムダを小さくし、(0.8倍するなど)

コスト関数が悪化した(大きくなった)場合、

ラムダを大きくします。(2倍するなど)

これにより、最適化が高速化し、かつ、安定化します。

また、非線形最小二乗法をより効率的に解く方法として、

Warm start: 過去に同じように最適化を解いた解を、初期値に使う方法
Multiple runs: 複数の初期値を使って、最適化を解き、最もf(x)が小さいものを選ぶ方法

なども併用されます。

Juliaによるレーベンバーグ・マーカート法による非線形最小二乗法

下記は、レーベンバーグ・マーカート法による

非線形最小二乗法のJuliaコードです。

下記の例では、非線形関数であるシグモイド関数を使って

最適化をしています。

"""
  solve nonlinear least square with levenberg marquardt
  xhat = argmin(|f(x)|^2)
"""
function solve_nonlinear_least_square_with_levenberg_marquardt(
        f, Df, x0, lambda0;
        kmax = 20, tol = 1e-6, lamr_n = 0.8, lamr_p = 2.0)

    n = length(x0)
    x=x0
    lambda = lambda0
    obj = zeros(0,1)
    residuals = zeros(0,1)
    xhist = x0'
    for k = 1:kmax
        fx = f(x)
        Dfk = Df(x)
        obj = [obj; norm(fx)^2]
        res = norm(2*Dfk'*fx)
        residuals = [residuals; res]
        if res < tol break end;
        tmp = Dfk'*fx
        if isa(tmp, Array)
            dx = (Dfk'*Dfk+lambda*eye(n)) \ tmp
        else
            dx = (Dfk'*Dfk+lambda*eye(n)) \ [tmp]
        end
        xt = x - dx
        if norm(f(xt)) < norm(fx)
            lambda = lamr_n*lambda
            x = xt
        else
            lambda = lamr_p*lambda
        end
        xhist = [xhist; x']
    end

    return x, Dict([ ("objectives", obj), ("residuals", residuals), ("x_history", xhist)])
end

function test_solve_nonlinear_least_square_with_levenberg_marquardt()
    println("test_solve_nonlinear_least_square_with_levenberg_marquardt")
    f(x) = (exp.(x) - exp.(-x)) ./ (exp.(x) + exp.(-x))
    Df(x) = 4.0 ./ (exp.(x) + exp.(-x)).^2
    xhat = solve_nonlinear_least_square_with_levenberg_marquardt(f,Df,[0.95], 1.0)
end


test_solve_nonlinear_least_square_with_levenberg_marquardt()

応用例1: 非線形最小二乗法による等距離地点探索

下記のJuliaのコードは、

二次元平面上の複数の既知の点から

等距離にある点を非線形最小二乗法で探索するサンプルコードです。

前述のレーベンバーグ・マーカート法の関数を使っています。

#
# location from range measurement code 
#
# author: Atsushi Sakai
#

using Base.Test
using PyPlot

include("../argmin.jl")
using argmin

function main()
    println(PROGRAM_FILE," start!!")
    A = [ 1.8  2.5;
          2.0  1.7;
          1.5  1.5;
          1.5  2.0;
          2.5  1.5]

    x0 = [3.0, 1.0]

    dist(x) = sqrt.( (x[1] .- A[:,1]).^2 + (x[2] .- A[:,2]).^2 )
    f(x) = dist(x)
    Df(x) = diagm(1 ./ dist(x)) * [ (x[1] .- A[:,1])  (x[2] .- A[:,2]) ]

    xhat, status = solve_nonlinear_least_square_with_levenberg_marquardt(f,Df,x0, 1.0)
    println(xhat)

    plot(A[:,1], A[:,2], "*r")
    plot(status["x_history"][:,1], status["x_history"][:,2], "-k")
    plot(xhat[1], xhat[2], "ob")
    axis("equal")
    show()

    println(PROGRAM_FILE," Done!!")
end

@time main()

上記を実行すると、

下記のようにそれぞれの点から等距離の点を探索してくれます。

f:id:meison_amsl:20181021090439p:plain

応用例2: 非線形フィッティング

続いて、レーベンバーグ・マーカート法を利用した

非線形フィッティング応用例を紹介します。

今回は、下記のような非線形モデルにおいて、

f:id:meison_amsl:20181021101229p:plain

4つのパラメータθを非線形最適化で推定します。

元の曲線線の一部の点を、入力データとして、

レーベンバーグ・マーカート法で最適化を実施するのが

下記のコードです。

#
# Non linear model fitting 
#
# author: Atsushi Sakai
#

include("../../argmin.jl")

using PyPlot
using argmin


function calc_y(theta, x)
    y = theta[1]*exp.(theta[2]*x) .* cos.(theta[3]*x .+ theta[4]);
    return y
end

function main()
    println(PROGRAM_FILE," start!!")
    theta_true = [1, -0.2, 2*pi/5, pi/3] # true parameters

    # Input data 
    M = 30
    xd = [5*rand(M); 5 .+ 15*rand(M)]
    yd = calc_y(theta_true, xd)

    f(theta) =  theta[1] * exp.(theta[2]*xd) .* cos.(theta[3] * xd .+ theta[4]) - yd;
    Df(theta) = hcat( exp.(theta[2]*xd) .* cos.(theta[3] * xd .+ theta[4]),
                            theta[1] * ( xd .* exp.(theta[2]*xd) .* cos.(theta[3] * xd .+ theta[4])),
                            -theta[1] * ( exp.(theta[2]*xd) .* xd .* sin.(theta[3] * xd .+ theta[4])),
                            -theta[1] * ( exp.(theta[2]*xd) .* sin.(theta[3] * xd .+ theta[4])) );
    theta0 = [1, 0, 1, 0]; # initial parameters
    x = linspace(0, 20, 500);
    y1 = calc_y(theta0, x)
    
    theta, history = solve_nonlinear_least_square_with_levenberg_marquardt(f, Df, theta0, 1.0)

    println("theta_true:", theta_true)
    println("theta:",theta)
    y = calc_y(theta, x)
    plot(xd, yd, "ob", label="Input data")
    plot(x, y1, "--r", label="Initial param")
    plot(x, y, "-r", label="Optimized param")
    legend()
    show()

    println(PROGRAM_FILE," Done!!")
end

@time main()