4 ms·

I have worked on a performance-portable math library. We implement BLAS, sparse matrix operations, a variety of solvers, some ODE stuff, and various utilities f

I have worked on a performance-portable math library. We implement BLAS, sparse matrix operations, a variety of solvers, some ODE stuff, and various utilities for a variety of serial and parallel execution modes on x86, ARM, and the three major GPU vendors. The simplest and highest-impact tests are all the edge cases - if an input matrix/vector/scalar is 0/1/-1/NaN, that usually tells you a lot about what the outputs should be. It can be difficult to determine sensible numerical limit for error in the algorithms. The simplest example is a dot product - summing floats is not associative, so doing it in parallel is not bitwise the same as serial. For dot in particular it's relatively easy to come up with an error bound, but for anything more complicated it takes a particular expertise that is not always available. This has been a work in progress, and sometimes (usually) we just picked a magic tolerance out of thin air that seems to work. Solvers are tested using analytical solutions and by inverting them, e.g. if we're solving Ax = y, for x, then Ax should come out "close" to the original y (see error tolerance discussion above). One of the most surprising things to me is that the suite has identified many bugs in vendor math libraries (OpenBLAS, MKL, cuSparse, rocSparse, etc.) - a major component of what we do is wrap up these vendor libraries in a common interface so our users don't have to do any work when they switch supercomputers, so in practice we test them all pretty thoroughly as well. Maybe I can let OpenBLAS off the hook due to the wide variety of systems they support, but I expected the other vendors would do a better job since they're better-resourced. For this reason we find regression tests to be useful as well.

essexedwards 2y ago
I've also been surprised many times by issues in numerical libraries. In addition to matrices with simple entries, I've found plenty of bugs just testing small matrices, with dimensions in {0,1,2,3,4}. Many libraries/routines fall over when the matrix is small, especially when one dimension is 0 or 1. Presently, I am working on cuSPARSE and I'm very keen to improve its testing and correctness. I would appreciate anything more you can share about bugs you've seen in cuSPARSE. Feel free to email me, eedwards at nvidia.com
- AlotOfReading 2y ago
  This is one of the reasons I argue that it's almost always better to prioritize speed and stability than accuracy specifically. No one actually knows what their thresholds are (including library authors), but the sky isn't falling despite that. Instabilities and nondeterminism will blow up a test suite pretty quickly though.
  - phdelightful 2y ago
    Yeah, we really just try to come up with very loose bounds since the analysis is hard. Even so, it does occasionally stop us from getting things way way wrong.
  - essexedwards 2y ago
    > No one actually knows what their thresholds are (including library authors) If low-level numerical libraries provided documentation for their accuracy guarantees, it would make it easier to develop software on top of those libraries. I think numerical libraries should be doing this, when possible. It's already common for special-function (e.g. sin, cos, sqrt) libraries to specify their accuracy in ULPs. It's less common for linear algebra libraries to specify their accuracy, but it's still quite doable for BLAS-like operations.
    - AlotOfReading 2y ago
      What I'm trying to convey is that the required accuracies for the application are what's unclear. To give an example of a case where accuracy matters, I regularly catch computational geometry folks writing code that branches differently on positive, negative, and 0 results. That application implies 0.5 ulp, which obviously doesn't match the actual implementation accuracy even if it's properly specified, so there's usually a follow-up conversation trying to understand what they really need and helping them achieve it.
- RhysU 2y ago
  Nx0, 0xN, and 0x0 matrices are great edge cases. 0-length vectors, too. Then, do the same with Nx1, 1xN, and 1x1 matrices.
