粗大メモ置き場

個人用,たまーに来訪者を意識する雑記メモ

MATLABで非線形最小二乗フィッティングする手順メモ

多数のデータ群から関数フィッティングを行う場合に非線形最小二乗法を用います。

MATLABにはOptimization Toolboxなどがあり,それを用いると簡単にフィッティングできます。
以下のブログなんかに解説がありますね。
MATLABで非線形最小二乗問題を解く - Py


このようなToolboxなしでも,
地道にこれを解く手順がありますので勘所をメモしておきます。
なお,数式とアルゴリズムの基礎がわかっている前提で書きます。

ダンピングを入れたVerはこちら。
ossyaritoori.hatenablog.com


アウトライン

matlabのバージョンは生協で買ったStudent version(2015a準拠)です。
本記事ではsymbolic math Toolboxを使って「偏微分」と「変数代入」を行います。
(頑張ればここは手作業でいけます。しないけど。)

ヤコビアンを求めて,Gauss-Newton法の枠組みでこれを計算します。

こんな風にフィッティング出来ます。
f:id:ossyaritoori:20170131220211p:plain

なお,Gauss-Newton法は案外簡単に発散するのでその場合は以下のLevenberg-Marquardt法を試してみてください。
ossyaritoori.hatenablog.com


問題設定

例題としてwikipedeiaの問題を考えます。
ガウス・ニュートン法 - Wikipedia


Fitting対象の関数は以下のようになります。

 y = \dfrac{ax}{x+b}

フィッティングするのはaとbの2つの変数です。

フィッティングするデータは次の様に表せます。

% data for fitting
xhat = [0.038, 0.194, 0.425, 0.626, 1.253, 2.500, 3.740].';
yhat = [0.050, 0.127, 0.094, 0.2122, 0.2729, 0.2665, 0.3317].';

通常こういう値のベクトルは縦向きですので後ろで転置してます。
今回は7×1のベクトル,すなわち7個のデータセットですね。

問題の定式化

残差と最小二乗和を以下の様に定義します。iはデータ群の各成分ってこと。

残差:  r_i = y_i - \dfrac{ax_i}{x_i+b}

残差の二乗和:  S = \sum_{i=1}^n (r_i)^2

Symbolic math Toolboxを使ってこれは次の様にかけます。

%% making evaluate funciton
syms a b x y

% evaluate fuction
func = a*x/(x+b);
R = y - func;
S = R^2;

ヤコビ行列の作成

ヤコビ行列はシンボリックの微分diffで簡単に作成できます。
(もっと簡単に作成する手法は一番下を参照)

実装上工夫する点は,

  • 代入できる値は早めに代入して,計算コストを減らす,(シンボリック計算は重い)


% partial derivative
Ra = diff(R,a);
Rb = diff(R,b);

% make jacobian matrix
Ja = subs(Ra,{x,y},{xhat,yhat});
Jb = subs(Rb,{x,y},{xhat,yhat});

Jacob = [Ja,Jb]; % jacobian

Ja,Jbの計算では,ヤコビ行列に先にフィッティングデータの値を代入して後の計算を楽にします。

補足:ヤコビ行列のサイズ

フィッティング係数がm個,データがn個あるとします。

この際,ヤコビ行列はn×mになることに注意。
f:id:ossyaritoori:20170131214013p:plain:w200
なんて親切な図だ←

注意:subsを用いた複数の配列の代入

subs関数を用いてシンボリックに配列を代入する時の書き方ですが,ややコツがいるかも?
公式を読んでいると以下のように代入しがちですがアウトです。

% 失敗例
subs(fa,[x,y],[xhat,yhat]); 

セル配列として受け渡すのが正解です。

% 成功例
subs(fa,{x,y},{xhat,yhat}); 

セル配列誰かに教わりたい。

反復計算による求解

ヤコビ行列が求まれば,その時々のパラメータの推定値 a_n,b_nをヤコビ行列Jと残差Rに代入して更新量がわかります。
パラメータをまとめてvとし,残差への代入値をrとすると,

