CPU和GPU中SVD的速度

Question

我正在测试Matlab R2014a中的svd，但似乎没有CPU与GPU的加速比。我使用了一张GTX 460显卡和一颗Core 2 duo E8500处理器。

以下是我的代码:

%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

此外，运行时间因不同运行而异，但 CPU 和 GPU 版本大致相同。为什么没有加速呢？

以下是一些测试：

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.

- mrgloom

---

你使用了哪些运行时？你的CPU型号是什么？在记录运行时间之前，你是否运行了几次基准测试代码以预热GPU？ - Divakar

此外，为了进行公正的基准测试，我认为您需要使用gather即gather(svd(Md))。 - Divakar

svd(Md) 正在计算手头矩阵的唯一奇异值，即 1000。当使用 [Ud,Sd,Vd]= svd(Md) 时，您仍然没有加速吗？ - Vitality

请注意，Matlab函数的CPU版本也经过了良好的优化，如果您拥有并行计算工具箱，则可能在幕后使用多个核心。因此，Matlab通常不是展示GPU加速的好软件。 - Tae-Sung Shin

Answer 1

通常，SVD是一种难以并行化的例程。您可以在这里检查到，即使使用高端的Tesla卡片，速度提升也不是很显著。

您有一张GTX460卡片 - Fermi架构。该卡片针对游戏优化（单精度计算），而不是HPC（双精度计算）。单精度/双精度吞吐量比为12。因此，该卡片具有873 GFLOPS SP / 72 GFLOPS DP。请在这里查看。

如果Md数组使用双精度元素，则对其进行的计算速度会相当慢。此外，在调用CPU例程时，有很高的机会使用所有CPU核心，从而降低在GPU上运行例程的可能收益。此外，在GPU运行中，您需要为传输缓冲区到设备付出时间。根据Divakar的建议，您可以使用Md = single（Md）将数组转换为单精度并再次运行基准测试。您可以尝试使用更大的数据集大小来查看是否有所改变。我不指望在您的GPU上获得太多收益。更新1：在您发布结果后，我看到DP / SP时间比为2。在CPU方面，这是正常的，因为您可以将2倍少的双精度值放入SSE寄存器中。然而，在GPU方面，仅有2的比率意味着gpu代码没有充分利用SM核心-因为理论比率为12。换句话说，我预计优化代码的SP性能比DP要好得多。看来情况并非如此。

- VAndrei

---

为了扩展有关单/双精度的声明，如果 OP 可以接受较小的精度损失，则可以建议在计时 GPU 代码之前使用 Md = single(Md)。 - Divakar

@Divakar 谢谢！我已经在帖子中添加了评论。 - VAndrei

Dongarra和他的合著者在《使用GPU加速数值密集线性代数计算》一文中概述了在GPU上加速SVD所面临的一些挑战。引用文章中的话：“有很多数学公式可以用来解决这些问题，但在所有情况下，由于算法的本质，设计一个高效的计算都是具有挑战性的。” - njuffa

Answer 2

正如VAndrei所说，SVD算法很难并行化。

你的主要问题是矩阵的大小。随着矩阵大小的增加，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

然而，由于奇异值分解(SVD)通常非常昂贵，尤其是对于大型矩阵。如果矩阵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

- Baiz

---

这个 MhReduced = Mh' * Mh 的技巧叫什么？它类似于协方差矩阵吗？ - mrgloom

我只知道它的名称叫“正规方程”。我对矩阵和线性方程只有很少的了解，但这可能会有所帮助：http://www.math.wsu.edu/faculty/genz/448/lessons/l401.pdf - Baiz

这是A'*A的一个效果（有时会发生），但并不重要，因为如果乘以-1，奇异向量仍然有效（现在找不到来源）。然而，我已经阅读了关于PCA分析的内容，并且有点困惑，特别是SVD的作用。无论如何，如果您想计算主成分，有两种经典方法：

1. eigenVectors = princomp(A);
2. [eigenVectors, eigenValues] = eig(cov(A)); 目前，当使用SVD时，我得到了不同的结果。也许我做错了什么。

- Baiz

http://math.stackexchange.com/questions/3869/what-is-the-intuitive-relationship-between-svd-and-pca - mrgloom

我已经阅读了几篇关于这个主题的文章，包括你发布的链接。请注意，那里使用的示例假定数据集具有零均值。当您使用 M = rand(n,m) 创建矩阵时，情况并非如此；您必须先对数据集进行去均值处理（请参见我的最后一条评论）。还要注意，如果列是变量，行是观测值，则协方差矩阵由 matCovar = (X'*X)/(nbRows-1) 给出。不确定为什么在您发布的链接中处理方式不同。 - Baiz

Answer 3

问题

首先，我使用以下代码在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代码的不同行为。

以下表格报告了时间（计时单位为秒; Intel Core i7-6700K CPU @ 4.00GHz, 16288 MB, 最大线程(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列是关于CPU和GPU Matlab版本的svd例程进行比较，当它用于计算奇异值或完整SVD时。可以看到，GPU版本可能比GPU版本慢得多。如一些答案中已经指出的，其动机在于并行化SVD计算具有固有的困难。

使用cuSOLVER？

在这一点上，显而易见的问题是：我们能否使用cuSOLVER来加速一些操作？确实，我们可以使用mexFiles使cuSOLVER例程在Matlab下运行。不幸的是，从上表的最后两列可以推断出，cuSOLVER的情况甚至更糟。这些列报告了使用cusolverDnSgesvd进行奇异值计算和完整SVD计算的时间，分别用于使用CUDA的仅奇异值计算和并行实现多个SVD的CUDA。正如可以看到的那样，cuSOLVER的cusolverDnSgesvd执行甚至比Matlab更差，如果考虑到它处理单精度，而Matlab处理双精度。

这种行为的动机在cusolverDnCgesvd性能与MKL比较中得到了进一步解释。在这篇文章中，cuSOLVER库的经理Joe Eaton表示。

我理解这里的困惑。我们为LU、QR和LDL^t因式分解提供了相当快的加速，这也是我们想要为SVD提供的。我们的cuSOLVER旨在首次提供作为CUDA工具包一部分的稠密和稀疏直接求解器；我们必须从某个地方开始。由于CULA不再受支持，我们觉得有必要在CUDA 7.0中将一些功能交到开发者手中。由于CUDA现在可以运行在更多的x86主机CPU上，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

如您所见，时间略微较长，但现在与CPU相当。

以下是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;
}

- Vitality

Answer 4

我曾经尝试在我的装备有GTX 460的笔记本电脑上并行化SVD超过一个月，这也是我本科论文的一部分。我做了很多实验后才发现MATLAB非常快，而且比我的代码表现更好。顺便说一下，我使用了单边Jacobi算法，但我还没有看到任何比MATLAB的svd更快的算法。在GPU上，如果您没有使用优雅的模型，内存复制的时间成本可能会非常高，我建议您阅读更多关于CUDA的内容。如果您需要任何帮助，请联系我。

- kzm11

---

这应该是一条评论而不是答案。 :) - pix