matlab ベクトル化入門


matlab高速化の常套手段「ベクトル化」をまとめてみました。
あくまで入門レベルです。

ベクトル化とは?

簡単に言うと、スカラ値を行列に置き換えることです。
これにより、計算速度を劇的に上げることができます。

やってみる

簡単なベクトル化。まずは配列で。

A = [1:100] ; % 1~100の配列
f = @(x) 2*x + 1 ; % テストのための計算式(関数ハンドルで生成。)
B = f(A') ; % ベクトル化計算(なんとなく列ベクトルのほうがイメージしやすいと思うので、「'」でAを転置。)

式で書くと、結果の行列Bはこんな感じで生成されます。

B = 2*
\begin{pmatrix}
1 \\
2 \\
... \\
100
\end{pmatrix}
+1=
\begin{pmatrix}
3 \\
4 \\
... \\
201
\end{pmatrix}

要素ごとに計算してくれます。ループで一要素ずつ計算するよりずっと早いです。
(いわゆるGPUの並列処理のような感じなのかな?)
結果も要素ごと、つまりAと同じサイズの配列でBが生成されます。

ちなみに、これをループで書くと、こんな感じ。

B = size(size(A)) ;
for i = 1:numel(A)
    B(i) = f(A(i)) ;
end

もちろん、配列だけでなく行列もできます。

C = randi(10, 5, 5) ; % 1~10の整数のランダムな5x5行列を生成
D = f(C) ;
D ;

DはCと同じサイズの行列として生成されます。
なんて素敵なんだ・・・!

行列と行列で

行列サイズが同じなら、複数の変数をベクトル化できます。
例えば、ランダムな値の積をとろうとしたとき、

E=randi(10,5,5);
F=randi(10,5,5);
G=E.*F;

これで、行列EとFの対応する要素の積を一気に計算できます。
この「.*」とは、要素ごとの積を意味します。「*」だけなら、行列の積です。
(「.*」を使う際は、EとFの行列サイズが同じであることが前提です)
で、これをループで書くと、

G = size(5, 5) ;
for i = 1:5
    for j = 1:5
        G(i, j) = E(i, j)*F(i, j) ;
    end
end

長い!
計算時間以前に、コードが長くなりますね。

ちょっとステップアップ

meshgridを組み合わせてみます。
X, Yの全ての組み合わせを一気に計算したいとき、

X = [1:50];
Y = [20:30] ;
[MX, MY] = meshgrid(X, Y) ;
g = @(x, y) x + y ;
M = g(MX, MY) ;

これはつまり、
MX: 10x50, MY: 10x50

MX = 
\begin{pmatrix}
1 & 2 & .. & 50\\
1 & 2 & .. & 50\\
... \\
1 & 2 & .. & 50\\
\end{pmatrix}
, MY=
\begin{pmatrix}
20 & 20 & ... & 20\\
21 & 21 & ... & 20\\
... \\
30 & 30 & ... & 30\\
\end{pmatrix}\\
M=
\begin{pmatrix}
1+20 & 2+20 & ... & 50+20\\
1+21 & 1+21 & ... & 50+20\\
... \\
1+30 & 1+30 & ... & 50+30\\
\end{pmatrix}

repmatを使ってMXとMYをそれぞれ生成することもできますが、meshgridなら一発で組み合わせ用の行列を用意できます。

さいごに

ベクトル化の利点は、高速化もさることながら、スカラ値向けの計算式をほぼそのまま使えることです。
式がそのままなので、ロジックを理解する手助けになるのではと感じています。
もちろん、前述の通り、積は「.*」か「*」かをはっきりさせる必要はありますが。

では!