Interval arithmetic and Prolog

Arithmetic over intervals, known as interval arithmetic is useful. Its a way to include uncertainty in calculating arbitrary formulae without probability and with support for a complete set of linear and non-linear operations. For example, in vanilla interval arithmetic: A + B = C where A=[0,5], B=[5,10] then C = [5,15].

A fancier version which automates the treatment of correlations between variables is called affine arithmetic. Same kind of thing – computation over intervals for arbitrary functions – but takes variable dependency into account.

think that a well written extension like this could greatly improve the range of problems which Prolog is apt for. @ridgeworks I note this paper on your github, would you agree?

1 Like

I think the Julia package provides many strong motivating examples:

I have written a pack for this, because my students round intermediate results in their calculations, e.g. with

t = frac(M - mu, s / sqrt(n))

Example: numerator is 1.368, my students, “round” it to 1.36, denominator something similar, result is then also an interval. I use it for an e-learning system.

The pack needs some optimization, as well as some default behavior for functions that are know to be monotonically increasing or decreasing.

So, interval arithmetic (IA) has been available in Prolog systems for well over 30 years. CLP(BNR) referenced in that paper was one of the earliest but, e.g., there’s Eclipse and CLIP (Hickey et/ al.) and near-Prolog systems such as Newton and Numerica (which I believe got absorbed into IBM’s ILOG solver). The current clpbnr pack for SWIP is an extended reimplementation (entirely in Prolog) of earlier versions of CLP(BNR), but also available with SWIP is the inclpr package, a CLP over reals supporting nonlinear polynomials, which also is based on IA.

To me the main rationale for interval arithmetic is to support a CLP instance over the continuous domain of reals using IEEE floating point arithmetic (which includes infinities and rounding modes) in a mathematically correct way, i.e. free of rounding errors.

It would seem that there are already “well written extensions” for Prolog so I’m not quite sure what you’re looking for. The main issue (IMHO) is connecting the folks with these kinds of applications to the enabling technology, which would also seem to apply to Prolog/logic programming in general.

1 Like

Or if you’re just looking for IA, using clpBNR:

?- A::integer(0,5), B::integer(5,10), {A+B==C}.
A::integer(0, 5),
B::integer(5, 10),
C::real(5, 15).

but also, because it’s CLP:

?- {A+B==C}, A::integer(0,5), B::integer(5,10).
A::integer(0, 5),
B::integer(5, 10),
C::real(5, 15).
1 Like

Will it be able to calculate intervals for non trivial examples? E.g. cos(A) + B^2 / C = D?

EDIT. Yes, but it appears to be wide.

?- {A*A/B + cos(B)*A ==C}, A::integer(1,100), B::integer(1,100).                                                               
A::integer(1, 100),                                                                                                            
B::integer(1, 100),                                                                                                            
C::real(-99.99000000000001, 10100.0). 

Verifying:

?- aggregate_all((min(C),max(C)),
                (between(1,100,A), between(1,100,B),
                C is A*A/B + cos(B)*A), Result).

Result = (-22.493788489021814, 10054.030230586814).

A side point, I’m also confused why in your A+B example, adding two explicitly declared integers results in a real range? This seem to be false. There is no selection of A and B which will produce 10.5 for example.

Me too. So I tried

?- {A+B==C}, A::integer(0,5), B::integer(5,10), C::integer.
A::integer(0, 5),
B::integer(5, 10),
C::integer(5, 15).

I think it’s perfectly reasonable that an ‘undeclared’ variable assumes the most general domain, like the arithmetic operators do.

I easily confuse intervals with domains since I haven’t used this stuff much before . I think intervals just claim to bracket the variables , and the game is to shrink the bracket as much as possible. Meanwhile domains imply that there is a circumstance in which the variable could take every possible value in a domain.

There is a limited set of circumstances in which the two are equivalent.

