Computational Divided Differencing

1
ComputationalDivided Differencing

• Thomas Reps
• University of Wisconsin

Joint work with Louis B. Rall
2
Program Transformation

• Transform F ? F
• For all x, F(x) and F (x) perform related
  computations
• Sometimes F needs x to be preprocessed
• F (h(x))

3
Example 1 Program Slicing
int Sum() int sum 0 int i 1 while (i lt
11) sum sum i i i
1 printf(d\n,sum) printf(d\n,i)
4
Example 1 Program Slicing
int Sum() int sum 0 int i 1 while (i lt
11) sum sum i i i
1 printf(d\n,sum) printf(d\n,i)
Backward slice with respect to printf(d\n,i)
5
Example 1 Program Slicing
int Sum_Select_i() int i 1 while (i lt
11) i i 1 printf(d\n,i)
Backward slice with respect to printf(d\n,i)
6
Example 1 Program Slicing
int Sum() int sum 0 int i 1 while (i lt
11) sum sum i i i
1 printf(d\n,sum) printf(d\n,i)
7
Example 1 Program Slicing
int Sum() int sum 0 int i 1 while (i lt
11) sum sum i i i
1 printf(d\n,sum) printf(d\n,i)
Backward slice with respect to printf(d\n,sum)

8
Example 1 Program Slicing
int Sum_Select_sum() int sum 0 int i
1 while (i lt 11) sum sum i i i
1 printf(d\n,sum)
Backward slice with respect to printf(d\n,sum)

9
Example 1 Program Slicing

• Given
• F program
• ? projection function
• Create F ?, such that
• F ?(x) (? ? F) (x)
• SumSelecti(x) (Selecti ? Sum)(x)
• SumSelectsum(x) (Selectsum ? Sum)(x)

10
Forward Slice
int Sum() int sum 0 int i 1 while (i lt
11) sum sum i i i
1 printf(d\n,sum) printf(d\n,i)
11
Forward Slice
int Sum() int sum 0 int i 1 while (i lt
11) sum sum i i i
1 printf(d\n,sum) printf(d\n,i)
Forward slice with respect to sum 0
12
Who Cares About Slices?

• Understanding programs
• Restructuring Programs
• Program Specialization and Reuse
• Program Differencing
• Testing (and Retesting)
• Year 2000 Problem
• Automatic Differentiation

13
What Are Slices Useful For?

• Understanding Programs
• What is affected by what?
• Restructuring Programs
• Isolation of separate computational threads
• Program Specialization and Reuse
• Slices specialized programs
• Only reuse needed slices
• Program Differencing
• Compare slices to identify changes
• Testing
• What new test cases would improve coverage?
• What regression tests must be rerun after a
  change?

14
Example 2 Partial Evaluation

• Given
• F D1 ? D2 ? D3
• s D1
• Create F s D2 ? D3, such that,
• F s(d) F(?s,d?)
• F s(d) F s(second(?s,d?)) F(?s,d?)

