matlab - Velocidad de SVD en CPU y GPU

matrix cuda

Question

Estoy probando svd en Matlab R2014a y parece que no hay aceleración CPU vs. GPU . Estoy usando una tarjeta GTX 460 y un Core 2 duo E8500 .

Aquí está mi código:

%test SVD
n=10000;
%host
Mh= rand(n,1000);
tic
%[Uh,Sh,Vh]= svd(Mh);
svd(Mh);
toc
%device
Md = gpuArray.rand(n,1000);
tic
%[Ud,Sd,Vd]= svd(Md);
svd(Md);
toc

Además, los tiempos de ejecución son diferentes de ejecución a ejecución, pero las versiones de CPU y GPU son más o menos iguales. ¿Por qué no hay aceleración?

Aquí hay algunas pruebas

for i=1:10
    clear;
    m= 10000;
    n= 100;
    %host
    Mh= rand(m,n);
    tic
    [Uh,Sh,Vh]= svd(Mh);
    toc
    %device
    Md = gpuArray.rand(m,n);
    tic
    [Ud,Sd,Vd]= svd(Md);
    toc
end

>> test_gpu_svd
Elapsed time is 43.124130 seconds.
Elapsed time is 43.842277 seconds.
Elapsed time is 42.993283 seconds.
Elapsed time is 44.293410 seconds.
Elapsed time is 42.924541 seconds.
Elapsed time is 43.730343 seconds.
Elapsed time is 43.125938 seconds.
Elapsed time is 43.645095 seconds.
Elapsed time is 43.492129 seconds.
Elapsed time is 43.459277 seconds.
Elapsed time is 43.327012 seconds.
Elapsed time is 44.040959 seconds.
Elapsed time is 43.242291 seconds.
Elapsed time is 43.390881 seconds.
Elapsed time is 43.275379 seconds.
Elapsed time is 43.408705 seconds.
Elapsed time is 43.320387 seconds.
Elapsed time is 44.232156 seconds.
Elapsed time is 42.984002 seconds.
Elapsed time is 43.702430 seconds.


for i=1:10
    clear;
    m= 10000;
    n= 100;
    %host
    Mh= rand(m,n,''single'');
    tic
    [Uh,Sh,Vh]= svd(Mh);
    toc
    %device
    Md = gpuArray.rand(m,n,''single'');
    tic
    [Ud,Sd,Vd]= svd(Md);
    toc
end

>> test_gpu_svd
Elapsed time is 21.140301 seconds.
Elapsed time is 21.334361 seconds.
Elapsed time is 21.275991 seconds.
Elapsed time is 21.582602 seconds.
Elapsed time is 21.093408 seconds.
Elapsed time is 21.305413 seconds.
Elapsed time is 21.482931 seconds.
Elapsed time is 21.327842 seconds.
Elapsed time is 21.120969 seconds.
Elapsed time is 21.701752 seconds.
Elapsed time is 21.117268 seconds.
Elapsed time is 21.384318 seconds.
Elapsed time is 21.359225 seconds.
Elapsed time is 21.911570 seconds.
Elapsed time is 21.086259 seconds.
Elapsed time is 21.263040 seconds.
Elapsed time is 21.472175 seconds.
Elapsed time is 21.561370 seconds.
Elapsed time is 21.330314 seconds.
Elapsed time is 21.546260 seconds.

Answer 1

En general, la SVD es una rutina difícil de paralelizar. Puede comprobar aquí que con una tarjeta Tesla de gama alta, la aceleración no es muy impresionante.

Usted tiene una tarjeta GTX460 - arquitectura Fermi . La tarjeta está optimizada para juegos (cálculos de precisión simple), no HPC (computación de precisión doble). La relación de rendimiento de Precisión Simple / Doble Precisión es 12. Por lo tanto, la tarjeta tiene 873 GFLOPS SP / 72 GFLOPS DP . Mira aquí .

Entonces, si la matriz de Md usa elementos de doble precisión, entonces el cálculo en ella sería bastante lento. También hay una alta probabilidad de que al llamar a la rutina de la CPU, se utilicen todos los núcleos de la CPU, reduciendo la posible ganancia de ejecutar la rutina en la GPU. Además, en la ejecución de la GPU, usted paga tiempo por transferir el búfer al dispositivo .

