`pagefun` による GPU 上の小さな行列問題のパフォーマンス改善

ライブスクリプトを開く

この例では、多くの独立した回転と平行移動を 3 次元環境のオブジェクトに適用するときに pagefun を使用してパフォーマンスを改善する方法を説明します。これは、小さな配列で大規模な計算バッチを実行する場合に関係するさまざまな問題の典型例です。

GPU は、大規模な行列で計算を実行する場合に最も効果的です。MATLAB® で通常これを実現するには、コードをベクトル化して、各命令で実行される作業を最大化します。大きなデータセットを多くの小さな行列演算に分かれた計算で処理しようとするとき、何百個もの GPU コアで同時に計算を実行してパフォーマンスを最大化するのは難しい場合があります。

関数 arrayfun を使用することで、スカラー演算を並列に GPU で実行できます。関数 pagefun は、同様に行列演算をバッチで実行する機能を追加します。関数 pagefun は、Parallel Computing Toolbox™ で gpuArray と共に使用することができます。

この例では、ロボットはセンサーを使って識別できる多数の対象物がある既知のマップを移動しています。ロボットは、オブジェクトの相対位置と向きを測定し、それらをマップ位置と照らし合わせることでマップ上のロボット自身の位置を特定します。ロボットが完全には位置を見失っていないと仮定すれば、両者の間のいかなる差異も利用し、カルマンフィルターなどを使用することによって位置を修正できます。アルゴリズムの初めの部分に焦点を当てて説明します。

マップの設定

大きな部屋で、位置と向きをランダムに配置したオブジェクトのマップを作成しましょう。

numObjects = 1000;
roomDimensions = [50 50 5]; % Length * breadth * height in meters

位置と向きを、3 行 1 列のベクトル T と 3 行 3 列の回転行列 R を使用して表します。これらの "変換" が N 個ある場合は、平行移動を 3 行 N 列の行列にまとめ、回転を 3 x 3 x N の配列にまとめます。この例の最後に示しているサポート関数 randomTransforms は、ランダムな値を用いて N 個の変換を初期化し、構造体を出力として提供します。

randomTransforms を使用してオブジェクト変換のマップとロボットの開始位置を設定します。

Map = randomTransforms(numObjects,roomDimensions);
Robot = randomTransforms(1,roomDimensions);

方程式の定義

マップ上に存在する対象物を正確に識別するために、ロボットがマップを変換してセンサーを起点に配置する必要があります。こうすることによってロボットは、現在検知しているものとこれから検知するものとを比較して、マップオブジェクトを見つけられるようになります。

マップオブジェクト $i$ のグローバルマップでの位置情報を変換することで、ロボットを基準とするそのマップオブジェクトの相対位置 $T_{r e l} (i)$ と向き $R_{r e l} (i)$ を求めることができます。

$\begin{array}{l} R_{r e l} (i) = R_{b o t}^{⊤} R_{m a p} (i) \\ T_{r e l} (i) = R_{b o t}^{⊤} (T_{m a p} (i) - T_{b o t}) \end{array}$

ここで、 $T_{b o t}$ と $R_{b o t}$ はロボットの位置と向きであり、 $T_{m a p} (i)$ と $R_{m a p} (i)$ はマップデータを表します。これに相当する MATLAB コードは次のようになります。

Rrel(:,:,i) = Rbot' * Rmap(:,:,i)
Trel(:,i) = Rbot' * (Tmap(:,i) - Tbot)

`for` ループを使用して多くの行列変換を CPU 上で実行

マップオブジェクトごとに、ロボットに対する相対的な位置情報を変換する必要があります。この例の最後に示しているサポート関数 loopingTransform は、すべての変換を順番にループ処理します。zeros の 'like' 構文に注目してください。この構文により、次のセクションで同じコードを GPU で使用できます。

