Скорость SVD в CPU и GPU

Как уже говорил В.Андрей, SVD — алгоритм, который сложно распараллелить.

Ваша основная проблема заключается в размере вашей матрицы. Производительность СВД быстро падает с ростом размера матрицы. Поэтому вашей главной целью должно быть уменьшение размера матрицы. Это можно сделать с помощью нормальных уравнений Гаусса (которые в основном представляют собой редукцию переопределенной линейной системы в смысле метода наименьших квадратов).

Это можно сделать, просто умножив транспонирование на матрицу:

MhReduced = Mh' * Mh;

Это уменьшает вашу матрицу до размера cols*cols (если cols — это количество столбцов Mh). Тогда вы просто позвоните [U,S,V] = svd(MhReduced);

Примечание. Использование этого метода может привести к сингулярным векторам с противоположным знаком (это важно, если вы сравниваете эти методы).

Если ваш матикс хорошо кондиционирован, это должно работать без проблем. Однако в случае плохо обусловленной матрицы этот метод может не дать полезного результата, тогда как непосредственное применение SVD все же может дать полезный результат из-за надежности SVD.

Это должно значительно увеличить вашу производительность, по крайней мере, с достаточно большими матрицами. Еще одним преимуществом является то, что вы можете использовать гораздо большие матрицы. Вам, вероятно, вообще не придется использовать GPU (поскольку либо матрицы настолько велики, что копирование на GPU стоит слишком дорого, либо после сокращения матрица настолько мала, что ускорение GPU не будет достаточно большим).

Также обратите внимание, что большая часть производительности теряется, если вы используете возвращаемые значения. Если вас интересует только производительность вычисления SVD, не принимайте никаких возвращаемых значений. Если вас интересует только «вектор решения», просто получите V (и получите доступ к последнему столбцу): [~,~, V] = svd(Mh);.

РЕДАКТИРОВАТЬ:

Я посмотрел на ваш пример кода, но я не уверен, что это такое, вы рассчитываете. Также я понял, что довольно сложно понять, что я сделал с A'*A, поэтому объясню подробно.

Для линейной системы с A*x=b, где A обозначает матрицу коэффициентов с m строками и n столбцами, x - вектор решения, а b - постоянный вектор (оба с m строками), решение можно вычислить следующим образом:

• если A квадрат (m=n): x = A^-1 * b,