15
Example 2 Partial Evaluation
for(each pixel p) Ray r Ray(eye,p)
Trace(object,r)
16
Example 2 Partial Evaluation
for(each pixel p) Ray r Ray(eye,p)
Trace(object,r)
17
Example 2 Partial Evaluation
TraceObj PEval(Trace,object) for(each pixel
p) Ray r Ray(eye,p) TraceObj(r)
18
Example 2 Partial Evaluation
TraceObj PEval(Trace,object) for(each pixel
p) Ray r Ray(eye,p) TraceObj(r)
19
Line-Character-Count Program
void line_char_count(FILE f) int lines
0 int chars BOOL eof_flag FALSE int
n extern void scan_line(FILE f, BOOL bptr,
int iptr) scan_line(f, eof_flag, n) chars
n while(eof_flag FALSE) lines lines
1 scan_line(f, eof_flag, n) chars chars
n printf(lines d\n,
lines) printf(chars d\n, chars)
20
Character-Count Program
void char_count(FILE f) int lines 0 int
chars BOOL eof_flag FALSE int n extern
void scan_line(FILE f, BOOL bptr, int
iptr) scan_line(f, eof_flag, n) chars
n while(eof_flag FALSE) lines lines
1 scan_line(f, eof_flag, n) chars chars
n printf(lines d\n,
lines) printf(chars d\n, chars)
21
Line-Character-Count Program
void line_char_count(FILE f) int lines
0 int chars BOOL eof_flag FALSE int
n extern void scan_line(FILE f, BOOL bptr,
int iptr) scan_line(f, eof_flag, n) chars
n while(eof_flag FALSE) lines lines
1 scan_line(f, eof_flag, n) chars chars
n printf(lines d\n,
lines) printf(chars d\n, chars)
22
Line-Count Program
void line_count(FILE f) int lines 0 int
chars BOOL eof_flag FALSE int n extern
void scan_line2(FILE f, BOOL bptr, int
iptr) scan_line2(f, eof_flag, n) chars
n while(eof_flag FALSE) lines lines
1 scan_line2(f, eof_flag, n) chars
chars n printf(lines d\n,
lines) printf(chars d\n, chars)
23
Specialization Via Slicing
wc -lc
Not partial evaluation!
void line_char_count(FILE f)
24
Ex. 3 Computational DifferentiationAutomatic
Differentiation
Transform F ? F
F(x0)
25
Computational Differentiation
26
Computational Differentiation
float F(float x) int i float ans
1.0 for(i 1 i lt n i) ans ans
fi(x) return ans
float delta . . . / small constant / float
F(float x) return (F(xdelta) - F(x)) /
delta
27
Computational Differentiation
float F (float x) int i float ans
1.0 for(i 1 i lt n i) ans ans
fi(x) return ans
28
Computational Differentiation
float F(float x) int i float ans
0.0 float ans 1.0 for(i 1 i lt n i)
ans ans fi(x) ans fi(x) ans
ans fi(x) return ans
29
Computational Differentiation
ans ans fi(x) ans fi(x) ans ans
fi(x)
0 0 1 1 f1(x) f1(x) 2
f1f2 f1f2 f1f2 3 f1f2f3
f1f2f3 f1f2f3 f1f2f3 . . .
. . .
30
Differentiation Arithmetic
class FloatD public float val'
float val FloatD(float)
FloatD operator(FloatD a, FloatD b) FloatD
ans ans.val' a.val' b.val' ans.val
a.val b.val return ans
31
Differentiation Arithmetic
class FloatD public float val'
float val FloatD(float)
FloatD operator(FloatD a, FloatD b) FloatD
ans ans.val' a.val b.val' a.val'
b.val ans.val a.val b.val return ans
32
Differentiation Arithmetic
float F(float x) int i float ans
1.0 for(i 1 i lt n i) ans ans
fi(x) return ans
33
Differentiation Arithmetic
FloatD F(FloatD x) int i FloatD ans
1.0 for(i 1 i lt n i) ans ans
fi(x) return ans
34
Differentiation Arithmetic
FloatD F(FloatD x) int i FloatD ans 1.0
// ans.val 0.0
// ans.val 1.0 for(i 1 i lt n i)
ans ans fi(x) return ans
35
Differentiation Arithmetic
FloatD F(FloatD x) int i FloatD ans
1.0 for(i 1 i lt n i) ans ans
fi(x) return ans
Similarly transformed version of fi
36
Differentiation Arithmetic
FloatD F(FloatD x) int i FloatD ans
1.0 for(i 1 i lt n i) ans ans
fi(x) return ans
Overloaded operator
37
Differentiation Arithmetic
FloatD F(FloatD x) int i FloatD ans
1.0 for(i 1 i lt n i) ans ans
fi(x) return ans
Overloaded operator
38
Laws for F
c 0 x 1 (c F)(x) F(x) (c F)(x) c
F(x) (F G)(x) F(x) G(x) (F G)(x)
F(x) G(x) F(x) G(x)
39
Differentiation Arithmetic
40
Chain Rule
41
Computational Differentiation
x1
y1
xi
yj1
xm
y F(x)
yn
42
The If Problem
float F(float x) float ans if (x 1.0)
ans 1.0 else ans xx
return ans
43
The If Problem
float F(float x) float ans' float ans
if (x 1.0) ans' 0.0 ans 1.0
else ans' xx ans xx
return ans'
F(1.0) 0.0 versus F(1.0) 2.0
44
Outline

