[rogy Advent Calendar 2017]「うわっ…わたしのMATLABコード、遅すぎ…？」

[rogy Advent Calendar 2017]「うわっ…わたしのMATLABコード、遅すぎ…？」

MATLAB
MATLAB, Octave, 遅い, 高速化

rogy Advent Calendar 2017 17日目の記事です． 前の記事はもりゅーのV-REPとC/C++で自作ロボットを動かすでした．

「MATLABは遅い」，ほぅ…

MATLABに対する愚痴として有名なのが，「MATLABは計算が遅い」というものです． 曰く，「MATLABはインタプリタ言語だから」とか，「大規模計算向けに作られていない」とか． 実際の話をすると，MATLABはインタプリタ言語ですが最近のバージョン(R2015b〜)ではJITコンパイラを使った最適化が常時走っていますし，並列計算やGPU計算によって大規模な計算も難なくこなせるように作られています．

では，コードが遅いのはなぜなのでしょうか？MATLABの問題なのでしょうか？コードを改善するだけで実行速度を速くできないのでしょうか？ 今回はそんな問いに答えていきたいと思います．

この記事の内容の多くはOctaveにも使える話かと思いますので，フリーソフトウェア狂のOctaveユーザの方も参考にしてみてください．

また，大抵のテクニックは公式のドキュメンテーション[1]にも書いてありますが，この記事ではさらに掘り下げてマニアックな方法についても説明していきます．

なぜ遅いのか

なぜ手元のコードが遅いのか． それには大きく分けて以下の3つが考えられます．

  1. 反復回数の多いforループを多用している
  2. 配列の大きさを頻繁に変更している
  3. アルゴリズムが悪い

1. 反復回数の多いforループを多用している

MATLABにはループ処理が遅いという欠点があります． したがって，反復回数の多いforループをなるべく削減することが高速化への第一歩になります[2]．

2. 配列の大きさを頻繁に変更している

MATLABの配列(行列やベクトル)は，動的にメモリを確保できるように作られています． したがって，たとえば

x = [0 1 2 3];  % x = [0 1 2 3]
x = [x x];  % x = [0 1 2 3 0 1 2 3]

というように，配列の大きさも動的に変更することができます．

しかし，これには変数にメモリを新たに割り当て直す処理が必要となるので，通常の演算に比べて時間がかかります． 2, 3回なら問題ありませんが，これを1000回，1万回と繰り返してしまうと，全体としてパフォーマンスが低下してしまいます．

これを解消するために，zeros や ones 等のコマンドを使用して，予め必要なサイズのメモリを確保しておくことが推奨されています[3]．

3. アルゴリズムが悪い

