for文を使わない！ベクトル化による高速計算テクニック
はじめに
MATLABを使い始めたばかりの方は、他のプログラミング言語の習慣でfor文を多用してしまいがちです。しかし、MATLABの真の力は行列演算にあります。この記事では、for文を使わずにベクトル化によって計算を高速化する方法を、実例を交えて解説します。
なぜベクトル化が重要なのか
MATLABは行列演算に最適化されており、要素ごとの処理をfor文で書くよりも、行列演算として一括処理する方が圧倒的に高速です。実際の計測結果を見てみましょう。
% データサイズ
n = 1000000;
x = rand(n, 1);
y = rand(n, 1);
% for文を使った場合
tic;
result1 = zeros(n, 1);
for i = 1:n
result1(i) = x(i) * y(i) + sin(x(i));
end
time_for = toc;
% ベクトル化した場合
tic;
result2 = x .* y + sin(x);
time_vec = toc;
fprintf('for文: %.4f秒\n', time_for);
fprintf('ベクトル化: %.4f秒\n', time_vec);
fprintf('速度向上: %.1f倍\n', time_for/time_vec);
基本的なベクトル化のテクニック
1. 要素ごとの演算
配列の各要素に対して同じ演算を行う場合、ドット演算子を使用します。
% 配列Aを定義
A = [1, 2, 3, 4, 5;
6, 7, 8, 9, 10;
11, 12, 13, 14, 15];
% NG: for文を使った例
result = zeros(size(A));
for i = 1:numel(A)
result(i) = A(i)^2 + 3*A(i) - 5;
end
% OK: ベクトル化
result = A.^2 + 3*A - 5;
2. 条件に基づく処理
論理インデックスを使うことで、if文を含むfor文を回避できます。
- 例.dataの中で15より大きい値は15に置換
% データ配列を定義
data = randn(1, 1000) * 10; % 平均0、標準偏差10の正規分布
threshold = 15;
% NG: for文とif文
for i = 1:length(data)
if data(i) > threshold
data(i) = threshold;
end
end
% OK: 論理インデックス
data(data > threshold) = threshold;
3. 累積和や累積積
MATLABには累積演算用の関数が用意されています。
% サンプルデータ
x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10];
% NG: for文で累積和
cumsum_manual = zeros(size(x));
cumsum_manual(1) = x(1);
for i = 2:length(x)
cumsum_manual(i) = cumsum_manual(i-1) + x(i);
end
% OK: 組み込み関数
cumsum_result = cumsum(x);
実践例：画像処理での応用
画像の各ピクセルに対してフィルタを適用する例を見てみましょう。
% サンプル画像の作成（実際の画像の代わり）
img = rand(480, 640, 3); % 480x640のRGB画像をシミュレート
% NG: for文を使った明度調整
[rows, cols, channels] = size(img);
img_bright1 = zeros(size(img));
for i = 1:rows
for j = 1:cols
for k = 1:channels
img_bright1(i,j,k) = img(i,j,k) * 1.5;
if img_bright1(i,j,k) > 1
img_bright1(i,j,k) = 1;
end
end
end
end
% OK: ベクトル化
img_bright2 = img * 1.5;
img_bright2(img_bright2 > 1) = 1;
% または min関数を使って
img_bright3 = min(img * 1.5, 1);
高度なベクトル化テクニック
1. meshgridを使った2次元演算
2次元の格子点での計算を効率化できます。
meshgrid : １次元ベクトルを２次元,３次元座標行列に変換
% グリッドサイズを設定
grid_size = 100; % 100x100のグリッド
range = linspace(-2, 2, grid_size);
% NG: 二重for文での実装
tic;
[X, Y] = meshgrid(range, range);
Z_for = zeros(size(X));
for i = 1:size(X,1)
for j = 1:size(X,2)
Z_for(i,j) = exp(-(X(i,j)^2 + Y(i,j)^2));
end
end
time_for = toc;
% OK: ベクトル化
tic;
[X, Y] = meshgrid(range, range);
Z_vec = exp(-(X.^2 + Y.^2));
time_vec = toc;
fprintf('meshgrid演算の比較:\n');
fprintf(' for文: %.4f秒\n', time_for);
fprintf(' ベクトル化: %.4f秒\n', time_vec);
fprintf(' 速度向上: %.1f倍\n', time_for/time_vec);
% 結果が同じであることを確認
fprintf(' 結果の差: %e\n', max(abs(Z_for(:) - Z_vec(:))));
2. bsxfunやimplicit expansionの活用
異なるサイズの配列間での演算を効率的に行えます。
% 各列から列平均を引く処理の比較
A = rand(1000, 500); % 1000行500列の行列
% NG: for文での実装
tic;
A_centered1 = zeros(size(A));
for j = 1:size(A, 2)
col_mean = mean(A(:, j));
for i = 1:size(A, 1)
A_centered1(i, j) = A(i, j) - col_mean;
end
end
time_for = toc;
% OK: ベクトル化（MATLAB R2016b以降）
tic;
A_centered2 = A - mean(A, 1);
time_implicit = toc;
% OK: bsxfun（古いバージョン）
tic;
A_centered3 = bsxfun(@minus, A, mean(A, 1));
time_bsxfun = toc;
fprintf('\n列平均を引く処理の比較:\n');
fprintf(' for文: %.4f秒\n', time_for);
fprintf(' implicit expansion: %.4f秒 (%.1f倍高速)\n', time_implicit, time_for/time_implicit);
fprintf(' bsxfun: %.4f秒 (%.1f倍高速)\n', time_bsxfun, time_for/time_bsxfun);
3. 行列演算を活用した距離計算
点群間の距離計算も効率化できます。
% N個の3次元点
N = 1000;
points = rand(N, 3);
% すべての点間のユークリッド距離を計算
% NG: 二重for文
tic;
dist_for = zeros(N, N);
for i = 1:N
for j = 1:N
diff = points(i,:) - points(j,:);
dist_for(i,j) = sqrt(sum(diff.^2));
end
end
time_for = toc;
% OK: ベクトル化（pdist2関数を使用）
tic;
dist_vec1 = pdist2(points, points);
time_pdist2 = toc;
% OK: 手動でのベクトル化
tic;
% 3次元に拡張して一度に計算
points_i = reshape(points, [N, 1, 3]);
points_j = reshape(points, [1, N, 3]);
diff = points_i - points_j;
dist_vec2 = sqrt(sum(diff.^2, 3));
time_manual = toc;
fprintf('\n距離行列計算の比較:\n');
fprintf(' for文: %.4f秒\n', time_for);
fprintf(' pdist2関数: %.4f秒 (%.1f倍高速)\n', time_pdist2, time_for/time_pdist2);
fprintf(' 手動ベクトル化: %.4f秒 (%.1f倍高速)\n', time_manual, time_for/time_manual);
pdist2の意外な遅さについて
上記の結果で、pdist2がfor文より遅い場合があることに驚かれるかもしれません。この現象の主な原因は以下の通りです。
-
汎用性のコスト：pdist2は多様な距離メトリクス（ユークリッド、マンハッタン、コサイン類似度など）に対応するため、内部で距離メトリクスの判定と適切な計算ルーチンの選択を行います。この柔軟性がオーバーヘッドとなります。
-
入力検証の負荷：pdist2は入力データの妥当性チェック、NaN値の処理、次元の整合性確認など、堅牢性のための処理を含んでいます。これらの処理は、単純な距離計算には不要な場合が多いです。
-
メモリ管理の違い：pdist2は内部で最適化されたメモリ管理を行いますが、特定のデータサイズでは、この最適化が逆にオーバーヘッドになることがあります。
このため、単純なユークリッド距離の計算で、データが整形済みの場合は、手動でベクトル化した方が高速になることがあります。ただし、pdist2は複雑な距離メトリクスや大規模データでは最適化が効いてくることもあるため、実際の使用場面でベンチマークを取ることが重要です。
パフォーマンス比較：実際のベンチマーク
様々な処理でのベクトル化の効果を測定してみましょう。
% ベンチマーク関数
function benchmark_vectorization()
sizes = [1000, 10000, 100000, 1000000];
fprintf('=== MATLABベクトル化の効果測定 ===\n\n');
for n = sizes
fprintf('データサイズ: %d要素\n', n);
fprintf('----------------------------------------\n');
% テストデータの準備
A = rand(n, 1);
B = rand(n, 1);
C = rand(n, 1);
threshold = 0.5;
% 1. 複合的な数式計算
fprintf('1. 複合的な数式: A.^2 + B.*sin(C) + sqrt(abs(A-B))\n');
% for文版
tic;
result1 = zeros(n, 1);
for i = 1:n
result1(i) = A(i)^2 + B(i)*sin(C(i)) + sqrt(abs(A(i)-B(i)));
end
time_for = toc;
% ベクトル化版
tic;
result2 = A.^2 + B.*sin(C) + sqrt(abs(A-B));
time_vec = toc;
fprintf(' for文: %.4f秒, ベクトル化: %.4f秒 (%.1f倍高速)\n', ...
time_for, time_vec, time_for/time_vec);
% 2. 条件付き処理
fprintf('2. 条件付き処理: 閾値を超える要素を処理\n');
% for文版
data = rand(n, 1);
tic;
result1 = data;
for i = 1:n
if data(i) > threshold
result1(i) = threshold + (data(i) - threshold) * 0.1;
end
end
time_for = toc;
% ベクトル化版
data = rand(n, 1);
tic;
result2 = data;
mask = data > threshold;
result2(mask) = threshold + (data(mask) - threshold) * 0.1;
time_vec = toc;
fprintf(' for文: %.4f秒, ベクトル化: %.4f秒 (%.1f倍高速)\n', ...
time_for, time_vec, time_for/time_vec);
% 3. 行列の各要素に対する複雑な演算
fprintf('3. 行列演算: 各要素の累乗と対数の組み合わせ\n');
% for文版
M = rand(round(sqrt(n)), round(sqrt(n)));
tic;
result1 = zeros(size(M));
for i = 1:size(M,1)
for j = 1:size(M,2)
if M(i,j) > 0.1
result1(i,j) = M(i,j)^3 + log(M(i,j));
else
result1(i,j) = M(i,j)^2;
end
end
end
time_for = toc;
% ベクトル化版
tic;
result2 = zeros(size(M));
mask = M > 0.1;
result2(mask) = M(mask).^3 + log(M(mask));
result2(~mask) = M(~mask).^2;
time_vec = toc;
fprintf(' for文: %.4f秒, ベクトル化: %.4f秒 (%.1f倍高速)\n', ...
time_for, time_vec, time_for/time_vec);
% 4. 移動平均の計算
fprintf('4. 移動平均: 窓サイズ10の移動平均\n');
window_size = 10;
data = rand(n, 1);
% for文版
tic;
result1 = zeros(n-window_size+1, 1);
for i = 1:(n-window_size+1)
result1(i) = mean(data(i:i+window_size-1));
end
time_for = toc;
% ベクトル化版（convolution使用）
tic;
kernel = ones(window_size, 1) / window_size;
result2 = conv(data, kernel, 'valid');
time_vec = toc;
fprintf(' for文: %.4f秒, ベクトル化: %.4f秒 (%.1f倍高速)\n', ...
time_for, time_vec, time_for/time_vec);
fprintf('\n');
end
% 特別なケース：画像フィルタリング
fprintf('=== 実用例: 画像フィルタリング ===\n');
img_size = 500;
img = rand(img_size, img_size);
kernel = [-1 -1 -1; -1 8 -1; -1 -1 -1]; % エッジ検出フィルタ
% for文版
tic;
result1 = zeros(img_size-2, img_size-2);
for i = 2:img_size-1
for j = 2:img_size-1
temp = 0;
for ki = -1:1
for kj = -1:1
temp = temp + img(i+ki, j+kj) * kernel(ki+2, kj+2);
end
end
result1(i-1, j-1) = temp;
end
end
time_for = toc;
% ベクトル化版（conv2使用）
tic;
result2 = conv2(img, kernel, 'valid');
time_vec = toc;
fprintf('画像サイズ %dx%d:\n', img_size, img_size);
fprintf(' for文: %.4f秒, ベクトル化: %.4f秒 (%.1f倍高速)\n', ...
time_for, time_vec, time_for/time_vec);
end
ベクトル化が難しい場合の対処法
すべての処理がベクトル化できるわけではありません。以下のような場合は工夫が必要です。
- 再帰的な処理: 前の結果に依存する計算
- 不規則なアクセスパターン: ランダムなインデックスへのアクセス
- 複雑な条件分岐: 多段階の条件判定
これらの場合でも、部分的なベクトル化や、arrayfun、cellfunなどの関数を使って改善できることがあります。
% 部分的なベクトル化の例
% フィボナッチ数列（完全なベクトル化は困難）
n = 100;
fib = zeros(n, 1);
fib(1:2) = 1;
% 最初の2項以外はfor文が必要
for i = 3:n
fib(i) = fib(i-1) + fib(i-2);
end
% ただし、行列の累乗を使った高速化は可能
F = [1 1; 1 0];
fib_matrix = F^n;
fib_n = fib_matrix(1, 2);
まとめ
MATLABでのベクトル化は、コードの実行速度を劇的に向上させる重要なテクニックです。主なポイントは以下の通りです。
- for文の代わりに行列演算を使用する
- 論理インデックスを活用して条件分岐を回避する
- 組み込み関数（sum、mean、cumsum等）を積極的に利用する
- meshgridやbsxfunで多次元演算を効率化する
ベクトル化されたコードは、高速であるだけでなく、より簡潔で読みやすくなることが多いです。MATLABを使う際は、常に「この処理はベクトル化できないか？」と考える習慣をつけましょう。