英文:
Vectorising logsum using numpy
问题
我需要计算np.sum(np.exp(L))，其中L是一个很长的数字列表，可能高达3000。np.float64会在exp(700)左右溢出。然而，np.logaddexp不受溢出影响，但仅适用于两个元素。
我编写了以下递归函数来计算一个列表的logsum。
import numpy as np
from smartprint import smartprint as sprint
def logsum(l):
if len(l) == 1:
return l[0]
return logsum([np.logaddexp(l[0], l[1])] + l[2:])
# 奇数长度
l = [-10, 1, 2, 45, 100]
assert (logsum(l) == np.log(np.sum(np.exp(l))))
sprint (logsum(l))
# 偶数长度
l = [-10, 1, 2, 43]
assert (logsum(l) == np.log(np.sum(np.exp(l))))
sprint (logsum(l))
输出：
logsum(l) : 100.0
logsum(l) : 43.0
另外，
logsum([1,2,3,4,50000]) 输出预期的 50000，而np.exp会溢出。现在，问题是我的列表L很大，有多达1000万个元素，我需要至少计算1000次logsum。因此，我想知道是否有一种方法可以使用np.logaddexp向量化，以便我可以用它来处理更长的列表，而不仅仅是两个元素。
英文:
I need to compute the np.sum(np.exp(L)) where L is a long list of numbers which can be as high as 3000. np.float64 will overflow around exp(700). However, np.logaddexp is immune to overflow but only works for two elements.
I wrote the following recursive function to compute it for a list.
import numpy as np
from smartprint import smartprint as sprint
def logsum(l):
if len(l) == 1:
return l[0]
return logsum([np.logaddexp(l[0], l[1])] + l[2:])
# odd length
l = [-10, 1, 2, 45, 100]
assert (logsum(l) == np.log(np.sum(np.exp(l))))
sprint (logsum(l))
# even length
l = [-10, 1, 2, 43]
assert (logsum(l) == np.log(np.sum(np.exp(l))))
sprint (logsum(l))
Output:
logsum(l) : 100.0
logsum(l) : 43.0
Additionally,
logsum([1,2,3,4,50000]) outputs 50000 as expected whereas the np.exp would overflow. Now, the problem is that my list L is huge and has upto 10 million elements and I need to compute the logsum at least 1000 times. So, I wonder if there is a way to vectorize using the np.logaddexp so that I can use it for a longer list instead of just two elements.
答案1
得分: 1
根据用户slothrop的评论，scipy中的logsumexp函数运行良好，已按预期进行了向量化。我之前的递归解决方案不再需要。在此处发布以保持完整性。
from scipy.special import logsumexp
logsumexp(l)
英文:
Based on user: slothrop's comment, the scipy function for logsumexp works well and is already vectorised as expected. My recursive solution is no longer needed. Posting here for completeness.
from scipy.special import logsumexp
logsumexp(l)
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论