8

I need to compute a normalized exponential of a vector in Matlab.

Simply writing 

res = exp(V)/sum(exp(V))

overflows in an element of V is greater than log(realmax) = 709.7827. (I am not sure about underflow conditions.)

How should I implement it to avoid numerical instability?

Update: I received excellent responses about how to avoid overflow. However, I am still happy to hear your thoughts about the possibility of underflow in the code.

user25004
  • 1,868
  • 1
  • 22
  • 47

2 Answers2

9

The following approach avoids the overflow by subtracting the exponents and then taking the exponential, instead of dividing the exponentials:

res = 1./sum(exp(bsxfun(@minus, V(:), V(:).')))

As a general rule, overflow can be avoided by working in the log domain for as long as possible, and taking the exponential only at the end.

Luis Mendo
  • 110,752
  • 13
  • 76
  • 147
  • Thanks. I am new to bsxfun. What if V is a matrix, and we want this to apply along dimension dim? – user25004 May 14 '14 at 22:00
  • @user25004 How would that be? `exp(V)` would then be a matrix and `sum(exp(V),dim)` would be a vector. How do you define `exp(V)/sum(exp(V,dim))` in that case? – Luis Mendo May 14 '14 at 22:03
  • I mean if the variable called dim is 1, I want your previous code to be applied to each column. If dim is 2, the code is applied row-wise. – user25004 May 14 '14 at 22:14
  • @user25004 That would be harder. And it changes the question completely. The simple answer would be: loop over each row or column. – Luis Mendo May 14 '14 at 22:15
  • @user25004 Welcome! Feel free to ask that new requirement as a new question, if the loop solution is not good enough – Luis Mendo May 14 '14 at 22:17
  • Nice solution, should have a better precision than my solution. – Daniel May 14 '14 at 23:00
4

The answer is pretty similar to your previous question. Use Math!

exp(V)=exp(V-max(V))*exp(max(V))
sum(exp(V))=sum(exp(V-max(V))*exp(max(V)))=exp(max(V)*sum(exp(V-max(V))))

Putting both together:

res=exp(V-max(V))*exp(max(V))/exp(max(V)*sum(exp(V-max(V)))=exp(V-max(V))/sum(exp(V-max(V)))

A code which is robust to the input range:

res=exp(V-max(V))/sum(exp(V-max(V)))
Daniel
  • 36,610
  • 3
  • 36
  • 69