8. ブロードキャスト

多重ループを配列演算によってベクトル化するにはブロードキャストという技術を使います。
リスト8-1の(4)にブロードキャストの一例を示します。
i,jの2重ループにブロードキャストを適用するには、 iに関する配列aは reshape(a, n, 1) とし、 jに関する配列bは reshape(b, 1, n) とします（n : 配列の大きさ）。

リスト8-1 ブロードキャストによる多重ループのベクトル化 (bcast.m)

% test of broadcast

clearvars;

N = 20000;
fn = 1;   % = 1/2/3/4
ntest = 10;

% setup
a = zeros(N, 1, 'double');
b = zeros(N, 1, 'double');
c = zeros(N, N, 'double');
for i = 1 : N
	a(i) = double(i);
	b(i) = double(i);
end

tic

% calculation
for i = 1 : ntest
	if     fn == 1
		c = fn1(a, b);
	elseif fn == 2
		c = fn2(a, b);
	elseif fn == 3
		c = fn3(a, b);
	elseif fn == 4
		c = fn4(a, b);
	end
end

dt = toc / ntest;
fmt = '(%d) N=%d, %.3f[sec], %e, %e\n';
exact = N^2 * (N + 1);
fprintf(fmt, fn, N, dt, sum(c(:)), exact);

% (1) j, i loop
function c = fn1(a, b)
	n = length(a);
	for j = 1 : n
		for i = 1 : n
			c(i, j) = a(i) + b(j);
		end
	end
end

% (2) j loop
function c = fn2(a, b)
	n = length(a);
	for j = 1 : n
		c(:, j) = a + b(j);
	end
end

% (3) i loop
function c = fn3(a, b)
	n = length(a);
	for i = 1 : n
		c(i, :) = a(i) + b;
	end
end

% (4) broadcast
function c = fn4(a, b)
	n = length(a);
	c = reshape(a, n, 1) + reshape(b, 1, n);
end

表8-1にMATLABでの計算時間を示します。"ntest"回行い平均を取っています。
(1) MATLABは"column major"であるためにループの順序は外側をjとすることが大切です。
(2) jについてfor文のループを取り、iについて配列演算を行うと(1)より少し速くなります。
(3) iについてfor文のループを取り、jについて配列演算を行うとメモリーアクセスが非常に遅くなり、 計算時間が非常に増えます。
(4) ブロードキャストを行うとfor文がなくなってコードが単純になるだけでなく、 (2)と比べても約6倍速くなります。 ただし、プログラムが複雑になると2次元の中間変数が増え使用メモリーが増えます。 中間変数を減らすと読みにくいコードになります。

以上から多重ループでは問題に応じて、 (1)(2)(4)のいずれかを選択することがよいと思われます。