matlab-performance-optimizer

MATLAB パフォーマンス最適化ツール

ベクトル化、メモリ管理、プロファイリングツールを使って MATLAB コードのパフォーマンスを最適化します。

このスキルを使う場面

• 遅い、または非効率な MATLAB コードを最適化する
• ループをベクトル化された演算に変換する
• メモリ使用量を削減する
• アルゴリズムのパフォーマンスを向上させる
• ユーザーが「遅い」「パフォーマンス」「最適化」「高速化」「効率的」「メモリ」と言及した場合
• ボトルネックを特定するためにコードをプロファイリングする
• 計算を並列化する

最適化の基本原則

1. ベクトル化（最も重要）

可能な限り、ループをベクトル化された演算に置き換えます。

低速 — ループを使用：

% Slow approach
n = 1000000;
result = zeros(n, 1);
for i = 1:n
    result(i) = sin(i) * cos(i);
end

高速 — ベクトル化：

% Fast approach
n = 1000000;
i = (1:n).';
result = sin(i) .* cos(i);

2. 配列を事前割り当てする

ループの前に必ず配列を事前割り当てします。

低速 — 配列を拡張：

% Very slow - array grows each iteration
result = [];
for i = 1:10000
    result(end+1) = i^2;
end

高速 — 事前割り当て：

% Fast - preallocated array
n = 10000;
result = zeros(n, 1);
for i = 1:n
    result(i) = i^2;
end

3. 組み込み関数を使う

MATLAB の組み込み関数は高度に最適化されています。

低速 — 手動で実装：

% Slow
sum_val = 0;
for i = 1:length(x)
    sum_val = sum_val + x(i);
end

高速 — 組み込み関数：

% Fast
sum_val = sum(x);

ベクトル化の手法

要素ごとの演算

要素ごとの演算には .*、./、.^ を使います。

% Instead of this:
for i = 1:length(x)
    y(i) = x(i)^2 + 2*x(i) + 1;
end

% Do this:
y = x.^2 + 2*x + 1;

論理インデックス

条件付きループを論理インデックスに置き換えます。

% Instead of this:
count = 0;
for i = 1:length(data)
    if data(i) > threshold
        count = count + 1;
        filtered(count) = data(i);
    end
end
filtered = filtered(1:count);

% Do this:
filtered = data(data > threshold);

行列演算

入れ子のループではなく、行列乗算を使います。

% Instead of this:
C = zeros(size(A, 1), size(B, 2));
for i = 1:size(A, 1)
    for j = 1:size(B, 2)
        for k = 1:size(A, 2)
            C(i,j) = C(i,j) + A(i,k) * B(k,j);
        end
    end
end

% Do this:
C = A * B;

累積演算

cumsum、cumprod、cummax、cummin を使います。

% Instead of this:
running_sum = zeros(size(data));
running_sum(1) = data(1);
for i = 2:length(data)
    running_sum(i) = running_sum(i-1) + data(i);
end

% Do this:
running_sum = cumsum(data);

メモリの最適化

適切なデータ型を使う

% Instead of default double (8 bytes)
data = rand(1000, 1000);  % 8 MB

% Use single precision when appropriate (4 bytes)
data = single(rand(1000, 1000));  % 4 MB

% Use integers when applicable
indices = uint32(1:1000000);  % 4 MB instead of 8 MB

スパース行列

ほとんどの要素がゼロの行列の場合：

% Dense matrix (wastes memory)
A = zeros(10000, 10000);
A(1:100, 1:100) = rand(100);  % 800 MB

% Sparse matrix (efficient)
A = sparse(10000, 10000);
A(1:100, 1:100) = rand(100);  % Only stores non-zeros

使っていない変数をクリアする

% Process large data
largeData = loadData();
processedData = processData(largeData);

% Clear when no longer needed
clear largeData;

% Continue with processed data
results = analyze(processedData);

インプレース演算

% Instead of creating copies
A = A + 5;  % In-place when possible

% Avoid unnecessary copies
B = A;      % Creates copy if A is modified later
B = A + 0;  % Forces copy

プロファイリングとベンチマーク

プロファイラーの使用

% Profile code execution
profile on
myFunction(inputs);
profile viewer
profile off

プロファイラーでは、次の情報を確認できます。

• 各関数の実行時間
• 各関数の呼び出し回数
• 最も時間がかかる行

実行時間の比較

% Time single execution
tic;
result = myFunction(data);
elapsedTime = toc;

% Benchmark with timeit (more accurate)
timeit(@() myFunction(data))