• Programs as data objects
• Partial evaluation
• Slicing
• Computational differentiation (CD)
• Computational divided differencing
• CDD generalizes CD
• Higher-order CDD
• Multi-variate CDD

45
Accurate Computation of Divided Differences
46
Divided Differences Slopes
47
Who Cares About Divided Differences?
Scientific and graphics programs Predict and
render modeled situations

• Interpolation/extrapolation
• Numerical integration
• Solving differential equations

48
Computing Divided Differences
float Square(float x) return x x
float Square_1DD_naive(float x0,float x1)
return (Square(x0) - Square(x1))/(x0 - x1)
49
But
Consider F(x) x2, when x0 ? x1
x0 x1
50
Computing Divided Differences
float Square(float x) return x x
float Square_1DD_naive(float x0,float x1)
return (Square(x0) - Square(x1))/(x0 - x1)
float Square_1DD(float x0,float x1) return x0
x1
51
Example Evaluation via Horner's rule
class Poly public Poly(int N, float
x) float Eval(float) private int
degree float coeff // Array
coeff0..degree
P(x) 2.1 x3 1.4 x2 .6 x 1.1
Poly P new Poly(4,2.1,-1.4,-0.6,1.1) float
val P-gtEval(3.005)
P(x) ((((0.0 x 2.1) x 1.4) x .6)
x 1.1)
52
Example Evaluation via Horner's rule
P(x) ((((0.0 x 2.1) x 1.4) x .6)
x 1.1)
float PolyEval(float x) int i float ans
0.0 for (i degree i gt 0 i--)
ans ans x coeffi return ans
53
Example Evaluation via Horner's rule
P(x) ((((0.0 x 2.1) x 1.4) x .6)
x 1.1)
float PolyEval_1DD(float x0,float x1) int
i float ans_1DD 0.0 float ans 0.0
for (i degree i gt 0 i--) ans_1DD
ans_1DD x1 ans ans ans x0
coeffi return ans_1DD
54
Laws for Fx0,x1
c 0 x 1 (cF)(x) F(x) (cF)(x)
cF(x) (FG)(x) F(x)G(x) (FG)(x)
F(x)G(x) F(x)G(x)
cx0,x1 0 xx0,x1 1 (cF)x0,x1 F
x0,x1 (cF)x0,x1 cFx0,x1 (FG)x0,x1
Fx0,x1Gx0,x1 (FG) x0,x1 F
x0,x1G(x1)
F(x0)Gx0,x1
55
Outline

• Programs as data objects
• Partial evaluation
• Slicing
• Computational differentiation (CD)
• Computational divided differencing
• CDD generalizes CD
• Higher-order CDD
• Multi-variate CDD

