Is Julia reliable for solving ordinary and stochastic differential equations?

Palli · May 6, 2023, 3:50pm

I see from last month:
2023: MIT Center for Computational Science & Engineering MathWorks Prize for Outstanding Master’s Research in Computational Science & Engineering. For Songchen Tan (MIT Julia Lab) for his work on TaylorDiff.jl for use NeuralPDE.jl physics-informed neural networks (PINNs)

So it’s a Prize (also) in the name of MathWorks, i.e. the maker of MATLAB. The SciML ecosystem if for sure doing something (well man things) right. The only language (of those) potentially competing on speed with Julia is C++. The others not unless you can wrap fast libraries from other fast languages, and that can be done for e.g. BLAS, but I understand more challenging for solving equations.

I believe all the issues Yuri brought up, are solvable, if not already, and the only really interesting intriguing issue he brought up that time was OffsetArrays.jl, that I very much doubt is an issue for solving your equations:

There is a lot of problems with OffsetArrays.jl but most other languages do not even something like OffsetArrays.jl or the expectation that most code written would automatically work with custom indexes.

I (still) have full confidence in Julia core developers and e.g. Chris, for math/technical computing, such as you bring up. Julia however isn’t bug-free (no software ever will be), but issues are handled well.

Yuri no longer has any open PR in Julia, all of his 11 have been merged or closed. But he also has open issues, including 2 from 2 weeks ago, rather interesting:

Issues · JuliaLang/julia · GitHub [after clicking the link the + needs to be changed to a space; I’m not sure how to put in a correct link for Discource, it breaks if I do it, so I’m not sure it’s possible…]

Note, I think all integer and floating-point division is correct in Julia for all types (handled by CPU instructions, i.e. assuming correctly), from his issue, it’s about floored division (e.g. div)

github.com/JuliaLang/julia

Incorrect results from floating-point division (`fld`, `cld`, `div`)

opened 12:38AM - 21 Apr 23 UTC

yurivish

bug maths float16 correctness bug ⚠

Julia has a [`div`](https://docs.julialang.org/en/v1/base/math/#Base.div) functi…on that is used to implement [floor division](https://docs.julialang.org/en/v1/base/math/#Base.fld) (`fld`) and [ceil division](https://docs.julialang.org/en/v1/base/math/#Base.cld) (`cld`). I think I’ve found a bug that causes all of these division functions to return incorrect results. Floor division is documented to return the largest integer less than or equal to `x / y`. It should never return a number that is *greater* than `x / y`. But it does: ```julia julia> 514 / Float16(0.75) Float16(685.5) julia> div(514, Float16(0.75)) # should round down, but rounds up instead Float16(686.0) julia> fld(514, Float16(0.75)) # likewise Float16(686.0) julia> fld(514, Float16(0.75)) ≤ 514 / Float16(0.75) false ``` Similarly, ceil division should never return a number that is *smaller* than regular division, but it does: ```julia julia> 515 / Float16(0.75) Float16(686.5) julia> cld(515, Float16(0.75)) # should round up, but rounds down instead Float16(686.0) julia> cld(515, Float16(0.75)) ≥ 515 / Float16(0.75) false ``` This behavior is not limited to 16-bit floats. Here’s a case where `fld` produces an incorrect result for `Float32` inputs: ```julia julia> 4_194_307 / Float32(0.75) # = 5592409.5 5.5924095f6 julia> fld(4_194_307, Float32(0.75)) # = 5592410, incorrectly rounded up 5.59241f6 julia> fld(4_194_307, Float32(0.75)) ≤ 4_194_307 / Float32(0.75) false ``` And here’s the same for `cld`: ```julia julia> 4_194_308 / Float32(0.75) # = 5592410.5 julia> cld(4_194_308, Float32(0.75)) # = 5592410, incorrectly rounded down 5.59241f6 julia> cld(4_194_308, Float32(0.75)) ≥ 4_194_308 / Float32(0.75) false ``` The equivalent operations in Python produce the correct results: ```python # For 16-bit floats: >>> np.float16(514) / np.float16(0.75) # regular division 685.5 >>> np.float16(514) // np.float16(0.75) # floor division 685.0 # For 32-bit floats: >>> np.float32(4_194_307) / np.float32(0.75) # regular division 5592409.5 >>> np.float32(4_194_307) // np.float32(0.75) # floor division 5592409.0 ``` Examples of this incorrect behavior are not hard to find – for most floats, you can find a divisor that will make either `fld` or `cld` return the wrong answer. Here are some examples for `Float16` where either `fld` or `cld` is incorrect: - `cld(1, Float16(0.000999)) < 1 / Float16(0.000999)` - `cld(2, Float16(0.001999)) < 2 / Float16(0.001999)` - `cld(3, Float16(0.002934)) < 3 / Float16(0.002934)` - `cld(4, Float16(0.003998)) < 4 / Float16(0.003998)` - `fld(5, Float16(0.004925)) > 5 / Float16(0.004925)` And here are some for `Float32`: - `fld(5, Float32(6.556511e-7)) > 5 / Float32(6.556511e-7)` - `fld(10, Float32(1.3113022e-6)) > 10 / Float32(1.3113022e-6)` - `fld(11, Float32(1.4305115e-6)) > 11 / Float32(1.4305115e-6)` - `cld(16, Float32(2.8014183e-6)) < 16 / Float32(2.8014183e-6)` - `cld(17, Float32(2.2053719e-6)) < 17 / Float32(2.2053719e-6)` For simplicity I’ve presented examples where the first argument is an integer; this bug also occurs for non-integral inputs. A divisor producing the wrong result can be found for over 51% of all possible 16-bit floats. I have not evaluated how widespread this is for `Float32`, but the results above suggest that it is similarly easy to find failures there too. I’ve tracked the invocations down to [this definition](https://github.com/JuliaLang/julia/blob/72aec423c2ab9f80c249d63fdd68b35833cfd7ed/base/div.jl#L370) of `div`, which has existed at least as far back as [Julia 1.6](https://github.com/JuliaLang/julia/blob/f9720dc2ebd6cd9e3086365f281e62506444ef37/base/div.jl#L279): ```julia div(x::T, y::T, r::RoundingMode) where {T<:AbstractFloat} = convert(T, round((x - rem(x, y, r)) / y)) ```

That error rate suggests that this bug manifests once out of every ~1,700 randomly chosen Float16 divisions
[…]
Based on a billion samples, the bug manifests once out of every ~10,000 randomly chosen Float32 divisions and once out of every ~74,000 Float64 divisions.

Division (and e.g. square root) CAN be correct, in general, can only be correctly rounded, so you always have a small error. According to chaos theory, small errors can blow up. I don’t think this is too worrying, i.e. to have a tiny bit smaller error, than the correctly rounded result. If you’re sensitive to it, then most likely the small error in the correctly rounded error too?

I’m curious, for differential equations, you can have arbitrary operators, e.g. division common, but would you have floored division often (or ever), or div? If not, then you don’t need to worry about that issue at least.

Since Yuri was posting new issues recently, he at least cares about Julia, and likely still uses, or why would he even have know or posted those two weeks ago.