Saturday, December 18, 2010

SICP 2.12 - 2.16: Extended Exercise: Interval Arithmetic (Part 2)

From SICP section 2.1.4 Extended Exercise: Interval Arithmetic

The solutions in the first half of this extended exercise, SICP 2.7 - 2.11: Extended Exercise: Interval Arithmetic (Part 1), are only one possible way to approach the problem of implementing interval arithmetic. Another way to represent intervals is as a center value and an added tolerance. The following alternate constructor and selectors are supplied in the text:
; Represent intervals as a center value and a width
(define (make-center-width c w)
(make-interval (- c w) (+ c w)))

(define (center i)
(/ (+ (lower-bound i) (upper-bound i)) 2))

(define (width i)
(/ (- (upper-bound i) (lower-bound i)) 2))

Here's a test case to show how it works:
> (define a (make-center-width 5 1))
> a
(4 . 6)
> (center a)
5
> (width a)
1

A third way of representing intervals, often used by engineers, is as a center value and a percent tolerance measured as the ratio of the width of the interval to its center value.



Exercise 2.12 asks us to define a constructor make-center-percent that takes a center and a percent tolerance and creates the desired interval. We'll also need to define a selector percent that extracts the percentage tolerance for a given interval.

We can define our new constructor in terms of make-center-width. We'll make our procedure take a whole number percentage p.
(define (make-center-percent c p)
(make-center-width c (* c (/ p 100.0))))

Since the percent tolerance of an interval is the ratio of the interval's width to its center, we can implement percent in terms of the provided width and center procedures.
(define (percent i)
(* 100.0 (/ (width i) (center i))))

Here's how we'd define the same interval as before with the new constructor:
> (define a (make-center-percent 5 20))
> a
(4.0 . 6.0)
> (center a)
5.0
> (width a)
1.0
> (percent a)
20.0



Exercise 2.13 asks us to show that there is, for small tolerances, a simple formula for the approximate percentage tolerance of the product of two intervals in terms of the tolerances of the factors. We can assume that all numbers are positive.

Recall the old definition of interval multiplication that we used was:

[a,b] × [c,d] = [min (ac, ad, bc, bd), max (ac, ad, bc, bd)]

Since we're assuming all positive numbers for this exercise, we can change the definition to a simpler form:

[a,b] × [c,d] = [ac, bd]

Now we need to complicate things again and look at what it means to represent an interval in terms of its center (ci) and percent tolerance (pi).

i = [ci - ci(pi/100), ci + ci(pi/100)]
i = [ci(1 - pi/100), ci(1 + pi/100)]

Multiplying two intervals x and y together would be:

xy = [cxcy(1 - px/100)(1 - py/100), cxcy(1 + px/100)(1 + py/100)]

We can use the FOIL method to combine some of the terms above.

(1 - px/100)(1 - py/100) = (1 - py/100 - px/100 + pxpy/10,000)
= 1 - (px + py)/100 + pxpy/10,000

(1 + px/100)(1 + py/100) = 1 + py/100 + px/100 + pxpy/10,000
= 1 + (px + py)/100 + pxpy/10,000

Inserting those back into the interval notation we get:

xy = [cxcy(1 - (px + py)/100 + pxpy/10,000), cxcy(1 + (px + py)/100 + pxpy/10,000)]

Since the exercise lets us assume that both px and py are small, and because we're only looking for an approximation, we can ignore the pxpy/10,000 terms because they'll be very small.

xy = [cxcy(1 - (px + py)/100), cxcy(1 + (px + py)/100)]

Now we have things back in terms of an interval's center and percent tolerance. The center of the product of the two intervals is cxcy, and the percent tolerance is (px + py). The approximate percentage tolerance of the product of the two intervals is the sum of the tolerances of the two factors.

We can check this with a quick example.
> (define a (make-center-percent 5 2))
> (define b (make-center-percent 10 3))
> (define c (mul-interval a b))
> (percent c)
4.997001798920647



Electrical resistor values are often expressed as a center value and a percent tolerance, so an engineer writing programs that work with resistances would be likely to make use of an interval arithmetic library. The formula for parallel resistors can be written in two algebraically equivalent ways:

and
The following two programs implement the two different computations:
(define (par1 r1 r2)
(div-interval (mul-interval r1 r2)
(add-interval r1 r2)))

(define (par2 r1 r2)
(let ((one (make-interval 1 1)))
(div-interval one
(add-interval (div-interval one r1)
(div-interval one r2)))))


Exercise 2.14 asks us to demonstrate that the two programs above give different answers for the same inputs by investigating the behavior of the system on a variety of arithmetic expressions.

First we'll make some intervals and look at the results of interval division in terms of center and percent.
> (define a (make-center-percent 100 5))
> (define b (make-center-percent 200 2))
> (define aa (div-interval a a))
> aa
(0.9047619047619049 . 1.1052631578947367)
> (define ab (div-interval a b))
> ab
(0.46568627450980393 . 0.5357142857142857)
> (center aa)
1.0050125313283207
> (center ab)
0.5007002801120448
> (percent aa)
9.97506234413964
> (percent ab)
6.993006993006991

