Avoiding NaN and Inf when computing exponential in Matlab-CodePudding

Suppose I have a vector S and am interested in computing the value ln(1 exp(S)) or ln(1 exp(-S)).

However, if an entry is too large or too small, I will encounter Inf problems or other numerical issues.

One solution is to notice that S=max(S,0) min(S,0) and manipulate the expression ln(1 exp(±S)) so that it is more convenient for the computer. For exemple ln(1 exp(S))=max(S,0) log(exp(-max(S,0)) exp(min(S,0)))

It in fact does not give Inf values like ln(1 exp(S)) would have given. However, I need to do another circus for ln(1 exp(-S)). So I am looking for the most elegant, efficient way to do it (for both ln(1 exp(S)) and ln(1 exp(-S)) and without loops of course). Here is a code with the solution I presented, feel free to come up with a better one or improve this one

S = [1000 -1000];
x = log(1 exp(S))
y = max(S,0) log(exp(-max(S,0)) exp(min(S,0)))

The vector y is the one we want and that we expect mathematically.

CodePudding user response：

There are two issues to correct. When S is too large in magnitude and positive, the exp(S) overflows into inf before the log( ) takes place, so the raw calculation is invalid. When S is too large in magnitude and negative, exp(S) is too small to have an effect when added to 1, so you end up with log(1) in those cases which is also invalid. I would propose writing a short function to correct the extreme cases. A simplistic approach would suggest something like this for starters:

% Calculates log(1 exp(S)) with overflow/underflow corrections
function result = log1pex(S)
expS = exp(S); % need this later for small number corrections
result = log(1 expS); % raw function evaluations
x = S >= 37; % detect points about where exp(S) 1 == exp(S) , doesn't need to be exact
result(x) = S(x); % correct those points
x = result == 0; % detect points where 1 exp(S) == 1
result(x) = expS(x); % correct those points
end

The check for S >= 37 simply detects when exp(S) will be too large and swamp the 1 in double precision, so the result can be taken to be log(exp(S)) = S. This check doesn't need to be precise ... it simply needs to cover the cases that cause inf problems.

The check for result == 0 is simply detecting the case where exp(S) is too small and gets swamped by the 1 in double precision. In that case simply apply the Taylor series for log(1 y) ~ y, where in this case y is exp(S).

I will state I haven't made a detailed examination of this latter case for best detection and formula to use. Maybe a slightly positive number would work better for detection and a couple of more Taylor series terms might be better in some cases to give a more accurate final result. E.g.,

tol = some small positive number
x = result < tol; % detect where 1 exp(S) ~ 1 , or could test S directly instead of result
result(x) = expS(x) - (expS(x).^2)/2   (expS(x).^3)/3; % correct those points

Or you could write a loop for the log(1 y) expansion to accommodate all terms that affect the result in double precision. If I get some time I may look at this later. But you get the idea.

EDIT

I took a quick look at log(1 exp(S)) vs log1pex(S) for the small corrections. When exp(S) is smaller than about 1e-8, two terms of the Taylor series are enough. And when exp(S) is smaller than about 1e-12 or so the log(1 exp(S)) raw calculation starts to become noticeably jagged when plotting, while the Taylor series calculation is smooth. So maybe this:

x = expS < 1e-8; % detect where 1 exp(S) ~ 1
result(x) = expS(x) - (expS(x).^2)/2; % correct those points