问CPU和GPU中的SVD速度

EN

一般来说，SVD是一个很难并行化的例程。你可以查看，使用高端的特斯拉显卡，加速效果并不是很令人印象深刻。

您有一张GTX460卡-- 。该卡针对游戏(单精度计算)进行了优化，而不是HPC (双精度计算)。单精度/双精度吞吐比为12，因此该卡具有873 GFLOPS SP / 72 GFLOPS DP。检查。

因此，如果Md数组使用双精度元素，那么它的计算将相当慢。此外，当调用CPU例程时，很可能会利用所有CPU核心，从而降低在GPU上运行例程的可能增益。此外，在GPU运行中，您需要为将缓冲区传输到设备而花费时间。

根据Divakar的建议，您可以使用Md = single(Md)将数组转换为单精度，然后再次运行基准测试。您可以尝试使用更大的数据集大小，以查看是否发生了变化。我不期望在你的GPU上从这个例程中获得太多好处。

更新1:

在您发布结果后，我看到DP/SP时间比率为2。在CPU端，这是正常的，因为您可以在SSE值寄存器中容纳的double值少2倍。然而，GPU端的比率仅为2意味着GPU代码没有充分利用SM核心-因为理论比率是12。换句话说，与DP相比，I期望优化代码的SP性能要好得多。情况似乎并非如此。

正如VAndrei已经指出的，奇异值分解是一种难以并行化的算法。

您的主要问题是矩阵的大小。随着矩阵尺寸的增大，奇异值分解的性能迅速下降。所以你的主要目标应该是减小矩阵的大小。这可以使用高斯法方程(这基本上是在最小二乘意义上的超定线性系统的缩减)来完成。

这可以通过简单地将转置乘以矩阵来完成：

MhReduced = Mh' * Mh;

这会将你的矩阵减少到cols*cols的大小(如果cols是Mh的列数)。然后您只需调用[U,S,V] = svd(MhReduced);

注意:使用这种方法可能会产生符号相反的奇异向量(如果您正在比较这些方法，这一点很重要)。

如果您的matix状态良好，这应该可以正常工作。然而，在矩阵病态的情况下，这种方法可能无法产生有用的结果，而由于奇异值分解的鲁棒性，直接应用奇异值分解仍然可以产生有用的结果。

这应该会立即提高您的性能，至少在矩阵足够大的情况下是这样。另一个优点是你可以使用更大的矩阵。你可能根本不需要使用GPU (因为矩阵太大，复制到GPU的成本太高，或者在减少矩阵后，矩阵太小，GPU的加速比不够大)。

另请注意，如果您使用返回值，则会损失很大一部分性能。如果您只对SVD计算的性能感兴趣，请不要采用任何返回值。如果您只对“解决方案向量”感兴趣，只需获取V(并访问最后一列)：[~,~, V] = svd(Mh);。

编辑：

我看过你的示例代码，但我不确定它是什么，你在计算。我也意识到很难理解我对A'*A做了什么，所以我会详细解释一下。

给定一个具有A*x=b的线性系统，A表示m行的系数矩阵，n个变量，x表示解向量，b表示常量向量(两者都有m行)，一个解可以计算如下：

• if A is square (m=n)：x = A^-1 * b，
• if A is not square (m!=n, m > n)：如果A不是正方形，则执行以下操作：

A*x=b

A'* A*x= A‘*b

X= (A‘* A)^-1 * A'*b

A" = (A'*A)^-1 * A'通常被称为伪逆。然而，这种计算确实会对矩阵的条件数产生负面影响。这个问题的一个解决方案是使用奇异值分解(SVD)。如果USV = svd(A)表示奇异值分解的结果，则伪逆由VS"U'给出，其中S"是通过取S.的非零元素的逆而形成的。所以A" = VS"U'。

x = A"*b

然而，由于奇异值分解是相当昂贵的，特别是对于大矩阵。如果矩阵A条件良好，并且不一定需要非常精确的结果(我们谈论的是1e-13或1e-14)，则可以使用通过(A'*A)^-1 * A计算伪逆的更快的方法。

如果您的情况实际上是A*x=0，只需使用奇异值分解并从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

issue

首先，我已经使用以下代码在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。还可以比较SVD代码的CPU和GPU版本的不同行为。

下表报告了时序( 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

当svd例程仅用于计算奇异值或完整的4时，第一个SVD列表示CPU和GPU Matlab版本的SVD之间的比较。可以看出，GPU版本比GPU版本慢得多。在上面的一些答案中已经指出了动机:将SVD计算并行化存在固有的困难。

使用cuSOLVER?的

在这一点上，显而易见的问题是:我们能用cuSOLVER获得一些加速吗？实际上，我们可以使用mexFiles来使cuSOLVER例程在Matlab下运行。不幸的是，cuSOLVER的情况甚至更糟，因为它可以从上表的最后两列中推断出来。这些列报告分别使用cusolverDnSgesvd进行仅奇异值计算和全奇异值分解计算的Singular values calculation only with CUDA和Parallel implementation for multiple SVDs using CUDA处的码的定时。可以看出，如果考虑到cuSOLVER的cusolverDnSgesvd是单精度的，而Matlab是双精度的，那么它的性能甚至比Matlab更差。

这种行为的动机在cusolverDnCgesvd performance vs MKL得到了进一步的解释，cuSOLVER库的经理Joe Eaton说

我理解这里的困惑。我们确实为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

可以看出，时间略高于CPU的情况，但现在与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;
}