% Compare multiple approaches
time1 = timeit(@() approach1(data));
time2 = timeit(@() approach2(data));
fprintf('Approach 1: %.6f s\nApproach 2: %.6f s\n', time1, time2);

よく使われる最適化パターン

パターン 1：find を論理インデックスに置き換える

% SLOW
indices = find(x > 5);
y = x(indices);

% FAST
y = x(x > 5);

パターン 2：repmat ではなく暗黙的拡張を使う

% SLOW - repmat to match dimensions
A = rand(1000, 5);
B = rand(1, 5);
C = A - repmat(B, size(A, 1), 1);

% FAST - implicit expansion (R2016b+)
C = A - B;

パターン 3：計算の繰り返しを避ける

% SLOW - recalculates each iteration
for i = 1:n
    result(i) = data(i) / sqrt(sum(data.^2));
end

% FAST - calculate once
norm_factor = sqrt(sum(data.^2));
for i = 1:n
    result(i) = data(i) / norm_factor;
end

% EVEN FASTER - vectorize
result = data / sqrt(sum(data.^2));

パターン 4：効率的な文字列演算

% SLOW - concatenating in loop
str = '';
for i = 1:1000
    str = [str, sprintf('Line %d\n', i)];
end

% FAST - cell array + join
lines = cell(1000, 1);
for i = 1:1000
    lines{i} = sprintf('Line %d', i);
end
str = strjoin(lines, '\n');

% FASTEST - vectorized sprintf
str = sprintf('Line %d\n', 1:1000);

パターン 5：異なるデータ型には table を使う

% Instead of separate arrays
names = cell(1000, 1);
ages = zeros(1000, 1);
scores = zeros(1000, 1);

% Use table
data = table(names, ages, scores);
% Faster access and better organization

アルゴリズム固有の最適化

畳み込みとフィルタリング

% Use built-in functions
filtered = conv(signal, kernel, 'same');
filtered = filter(b, a, signal);

% For 2D
filtered = conv2(image, kernel, 'same');
filtered = imfilter(image, kernel);

% FFT-based for large kernels (zero-pad for linear convolution)
nfft = length(signal) + length(kernel) - 1;
filtered = ifft(fft(signal, nfft) .* fft(kernel, nfft));

距離の計算

% Instead of nested loops for pairwise distances
% SLOW
n = size(points, 1);
distances = zeros(n, n);
for i = 1:n
    for j = 1:n
        distances(i,j) = norm(points(i,:) - points(j,:));
    end
end

% FAST - vectorized
distances = pdist2(points, points);

並べ替えと検索

% Presort for multiple searches
sortedData = sort(data);

% Binary search on sorted data
idx = find(sortedData >= value, 1, 'first');

% Use ismember for set operations
[isPresent, locations] = ismember(searchValues, data);

% Use unique for removing duplicates
uniqueData = unique(data);

並列計算

単純な並列ループ（parfor）

% Convert for to parfor for independent iterations
parfor i = 1:n
    results(i) = expensiveFunction(data(i));
end

parfor の要件：

• 各反復は独立している必要があります
• ループ変数は連続した整数である必要があります
• 変数は、ループ変数、スライス変数、ブロードキャスト変数、またはリダクション変数に分類できる必要があります

配列の並列演算

% Create parallel pool
parpool('local', 4);  % 4 workers

% Use parfeval for asynchronous parallel execution
futures = parfeval(@expensiveFunction, 1, data);
result = fetchOutputs(futures);

% GPU arrays for massive parallelization
gpuData = gpuArray(data);
result = arrayfun(@myFunction, gpuData);
result = gather(result);  % Bring back to CPU

高度な最適化

クリティカルな処理に MEX 関数を使う

パフォーマンスが重要なコードを C/C++ に変換します。

% Create MEX file for bottleneck function
% Write myFunction.c, then compile:
% mex myFunction.c

% Call like regular MATLAB function
result = myFunction(inputs);

キャッシュされた結果に永続変数を使う

function result = expensiveComputation(input)
    persistent cachedData cachedInput

    if isequal(input, cachedInput)
        % Return cached result
        result = cachedData;
        return;
    end

    % Compute and cache
    result = computeExpensiveOperation(input);
    cachedData = result;
    cachedInput = input;
end

JIT アクセラレーションのベストプラクティス

MATLAB の JIT（ジャストインタイム）コンパイラーは、次のコードを最適化します。

• スカラー演算を行う単純な for ループ
• 動的な機能を使わない関数

JIT に適したコード：

function result = jitFriendly(n)
    result = 0;
    for i = 1:n
        result = result + i;
    end
