MyEnigma

とあるエンジニアのブログです。#Robotics #Programing #C++ #Python #MATLAB #Vim #Mathematics #Book #Movie #Traveling #Mac #iPhone

スタンフォード大学の学生が学ぶ、非線形最小二乗法とその応用1:ガウス・ニュートン法編

目次

 

はじめに

スタンフォード大学には

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

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

 

今回の記事では、

この授業の教科書である

Introduction to Applied Linear Algebraを

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

 

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

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

 

本記事は、

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

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

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

 

非線形最小二乗法

非線形最小二乗法は下記のように定式化できます。

f:id:meison_amsl:20181014083447p:plain:w300

以前の記事で取り扱っていた線形最小二乗法は、

上記の目的関数f(x)が Ax=bのような線形関数の場合です。

つまり、通常の線形最小二乗法は非線形最小二乗法の一部になります。

しかし、非線形最小二乗法では、

f(x)はアフィン変換ではなくなり、

これまでのような直接的な解法は使えなくなります。

 

一般的に、

非線形最小二乗法は、計算時間が長く、

安定して良い解を計算することが非常に難しいですが、

ヒューリスティックな方法を使うことで解くことができます。

 

非線形最小二乗法の応用例としては、

  • 非線形データフィッティング

などがあります。  

非線形最小二乗法の解法1: ガウス・ニュートン法

ガウス・ニュートン法は非線形最小二乗法を解く代表的な方法です。

 

基本的には、ある初期値からスタートして、

現在のxの周りでテイラー展開をして、線形最小二乗法を解き、

その解を使って、xを逐次的に更新して非線形最小二乗問題を解くアルゴリズムです。

 

もう少し詳しく説明すると、

ガウス・ニュートン法では、

関数f(x)の最小値を計算するために、

テイラー展開で線形化した下記の式の二乗を最小化するように

最適値を計算します。

f:id:meison_amsl:20181015112618p:plain

 

では具体的に、ガウス・ニュートン法を導出してみたいと思います。

まず、関数f(x)をテイラー展開の二次項までで近似すると、

ここで⊿xで微分するし、式を整理すると、

となります。これをベクトルで拡張すると、

となり、Hはヘッセ行列で、gは勾配ベクトルになります。

このδxの分だけ、xを更新していけば、f(x)を最小化することができます。

 

ちなみこのように、ヘッセ行列を使って解く手法が、

ニュートン法(ガウスが付かない)です。(複雑ですね。。。笑)

myenigma.hatenablog.com

 

一方、ガウス・ニュートン法では、

ヘッセ行列Hをヤコビ行列の二乗で近似し、

勾配ベクトルはヤコビ行列と関数f(x)の出力のベクトルの積で計算できるとすると、

xの更新式は下記の通りになります。

f:id:meison_amsl:20181014092950p:plain

Df(x)がテイラー展開によって近似されたヤコビ行列です。

 

つまり、ニュートン法のヘッセ行列を

ヤコビ行列の二乗で近似したのが、

ガウス・ニュートン法であるといえます。

多くの実問題の場合、

ヘッセ行列を計算することは難しかったり、

大規模な逆行列を計算することが難しいことが多いため、

その代わりとしてガウス・ニュートン法が使われることが多いようです。

 

上記の更新式でxを更新しつつ、下記の3つ条件を元に停止判定を実施します。

  • 1 ヤコビ行列のノルムがが小さくなった場合

  • 2 xの更新量が小さくなった場合、

  • 3 最大イタレーション数を超えた場合

 

また、線形方程式の数m と変数の数nが同じ場合、

つまりヤコビ行列が正方行列の場合、

前述の更新式はもっと簡素化することができます。

f:id:meison_amsl:20181014100153p:plain

これはヤコビ行列の二乗の逆行列とヤコビ行列の掛け算で

正方行列では互いに相殺させることができるからです。

この特別な形のガウス・ニュートン法は

ニュートン・ラフソン法と呼ばれます。(またニュートン出てきた。。)

 

Juliaによるガウス・ニュートン法のサンプルコード

上記のガウスニュートン法とニュートン・ラフソン法を

Juliaで実装したのが下記です。

"""
  solve nonlinear least square with newton-raphson method

  The inputs have to be length(x) == length(f(x)).
  If it is not, you cant use solve_nonlinear_least_square_with_gauss_newton

  xhat = argmin(|f(x)|^2)
"""
function solve_nonlinear_least_square_with_newton_raphson(
        f, Df, x1; kmax = 20, tol = 1e-6)

    x=x1
    @assert length(x) == length(f(x))
    for k = 1:kmax
        fk = f(x)
        dx = fk / Df(x)
        if norm(dx) < tol break end;
        x = x - dx
    end

    return x
end

"""
  solve nonlinear least square with gauss-newton method

  xhat = argmin(|f(x)|^2)
"""
function solve_nonlinear_least_square_with_gauss_newton(
        f, Df, x1; kmax = 20, tol = 1e-6)
    x=x1
    for k=1:kmax
        fk = f(x)
        dfx = Df(x)
        dx = (dfx'*dfx) \ (dfx'*fk)
        if norm(dx) < tol break end;
        x = x - dx'
    end

    return x
end


function test_solve_nonlinear_least_square_with_newton_raphson()
    println("test_solve_nonlinear_least_square_with_newton_raphson")

    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_newton_raphson(f,Df,0.95)
end


function test_solve_nonlinear_least_square_with_gauss_newton()
    println("test_solve_nonlinear_least_square_with_gauss_newton")

    f(x) = (exp.(x) - exp.(-x)) / (exp.(x) + exp.(-x))
    Df(x) = 4 ./ (exp.(x) + exp.(-x)).^2
    xhat = solve_nonlinear_least_square_with_gauss_newton(f,Df,[0.95 1.15]')
end

test_solve_nonlinear_least_square_with_newton_raphson()
test_solve_nonlinear_least_square_with_gauss_newton()

github.com

  

参考資料

www.slideshare.net

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

 

MyEnigma Supporters

もしこの記事が参考になり、

ブログをサポートしたいと思われた方は、

こちらからよろしくお願いします。

myenigma.hatenablog.com