Стандартное отклонение с использованием CUDA
Я пишу код, чтобы найти среднее и стандартное отклонение 6 векторов с 8000 элементов в каждом. Мне было интересно, смогу ли я сделать это с помощью CUDA и ускорить операцию. Я мог бы подумать, как найти среднее с помощью CUDA, но я не могу понять, как рассчитать стандартное отклонение с помощью CUDA. Может кто-нибудь помочь мне здесь, пожалуйста?
5 ответов
Я решил эту проблему в CUDA для интеллектуального анализа данных. Я не использовал никаких библиотек. Но это дало мне хорошие результаты. Задача состоит в том, чтобы найти стандартное отклонение и среднее значение 128*1 миллиона образцов. Это то, что я сделал.
Мое устройство имеет общую память 16 КБ. И я использую поплавки. Таким образом, общая память может вместить 4000 элементов. Максимальный поток на блок для моего устройства составляет 512. Таким образом, я могу иметь 8 блоков. Если я разделю 16KB на 8 блоков, я получу 2000KB (то есть 1 float для 1 потока). Как правило, это не совпадает. Если у вас есть лучшее устройство, вам нужно сделать это снова.
Чтобы найти стандартное отклонение, каждый блок имеет 512 элементов. Вы можете найти квадрат (элемент-среднее), используя одну нить.
Следующая задача - добавить это и найти сумму этих элементов. Попробуйте тот же метод, который вы использовали для нахождения среднего значения. Для 512 элементов. Скопируйте результат в глобальную память.
Итерация. Найдите квадратный корень из результата.
PS: Планируйте соответственно так, чтобы глобальные вызовы памяти были минимальными. Среднее и стандартное отклонение часто вызывает данные из памяти.
Вот пример Thrust, который вычисляет количество сводных статистических данных за один проход, включая среднее и стандартное. отклонение.
#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;
}
Это выходит за рамки моей компетенции, но существуют однопроходные итерационные алгоритмы для вычисления стандартного отклонения, которые могут быть преобразованы в сокращение. В частности, я имею в виду алгоритм Уэлфорда, описанный в Knuth, TAOCP, vol. 2. Один недостаток состоит в том, что для каждого шага требуется разделение, но это, вероятно, будет хорошо сбалансировано с необходимыми обращениями к памяти. Полезная онлайн-ссылка на алгоритм выглядит так:
Вот решение для нахождения среднего значения и дисперсии вектора с использованием 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;
}
Вот это решение в автономном исходном файле вместе с некоторой информацией о профилировании, сравнивающей ЦП и ГП (обратите внимание, чтобы ответить на вопрос точно так, как задано, необходимо позвонить в
mean_and_var6 раз на 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, пытающийся приспособиться к C++ для целей Thrust) и отсутствие комментариев.
Во всех случаях ответы согласуются со средним значением и дисперсией населения, рассчитанными с использованием 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:
Выводы из профилирования:
- Если вы собираетесь использовать Thrust, убедитесь, что вы используете память устройства!
- Для достаточно большого размера ввода GPU работает быстрее, чем CPU.
- Для меньших размеров ввода, в зависимости от платформы, может быть быстрее использование ЦП, чем GPU (если кто-нибудь может сказать мне, почему и как улучшить мой код Thrust, чтобы сделать его быстрее для меньших размеров ввода, это было бы здорово! )
Поздний ответ, но я решил эту проблему, используя thrust::transform_reduce в моем коде (проверено с 100k плавает на GTX 1070):
#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;
}