Yes, this is a well known property of interval arithmetic - the narrowest bounds for an expression are only generated when each variable only occurs once in that expression. It’s often referred to as the dependency issue/problem and is fairly widely discussed in the academic literature. There’s a simple example in the clpbnr User Guide (courtesy F. Benhamou): the equation Y is X*X+1 has a mathematical equivalent (using rationals to avoid floating point nonsense) Y is (X + 1r2)**2 - 1r4. If X is an interval in the range (-1,1), evaluation of the two expressions yields different interval values - the second one is narrower because there is only one ocurrance of X:

?- X::real(-1,1),{Y1 is X*(X+1), Y2 is (X+1/2)**2 - 1/4}.
X::real(-1, 1),
Y1::real(-2, 2),
Y2::real(-0.25, 2).

But there’s a second reason why using interval arithmetic on this particular example is potentially sub-optimal; consider cos(B) where B::integer(1,100). If B was a continuous (real) interval, the bounds of cos(B) would be (-1.0,1.0). But because B is an integer in that same range it’s possible (probable?) that the domain is narrower, i.e., no integral value of B in that range would actually produce a -1.0 or 1.0. However there’s no way of discovering that without an exhaustive search through all the possible values of B (as in your use of aggregate_all), and that’s something no interval arithmetic library will likely do.

There are other practical reasons why domains in interval constraint system don’t narrow to optimal values (weak constraint propagators, accumulated rounding errors, etc.) but that’s a whole other discussion.

So when you’ve got a variable with an interval domain, it’s important to remember that not all values in the domain are solutions, or even that there are any solutions. The only guarantee is that there are no solutions outside the domain. Finding solutions is a separate step (searching or labelling), unless the domain has narrowed to a single point value.

The set of integers is a subset of reals so it’s not false, but perhaps not as “narrow” as you’d like. As noted above, because a value is an element of an interval domain does not mean it’s a solution. (If you tried to unify C with 10.5 it would fail.)

The operational semantics of type integer in clpBNR is that the bounds are always integers (calculated non-integer bounds will contract inward to the closest integer) and integer variables can only be unified with integer values in their domain.

That’s the conclusion I came to. I suppose it’s possible to infer a more constrained type of a variable from the type of the variables it’s constrained by, but I don’t see much point; the effect on the solution space is the same (i.e., solutions for C can only be integers).

3 Likes

Put this way, a domain is a set of actual solutions which doesn’t really exist until the set is known, i.e., it’s a somewhat abstract notion. Instead, I would say a domain is a set of possible solutions, and an interval is one way of specifying that set in a concrete way. But there are other ways of specifying a domain, e.g., symbolically (set expressions) or a union of intervals.

3 Likes

I’ve put the intervals on github, if someone wishes to contribute. A student will also add some stuff.

It’s just “functional”, i.e. cannot solve an equation for an unknown as clpbnr does. And it’s really just a stub, with only +, -, *, and /.

But it illustrates the three ways the problems are solved, + and - being monotonical in their arguments, * using combinations of the bounds of the intervals, and / non-deterministic.

1 Like

It would be interesting to gauge the demand for a functional interval arithmetic library. Ideally (IMO) it would be nice to integrate it into the standard Prolog functional arithmetic (is/3 et. al.) so you could write things like R is 0..5 + 2..pi resulting in R = 2..8.141592653589795.

I experimented with the idea of extending arithmetic types using booleans, strings, and arrays a while ago - see GitHub - ridgeworks/arithmetic_types: Extend SWIP arithmetic with user defined types. Doing something similar for intervals would entail a fair amount of work to support comparison (for ordering), additional functions, rounding modes(?), … .

The inclpr package includes an interval arithmetic module (inclpr_interval_arithmetic) but I think it’s intended for use with the package - bounds are aways floats and I wouldn’t call the API user friendly.

I’ll circle back around to it some time soon: I have a fetish for these lesser known high leverage concepts so far as they can be applied to my sort of garden variety problems :slight_smile:

Very great explanation.

This reminds me of an exercise from SICP, see 2.1.4 Extended Exercise: Interval Arithmetic.

The Exercise 2.13 - 2.16 asks readers to explain why this happens.

My answer was that, for Interval Arithmetic, addition preserves the tolerance, but multiplication and division add their tolerance together.

So given two variables, if we just do addition and subtraction, no problem.

