NaNs when computing exponential and normalizing

25 views (last 30 days)
f10w
f10w on 17 Jul 2014
Commented: f10w on 21 Jul 2014
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.

Answers (1)

John D'Errico
John D'Errico on 17 Jul 2014
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.
  3 Comments
Roger Stafford
Roger Stafford on 18 Jul 2014
@Khue. Yes, I was about to suggest something of the sort until I read your comment. It is a perfectly valid method.

Sign in to comment.

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!