Я физик, пишу программу, которая включает генерацию нескольких (порядка нескольких миллиардов) случайных чисел, взятых из гауссовского распределения. Я пытаюсь использовать 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;
}
}
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 где эта самая техника обсуждается и допускается для:
Распределения иногда хранят значения из своего связанного источника
случайные числа через звонки их оператору (). Например, общий
Метод генерации нормально распределенных случайных чисел заключается в
получить два равномерно распределенных случайных числа и вычислить два
нормально распределенные случайные числа из них. Для сброса
кэш случайных чисел распределения в определенное состояние, каждый
Распределение имеет функцию-член сброса. Это должно быть вызвано на
распределение всякий раз, когда связанный с ним двигатель заменяется или восстанавливается.
Некоторые мысли поверх превосходного ответа HH
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;
}
Если скорость имеет первостепенное значение, я бы попытался вычислить все, что смогу за пределами основного цикла 10 ^ 8. Можно ли предварительно вычислить sqrt (vec1 [j]), чтобы сэкономить на вызове sqrt ()? Это возможно
иметь vec2 [j] / sqrt (vec1 [j]) как один вектор?
Если невозможно предварительно вычислить эти векторы, я бы попытался сэкономить на доступе к памяти. Хранение частей vec2 [j] и vec1 [j] вместе может помочь с извлечением одной строки кэша вместо двух. Так объявляю vector<pair<double,double>> vec12(dim_vec);
и использовать в выборке y+=l*vec12[j].first/sqrt(vec12[j].second)