• если A не квадрат (m!=n, m > n):

  A * x = b

  A'* A * x = A' * b

  x = (A' * A)^-1 * A'*b

A" = (A'*A)^-1 * A' обычно называют псевдоинверсией. Однако этот расчет отрицательно влияет на число обусловленности матрицы. Решением этой проблемы является использование разложения по сингулярным числам (SVD). Если USV = svd(A) обозначает результаты SVD, псевдообратное значение задается как VS"U', при этом S" формируется путем взятия обратного значения ненулевых элементов S. Итак, A" = VS"U'.

x = A"*b

Однако поскольку СВД довольно затратна, особенно с большими матрицами. Если матрица A хорошо обусловлена ​​и очень точные результаты не требуются (мы говорим о 1e-13 или 1e-14), можно использовать гораздо более быстрый подход, вычисляя псевдоинверсию через (A'*A)^-1 * A.

Если ваш случай на самом деле A*x=0, просто используйте SVD и прочитайте последний вектор-столбец из V, это решение.

Если вы используете SVD не для решения линейной системы, а для результатов U и S (как предполагает ваш пример), я не уверен, что то, что я опубликовал, вам поможет.

Источники: 1, 2, 3

Вот пример кода для тестирования. Протестируйте его с большими матрицами, вы увидите, что использование (A'*A)^-1 * A' намного быстрее, чем альтернативы.

clear all

nbRows = 30000;
nbCols = 100;
% Matrix A
A = rand(nbRows,nbCols);

% Vector b
b = rand(nbRows,1);

% A*x=b

% Solve for x, using SVD
% [U,S,V]=svd(A,0);
% x= V*((U'*b)./diag(S))
tic
[U1,S1,V1]=svd(A,0);
x1= V1*((U1'*b)./diag(S1));
toc

tic
[U1,S1,V1]=svd(A,0);
x2 = V1*inv(S1)*U1'*b;
toc

% Solve for x, using manual pseudo-inverse
% A*x=b
% A'*A*x = A'*b
% x = (A'*A)^-1 * A'*b
tic
x3 = inv(A'*A) * A'*b;
toc

% Solve for x, let Matlab decide how (most likely SVD)
tic
x4 = A\b;
toc

Проблема

Прежде всего, я воспроизвел вашу проблему в Matlab2016b, используя следующий код:

clear all
close all
clc

Nrows = 2500;
Ncols = 2500;

NumTests = 10;

h_A = rand(Nrows, Ncols);
d_A = gpuArray.rand(Nrows, Ncols);

timingCPU = 0;
timingGPU = 0;

for k = 1 : NumTests
    % --- Host
    tic
    [h_U, h_S, h_V] = svd(h_A);
%     h_S = svd(h_A);
    timingCPU = timingCPU + toc;

    % --- Device
    tic
    [d_U, d_S, d_V] = svd(d_A);
%     d_S = svd(d_A);
    timingGPU = timingGPU + toc;
end

fprintf('Timing CPU = %f; Timing GPU = %f\n', timingCPU / NumTests, timingGPU / NumTests);

С помощью приведенного выше кода можно либо вычислить только сингулярные значения, либо вычислить полный SVD, включая сингулярные векторы. Также можно сравнить различное поведение CPU и GPU версий кода SVD.

Время указано в следующей таблице (время в s; Intel Core i7-6700K CPU @ 4.00GHz, 16288 MB, Max threads(8), GTX 960):

              Sing. values only | Full SVD         | Sing. val. only | Full
                                |                  |                 |
Matrix size   CPU      GPU      | CPU       GPU    |                 |
                                |                  |                 |
 200 x  200   0.0021    0.043   |  0.0051    0.024 |   0.098         |  0.15
1000 x 1000   0.0915    0.3     |  0.169     0.458 |   0.5           |  2.3
2500 x 2500   3.35      2.13    |  4.62      3.97  |   2.9           |  23
5000 x 5000   5.2      13.1     | 26.6      73.8   |  16.1           | 161

Первые столбцы 4 относятся к сравнению между версиями Matlab CPU и GPU подпрограммы svd, когда она используется для вычисления только сингулярных значений или полного SVD. Как видно, версия на GPU может быть значительно медленнее версии на GPU. Мотивация уже была указана в некоторых ответах выше: существует неотъемлемая сложность распараллеливания вычислений SVD.

Используете cuSOLVER?

На данный момент возникает очевидный вопрос: можем ли мы получить некоторое ускорение с помощью cuSOLVER? Действительно, мы могли бы использовать mexFiles, чтобы подпрограммы cuSOLVER работали под Matlab. К сожалению, ситуация с cuSOLVER еще хуже, о чем можно судить по последним двум столбцам вышеприведенной таблицы. Такие столбцы сообщают о времени кодов в Расчет сингулярных значений только с CUDA и Параллельная реализация для нескольких SVD с использованием CUDA с использованием cusolverDnSgesvd для расчета только сингулярных значений и полного расчета SVD соответственно. Как видно, cusolverDnSgesvd cuSOLVER работает даже хуже, чем Matlab, если учесть, что он работает с одинарной точностью, а Matlab с двойной точностью.

Мотивация такого поведения дополнительно объясняется в производительность cusolverDnCgesvd по сравнению с MKL, где говорит Джо Итон, менеджер библиотеки cuSOLVER

Я понимаю путаницу здесь. Мы обеспечиваем приличное ускорение для факторизаций LU, QR и LDL^t, то же самое мы хотели бы сказать и для SVD. Наша цель с cuSOLVER состоит в том, чтобы предоставить плотные и разреженные прямые решатели как часть инструментария CUDA в первый раз; мы должны с чего-то начать. Поскольку CULA больше не поддерживается, мы сочли необходимым передать некоторые функции разработчикам CUDA 7.0. Поскольку в настоящее время CUDA работает на более чем x86 хостах CPUs, cuSOLVER заполняет потребность там, где нет MKL. При этом мы можем добиться большего успеха с SVD, но придется подождать следующего релиза CUDA, поскольку приоритеты и сроки уже сжаты.

Использование других библиотек

На данный момент другими возможностями являются использование других библиотек, таких как

1. CULA;
2. MAGMA;
3. ArrayFire.

CULA не предлагается бесплатно, поэтому я не пробовал.

У меня были некоторые проблемы с установкой с зависимостями MAGMA, поэтому я не исследовал этот момент дальше (отказ от ответственности: я ожидаю, что со временем я смогу решить такие проблемы).

Затем я, наконец, остановился на использовании ArrayFire.

Используя ArrayFire, у меня было следующее время для полного вычисления SVD:

 200 x  200      0.036
1000 x 1000      0.2
2500 x 2500      4.5
5000 x 5000     29

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

Вот код ArrayFire:

#include <arrayfire.h>
#include <cstdio>
#include <cstdlib>
#include <fstream>

using namespace af;

int main(int argc, char *argv[])
{
    const int N = 1000;

    try {

        // --- Select a device and display arrayfire info
        int device = argc > 1 ? atoi(argv[1]) : 0;
        af::setDevice(device);
        af::info();

        array A = randu(N, N, f64);
        af::array U, S, Vt;

        // --- Warning up
        timer time_last = timer::start();
        af::svd(U, S, Vt, A);
        S.eval();
        af::sync();
        double elapsed = timer::stop(time_last);
        printf("elapsed time using start and stop = %g ms \n", 1000.*elapsed);

        time_last = timer::start();
        af::svd(U, S, Vt, A);
        S.eval();
        af::sync();
        elapsed = timer::stop(time_last);
        printf("elapsed time using start and stop = %g ms \n", 1000.*elapsed);

    }
    catch (af::exception& e) {

        fprintf(stderr, "%s\n", e.what());
        throw;
    }

    return 0;
}

Я пытался распараллелить SVD на своем ноутбуке с GTX 460 в течение более одного месяца, что также было частью моей дипломной работы, я провел так много экспериментов, что позже обнаружил, что MATLAB чрезвычайно быстр и, кстати, превосходит мой код. , я использовал одну сторону Якоби, и я еще не видел ни одной статьи, раскрывающей алгоритм быстрее, чем svd в MATLAB. На графическом процессоре затраты времени на копирование памяти могут быть очень высокими, если вы не используете элегантную модель. Я отсылаю вас к дополнительной информации о CUDA. Если вам нужна помощь, пожалуйста, свяжитесь со мной.