Slide 1

Slide 1 text

Interval Arithmetic Affine Arithmetic Calculation method for estimating the error

Slide 2

Slide 2 text

No content

Slide 3

Slide 3 text

No content

Slide 4

Slide 4 text

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

Slide 5

Slide 5 text

Error in function How can we describe the error? float_with_error intersect( vec3 org, vec3 dir);

Slide 6

Slide 6 text

Interval Arithmetic The naive way

Slide 7

Slide 7 text

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

Slide 8

Slide 8 text

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:: is_iec559);

Slide 9

Slide 9 text

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

Slide 10

Slide 10 text

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

Slide 11

Slide 11 text

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

Slide 12

Slide 12 text

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

Slide 13

Slide 13 text

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

Slide 14

Slide 14 text

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

Slide 15

Slide 15 text

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.

Slide 16

Slide 16 text

Interval Arithmetic Plot a function like the following. ℎ 𝑥 = 𝑔 𝑔 𝑥 The box is stretched too vertically, and the result looks too conservative.

Slide 17

Slide 17 text

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

Slide 18

Slide 18 text

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.

Slide 19

Slide 19 text

Affine arithmetic Keep tracking where the error came from.

Slide 20

Slide 20 text

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.

Slide 21

Slide 21 text

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

Slide 22

Slide 22 text

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

Slide 23

Slide 23 text

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.

Slide 24

Slide 24 text

Affine arithmetic “+” “−” There are two affine quantities. ො 𝑥 = 𝑥0 + ෍ 𝑖=1 𝑛 𝑥𝑖 𝜀𝑖 ො 𝑦 = 𝑦0 + ෍ 𝑖=1 𝑛 𝑦𝑖 𝜀𝑖 Additions and subtractions simply add up each coefficient. ො 𝑥 ± ො 𝑦 = (𝑥0 ±𝑦0 ) + ෍ 𝑖=1 𝑛 (𝑥𝑖 ± 𝑦𝑖 ) 𝜀𝑖

Slide 25

Slide 25 text

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?

Slide 26

Slide 26 text

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.

Slide 27

Slide 27 text

Affine arithmetic “×” 𝑄 = ෍ 𝑖=1 𝑛 ෍ 𝑗=1 𝑛 𝑥𝑖 𝑦𝑗 𝜀𝑖 𝜀𝑗 The conservative estimator is as follows ത 𝑄 = −𝑢𝑣, 𝑢𝑣 𝑢 = ෍ 𝑖=1 𝑛 𝑥𝑖 𝑣 = ෍ 𝑖=1 𝑛 𝑦𝑖 The computational complexity is 𝑂(𝑛).

Slide 28

Slide 28 text

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

Slide 29

Slide 29 text

Affine arithmetic “×” Summarize the calculations so far. ො 𝑥 ∙ ො 𝑦 = 𝑥0 ∙ 𝑦0 + ෍ 𝑖=1 𝑛 𝑥0 𝑦𝑖 + 𝑥𝑖 𝑦0 𝜀𝑖 + ෍ 𝑖=1 𝑛 ෍ 𝑗=1 𝑛 𝑥𝑖 𝑦𝑗 𝜀𝑖 𝜀𝑗 ≈ 𝑥0 ∙ 𝑦0 + ෍ 𝑖=1 𝑛 𝑥0 𝑦𝑖 + 𝑥𝑖 𝑦0 𝜀𝑖 + ෍ 𝑖=1 𝑛 𝑥𝑖 ෍ 𝑖=1 𝑛 𝑦𝑖 𝜀𝑘

Slide 30

Slide 30 text

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 𝑧𝑘.

Slide 31

Slide 31 text

Affine arithmetic “/” In AA, instead of directly implementing division, implement 1 𝑦 first, by using 𝑥 𝑦 = 𝑥 1 𝑦 .

Slide 32

Slide 32 text

Affine approximation

Slide 33

Slide 33 text

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 𝑛 𝑥𝑖

Slide 34

Slide 34 text

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.

Slide 35

Slide 35 text

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.

Slide 36

Slide 36 text

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.

Slide 37

Slide 37 text

Affine arithmetic “√” p, q, and δ can also be found for square roots

Slide 38

Slide 38 text

No content

Slide 39

Slide 39 text

Result

Slide 40

Slide 40 text

Affine arithmetic Plot g(x) using AA. The result is conservative as in IA but even tighter.

Slide 41

Slide 41 text

Affine arithmetic For h(x), the result is significantly tighter than IA!!

Slide 42

Slide 42 text

Affine arithmetic

Slide 43

Slide 43 text

Wrap-up

Slide 44

Slide 44 text

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.

Slide 45

Slide 45 text

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.