WIP: Remove `add_tiny` and add `estimate_magitude` #125

wesselb · 2020-11-30T13:50:13Z

The function add_tiny was a bit of a hack to deal with zeros. The new addition estimate_magnitude deals with this in a better way, at the cost of more function evaluations.

This should address #124:

julia> (p -> FiniteDifferences.central_fdm(p, 1)(cosc, 0)).(2:10) .+ (pi ^ 2) / 3
9-element Array{Float64,1}:
 -2.487109898532124
 -3.1522345644852123e-6
 -4.327381120106111e-9
 -8.283298491562618e-10
  1.3939516207983615e-11
  3.93636234718997e-11
  9.558132063602898e-12
 -2.1227464230832993e-13
  1.213251721310371e-12

The package can be pushed to estimate at high accuracies:

julia> FiniteDifferences.central_fdm(10, 1, adapt=3)(cosc, 0) + (pi ^ 2) / 3
-3.019806626980426e-14

This PR also adds a test that checks that above two cases.

willtebbutt

Thanks for sorting this out promptly!

Would it make sense to construct some test cases where the conditions present in Accuracy is a bit brittle #124 are nearly satisfied, for the sake of robustness testing? i.e. what if M is only very slightly greater than 0 or something?
Some unit testing for estimate_magnitude would be a good idea, particularly around input types that are not Float64s.
Integration testing for Float32s would be a good idea if you think we know how to construct them properly.

src/methods.jl

wesselb · 2020-11-30T14:41:37Z

To answer your questions,

I'm not sure. We could do that, but then we would have to determine an arbitrary threshold that determines what a pathological output of f(x) is. Admittedly, in perturbing x, we already have such an arbitrary threshold. I'm on the fence. What would you propose?
I've added type constraints and tests for estimate_magnitude.
Yeah, I fully agree that integration testing for Float32s is something that we should do properly, but my sense is that that is bigger than this PR. A whole test suite of hard cases for Float64s and Float32s and expected accuracies would be great.

willtebbutt · 2020-11-30T14:43:35Z

Hmmm good point. I think we're fine as-is (in this PR).
Thanks.
Fair point, let's leave it for the future.

willtebbutt

Approved subject to patch-version bump.

edit: we need CI back up and running. This cannot be merged until that happens.

wesselb · 2020-11-30T14:48:04Z

Yes, ugh, that's annoying. Moreover, @oxinabox will check tomorrow if these fixes help with some other accuracy issues that he was encountering, so perhaps good to also wait for that.

willtebbutt · 2020-11-30T14:51:17Z

In the mean time @giordano could you confirm if this branch resolves the issues that you've been encountering? In particular, are there any test cases over and above the one you provided in #124 that you've been having particular problems with?

giordano · 2020-11-30T14:56:15Z

Ok, I can test this later today. There were indeed a couple of other functions that were failing to achieve the same accuracy as Calculus, cosc was the first one. Thanks!

giordano · 2020-11-30T23:04:31Z

Another similar case, expected derivative in 0 is 0:

julia> Calculus.derivative(sinc, 0)
0.0

# v0.11.3
julia> (p -> FiniteDifferences.central_fdm(p, 1)(sinc, 0)).(2:10)
9-element Array{Float64,1}:
  0.0
  0.0
  5.444645613907266e-13
  2.7755575615628914e-16
 -9.538516174483089e-15
  2.42861286636753e-16
 -3.7404985543720125e-15
 -4.822111254688986e-16
  2.809173270009664e-15

# this PR
julia> (p -> FiniteDifferences.central_fdm(p, 1)(sinc, 0)).(2:10)
9-element Array{Float64,1}:
  0.0
  0.0
  5.444645613907266e-13
  4.312032877555794e-14
 -9.538516174483089e-15
  2.3967209127499836e-15
 -3.7404985543720125e-15
 -4.822111254688986e-16
  2.809173270009664e-15

I wouldn't say the accuracy in this case is bad, it's just that comparing with 0 is complicated.

Also the derivative of this function f(t) should be 0 (OK, I must admit that I have some very sick test cases)

julia> f(t) = (x -> QuadGK.quadgk(sin, -x, x)[1])(t)
f (generic function with 1 method)

julia> Calculus.derivative(f, 4)
7.065611032273702e-11

# v0.11.3
julia> (p -> FiniteDifferences.central_fdm(p, 1)(f, 4)).(2:10)
9-element Array{Float64,1}:
 -0.015243579281754268
  0.0
  4.131607608911521e-9
 -1.7704270075193033e-9
  2.7545394016191937e-9
 -8.188893009030013e-12
 -6.983545129351427e-11
 -3.0697586011503785e-11
 -2.164732008407795e-11