更新式:  v_{new} = v + \Delta v
に対して

 \Delta v =  - (J^TJ)^{-1}J^T \  r

が求める計算です。
これまたsubsでの代入に注意。

% 残差に事前にフィッティングするデータを代入。
r = subs(R,{x,y},{xhat,yhat});

%ヤコビ行列にデータvを代入
J = double(subs(Jacob,[a,b],v'));
% 残差にも代入
df = double(subs(r,[a,b],v'));

% 更新量Δを計算して更新
delta = - ( J.' * J )\J.' *  df;
v = v + delta;

symsだと計算した値がシンボリックのままになるので,
doubleで一度数値に直しています。

補足:行列のサイズ

フィッティング係数がm個,データがn個あるとします。
この例題では,m=2,n=7ですね。

図解すると次のとおりになります。
f:id:ossyaritoori:20170131215154p:plain

な,なんて親切な図なんだ←


計算結果

今回は固定ステップ itaration=10 でやりました。
一般的には,「誤差が閾値を下回るまで」という条件付けをします。
正確さを気にしなければ誤差の微分見て鈍くなった所で止めてもいいかと。

フィッティングの結果を図示します。
f:id:ossyaritoori:20170131220211p:plain


各ステップ毎の誤差の二乗和Sは次のようになっています。
f:id:ossyaritoori:20170131220401p:plain

順調に減って居ますが途中から鈍くなってますね。今回最適化した結果は

a = 0.3618
b = 0.5563

です。wikiの値とほぼ同じですね。

まとめと今後

simbolicでGauss-newton法のフィッティングを行いました。
今後やるとしたら以下についてまとめる予定です。

サンプルコード

終了条件が誤差の大きさではなく,反復回数になっているので気に食わない方は変更してね。

%% get data
clear all;
% data for fitting
xhat = [0.038, 0.194, 0.425, 0.626, 1.253, 2.500, 3.740].';
yhat = [0.050, 0.127, 0.094, 0.2122, 0.2729, 0.2665, 0.3317].';

%% making evaluate funciton
syms a b x y

% evaluate fuction
func = a*x/(x+b);
R = y - func;
S = R^2;

% partial derivative
Ra = diff(R,a);
Rb = diff(R,b);

% make jacobian matrix
Ja = subs(Ra,{x,y},{xhat,yhat});
% subs(fa,[x,y],[xhat,yhat]); ではダメ
Jb = subs(Rb,{x,y},{xhat,yhat});

Jacob = [Ja,Jb]; % jacobian


%% initial value
vinit = [0.9;0.2];
v = vinit;
iteration = 10;

% evaluate function
SumOfError = zeros(iteration,1);
SEF = subs(S,{x,y},{xhat,yhat});

%% iterative optimization
% error function
r = subs(R,{x,y},{xhat,yhat});

for i = 1:iteration
   J = double(subs(Jacob,[a,b],v'));
   df = double(subs(r,[a,b],v'));
   delta = - ( J.' * J )\J.' *  df;
   v = v + delta;

   SOE = double(subs(SEF,[a,b],v'));
   SumOfError(i) = SOE.' * SOE;
end

%% showing
fit = subs(func,x,xhat);
yfit = double( subs(fit,[a,b],v') );

figure(1);
plot(xhat,yhat,'+',xhat,yfit);
xlabel('x');
ylabel('y');
legend('data','fitted value','Location','Best');
grid on;

figure(2);
plot(SumOfError);
xlabel('iteration');
ylabel('Sum Of Squared Error');
legend('Squared Error value')
grid on;

参考資料

基本的には教科書を読もう。

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

Python NumPy サンプルコード: Gauss-Newton 法で非線形最小二乗問題を解く | org-技術

jacobian 関数の使用 :2017/3/7 更新

「ヤコビ行列はシンボリックの微分diffで簡単に作成できます。」
matlabにはjacobian関数が存在し,ヤコビ行列を一発で作成してくれます。

% partial derivative
Jacob = jacobian(R,[a b]);

Referenceをちゃんと読もうな。