【MATLAB】polyfit 関数に手を加えて実行速度アップ（注：リスクあり）

やったこと

関数の実行速度を速くするために、MATLAB の関数 polyfit を書き換えて（＋別名保存して）
、実行速度が約5倍になったという話。R2019a で実施しています。

とはいえ、ルートディレクトリ（下記）にある MATLAB 関数を直接いじってしまうのはあらゆる問題の原因になりかねないので、書き換えたものは別名（例 my_polyfit）の関数として保存します。

>> matlabroot
ans =
    'C:\Program Files\MATLAB\R2019a'

保存場所は適宜関数を呼び出すスクリプトと同じフォルダに入れるなど、影響範囲は小さい方が安全ですね。

注１：MATLAB ルートディレクトリにある関数は編集しません。コピーします。
注２：掲載内容をよく読み効果と副作用などのリスクを理解したうえで活用ください。

きっかけ

MATLAB Answers にこんな質問が寄せられていました。

引用元：多次元配列に対するPolyfitの高速化 (MATLAB Answers より)

　私は測定水準x, 測定値yに対し y=a0 + a1x + a2x^2 の形でカーブフィッティングを行いたいと考えています。例えばxの次元は100万×3, yの次元は100万×3であり、１行１行について上記カーブフィッティングを行いたいです。
arrayfunを用いて下記のような関数を作成すると計算は可能でしたが実行速度が大変遅いです。実行速度を改善するために、よりよい実装方法はないでしょうか。

関数の呼び出し方で何とかなればよいのですが、残念ながら polyfit はベクトルデータに対してのみ対応となっていて、行列を与えて複数のフィッティング結果を返すことができません（R2019b 時点）。上の質問では arrayfun 関数を用いて forループは避けていますが実行速度がいまいちとのこと。

そんな質問に対して以下の解決策が見つかりました。

引用元：Multiple use of polyfit - could I get it faster? (MATLAB Answers より)

you can omit the nice and secure error checks of POLYFIT:

そう、リスクを承知の上で エラーチェックなどを消してしまえばよいと。この方法で上の質問主さんの場合は、元の実装に比べて50倍速くなったとのこと。

そもそも時間がかかるのは何が原因なのか？

実行速度が気になった場合は、まずプロファイラーを使って何がボトルネックになっているのかをチェックします。参照：プロファイラーの使い方

コードは以下（sample1.m）に書きますが、polyfit.m の内部に関しては以下のような結果になりました。
　

警告チェック部分 warnIfLargeConditionNum.. や warning 関数に、大きく時間がとられていることが分かります。

>> edit polyfit

と実行して MATLAB の提供する polyfit.m の中身も確認してみてください。

行列の条件数が悪い場合など、確かにチェックは必要ですが、その辺はちゃんと理解したうえで警告チェック部分を消しちゃっても問題ないよー、というケースもあるでしょう。多分。

次は、警告チェック部分を削除した関数を作って効果をみてみます。

プロファラーで速度検証したスクリプトです。

sample1.m

% 例として
% 乱数100組の x,y を使って
%  y=a0 + a1*x + a2*x^2 の形のカーブフィッティングを1000回繰り返します
x = rand(100,1e3);
y = rand(100,1e3);
dim = 3; 
m = size(x,2); % m = 1000;

c1 = zeros(m,dim+1); % 係数保存用
tic
for ii=1:m
   c1(ii,:) = polyfit(x(:,ii),y(:,ii),dim);
end
toc

polyfit をベースに my_polyfit 関数作成

polyfit.m　が保存されている場所は

>> which polyfit
C:\Program Files\MATLAB\R2019a\toolbox\matlab\polyfun\polyfit.m

で確認できます。この polyfit.m をコピーして、プロファイラを使って判明した時間のかかる警告チェック部分をコメントアウトし、my_polyfit.m とでも作業ディレクトリに別名保存します。mファイル内の関数名の変更も忘れずに。

さて、速くなったかな？

以下 sample2.m のように同じデータに対して処理速度を比較すると

>> sample2
経過時間は 0.294819 秒です。
経過時間は 0.067098 秒です。

となりました。速くなってますね。

速度チェック用スクリプトです。

sample2.m

% 乱数100組の x,y を使って
%  y=a0 + a1*x + a2*x^2 の形のカーブフィッティングを1000回繰り返します
x = rand(100,1e3);
y = rand(100,1e3);
dim = 3; 
m = size(x,2); % m = 1000;

c1 = zeros(m,dim+1); % 係数保存用
tic
for ii=1:m
   c1(ii,:) = polyfit(x(:,ii),y(:,ii),dim);
end
toc

c2 = zeros(m,dim+1); % 係数保存用
tic
for ii=1:m
   c2(ii,:) = my_polyfit(x(:,ii),y(:,ii),dim);
end
toc

% 一応結果の等価性を確認
isequal(c1,c2)

まとめ

リスクや副作用を気にしなくてよい状況であれば、MATLAB 側で用意されている関数もコピー＆編集することで目的が果たせる場合もあるという例でした。

polyfit 関数の場合、誤差推定も求めるために、係数を計算する際に

[Q,R] = qr(V,0);
p = R\(Q'*y);

と、いったん QR 分解をしています。
ここで、誤差推定もいらず係数だけが求まればよいということであれば

p = V\y;

でもOK。こちらはかなり速くなります。興味ある方は試してみてくださいね