What's notable here is that the center of the result of dividing an interval by itself is not 1, but just an approximation to it. This will be important when we explain what's wrong with our interval system in the next exercise. For now we're just trying to show that something isn't right.

We can use the same intervals that we defined above to illustrate the difference between the two parallel resistance procedures.
> (define apb1 (par1 a b))
> (define apb2 (par2 a b))
> apb1
(60.25889967637541 . 73.6082474226804)
> apb2
(63.986254295532646 . 69.32038834951456)

> (define apa1 (par1 a a))
> (define apa2 (par2 a a))
> apa1
(42.97619047619048 . 58.026315789473685)
> apa2
(47.5 . 52.49999999999999)

This verifies that there is a significant difference between the two procedures.


Exercise 2.15 points out that a formula to compute with intervals using the library we've been developing will produce tighter error bounds if it can be written in such a form that no variable that represents an uncertain number (an interval) is repeated. The conclusion is that par2 is a "better" program for parallel resistances than par1. We need to explain if this is correct, and why.

First, this is a good time to show that the two formulas used to develop par1 and par2 are algebraically equivalent. We're trying to show that


We'll start with the formula on the right hand side and derive the formula on the left. To do so, all we really need to do is multiply by R1/R1 and R2/R2. Since both of these fractions equal 1, these are valid transformations.


So the formulas are algebraically equivalent, but they don't give the same answer. Why could that be? The answer lies in the trick we used just now to show equivalence. We used the ratios R1/R1 and R2/R2 to change the formula and said that it was okay because that's just like multiplying by 1. But R1 and R2 represent resistor values, which are intervals, and we saw in exercise 2.14 that dividing an interval by itself doesn't equal 1, it just approximates it. Transforming the equation in this way introduces error. That's why the observation that we can get tighter error bounds if we avoid repeating variables that represent uncertain numbers is correct.



Exercise 2.16 asks us to explain why equivalent algebraic expressions may lead to different answers. It goes on to ask if we can devise an interval-arithmetic package that does not have this shortcoming, or if the task is impossible.

As the text points out, this is a very difficult problem. Wikipedia to the rescue.

The so-called dependency problem is a major obstacle to the application of interval arithmetic. Although interval methods can determine the range of elementary arithmetic operations and functions very accurately, this is not always true with more complicated functions. If an interval occurs several times in a calculation using parameters, and each occurrence is taken independently then this can lead to an unwanted expansion of the resulting intervals.

...

In general, it can be shown that the exact range of values can be achieved, if each variable appears only once. However, not every function can be rewritten this way.


In short, no we cannot design an interval arithmetic package that does not have this shortcoming in the general case. The best we can do, as was indicated in the previous exercise, is to try and write formulas that avoid repeating variables that represent intervals. This is not always possible.


Related:

For links to all of the SICP lecture notes and exercises that I've done so far, see The SICP Challenge.

3 comments:

gavri said...

I'm sure I'm missing something here, but I can't make the connection between the floating point arithmetic explanations for exercises 2.14, 21.5 and the Wikpedia explanation for exercise 2.16 which seems to be based on how once intervals don't get smaller in the progression of a calculation :-/

For example, with
(def x (make-interval 3 4)

(sub-interval (add-interval x x) x)

gives

(cons 2 5)

instead of (cons 3 4), which is what we'd get if the terms were cancelled before the values were substituted.

freyus said...

Bill, I initially arrived at the same answer for 2.15, but I don't think it's correct. Here's why...

Consider if you start with the first expression:

(1) R1*R2/(R1+R2)

and manipulate it to get the second expression:

(2) 1/(1/R1 + 1/R2)

We get the second expression by multiplying the first expression by (1/(R1*R2))/(1/(R1*R2)). Although this quantity is, algebraically, equivalent to 1, it has an associated error tolerance. Taking the product will *increase* in the error tolerance of the result. By the same logic expressed in your post, (2) must have a higher error tolerance than (1)!

I think the answer relates to the fact that any multiplication or division of two intervals causes an increase in total error tolerances. A reduction in how many times a variable appears should be accompanied by a reduction in multiplies and divides, and thus a lower error tolerance.

Unknown said...

For 2.16, we figured out that the general case of the problem is equivalent to asking the following question:

For a given expression, can it be reduced so that every interval term appears only once?

Which can be changed to the statement:

This expression is in its 'most-reduced' (or canonical, or normal) form.

Which, seems to be related to the Entscheidungsproblem (https://en.wikipedia.org/wiki/Entscheidungsproblem), as we're making a formal statement that we must prove true or false, and this was formally proved undecidable by Alonzo Church and Alan Turing. (see https://en.wikipedia.org/wiki/Church%E2%80%93Turing_thesis)

It is therefore impossible to write an interval arithmetic package that is generally capable of resolving the dependency problem.

I was really surprised how deep that question was.