Arjen Markus (9 july 2003) Having read a book about numerical accuracy and stability, as well as several technical descriptions about interval arithmetic (certain Fortran compilers support it), I wanted to have this facility in Tcl too.
The idea: a number is not completely accurate, but there is some error margin. Then summing two such numbers becomes a matter of manipulating intervals:
x = [[a,b]] y = [[c,d]]
then:
x+y = [[a+c,b+d]]
for instance.
I know that it has severe drawbacks (and for these I want to provide some counter measures, like an alternative technique, called Running error analysis):
Well, creating the procs was easy and fast, so we can focus on obtaining experience :) Here is the script with a simple example.
Note: the functional or prefix style of arithmetic is really useful here - it will allow us to overload the procedures without changing the code everywhere.
# interval.tcl --
# Rudimentary support for interval arithmetic in Tcl
#
# To do:
# Elementary functions - cos, sin, tan, abs, exp, log, pow
# General function - function
#
namespace eval ::math::interval {
namespace export + - * / \
contains interval begin-end \
square abs
}
# + --
# Add two interval numbers
#
# Arguments:
# operand1 First operand (interval number or regular number)
# operand2 Second operand (interval number or regular number)
# Result:
# Sum of the two intervals
#
proc ::math::interval::+ {operand1 operand2} {
foreach {begin1 end1} [concat $operand1 $operand1] {break}
foreach {begin2 end2} [concat $operand2 $operand2] {break}
return [list [expr {$begin1+$begin2}] [expr {$end1+$end2}]]
}
# - --
# Subtract two interval numbers
#
# Arguments:
# operand1 First operand (interval number or regular number)
# operand2 Second operand (interval number or regular number)
# Result:
# Difference of the two intervals
#
proc ::math::interval::- {operand1 operand2} {
foreach {begin1 end1} [concat $operand1 $operand1] {break}
foreach {begin2 end2} [concat $operand2 $operand2] {break}
return [list [expr {$begin1-$begin2}] [expr {$end1-$end2}]]
}
# * --
# Multiply two interval numbers
#
# Arguments:
# operand1 First operand (interval number or regular number)
# operand2 Second operand (interval number or regular number)
# Result:
# Difference of the two intervals
#
proc ::math::interval::* {operand1 operand2} {
foreach {begin1 end1} [concat $operand1 $operand1] {break}
foreach {begin2 end2} [concat $operand2 $operand2] {break}
set v1 [expr {$begin1*$begin2}]
set v2 [expr {$begin1*$end2 }]
set v3 [expr {$end1 *$begin2}]
set v4 [expr {$end1 *$end2 }]
set min $v1
set max $v1
foreach v [list $v2 $v3 $v4] {
set min [expr { $min<$v? $min : $v }]
set max [expr { $max>$v? $max : $v }]
}
return [list $min $max]
}
# / --
# Divide two interval numbers
#
# Arguments:
# operand1 First operand (interval number or regular number)
# operand2 Second operand (interval number or regular number)
# Result:
# Difference of the two intervals
# Note:
# The second interval may include 0!
#
proc ::math::interval::/ {operand1 operand2} {
foreach {begin1 end1} [concat $operand1 $operand1] {break}
foreach {begin2 end2} [concat $operand2 $operand2] {break}
if { $begin2 <= 0.0 && $end2 >= 0.0 } { error "Division by interval containing zero" }
* $operand1 [list [expr {1.0/$begin2}] [expr {1.0/$end2}]]
}
# contains --
# Check if the interval contains a given value
#
# Arguments:
# interval Interval number in question
# value Value that is or is not contained
# Result:
# 0 if not, 1 if it is contained
#
proc ::math::interval::contains {interval value} {
foreach {begin end} [concat $interval $interval] {break}
if { $begin <= $value && $end >= $value } {
return 1
} else {
return 0
}
}
# interval --
# Create an interval number from a centre and a tolerance
#
# Arguments:
# centre Central number for the interval
# tolerance Tolerance (relative error, defaults to 1.0e-10)
# Result:
# Interval number
# Note:
# For centre = 0.0, the tolerance is chosen to be absolute
# - a hack in fact.
#
proc ::math::interval::interval { centre {tolerance 1.0e-10} } {
if { $centre != 0.0 } {
set begin [expr {(1.0-$tolerance)*$centre}]
set end [expr {(1.0+$tolerance)*$centre}]
} else {
set begin [expr {-$tolerance}]
set end $tolerance
}
return [list $begin $end]
}
# begin-end --
# Create an interval number from two extremes
#
# Arguments:
# low Lower bound for the interval
# high Higher bound for the interval
# Result:
# Interval number
#
proc ::math::interval::begin-end { low high } {
if { $low < $high } {
return [list $low $high]
} else {
return [list $high $low]
}
}
# square --
# Return the interval containing the square of an interval
#
# Arguments:
# interval Interval to be squared
# Result:
# Encompassing interval
# Note:
# The trick is that the square of a number can not be negative!
# And the operands are correlated so that a direct multiplication
# gives a very pessimistic answer
#
proc ::math::interval::square { interval } {
foreach {begin end} [concat $interval $interval] {break}
set min2 [expr {$begin*$begin}]
set max2 [expr {$end*$end}]
if { $begin < 0.0 && $end > 0.0 } {
return [list 0.0 [expr {$min2>$max2? $min2 : $max2}]]
} else {
return [list [expr {$min2<$max2? $min2 : $max2}] \
[expr {$min2>$max2? $min2 : $max2}] ]
}
}
# abs --
# Return the interval containing the absolute value of an interval
#
# Arguments:
# interval Interval to be treated
# Result:
# Encompassing interval
#
proc ::math::interval::abs { interval } {
foreach {begin end} [concat $interval $interval] {break}
set min [expr {abs($begin)}]
set max [expr {abs($end)}]
if { $begin < 0.0 && $end > 0.0 } {
return [list 0.0 [expr {$min>$max? $min : $max}]]
} else {
return [list [expr {$min<$max? $min : $max}] \
[expr {$min>$max? $min : $max}] ]
}
}
# Test code
#
namespace import ::math::interval::*
puts [+ {0.9 1.1} {2.1 2.2}]
puts [+ {0.9 1.1} 1]
puts [+ 1 1]
puts [- {0.9 1.1} {2.1 2.2}]
puts [* {0.9 1.1} {2.1 2.2}]
puts [/ {0.9 1.1} {2.1 2.2}]
puts "Square:"
puts [square {1 2}]
puts [square {-1 2}]
puts [square {-1 -2}]
puts "Abs:"
puts [abs {1 2}]
puts [abs {-1 2}]
puts [abs {-1 -2}]
#
# Now a more practical example
#
puts "Integrating a parabolic function f(x) = alpha * x**2"
puts "Assumes that all steps have independent uncertainties!"
set alpha [interval 1 0.001]
set alphapl 1
set dx 0.1
set sum 0.0
set sumpl 0.0
for { set i 0 } { $i < 300 } { incr i } {
set x [* $i $dx]
set sum [+ $sum [* $alpha [* $x $x]]]
set sumpl [+ $sumpl [* $alphapl [* $x $x]]]
}
puts "(Interval) estimate of integral: $sum"
puts "(Plain) estimate of integral: $sumpl"The output:
3.0 3.3 1.9 2.1 2 2 -1.2 -1.1 1.89 2.42 0.409090909091 0.52380952381 Square: 1 4 0.0 4 1 4 Abs: 1 2 0.0 2 1 2 Integrating a parabolic function f(x) = alpha * x**2 Assumes that all steps have independent uncertainties! (Interval) estimate of integral: 89460.9495 89640.0505 (Plain) estimate of integral: 89550.5 89550.5
TV I remember the idea from 1st years physics practicum: propagation of measurement errors. Maybe give the whole idea a standard deviation functional integral interpretation and we're in interesting physicist area..
A main point to think about which is very relevant in digital signal processing as opposed to electronics and other continuous equations based computations and approximation methods is the correlation measure of the error, also expressable in bits (the information unit), which comes from rounding errors in the computations themselves, also in the inevitable limited accuracy computations of the intervals, so you make sure they are real outer bounds.
For instance
(bin) 4 % set a [expr 1.0/6] 0.166666666667 (bin) 5 % expr 1.0/$a 5.99999999999
is simplistic example of what in senstive boundary computations (as in electrical engineering is well known rule in differential equations related problems: most problems arise in the boundary conditions) can easily make your 'certain' interval just wrong.
I just reread and found you mention correlation yourself. The accumulation of error does however not take all information away, the uncertainty interval becomes large, as does the product, their ratio doesn't change to dramatically I'd think. The worst case does, though, 10 . 10 =100, (10-1) . (10-1) = 81, so 10% * 10% --> ca. 20 %. It's of course just the expectation value which still carries information.
AM Yes, that is the main drawback of this type of calculations. You can do something about it, but I still have to digest the details (from the documentation of one of the compilers that support interval arithmetic).
t. In limited bit width audio and image dsp, some people noise up the data, to make the overall quantisation more based on randomness, so the effect of correlated repeated quantisation noise is made less correlates at the expense of overall decreased signal to noise ration. That can help making the overall energy spectrum errors smoothed out, so less pronounced peaks occur at certain points. In case it is possible to propagate the errors of each part of the computation, a bit accurate analysis is of course possible, and capable of exactly predicting the overall error bounds. With floating point computations, that can be quite hard, though overestimation should be possible.
If you'd have something like your own example, the multiplies should be well behaved enough to do actual error bound prediction, but when you'd have something like :
# a=10; b=10 expr 1/($a*$b-90)
things can get pretty tricky.