There are a couple of different contexts at work here.
Maybe you are doing essentially a one-off sanity check, like interactively eyeballing an answer to see if a result is not absurd. But using isapprox() instead of your eyes.
Maybe you are about to cement-in a test in a routine in the midst of a lengthy, critical, and difficult-to-read program typically used by others.
Maybe like that Bessel program? ]Oddly enough I’ve written
a (too long- 21 page) paper on evaluating this function… https://people.eecs.berkeley.edu/~fateman/papers/hermite.pdf ]
I’m reminded of the ChatGPT (etc) programs that offer to write you a program from a natural language description. The implicit (or
sometimes explicit) warning is that you should confirm that the program is correct. My thought is that if you use isapprox() in that second context, you should be warned.
A historical note… I seem to recall that Iverson’s APL, circa 1957) had a symbol for approximate equality; I haven’t tracked it down though. If no one else has mentioned it here, it occurs to me that it may be some sort of historical signpost. APL was (is?) worth knowing about.
Do you mean for program validation? Why is calling isapprox (ideally with custom tolerances, at least in tests to validate the precise accuracy as opposed to tests for coarse bugs) worse than writing your own similar test over and over (often done badly)?
That’s what I don’t get here. Definitely, you should know what any function (not just isapprox) does if you are going to call it! But what do you have against code re-use? All of the examples you have given of approximate comparisons were special cases of isapprox.
In brief, anyone reading/ perhaps modifying/ your code will have to know (or ignore) the meaning of isapprox. Also maybe
find out where you have set or defined atol, rtol, norm. If you are repeatedly comparing floating-point values for
approximate equality with similar but not identical code segments, written badly, then you have other problems.
Testing for abs(a-b)<test is really pretty clear at least for a,b, floating-point values, if that’s what you mean.
Offhand, I have not much uniform intuition as to how I might be comparing structures, matrices, vectors approximately
Say a,b are polynomials or symbolic expressions.
. Do they have similar coefficients? similar location of zeros? similar limiting values? similar some norm?
I agree, you should understand the functions you call. But the definition of isapprox is not especially complicated (if you know what relative and absolute error mean), and it is widely used for validation suites in the Julia ecosystem so it is not an obscure function. The hard part is understanding floating-point comparisons in the first place.
But most validation tests are not like this, they are relative-error tests, because for most floating-point comparisons it is more natural to set a relative tolerance than an absolute one. And those are much more cumbersome to write: norm(a - b) <= rtol * max(norm(a), norm(b)) or similar.
For one thing, this involves writing a and/or b twice, which means that you have to write a function or use temporary variables if they are expressions and you don’t want to write/eval the expressions twice. And if the arguments are arrays, many people get it wrong on their first try — they check approximate equality elementwise rather than using a norm to set the overall scale.
The Julia stdlib and packages have thousands of validation tests that check relative errors; I don’t understand why they should all have to rewrite this code, even though it is short. (Many of the existing tests could benefit from tightened tolerances, however; @test supports compact syntax such as @test x ≈ y rtol=1e-13 for this.)
Just sharing a comment to say that I appreciate the opinion of various experts in the subject. And to say that there are different objectives at play in the arguments.
From one side @stevengj seems to promote a deeper understanding of floating point arithmetic. From the other side @Richard_Fateman seems to value a friendlier end-user experience that escapes the floating point difficulties.
As someone not trained in the deepest details of floating points, I would really appreciate if a scientific language could hide these counter-intuitive results from me. I bet that the isapprox behavior near 0 is surprising for >95% of users of programming languages, including Julia users.
It seems like he’s advocating the opposite? He wants to remove/discourage the isapprox function and tell people to re-implement the desired comparison manually whenever they need it.
Yes, 95% of people think that 10^{-8} is approximately zero, because they’ve not internalized that to say whether something is big or small you first have to ask “compared to what?”
I think most people would compare it with what they know about floating points. If they know that Float64 has n digits of precision, they will create a threshold in their minds that discards anything beyond k digits. I know that k isn’t well-defined for all cases based on the comments above, but people still will apply similar criteria to pick “atol” implicitly in their applications. This “fallback” value that most people pick without thinking too much (assuming order of magnitude 1 for example) is something that could be fine tuned with ScopedValues.jl as suggested above.
My personal opinion is that this behavior would be more intuitive for end-users than simply returning false whenever one of the arguments is 0. These are all design choices after all, and if we know our audience is well served with the current choice, that is fine.
I think “most people” know next to nothing about floating points. If however we speak about (potential) Julia users, they come from very different backgrounds and will have very different notion about approximate equality.
For me, Unitful.jl is one of the largest Julia blessings. I’d let it convert all units to default SI units - meaning, even in our macroworld there would be numbers spanning a dozen of orders of magnitude both up and down. You cannot have the same atol for the Young’s modulus of steel (1e11 Pa) and the diffusion coefficient of hydrogen in steel (1e-15 m²/s).
I’m not convinced that the default atol=0 is such a problem for the average user? They might get surprised by it once, see the docstring where it elaborates on it (though maybe it could be given as an admonition to highlight it better), and then understand to choose a better tolerance from then on (which I believe is a good thing - you should be thinking about tolerances). As @stevengj already mentions, I don’t see what other good default could be chosen. Not everything can be so simple.
I also don’t think the note in isapprox’s docstring (or the main point behind this discussion) really requires anyone to know anything about floating point arithmetic anyway. It would be obvious (upon reflection, at least) that comparing zero to other numbers is hard since zero loses all information about units in the first place as well (among the other reasons given already).
Right. As I’ve said many times, many people’s mental model of floating-point arithmetic is closer to fixed-point arithmetic: they think of it as a fixed number of digits past the decimal point, rather than a fixed number of significant digits.
which are satisfied by x ≈ y (i.e. isapprox) in Julia (until you scale to overflow/underflow, of course, and NaN is not reflexive due to IEEE equality rules unless you pass nans=true). They aren’t satisfied if you had a nonzero default atol, because that would violate scale-invariance (3). (In contrast, numpy.isclose violates both 2 and 3.)
In particular, the APL tolerant comparison is defined by:
Tolerant comparison considers two numbers to be equal if they are within some neighborhood. The neighborhood has a radius of ⎕ct times the larger of the two in absolute value.
where ⎕ct is a “system variable”. This definition corresponds exactly to Julia’s isapprox(x,y, rtol=⎕ct), i.e. abs(x - y) ≤ rtol*max(abs(x),abs(y)).
I couldn’t find any clear description of their default system rtol value, and tastes could certainly vary on what choice to make here.@jar1 points out this paper which suggests that the default tolerance was 2^-32 ≈ 2.3e-10, a bit smaller than Julia’s √ε ≈ 1.5e-8 (in Float64 precision).
So, it seems the only difference between APL’s “tolerant comparison” and Julia’s x ≈ y might be the choice of the default rtol.
(As @juliohm noted above, Python 3.5’s math.isclose is also equivalent to ≈, except with a default rtol of 1e-9 … very similar to Julia’s choice except that it effectively assumes the precision is always Float64.)
The way I see it, the default atol of 0 means that you might get some unexpected test failures that you need to investigate and find out that you need to modify your tests or set an appropriate atol. I rather have that than having tests pass that shouldn’t because your data was scaled in a way which wasn’t appropriate for a non-zero default atol.