Computational Divided Differencing

# Computational Divided Differencing

## Computational Divided Differencing

- - - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - - -
##### Presentation Transcript

1. ComputationalDivided Differencing ThomasReps UniversityofWisconsin 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. Backward slicewith respect to “printf(“%d\n”,i)” Example 1: Program Slicing int Sum() { int sum = 0; int i = 1; while (i < 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 < 11) { sum = sum + i; i = i + 1; } printf(“%d\n”,sum); printf(“%d\n”,i); } Backward slicewith respect to “printf(“%d\n”,i)”

5. Example 1: Program Slicing int Sum_Select_i() { int i = 1; while (i < 11) { i = i + 1; } printf(“%d\n”,i); } Backward slicewith respect to “printf(“%d\n”,i)”

6. Backward slicewith respect to “printf(“%d\n”,sum)” Example 1: Program Slicing int Sum() { int sum = 0; int i = 1; while (i < 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 < 11) { sum = sum + i; i = i + 1; } printf(“%d\n”,sum); printf(“%d\n”,i); } Backward slicewith respect to “printf(“%d\n”,sum)”

8. Example 1: Program Slicing int Sum_Select_sum() { int sum = 0; int i = 1; while (i < 11) { sum = sum + i; i = i + 1; } printf(“%d\n”,sum); } Backward slicewith respect to “printf(“%d\n”,sum)”

9. Example 1: Program Slicing • Given • F: program • : projection function • CreateF, such that • F(x) = (F) (x) • SumSelect[i](x) = (Select[i] Sum)(x) • SumSelect[sum](x) = (Select[sum] Sum)(x)

10. 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)

11. Trace( ,r); – Example 2: Partial Evaluation for(each pixel p){ Ray r = Ray(eye,p); Trace(object,r); }

12. Trace( ,r); – Example 2: Partial Evaluation for(each pixel p){ Ray r = Ray(eye,p); Trace(object,r); }

13. Example 2: Partial Evaluation TraceObj = PEval(Trace,object); for(each pixel p){ Ray r = Ray(eye,p); TraceObj(r); }

14. Example 2: Partial Evaluation TraceObj = PEval(Trace,object); for(each pixel p){ Ray r = Ray(eye,p); TraceObj(r); }

15. F(x0) Comp. Diff. d dx  F’(x0) F’(x0) Ex. 3: Computational Differentiation[“Automatic Differentiation”] TransformFF’ F(x0)

16. Computational Differentiation

17. Computational Differentiation float F(float x) { int i; float ans = 1.0; for(i = 1; i <= n; i++) { ans = ans * f[i](x); } return ans; } float delta = . . .; /* small constant */ float F’(float x) { return (F(x+delta) - F(x)) / delta; }

18. Computational Differentiation float F (float x) { int i; float ans = 1.0; for(i = 1; i <= n; i++) { ans = ans * f[i](x); } return ans ; }

19. Computational Differentiation floatF’(float x) { int i; float ans’ = 0.0; float ans = 1.0; for(i = 1; i <= n; i++) { ans’ = ans * f’[i](x) + ans’ * f[i](x); ans = ans * f[i](x); } returnans’; }

20. Iter. ans’ ans Computational Differentiation ans’ = ans * f’[i](x) + ans’ * f[i](x); ans = ans * f[i](x); 0 0 1 1 f1’(x) f1(x) 2 f1*f2’ + f1’*f2 f1*f2 3 f1*f2*f3’ + f1*f2’*f3 + f1’*f2*f3 f1*f2*f3 … . . . . . .

21. 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; }

22. 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; }

23. Differentiation Arithmetic float F(float x) { int i; float ans = 1.0; for(i = 1; i <= n; i++) { ans = ans * f[i](x); } return ans; }

24. Differentiation Arithmetic FloatD F(FloatD x) { int i; FloatD ans = 1.0; for(i = 1; i <= n; i++) { ans = ans * f[i](x); } return ans; }

25. Differentiation Arithmetic FloatD F(FloatD x) { int i; FloatD ans = 1.0; // ans.val’ = 0.0 // ans.val = 1.0 for(i = 1; i <= n; i++) { ans = ans * f[i](x); } return ans; }

26. Differentiation Arithmetic FloatD F(FloatD x) { int i; FloatD ans = 1.0; for(i = 1; i <= n; i++) { ans = ans * f[i](x); } return ans; } Similarly transformed version off[i]

27. Differentiation Arithmetic FloatD F(FloatD x) { int i; FloatD ans = 1.0; for(i = 1; i <= n; i++) { ans = ans * f[i](x); } return ans; } Overloaded operator*

28. Differentiation Arithmetic FloatD F(FloatD x) { int i; FloatD ans = 1.0; for(i = 1; i <= n; i++) { ans = ans * f[i](x); } return ans; } Overloaded operator=

29. 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)