end

JIT に適さないコード（避ける）：

function result = jitUnfriendly(n)
    result = 0;
    for i = 1:n
        eval(['x' num2str(i) ' = i;']);  % Dynamic code
    end
end

パフォーマンスのチェックリスト

最適化したコードを仕上げる前に、次の点を確認します。

• 可能な箇所でループをベクトル化している
• ループの前に配列を事前割り当てしている
• 手動実装ではなく組み込み関数を使っている
• find とインデックス指定の組み合わせではなく、論理インデックスを使っている
• 適切なデータ型（single と double、整数）を使っている
• スパースデータにスパース行列を使っている
• 繰り返しの計算をループの外に移している
• 効率的な方法で文字列を連結している
• 実際のボトルネックを特定するためにコードをプロファイリングしている
• 要素ごとのループではなく行列演算を使っている
• 独立した演算に並列計算を検討している
• メモリを大量に使う演算を最適化している
• コストの高い処理を繰り返し呼び出す場合にキャッシュを実装している

プロファイリングのワークフロー

1. まず測定する：最適化の前にプロファイリングを行います

   profile on
   myScript;
   profile viewer
   

2. ボトルネックを特定する：最も時間のかかる関数に注目します

3. 最適化する：適切な手法を適用します

4. 再度測定する：改善を確認します

   % Before
   time_before = timeit(@() myFunction(data));
   
   % After optimization
   time_after = timeit(@() myFunctionOptimized(data));
   
   fprintf('Speedup: %.2fx\n', time_before/time_after);
   

5. 繰り返す：残りのボトルネックについて繰り返します

よくあるパフォーマンス上の落とし穴

落とし穴 1：時期尚早な最適化

• まずプロファイリングを行い、その後で最適化します
• 思い込みではなく、実際のボトルネックに注目します

落とし穴 2：過剰なベクトル化

• ループのほうが分かりやすく、十分に高速な場合もあります
• 読みやすさとパフォーマンスのバランスを取ります

落とし穴 3：メモリアクセスパターンを無視する

% SLOW - inner loop over columns (row-major traversal in column-major MATLAB)
for i = 1:rows
    for j = 1:cols
        A(i,j) = process(i, j);
    end
end

% FAST - inner loop over rows (column-major traversal, contiguous memory)
for j = 1:cols
    for i = 1:rows
        A(i,j) = process(i, j);
    end
end

% FASTEST - vectorized
[I, J] = ndgrid(1:rows, 1:cols);
A = process(I, J);

落とし穴 4：不要なデータ型変換

% SLOW - repeated conversions
for i = 1:n
    x = double(data(i));
    result(i) = sin(x);
end

% FAST - convert once
x = double(data);
result = sin(x);

最適化の例

例 1：画像処理

% SLOW
[rows, cols] = size(image);
output = zeros(rows, cols);
for i = 2:rows-1
    for j = 2:cols-1
        output(i,j) = mean(image(i-1:i+1, j-1:j+1), 'all');
    end
end

% FAST
kernel = ones(3,3) / 9;
output = conv2(image, kernel, 'same');

例 2：統計解析

% SLOW
n = size(data, 1);
means = zeros(n, 1);
for i = 1:n
    means(i) = mean(data(i, :));
end

% FAST
means = mean(data, 2);

例 3：時系列処理

% SLOW
n = length(signal);
movingAvg = zeros(size(signal));
window = 10;
for i = window:n
    movingAvg(i) = mean(signal(i-window+1:i));
end

% FAST - trailing window: [window-1 past samples, 0 future samples]
movingAvg = movmean(signal, [window-1 0]);

パフォーマンスのトラブルシューティング

問題：ベクトル化した後もコードが遅い

• 解決策：プロファイリングで新たなボトルネックを特定します。アルゴリズムの計算量も検討します

問題：メモリ不足エラーが発生する

• 解決策：より小さいデータ型を使用し、データを分割して処理し、疎行列を使用します

問題：parfor が for ループより遅い

• 解決策：オーバーヘッドが利点を上回っていないか確認し、各反復に十分な処理コストがあることを確かめます

問題：GPU 計算が CPU より遅い

• 解決策：データ転送のオーバーヘッドが計算時間を上回っている可能性があります。大規模な配列で使用します

追加リソース

• profile viewer を使用してパフォーマンスを分析します
• memory を使用してメモリ使用量を確認します
• doc で次を参照します：timeit、tic/toc、parfor、gpuArray、sparse
• MATLAB のパフォーマンスおよびメモリに関するドキュメントを確認します