Повышение эффективности малых матричных задач на графическом процессоре с помощью PAGEFUN

В этом примере показано, как использовать pagefun улучшить эффективность применения большого количества независимых вращений и перемещений к объектам в 3-D окружении. Это типично для области значений задач, связанных с большим количеством вычислений на небольших массивах.

Графические процессоры наиболее эффективны при выполнении вычислений на очень больших матрицах. В MATLAB ® это обычно достигается путем векторизации кода, чтобы максимизировать работу, проделанную в каждой инструкции. Когда у вас есть большой набор данных, но вычисления разделены на много небольших матричных операций, может быть трудно максимизировать эффективность, одновременно работая на многих сотнях ядер графический процессор.

The arrayfun и bsxfun функции позволяют выполнять скалярные операции параллельно на графическом процессоре. The pagefun функция добавляет возможность выполнения матричных операций в пакете аналогичным образом. The pagefun функция доступна в Parallel Computing Toolbox™ для использования с gpuArrays.

В этом примере робот перемещается по известной карте, содержащей большое количество функций, которые робот может идентифицировать с помощью своих датчиков. Робот находит себя в карте путем измерения относительного положения и ориентации этих объектов и сравнения их с местоположениями карт. Если предположить, что робот не полностью потерян, он может использовать любое различие между ними, чтобы исправить свое положение, например, с помощью фильтра Калмана. Мы сосредоточимся на первой части алгоритма.

Этот пример является функцией, так что вспомогательные функции могут быть вложены в него.

function paralleldemo_gpu_pagefun

Настройте карту

Создадим карту объектов с рандомизированными положениями и ориентациями в большом помещении.

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

Мы представляем положения и ориентации, используя векторы 3 на 1 T и матрицы вращения 3 на 3 R. Когда у нас есть N этих преобразований, мы упаковываем переводы в 3-by-N матрицу и повороты в 3-by-3-by-N массив. Следующая функция инициализирует N преобразований со случайными значениями, обеспечивая структуру как выход:

    function Tform = randomTransforms(N)
        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

Теперь используйте это, чтобы настроить карту преобразований объектов и начальное место для робота.

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

Задайте уравнения

Чтобы правильно идентифицировать характеристики карты, робот должен преобразовать карту, чтобы поместить ее датчики в источник. Затем он может найти объекты карты, сравнивая то, что он видит с тем, что он рассчитывает видеть.

Для объекта карты$i$ мы можем найти его положение относительно робота${T_\mathrm{rel}}(i)$ и ориентации путем${\mathbf{R}_\mathrm{rel}}(i)$ преобразования его местоположения глобальной карты:

$$\begin{array}{l}
 \mathbf{R}_\mathrm{rel}(i) \;=\;
 \mathbf{R}_\mathrm{bot}^\top \mathbf{R}_\mathrm{map}(i) \\
 T_\mathrm{rel}(i) \;=\;
 \mathbf{R}_\mathrm{bot}^\top ({T_\mathrm{map}}(i) - T_\mathrm{bot})
\end{array}$$

где$T_\mathrm{bot}$ и$\mathbf{R}_\mathrm{bot}$ являются положением и ориентацией робота, и${T_\mathrm{map}}(i)$${\mathbf{R}_\mathrm{map}}(i)$ представляют данные карты. Эквивалентный код MATLAB выглядит следующим образом:

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

Выполните много матричных преобразований на центральном процессоре с помощью for цикл

Мы должны преобразовать каждый объект карты в его местоположение относительно робота. Мы можем сделать это последовательно, закольцовывая все преобразования в свою очередь. Обратите внимание на синтаксис 'like' для zeros который позволит использовать тот же код на графическом процессоре в следующем разделе.

    function Rel = loopingTransform(Robot, Map)
        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

К времени вычисления мы используем timeit функция, которая будет вызывать loopingTransform несколько раз, чтобы получить среднее значение времени. Поскольку это требует функции без аргументов, мы используем @() синтаксис для создания анонимной функции правой формы.

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

Попробуйте тот же код на графическом процессоре

Запустить этот код на графическом процессоре - это просто вопрос копирования данных в gpuArray. Когда MATLAB встречается с данными, хранящимися на графическом процессоре, он будет запускать любой код, использующий его на графическом процессоре, пока он поддерживается.

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

Теперь мы звоним gputimeit, что эквивалентно timeit для кода, который включает в себя расчет графический процессор. Перед записью времени убедитесь, что все операции GPU завершены.

fprintf('Computing...\n');
gpuTime = gputimeit(@()loopingTransform(gRobot, gMap));

fprintf('It takes %3.4f seconds on the GPU to execute %d transforms.\n', ...
        gpuTime, numObjects);
fprintf(['Unvectorized GPU code is %3.2f times slower ',...
    'than the CPU version.\n'], gpuTime/cpuTime);
Computing...
It takes 0.5588 seconds on the GPU to execute 1000 transforms.
Unvectorized GPU code is 53.90 times slower than the CPU version.

Пакетная обработка с использованием pagefun

Версия графический процессор выше была очень медленной, потому что, хотя все вычисления были независимыми, они выполнялись последовательно. Использование pagefun мы можем запустить все расчеты параллельно. Мы также нанимаем bsxfun вычислить переводы, поскольку они являются поэлементными операциями.

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

gpuPagefunTime = gputimeit(@()pagefunTransform(gRobot, gMap));
fprintf(['It takes %3.4f seconds on the GPU using pagefun ',...
    'to execute %d transforms.\n'], gpuPagefunTime, numObjects);
