[ Marquette | MSCS | Corliss | Vladik's Intervals ]
Interval FAQ: [ Entry page | Contents | Search ]
About Sun's f95 with Interval Support

# Interval FAQ How do I convert a real/double code to intervals?

## How do I convert a real/double code to intervals?

OK, I've heard that validated computations using interval arithmetic are useful. I have this old [ Fortran | C | ... ] [real | double precision | float | double | ...] code. How can I easily convert it to intervals?

See related FAQ: Can I replace real arithmetic with intervals?

### Response by George Corliss:

(georgec@mscs.mu.edu), Department of Mathematics, Statistics, and Computer Science, Marquette University.

### Introduction

Welcome. I'm glad you are interested. I find the guarantees of validated computations useful for computations with

• Uncertain parameter values,
• Questionable accuracy,
• VERY important consequences,
• Existence of a solution is in doubt,
• Uniqueness of a solution is in doubt.

The advice offered here may put you off. I am attempting to be honest. Converting legacy code to interval arithmetic is somewhat like converting it for parallelism or for automatic differentiation. It can be done. Tools and compilers can help, but there are pitfalls for the unwary. My purpose is to offer guidance through some potential pitfalls.

Do you want to learn to write interval code, or do you want to get a validated result? In the first case, start with something very small to learn the basics. In the second case, look for existing software to see how you express your problem to it. Then, start with something very small, and grow from there. I do not recommend that your first interval code be hundreds of lines long in any case.

Outline: [ Introduction | Mechanics | Containment | Performance | Examples | Disclaimer | Readings ]

### Mechanics

#### Environment

You will need some environment or processor supporting interval arithmetic. See IFAQ - "Interval compilers and programming environments?"

#### Standard Code

Some questions to ask:

• Do you have the source? If not, you are out of luck.
• Is the code 100% standard-conforming? If not, fix it before attempting conversion.
• How old is it? Older codes tend to be more challenging.
• How long is it? Longer codes tend to be more challenging. I recommend your first attempt not more than a few tens of lines.
• Did you write it? Revising someone else's code is more challenging.
But you already knew all those hints, didn't you?

#### Change Types

In general, you want to change types real, double precision, float, double, etc. to INTERVAL as supported by the programming environment you select. However, in most programs, not all variables are re-declared. In concept:

Before:

DOUBLE PRECISION A, X, Y(20)

After:

DOUBLE PRECISION A
INTERVAL X, Y(20)

Of course, the exact syntax depends on the programming environment.

Advice: Consider carefully which variable should be intervals.

Outline: [ Introduction | Mechanics | Containment | Performance | Examples | Disclaimer | Readings ]

### Containment

The fundamental demand of interval arithmetic is containment: "Thou shalt not lie."

In order to guarantee containment in legacy code, there are a few aspects that you must consider. Mechanical transformations rarely give exactly what you intend. Each of these pitfalls can yield to inaccurate floating-point results, too.

#### Literals

Danger: Most real literals are not exact floating point numbers.

Consider the intent of every floating point literal. Most real numbers are not exactly representable, so you should consider explicit conversion to interval. In concept:

Before:

A = 0.1

After:

A = interval (0.1, 0.1)

Alternatives:

1. Literal is an approximation. An interval literal can capture the uncertainty. You want a relatively thick interval, e.g. interval (0.092, 0.16).
2. Literal is exact, but it is rounded to an internal floating point representation. An interval encloses the true value. You want an interval only one or two ULPs (Unit in the Last Place) wide, e.g. interval(0.1, 0.1). Beware: The literal "0.1" may be rounded by the compiler before the run-time function interval() gets a chance.
3. The floating point representation of the literal is exact. Beware: This is rarely what you really intend, but it is usually what you get if you leave the floating point literal in the code.
4. The literal is exact. Beware: That is probably not true in floating point arithmetic, either.

Of course, it is easier to convert literals when they are declared as constants, rather than appearing as "magic numbers" in the program text. Read carefully how your programming environment handles floating point and interval literals.

Advice: Look at, and think about, every floating point literal in your program text.

#### Input/Output of interval and floating point values

Danger: Real numbers entered or printed are rarely exact.

The fundamental demand of interval arithmetic, "Thou shalt not lie," is often at odds with our expectation that READ and WRITE are inverses (which is rarely true for floating point numbers, either).

Containment requires that READing yields an interval containing the interval represented by the character string entered.

Similarly, the characters produced by a WRITE must represent an interval containing the internal representation.

Most interval programming environments provide special facilities for I/O. Read the instructions carefully.

Advice: Look at, and think about, every I/O READ, WRITE, print, scan, etc. of floating point data in your program text and the input you provide.

#### IF and similar branching statements

Danger: The direct interval analogue of IF is ambiguous.

Suppose

```DOUBLE PRECISION X
IF (X <= 0.0) THEN
F = -1.0
ELSE
F = 1.0
END IF```

were changed into (it should not be)

```INTERVAL X
IF (X <= 0.0) THEN
F = -1.0
ELSE
F = 1.0
END IF```

When X = [-1.0, -0.8] or X = [0.8, 1.0], the meaning is clear. F = -1.0 or F = 1.0, respectively.

What should we do when X = [-0.1, 0.2]? We would need a three-valued logic: True, False, Maybe. The requirement of containment specifies F must contain [-1.0, 1.0]. For this simple example, we could imagine executing both branches and taking the interval hull, but any such automatic transformation is extremely difficult for realistic examples. The result is probably not what you intend, anyway.

Instead, several interval programming environments provide "certainly" and "possibly" logical operators. Other interval environments interpret logical operators as either "certainly" or "possibly."

"Certainly" means that the operator is True for each pair of elements, e.g.

• X is certainly less than Y if x < y, for every x in X and Y in Y, ie. Supremum(X) < Infimum(Y)
• X is certainly equal to Y if x = y, for every x in X and Y in Y, ie. Infimum(X) = Supremum(X) = Infimum(Y) = Supremum(Y)

"Possibly" means that the operator is True for some pair of elements, e.g.

• X is possibly less than Y if x < y, for some x in X and Y in Y, ie. Infimum(X) < Supremum(Y)
• X is possibly equal to Y if x = y, for some x in X and Y in Y, ie. X intersect Y is not empty

The distinction of "certainly" and "possibly" applies only to logical operations involving intervals.

Determine how your interval programming environment handles logical operations. You MUST consider every IF or SWITCH statement involving floating point quantities. Alternatives:

1. "Certainly" logical operator
2. "Possibly" logical operator
3. Do-it-yourself by explicitly using endpoints
4. Leave, if operands remain as floating point types

In general, Alternative 3 is not as attractive as it might appear, and you cannot get away with Alternative 4 as often as you might hope.

Advice: Express each logical operation on floating point objects as "certainly" or "possibly" logical operators.

#### What is contained?

Danger: A guaranteed enclosure may not enclose what you think it does.

We interval folks are fond of promoting the validated enclosures provided by interval algorithms. Automatic differentiation provides a derivative of what your code does, not what you may think it does. Similarly, an intervalization of an algorithm computes an enclosure of what the original code should compute, but that may not be what you think it computes.

For example, a Runge-Kutta code evaluated in interval arithmetic computes a guaranteed enclosure of the approximation the original code should compute. That may be useful for tracking possible lost accuracy in the algorithm. However, it is not an enclosure of the solution to the differential equation being solved, and certainly not an enclosure of the original problem being modeled, because it does not enclose truncation errors.

In general, good floating-point algorithms make poor interval algorithms. If you are solving a standard mathematical problem, consider looking for suitable existing interval software.

Advice: Think carefully about what you have enclosed. Consider existing interval software.

Outline: [ Introduction | Mechanics | Containment | Performance | Examples | Disclaimer | Readings ]

### Performance -- Tightness

The requirement for an interval environment is containment; they compete on the basis of tightness, speed, and ease of use. It is in the area of achieving tightness where the most subtle pitfalls await.

#### Expression rearrangement

Danger: Mathematically equivalent expressions are nor necessarily equivalent in interval arithmetic.

Hint: Mathematically equivalent expressions are not necessarily equivalent in floating-point arithmetic, either,

The fundamental theorem of interval arithmetic asserts that an evaluation in interval arithmetic of expressions or functions most commonly used in scientific computing yields an enclosure of the true value. However, sometimes a mathematically equivalent expression can yield a much tighter enclosure.

For example, f(x) = 2x3 - 3 x2 + 1 evaluated in on X = [1, 2] yields [-4, 14], while the equivalent expression x2 (2x - 3) + 1 yields [-3, 5]. Tight bounds are [0 = f(1), 5 = f(2)], since f is monotonic.

In general, expressions evaluate more tightly when each variable appears as few times as possible.

If tighter evaluations are necessary, consider centered forms, Taylor forms, or other more sophisticated tools. It is often the case that relatively costly strategies for computing tight bounds pay off in overall algorithmic efficiency because the underlying algorithm works better. However, such more sophisticated forms of function bounding are not interval versions of real code.

Advice: A little mathematics beats a lot of computing. Simplify expressions.

#### Good algorithms

Danger: Interval versions of good floating-point algorithms are rarely good interval algorithms.

An interval version of a floating-point algorithm does yield enclosures of the results computed by the floating-point algorithm. That is often useful. However, the enclosing intervals are sometimes too wide to provide insight. In that case, algorithms designed for validation are called for.

For example, variants of Gaussian elimination work well in practice for the approximate solution of dense linear systems. However, Wilkinson showed that the worst-case accumulation of round-off errors could be very bad. Interval arithmetic must enclose the worst-case round-off error, so an interval arithmetic version of Gaussian elimination sometimes yields very pessimistic bounds. Good interval linear equation solvers use residual correction or another contractive iteration.

Good interval algorithms usually use one or more of these general strategies:

• Contractive iteration
• Intersection of several answers
• Approximation, followed by validation
• Enclosure of truncation errors

Also, as pointed out above, the result of an interval version of a floating-point algorithm probably does not enclose the true answer to the mathematical problem posed.

Advice: Think carefully about what you are enclosing. Use existing quality software whenever possible.

Outline: [ Introduction | Mechanics | Containment | Performance | Examples | Disclaimer | Readings ]

### Examples

Yes, I should offer some examples. In truth, I have only trivial examples of real code converted to interval code because I am not in the habit of converting legacy codes. I write my interval codes from scratch. However, if you have an example to share, please send it to me (georgec@mscs.mu.edu), and I'll put it here.

Outline: [ Introduction | Mechanics | Containment | Performance | Examples | Disclaimer | Readings ]

### Disclaimer

In this brief treatment, I cannot address all possible pitfalls. Conversion of legacy code to intervals is probably more challenging than converting the same code to parallel processing or for automatic differentiation. It is feasible, and there are rewards, but it is not a mechanical process.

Outline: [ Introduction | Mechanics | Containment | Performance | Examples | Disclaimer | Readings ]

See IFAQ - "What books should I start with?"

Outline: [ Introduction | Mechanics | Containment | Performance | Examples | Disclaimer | Readings ]

[ Marquette | MSCS | Corliss | Vladik's Intervals ]
Interval FAQ: [ Entry page | Contents | Search ]
About Sun's f95 with Interval Support

If you have a question related to validated computing, interval analysis, or related matters, I recommend

Access count since 20 Apr 1999 : 