Numerical Stability Analysis and Solutions for RuntimeWarning: invalid value encountered in double_scalars in NumPy

Problem Background and Phenomenon Analysis

In scientific computing, handling extreme numerical values is a common challenge. When using NumPy for computations involving exponential functions, developers frequently encounter the RuntimeWarning: invalid value encountered in double_scalars warning. This warning typically indicates invalid numerical operations during computation, such as division by zero or NaN values.

Consider the following typical scenario:

import numpy

d = numpy.array([[1089, 1093]])
e = numpy.array([[1000, 4443]])
answer = numpy.exp(-3 * d)
answer1 = numpy.exp(-3 * e)
res = answer.sum()/answer1.sum()

Executing this code produces a nan result accompanied by warning messages. The root cause lies in the answer1.sum() calculation resulting in 0, leading to division-by-zero errors.

Numerical Underflow Mechanism Analysis

When the parameters in the exponential function numpy.exp(-3 * e) are large negative numbers, the computation results approach 0. In double-precision floating-point representation, underflow occurs when values are smaller than approximately 1e-308, causing them to be approximated as 0. Detailed analysis:

For denominator calculation: numpy.exp(-3 * 1000) + numpy.exp(-3 * 4443)

numpy.exp(-3000) ≈ 0 (actual value is approximately 5.148e-1304, far below double-precision minimum)

numpy.exp(-13329) ≈ 0 (actual value is even smaller)

Therefore, the denominator sum results in 0, causing the division operation to produce nan.

Mathematical Derivation and Stability Optimization

Numerical underflow issues can be avoided through mathematical transformations. Consider the general form of computation:

result = (exp(a) + exp(b)) / (exp(c) + exp(d))

Apply logarithmic identities for transformation:

log(result) = log(exp(a) + exp(b)) - log(exp(c) + exp(d))

Further derivation for the numerator part:

log(exp(a) + exp(b)) = a + log(1 + exp(b - a))

This transformation converts large-value exponential operations into relatively smaller logarithmic operations, effectively avoiding numerical underflow. In the original problem, let x = 3*1089, y = 3*1093, z = 3*1000, k = 3*4443, then:

log(numerator) = -x + log(1 + exp(-y + x))
log(denominator) = -z + log(1 + exp(-k + z))

Since exp(-k + z) is extremely small, log(1 + exp(-k + z)) ≈ 0, therefore:

log(result) ≈ -x + log(1 + exp(-y + x)) + z

NumPy Built-in Function Solution

NumPy provides the np.logaddexp function specifically designed for such problems:

import numpy as np

d = np.array([[1089, 1093]])
e = np.array([[1000, 4443]])

log_res = np.logaddexp(-3*d[0,0], -3*d[0,1]) - np.logaddexp(-3*e[0,0], -3*e[0,1])
res = np.exp(log_res)

print(f"Computation result: {res}")
print(f"Logarithmic result: {log_res}")

The np.logaddexp function implements numerically stable computation of log(exp(a) + exp(b)), avoiding intermediate exponential operations.

SciPy Extension Solution

For more complex multi-element summations, SciPy's logsumexp function can be used:

import numpy as np
from scipy.special import logsumexp

d = np.array([[1089, 1093]])
e = np.array([[1000, 4443]])

res = np.exp(logsumexp(-3*d) - logsumexp(-3*e))

print(f"SciPy computation result: {res}")

The logsumexp function can handle input arrays of any dimension, providing better versatility.

Numerical Accuracy Verification

Comparison between approximate solutions obtained through mathematical derivation and exact solutions:

# Approximate solution
approx_result = np.exp(-266.99999385580668)

# Exact solution (via high-precision computation)
# Wolfram Alpha: 1.105034914720621496 × 10^-116

error = abs(res - approx_result)
relative_error = error / res

print(f"Relative error: {relative_error:.2e}")

Actual computations show that the relative error between results obtained using numerically stable methods and exact solutions is on the order of 1e-13, fully meeting engineering computation requirements.

Preventive Programming Practices

In addition to using specialized numerically stable functions, the following preventive measures can be adopted:

Input Validation: Check input data reasonability before computation

def safe_division(numerator, denominator):
    if denominator == 0:
        # Return appropriate value based on business logic
        return 0.0  # or np.nan
    return numerator / denominator

Numerical Range Checking: Identify parameter ranges that may cause numerical issues

def check_exponent_safety(exponent, threshold=700):
    """Check if exponent parameters are safe"""
    if np.any(np.abs(exponent) > threshold):
        print("Warning: Exponent parameters may cause numerical issues")
        return False
    return True

Summary and Best Practices

Addressing numerical stability issues in NumPy requires comprehensive consideration of mathematical principles and programming practices. Key points include:

1. Understanding the limitations of floating-point representation, particularly underflow and overflow boundaries

2. Mastering the mathematical foundations of numerical stability techniques like log-sum-exp

3. Proficiently using specialized functions provided by NumPy and SciPy

4. Adding appropriate validation and fault-tolerant mechanisms at critical computation stages

Through the methods introduced in this paper, developers can effectively avoid the RuntimeWarning: invalid value encountered in double_scalars warning, ensuring accuracy and stability in scientific computations.