MATLABで非線形最小二乗フィッティングする手順メモ
多数のデータ群から関数フィッティングを行う場合に非線形最小二乗法を用います。
MATLABにはOptimization Toolboxなどがあり，それを用いると簡単にフィッティングできます。
以下のブログなんかに解説がありますね。
MATLABで非線形最小二乗問題を解く - Py
このようなToolboxなしでも，
地道にこれを解く手順がありますので勘所をメモしておきます。
なお，数式とアルゴリズムの基礎がわかっている前提で書きます。
ダンピングを入れたVerはこちら。
ossyaritoori.hatenablog.com
アウトライン
matlabのバージョンは生協で買ったStudent version（2015a準拠)です。
本記事ではsymbolic math Toolboxを使って「偏微分」と「変数代入」を行います。
（頑張ればここは手作業でいけます。しないけど。）
ヤコビアンを求めて，Gauss-Newton法の枠組みでこれを計算します。
こんな風にフィッティング出来ます。
なお，Gauss-Newton法は案外簡単に発散するのでその場合は以下のLevenberg-Marquardt法を試してみてください。
ossyaritoori.hatenablog.com
問題設定
例題としてwikipedeiaの問題を考えます。
ガウス・ニュートン法 - Wikipedia
Fitting対象の関数は以下のようになります。
フィッティングするのは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はデータ群の各成分ってこと。
残差：
残差の二乗和：
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になることに注意。
注意：subsを用いた複数の配列の代入
subs関数を用いてシンボリックに配列を代入する時の書き方ですが，ややコツがいるかも？
公式を読んでいると以下のように代入しがちですがアウトです。
% 失敗例 subs(fa,[x,y],[xhat,yhat]);
セル配列として受け渡すのが正解です。
% 成功例 subs(fa,{x,y},{xhat,yhat});
セル配列誰かに教わりたい。
反復計算による求解
ヤコビ行列が求まれば，その時々のパラメータの推定値 をヤコビ行列Jと残差Rに代入して更新量がわかります。
パラメータをまとめてvとし，残差への代入値を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ですね。
図解すると次のとおりになります。
計算結果
今回は固定ステップ itaration=10 でやりました。
一般的には，「誤差が閾値を下回るまで」という条件付けをします。
正確さを気にしなければ誤差の微分見て鈍くなった所で止めてもいいかと。
フィッティングの結果を図示します。
各ステップ毎の誤差の二乗和Sは次のようになっています。
順調に減って居ますが途中から鈍くなってますね。今回最適化した結果は
a = 0.3618
b = 0.5563
です。wikiの値とほぼ同じですね。
まとめと今後
simbolicでGauss-newton法のフィッティングを行いました。
今後やるとしたら以下についてまとめる予定です。
- Levenberg-Marquardtアルゴリズムの実装
- 有限差分ヤコビアンを用いる -> symbolic抜きに擬似(数値)微分
サンプルコード
終了条件が誤差の大きさではなく，反復回数になっているので気に食わない方は変更してね。
%% 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;
