Verified interval arithmetic in pure R

What this package guarantees, and what it does not

Every operation in this package returns an enclosure: an interval that provably contains the exact real result. Not an estimate of it, not an error bar around it. If the true answer is a real number y, then y lies inside the returned interval, and that statement is a theorem about the code rather than an observation about a sample of inputs.

The guarantee is deliberately narrow. An enclosure that contains the answer may still be wide, and a wide enclosure is honest but not useful. Throughout the package the two properties are kept apart and reported apart: validity is never traded for tightness, and where a result cannot be tightened the package says so instead of narrowing it.

Conformance with IEEE Std 1788.1-2017 is not claimed, for the reason the standard itself gives: its subclause 1.5 makes conformance a list of requirements an implementation shall satisfy, with no partial grade available. What is followed, measured one requirement at a time, is documented in ?RobustArithmetic and summarised at the end of this vignette.

The kernel: intervals and outward rounding

An interval is built from its endpoints and carries a decoration that records what is known about the operation that produced it.

x <- ra_interval(1, 2)
y <- ra_interval(3, 4)
ra_add(x, y)
#> <ra_ivl[1]>
#> [3.999999999999999, 6.000000000000001]_com
ra_mul(x, y)
#> <ra_ivl[1]>
#> [2.999999999999999, 8.000000000000002]_com

The endpoints of a computed interval are rounded outward, so the enclosure can only grow and never silently shrink past the true result. R offers no way to change the processor’s rounding mode, and changing it would not be a local act in any case: it is per-thread state of the processor, so it would govern every floating-point operation executed afterwards on that thread, inside this package or anywhere else. Outward rounding is therefore obtained arithmetically, from the predecessor and successor formulas of Rump, Zimmermann, Boldo and Melquiond (2009), which are valid under round-to-nearest.

ra_pred(1)
#> [1] 1
ra_succ(1)
#> [1] 1

Division by an interval containing zero is where a naive implementation quietly returns nonsense. Here it returns the mathematically correct answer, which is a pair of intervals, through the extended division of the interval Newton literature.

ra_div_extended(ra_interval(1, 2), ra_interval(-1, 1))
#> $first
#> <ra_ivl[1]>
#> [-Inf, -0.9999999999999998]_trv
#> 
#> $second
#> <ra_ivl[1]>
#> [0.9999999999999998, Inf]_trv

Elementary functions at two levels

Elementary functions are provided at two levels, and the difference between them is a difference in what can be proved rather than a difference in speed alone.

The fast level evaluates over the system math library and widens the result by a declared slack derived from published accuracy measurements. The rigorous level evaluates over Rmpfr with a directed-rounding bridge. When a verdict would otherwise fall inside the slack — that is, when the fast level cannot decide the question being asked — the package escalates through a ladder of precisions rather than returning a confident answer it cannot defend.

ra_elem("exp", ra_interval(0, 1))
#> <ra_ivl[1]>
#> [0.9999999999999996, 2.718281828459047]_com*
#> * measured provenance: fast-level slack, a declared convention over a measured error, not a theorem
ra_has_mpfr()

The slack is not a tuning knob chosen for convenience. It is measured against the library actually installed, which matters more than it sounds: R does not call the system libm along a single fixed path, so a bound published for one library does not certify what R computes. ra_measure_library_error() is the instrument that measures it.

Expressions, and why the centered form exists

Evaluating an expression by substituting an interval for each occurrence of the variable — the natural extension — is valid but often loose, because each occurrence is treated as if it varied independently of the others. The centered extension trades that looseness for a derivative bound.

ra_eval_natural(quote(x^2 - x), list(x = ra_interval(0, 1)))
#> <ra_ivl[1]>
#> [-1.0000000000000003, 1.0000000000000005]_com
ra_enclose_expr(quote(x^2 - x), ra_interval(0, 1))
#> <ra_ivl[1]>
#> [-0.7500000000000008, 0.2500000000000008]_com

The natural extension gives [-1, 1]; ra_enclose_expr() gives [-0.75, 0.25], which is the exact range. Neither form is uniformly better and the package does not pretend otherwise, so ra_enclose_expr() takes the best of the methods it is given —- natural, mean value and a monotonicity test -— rather than committing to one.

Roots: the interval Newton operator

ra_solve() finds roots of a univariate expression over a box and returns, for each one, a word about what was established. The engine is the Hansen-Sengupta operator with extended division and epsilon-inflated candidate verification, driven by a subdivision (paving) loop.

