如何在CUDA中以最大性能归一化矩阵列？

Question

如何在CUDA中有效地对矩阵列进行归一化？

我的矩阵是以列为主存储的，通常大小为2000x200。

该操作可以用以下matlab代码表示。

A = rand(2000,200);

A = exp(A);
A = A./repmat(sum(A,1), [size(A,1) 1]);

这些操作可以通过Thrust、cuBLAS和/或cuNPP有效地完成吗？

以下是包括4个核心的快速实现。

想知道是否可以将这些操作合并为1或2个核心以提高性能，尤其是由cublasDgemv()执行的列求和步骤。

#include <cuda.h>
#include <curand.h>
#include <cublas_v2.h>
#include <thrust/device_vector.h>
#include <thrust/device_ptr.h>
#include <thrust/transform.h>
#include <thrust/iterator/constant_iterator.h>
#include <math.h>

struct Exp
{
    __host__ __device__ void operator()(double& x)
    {
        x = exp(x);
    }
};

struct Inv
{
    __host__ __device__ void operator()(double& x)
    {
        x = (double) 1.0 / x;
    }
};

int main()
{
    cudaDeviceSetCacheConfig(cudaFuncCachePreferShared);
    cublasHandle_t hd;
    curandGenerator_t rng;
    cublasCreate(&hd);
    curandCreateGenerator(&rng, CURAND_RNG_PSEUDO_DEFAULT);

    const size_t m = 2000, n = 200;
    const double c1 = 1.0;
    const double c0 = 0.0;

    thrust::device_vector<double> A(m * n);
    thrust::device_vector<double> sum(1 * n);
    thrust::device_vector<double> one(m * n, 1.0);

    double* pA = thrust::raw_pointer_cast(&A[0]);
    double* pSum = thrust::raw_pointer_cast(&sum[0]);
    double* pOne = thrust::raw_pointer_cast(&one[0]);

    for (int i = 0; i < 100; i++)
    {
        curandGenerateUniformDouble(rng, pA, A.size());


        thrust::for_each(A.begin(), A.end(), Exp());

        cublasDgemv(hd, CUBLAS_OP_T, m, n,
                &c1, pA, m, pOne, 1, &c0, pSum, 1);

        thrust::for_each(sum.begin(), sum.end(), Inv());

        cublasDdgmm(hd, CUBLAS_SIDE_RIGHT, m, n, pA, m, pSum, 1, pA, m);
    }

    curandDestroyGenerator(rng);
    cublasDestroy(hd);

    return 0;
}

- kangshiyin

---

可以使用CUDA来有效地完成。展示一些你编写的CUDA代码，以实现你想要的功能。 - sgarizvi

代码已添加。寻求性能提升。 - kangshiyin

Answer 1

我使用CUDA 5.0在M2090上比较了3种方法的性能。

1. [173.179 us] 在问题中展示的cublas实现
2. [733.734 us] 使用thrust::reduce_by_key的纯Thrust实现，来自@talonmies
3. [1.508 ms] 使用thrust::inclusive_scan_by_key的纯Thrust实现

可以看出：

1. 在这种情况下，cublas具有最高的性能；
2. 无论是thrust::reduce_by_key还是thrust::inclusive_scan_by_key都会启动多个内核，从而导致额外的开销；
3. thrust::inclusive_scan_by_key相对于thrust::reduce_by_key会向DRAM写入更多的数据，这可能是内核时间较长的原因之一；
4. cublas和thrust方法之间的主要性能差异在于矩阵列求和。thrust较慢可能是因为thrust::reduce_by_key被设计用于对长度不同的段进行归约，但cublas_gemv()只能应用于固定长度的段（行/列）。

当矩阵A足够大以忽略内核启动开销时，cublas方法仍然表现最佳。在A_{20,000 x 2,000}上的分析结果如下所示。

将第一个for_each操作与cublasSgemv调用合并，如@talonmies所示，可能进一步提高性能，但我认为应使用手写内核而不是thrust::reduce_by_key。

以下是3种方法的代码。

