TopolearnProbabilistic and scientific code fails silently with NaN, inf or wrong digits when it ignores floating-point limits.
Floating point has about 16 significant digits and a range up to roughly 1e308. Algorithms that are correct on paper fail in three classic ways.
Overflow and underflow. np.exp(1000) is inf; a product of 500 small probabilities is 0.0. Work in log space: sum log-probabilities instead of multiplying probabilities, and use scipy.special.logsumexp(a) to compute log(sum(exp(a))) safely by factoring out the maximum. special.softmax and special.log_softmax do the same for normalisation, and distributions offer logpdf, logcdf, logsf.
Cancellation. Subtracting nearly equal numbers loses digits: 1 - cos(x) or log(1 + x) for tiny x. Use np.log1p, np.expm1, special.cosm1, special.xlogy (handles 0 * log 0), and special.expit / special.logit for the sigmoid.
Large factorials and binomials. Use special.gammaln, special.betaln and special.comb(n, k, exact=True) rather than computing factorials directly.
Ill-conditioning. A problem with a large condition number amplifies input errors whatever the algorithm. Signs include LinAlgWarning: Ill-conditioned matrix, and optimisers that stall. Rescale variables, centre data, use orthogonal factorisations (QR, SVD) instead of normal equations, and regularise.
Also, tolerances such as rtol, xtol and atol are relative to scale; set them consciously rather than trusting defaults.
Example
import numpy as np
from scipy import special
logp = np.array([-1000.0, -1001.0, -1002.0]) # log-likelihoods
print(np.log(np.sum(np.exp(logp)))) # -inf + RuntimeWarning (exp underflows to 0)
print(round(special.logsumexp(logp), 4)) # -999.5924
x = 1e-12
print(np.log(1 + x), np.log1p(x)) # 1.000088900581841e-12 9.999999999995e-13
print(special.expit(-800.0), special.log_expit(-800.0)) # 0.0 -800.0
print(special.xlogy(0, 0)) # 0.0 (0 * log 0 treated as 0)
print(special.comb(1000, 500, exact=True) % 1000) # 320 (exact big integer)
print(round(special.gammaln(1001), 2)) # 5912.13 (log of 1000!)Remember
- Work in log space; use
logsumexp,log_softmax,logpdf - Use
log1p,expm1,xlogy,expitto avoid cancellation and overflow - Rescale and use orthogonal factorisations for ill-conditioned problems
Watch out for
- Multiplying many probabilities directly and getting 0.0 or NaN.