Эффективный способ вычисления среднего геометрического числа

Я думаю, что я нашел способ сделать это, это объединило две процедуры в вопросе, подобно идее Питера. Вот пример кода.

double geometric_mean(std::vector<double> const&data)
{
const double too_large = 1.e64;
const double too_small = 1.e-64;
double sum_log = 0.0;
double product = 1.0;
for(auto x:data) {
product *= x;
if(product > too_large || product < too_small) {
sum_log+= std::log(product);
product = 1;
}
}
return std::exp((sum_log + std::log(product))/data.size());
}

Плохая новость: это идет с веткой. Хорошая новость: предсказатель ветвления, вероятно, понимает это почти всегда правильно (ветвь должна срабатывать очень редко).

Ветви можно было бы избежать, используя идею Питера о постоянном количестве терминов в продукте. Проблема в том, что переполнение / недополнение может происходить в течение всего лишь нескольких терминов, в зависимости от значений.

Вы можете ускорить это путем умножения чисел, как в исходном решении, и преобразования только в логарифмы для каждого определенного числа умножений (в зависимости от размера ваших начальных чисел).

Другой подход, который дал бы лучшую точность и производительность, чем метод логарифма, состоял бы в том, чтобы компенсировать показатели вне диапазона на фиксированную величину, поддерживая точный логарифм аннулированного превышения. Вот так:

const int EXP = 64; // maximal/minimal exponent
const double BIG = pow(2, EXP); // overflow threshold
const double SMALL = pow(2, -EXP); // underflow threshold

double product = 1;
int excess = 0; // number of times BIG has been divided out of product

for(int i=0; i<n; i++)
{
product *= A[i];
while(product > BIG)
{
product *= SMALL;
excess++;
}
while(product < SMALL)
{
product *= BIG;
excess--;
}
}

double mean = pow(product, 1.0/n) * pow(BIG, double(excess)/n);

Все умножения на BIG а также SMALL точны, и нет звонков log (трансцендентная и, следовательно, особенно неточная функция).

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

log(abcde) = 5*log(K)

log(ab) + log(cde)  = 5*log(k)

Суммирование журналов для стабильного вычисления продуктов совершенно нормально и довольно эффективно (если этого недостаточно: есть способы получить векторизованные логарифмы с несколькими операциями SSE — есть также векторные операции Intel MKL).

Чтобы избежать переполнения, общепринятым методом является предварительное деление каждого числа на максимальную или минимальную запись величины (или суммирование логарифмических разностей до логарифмического максимума или логарифмического минимума). Вы также можете использовать сегменты, если числа сильно различаются (например, суммируйте журнал малых и больших чисел отдельно). Обратите внимание, что, как правило, ничего из этого не требуется, за исключением очень больших наборов, так как журнал double никогда не бывает огромным (между -700 и 700).

Кроме того, вы должны отслеживать знаки отдельно.

вычисления log x обычно содержит то же количество значащих цифр, что и xкроме случаев, когда x близко к 1: вы хотите использовать std::log1p если вам нужно вычислить prod(1 + x_n) с небольшим x_n,

Наконец, если у вас есть проблемы с ошибками округления при суммировании, вы можете использовать Суммирование Кахана или варианты.

Вместо использования логарифмов, которые очень дороги, вы можете напрямую масштабировать результаты по степеням двух.

double geometric_mean(std::vector<double> const&data) {
double huge = scalbn(1,512);
double tiny = scalbn(1,-512);
int scale = 0;
double product = 1.0;
for(auto x:data) {
if (x >= huge) {
x = scalbn(x, -512);
scale++;
} else if (x <= tiny) {
x = scalbn(x, 512);
scale--;
}
product *= x;
if (product >= huge) {
product = scalbn(product, -512);
scale++;
} else if (product <= tiny) {
product = scalbn(product, 512);
scale--;
}
}
return exp2((512.0*scale + log2(product)) / data.size());
}