Según la sugerencia de Divakar, podría usar Md = single(Md) para convertir su matriz a precisión simple y ejecutar nuevamente el benchmark. Puede probar e ir con un tamaño de conjunto de datos mayor para ver si algo cambia. No espero ganar demasiado para esta rutina en su GPU.

Actualización 1:

Después de publicar los resultados, vi que la relación de tiempo DP / SP es 2. En el lado de la CPU, esto es normal, porque puede ajustar dos veces menos valores double en los registros SSE. Sin embargo, una relación de solo 2 en el lado de la GPU significa que el código de la GPU no hace el mejor uso de los núcleos SM, porque la relación teórica es 12. En otras palabras, habría esperado un rendimiento de SP mucho mejor para un código optimizado, en comparación con DP . Parece que este no es el caso.

Answer 2

Como VAndrei ya ha declarado, la SVD es un algoritmo que es difícil de paralelizar.

Tu problema principal es el tamaño de tu matriz. El rendimiento de la SVD cae rápidamente con un tamaño de matriz creciente. Entonces su objetivo principal debe ser reducir el tamaño de la matriz. Esto se puede lograr utilizando ecuaciones normales de Gauss (que es básicamente una reducción de un sistema lineal sobredeterminado en el sentido de mínimos cuadrados).

Esto se puede hacer simplemente multiplicando la transposición en la matriz:

MhReduced = Mh'' * Mh;

Esto reduce su matriz al tamaño de cols * cols (si cols es el número de columnas de Mh). Luego solo llama a [U,S,V] = svd(MhReduced);

Nota: El uso de este método puede generar vectores singulares con signo opuesto (solo importante si está comparando estos métodos).

Si tu matix está bien acondicionado, debería funcionar sin problemas. Sin embargo, en el caso de una matriz mal acondicionada, este método puede no producir un resultado utilizable, mientras que la aplicación directa de SVD aún podría producir un resultado utilizable debido a la solidez de la SVD.

Esto debería aumentar tu rendimiento inmensamente, al menos con matrices lo suficientemente grandes. Otra ventaja es que puedes usar matrices mucho más grandes. Probablemente no tenga que usar la GPU en absoluto (ya que las matrices son tan grandes que la copia a la GPU cuesta demasiado o después de la reducción, la matriz es tan pequeña que la aceleración de la GPU no será lo suficientemente grande).

También tenga en cuenta que se pierde un gran trozo de rendimiento si usa valores devueltos. Si solo está interesado en el rendimiento de la calumnia SVD, no tome ningún valor de retorno. Si solo está interesado en el "vector de solución", simplemente obtenga V (y acceda a la última columna): [~,~, V] = svd(Mh); .

EDITAR:

He revisado tu código de muestra, pero no estoy seguro de qué se trata, estás calculando. También me di cuenta de que es bastante difícil entender lo que hice con A''*A , así que lo explicaré en detalle.

Dado un sistema lineal con A*x=b , A que denota la matriz de coeficientes con m filas y n cols, x el vector solución yb el vector constante (ambos con m filas), una solución puede calcularse de la siguiente manera:

• si A es cuadrado ( m=n ): x = A^-1 * b ,

• si A no es cuadrado ( m!=n, m > n ):

  A * x = b

  A ''* A * x = A'' * b

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

A" = (A''*A)^-1 * A'' se llama típicamente pseudoinverso. Sin embargo, este cálculo influye negativamente en el número de condición de la matriz. Una solución a este problema es usar una descomposición de valores singulares (SVD). Si USV = svd (A) denota los resultados de la SVD, el pseudoinverso es dado por VS"U'' , con S" se forma tomando el inverso de los elementos distintos de cero de S. Entonces A" = VS"U'' .

x = A"*b

Sin embargo, dado que una SVD es bastante costosa, especialmente con matrices grandes. Si la matriz A está bien condicionada y los resultados muy precicse no son necesariamente necesarios (estamos hablando de 1e-13 o 1e-14), el enfoque mucho más rápido al calcular la vía peseudo-inversa (A''*A)^-1 * A puede ser usado.

Si su caso en realidad es A*x=0 , simplemente use un SVD y lea el último vector de columna de V, es la solución.

Si usas la SVD para no resolver un sistema lineal, sino los resultados de U y S (como sugiere tu ejemplo), no estoy seguro de que lo que he publicado te ayude.

