Rigorous approximation of
Download the notebook for this tutorial.
Introduction
Following the book Validated Numerics (Princeton, 2011) by Warwick Tucker, we find rigorous (i.e., guaranteed, or validated) bounds on using interval arithmetic, via the IntervalArithmetic.jl
package. First, let's import the package
using IntervalArithmetic
Calculating via a simple sum
There are many ways to calculate . For illustrative purposes, we will use the following sum
It is known that the exact value is . Thus, if we can calculate rigorous bounds on , then we can find rigorous bounds on .
The idea is to split up into two parts, , with and .
By bounding using integrals from below and above, we can see that . Rigorous bounds on are found by calculating it using interval arithmetic.
may be calculated by summing either forwards or backwards. A naive (non-interval) version could be the following:
function forward_sum_naive(N)
S_N = 0.0
for i in 1:N
S_N += 1. / (i^2)
end
S_N
end
S = forward_sum_naive(10000)
err = abs(S - pi^2/6.) # error
@show S, err
(S, err) = (1.6448340718480652, 9.999500016122376e-5)
Interval method
To find rigorous bounds for , we use interval arithmetic: each term is enclosed in an interval that is guaranteed to contain the true real value. A first idea is simply to wrap each term in the @interval macro, which converts its arguments into containing intervals:
function forward_S_N(N)
S_N = @interval(0.0)
for i in 1:N
S_N += @interval( 1. /(i^2) )
end
S_N
end
N = 10^5
@time rigorous_approx_S_N = forward_S_N(N)
2.710170 seconds (16.48 M allocations: 705.208 MiB, 5.64% gc time)
[1.64492, 1.64493]
We incorporate the respective bound on to obtain the bounds on , and hence on . We can also optimize the code by creating the interval @interval(1)
only once:
function forward_sum(N)
S_N = @interval(0.0)
interval_one = @interval(1.)
for i in 1:N
S_N += interval_one / (i^2)
end
T_N = interval_one / @interval(N, N+1)
S = S_N + T_N
sqrt(6*S)
end
N = 10^6
@time S = forward_sum(N)
@show S, diam(S)
0.657669 seconds (5.45 M allocations: 178.428 MiB, 7.03% gc time)
(S, diam(S)) = ([3.14159, 3.1416], 2.1298829366855898e-10)
We can ask for the midpoint–radius representation, which shows that the calculated bounds are correct to around 10 decimal places:
midpoint_radius(S)
(3.14159265358984, 1.0649436887888442e-10)
We may check that the true value of is indeed contained in the interval:
big(π) ∈ S
true
We may repeat the calculation, but now summing in the opposite direction. Due to the way that floating-point arithmetic works, this gives a more accurate answer.
function reverse_sum(N)
S_N = @interval(0.0)
interval_one = @interval(1.)
for i in N:-1:1
S_N += interval_one / (i^2)
end
T_N = interval_one / @interval(N, N+1)
S = S_N + T_N
sqrt(6*S)
end
N = 10^6
@time S = reverse_sum(N)
@show S, diam(S)
0.169420 seconds (3.03 M allocations: 77.610 MiB, 11.94% gc time)
(S, diam(S)) = ([3.14159, 3.1416], 9.57456336436735e-13)
Note that the sqrt function is guaranteed (by the IEEE 754 standard) to give correctly-rounded results, so the resulting bounds on are guaranteed to be correct.
Note also that due to the way that the IntervalArithmetic package works, we can make the code simpler, at the expense of some performance. Only two explicit calls to the @interval
macro are now required:
function reverse_sum2(N)
S_N = @interval(0.0)
for i in N:-1:1
S_N += 1 / (i^2)
end
T_N = 1 / @interval(N, N+1)
S = S_N + T_N
sqrt(6*S)
end
N = 10^6
@time S = reverse_sum2(N)
@show S, diam(S)
0.137809 seconds (2.06 M allocations: 48.472 MiB, 15.39% gc time)
(S, diam(S)) = ([3.14159, 3.1416], 9.57456336436735e-13)