ra_solve(quote(x^2 - 2), ra_interval(0, 2))
#> <ra_paving> roots of x^2 - 2 over [0, 2]_com
#>   unique roots: 1  (provenance: theorem)
#>     [1.4142135623730947, 1.4142135623730956]_trv  unique
#>   not excludable: 0
#>   absence demonstrated over the rest of the window (0 excluded boxes, re-verified rigorously)
#>   budget: 2 of 10000 boxes

The important property is the failure mode. The loop has a budget, and when the budget runs out it abstains by name and prints what it spent, rather than returning a smaller set of roots as though the search had been exhaustive. An abstention is information; a silently truncated answer is not.

Fixed points in higher dimensions, without interval matrices

The solver above is univariate by contract. For fixed points of maps on R^n the package takes a different route, and the route is the point of the design: the whole certificate is a scalar inequality, so a fixed point in R^n is certified without interval matrices, without factorisations and without an interval spectral decomposition. The norm collapses the dimension into three numbers, and the caller may bound those three numbers by any means they can defend.

ra_ball_certificate(rho = 0.2, lipschitz = 0.5)
#> fixed point certificate over a ball in the unspecified norm
#>   confinement : certified by the lipschitz route, radius 0.40000000000000013
#>   contraction : certified, the fixed point is unique in the ball
#>   the centre is within 0.4000000000000002 of it
#>   bounds used : rho 0.20000000000000002, L 0.5, image not given
#>   provenance  : theorem

Two theorems are applied and the object says which one decided the verdict. Confinement — the map carries the ball into itself — gives existence by Brouwer and says nothing about what happens inside. Contraction adds uniqueness and convergence by Banach. They are certified separately because they are not the same statement, and a package that reported one word for both would be reporting less than it knows.

ra_ball_certificate(rho = 0.24, lipschitz = 1.6, radius = 0.4, image_radius = 0.4)
#> fixed point certificate over a ball in the unspecified norm
#>   confinement : certified by the image route, radius 0.4
#>   contraction : not certified, so what happens inside is not claimed
#>   reason      : L = 1.6000000000000001 is not below one, so there is no contraction.
#>   bounds used : rho 0.24, L 1.6000000000000001, image 0.4000000000000001
#>   provenance  : theorem

That second call is confined but not contractive: the image radius route needs no Lipschitz constant, which is why it is not redundant.

Two assemblers help build the three inputs in the matrix case without factorising anything: ra_gershgorin() encloses eigenvalues from the matrix entries, and ra_spectral_sum() bounds a spectral quantity of a sum from bounds on its terms.

ra_gershgorin(matrix(c(2, 0.1, 0.1, 3), 2, 2))
#> <ra_ivl[1]>
#> [1.8999999999999996, 3.100000000000001]_com

The responsibility that cannot be delegated: the three bounds must be measured in the same norm, and the package cannot check that. The label travels with the object. Mixing norms silently would be a failure mode indistinguishable, at run time, from a valid certificate. Likewise the Lipschitz constant must hold over the whole ball and not merely at its centre; a caller who bounds the derivative at the centre receives a certificate that means nothing, and no arithmetic can detect it. In one dimension the package closes that loop itself, in ra_certify_fixed_point().

A refusal certifies nothing

This is worth stating on its own, because it is the most common way a verified result is misread. The hypotheses behind every certificate here are sufficient and not necessary. There is no test of exclusion in this package and no word of demonstrated absence. A map that this package declines to certify may still have a fixed point in the ball; a root that the solver does not report may still exist inside a box the budget never reached. What the package proves, it proves; what it declines, it declines without claiming the opposite.

What is covered, and what is not

Measured against IEEE Std 1788.1-2017, one requirement at a time: the interval type and the decoration system of clause 5; 22 of the 39 arithmetic operations of Table 4.1; and the seven numeric functions of Table 4.3. What is not provided is the cancellative operations, the interval comparison relations, the text input of subclause 6.8.2, the interchange representation of subclause 7.3, and the tightest accuracy that subclause 6.5.2 requires of the basic operations — which here are one unit in the last place wider at each end. The text output of subclause 6.8.3 is provided, by ra_interval_to_text() and by the printing of every interval: the pair of decimals on the page encloses the pair of doubles underneath, which is what that subclause asks and what rounding each endpoint to nearest would break.

ra_operator_table() prints the coverage as data rather than as prose, and ra_verify_operator_table() checks the claim against the code.