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

Ответ 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 может переполняться, вы можете определить его как double вместо long long и умножить на invN на каждом шаге, но вы можете потерять много точности с этим подход.

РЕДАКТИРОВАТЬ. Для больших входов мы можем разбить вычисление в нескольких ведрах:

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

Ответ 2

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

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

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

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

Ответ 3

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

Ответ 4

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

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

Ответ 5

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

log(abcde) = 5*log(K)

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

Ответ 6

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

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

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

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

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

Ответ 7

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

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