#include <cuda.h>
#include <curand.h>
#include <cublas_v2.h>
#include <thrust/device_vector.h>
#include <thrust/device_ptr.h>
#include <thrust/transform.h>
#include <thrust/reduce.h>
#include <thrust/scan.h>
#include <thrust/iterator/counting_iterator.h>
#include <thrust/iterator/transform_iterator.h>
#include <thrust/iterator/discard_iterator.h>
#include <thrust/iterator/permutation_iterator.h>
#include <math.h>

struct Exp: public thrust::unary_function<double, double>
{
    __host__ __device__ double operator()(double x)
    {
        return exp(x);
    }
};

struct Inv: public thrust::unary_function<double, double>
{
    __host__ __device__ double operator()(double x)
    {
        return (double) 1.0 / x;
    }
};

template<typename T>
struct MulC: public thrust::unary_function<T, T>
{
    T C;
    __host__ __device__ MulC(T c) :
        C(c)
    {
    }
    __host__ __device__ T operator()(T x)
    {
        return x * C;
    }
};

template<typename T>
struct line2col: public thrust::unary_function<T, T>
{
    T C;
    __host__ __device__ line2col(T C) :
            C(C)
    {
    }

    __host__ __device__ T operator()(T i)
    {
        return i / C;
    }
};

int main()
{
    cudaDeviceSetCacheConfig(cudaFuncCachePreferShared);
    cublasHandle_t hd;
    curandGenerator_t rng;
    cublasCreate(&hd);
    curandCreateGenerator(&rng, CURAND_RNG_PSEUDO_DEFAULT);

    const size_t m = 2000, n = 200;
    const double c1 = 1.0;
    const double c0 = 0.0;

    thrust::device_vector<double> A(m * n);
    thrust::device_vector<double> B(m * n);
    thrust::device_vector<double> C(m * n);
    thrust::device_vector<double> sum1(1 * n);
    thrust::device_vector<double> sum2(1 * n);
    thrust::device_vector<double> one(m * n, 1);

    double* pA = thrust::raw_pointer_cast(&A[0]);
    double* pB = thrust::raw_pointer_cast(&B[0]);
    double* pSum1 = thrust::raw_pointer_cast(&sum1[0]);
    double* pSum2 = thrust::raw_pointer_cast(&sum2[0]);
    double* pOne = thrust::raw_pointer_cast(&one[0]);

    curandGenerateUniformDouble(rng, pA, A.size());

    const int count = 2;

    for (int i = 0; i < count; i++)
    {
        thrust::transform(A.begin(), A.end(), B.begin(), Exp());
        cublasDgemv(hd, CUBLAS_OP_T, m, n, &c1, pB, m, pOne, 1, &c0, pSum1, 1);
        thrust::transform(sum1.begin(), sum1.end(), sum1.begin(), Inv());
        cublasDdgmm(hd, CUBLAS_SIDE_RIGHT, m, n, pB, m, pSum2, 1, pB, m);
    }

    for (int i = 0; i < count; i++)
    {
        thrust::reduce_by_key(
                thrust::make_transform_iterator(thrust::make_counting_iterator(0), line2col<int>(m)),
                thrust::make_transform_iterator(thrust::make_counting_iterator(0), line2col<int>(m)) + A.size(),
                thrust::make_transform_iterator(A.begin(), Exp()),
                thrust::make_discard_iterator(),
                sum2.begin());
        thrust::transform(
                A.begin(), A.end(),
                thrust::make_permutation_iterator(
                        sum2.begin(),
                        thrust::make_transform_iterator(thrust::make_counting_iterator(0), line2col<int>(m))),
                C.begin(),
                thrust::divides<double>());
    }

    for (int i = 0; i < count; i++)
    {
        thrust::inclusive_scan_by_key(
                thrust::make_transform_iterator(thrust::make_counting_iterator(0), line2col<int>(m)),
                thrust::make_transform_iterator(thrust::make_counting_iterator(0), line2col<int>(m)) + A.size(),
                thrust::make_transform_iterator(A.begin(), Exp()),
                C.begin());
        thrust::copy(
                thrust::make_permutation_iterator(
                        C.begin() + m - 1,
                        thrust::make_transform_iterator(thrust::make_counting_iterator(0), MulC<int>(m))),
                thrust::make_permutation_iterator(
                        C.begin() + m - 1,
                        thrust::make_transform_iterator(thrust::make_counting_iterator(0), MulC<int>(m))) + n,
                sum2.begin());
        thrust::transform(
                A.begin(), A.end(),
                thrust::make_permutation_iterator(
                        sum2.begin(),
                        thrust::make_transform_iterator(thrust::make_counting_iterator(0), line2col<int>(m))),
                C.begin(),
                thrust::divides<double>());
    }

    curandDestroyGenerator(rng);
    cublasDestroy(hd);

    return 0;
}