fprintf(['Vectorized GPU code is %3.2f times faster ',...
    'than the CPU version.\n'], cpuTime/gpuPagefunTime);
fprintf(['Vectorized GPU code is %3.2f times faster ',...
    'than the unvectorized GPU version.\n'], gpuTime/gpuPagefunTime);
It takes 0.0008 seconds on the GPU using pagefun to execute 1000 transforms.
Vectorized GPU code is 13.55 times faster than the CPU version.
Vectorized GPU code is 730.18 times faster than the unvectorized GPU version.

Объяснение

Первым вычислением стал расчет оборотов. Это включало умножение матрицы, которое переводится в функцию mtimes (*). Передадим это в pagefun наряду с двумя наборами вращений, которые нужно умножить:

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

Robot.R' является матрицей 3 на 3 и Map.R является массивом 3-by-3-by-N. The pagefun функция соответствует каждой независимой матрице от карты до того же вращения робота и дает нам необходимую 3-by-3-by-N выход.

Вычисление преобразования также включает умножение матрицы, но нормальные правила умножения матрицы позволяют этому выйти за пределы цикла без каких-либо изменений. Однако это также включает вычитание Robot.T от Map.T, которые являются различными размерами. Поскольку это операция элемента за элементом, мы можем использовать bsxfun чтобы соответствовать размерностям так же, как и pagefun сделал для вращений:

Rel.T = Robot.R' * bsxfun(@minus, Map.T, Robot.T);

На этот раз нам нужно было использовать поэлементный оператор, который преобразуется в функцию minus (-).

Более совершенная векторизация графический процессор - решение задачи «потерянного робота»

Если бы наш робот находился в неизвестной части карты, он мог бы использовать глобальный алгоритм поиска, чтобы найти себя. Алгоритм будет тестировать несколько возможных местоположений путем выполнения вышеописанных расчетов и поиска хорошего соответствия между объектами, наблюдаемыми датчиками робота, и тем, что он ожидал бы видеть в этом положении.

Теперь у нас есть несколько роботов, а также несколько объектов. N объектов и M роботов должны выдать нам N * M преобразований. Чтобы отличить 'пространство роботов' от 'пространство объектов', мы используем 4-ую размерность для вращений и 3-е для перемещений. Это означает, что вращение нашего робота будет 3-by-3-by-1-by-M, и перемещения будут 3-by-1-by-M.

Мы инициализируем наш поиск со случайными местоположениями роботов. Хороший алгоритм поиска использовал бы топологические или другие подсказки, чтобы заполнить поиск более интеллектуально.

numRobots = 10;
Robot = randomTransforms(numRobots);
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);

Наша новая функция циклического преобразования требует двух вложенных циклов, чтобы закольцовывать роботов, а также объекты.

    function Rel = loopingTransform2(Robot, Map)
        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

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

Для наших GPU timings мы используем tic и toc на этот раз, потому что в противном случае вычисление займет слишком много времени. Это будет достаточно точно для наших целей. Чтобы убедиться, что любые затраты, связанные с созданием выходных данных, включены, мы вызываем loopingTransform2 с одной выходной переменной, так же как timeit и gputimeit делать по умолчанию.

fprintf('Computing...\n');
tic;
gRel = loopingTransform2(gRobot, gMap); %#ok<NASGU> Suppress unused variable warning
gpuTime = toc;

fprintf('It takes %3.4f seconds on the GPU to execute %d transforms.\n', ...
        gpuTime, numObjects*numRobots);
fprintf(['Unvectorized GPU code is %3.2f times slower ',...
    'than the CPU version.\n'], gpuTime/cpuTime);
Computing...
It takes 7.0564 seconds on the GPU to execute 10000 transforms.
Unvectorized GPU code is 47.26 times slower than the CPU version.

Как и прежде, циклическая версия работает на графическом процессоре намного медленнее, потому что она не делает вычислений параллельно.

Новая pagefun версия должна включать в себя transpose оператор, а также mtimes в вызов на pagefun. Нам также нужно squeeze транспонированные ориентации робота, чтобы поместить разброс по роботам в 3-ью размерность, чтобы соответствовать перемещениям. Несмотря на это, получившийся код значительно компактнее.

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

Еще раз, pagefun и bsxfun расширение размерностей соответствующим образом. Так где мы умножаем 3-by-3-by-1-by-M матрицу Rt с 3-by-3-by-N-by-1 матрицей Map.RМы получаем 3-by-3-by-N-by-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);
fprintf(['Vectorized GPU code is %3.2f times faster ',...
    'than the CPU version.\n'], cpuTime/gpuPagefunTime);
fprintf(['Vectorized GPU code is %3.2f times faster ',...
    'than the unvectorized GPU version.\n'], gpuTime/gpuPagefunTime);
It takes 0.0025 seconds on the GPU using pagefun to execute 10000 transforms.
Vectorized GPU code is 59.97 times faster than the CPU version.
Vectorized GPU code is 2834.45 times faster than the unvectorized GPU version.

Заключение

The pagefun функция поддерживает ряд 2-D операций, а также большинство скалярных операций, поддерживаемых arrayfun и bsxfun. Вместе эти функции позволяют векторизировать область значений расчетов с помощью матричной алгебры и манипуляции с массивами, устраняя необходимость в циклах и делая огромные увеличения эффективности.

Всякий раз, когда вы делаете небольшие вычисления на данных графический процессор в цикле, вы должны рассмотреть преобразование в пакетную реализацию таким образом. Это также может быть возможностью использовать графический процессор для повышения эффективности в тех случаях, когда ранее он не повышал эффективность.

end