3

I want to use numpy's logsumexp() in python 2.7. The formula I need to solve looks like this:

log ( 1 + e^a1 + e^a2 + e^a3 + ... e^an - e^ax )

The last term which is a negative number just has to be appended on. Excluding this last term, I would do the following:

myarray = numpy.array([0, a1, a2, a3, ..., an])

That way, with the first element being 0, then e^0 = 1 and so I have my first term, which is 1. Then I would just use

  result = numpy.logsumexp(myarray)

and I would get the correct result.

But now I have to append a -e^ax, and because it's negative, I can't simply append ax to the end of myarray. I also can't append -ax because that's just wrong, it would mean that I'm adding 1/e^ax, instead of -e^ax.

Is there any direct way to append this so that I can still use logsumexp()? The only reason I'm insisting on using logsumexp() rather than separately using numpy.exp() and numpy.sum() and numpy.log() is because I have the impression that logsumexp also incorporates stability within it in order to prevent underflows (correct me if I'm wrong). However if there's no other way around then I guess I have no choice.

user961627
  • 12,379
  • 42
  • 136
  • 210

1 Answers1

5

According to scipy.misc.logsumexp documentation:

scipy.misc.logsumexp(a, axis=None, b=None)

Parameters:   
    b: array-like, optional
        Scaling factor for exp(a).
        Must be of the same shape as a or broadcastable to a.
    New in version 0.12.0.

So, you could add list of factors like this:

In [2]: a = [0, 1, 3, 2]

In [3]: logsumexp(a, b=[1] * (len(a) - 1) + [-1])
Out[3]: 2.7981810916785101
awesoon
  • 32,469
  • 11
  • 74
  • 99
  • Thanks - but I don't exactly understand how to use it for my case? I need to subtract e^ax from the list. I tried to follow what you wrote using a calculator but got a different answer, or I guess my python isn't good enough to follow. – user961627 Aug 30 '14 at 09:35
  • Hmm... Line #23 should do the trick. Could you, please, add an example, producing unexpected result? – awesoon Aug 30 '14 at 09:39
  • Alright, example: I need to take the log of the sum of `e^0 + e^1 + e^3 - e^2`. So the sum is about `16.41476`. Then the log of this is about `2.79818...`. So the final result should be `2.79818`. (Not rounded exactly). – user961627 Aug 30 '14 at 09:55
  • I've updated the post with your example, does it solve the problem? – awesoon Aug 30 '14 at 10:15
  • Thank you! Yes this works. Could you please explain though? I understand that `b` equals to 1 * the last element of `a`, but why the `+[-1]` ? – user961627 Aug 30 '14 at 11:20
  • No, expression `[1] * (len(a) - 1) + [-1]` generates list with all `1` except last, which is `-1`. Just input this in python. – awesoon Aug 30 '14 at 16:34
  • Would you happen to know this? http://stackoverflow.com/questions/25651480/does-scipy-logsumexp-deal-with-the-underflow-challenge – user961627 Sep 03 '14 at 19:03