dahart 2y ago
> sometimes (usually) we just picked a magic tolerance out of thin air that seems to work. Probably worth mentioning that in general the tolerance should be relative error, not absolute, for floating point math. Absolute error tolerance should only be used when there’s a maximum limit on the magnitude of inputs, or the problem has been analyzed and understood. I know that doesn’t stop people from just throwing in 1e-6 all over the place, just like the article did. (Hey I do it too!) But if the problem hasn’t been analyzed then an absolute error tolerance is just a bug waiting to happen. It might seem to work at first, but then catch and confuse someone as soon as the tests use bigger numbers. Or maybe worse, fail to catch a bug when they start using smaller numbers.
- tomsmeding 2y ago
  But then relative error is also not a panacea. If I compute 1 + 1e9, then producing 1e9 - 1 instead would fall within a relative error bound of 1e-6 easily. More generally, relative error works only if your computation scales "multiplicatively" from zero; if there's any additive component, it's suspect. Of course, as you say, absolute error is also crap in general: it's overly restrictive for large inputs and overly permissive for small ones. I'm not a numerics person, but I do end up needing to decide on something sensible for error bounds on computations sometimes. How does one do this properly? Interval arithmetic or something?
  - dahart 2y ago
    > More generally, relative error works only if your computation scales “multiplicatively” from zero; if there’s any additive component, it’s suspect. IEEE floating point inherently scales from zero, the absolute error in any computation is proportional to the magnitude of the input numbers, whether you’re adding or multiplying or doing something else. It’s the reason that subtracting two large numbers has a higher relative error than subtracting two small numbers, c.f. catastrophic cancellation. > How does one do this properly? There’s a little bit of an art to it, but you can start by noting the actual result of any accurate operation has a maximum error of 0.5 LSB (least significant bit) simply as a byproduct of having to store the result in 32 or 64 bits; essentially just think about every single math instruction being required to round the result so it can fit into a register. Now write an expression for your operation in terms of perfect math. If I’m adding two numbers it will look like a[+]b = (a + b)*(1 + e), where e is your 0.5 LSB epsilon value. For 32 bit float, e == +/- 1e^-24. In this case I differentiate between digital addition with a finite precision result, and perfect addition, using [+] for digital addition and + for perfect addition. This gets hairy and you need more tricks for anything complicated, but multiplying each and every operation by (1+e) is the first step. It quickly becomes apparent that the maximum error is bounded by |e| * (|a|+|b|) for addition or |e| * (|a| * |b|) for multiply… substitute whatever your operation is. When doing more complicated order-dependent error analysis, it’s helpful to use bounds and to allow error estimates to grow slightly in order to simplify expressions. This way you can prove the error is less than a certain expression, but the expression might be conservative. A 3d dot product is a good example to work though using (1+e). Typically it’s reasonable to drop e^2 terms, even though it will technically compromise your error bound proof by some minuscule amount. a[*]x [+] b[*]y [+] c[*]z = ((ax(1+e) + bx(1+e)) + cx(1+e))(1+e) = ((ax+axe + by+bye)(1+e) + cz+cze)(1+e) = (ax+by+(ax+by)e + axe+bye+(ax+by)ee + cz+cze)(1+e) = (ax + bx + 2e(ax+by) + e^2(ax+by) + cz+cze)(1+e) = ax + by + 2e(ax+by) + e^2(ax+by) + cz + cze + axe + bye + 2e^2(ax+by) + e^3(ax+by) + cze + cze^2 = ax+by+cz + 3e(ax+by) + 2e(cz) + 3e^2(ax+by) + e^2cz + e^3(ax+by) Now drop all the higher order terms of e. = ax+by+cz + 3e(ax+by) + 2e(cz) Now also notice that 2e|cz| <= 3e|cz|, so we can say the total error bound: <= (ax + by + cz) + 3e( |a||x| + |b||y| + |c||z| ) And despite the intermediate mess, this suddenly looks very conceptually simple and doesn’t depend on the order of operations. If the input values are all positive, then we can say the error is proportional to 3 times the magnitude of the dot product. And it’s logical too because we stacked 3 math operations, one multiply for each element of the sum and two adds. Sorry if that was way too much detail… I got carried away. :P I glossed over some topics, and there could be mistakes but that’s the gist. I’ve had to do this for my work on a much more complicated example, and it took a few tries. There is a good linear algebra book about this, I think called Accuracy and Stability of Numerical Algorithms (Nicholas Higham). The famous PBR graphics book by Pharr et. al. also talks about error estimation techniques.
boscillator 2y ago
100% agree that picking numeric tolerances can be tricky. It is especially tricky when writing a generic math library like you are. If your doing something more applied it can help to take limits from your domain. For the example in the blog, if you're using GPS to determine you're position on earth, you probably know how precise physics allows that answer to be, and you only need to test to that tolerance (or an order of magnitude more strict, to give some wiggle room.)
deleted 2y ago
[deleted]

essexedwards 2y ago

I've also been surprised many times by issues in numerical libraries. In addition to matrices with simple entries, I've found plenty of bugs just testing small matrices, with dimensions in {0,1,2,3,4}. Many libraries/routines fall over when the matrix is small, especially when one dimension is 0 or 1. Presently, I am working on cuSPARSE and I'm very keen to improve its testing and correctness. I would appreciate anything more you can share about bugs you've seen in cuSPARSE. Feel free to email me, eedwards at nvidia.com

AlotOfReading 2y ago
This is one of the reasons I argue that it's almost always better to prioritize speed and stability than accuracy specifically. No one actually knows what their thresholds are (including library authors), but the sky isn't falling despite that. Instabilities and nondeterminism will blow up a test suite pretty quickly though.
- phdelightful 2y ago
  Yeah, we really just try to come up with very loose bounds since the analysis is hard. Even so, it does occasionally stop us from getting things way way wrong.
- essexedwards 2y ago
  > No one actually knows what their thresholds are (including library authors) If low-level numerical libraries provided documentation for their accuracy guarantees, it would make it easier to develop software on top of those libraries. I think numerical libraries should be doing this, when possible. It's already common for special-function (e.g. sin, cos, sqrt) libraries to specify their accuracy in ULPs. It's less common for linear algebra libraries to specify their accuracy, but it's still quite doable for BLAS-like operations.
  - AlotOfReading 2y ago
    What I'm trying to convey is that the required accuracies for the application are what's unclear. To give an example of a case where accuracy matters, I regularly catch computational geometry folks writing code that branches differently on positive, negative, and 0 results. That application implies 0.5 ulp, which obviously doesn't match the actual implementation accuracy even if it's properly specified, so there's usually a follow-up conversation trying to understand what they really need and helping them achieve it.
