MATLABで、for文・parfor文中で関数内のパラメータを様々に変えて計算をする方法(クロージャの利用)

たとえば以下のような状況があると思います。

  • ある微分方程式の解軌道を、for文中で異なるパラメータごとに描きたい
  • ある方程式の解を、for文中で異なるパラメータごとに求めたい

このようなとき、ode45やfsolveを使いますが、for文中で関数の中身のパラメータの値を直接いじることはできません。

あまりよくない方法(parfor内では使えない方法)

そのため、たとえばglobalを使ってグローバル変数を定義し、関数内から変えたりする方法を取ったりします。例として、関数odefunc内のpという変数を様々に変えてode45をする場合、次のようになります:

% odefunc.m
function ret = odefunc(t,x)
    global p;
    ret = [-p*x(2); x(1)];
end
% main.m
global p;
for p = 1:10
    [t, y] = ode45(@odefunc, linspace(0,10,100),[1; 1]); 
    figure;
    plot(t,y);
end

この方法は比較的簡単ですが、parfor内では使えないという重大な欠点があります。parforはfor文と同じ文法で簡単に並列計算・高速化を行うための方法(MATLABのparforの使い方メモ(高速化のために) - woodrush’s diary 参照)なのですが、globalな変数を用いていると並列計算が行えないため、この方法が使えません。

良い方法 - 「関数を返す関数」を使う方法(parfor内で使える方法)

これは、「関数を返す関数」を使うとかなり綺麗に解決することができます。まずはソースを見てみましょう:

% odefuncFactory.m
function retfunc = odefuncFactory(p)
    function ret = odefunc(t,x)
        ret = [-p*x(2); x(1)];
    end
    retfunc = @odefunc;
end
% main.m
parfor p = 1:10
    odefunc = odefuncFactory(p);
    [t, y] = ode45(odefunc, linspace(0,10,100), [1; 1]); 
    figure;
    plot(t,y);
end

ここでodefuncFactoryという怪しい関数が出てきました。これはその名の通り、「パラメータpを与えると、そのpを使ったodefuncを返す関数」つまり「関数を返す関数」になっています。main.mの3行目にて、odefuncFactoryに所望のpを入れ、そのpを使ったodefuncを受け取っています。「odefunc = 」と書くことで、返ってきた関数をodefuncという変数に格納しています。(このように、MATLABでは関数も変数のように扱うことができるのです。)受け取ったodefuncをode45に渡し、シミュレーションを行っています。

この方法を使うと、異なるpごとに異なる関数を毎回作成することになるので、parfor内でも使えるようになるのです。このことを強調するために、2つ目のmain.mではfor文をparfor文にしました。

おまけ - 更に簡潔に

odefuncの程度のような短い関数ならば、無名関数を使うと更に短く書くことができます:

% main.m
odefuncFactory = (@ (p) (@ (t,x) [-p*x(2); x(1)]));
for p = 1:10
    odefunc = odefuncFactory(p);
    [t, y] = ode45(odefunc, linspace(0,10,100), [1; 1]); 
    figure;
    plot(t,y);
end

% 更に短く
parfor p = 1:10
    odefunc = (@ (t,x) [-p*x(2); x(1)]);
    [t, y] = ode45(odefunc, linspace(0,10,100), [1; 1]); 
    figure;
    plot(t,y);
end

main.mだけになってしまいました。一つ目の例のように、無名関数も「関数を返す関数」の作成に使えます。更に、odefuncくらいの関数なら、odefuncFactoryを作るまでもなく、二つ目の例のように直接odefuncを無名関数で構成することもできます。
ただし、最初の例では無名関数odefuncFactoryをparfor内で定義しているという違いに注意して下さい。無名関数を使うと何かとparforの制約に引っかかりやすくなるので、通常は個別のmファイルを作ることをおすすめします。

もう少し詳しく

この一連方法では「関数を返す関数」よりも、クロージャという概念が本質的に活用されています。2つ目のmain.mを見てみましょう。main.m中には「odefuncFactory関数の中にodefunc関数が包まれている」という階層構造があることに注目して下さい。次に、odefuncFactory.mのソースをよく見ると、「function ret = odefunc(t,x)」の中には変数pの定義が含まれていない(p= ... と書いていない)ことに気づきます。普通ならばどこにもpの定義が書いていなければエラーが起こります。それでも動くのは、odefuncの外のodefuncFactoryでpが定義されているからです。odefuncは、その外のodefuncFactoryに存在する変数pの値が「覗ける」ので、このプログラムは動くわけです。つまり、odefuncFactoryは、odefuncと共に、odefuncの外側の変数たちの知識、つまり外側の環境へのアクセスも与えています。このように、関数とその外側の環境のペアをクロージャと言います。odefuncFactoryは「関数を返す関数」と言いましたが、更に正確に言うと「クロージャを返す関数」だったわけです。

このようにMATLABの無名関数はクロージャを返してくれるという性質から、色々とマニアックなことができます。たとえば、条件分岐や複数行実行などが行えたりします。詳しくはMATLABの無名関数内で条件分岐・複数行実行を行う方法 - woodrush’s diaryを御覧ください。