Floating point problems in asymptotic functions approaching zero - Python

Question

New to python coming from MATLAB.

I am using a hyperbolic tangent truncation of a magnitude-scale function. I encounter my problem when applying the 0.5 * math.tanh(r/rE-r0) + 0.5 function onto an array of range values r = np.arange(0.1,100.01,0.01). I get several 0.0 values for the function on the side approaching zero, which cause domain issues when I perform the logarithm:

P1 = [ (0.5*m.tanh(x / rE + r0 ) + 0.5) for x in r] # truncation function

I use this work-around:

P1 = [ -m.log10(x) if x!=0.0 else np.inf for x in P1 ]

which is sufficient for what I am doing but is a bit of a band-aid solution.

As requested for mathematical explicitness:

In astronomy, the magnitude scale works roughly as such:

mu = -2.5log(flux) + mzp # apparent magnitude

where mzp is the magnitude at which one would see 1 photon per second. Therefore, greater fluxes equate to smaller (or more negative) apparent magnitude. I am making models for sources which use multiple component functions. Ex. two sersic functions with different sersic indices with a P1 outer truncation on the inner component and a 1-P1 inner truncation on the outer component. This way, when adding the truncation function to each component, the magnitude as defined by radius, will become very large because of how small mu1-2.5*log(P1) gets as P1 asymptotically approaches zero.

TLDR: What I would like to know is if there is a way of preserving floating points whose accuracy is insufficient to be distinguishable from zero (in particular in the results of functions that asymptotically approach zero). This important because when taking the logarithm of such numbers a domain error is the result.

The last number before the output in the non-logarithmic P1 starts reading zero is 5.551115123125783e-17, which is a common floating point arithmetic rounding error result where the desired result should be zero.

Any input would be greatly appreciated.

@user:Dan without putting my whole script:

xc1,yc1 = 103.5150,102.5461;
Ee1 = 23.6781;
re1 = 10.0728*0.187;
n1 = 4.0234;
# radial brightness profile (magnitudes -- really surface brightness but fine in ex.)
mu1 = [ Ee1 + 2.5/m.log(10)*bn(n1)*((x/re1)**(1.0/n1) - 1) for x in r];

# outer truncation
rb1 = 8.0121
drs1 = 11.4792

P1 = [ (0.5*m.tanh( (2.0 - B(rb1,drs1) ) * x / rb1 + B(rb1,drs1) ) + 0.5) for x in r]

P1 = [ -2.5*m.log10(x) if x!=0.0 else np.inf for x in P1 ] # band-aid for problem

mu1t = [x+y for x,y in zip(P1,mu1)] # m1 truncated by P1 

where bn(n1)=7.72 and B(rb1,drs1) = 2.65 - 4.98 * ( r_b1 / (-drs1) );

mu1 is the magnitude profile of the component to be truncated. P1 is the truncation function. Many of the final entries for P1 are zero, which is due to the floating points being undistinguished from zero due to the floating point accuracy.

An easy way to see the problem:

>>> r = np.arange(0,101,1)
>>> P1 = [0.5*m.tanh(-x)+0.5 for x in r]
>>> P1
[0.5, 0.11920292202211757, 0.01798620996209155, 0.002472623156634768, 0.000335350130466483, 4.539786870244589e-05, 6.144174602207286e-06, 8.315280276560699e-07, 1.1253516207787584e-07, 1.5229979499764568e-08, 2.0611536366565986e-09, 2.789468100949932e-10, 3.775135759553905e-11, 5.109079825871277e-12, 6.914468997365475e-13, 9.35918009759007e-14, 1.2656542480726785e-14, 1.7208456881689926e-15, 2.220446049250313e-16, 5.551115123125783e-17, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

Note also the floats before zeros.

Answer 1

Recall that the hyperbolic tangent can be expressed as (1-e^{-2x})/(1+e^{-2x}). With a bit of algebra, we can get that 0.5*tanh(x)-0.5 (the negative of your function) is equal to e^{-2x}/(1+e^{-2x}). The logarithm of this would be -2*x-log(1+exp(-2*x)), which would work and be stable everywhere.

That is, I recommend you replace:

P1 = [ (0.5*m.tanh( (2.0 - B(rb1,drs1) ) * x / rb1 + B(rb1,drs1) ) + 0.5) for x in r]

P1 = [ -2.5*m.log10(x) if x!=0.0 else np.inf for x in P1 ] # band-aid for problem

With this simpler and more stable way of doing it:

r = np.arange(0.1,100.01,0.01)
#r and xvals are numpy arrays, so numpy functions can be applied in one step
xvals=(2.0 - B(rb1,drs1) ) * r / rb1 + B(rb1,drs1)
P1=2*xvals+np.log1p(np.exp(-2*xvals))

Answer 2

Two things you can try.

(1) brute force approach: find a variable-precision float arithmetic package and use that instead of built-in fixed precision. I am playing with your problem in Maxima [1] and I find that I have to increase the float precision quite a lot in order to avoid underflow, but it is possible. I can post the Maxima code if you want. I would imagine that there is some suitable variable-precision float library for Python.

(2) approximate log((1/2)(1 + tanh(-x)) with a Taylor series or some other kind of approximation in order to avoid the log(tanh(...)) altogether.

[1] http://maxima.sourceforge.net