RhysU 2y ago
Nx0, 0xN, and 0x0 matrices are great edge cases. 0-length vectors, too. Then, do the same with Nx1, 1xN, and 1x1 matrices.

dahart 2y ago

> sometimes (usually) we just picked a magic tolerance out of thin air that seems to work. Probably worth mentioning that in general the tolerance should be relative error, not absolute, for floating point math. Absolute error tolerance should only be used when there’s a maximum limit on the magnitude of inputs, or the problem has been analyzed and understood. I know that doesn’t stop people from just throwing in 1e-6 all over the place, just like the article did. (Hey I do it too!) But if the problem hasn’t been analyzed then an absolute error tolerance is just a bug waiting to happen. It might seem to work at first, but then catch and confuse someone as soon as the tests use bigger numbers. Or maybe worse, fail to catch a bug when they start using smaller numbers.

tomsmeding 2y ago
But then relative error is also not a panacea. If I compute 1 + 1e9, then producing 1e9 - 1 instead would fall within a relative error bound of 1e-6 easily. More generally, relative error works only if your computation scales "multiplicatively" from zero; if there's any additive component, it's suspect. Of course, as you say, absolute error is also crap in general: it's overly restrictive for large inputs and overly permissive for small ones. I'm not a numerics person, but I do end up needing to decide on something sensible for error bounds on computations sometimes. How does one do this properly? Interval arithmetic or something?
- dahart 2y ago
  > More generally, relative error works only if your computation scales “multiplicatively” from zero; if there’s any additive component, it’s suspect. IEEE floating point inherently scales from zero, the absolute error in any computation is proportional to the magnitude of the input numbers, whether you’re adding or multiplying or doing something else. It’s the reason that subtracting two large numbers has a higher relative error than subtracting two small numbers, c.f. catastrophic cancellation. > How does one do this properly? There’s a little bit of an art to it, but you can start by noting the actual result of any accurate operation has a maximum error of 0.5 LSB (least significant bit) simply as a byproduct of having to store the result in 32 or 64 bits; essentially just think about every single math instruction being required to round the result so it can fit into a register. Now write an expression for your operation in terms of perfect math. If I’m adding two numbers it will look like a[+]b = (a + b)*(1 + e), where e is your 0.5 LSB epsilon value. For 32 bit float, e == +/- 1e^-24. In this case I differentiate between digital addition with a finite precision result, and perfect addition, using [+] for digital addition and + for perfect addition. This gets hairy and you need more tricks for anything complicated, but multiplying each and every operation by (1+e) is the first step. It quickly becomes apparent that the maximum error is bounded by |e| * (|a|+|b|) for addition or |e| * (|a| * |b|) for multiply… substitute whatever your operation is. When doing more complicated order-dependent error analysis, it’s helpful to use bounds and to allow error estimates to grow slightly in order to simplify expressions. This way you can prove the error is less than a certain expression, but the expression might be conservative. A 3d dot product is a good example to work though using (1+e). Typically it’s reasonable to drop e^2 terms, even though it will technically compromise your error bound proof by some minuscule amount. a[*]x [+] b[*]y [+] c[*]z = ((ax(1+e) + bx(1+e)) + cx(1+e))(1+e) = ((ax+axe + by+bye)(1+e) + cz+cze)(1+e) = (ax+by+(ax+by)e + axe+bye+(ax+by)ee + cz+cze)(1+e) = (ax + bx + 2e(ax+by) + e^2(ax+by) + cz+cze)(1+e) = ax + by + 2e(ax+by) + e^2(ax+by) + cz + cze + axe + bye + 2e^2(ax+by) + e^3(ax+by) + cze + cze^2 = ax+by+cz + 3e(ax+by) + 2e(cz) + 3e^2(ax+by) + e^2cz + e^3(ax+by) Now drop all the higher order terms of e. = ax+by+cz + 3e(ax+by) + 2e(cz) Now also notice that 2e|cz| <= 3e|cz|, so we can say the total error bound: <= (ax + by + cz) + 3e( |a||x| + |b||y| + |c||z| ) And despite the intermediate mess, this suddenly looks very conceptually simple and doesn’t depend on the order of operations. If the input values are all positive, then we can say the error is proportional to 3 times the magnitude of the dot product. And it’s logical too because we stacked 3 math operations, one multiply for each element of the sum and two adds. Sorry if that was way too much detail… I got carried away. :P I glossed over some topics, and there could be mistakes but that’s the gist. I’ve had to do this for my work on a much more complicated example, and it took a few tries. There is a good linear algebra book about this, I think called Accuracy and Stability of Numerical Algorithms (Nicholas Higham). The famous PBR graphics book by Pharr et. al. also talks about error estimation techniques.

boscillator 2y ago

100% agree that picking numeric tolerances can be tricky. It is especially tricky when writing a generic math library like you are. If your doing something more applied it can help to take limits from your domain. For the example in the blog, if you're using GPS to determine you're position on earth, you probably know how precise physics allows that answer to be, and you only need to test to that tolerance (or an order of magnitude more strict, to give some wiggle room.)

deleted 2y ago

[deleted]