使用 CUDA 的标准差

Question

我正在编写一个代码来查找 6 个向量的均值和标准差，每个向量有 8000 个元素。我想知道我是否可以使用 CUDA 来加快操作速度。我可以想到如何使用 CUDA 找到平均值，但我无法理解如何使用 CUDA 计算标准差。有人可以帮我吗？

Answer 1

Here is a Thrust example 一次计算多个汇总统计信息，包括平均值和标准差。偏差。

#include <thrust/device_vector.h>
#include <thrust/host_vector.h>
#include <thrust/transform_reduce.h>
#include <thrust/functional.h>
#include <thrust/extrema.h>
#include <cmath>
#include <limits>

// This example computes several statistical properties of a data
// series in a single reduction.  The algorithm is described in detail here:
// http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Parallel_algorithm
//
// Thanks to Joseph Rhoads for contributing this example


// structure used to accumulate the moments and other 
// statistical properties encountered so far.
template <typename T>
struct summary_stats_data
{
    T n;
    T min;
    T max;
    T mean;
    T M2;
    T M3;
    T M4;

    // initialize to the identity element
    void initialize()
    {
      n = mean = M2 = M3 = M4 = 0;
      min = std::numeric_limits<T>::max();
      max = std::numeric_limits<T>::min();
    }

    T variance()   { return M2 / (n - 1); }
    T variance_n() { return M2 / n; }
    T skewness()   { return std::sqrt(n) * M3 / std::pow(M2, (T) 1.5); }
    T kurtosis()   { return n * M4 / (M2 * M2); }
};

// stats_unary_op is a functor that takes in a value x and
// returns a variace_data whose mean value is initialized to x.
template <typename T>
struct summary_stats_unary_op
{
    __host__ __device__
    summary_stats_data<T> operator()(const T& x) const
    {
         summary_stats_data<T> result;
         result.n    = 1;
         result.min  = x;
         result.max  = x;
         result.mean = x;
         result.M2   = 0;
         result.M3   = 0;
         result.M4   = 0;

         return result;
    }
};

// summary_stats_binary_op is a functor that accepts two summary_stats_data 
// structs and returns a new summary_stats_data which are an
// approximation to the summary_stats for 
// all values that have been agregated so far
template <typename T>
struct summary_stats_binary_op 
    : public thrust::binary_function<const summary_stats_data<T>&, 
                                     const summary_stats_data<T>&,
                                           summary_stats_data<T> >
{
    __host__ __device__
    summary_stats_data<T> operator()(const summary_stats_data<T>& x, const summary_stats_data <T>& y) const
    {
        summary_stats_data<T> result;

        // precompute some common subexpressions
        T n  = x.n + y.n;
        T n2 = n  * n;
        T n3 = n2 * n;

        T delta  = y.mean - x.mean;
        T delta2 = delta  * delta;
        T delta3 = delta2 * delta;
        T delta4 = delta3 * delta;

        //Basic number of samples (n), min, and max
        result.n   = n;
        result.min = thrust::min(x.min, y.min);
        result.max = thrust::max(x.max, y.max);

        result.mean = x.mean + delta * y.n / n;

        result.M2  = x.M2 + y.M2;
        result.M2 += delta2 * x.n * y.n / n;

        result.M3  = x.M3 + y.M3;
        result.M3 += delta3 * x.n * y.n * (x.n - y.n) / n2; 
        result.M3 += (T) 3.0 * delta * (x.n * y.M2 - y.n * x.M2) / n;

        result.M4  = x.M4 + y.M4;
        result.M4 += delta4 * x.n * y.n * (x.n * x.n - x.n * y.n + y.n * y.n) / n3;
        result.M4 += (T) 6.0 * delta2 * (x.n * x.n * y.M2 + y.n * y.n * x.M2) / n2;
        result.M4 += (T) 4.0 * delta * (x.n * y.M3 - y.n * x.M3) / n;

        return result;
    }
};

template <typename Iterator>
void print_range(const std::string& name, Iterator first, Iterator last)
{
    typedef typename std::iterator_traits<Iterator>::value_type T;

    std::cout << name << ": ";
    thrust::copy(first, last, std::ostream_iterator<T>(std::cout, " "));  
    std::cout << "\n";
}