計算の遅いアルゴリズムを使用してしまうと，実行時間は長くなります(#それはそう)． これはMATLAB側ではどうしようもありませんが，MATLABに予め実装されている高速なアルゴリズムを使って，自作の遅いアルゴリズムを置き換えることができる場合があります．

それではどのようにすればこれらの問題を解決できるのか，具体的な例を使って見ていきましょう． また，これらの問題を解決した上でさらに高速化するにはどのようなテクニックが有効かについても解説します．

方針1. forループバスターになる

アルゴリズムを実装する際にforは欠かせない構文ですから，そんなに都合よくfor文を消せないのではないか，と思うのではないかと思います． 実はその考えが落とし穴です． というのも，MATLABは行列計算をベースにした言語で，forループを使わずに様々な処理が書けるようになっています． また，MATLABは単なるプログラミング言語ではなく，豊富な標準関数も含めて1つのソフトウェアになっていて，それを活用することこそが高速化の近道だったりします．

例1-1. ベクトル・行列の要素ごとの演算

ベクトルや行列の要素ごとの演算は，要素ごとの処理を行う演算子を使用することで，forループよりも高速に処理できます[2]． 演算子以外にも，sqrt や sin などの主要な関数はベクトル・行列の入力に対応しています．

演算子
.* 要素ごとの積
./ 要素ごとの商
.^ 要素ごとの累乗

また多少気持ち悪いですが，ベクトルや行列にスカラを足すことで全要素へ足すことができるので，(1:5)+ones(1,5)*8 ではなく (1:5)+8 のように書けます．

Before

clear

n = 5000000;

Y = [];
X = linspace(0,5,n);
for k=1:n
    Y = X(k)^2-3*X(k)+4;
end

After

clear

n = 5000000;

X = linspace(0,5,n);
Y = X.^2-3*X+4;

例1-2. ディジタルフィルタ処理

ハイパスフィルタやローパスフィルタを時系列データにかける際に，forループを使って実装してしまいがちですが，素直に filter 関数[4]を使いましょう． つぎの例では，カットオフ周波数5Hzのハイパスフィルタを適当な乱数データに対して適用しています．

Before

clear
close all

n = 101;
t = linspace(0,1,n);
x = rand(n,1);
y = [];

dt = t(2)-t(1);
f_c = 5;
alpha = 1/(1+2*pi*f_c*dt);

y(1) = x(1);
for k=2:length(x)
    y = [y; alpha*(y(k-1)+x(k)-x(k-1)) ];
end

plot(t,x,t,y)

After

clear
close all

n = 101;
t = linspace(0,1,n);
x = rand(n,1);
y = [];

dt = t(2)-t(1);
f_c = 5;
alpha = 1/(1+2*pi*f_c*dt);

y = filter([alpha -alpha], [1 -alpha], x);

plot(t,x,t,y)

たとえば n=100000 程度にすると違いが大きくあらわれます．

>> tic; example_filter; toc
Elapsed time is 3.986036 seconds.
>> tic; example_filter_improved; toc
Elapsed time is 0.065683 seconds.

例1-3. 条件にマッチした要素を取り出す

ベクトルや行列から，特定の条件にマッチした要素(とそのインデックス)だけ取り出したい場合があります． 論理インデックス[5]と find 関数[6]を使用すると，forループを使わずに高速に処理できます．

Before

clear

x = randi(10,1000,1);
idx = [];
val = [];

for k=1:length(x)
    if mod(x(k), 3) == 0
        idx = [idx; k];
        val = [val; x(k)];
    end
end

After

clear

x = randi(10,1000,1);

mask = mod(x, 3) == 0; % 各要素に対して条件を満たすかどうかの論理値を計算
idx = find(mask);
val = x(mask);

例1-4. 配列要素の重複をなくす

配列要素の重複をなくしたいとき，ついついアルゴリズムを書いてしまいがちですが，unique 関数[7]を使いましょう．

Before

clear

x = [1 2 3 4 2 5 3 5 7 2 4 6 9];

y = [];
for k=1:length(x)
    if ~ismember(x(k), y)
        y = [y x(k)];
    end
end

y

After

clear

x = [1 2 3 4 2 5 3 5 7 2 4 6 9];

y = unique(x, 'stable')

実行結果はいずれも

y = 
     1     2     3     4     5     7     6     9

です． 'stable' オプションをつけずに実行すると，ソートされた状態で結果が出力されます．

y = 
     1     2     3     4     5     6     7     9

たとえば要素数を100000にして実行速度を比較すると，つぎのようになります．

>> tic; example_unique; toc
Elapsed time is 0.232489 seconds.
>> tic; example_unique_improved; toc
Elapsed time is 0.014190 seconds.

標準で用意されている関数を使おう！

上記の例で少し触れているように，MATLABには便利な関数が豊富に実装されています． 大抵の場合，用意されている関数を使用したほうが自前で実装するよりも高速に処理できます．

行列の生成と変形

使い始めたばかりの方は，単位行列や零行列，対角行列などを生成する eye や zeros, ones, diag 等[8]を確認してみましょう．

また，つぎのような関数を使いこなせると，forループを削減できる機会が多くなります．

関数
fliplr 左右反転
flipud 上下反転
reshape 行列・ベクトルのサイズの変更
repmat 行列・ベクトルの繰り返し
kron クロネッカー積

アルゴリズム

ソートや最大値・最小値，非ゼロ要素のカウントなど，様々な処理が標準で用意されています．

関数
sort ソート
sortrow 行を保存したソート
nnz 非ゼロ要素の個数
nzeros 非ゼロ要素のリスト
min, max 最小値，最大値

arrayfun, cellfun

arrayfun[9], cellfun[10] を使うと，配列の各要素に関数や無名関数(いわゆるラムダ)を適用できます． arrayfun, cellfun 自体は高速化には貢献しません[11]が，GPU配列と併用することでforループからの大幅な高速化が実現できます(後述)．

方針2. プロファイラを活用する

重い処理はないか調べる際に便利なのがプロファイラです． forループのようにあからさまなボトルネック以外でも，これを使うことで発見できる場合があります．

使用手順は以下のとおりです．

  1. プロファイリングの開始
    まず，profile on を実行してプロファイリングを開始します．
  2. プログラムの実行
  3. プロファイリングの終了，結果の閲覧
    結果を閲覧するには，profile off でプロファイリングを停止してから，profile viewer でビューワを起動します．

ビューワが起動すると，下の画像のような画面が出てきます．

関数名をクリックすると，行ごとの実行時間も確認することができます．

上記の例では，配列のサイズ変更に時間がかかっていることがわかります．

ただし，プロファイラを有効にすると通常時よりも遅く動作しますし，当然メモリも消費しますので，大きなプログラムのプロファイリングをする際には注意してください．

方針3. メモ化を活用する

何度も同じ値で関数が呼ばれる際は，memoize 関数[12]を使用してメモ化することで，処理を高速化できます． 計算結果をためておいて，以前計算した値が入ってきたら結果をそのまま出力するというイメージです．

memoize の使用例については，公式のサンプル(組み合わせの計算)が非常に良い例になっています．

ただし，memoize を経由して呼ばれた関数の処理結果はメモリ上にキャッシュされるので，十分なメモリがあるか・結果が大きな配列にならないか等を予め確認しておく必要があることに注意が必要です． キャッシュされたデータを消去するには，clearCache 関数や，clearAllMemoizedCache 関数を使いましょう．

方針4. スパース行列を使えないか検討する

MATLABの変数は既定では密行列を仮定していますが，使用する行列がスパース(つまりほとんどの要素が0)の場合は，スパース行列の機能[13]を使うことでメモリの使用量や計算時間を節約することができます． メモリの容量の都合上，通常のMATLAB変数では数万✕数万のようなサイズの行列は実用にたえませんが，スパース行列を使用することで文法はそのままに問題なく処理できるようになります．

たとえば，大規模なネットワークや階層化された構造を考える際に，グラフ理論を使用することがあります． 大抵の用途では1つのノードから伸びるエッジは少ないため，グラフを表現する隣接行列はスパースな行列になります． このような場合には，スパース行列を使用したほうがメモリ消費量も計算時間も削減できるというわけです．

方針5. 並列化する

それでもforを取り除けない場合は，処理を並列化して高速化できないかを検討すると良いでしょう． MATLABでは，forと同じ構文で使える並列化バージョンのfor文として，parfor [14]が用意されています． parfor を使用するには，Parallel Computing Toolboxが必要です．

出力結果が互いに干渉せず，かつ出力結果をつぎのループで使用しない場合であれば，基本的にそのままforをparforに置き換えるだけで並列化のご利益を受けることができます．

方針6. GPU配列を使用する

対応しているGPUと，Parallel Computing Toolboxが必要です．

わたしはつよいGPUを積んだPCを持っておらず試したことがないので，詳しくはこちらの記事をご覧ください:

大抵の場合は上記方針1〜3のテクニックで高速化できるのですが，gpuArrayを使用する場合は arrayfun や cellfun を高速化できるという恩恵が受けられます．

例6-1. ベクトル入力に対応していない関数をベクトルに適用する

配列入力に対応していない関数を配列の各要素に適用する場合，forループをどうしても使用しなければならないことがあります． 以下のコードは通常の場合はforループよりも遅くなりますが，gpuArray を使用した場合には高速になります．

Before

clear

Y = cell(1000000,1);
X = gpuArray(1:1000000);

for k=1:length(X)
    Y{k} = factor(X(k));
end

After

clear

X = gpuArray(1:1000000);
Y = arrayfun(@factor, X, 'UniformOutput', false);

方針7. JITコンパイラの特性を利用する

冒頭でも述べたように，MATLABはインタプリタ言語ですが，JITコンパイラ[15]が常時走っています． 細かい挙動に関しては公開されていないため不明ですが，どうやらサーチパスやスクリプトファイル上にある関数をロードする際にコンパイルしているようです．

以下の例は，JITの特性を利用(?)した，関数化による高速化テクニックを使用したものです．

例7-1. 1つのarrayfun, cellfunから2つ以上の結果を取り出す

たとえば，同じ数に対して factor と primes の結果を一緒に取ってきたい場合があるとします． ふつうに arrayfun でそれぞれ処理するよりも，1つの関数にした方がわずかに速くなります．

Before

clear

X = int16(rand(1000,1)*200);
Y = arrayfun(@factor, X, 'UniformOutput', false);
Z = arrayfun(@primes, X, 'UniformOutput', false);

After

clear

X = int16(rand(1000,1)*200);

[Y,Z] = arrayfun(@f,X,'UniformOutput',false);

function [out1,out2] = f(x)
out1 = factor(x); out2 = primes(x);
end

500000個程度の配列では以下のように実行時間はほとんど変わりませんが，割合としてはおよそ10%の削減になっています．

>> tic; example_arrayfun_2out; toc
Elapsed time is 10.625145 seconds.
>> tic; example_arrayfun_2out_improved; toc
Elapsed time is 9.404639 seconds.

3つ，4つと出力の数が増えていくと，この差は大きく開いていきます． ちなみに4出力にした場合はつぎのようになり，約20%の削減になっていることがわかります．

>> tic; example_arrayfun_4out; toc
Elapsed time is 21.171779 seconds.
>> tic; example_arrayfun_4out_improved; toc
Elapsed time is 17.079610 seconds.

注意として，clear all ではなく clear を使用する必要があることを補足しておきます． clear all ではコンパイル済みの関数もクリアされてしまうため，パフォーマンスが低下します([16])．

方針8. (おまけ)Simulinkの高速化

Simulinkには，実行モードという概念[17]が存在します． モードはシミュレーション時間の隣りにあるドロップダウンリストから選ぶことができます．

モード
ノーマル 通常のモード
エクスターナル ハードウェアやDesktop Real-Time等と接続して実行
アクセラレータ JITコンパイラを使用して高速化
ラピッドアクセラレータ Simulink Coderを使用して実行ファイルを作成して高速化

これらのうち，(ラピッド)アクセラレータモード[18]を使用することで，重いモデルを高速動作させられる場合があります． 実行時間がさほどかからないシミュレーションの場合は，起動にかかるオーバヘッドのほうが大きくなる可能性があることに注意が必要です．

まとめ

MATLABは一般的なプログラミング言語とは異なる特性を持っているため，高速化には特殊なテクニックが必要になることがあります．

MATLABでプログラムを書く際の基本的な考え方として，

  • forループを使わずに処理する
  • 自前実装せずに標準関数を使用する
  • プロファイラでボトルネックを探す
  • 並列化機能が使えないか試してみる というのが重要で，大抵の場合はこれらを検討するだけで大幅な高速化が見込めます．

みなさんもぜひお試しください．

  1. パフォーマンス向上の手法 - MATLAB & Simulink - MathWorks 日本 ↩

  2. ベクトル化 - MATLAB & Simulink - MathWorks 日本 ↩ ↩

  3. 事前割り当て - MATLAB & Simulink - MathWorks 日本 ↩

  4. 1 次元のデジタル フィルター - MATLAB filter - MathWorks 日本 ↩

  5. 行列のインデックス - MATLAB & Simulink - MathWorks 日本 ↩

  6. 非ゼロ要素のインデックスと値を見つける - MATLAB find - MathWorks 日本 ↩

  7. 配列の一意の値 - MATLAB unique - MathWorks 日本 ↩

  8. 行列および配列 - MATLAB & Simulink - MathWorks 日本 ↩

  9. 配列の個々の要素に関数を適用 - MATLAB arrayfun - MathWorks 日本 ↩

  10. セル配列の各セルに関数を適用 - MATLAB cellfun - MathWorks 日本 ↩

  11. array/cellfun vs. for loop - MATLAB Answers - MATLAB Central ↩

  12. メモ化のセマンティクスを関数ハンドルに追加 - MATLAB memoize - MathWorks 日本 ↩

  13. メモ化のセマンティクスを関数ハンドルに追加 - MATLAB memoize - MathWorks 日本 ↩

  14. 並列 for ループ - MATLAB parfor - MathWorks 日本 ↩

  15. JITコンパイラを使った最適化 ↩

  16. clear - MATLAB & Simulink - MathWorks 日本 ↩

  17. シミュレーション モードの選択 - MATLAB & Simulink - MathWorks 日本 ↩

  18. アクセラレータ モードの動作 - MATLAB & Simulink - MathWorks 日本 ↩

Share Comments