NaNs when computing exponential and normalizing

I have an mxn matrix A, with very small or very large elements. For example:

    A = -1e4*randn(10,20);

I would like to create a new matrix C of the same size, as follows:

First, define a matrix B whose elements are the exponential of the elements of A:

    B = exp(A)

Then, the matrix C is defined such that each column of C is proportional to the corresponding column of B and the sum of each column of C is equal to 1. In other words, if we take an element of B and divide it by the sum of all the elements of the same columns, then we obtain the corresponding element of C:

    C = bsxfun(@rdivide,B,sum(B));

Mathematically:

Clearly, all the elements of C are between 0 and 1. However, when computing using the following code, the obtained matrix contains many NaNs:

    A = -1e4*randn(10,20);

    B = exp(A);

    C = bsxfun(@rdivide,B,sum(B));

    sum(sum(isnan(ee)))

I hope somebody can suggest a way to overcome this issue. Thank you very much in advance.

Answer 1

[John D'Errico]

Sigh. You have won the award for being the 10 millionth person to trip over the pitfalls of floating point arithmetic and the issue of a limited dynamic range. No matter what limits were imposed of course, SOMEBODY will try to exceed them. It is in the nature of computing that people assume their computers can do anything, that they are all powerful. After all, the computers we see in the movies are exactly that.

As well, despite what you may think, floating point arithmetic is NOT mathematics, although it is vaguely related. So while you may be able to write some expression down that makes sense mathematically, in terms of floating point arithmetic, you may see a meaningless result if you exceed the limits of the precision in which you work.

A double has a limited exponent range, with an exponent that will not allow you to form numbers larger then 1e308 or so, and about 1e-308 in the other direction. So when you try to form a number like exp(1e4), of course it fails. The number is simply not representable as a double.

So if you intend to compute with numbers of those magnitudes, then you need to work with a tool that can do those computations, OR reformulate things so that you stay in the workable range of a double. Remember that in in terms of a double,

1 + 1e-17 == 1
ans =
     1

and even if the difference in magnitude is smaller than that, you will still see a massive loss of information when the numbers you add and subtract vary by many orders of magnitude. So for example, we have roughly a 10% error in this simple result:

1 + 1e-15 - 1
ans =
      1.11022302462516e-15

So the computations you show will surely yield meaningless results even if you did not have NaNs confounding the problem.

Again, you will need to work with tools that can handle these numbers. So either the symbolic toolbox with vpa, or my own HPF toolbox. Remember though that numbers as large as exp(1e4) are on the order of magnitude of 1e4300.

DefaultNumberOfDigits 50
exp(hpf(1e4))
ans =
8.8068182256629215872614960076445610035200040855913e4342

So any arithmetic on numbers that vary by 4000+ decimal digits must work in more precision than that, or risk too much loos in the lower bits.

[f10w]

Hi @John. Thanks for the answer.

When I posted the question, I thought that there should be some analytical tricks that help to overcome that numerical issue. Such a well-known trick is the log-sum-exp trick.

With the same idea, we can compute C using:

Thus, the new code is:

A = bsxfun(@minus,A,max(A));
B = exp(A);
C = bsxfun(@rdivide,B,sum(B));

I have the feeling that this might give a very good approximation, but I don't know how good it is. What do you think about that?

Maybe we can try using the toolboxes that you suggested and then comparing the results with the above method.

[Roger Stafford]

@Khue. Yes, I was about to suggest something of the sort until I read your comment. It is a perfectly valid method.

[f10w]

Thanks, @Roger.