NumPy QuadDType: Quadruple Precision for Everyone
labs.quansight.org
labs.quansight.org
In general, using a well-placed compensated (or exact) sum/dot product algorithm will solve most of your problems at a low cost. Nowadays, most people are familiar with Kahan summation (and related tricks like sorting numbers), but they do not know that it does not help much (often less than 128-bit floating points) and that the state of the art for those algorithms has moved quite a bit since. Nowadays, you can ensure that your sum/BLAS operation produces a result within one ulp of the exact result!
I highly recommend perusing the Handbook of Floating-Point Arithmetic (the Accurate Algorithm online page also seems fine [0]) to get an idea of what can be done. Implementations of these ideas can be found in most languages used for numerical computation (a quick Google search gave me options for C++/Python [1], Rust [2], etc.).
[0]: https://accurate-algorithms.readthedocs.io/en/latest/index.h...
It's interesting that 64 bits are often necessary (as opposed to 32) but 128 bits so rarely are.
1. What’s the implications when parallelization is considered? I.e. is those accurate sum algorithms parallelizable, and how would it compares to traditional parallelized sum? 2. What is the performance penalty I should be expecting?
2. It varies wildly depending on the algorithm you are using, the implementation of said algorithm, and the operations you are performing outside of the summation (if you sum the results of a loop performing non-trivial operations, that can easily overlap with the summation overhead making it essentially free). The only rule of thumb I can give is that runtime cost increases with desired accuracy. My recommendation would be to switch from a naive sum to using an exact accumulator and observe the difference in output value (letting you know how impactful numerical errors in that sum were) and runtime (letting you know how expensive the switch is for you), then dial it down if needed.
You will indeed need higher precision if you need a very large number of digits (more than what double precision gives you).
However, even with high precision, the algorithms you use to compute your moments matter: you might not have as many digits as you hope for. Variance, for example, is famously numerically unstable (sometimes giving a negative variance) when computed naively. A stable summation algorithm plugged into a variance computation (even a carefully written one) can sometimes give you six additional digits of precision in the variance that you might not get when doubling the precision (it depends on how big of a cancellation you had in your operations; too big, and doubling the precision will not be enough to recover, while a compensated summation would do fine).
I have seen this play out once in a clustering algorithm that would only work on some data when using a compensated summation to compute the variance.
You are also right, that if you can solve your problem with high precision, it's unlikely that just quad precision will be enough. For example the inverse Laplace transform. In some cases you need thousands of digits. The cases where double precision is not enough, but quad precision is are probably rare. But it's worth a try.
The article also gives an example from quantum physics. Some chaotic or highly numerically unstable problems just need quad precision to not blow up; primarily problems/algorithms (can be the problem itself or the algorithm to solve it that causes the issues) that tend to catastrophically amplify small errors. Matrix inversion is another example, to the point that most algorithms try to do almost anything to avoid explicitly inverting a matrix. But sometimes you just have to.
Tricks like clever factorization (which a lot of factorization algorithms have their own severe numerical issues e.g. some of the ways to compute QR or SVD), preconditioners, sparse and/or iterative algorithms like GMRES, randomized algorithms (my favorite) etc are all workarounds you wouldn't need if there was a fast and stable way to exactly invert any arbitrary non-singular matrix. Well, you would have less need for, there are other benefits to those methods but improving numerical stability by avoiding inverses is one of the main ones.
Do stats often deal with distict probabilities below 10^-300? And do they need to store them all over, or just in a couple variables that could do something clever?
Yes and with very wide dynamic range, though you really try to avoid it using other tricks. A lot of methods involve something resembling optimization of a likelihood function, which is often [naively] in the form of a product of a lot of probabilities (potentially hundreds or more) or might also involve the ratio of two very small numbers. Starting out far away from the optimum those probabilities are often very small, and even close to the optimum they can still be unreasonably small while still having extremely wide dynamic range. Usually when you really can't avoid it there's only a few operations where increased precision helps, but again even then it's usually a bandaid in lieu of a better algorithm or trick to avoid explicitly computing such small values. Still, I've had a few cases where it helped a bit.
Once you decide you need to go over 64 bits, the exact data type is largely a matter of convenience, and I can easily see two doubles being more convenient than either a quad or a 128-bit integer. You have more than 100 bits of precision in such a scheme, and it runs pretty fast.
If for e.g. you see bare summation rather than some sort of sorting by magnitude then expect bad results as numbers get bigger.
I looked up my consumer card on a whim (playing with sprinkling some GPGPU into my code—cupy is actually really easy to use and ergonomic!), and apparently doubles are 1/32 as fast as singles. Yikes!
For example: simulating the three body problem (three celestial bodies of equalish size orbiting eachother). Very small differences get amplified, and any difference at all in how the math is implemented on different systems will result in divergent solutions.
For some context, an electron "has structure" roughly on a length scale of about 10^-15 meters (details here unimportant). If you're representing a quantity in O(1) meters, then doing a handful of floating-point operations in double will likely accumulate error on the order of O(1) electrons. Incredible accuracy!
I think the key here is to be mindful of units and make sure you're combining like with like and scaling things appropriately as you go. If you're doing things the wrong way, quad will punish you eventually, too.
An underappreciated tool for working more confidently with floats is interval arithmetic. Worth checking out if you actually (truly) need to do things robustly (most people don't).
Proposals like Posits aim to improve on this :)
A great thing in the Posit proposal is that they also consider a distinct accumulator type for things like sums and dot product. In my experience (I did a PhD on measuring floating point error in large numerical simulations) those operations are the ones that are most likely to sneak on you and need fixing. But most people don't realize that those accumulator type are actually easy to use with IEEE floats, they are entirely orthogonal to using Posits.
Why introduce a new alias that might be 128 bits but also 80 ? IMO the world should focus on well defined types (f8, f16, f32, f64, f80, f128), then maybe add aliases.
If you have written a linear system solver, you might prefer to express yourself in terms of single/double precision. The user is responsible for knowing if their matrix can be represented in single precision (whatever that means, it is their matrix and their hardware after all), how well conditioned it is, all that stuff. You might rather care that you are working in single precision, and that there exists a double precision (with, it is assumed, hardware support, because GPUs can’t hurt us if we pretend they don’t exist) to do iterative refinement in if you need to.
Math functions are provided via libquadmath. Newer glibc also provides quad precision functions with slightly different names, IIRC using a f128 suffix rather than "q" like libquadmath.
Notably lacking support for ARM.
> Significand: 113 bits
> Exponent: 15 bits
Seems like a worse tradeoff than Dec128 [1] (110 bits, 17 bits) unless it's way faster.
[1] https://www.intel.com/content/www/us/en/developer/articles/t...
I was scratching my head for a bit there... It's also good to see that custom dtypes have come a long way. This is a great application of a custom dtype, and that wouldn't have been possible a few years back (well, okay, a "few" to me is probably like 8 years, but still).
Can the comparison also be done using the max precision unsigned integers instead of just float64?
I was just complaining about long double on Windows working on an algorithm that can benefit from higher precision, this post is serendipitous.