However, if we do multiplication and division, then

  1. If the two variables are independent, no problem. Because adding tolerance together is necessary.
  2. If the two variables are dependent, then there is a problem. Because, for example A/A, it is definitely a single point 1, but for interval arithmetic would add their tolerance together.

Also, in Exercise 2.16, the authors asks an open question (I’m not sure if it’s currently solved today though):

Exercise 2.16: Explain, in general, why equivalent algebraic expressions may lead to different answers. Can you devise an interval-arithmetic package that does not have this shortcoming, or is this task impossible? (Warning: This problem is very difficult.)

My go-to reference for questions like this: “Introduction to INTERVAL ANALYSIS, Moore et. al, 2009” says:

We stress that two expressions can be equivalent in real arithmetic but not equivalent in interval arithmetic. This is due to the lack of distributivity and additive and multiplicative inverses in interval arithmetic.

(Explanations earlier in the text explain why this is so.)

A subsequent quote:

It turns out that any natural interval extension of a rational function in which each variable occurs only once (if at all) and to the first power only will compute the exact range of values provided that no division by an interval containing zero occurs.

This looks pretty fundamental to me (a non-expert).

1 Like

Thanks for sharing the great information!

Just want to clarify, do you think “to the first power” here refers to the first power for each variable, or does it refer to the 1 degree multivariate rational function?

For example, the condition:

each variable occurs only once (if at all) and to the first power

Do you think this condition hold for the expression X*Y*Z/W? (Assuming W interval containing no zero).

I think it holds, although the total degree is 3, but each variable is the first power. (Correct me if I am wrong).

Thanks.

I think it’s the former, i.e., X*Y*Z/W should optimally narrow. Another quote from the book:

The overestimation when we compute a bound on the range of X^2 as X*X is due to the phenomenon of interval dependency. Namely, if we assume x is an unknown number known to lie in the interval X, then, when we form the product x*x, the x in the second factor, although known only to lie in X must be the same as the x in the first factor, whereas, in the definition of the interval product X*X, it is assumed that the values in the first factor and the values in the second factor vary independently.
Interval dependency is a crucial consideration when using interval computations.

It’s a bit difficult to show this using clpbnr because constraints undergo a symbolic simplification which converts multiplications to powers if possible, since powers are supported as a separate primitive. But using a non-exported predicate:

?- X::real(-1,1), clpBNR:constrain_(Y1==X*X), clpBNR:constrain_(Y2==X**2).
X::real(-1, 1),
Y1::real(-1, 1),
Y2::real(0, 1).

So I think if interval powers are not decomposed to (repeated) multiplications of the same variable, the power restriction from the book doesn’t apply, i.e., X**2*Y*Z/W should also optimally narrow. (Note that the power function is anything but simple given mixed mode arithmetic.)

1 Like

For the record:

Here’s a small demo implementing interval addition using the clpfd range syntax:

:- module(type_interval,
	[
	 op(450, xfx, '..')  % same as clpfd def.
	]).

:- use_module(library(arithmetic_types)).

:- arithmetic_function('..'/2).
:- arithmetic_function('+'/2).

to_interval(N,N,N)    :- number(N), !.
to_interval(L..U,L,U) :- number(L), number(U), L =< U.

..(L,U,R) :- to_interval(R,L,U).  % point intervals will evaluate to number

+(L1..U1,I2,R) :- !,              % first or both interval arguments
	to_interval(I2,L2,U2),
	interval_add(L1,U1, L2,U2, R).
+(N1,L2..U2,R) :-                 % second interval argument 
	number(N1),                   % first must eval to a number
	interval_add(N1,N1, L2,U2, R).

interval_add(L1,U1, L2,U2, R) :-  % interval result
	UL is L1+L2, UR is U1+U2,
	to_interval(R,UL,UR).

so:

?- R is 1..2 + 3..4.
R = 4..6.

?- R is 1..2 + 3.5.
R = 4.5..5.5.

?- R is cos(0)..(pi+1) + 1r3..1r2.
R = 1.3333333333333333..4.641592653589793.

?- R is 1..(3-2).
R = 1.
1 Like