- kangshiyin

Answer 2

您应该能够将第一个 for_each 操作与 cublasSgemv 调用融合为单个 reduce_by_key 调用。如果您定义/重新定义函数对象如下：

struct Accessor : public thrust::unary_function<int,int>
{
    int lda;
    __host__ __device__ Accessor(int _lda) : lda(_lda) {};
    __host__ __device__ int operator()(const int& idx)
    {
        return idx/lda;
    }
};

struct Exp : public thrust::unary_function<double,double>
{
    __host__ __device__ double operator()(const double& x)
    {
        return exp(x);
    }
};

struct Inv : public thrust::unary_function<double,double>
{
    __host__ __device__ double operator()(const double& x)
    {
        return double(1.0) / x;
    }
};

您可以将标准化输出计算为：

Accessor columns(m);
thrust::reduce_by_key(
        thrust::make_transform_iterator(thrust::make_counting_iterator(int(0)), columns),
        thrust::make_transform_iterator(thrust::make_counting_iterator(int(m*n)), columns),
        thrust::make_transform_iterator(A.begin(), Exp()),
        thrust::make_discard_iterator(),
        sum.begin());

thrust::for_each(sum.begin(), sum.end(), Inv());

cublasDdgmm(hd, CUBLAS_SIDE_RIGHT, m, n, pA, m, pSum, 1, pA, m);

[免责声明：所有代码都是在浏览器中编写的，未经测试，请自行承担风险]

除了减少核心调用次数外，使用高级迭代器还可以消除对大型单位矩阵的需求，这应该会减少内存占用量以及执行求和和指数运算所需的总内存交换次数。

- talonmies

---

迭代器确实很“花哨”。我比较了cublas和thrust的方法。虽然thrust::reduce_by_key可能需要更低的内存带宽，但与cublasDgemv相比仍然较慢。有什么想法吗？ - kangshiyin

我怀疑相对性能会在很大程度上取决于使用的GPU和类型。在使用32位类型的不同GPU上，您可能会发现缩减方法比纯CUBLAS实现更接近性能。Thrust开发人员已经承认，自他们在Thrust中进行当前实现以来，先进的缩减技术已经有所发展，但总体而言，树状缩减模式始终不如作为FMAD流表达的最佳方案高效，就像这种情况一样。 - talonmies

我还建议尝试使用thrust :: transform而不是thrust_for_each。在某些情况下（诚然，一段时间以前），我发现它比for_each快一点。但它可能不会改变性能太多。 - talonmies

Answer 3

您可以按照以下方式使用ArrayFire：

array A = randu(2000, 2000);
A = exp(A);
A /= tile(sum(A, 0), A.dims(0), 1);

你也可以用thrust来实现这个功能。但是如果你要处理矩阵（而不是普通向量），就必须使用for循环，这样效率就会降低。

免责声明：我是Accelereyes的开发人员，正在开发arrayfire。

编辑：根据要求，我正在生成新的基准测试。

编辑：由于这个基准测试，我们在代码中发现了exp的性能问题。我们正在审查和修复它。

- Pavan Yalamanchili

---

谢谢！令人印象深刻的是，代码可以像Matlab一样简单。您能否与我的代码进行性能比较？因为我手头没有ArrayFire库。 - kangshiyin

@EricShiyinKang 更新了结果。 - Pavan Yalamanchili

我认为你的基准代码存在问题，这导致了cublas/thrust方法的池计时结果。这是修改后的bench.cu。 - kangshiyin

@EricShiyinKang，你为什么要在循环外和循环内都生成随机数呢？同时，我意识到在计时器::stop之前没有使用设备同步，导致了对thrust和arrayfire的结果产生偏差。我正在重新修改代码。 - Pavan Yalamanchili

在CURAND参考手册的性能注释中提到，curandCreateGenerator()后第一次调用curandGenerateUniformDouble()需要额外的时间。 - kangshiyin