Реализация logsumexp в C? - PullRequest
Новая
       13

Реализация logsumexp в C?

11 голосов
/ 13 ноября 2010

Кто-нибудь знает о числовой библиотеке C с открытым исходным кодом, которая обеспечивает функцию logsumexp?

Функция logsumexp(a) вычисляет сумму журнала экспонент (e ^ {a_1} + ... e ^ {a_n}) компонентов массива a, избегая числового переполнения.

1 Ответ

10 голосов
/ 13 ноября 2010

Вот очень простая реализация с нуля (проверено, хотя бы минимально): 

double logsumexp(double nums[], size_t ct) {
  double max_exp = nums[0], sum = 0.0;
  size_t i;

  for (i = 1 ; i < ct ; i++)
    if (nums[i] > max_exp)
      max_exp = nums[i];

  for (i = 0; i < ct ; i++)
    sum += exp(nums[i] - max_exp);

  return log(sum) + max_exp;
}

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

Если вы хотите, чтобы он работал без сбоев при задании 0 аргументов, вам придется добавить регистр для этого:)

...