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

Question

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

double geometric_mean(std::vector<double> const&data) // failure 
{ 
    auto product = 1.0; 
    for(auto x:data) product *= x; 
    return std::pow(product,1.0/data.size()); 
} 

Однако, это вполне может потерпеть неудачу из сгущенного или переполнения в накопленном product (примечание: long double действительно не избежать этой проблемы). Таким образом, следующий вариант подвести вверх логарифмы:

double geometric_mean(std::vector<double> const&data) 
{ 
    auto sumlog = 0.0; 
    for(auto x:data) sum_log += std::log(x); 
    return std::exp(sum_log/data.size()); 
} 

Это работает, но требует std::log() для каждого элемента, который является потенциально медленно. Могу ли я избежать этого? Например, отслеживая (эквивалент) экспонента и мантиссу накопленного product отдельно?

Answer 1

«расколом показатель и мантисса» решение:

double geometric_mean(std::vector<double> const & data) 
{ 
    double m = 1.0; 
    long long ex = 0; 
    double invN = 1.0/data.size(); 

    for (double x : data) 
    { 
     int i; 
     double f1 = std::frexp(x,&i); 
     m*=f1; 
     ex+=i; 
    } 

    return std::pow(std::numeric_limits<double>::radix,ex * invN) * std::pow(m,invN); 
} 

Если вы обеспокоены тем, что ex может переполнить вы можете определить его как двойной вместо long long, и умножьте на invN на каждом шагу, но вы можете потерять много точности при таком подходе.

EDIT Для больших входов, мы можем разбить вычисления в несколько ведер:

double geometric_mean(std::vector<double> const & data) 
{ 
    long long ex = 0; 
    auto do_bucket = [&data,&ex](int first,int last) -> double 
    { 
     double ans = 1.0; 
     for (;first != last;++first) 
     { 
      int i; 
      ans *= std::frexp(data[first],&i); 
      ex+=i; 
     } 
     return ans; 
    }; 

    const int bucket_size = -std::log2(std::numeric_limits<double>::min()); 
    std::size_t buckets = data.size()/bucket_size; 

    double invN = 1.0/data.size(); 
    double m = 1.0; 

    for (std::size_t i = 0;i < buckets;++i) 
     m *= std::pow(do_bucket(i * bucket_size,(i+1) * bucket_size),invN); 

    m*= std::pow(do_bucket(buckets * bucket_size, data.size()),invN); 

    return std::pow(std::numeric_limits<double>::radix,ex * invN) * m; 
} 

Answer 2

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

Answer 3

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

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()); 
} 

К сожалению, это связано с веткой.

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

Answer 4

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

log(abcde) = 5*log(K) 

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

Answer 5

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

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 (трансцендентная, и поэтому особенно неточные, функции).

Answer 6

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

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

Answer 7

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

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

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

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

И наконец, если у вас возникли проблемы с округлением при суммировании, вы можете использовать Kahan summation или варианты.