Числовая функция для записи суммы в Python - PullRequest
Новая
       39

Числовая функция для записи суммы в Python

7 голосов
/ 20 сентября 2011

Учитывая log(a) и log(b), я хочу вычислить log(a+b) (численно устойчивым образом).

Я написал для этого небольшую функцию: 

def log_add(logA,logB):
    if logA == log(0):
        return logB
    if logA<logB:
        return log_add(logB,logA)
    return log( 1 + math.exp(logB-logA) ) + logA

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

Вам известна стандартная функция math или numpy для вычисления log(a+b) из log(a) и log(b)?

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

Заранее спасибо, численные методы, ниндзя!

Ответы [ 2 ]

9 голосов
/ 20 сентября 2011

Примечание: лучший ответ до сих пор - просто использовать numpy.logaddexp(logA,logB).

Почему именно вы сравниваете с log(0)? Это равно -numpy.inf, в этом случае вы получите log(1 + math.exp(-inf-logB) ) + logB, который сводится к logB. Этот вызов всегда выдаст предупреждающее сообщение, которое будет очень медленным.

Я мог бы придумать эту строчку. Однако вам нужно действительно измерить, чтобы увидеть, если это на самом деле быстрее. Он использует только одну «сложную» функцию вычисления вместо двух, которые вы используете, и рекурсия не происходит, if все еще там, но скрыто (и может быть оптимизировано) в fabs / maximum. 

def log_add(logA,logB):
    return numpy.logaddexp(0,-numpy.fabs(logB-logA)) + numpy.maximum(logA,logB)

редактирование:

Я быстро сделал timeit () со следующими результатами:

  1. Ваша оригинальная версия заняла около 120 с
  2. Моя версия заняла около 30 с
  3. Я удалил сравнение с журналом (0) из вашей версии, и оно уменьшилось до 20 с
  4. Я отредактировал свой код, чтобы сохранить logaddexp, но также работал с вашим рекурсивом, если и он уменьшился до 18.

Обновленный код, вы также можете переключать рекурсивный вызов с помощью встроенной обновленной формулы, но это мало изменило мои временные тесты: 

def log_add2(logA, logB):
    if logA < logB:
        return log_add2(logB, logA)
    return numpy.logaddexp(0,logB-logA)+logA

Редактировать 2:

Как отметил в комментариях pv, вы могли бы просто сделать numpy.logaddexp(logA, logB), что сводится к вычислению log(exp(logA)+exp(logB)), что, конечно, равно log(A+B). Я рассчитал его (на той же машине, что и выше), и он пошел еще дальше до 10 секунд. Итак, мы снизились до 1/12, неплохо;).

0 голосов
/ 20 сентября 2011 
def log_add(logA, logB): 
    return math.log(math.exp(logA) + math.exp(logB))

слишком медленно? или неточно?

...