はじめに

スタンフォード大学には

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

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

今回の記事では、

この授業の教科書である

Introduction to Applied Linear Algebraを

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

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

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

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

また、2021年に日本語訳も出版されています。

スタンフォード　ベクトル・行列からはじめる最適化数学 (ＫＳ情報科学専門書)

本記事は、

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

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

myenigma.hatenablog.com

非線形最小二乗法

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

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

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

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

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

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

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

一般的に、

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

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

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

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

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

や

GPSによる位置計測: 理解するためのＧＰＳ測位計算プログラム入門

などがあります。　

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

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

ガウス・ニュートン法 - Wikipedia

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

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

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

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

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

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

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

最適値を計算します。

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

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

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

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

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

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

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

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

myenigma.hatenablog.com

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

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

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

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

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

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

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

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

多くの実問題の場合、

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

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

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

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

1 ヤコビ行列のノルムがが小さくなった場合
2 xの更新量が小さくなった場合、
3 最大イタレーション数を超えた場合

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

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

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

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

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

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

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

ニュートン法 - Wikipedia

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()