int main(void)
{
    typedef float T;

    // initialize host array
    T h_x[] = {4, 7, 13, 16};

    // transfer to device
    thrust::device_vector<T> d_x(h_x, h_x + sizeof(h_x) / sizeof(T));

    // setup arguments
    summary_stats_unary_op<T>  unary_op;
    summary_stats_binary_op<T> binary_op;
    summary_stats_data<T>      init;

    init.initialize();

    // compute summary statistics
    summary_stats_data<T> result = thrust::transform_reduce(d_x.begin(), d_x.end(), unary_op, init, binary_op);

    std::cout <<"******Summary Statistics Example*****"<<std::endl;
    print_range("The data", d_x.begin(), d_x.end());

    std::cout <<"Count              : "<< result.n << std::endl;
    std::cout <<"Minimum            : "<< result.min <<std::endl;
    std::cout <<"Maximum            : "<< result.max <<std::endl;
    std::cout <<"Mean               : "<< result.mean << std::endl;
    std::cout <<"Variance           : "<< result.variance() << std::endl;
    std::cout <<"Standard Deviation : "<< std::sqrt(result.variance_n()) << std::endl;
    std::cout <<"Skewness           : "<< result.skewness() << std::endl;
    std::cout <<"Kurtosis           : "<< result.kurtosis() << std::endl;

    return 0;
}

Answer 2

这超出了我的专业领域，但是有用于计算标准差的单次迭代算法，可以转换为减少。特别是，我正在考虑 Welford 算法，如 Knuth, TAOCP, vol. 中所述。 2. 一个缺点是每一步都需要除法，但这很可能与必要的内存访问很好地平衡。该算法的可用在线参考似乎是：

http://www.johndcook.com/standard_deviation.html

Answer 3

我已经在 CUDA 中为数据挖掘解决了这个问题。我没有使用任何库。但是，它给了我很好的结果。问题是求128*100万个样本的标准差和均值。这就是我所做的。

1. 我的设备有 16KB 的共享内存。而且，我正在使用花车。因此，共享内存可以容纳 4,000 个元素。我的设备每个块的最大线程数是 512。所以，我可以有 8 个块。如果我将 16KB 分成 8 个块，我将得到 2,000KB（即 1 个浮点数用于 1 个线程）。一般来说，这不会匹配。如果您有更好的设备，则需要重新计算。

2. 要找到标准差，每个块有 512 个元素。您可以使用单个线程找到平方（元素均值）。

3. 下一个挑战是添加它并找到这些元素的总和。尝试使用与求均值相同的方法。对于 512 个元素。将结果复制到全局内存。

4. 迭代。求结果的平方根。

PS：相应地进行计划，以使全局内存调用最少。均值和标准差经常从内存中调用数据。

Answer 4

一个迟到的答案，但我已经在我的代码中使用thrust::transform_reduce 解决了这个问题（在 GTX 1070 上使用 100k 浮点数进行了测试）：

#include <thrust/transform_reduce.h>
#include <thrust/device_vector.h>
#include <thrust/functional.h>

#include <functional>
#include <cmath>

/*
 * @struct varianceshifteop
 * @brief a unary function that shifts input data
 * by their mean and computes the squares of them
 */
struct varianceshifteop
    : std::unary_function<float, float>
{
    varianceshifteop(float m)
        : mean(m)
    { /* no-op */ }

    const float mean;

    __device__ float operator()(float data) const
    {
        return ::pow(data - mean, 2.0f);
    }
};

int main(int argc, char** argv)
{
    thrust::device_vector<float> data{ ... };

    // sum elements and divide by the number of elements
    float mean = thrust::reduce(
        data.cbegin(),
        data.cend(),
        0.0f,
        thrust::plus<float>()) / data.size();

    // shift elements by mean, square, and add them
    float variance = thrust::transform_reduce(
            data.cbegin(),
            data.cend(),
            varianceshifteop(mean),
            0.0f,
            thrust::plus<float>()) / (data.size() - 1);

    // standard dev is just a sqrt away
    float stdv = std::sqrtf(variance);

    return 0;
}

Answer 5

这是使用Thrust 求向量均值和方差的解决方案：

template <typename T> struct square {
    __host__ __device__ T operator()(const T& x) const {
        return x * x;
    }
};