# this PR
julia> (p -> FiniteDifferences.central_fdm(p, 1)(f, 4)).(2:10)
9-element Array{Float64,1}:
  0.01641638262171357
 -1.5410326104441793e-5
  2.1935477946857812e-8
 -2.7022157015905038e-8
  9.047605793628396e-11
 -1.3604658271651211e-9
 -3.154690991530035e-11
  4.959996079211705e-11
  7.027698887063365e-13

giordano · 2020-11-30T23:18:00Z

In the end I think I'll default to 8 grid points in Measurements.jl and call it a day

wesselb · 2020-11-30T23:52:48Z

Ah, there was one other edge case that wasn't taken care of. The results look now as follows:

julia> (p -> FiniteDifferences.central_fdm(p, 1)(cosc, 0)).(2:10) .+ (pi ^ 2) / 3
9-element Array{Float64,1}:
 -2.487109898532124
 -3.1522345644852123e-6
 -4.327381120106111e-9
 -8.283298491562618e-10
  1.3939516207983615e-11
  3.93636234718997e-11
  9.558132063602898e-12
 -2.1227464230832993e-13
  1.213251721310371e-12

julia> (p -> FiniteDifferences.central_fdm(p, 1)(sinc, 0)).(2:10)
9-element Array{Float64,1}:
  0.0
  0.0
  5.444645613907266e-13
  4.312032877555794e-14
 -9.538516174483089e-15
  2.3967209127499836e-15
 -3.7404985543720125e-15
 -4.822111254688986e-16
  2.809173270009664e-15

julia> (p -> FiniteDifferences.central_fdm(p, 1)(f, 4)).(2:10)
9-element Array{Float64,1}:
 -2.3988103001624555e-13
  1.430803955495777e-12
 -2.4073864725037792e-14
 -2.2694686733052566e-15
  1.52180247233912e-13
 -1.1700430263506674e-15
  2.5255867408048164e-13
  7.01087210768824e-13
 -1.0381492304012043e-14

The latter two are near machine epsilon. I don't think you can get much better.

The hardest case is still cosc, but you can push it:

julia> FiniteDifferences.central_fdm(10, 1, adapt=3)(cosc, 0) + (pi ^ 2) / 3
-3.019806626980426e-14

You can also push Richardson extrapolation:

julia> FiniteDifferences.extrapolate_fdm(FiniteDifferences.central_fdm(2, 1), cosc, 0.0, contract=0.5)[1] + (pi ^ 2) / 3
-5.702105454474804e-13

giordano · 2020-12-01T00:07:13Z

With the latest version all current tests pass already with FiniteDifferences.central_fdm(4, 1), that's amazing! I'll probably also play with the other extrapolation methods, but this is looking much better now, thanks!

wesselb · 2020-12-01T00:08:31Z

Nice! I'm very happy to hear that. :)

oxinabox · 2020-12-01T17:42:52Z

src/methods.jl

+    # pathological input for `f`. Perturb `x`. Assume that the pertubed value for `x` is
+    # highly unlikely also a pathological value for `f`.
+    Δ = convert(T, 0.1) * max(abs(x), one(x))
+    return float(maximum(abs, f(x + Δ)))


should we make this recursive?

Suggested change

return float(maximum(abs, f(x + Δ)))

return estimate_magitude(f, x + Δ)

That would recurse infinitely for the function x -> 0.0. Moreover, a function might just actually be 0 in a neighbour around x.

oxinabox · 2020-12-01T18:34:14Z

src/methods.jl

+    # Estimate the round-off error. It can happen that the function is zero around `x`, in
+    # which case we cannot take `eps(f(x))`. Therefore, we assume a lower bound that is
+    # equal to `eps(T) / 1000`, which gives `f` four orders of magnitude wiggle room.
+    ε = max(eps(estimate_magitude(f, x)), eps(T) / 1000) * factor


should we move this to a little function?
estimate_roundoff_error

Totally! Pushed the change.

oxinabox · 2020-12-01T18:56:16Z

Yep, I can confirm that this fixes the problems we were having with ChainRules

oxinabox · 2020-12-02T17:46:39Z

I think this needs rebasing so CI will run.
Other than that this LGTM, hopefully can merge soon?

wesselb · 2020-12-03T13:57:32Z

@oxinabox The tests fail on nightly. Are we okay with that?

oxinabox · 2020-12-03T16:31:11Z

@oxinabox The tests fail on nightly. Are we okay with that?

yes this is fine.