計算時間を測定するには、関数 timeit を使用します。この関数は loopingTransform を複数回呼び出して平均時間を取得します。これには引数をもたない関数が必要なため、@() 構文を使用して正しい形式の無名関数を作成します。

cpuTime = timeit(@()loopingTransform(Robot,Map,numObjects));
fprintf('It takes %3.4f seconds on the CPU to execute %d transforms.\n', ...
        cpuTime, numObjects);

It takes 0.0047 seconds on the CPU to execute 1000 transforms.

GPU での同じコードの試用

このコードを GPU 上で実行するには、単にデータを gpuArray にコピーするだけで済みます。MATLAB が GPU に格納されたデータを認識した場合、gpuArray がサポートされていれば、それを使用して任意のコードを実行します。

gMap.R = gpuArray(Map.R);
gMap.T = gpuArray(Map.T);
gRobot.R = gpuArray(Robot.R);
gRobot.T = gpuArray(Robot.T);

次に、gputimeit を呼び出します。これは GPU 計算を含むコードの timeit に相当します。時間を記録する前に、すべての GPU 演算が終了していることを確認します。

fprintf('Computing...\n');

Computing...

gpuTime = gputimeit(@()loopingTransform(gRobot,gMap,numObjects));

fprintf('It takes %3.4f seconds on the GPU to execute %d transforms.\n', ...
        gpuTime, numObjects);

It takes 0.2745 seconds on the GPU to execute 1000 transforms.

fprintf(['Unvectorized GPU code is %3.2f times slower ',...
    'than the CPU version.\n'], gpuTime/cpuTime);

Unvectorized GPU code is 58.19 times slower than the CPU version.

`pagefun` を使用したバッチ処理

GPU を使用した上記の方法では、すべての計算が独立しているにもかかわらず順番にしか実行されないため、非常に時間がかかりました。pagefun を使用すると、すべての計算を並列に実行することができます。この例の最後に示しているサポート関数 pagefunTransform は、for ループの代わりに pagefun を使用して関数 loopingTransform と同じ変換を適用します。

gpuPagefunTime = gputimeit(@()pagefunTransform(gRobot, gMap));
fprintf(['It takes %3.4f seconds on the GPU using pagefun ',...
    'to execute %d transforms.\n'], gpuPagefunTime, numObjects);

It takes 0.0003 seconds on the GPU using pagefun to execute 1000 transforms.

fprintf(['Vectorized GPU code is %3.2f times faster ',...
    'than the CPU version.\n'], cpuTime/gpuPagefunTime);

Vectorized GPU code is 17.88 times faster than the CPU version.

fprintf(['Vectorized GPU code is %3.2f times faster ',...
    'than the unvectorized GPU version.\n'], gpuTime/gpuPagefunTime);

Vectorized GPU code is 1040.61 times faster than the unvectorized GPU version.

説明

初めの計算は、回転の計算でした。これは行列乗算を含み、関数 mtimes (*) に変換することができました。これを、乗算する 2 組の回転と共に pagefun に渡します。