template <typename T> void mean_and_var(T a, int n, double* p_mean, double* p_var) {
    double sum = thrust::reduce(a, &a[n], 0.0, thrust::plus<double>());
    double sum_square = thrust::transform_reduce(
        a,
        &a[n],
        square<double>(),
        0.0,
        thrust::plus<double>()
    );
    double mean = sum / n;
    *p_mean = mean;
    *p_var = (sum_square / n) - mean*mean;
}

这是一个独立源文件中的解决方案以及一些分析信息，将 CPU 与 GPU 进行比较（注意完全按照要求回答问题，有必要在 6 个不同的位置上调用 mean_and_var 6 次长度为 8000 的向量，但这应该会给你要点）：

#include "stdio.h"
#include <thrust/reduce.h>
#include <thrust/device_vector.h>

#define PROFILING_INIT                                              \
    cudaEvent_t start, stop;                                        \
    float elapsedTime;

#define PROFILING_START                                             \
    cudaEventCreate(&start);                                        \
    cudaEventCreate(&stop);                                         \
    cudaEventRecord(start, 0);

#define PROFILING_STOP                                              \
    cudaEventRecord(stop, 0);                                       \
    cudaEventSynchronize(stop);                                     \
    cudaEventElapsedTime(&elapsedTime, start, stop);                \
    printf("Time elapsed:  %.3g ms\n", elapsedTime);

#define N (6*8000)
// #define N (6*8000*10)
// #define N (6*8000*100)
double a[N];

template <typename T> struct square {
    __host__ __device__ T operator()(const T& x) const {
        return x * x;
    }
};

void mean_and_var_cpu(double* a, int n, double* p_mean, double* p_var) {
    double sum = 0, sum_square = 0, mean;
    for (int i = 0; i < n; i++) {
        sum += a[i];
        sum_square += (a[i] * a[i]);
    }
    mean = sum / n;
    *p_mean = mean;
    *p_var = (sum_square / n) - mean*mean;
}

template <typename T> void mean_and_var(T a, int n, double* p_mean, double* p_var) {
    double sum = thrust::reduce(a, &a[n], 0.0, thrust::plus<double>());
    double sum_square = thrust::transform_reduce(a, &a[n], square<double>(), 0.0, thrust::plus<double>());
    double mean = sum / n;
    *p_mean = mean;
    *p_var = (sum_square / n) - mean*mean;
}

int main() {
    for (int i = 0; i < N; i++) {
        a[i] = i;
    }

    double mean, var;

    PROFILING_INIT;

    printf("With thrust:\n");
    PROFILING_START;
    mean_and_var<double*>(a, N, &mean, &var);
    PROFILING_STOP;
    printf("Mean = %f, var = %f\n", mean, var);

    printf("With thrust, using device memory:\n");
    thrust::device_vector<double> a_dev(N);
    thrust::copy(a, &a[N], a_dev.begin());
    PROFILING_START;
    mean_and_var<thrust::device_ptr<double>>(&a_dev[0], N, &mean, &var);
    PROFILING_STOP;
    printf("Mean = %f, var = %f\n", mean, var);

    printf("On CPU:\n");
    PROFILING_START;
    mean_and_var_cpu(a, N, &mean, &var);
    PROFILING_STOP;
    printf("Mean = %f, var = %f\n", mean, var);

}

对代码质量（我主要是一名 C 程序员，为了 Thrust 的目的而尝试适应 C++）和缺少 cmets 表示歉意。

在所有情况下，答案都与使用 Python statistics 模块计算的均值和总体方差一致：

import statistics
x = 8000*6
print(statistics.mean(list(range(x))))
# print(statistics.variance(list(range(x))))
print(statistics.pvariance(list(range(x))))

以下是一些分析信息，所有这些信息都是在 Jetson Nano 开发板上执行的：

N	Time taken for Thrust (ms)	Time taken for Thrust using device memory (ms)	Time taken for naive CPU implementation (ms)
6*8000	11.2	2.38	1.11
6*8000*10	108	7.62	11.7
6*8000*100	1.12e+03	41.7	109

分析的结论：

• 如果您要使用 Thrust，请确保使用设备内存！
• 对于足够大的输入大小，GPU 比 CPU 快
• 对于较小的输入尺寸，取决于平台，使用 CPU 可能比使用 GPU 更快（如果有人能告诉我为什么/如何改进我的 Thrust 代码以使其在较小的输入尺寸下更快，那会很棒！）