56
F(x0)
57
(No Transcript)
58
First Divided Difference
if x0 ? x1
59
Laws for Fx0,x1
c 0 x 1 (cF)(x) F(x) (cF)(x)
cF(x) (FG)(x) F(x)G(x) (FG)(x)
F(x)G(x) F(x)G(x)
cx0,x1 0 xx0,x1 1 (cF)x0,x1 F
x0,x1 (cF)x0,x1 cFx0,x1 (FG)x0,x1
Fx0,x1Gx0,x1 (FG) x0,x1 F
x0,x1G(x1)
F(x0)Gx0,x1
60
Laws for Fx0,x0
c 0 x 1 (cF)(x) F(x) (cF)(x)
cF(x) (FG)(x) F(x)G(x) (FG)(x)
F(x)G(x) F(x)G(x)
cx0,x0 0 xx0,x0 1 (cF)x0,x0 F
x0,x0 (cF)x0,x0 cFx0,x0 (FG)x0,x0
Fx0,x0Gx0,x0 (FG) x0,x0 F
x0,x0G(x0)
F(x0)Gx0,x0
61
Example Evaluation via Horner's rule
float PolyEval(float x) int i float ans
0.0 for (i degree i gt 0 i--)
ans ans x coeffi return ans
62
Example Evaluation via Horner's rule
float PolyEval_1DD(float x0,float x1) int
i float ans_1DD 0.0 float ans 0.0
for (i degree i gt 0 i--) ans_1DD
ans_1DD x1 ans ans ans x0
coeffi return ans_1DD
63
Example Evaluation via Horner's rule
float PolyEval_1DD(float x0,float x1) int
i float ans_1DD 0.0 float ans 0.0
for (i degree i gt 0 i--) ans_1DD
ans_1DD v ans ans ans v coeffi
return ans_1DD
64
Example Evaluation via Horner's rule
float PolyEval'(float x) int i float
ans' 0.0 float ans 0.0 for (i
degree i gt 0 i--) ans' ans' x
ans ans ans x coeffi return
ans'
65
Example Evaluation via Horner's rule
float PolyEval'(float x) int i float
ans' 0.0 float ans 0.0 for (i
degree i gt 0 i--) ans' ans' v
ans ans ans v coeffi return
ans'
Eval_1DD(v,v) Eval(v)
66
Outline

• Programs as data objects
• Partial evaluation
• Slicing
• Computational differentiation (CD)
• Computational divided differencing
• CDD generalizes CD
• Higher-order CDD
• Multi-variate CDD

67
Accurate Computation of Divided Differences
68
F x0, x1 , x2 Slopes of Slopes
69
Divided-Difference Table
F(x0) Fx0,x1 F(x1) Fx0,x1,x2 Fx1,x
2 Fx0,x1,x2,x3 F(x2)
Fx1,x2,x3 Fx2,x3 F(x3)
70
Divided-Difference Table
0 0 0 0 0 0
F?x0, x1, x2, x3
71
Laws for F?x0,,xn
c? cI x? Ax0,,xn (cF)?(x)
F?(x) (cF)?(x) cF?(x) (FG)?(x)
F?(x)G?(x) (FG)?(x) F?(x)G?(x)
cx0,x1 0 xx0,x1 1 (cF)x0,x1
Fx0,x1 (cF)x0,x1 cFx0,x1 (FG)x0,x1
Fx0,x1Gx0,x1 (FG) x0,x1 F
x0,x1G(x1)
F(x0)Gx0,x1
72
2.1? 2.1I
73
2.1? 2.1I
74
2.1x3 1.4x2 .6x 1.1
x0 3.0 x1 4.0 x2 5.0 x3 7.0
75
2.1x3 1.4x2 .6x 1.1
x0 3.0 x1 3.001 x2 3.002 x3 3.005
76
2.1x3 1.4x2 .6x 1.1
x0 3.0 x1 3.0001 x2 3.0002 x3
3.0005
77
2.1x3 1.4x2 .6x 1.1
x0 3.0 x1 3.0 x2 3.0 x3 3.0
78
Example Evaluation via Horner's rule
class Poly public Poly(int N, float
x) float Eval(float) private int
degree float coeff // Array
coeff0..degree
79
Example Evaluation via Horner's rule
class Poly public Poly(int N, float
x) float Eval(float) FloatDD
Eval(FloatDD) private int degree
float coeff // Array coeff0..degree
P(x) 2.1 x3 1.4 x2 .6 x 1.1
Poly P new Poly(4,2.1,-1.4,-0.6,1.1) FloatDD
xs(4, 3.0, 4.0, 5.0, 7.0) FloatDD ddt
P-gtEval(xs)
80
Example Evaluation via Horner's rule
float PolyEval(float x) int i float ans
0.0 for (i degree i gt 0 i--) ans
ans x coeffi return ans
81
Example Evaluation via Horner's rule
FloatDD PolyEval(FloatDD x) int i
FloatDD ans 0.0 for (i degree i gt 0
i--) ans ans x coeffi return
ans
82
2.1x3 1.4x2 .6x 1.1
83
2.1x3 1.4x2 .6x 1.1
84
2.1x3 1.4x2 .6x 1.1
85
2.1x3 1.4x2 .6x 1.1
86
2.1x3 1.4x2 .6x 1.1
87
Example Plotting a Function
88
Outline

