1 / 105

Computational Divided Differencing

Computational Divided Differencing. Thomas Reps University of Wisconsin. Joint work with Louis B. Rall. 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 ) ).

gustav
Download Presentation

Computational Divided Differencing

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

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

More Related