遅い部分を特定する

時間はticとtocを組み合わせて計る。tic/tocのヘルプに便利な例がある。
プロファイラーを使って、コードの遅い部分を特定する。

ベクトル計算の利用

ベクトル計算とは、簡単に行ってしまうとfor文を使わないやり方。例えば、meshgridで生成した格子データを関数に放り込むとベクトル計算で処理される。。

例：

>> [Xgrid, Ygrid] = meshgrid(linspace(0,1), linspace(0,1));
>> Zvalue = cos(Xgrid .* Ygrid)

meshgridで生成した格子データなど、配列用の関数を駆使すること。

例： cusum, bsxfun, reshape, accumarray, histc, diff, repmat, permute, sparse等

これらの情報は「help->ユーザーガイド->プログラミングの基礎->パフォーマンス->パフォーマンスを向上させる方法」に記載あり

列優先

メモリが列優先であることを意識する。

列優先の良い例：

x = zeros(N, M);
for j=1:M
  for i=1:N
    x(i,j)の処理
  end
end

線形インデックス

subscript index(添え字)ではなくLinear index(線形インデックス)を用いる。

添え字を使った悪い例：

x = zeros(N, M);
for j=1:M
  for i=1:N
    x(i,j)の処理
  end
end

線形インデックスを使った良い例：

x = zeros(1,N*M);
for k=1:N*M
    x(k)の処理
end
X = reshape(N,M);

補足すると、X(i,j) = x(i+(j-1)*N)という対応になっている。

速度を比較するプログラムを作ったので、以下のプログラムを実行してみると線形インデックスの恩恵が分かる。

clear all
close all

Nx=200;
Ny=100;
REPS = 1000;
[x, y] = meshgrid(linspace(0,1,Nx), linspace(0,1,Ny));

%--------------------------------------------------------------------------
fprintf('-----------------------------------------------------\n')
fprintf('  Recommended way\n')
tic;
for t=1:REPS
    z1 = cos(x.*y);
end
averageTime = toc/REPS;
fprintf('Average time for z1: %.3e\n', averageTime);

%--------------------------------------------------------------------------
fprintf('-----------------------------------------------------\n')
fprintf('  Using subscript index\n')
tic;
for t=1:REPS
    z2 = zeros(Ny, Nx);
    for j=1:Nx
        for i=1:Ny
            z2(i,j) = cos(x(i,j)*y(i,j));
        end
    end
end
averageTime = toc/REPS;
fprintf('Average time for z2: %.3e\n', averageTime);
% fprintf('isequal(z1,z2): %d\n', isequal(z1,z2));

%--------------------------------------------------------------------------
fprintf('-----------------------------------------------------\n')
fprintf('  Using linear index for matrix\n')
tic;
for t=1:REPS
    z3 = zeros(Ny, Nx);
    for j=1:Nx
        for i=1:Ny
            index = i+(j-1)*Ny;
            z3(index) = cos(x(index)*y(index));
        end
    end
end
averageTime = toc/REPS;
fprintf('Average time for z3: %.3e\n', averageTime);
% fprintf('isequal(z1,z3): %d\n', isequal(z1,z3));

%--------------------------------------------------------------------------
fprintf('-----------------------------------------------------\n')
fprintf('  Using vector and reshape\n')
tic;
for t=1:REPS
    z4 = zeros(1, Ny*Nx);
    x4 = reshape(x,1,Ny*Nx);
    y4 = reshape(y,1,Ny*Nx);
    for i=1:Ny*Nx
        z4(i) = cos(x(i)*y(i));
    end
    z4 = reshape(z4,Ny,Nx);
    x4 = reshape(x4,Ny,Nx);
    y4 = reshape(y4,Ny,Nx);
end
averageTime = toc/REPS;
fprintf('Average time for z4: %.3e\n', averageTime);
% fprintf('isequal(z1,z4): %d\n', isequal(z1,z4));

論理的インデックスの使用

例：配列aの中から、5より大きい要素を抽出して配列bにコピーしたいとする。これを論理インデックスを使用すると以下のように書ける。

>> a=1:10;
>> b = a(a>5);

JIT/Acceleratorを意識する

MATLABはインタープリター言語なので、以前はfor文とかが苦手だったが、現在はfor文を見つけるとそこだけコンパイルしてマシン語にして実行するため高速になった。コンパイルする機能はJIT/Acceleratorとか呼ばれていて、このJITの働きを阻害するようなコーディングはよくない。例えば

ループの中で変数サイズを変更する
ループの中で変数の型を変更する
等間隔でないループインデックス

レンダリング

レンダーをMatlab独自のものから、ハードウェアレンダリングであるOpenGLに変える

>> set(gcf, 'render', 'opengl')

imakov’s blog

MATLABの高速化