Rel.R = pagefun(@mtimes, Robot.R', Map.R);

Robot.R' は 3 行 3 列の行列であり、Map.R は 3 x 3 x N の配列です。関数 pagefun が、マップから得た個々の独立した行列をロボットの該当する回転に一致させ、必要な 3 x 3 x N の出力を提供します。

また平行移動計算は行列乗算も含みますが、行列乗算の通常のルールでは、これを変更せずにループ外に出すことができます。

Rel.T = Robot.R' * (Map.T - Robot.T);

より高度な GPU のベクトル化 - "位置特定できないロボット" 問題の解決

ロボットがマップ上の不明な場所にいる場合、グローバル検索アルゴリズムを使用してロボット自身の位置を特定することがあります。このアルゴリズムは上記の計算を実行し、またロボットのセンサーによって確認済みのオブジェクトと、その位置から確認できるであろうものとの間の有用な対応を検索することで考えられる位置情報を多数テストします。

今度は、複数のロボットと複数のオブジェクトがあります。N 個のオブジェクトと M 個のロボットは、N*M に変換されます。'ロボットの空間' と 'オブジェクトの空間' を区別するために、回転に対しては 4 次元、平行移動に対しては 3 次元を使用します。つまり、ロボットの回転は 3 x 3 x 1 x M、平行移動は 3 x 1 x M となります。

ロボットをランダムに配置し、検索を初期化します。優れた検索アルゴリズムは位相的な手がかりまたはその他の手がかりを使用して、より高度な検索を設定します。

numRobots = 10;
Robot = randomTransforms(numRobots,roomDimensions);
Robot.R = reshape(Robot.R, 3, 3, 1, []); % Spread along the 4th dimension
Robot.T = reshape(Robot.T, 3, 1, []); % Spread along the 3rd dimension
gRobot.R = gpuArray(Robot.R);
gRobot.T = gpuArray(Robot.T);

この例の最後で定義しているサポート関数 loopingTransform2 は、2 つの入れ子にされたループを使用してループ変換を実行することでロボットとオブジェクトをループ処理します。

cpuTime = timeit(@()loopingTransform2(Robot,Map,numObjects,numRobots));
fprintf('It takes %3.4f seconds on the CPU to execute %d transforms.\n', ...
        cpuTime, numObjects*numRobots);

It takes 0.0852 seconds on the CPU to execute 10000 transforms.

今回は GPU の時間測定に、tic と toc を使用します。これらを使用しないと、計算に時間がかかりすぎるためです。この目的にとっては、この方法でも十分に正確です。出力データの作成に関連するすべてのコストが確実に含まれるように、timeit と gputimeit が既定で行う場合と同様に、単一の出力変数で loopingTransform2 を呼び出します。

fprintf('Computing...\n');

Computing...

tic;
gRel = loopingTransform2(gRobot,gMap,numObjects,numRobots);
gpuTime = toc;

fprintf('It takes %3.4f seconds on the GPU to execute %d transforms.\n', ...
        gpuTime, numObjects*numRobots);

It takes 3.6775 seconds on the GPU to execute 10000 transforms.

fprintf(['Unvectorized GPU code is %3.2f times slower ',...
    'than the CPU version.\n'], gpuTime/cpuTime);

Unvectorized GPU code is 43.17 times slower than the CPU version.

既に述べたように、ループによる方法を GPU 上で実行した場合、計算を並列に行わないため速度は非常に遅くなります。

pagefun による新しい方法では、mtimes に加えて transpose 演算子も pagefun の呼び出しに組み込む必要があります。また、squeeze を使用して、転置されたロボットの向きから大きさが 1 の次元を削除し、拡張をロボットに適用して 3 次元にし、平行移動と一致させます。上記にもかかわらず、結果のコードはかなりコンパクトになります。この例の最後に示しているサポート関数 pagefunTransform2 は、入れ子にされた for ループの代わりに pagefun の呼び出しを 2 つ使用して関数 loopingTransform2 と同じ変換を適用します。

ここでも、pagefun が次元を適切に拡張します。したがって、3 x 3 x 1 x M の行列 Rt と 3 x 3 x N x 1 の行列 Map.R を乗算する場合には、3 x 3 x N x M の行列が出力されます。

gpuPagefunTime = gputimeit(@()pagefunTransform2(gRobot,gMap));
fprintf(['It takes %3.4f seconds on the GPU using pagefun ',...
    'to execute %d transforms.\n'], gpuPagefunTime, numObjects*numRobots);

It takes 0.0012 seconds on the GPU using pagefun to execute 10000 transforms.

fprintf(['Vectorized GPU code is %3.2f times faster ',...
    'than the CPU version.\n'], cpuTime/gpuPagefunTime);

Vectorized GPU code is 72.28 times faster than the CPU version.

fprintf(['Vectorized GPU code is %3.2f times faster ',...
    'than the unvectorized GPU version.\n'], gpuTime/gpuPagefunTime);

Vectorized GPU code is 3120.72 times faster than the unvectorized GPU version.

まとめ

関数 pagefun は、arrayfun がサポートする大部分のスカラー演算に加えて、多くの 2 次元演算をサポートします。これらの関数を一緒に使用することで、行列代数や配列操作などの幅広い計算をベクトル化できるようになります。そのためループの必要がなくなり、パフォーマンスが大幅に向上します。

いずれの場所で GPU データに対して小規模な計算をループで行っている場合でも、この方法によるバッチ実装への転換を検討するようお勧めします。このようにすることで、それまでパフォーマンスの改善が見られなかった場所で、GPU を最大限に利用してパフォーマンスを向上させる機会が得られることもあります。

サポート関数

ランダム変換関数

関数 randomTransforms は、複数のオブジェクトについて、指定された次元の部屋でオブジェクト変換のマップとロボットの開始位置を設定します。

function Tform = randomTransforms(N,roomDimensions)
Tform.T = zeros(3, N);
Tform.R = zeros(3, 3, N);

for i = 1:N
    Tform.T(:,i) = rand(3, 1) .* roomDimensions';
    % To get a random orientation, we can extract an orthonormal
    % basis for a random 3-by-3 matrix.
    Tform.R(:,:,i) = orth(rand(3, 3));
end

end

ループ変換関数

関数 loopingTransform は、変換を順番にループ処理することで、各オブジェクトの位置をロボットに対する相対的な位置に変換します。

function Rel = loopingTransform(Robot,Map,numObjects)
Rel.R = zeros(size(Map.R), 'like', Map.R); % Initialize memory
Rel.T = zeros(size(Map.T), 'like', Map.T); % Initialize memory

for i = 1:numObjects
    Rel.R(:,:,i) = Robot.R' * Map.R(:,:,i);
    Rel.T(:,i) = Robot.R' * (Map.T(:,i) - Robot.T);
end

end

`pagefun` 変換関数

関数 pagefunTransform は、関数 pagefun を使用して変換を適用することで、各オブジェクトの位置をロボットに対する相対的な位置に変換します。

function Rel = pagefunTransform(Robot,Map)
Rel.R = pagefun(@mtimes, Robot.R', Map.R);
Rel.T = Robot.R' * (Map.T - Robot.T);
end

入れ子にされたループ変換関数

関数 loopingTransform2 は、2 つの入れ子にされたループを使用してループ変換を実行することでロボットとオブジェクトをループ処理します。この変換で、各オブジェクトの位置を各ロボットに対する相対的な位置にマップします。

function Rel = loopingTransform2(Robot,Map,numObjects,numRobots)
Rel.R = zeros(3, 3, numObjects, numRobots, 'like', Map.R);
Rel.T = zeros(3, numObjects, numRobots, 'like', Map.T);

for i = 1:numObjects
    for j = 1:numRobots
        Rel.R(:,:,i,j) = Robot.R(:,:,1,j)' * Map.R(:,:,i);
        Rel.T(:,i,j) = ...
            Robot.R(:,:,1,j)' * (Map.T(:,i) - Robot.T(:,1,j));
    end
end

end

`pagefun` の呼び出しを 2 つ使用した変換関数

関数 pagefunTransform2 は、関数 pagefun の呼び出しを 2 つ使用することで、各オブジェクトの位置を各ロボットに対する相対的な位置にマップします。

function Rel = pagefunTransform2(Robot, Map)
Rt = pagefun(@transpose, Robot.R);
Rel.R = pagefun(@mtimes, Rt, Map.R);
Rel.T = pagefun(@mtimes, squeeze(Rt), ...
    (Map.T - Robot.T));
end

参考

pagefun | gpuArray | arrayfun | gputimeit

pagefun による GPU 上の小さな行列問題のパフォーマンス改善