wesselb added 3 commits November 30, 2020 13:13

Remove add_tiny and add estimate_magitude

2c1d492

Improve comment

c5a4b40

Add test for #124

795b2ef

wesselb requested review from willtebbutt and oxinabox November 30, 2020 13:50

willtebbutt requested changes Nov 30, 2020

View reviewed changes

src/methods.jl Outdated Show resolved Hide resolved

src/methods.jl Show resolved Hide resolved

src/methods.jl Outdated Show resolved Hide resolved

src/methods.jl Outdated Show resolved Hide resolved

Remove inlines and add tests

65c38a0

Fix equality checking in tests

15ef56c

willtebbutt approved these changes Nov 30, 2020

View reviewed changes

wesselb added 2 commits November 30, 2020 14:46

Fix test

f4b5dcb

Increase patch version

a464d5d

wesselb added 2 commits November 30, 2020 14:55

Actually test the edge case

be02916

Add missing change

c7a9b97

Convert asserts to tests

478f3ce

Fix and test one other edge case

740463a

oxinabox reviewed Dec 1, 2020

View reviewed changes

oxinabox mentioned this pull request Dec 1, 2020

Tighten up a bunch of test accurasy requirements now that have FiniteDifferences v0.11.4 JuliaDiff/ChainRules.jl#315

Merged

oxinabox approved these changes Dec 1, 2020

View reviewed changes

Add and test estimate_roundoff_error

370bb16

oxinabox mentioned this pull request Dec 2, 2020

Unrevert "Tighten up a bunch of test accurasy requirements now that have FiniteDifferences v0.11.4 ) JuliaDiff/ChainRules.jl#317

Merged

Merge branch 'master' into wb/rework_add_tiny

08334bc

Bump patch number

be7668e

wesselb merged commit cfef6f8 into master Dec 3, 2020

giordano mentioned this pull request Dec 16, 2020

Accuracy is a bit brittle #124

Closed

oxinabox mentioned this pull request Mar 23, 2021

WIP: Rework and Simplify Internals #109

Merged

1 task

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

WIP: Remove `add_tiny` and add `estimate_magitude` #125

WIP: Remove `add_tiny` and add `estimate_magitude` #125

wesselb commented Nov 30, 2020

willtebbutt left a comment

wesselb commented Nov 30, 2020

willtebbutt commented Nov 30, 2020

willtebbutt left a comment •

edited

Loading

wesselb commented Nov 30, 2020

willtebbutt commented Nov 30, 2020

giordano commented Nov 30, 2020

giordano commented Nov 30, 2020

giordano commented Nov 30, 2020

wesselb commented Nov 30, 2020

giordano commented Dec 1, 2020

wesselb commented Dec 1, 2020

oxinabox Dec 1, 2020

wesselb Dec 1, 2020

oxinabox Dec 1, 2020

wesselb Dec 1, 2020

oxinabox commented Dec 1, 2020

oxinabox commented Dec 2, 2020 •

edited

Loading

wesselb commented Dec 3, 2020

oxinabox commented Dec 3, 2020

	return float(maximum(abs, f(x + Δ)))
	return estimate_magitude(f, x + Δ)

WIP: Remove add_tiny and add estimate_magitude #125

WIP: Remove add_tiny and add estimate_magitude #125

Conversation

wesselb commented Nov 30, 2020

willtebbutt left a comment

Choose a reason for hiding this comment

wesselb commented Nov 30, 2020

willtebbutt commented Nov 30, 2020

willtebbutt left a comment • edited Loading

Choose a reason for hiding this comment

wesselb commented Nov 30, 2020

willtebbutt commented Nov 30, 2020

giordano commented Nov 30, 2020

giordano commented Nov 30, 2020

giordano commented Nov 30, 2020

wesselb commented Nov 30, 2020

giordano commented Dec 1, 2020

wesselb commented Dec 1, 2020

oxinabox Dec 1, 2020

Choose a reason for hiding this comment

wesselb Dec 1, 2020

Choose a reason for hiding this comment

oxinabox Dec 1, 2020

Choose a reason for hiding this comment

wesselb Dec 1, 2020

Choose a reason for hiding this comment

oxinabox commented Dec 1, 2020

oxinabox commented Dec 2, 2020 • edited Loading

wesselb commented Dec 3, 2020

oxinabox commented Dec 3, 2020

WIP: Remove `add_tiny` and add `estimate_magitude` #125

WIP: Remove `add_tiny` and add `estimate_magitude` #125

willtebbutt left a comment •

edited

Loading

oxinabox commented Dec 2, 2020 •

edited

Loading