Continuing the discussion from Formalizing Julia in Lean (and finding 400+ bugs in Base):
I see at the Java docs [see linked from the issue] even the sine function is guaranteed semi-monotonic… Semi-monotonic doesn’t imply correctly rounded I think i.e. it’s a stricter goal. I just discovered the non-monotonic bug and existence of semi- definition. Correctly rounding is a trade-off might be slower, not sure about semi-.
I think monotonic/semi-monotonic should be cheaper than correctly rounded (although it’s of course much harder for Float64 especially on multiple argument functions)
Perhaps a naive question, but is there any scenario where this non-monotonicity would be problematic? If you’re down to doing math at the ulp-scale, don’t you have much bigger problems than monotonicity?
I don’t have a specific experience. And indeed the error itself is tiny. But it’s an error in the order, which feels like it could cause silent damage somewhere. tanh(x) being monotonic is just something I’ve always assumed, so surely others would as well, and that broken assumption could lead to bugs.
For example say you have a grid and map a tanh over it, like tanh.(grid). Down the pipeline you then run an algorithm that assumes a sorted input. Well you had thought it was sorted, because your grid was sorted, and tanh is monotonic, but you never knew about this quirk. And thus your algorithm generates some weird bug because the result was not in fact sorted.
I mean maybe there’s some sort of derivative bug too. Like you could imagine that I might write code that assumes a finite-difference gradient is >= 0 for a monotonic function. But with a one-ulp step at x = 0.49214557f0 it comes out as -1.0. So anything that branches on that sign goes the wrong way, or throws a DomainError if you take a sqrt of it.
There’s just tons of things that I think could go wrong when you have commonsense correctness issues like this. Sure it could be documented, but most people will not read the docs for something as simple as tanh (my AI agent included). So any convention-breaking behavior might bite them.
Julia doesn’t guarantee commutativity and distributivity of simple arithmetic with floating point, the native floating point operations do not. I see no particular reason why it should guarantee monotonicity, or any other mathematical property involving reals. Floating point computations are inherently approximate.
But it’s not about mathematical properties of reals, it’s about ordering. Floating point numbers are ordered, so it would make sense that passing them through a monotone function like tanh would preserve the ordering. Monotinicity survives rounding as well, the type of computations you’re doing matter.
It’s a tall order:
julia> x = 1.0
1.0
julia> nextfloat(x) > x
true
julia> sqrt(nextfloat(x)) > sqrt(x)
false
sqrt is strictly monotonic. Would it be reasonable to preserve strict monotonicity as well?
With floating point I think the only thing that can be guaranteed is that x > y ⟹ f(x) > f(y) ± ϵ for f approximating an increasing function, and ϵ some preferably small number depending on x, y, and f. Everything concerning reals is only approximated by floating point.
One could make provisions in system libraries to ensure preservation of monotonicity, though it would have to be on a case by case basis. This boils down to how these functions are actually computed, with Taylor expansions or other approximations utilizing native operations which are only approximate. E.g. (a + x) - x is only approximately equal to a. It would be very hard to ensure monotonicity for every conceivable user-written function which is provably monotonic in a mathematical sense.
No, because strict monotonicity is not preserved under rounding.
But even simple native arithmetic is bad in this sense, and can reverse the order:
julia> f(x) = (1.0 + x)^2 - 1 - 2x # mathematically equal to x^2.
f (generic function with 1 method)
julia> x = 0.7385297827525479;
julia> y = 0.738529782752548:
julia> x < y
true
julia> f(y) < f(x)
true
julia> y^2 < x^2
false
sure, you’re not going to be able to make any such claims for user generated functions unless maybe in some special cases. But for common math functions in base like tanh it should be doable or at least highly desired.
That doesn’t look simple at all to me btw
That’s pretty much hopeless for most functions. Even a simple monotonic function like f(x) = x / 2 on the interval [0, 1], which for most inputs is computed exactly, must fail simply because there are fewer floating point numbers available in the output domain [0, 0.5] than in the input domain [0, 1].
I was thinking if you calculate gradients, e.g. with tanh, then they will be wrong. Tanh is or was used by neural networks, not other activation functions are more common. But I believe it’s still used in other places in networks, maybe gradients apply there too. The gradients are done numerically, I’m assuming, though this non-monotonicity is I guess only apparent if you use the nextfloat value or at best a very close one, and you likely have larger steps. It’s easy to check if this error (for e.g. tanh) also appears for larger steps.
You must test for:
julia> sqrt(nextfloat(x)) >= sqrt(x)
true
I.e. is it non-decreasing, since you compress your values into a smaller range. It’s impossible to avoid the equality, i.e. ϵ == 0, but that should be the only exception to increasing (the mathematical truth), because of finite floating point.
So, what’s the promise julia should make? Something like this: Whenever you have a differentiable mathematical function with f'(x) > 0 in x, its julia implementation should have f(nextfloat(x)) >= f(x)? Or should it avoid derivatives, and use monotonic in a neighbourhood instead?
Yes I think so. At first we should just aim for this, e.g. fix tanh (also sin?), then see if we can do for all relevant and then provide a guarantee. I think if we adopt correctly rounded we get this for free, but that can only be tested for up to Float32. Here’s their definition (can it be done for Float64 too?):
Therefore, most methods with more than 0.5 ulp errors are required to be semi-monotonic: whenever the mathematical function is non-decreasing, so is the floating-point approximation, likewise, whenever the mathematical function is non-increasing, so is the floating-point approximation. Not all approximations that have 1 ulp accuracy will automatically meet the monotonicity requirements.
Numerically yes, but usually with other methods than finite differences. And if you do use finite differences, the step anyway needs to be quite a bit larger than the scale at which you might see a non-monotonicity.
I, meaning AI, did a lot of analysis, and first I report verifiable parts:
Here it is — the first reversal in that search whose output jump exceeds 1 ulp. It’s tanpi on the increasing branch (0, 0.5):
x = 0.2356103819462645
xn = nextfloat(x)tanpi(x) = 0.9134413873825893
tanpi(xn) = 0.9134413873825891 # 2 ulps down, wrong for an increasing branch
julia> bitstring(tanpi(x))
“0011111111101101001110101110100101101110101100111101111111110100”julia> bitstring(tanpi(xn))
“0011111111101101001110101110100101101110101100111101111111110010”
The two bit patterns differ in the last two bits, confirming it is a 2-ulp drop, not a printing artifact.
It recovers immediately after that:
step 0: 0.9134413873825893 (start)
step 1: 0.9134413873825891 (-2 ulp) ← reversal
step 2: 0.9134413873825894 (+1 ulp) ← recovered
step 3: 0.9134413873825895 (+2 ulp)
step 4: 0.9134413873825897 (+4 ulp)So the value drop is 2 ulps, but the input gap is still 1 step — this remains a single adjacent-pair inversion, not a sustained reversed interval. Of the 10 hits printed by that search, the first six were exactly 1 ulp; this was the first 2-ulp one.
Some discussion I had regarding my concerns:
Great point on tanpi — correct rounding is not preserved by division, so correctly-rounded sinpi and cospi do not imply a correctly-rounded or monotone tanpi. Let me actually certify Float16 exhaustively, and also test your “fully correctly rounded” hypothesis against BigFloat.
It didn’t “certify Float16” rather found issues I filed (an I have PR upcoming):
Confirmed cospi reversal (mathematically decreasing on [0,1], but the approximation increases):
x = 0.23383972257699212
xn = nextfloat(x)cospi(x) = 0.7420793865711964
cospi(xn) = 0.7420793865711965 # wrong direction
tanpi is interesting because it literally computes sinpi_kernel_wide(rx) / cospi_kernel_wide(rx) — i.e. it consumes the same cospi kernel that we just showed is broken. So it’s a good candidate to inherit the reversal. tan uses its own tan_kernel with a double-double path, so it may be clean.
For tanh we found adjacent reversals, so I computed the maximum gap: 1 ulp. For sin, cos, sinpi, and cospi, the scans have found no adjacent reversals, so there cannot be a 2-ulp or larger gap inside any monotonic interval.
tanpi x=0.2356103819462645 tanpi(x)=0.9134413873825893 tanpi(next)=0.9134413873825891 delta=-2.220446049250313e-16 ulps=-2.0 BAD=true
tanpi x=0.24685571942890097 tanpi(x)=0.9804365152578822 tanpi(next)=0.980436515257882 delta=-2.220446049250313e-16 ulps=-2.0 BAD=true
tanpi x=0.22717042552698963 tanpi(x)=0.8659423889513675 tanpi(next)=0.8659423889513673 delta=-2.220446049250313e-16 ulps=-2.0 BAD=true
tanh x=-0.4843550140218486 tanh(x)=-0.4497246087273284 tanh(next)=-0.44972460872732845 delta=-5.551115123125783e-17 ulps=-1.0 BAD=true
cospi x=0.23383972257699212 cospi(x)=0.7420793865711964 cospi(next)=0.7420793865711965 delta=1.1102230246251565e-16 ulps=1.0 BAD=true
cospi x=0.2445573715344167 cospi(x)=0.7190933104883361 cospi(next)=0.7190933104883362 delta=1.1102230246251565e-16 ulps=1.0 BAD=true
cospi x=0.2396223272407172 cospi(x)=0.7297803251322194 cospi(next)=0.7297803251322195 delta=1.1102230246251565e-16 ulps=1.0 BAD=true
cospi x=0.22961747798032528 cospi(x)=0.7509052450302832 cospi(next)=0.7509052450302833 delta=1.1102230246251565e-16 ulps=1.0 BAD=true
On “only a single gap”: for cospi I have one confirmed example so far, not a claim that it is the only one — it’s rare in my random windows. For tanh the random scan found 152 such reversals in 20M adjacent pairs. All of them, including the cospi one, were only a 1-step gap: the function recovers on the next step, so I found no larger gap. Let me search specifically for more cospi examples rather than relying on the rare random-window hit.
One of the codes used in analysis to find reversals:
using Random
const MAXSTEP = 64
function survey(name, f, lo, hi, increasing, n, seed)
rng = MersenneTwister(seed)
reversals = 0
maxulp = 0.0 # largest |drop| in ulps
maxulp_ex = nothing
maxgap = 0 # longest recovery in steps (1 = only step 1 is wrong)
maxgap_ex = nothing
for _ in 1:n
x = lo + rand(rng) * (hi - lo)
xn = nextfloat(x)
y = f(x)
yn = f(xn)
bad = increasing ? (yn < y) : (yn > y)
bad || continue
reversals += 1
u = abs((yn - y) / eps(y))
if u > maxulp
maxulp = u
maxulp_ex = (x, y, xn, yn)
end
# First k >= 1 with f(x_k) on the correct side of f(x0) again.
z = x
k = 0
while k < MAXSTEP
z = nextfloat(z)
k += 1
v = f(z)
ok = increasing ? (v >= y) : (v <= y)
ok && break
end
gap = k - 1
if gap > maxgap
maxgap = gap
maxgap_ex = (x, y, k, gap)
end
end
println(name, " reversals=", reversals,
" maxdrop_ulps=", maxulp, " maxgap_steps=", maxgap)
maxulp_ex !== nothing &&
println(" maxdrop example: x=", maxulp_ex[1], " y=", maxulp_ex[2],
" yn=", maxulp_ex[4], " ulps=", (maxulp_ex[4]-maxulp_ex[2])/eps(maxulp_ex[2]))
maxgap_ex !== nothing &&
println(" maxgap example: x=", maxgap_ex[1], " recovers at step ", maxgap_ex[3])
return reversals, maxulp, maxgap
end
# cospi: decreasing on (0,1); reversals are in ~[0.1,0.25]
survey("cospi", cospi, 0.05, 0.30, false, 2_000_000_000, 31)
# tanpi: increasing on (0,0.5); reversals are in ~[0.1,0.25]
survey("tanpi", tanpi, 0.05, 0.30, true, 2_000_000_000, 32)
# tanh: increasing; reversals are in |x| in ~[0.125,0.5]
survey("tanh+", tanh, 0.12, 0.50, true, 300_000_000, 33)
survey("tanh-", tanh, -0.50, -0.12, true, 300_000_000, 34)
(post deleted by author)
Pain.