log and rsqrt epsilon defaults underflow to zero in fp16, so folding returns inf #2886

What

log and rsqrt cast their default epsilon to the input dtype:

# elementwise_unary.py
epsilon=nptype_from_builtin(self.x.dtype)(1e-45)   # log
epsilon=nptype_from_builtin(self.x.dtype)(1e-12)   # rsqrt

Both values are representable in fp32, but np.float16(1e-12) and np.float16(1e-45) are both exactly 0.0. So in fp16 the epsilon the docstrings describe as "a small constant added for stability" is removed entirely, and constant folding returns infinities:

                epsilon                  folded value for x = [0.0, 1e-8]
rsqrt  fp16     0.0                      [inf inf]
rsqrt  fp32     1e-12                    [1.0000e+06 9.9995e+03]
log    fp16     0.0                      [-inf -inf]
log    fp32     1e-45                    [-103.27893  -18.420681]

fp16 is the dtype most Core ML models are converted to, so the stabilizer is missing in the case it exists for.

Reproduction

import numpy as np
from coremltools.converters.mil.mil import Builder as mb

for op_name in ("rsqrt", "log"):
    for dt in (np.float16, np.float32):
        x = np.array([0.0, 1e-8], dtype=dt)

        @mb.program(input_specs=[mb.TensorSpec(shape=(1,))])
        def prog(unused):
            return getattr(mb, op_name)(x=x)

        f = list(prog.functions.values())[0]
        op = [o for o in f.operations if o.op_type == op_name][0]
        print(op_name, dt.__name__, np.array(op.epsilon.val), np.array(op.outputs[0].val))

coremltools 9.0, macOS arm64.

A second, related observation

While narrowing this down I found that constant folding and the runtime disagree for fp32 rsqrt, independently of the fp16 problem. The same op, with the same input, gives different answers depending on whether the input is known at convert time:

         x   const-folded        runtime     documented (1/sqrt(x+1e-12))
     0e+00   1000000.0000      4096.0000   1000000.0000
     1e-08      9999.5000      3790.3638      9999.5001
     1e-06       999.9995       971.4670       999.9995
     1e-04       100.0000        99.9702       100.0000

The folded path matches the documented y = 1 / sqrt(x + epsilon) with epsilon = 1e-12. The runtime behaves as though epsilon is roughly 2**-24 (~6e-08). That is a 244x difference at x = 0.

I have not filed that as a separate issue because I cannot tell which side is intended, and the runtime is not something coremltools can change. If the runtime's epsilon is deliberate, then value_inference is producing values the hardware will never reproduce, which seems worth deciding on explicitly. Happy to split it out if you'd prefer.

Interestingly, the fp16 fix in the linked PR makes fp16 rsqrt(0) fold to 4096.0 — the same value the runtime returns — so for fp16 the two paths end up agreeing.