30. F v F(v) Computational Differentiation v F(v) Differentiating version of F w F’(v)*w Differentiation Arithmetic

31. G F v G(v) F(G(v)) Computational Differentiation Computational Differentiation F(G(v)) Differentiating version of F v G(v) Differentiating version of G 1.0 G’(v) F’(G(v))*G’(v) Chain Rule

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

33. Accurate Computation of Divided Differences

34. F F(x1) F(x0) x0 x1 Divided Differences = Slopes

35. Interpolation/extrapolation Numerical integration Solving differential equations Who Cares About Divided Differences? Scientific and graphics programs: Predict and render modeled situations

36. round-off error! division by a small quantity No Way! 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); }

37. F [x0,x1] = F(x0) –F(x1) x0–x1 = x02–x12 x0–x1 But … Consider F(x) = x2, when x0 x1 = x0 + x1

38. 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; }

39. Example: Evaluation via Horner's rule class Poly { public: Poly(int N, float x[]); float Eval(float); private: int degree; float *coeff; // Array coeff[0..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->Eval(3.005); P(x) = ((((0.0 * x + 2.1) * x– 1.4) * x– .6) * x + 1.1)

40. Example: Evaluation via Horner's rule P(x) = ((((0.0 * x + 2.1) * x– 1.4) * x– .6) * x + 1.1) float Poly::Eval(float x) { int i; float ans = 0.0; for (i = degree; i >= 0; i--) { ans = ans * x + coeff[i]; } return ans; }

41. Example: Evaluation via Horner's rule P(x) = ((((0.0 * x + 2.1) * x– 1.4) * x– .6) * x + 1.1) float Poly::Eval_1DD(float x0,float x1) { int i; float ans_1DD = 0.0; float ans = 0.0; for (i = degree; i >= 0; i--) { ans_1DD = ans_1DD * x1 + ans; ans = ans * x0 + coeff[i]; } return ans_1DD; }

42. Laws for F[x0,x1] 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) c[x0,x1] = 0 x[x0,x1] = 1 (c+F)[x0,x1] = F [x0,x1] (c*F)[x0,x1] = c*F[x0,x1] (F+G)[x0,x1] = F[x0,x1]+G[x0,x1] (F*G) [x0,x1] = F [x0,x1]*G(x1) + F(x0)*G[x0,x1]

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

44. F(x0) F(x0) - F(x1) x0 - x1 Standard Div. Diff. Comp. Diff. d dx F [x0,x1]  F’(x0) F’(x0) x1x0 F(x0)-F(x1) x0-x1 x1x0 ? F(x0)

45. F(x0)  F(x0) F(x0) - F(x1) x0 - x1 Comp. Div. Diff. Comp. Diff. d dx F [x0,x1]  F’(x0) F’(x0) F[x0,x1] x1x0 x1x0

46. d dz F(z) F[x0, x0] = z = x0 First Divided Difference if x0x1

47. Laws for F[x0,x1] 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) c[x0,x1] = 0 x[x0,x1] = 1 (c+F)[x0,x1] = F [x0,x1] (c*F)[x0,x1] = c*F[x0,x1] (F+G)[x0,x1] = F[x0,x1]+G[x0,x1] (F*G) [x0,x1] = F [x0,x1]*G(x1) + F(x0)*G[x0,x1]

48. Laws for F[x0,x0] 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) c[x0,x0] = 0 x[x0,x0] = 1 (c+F)[x0,x0] = F [x0,x0] (c*F)[x0,x0] = c*F[x0,x0] (F+G)[x0,x0] = F[x0,x0]+G[x0,x0] (F*G) [x0,x0] = F [x0,x0]*G(x0) + F(x0)*G[x0,x0]

49. Example: Evaluation via Horner's rule float Poly::Eval(float x) { int i; float ans = 0.0; for (i = degree; i >= 0; i--) { ans = ans * x + coeff[i]; } return ans; }

50. Example: Evaluation via Horner's rule float Poly::Eval_1DD(float x0,float x1) { int i; float ans_1DD = 0.0; float ans = 0.0; for (i = degree; i >= 0; i--) { ans_1DD = ans_1DD * x1 + ans; ans = ans * x0 + coeff[i]; } return ans_1DD; }