アフィン演算

Transcript

1. Interval Arithmetic Affine Arithmetic Calculation method for estimating the error
4. Error in function Consider the function of intersection. A "float" is not a real number, but a decimal number expressed to a finite number of digits. Since the calculation is done with a finite number of digits, there is an error in the result. float intersect( vec3 org, vec3 dir);
5. Error in function How can we describe the error? float_with_error intersect( vec3 org, vec3 dir);
6. Interval Arithmetic The naive way
7. Interval Arithmetic Expressing a number as an interval. That is, try to hold from the lowest to the highest value. // Interval Arithmetic(IA) class IA { private: float low_; float high_; };
8. IEEE 754 The result of a float or double calculation contains an error. But it does not produce values that are completely meaningless. If IEEE 754 is satisfied, +,−,×,/,√ and FMA will return the closest value to the result calculated with real numbers(infinite precision). In C++, To check if it is supported, use is_iec559. float a,b,c; // d is the closest float to a+b float d = a + b; // And so on float e = a * b; float f = a / b; float g = sqrt(a); float h = std::fma(a,b,c); // No guarantee float i = std::sin(a); // Check for IEEE 754 support static_assert( std::numeric_limits<float>:: is_iec559);
9. Neighboring float Floats have a rounding error in their calculations. Of course, the result of interval calculation also has it. In C++, you can use std::fesetround() to set the rounding direction. But it is expensive to execute. A simple way to make it safe is to always round outside the interval. std::nextafter() can be used to get the float next to it. IA(float low, float high) { low_ = std::nextafter(low, -inf); high_ = std::nextafter(high, +inf); }
10. Interval Arithmetic “+” Addition is just adding low to low, high to high. Interval should be the most conservative. // addition IA operator+(const IA& rhs) const { return IA(low_ + rhs.low_, high_ + rhs.high_); }
11. Interval Arithmetic “−” Unlike addition, subtraction does not result in a combination of low to low. It is easy to understand if you think of it as a worst case scenario where the interval widens. // subtraction IA operator-(const IA& rhs) const { return IA(low_ - rhs.high_, high_ - rhs.low_); }
12. Interval Arithmetic “×” Multiplication is a slightly more complicated. It considers all the cases where the interval has a negative value or not, and whether it crosses 0 or not. The largest and smallest of all the combinations of low and high will be the new interval low and high. // multiplication IA operator * (const IA& rhs) const { const float x0 = high_ * rhs.high_; const float x1 = high_ * rhs.low_; const float x2 = low_ * rhs.high_; const float x3 = low_ * rhs.low_; return IA(std::min({x0,x1,x2,x3}), std::max({x0,x1,x2,x3})); }
13. Interval Arithmetic “/” Division is the same as multiplication, with the maximum and minimum of all combinations being the interval. Only when the denominator is likely to be zero, exception handling is necessary. // division IA operator / (const IA& rhs) const { if ((rhs.low_ <= 0.0f)&&(0.0f <= rhs.high_)) { return IA(-inf, +inf); } else { const float x0 = high_ / rhs.high_; const float x1 = high_ / rhs.low_; const float x2 = low_ / rhs.high_; const float x3 = low_ / rhs.low_; return IA(std::min({ x0, x1, x2, x3 }), std::max({ x0, x1, x2, x3 })); } }
14. Interval Arithmetic “√” Since the square root is a monotonically increasing function, simply store the low and high results // square root IA sqrt() const { return IA( sqrt(max(low_,0.0f)), sqrt(max(high_,0.0f))); }
15. Interval Arithmetic The graph plots the function 𝑔(𝑥) using IA. 𝑔 𝑥 = 𝑥2 − 𝑥 + 1 2 𝑥2 + 1 2 The enclosing box is a vertical representation of the IA result when its width is the error. The results are properly conservative and also seem to be fine.
16. Interval Arithmetic Plot a function like the following. ℎ 𝑥 = 𝑔 𝑔 𝑥 The box is stretched too vertically, and the result looks too conservative.
17. Interval Arithmetic IA has the problem that the error is estimated too conservatively. For example, if you calculate “X – X” for a interval [4,6], The answer is [-2,2]. Thinking semantically, it has to be [0,0]. As a more subtle example, “X(10-X)” would be [16,36], but semantically it's actually [24,25]!! // [4,6]-[4,6] = [-2,2] IA ia(4.0, 6.0); IA ans0 = ia - ia; // [4,6]*(10-[4,6]) // = [4,6]*[4,6] // = [16,36] IA ans1 = ia * (10.0 - ia);
18. Interval Arithmetic To solve the “X-X” problem, we may need to track the symbols. To solve the “X(10-X)” problem, we may need to track where the error came from as well.
19. Affine arithmetic Keep tracking where the error came from.
20. Affine form In Affine form, variables are represented as follows. ො 𝑥 = 𝑥0 + 𝑥1 𝜀1 + 𝑥2 𝜀2 + ⋯ + 𝑥𝑛 𝜀𝑛 = 𝑥0 + ෍ 𝑖=1 𝑛 𝑥𝑖 𝜀𝑖 The symbol 𝜺𝒊 is an error that is assumed to be shared with other variables and takes the value between [-1,1]. 𝒙𝒊 is a scalar value that represents how much influence 𝜺𝒊 gives to ෝ 𝒙. 𝑥0 means the middle of the range of values, and is called central value.
21. Affine form To convert from interval to affine form, calculate the center and width of the interval range. ҧ 𝑥 = 𝑙𝑜𝑤, ℎ𝑖𝑔ℎ ො 𝑥 = ℎ𝑖𝑔ℎ + 𝑙𝑜𝑤 2 + ℎ𝑖𝑔ℎ − 𝑙𝑜𝑤 2 𝜀1 Give one example. ҧ 𝑥 = −2, 6 ො 𝑥 = 6 + −2 2 + 6 − −2 2 𝜀1 = 2 + 4𝜀1
22. Affine form To convert from affine form to interval, sum up the absolute values of all 𝑥𝑖 to get 𝜉. 𝜉 = ෍ 𝑖=1 𝑛 𝑥𝑖 Since this is the half extension of interval, the range from central value to 𝜉 is interval. ො 𝑥 = 𝑥0 − 𝜉, 𝑥0 + 𝜉 Give one example. ො 𝑥 = 2 + 1𝜀1 − 3𝜀2 ො 𝑥 = 2 − (|1| + | − 3|), 2 + (|1| + | − 3|)| = −2, 6
23. Affine form There are two affine quantities. ො 𝑥 = 20 − 4𝜀1 + 2𝜀3 + 3𝜀4 ො 𝑦 = 10 − 2𝜀1 + 1𝜀2 − 1𝜀4 These two intervals are as follows. ҧ 𝑥 = 11,29 ത 𝑦 = 6,14 But the region with these two affine quantities is tighter than 11,31 × [6,14], as shown in the figure on the right.
24. Affine arithmetic “+” “−” There are two affine quantities. ො 𝑥 = 𝑥0 + ෍ 𝑖=1 𝑛 𝑥𝑖 𝜀𝑖 ො 𝑦 = 𝑦0 + ෍ 𝑖=1 𝑛 𝑦𝑖 𝜀𝑖 Additions and subtractions simply add up each coefficient. ො 𝑥 ± ො 𝑦 = (𝑥0 ±𝑦0 ) + ෍ 𝑖=1 𝑛 (𝑥𝑖 ± 𝑦𝑖 ) 𝜀𝑖
25. Affine arithmetic “×” Compute the ො 𝑥 ∙ ො 𝑦 . ො 𝑥 ∙ ො 𝑦 = 𝑥0 + ෍ 𝑖=1 𝑛 𝑥𝑖 𝜀𝑖 ∙ 𝑦0 + ෍ 𝑖=1 𝑛 𝑦𝑖 𝜀𝑖 = 𝑥0 ∙ 𝑦0 + ෍ 𝑖=1 𝑛 𝑥0 𝑦𝑖 + 𝑥𝑖 𝑦0 𝜀𝑖 + ෍ 𝑖=1 𝑛 𝑥𝑖 𝜀𝑖 ∙ ෍ 𝑖=1 𝑛 𝑦𝑖 𝜀𝑖 = 𝑥0 ∙ 𝑦0 + ෍ 𝑖=1 𝑛 𝑥0 𝑦𝑖 + 𝑥𝑖 𝑦0 𝜀𝑖 + ෍ 𝑖=1 𝑛 ෍ 𝑗=1 𝑛 𝑥𝑖 𝑦𝑗 𝜀𝑖 𝜀𝑗 The first half part is a linear formula for 𝜺𝒊 .(It matches AA) But how do we do the second half?
26. Affine arithmetic “×” Place the second half of the part as 𝑄. 𝑄 = ෍ 𝑖=1 𝑛 ෍ 𝑗=1 𝑛 𝑥𝑖 𝑦𝑗 𝜀𝑖 𝜀𝑗 Once the maximum and minimum values of this function are obtained, it can be expressed with new noise symbols 𝜀𝑘 as follows. 𝑄 = 𝑚𝑎𝑥 + min 2 + 𝑚𝑎𝑥 − min 2 𝜀𝑘 Intuitively, we could calculate all the combinations of 𝜀𝑖 𝜀𝑗 and get the minimum and maximum, but this computational complexity would be O(𝑛2). Using a smart algorithm, we can drop it to O 𝑚 log 𝑚 where m is the number of non-zero parameters. But still, it's expensive.
27. Affine arithmetic “×” 𝑄 = ෍ 𝑖=1 𝑛 ෍ 𝑗=1 𝑛 𝑥𝑖 𝑦𝑗 𝜀𝑖 𝜀𝑗 The conservative estimator is as follows ത 𝑄 = −𝑢𝑣, 𝑢𝑣 𝑢 = ෍ 𝑖=1 𝑛 𝑥𝑖 𝑣 = ෍ 𝑖=1 𝑛 𝑦𝑖 The computational complexity is 𝑂(𝑛).
28. Affine arithmetic “×” Show an example. 𝑄 = 2𝜀1 + 1𝜀2 + 0𝜀3 −2𝜀1 + 0𝜀2 + 1𝜀3 𝑢 = |2| + |1| + |0| = 3 𝑣 = −2 + |0| + |1| = 3 𝑢𝑣 = 3 × 3 = 9 ത 𝑞 = −9,9 The actual interval is 𝑄𝑚𝑎𝑥 = 1  𝑤ℎ𝑒𝑛 𝜀1 , 𝜀2 , 𝜀3 = 0,1,1 𝑄𝑚𝑖𝑛 = −9 𝑤ℎ𝑒𝑛 𝜀1 , 𝜀2 , 𝜀3 = 1,1, −1 ത 𝑞 = −9,1
29. Affine arithmetic “×” Summarize the calculations so far. ො 𝑥 ∙ ො 𝑦 = 𝑥0 ∙ 𝑦0 + ෍ 𝑖=1 𝑛 𝑥0 𝑦𝑖 + 𝑥𝑖 𝑦0 𝜀𝑖 + ෍ 𝑖=1 𝑛 ෍ 𝑗=1 𝑛 𝑥𝑖 𝑦𝑗 𝜀𝑖 𝜀𝑗 ≈ 𝑥0 ∙ 𝑦0 + ෍ 𝑖=1 𝑛 𝑥0 𝑦𝑖 + 𝑥𝑖 𝑦0 𝜀𝑖 + ෍ 𝑖=1 𝑛 𝑥𝑖 ෍ 𝑖=1 𝑛 𝑦𝑖 𝜀𝑘
30. Affine arithmetic round off Unlike interval, the number doesn’t always get larger or smaller. At the time of calculating each coefficient 𝑧𝑖, neither roundup nor round down should be done. The differences should be summed up and added up to 𝑧𝑘.
31. Affine arithmetic “/” In AA, instead of directly implementing division, implement 1 𝑦 first, by using 𝑥 𝑦 = 𝑥 1 𝑦 .
32. Affine approximation
33. Affine approximation Suppose there is an interval [𝒂, 𝒃], and consider making 𝒚 = 𝒑𝒙 + 𝒒 ± 𝜹, which encloses y = 𝑓(𝑥) in that range. This [𝒂, 𝒃] is calculated from the 𝑥0 and 𝜉 of the affine form. 𝑎 = 𝑥0 − ξ 𝑏 = 𝑥0 + 𝜉 𝜉 = ෍ 𝑖=1 𝑛 𝑥𝑖
34. Affine approximation: MinRange Consider using the slope q as the slope of 𝑓(𝑥) at b(the end of the interval). That is, f'(b). 𝑝 ≔ 𝑓′ 𝑏 Since the function is convex, the line ±δ ( dashed line) always passes through 𝑓(𝑎) and 𝑓(𝑏), and the yellow line is drawn in the center of the lines. 𝑞 ≔ 𝑓 𝑎 + 𝑓 𝑏 − 𝑝 𝑎 + 𝑏 2 𝛿 ≔ 𝑓 𝑏 − 𝑓 𝑎 − 𝑝 𝑏 − 𝑎 2 This is called the Min-Range approximation.
35. Affine approximation The width goes from 1 to − 1 4 (∵ 1.25,1.0 ) The center goes from 1.5 to 1.125. Transform the original Affine form so that it fits within this range. That is... 𝑥𝑖 : Multiply by slope p(=-1/4). 𝑥𝑘 : Assign half the width of the dashed lines. 𝑥0 : Multiply by p and add the y-intercept (=9/8) of the line.
36. Affine approximation: Chebyshev There is a way to make it tighter. Apply the slope as the change in the function. 𝑝 ≔ 𝑓 𝑏 − 𝑓 𝑎 𝑏 − 𝑎 𝑞 ≔ 𝑓 𝑎 + 𝑓 𝜉 − 𝑝 𝑎 + 𝜉 2 𝛿 ≔ 𝑓 𝜉 − 𝑓 𝑎 − 𝑝 𝜉 − 𝑎 2 𝜉 is such that 𝑓′ 𝜉 = 𝑝. This is called the Chebyshev approximation.
37. Affine arithmetic “√” p, q, and δ can also be found for square roots
39. Result
40. Affine arithmetic Plot g(x) using AA. The result is conservative as in IA but even tighter.
41. Affine arithmetic For h(x), the result is significantly tighter than IA!!
42. Affine arithmetic
43. Wrap-up
44. Wrap-up Talked about Interval arithmetic(IA) and Affine arithmetic(AA) as representations of real number with error. IA keeps only the interval of error. The mechanism is very simple, but the errors increase exponentially. AA keeps track of where the error came from. It is more computationally intensive, but the error is smaller than IA. There are two affine approximations, Min-Range and Chebyshev. A very simple implementation has been submitted to GitHub.
45. References [1] Rump, S. M., & Kashiwagi, M. (2015). Implementation and improvements of affine arithmetic. Nonlinear Theory and Its Applications, IEICE, 6(3), 341–359. https://doi.org/10.1587/nolta.6.341 [2] Stolfi, J., & de Figueiredo, L. H. (2003). An Introduction to Affine Arithmetic. Trends in Computational and Applied Mathematics, 4(3), 297–312. https://doi.org/10.5540/tema.2003.04.03.0297 [3] Comba, J. L. D., & Stolfi, J. (1993). Affine Arithmetic and its Applications to Computer Graphics.