Fuentes: 1 , 2 , 3

Aquí hay un código de muestra para probar. Pruébelo con matrices grandes, verá que usar (A''*A)^-1 * A'' es mucho más rápido que las alternativas.

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

Answer 3

He tratado de paralelizar SVD en mi computadora portátil equipada con GTX 460 durante más de un mes, que también fue parte de mi tesis de pregrado, hice tantos experimentos que más tarde descubrí que MATLAB es extremadamente rápido y supera mi código, por cierto , Utilicé un lado, Jacobi, y todavía no he visto ningún documento que revele un algoritmo más rápido que el svd de MATLAB. En GPU, el costo de tiempo de la copia de memoria puede ser muy alto si no está usando un modelo elegante, lo remito para que lea más sobre CUDA. Si necesita ayuda, contácteme.

Answer 4

La cuestión

Antes que nada, he replicado su problema en Matlab2016b usando el siguiente código:

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);

Con el código anterior, es posible calcular solo los valores singulares o calcular el SVD completo incluyendo los vectores singulares. También es posible comparar el comportamiento diferente de las versiones de CPU y GPU del código SVD.

El tiempo se informa en la siguiente tabla (tiempo en s : Intel Core i7-6700K CPU @ 4.00GHz , 16288 MB , Max threads(8) Intel Core i7-6700K CPU @ 4.00GHz 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

Las primeras 4 columnas se refieren a una comparación entre las versiones CPU y GPU Matlab de la rutina svd cuando se usa para calcular solo los valores singulares o la SVD completa. Como se puede ver, la versión de la GPU puede ser significativamente más lenta que la GPU. La motivación ya ha sido señalada en algunas respuestas anteriores: existe una dificultad inherente para paralelizar el cálculo SVD.

¿Usando cuSOLVER?

En este punto, la pregunta obvia es: ¿podemos cuSOLVER poco con cuSOLVER ? De hecho, podríamos usar mexFiles para hacer que las rutinas cuSOLVER ejecuten bajo Matlab. Desafortunadamente, la situación con cuSOLVER es aún peor, como se puede deducir de las dos últimas columnas de la tabla anterior. Tales columnas informan el tiempo de los códigos en el cálculo de valores singulares solo con la implementación de CUDA y Paralelo para múltiples SVD utilizando CUDA utilizando cusolverDnSgesvd para el cálculo de valores únicos y el cálculo completo de SVD, respectivamente. Como se puede ver, cuSOLVER de cusolverDnSgesvd funciona incluso peor que Matlab, si se tiene en cuenta que se trata de precisión simple, mientras que Matlab con precisión doble.

La motivación para este comportamiento se explica con más detalle en el rendimiento de cusolverDnCgesvd frente a MKL, donde Joe Eaton, gerente de la biblioteca cuSOLVER , dice:

Entiendo la confusión aquí. Proporcionamos una aceleración decente para las factorizaciones LU , QR y LDL^t , que es lo que también nos gustaría decir para SVD . Nuestro propósito con cuSOLVER es proporcionar soluciones directas densas y dispersas como parte del conjunto de herramientas de CUDA por primera vez; tenemos que comenzar en algún lado Como CULA ya no es compatible, sentimos que era urgente poner alguna funcionalidad en manos de los desarrolladores en CUDA 7.0 . Dado que CUDA ejecuta en más CPUs host x86 estos días, cuSOLVER una necesidad donde no hay MKL . Habiendo dicho eso, podemos hacer mejor con SVD , pero tendrá que esperar a que la próxima publicación de CUDA , las prioridades y los plazos ya sean estrictos.

Usando otras bibliotecas

En este punto, otras posibilidades están usando otras bibliotecas como

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

CULA no se ofrece de forma gratuita, así que no lo he probado.

Tuve algunos problemas de instalación con las dependencias de MAGMA , por lo que no he investigado más este punto (descargo de responsabilidad: espero que, con más tiempo, pueda resolver estos problemas).

Finalmente, finalmente terminé usando ArrayFire .

Usando ArrayFire , tuve el siguiente tiempo para el cálculo completo de SVD :

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

Como se puede ver, el tiempo es ligeramente más alto, pero ahora comparable, al caso de la CPU.

Aquí está el código 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;
}