c ++ 11 — несколько случайных чисел переполнение стека

Я физик, пишу программу, которая включает генерацию нескольких (порядка нескольких миллиардов) случайных чисел, взятых из гауссовского распределения. Я пытаюсь использовать C ++ 11. Генерация этих случайных чисел отделяется операцией, которая должна занимать очень мало времени. Больше всего меня беспокоит, может ли тот факт, что я генерирую так много случайных чисел с таким небольшим временным интервалом, потенциально привести к неоптимальной производительности. Я проверяю определенные статистические свойства, которые в значительной степени зависят от независимости случайности чисел, поэтому мой результат особенно чувствителен к этим вопросам. У меня такой вопрос, с типами чисел, которые я упоминаю ниже в коде (упрощенная версия моего фактического кода), я делаю что-то явно (или даже слегка) неправильно?

#include <random>

// Several other includes, etc.

int main () {

int dim_vec(400), nStats(1e8);
vector<double> vec1(dim_vec), vec2(dim_vec);

// Initialize the above vectors, which are order 1 numbers.

random_device rd;
mt19937 generator(rd());
double y(0.0);
double l(0.0);

for (int i(0);i<nStats;i++)
{
for (int j(0);j<dim_vec;j++)
{
normal_distribution<double> distribution(0.0,1/sqrt(vec1[j]));
l=distribution(generator);
y+=l*vec2[j];
}
cout << y << endl;
y=0.0;
}
}

2

Решение

normal_distribution разрешено иметь состояние. И в этом конкретном распределении обычно генерируют номера в паре с каждым другим вызовом, а при нечетных вызовах возвращают второе кэшированное число. Создавая новый дистрибутив при каждом вызове, вы выбрасываете этот кеш.

К счастью, вы можете «сформировать» один дистрибутив, вызвав его с помощью normal_distribution :: param_type’s:

 normal_distribution<double> distribution;
using P = normal_distribution<double>::param_type;
for (int i(0);i<nStats;i++)
{
for (int j(0);j<dim_vec;j++)
{
l=distribution(generator, P(0.0,1/sqrt(vec1[j])));
y+=l*vec2[j];
}
cout << y << endl;
y=0.0;
}

Я не знаком со всеми реализациями std::normal_distribution, Однако я написал один для Libc ++. Поэтому я могу с некоторой уверенностью сказать вам, что мое небольшое переписывание вашего кода окажет положительное влияние на производительность. Я не уверен, какое влияние это окажет на качество, за исключением того, что я знаю, что это не ухудшит его.

Обновить

относительно Северин ПаппадоНиже приведен комментарий о законности создания пар чисел за раз в распределении: см. N1452 где эта самая техника обсуждается и допускается для:

Распределения иногда хранят значения из своего связанного источника
случайные числа через звонки их оператору (). Например, общий
Метод генерации нормально распределенных случайных чисел заключается в
получить два равномерно распределенных случайных числа и вычислить два
нормально распределенные случайные числа из них. Для сброса
кэш случайных чисел распределения в определенное состояние, каждый
Распределение имеет функцию-член сброса. Это должно быть вызвано на
распределение всякий раз, когда связанный с ним двигатель заменяется или восстанавливается.

6

Другие решения

Некоторые мысли поверх превосходного ответа HH

  1. Нормальное распределение (mu, sigma) генерируется из нормального (0,1) смещения и масштаба:

N (мю, сигма) = мю + N (0,1) * сигма

если ваше среднее значение (mu) всегда равно нулю, вы можете упростить и ускорить (не добавляя 0.0) свой код, выполнив что-то вроде

normal_distribution<double> distribution;
for (int i(0);i<nStats;i++)
{
for (int j(0);j<dim_vec;j++)
{
l  = distribution(generator);
y += l*vec2[j]/sqrt(vec1[j]);
}
cout << y << endl;
y=0.0;
}
  1. Если скорость имеет первостепенное значение, я бы попытался вычислить все, что смогу за пределами основного цикла 10 ^ 8. Можно ли предварительно вычислить sqrt (vec1 [j]), чтобы сэкономить на вызове sqrt ()? Это возможно
    иметь vec2 [j] / sqrt (vec1 [j]) как один вектор?

  2. Если невозможно предварительно вычислить эти векторы, я бы попытался сэкономить на доступе к памяти. Хранение частей vec2 [j] и vec1 [j] вместе может помочь с извлечением одной строки кэша вместо двух. Так объявляю vector<pair<double,double>> vec12(dim_vec); и использовать в выборке y+=l*vec12[j].first/sqrt(vec12[j].second)

1

По вопросам рекламы [email protected]