• Programs as data objects
• Partial evaluation
• Slicing
• Computational differentiation (CD)
• Computational divided differencing
• CDD generalizes CD
• Higher-order CDD
• Multi-variate CDD

89
Divided-Difference Table A Polynomial
2.1x3 1.4x2 .6x 1.1
x0 3.0 x1 4.0 x2 5.0 x3 7.0
90
Newton Interpolating Polynomial
91
Divided-Difference Table A Polynomial
2.1x3 1.4x2 .6x 1.1
x0 3.0 x1 3.0 x2 3.0 x3 3.0
92
Taylor Series
(xx0)m Fx0,,x0
m1 x0s
93
Divided-Difference Table A Polynomial
94
Divided-Difference Table A Polynomial
95
Divided-Difference Table A Polynomial
96
Divided-Difference Table A Polynomial
97
Divided-Difference Table A Polynomial
98
Linear Interpolation
99
Linear Interpolation
100
Linear Interpolation
101
Linear Interpolation
102
First Divided Difference Slope
103
Divided-Difference Table A Polynomial
x1
x0
104
First Divided Difference Slope
105
First Divided Difference Slope
x0
x1
106
Divided-Difference Table A Surface
0
x0
x1
107
Avoid the Subtractions!
Use a divided-difference arithmetic Divided-diffe
rence tables of divided-difference tables
108
Quadratic Interpolation
109
Divided-Difference Table A Surface
110
Example Evaluation via Horner's rule
class BivariatePoly public float
Eval(float,float) private int
degree1,degree2 // Array coeff0..degree10.
.degree2 float coeff
float val P-gtEval(x,y)
111
Example Evaluation via Horner's rule
float BivariatePolyEval(float x,float y)
float ans 0.0 for (int i degree1 i gt 0
i--) float temp 0.0 for (int j
degree2 j gt 0 j--) temp temp y
coeffij ans ans x temp
return ans
112
Example Evaluation via Horner's rule
DDAlt2gt BivariatePolyEval(DDAlt2gt x, DDAlt2gty)

DDAlt2gt ans 0.0 for (int i degree1 i
gt 0 i--) DDAlt2gt temp 0.0 for (int
j degree2 j gt 0 j--) temp temp y
coeffij ans ans x temp
return ans
DDAlt2gt x(VAR, 1, x0, x1, x2, x3) DDAlt2gt y(VAR,
2, y0, y1, y2) DDAlt2gt surface P-gtEval(x,y)
113
Relationships Between CDD and CD
F?x0, x1, x2, x3
114
Relationships Between CDD and CD
F?x0, x0, x0, x0
115
Relationships Between CDD and CD
116
Relationships Between CDD and CD
117
Relationships Between CDD and CD
118
Related Work

• Computational differentiation
• Wengert, Rall, Iri, Griewank, Bischof, Carle,
• Accurate divided differencing
• Matrix operations Opitz 64 (McCurdy 80, Reps 98)
• Library functions McCurdy 80, Kahan Fateman
  85, Rall Reps 00
• Controlling round-off error
• Kulischs group at Karlsruhe
• Operator overloading
• partial evaluation, safe pointers, expression
  templates