IntervalRootFinding.jl
IntervalRootFinding.jl copied to clipboard
`ForwardDiff` cannot guarantee valid enclosures of derivatives
Unfortunately, computing the derivative with ForwardDiff.derivative
over an interval cannot guarantee a rigorous enclosure of the true derivative, because internally ForwardDiff
might still use some non-rigorous floating computations somewhere. Here is an example:
Consider the function $f(x) = 2^x$, whose derivative is $f'(x)=\ln(2)2^x$ and hence $f'(0) = \ln(2)$. Now with ForwardDiff
julia> using IntervalArithmetic, ForwardDiff
julia> f(x) = 2^x
f (generic function with 1 method)
julia> df = ForwardDiff.derivative(f, 0..0)
[0.693147, 0.693148]
julia> log(big"2")
0.6931471805599453094172321214581765680755001343602552541206800094933936219696955
julia> log(big"2") ∈ df
false
julia> diam(df)
0.0
The reason for this can be found in the source code of ForwardDiff
here, basically log(2)
is computed with floating point arithmetic and hence the true natural logarithm is not included
I am not sure how much this mines the newton method (or if it's already taken into account in this package)
this could be easily solved by
julia> ff(x) = Interval(2) ^ x
ff (generic function with 1 method)
julia> ForwardDiff.derivative(ff, 0..0)
[0.693147, 0.693148]
julia> diam(ans)
1.1102230246251565e-16
julia> log(big"2") ∈ ForwardDiff.derivative(ff, 0..0)
true
I guess documenting that to get verified results you need to have intervals everywhere in your function definition could be enough
Interesting; I guess there's a missing promotion in the ForwardDiff code?
Probably there should be a promotion here
Fixed on IntervalArithmetic v0.22.8:
julia> f(x) = interval(2)^x # Need to use interval to be guaranteed
f (generic function with 1 method)
julia> df = ForwardDiff.derivative(f, interval(0, 0))
[0.693147, 0.693148]_com
julia> in_interval(log(big"2"), df)
true
Of course it is only fixed for power and if you use interval for everything. Otherwise, the best we can do is to tell that the computation is not guaranteed, which also currently works.
julia> g(x) = 2^x # Mixing interval with numbers
g (generic function with 1 method)
julia> dg = ForwardDiff.derivative(g, interval(0, 0)) # This has the Not Guaranteed tag, telling we should not trust the result
[0.693147, 0.693148]_com_NG
julia> in_interval(log(big"2"), dg)
false
@OlivierHnt This is a good example of the usefulness of the NG tag !