Libertem: Numerical precision for high-resolution detectors

[uellue]

The test for numerical precision added in #226 uncovers that float32 may not be enough for 4096x4096 detectors to give sufficient numerical precision with standard reduce operations. We might have to use better algorithms like Kahan summation or calculate with float64.

For now, the test is "cheated" by using 1.0 as value, which doesn't give numerical errors.

[sk1p]

Relevant links about error boundaries in BLAS:

https://www.netlib.org/lapack/lug/node108.html
https://software.intel.com/en-us/mkl-developer-reference-c-error-analysis

[uellue]

Perhaps we can introduce a job option to perform the calculations with float64/complex128, i.e. convert the data in ApplyMasksTask.__call__()? That way the user can decide between fast and precise math. IMO this is better than Kahan summation because it should improve the situation for all reduction operations, not only the sum, and works with all libraries that support float64/complex128.

[uellue]

Another point is that the situation is particularly bad if all values are equal because all values have the same systematic error. If the values scatter, the situation should be much better.

Edit: Tested that, but couldn't get the test to pass. We should stick with the option for higher precision.

[ju-bar]

I propose to use butterfly summation. My tests show that it provides accuracy as well as better speed compared to straight summation.
Examples: fdncs2 (preserves array data) or fdncs2m (modifies array data)
in https://github.com/ju-bar/JMultiSliceLib/blob/master/src/JMultiSlice.cpp

[ju-bar]

I've used a small test program
https://github.com/ju-bar/TestSum
for checking different summation algorithms concerning accuracy and speed. Here is the result:
https://github.com/ju-bar/TestSum/blob/master/TestSum/output.txt

[uellue]

Hm, interesting! Thank you for testing this @ju-bar.

Since we use the summation only indirectly as part of the numpy or torch implementation of the dot product, it would have to be changed there. I'm a bit curious how that is implemented there and if the BLAS implementation or torch could actually be improved with a better sum.

Discussions with @sk1p: The most pragmatic solution for LiberTEM is probably passing an option to the dataset to convert the values to float32 or float64 (resp. complex64 or complex128) for users to decide between speed and accuracy. The dataset is the best place because, for example, K2IS has to be unpacked anyway and we can unpack right into the desired type, eliminating a copy step. Going to float64 rather than changing the sum algorithm will improve the accuracy of all calculations in all libraries, not only the sum in specific places.

[sk1p]

The most pragmatic solution for LiberTEM is probably passing an option to the dataset to convert the values to float32 or float64 (resp. complex64 or complex128) for users to decide between speed and accuracy.

In the low-level DataSet interfaces, this is now possible: you can specify your desired dtype for get_tiles(...). This parameter sets the dtype for the initial buffer inside of the dataset (conversion is done while reading the data); we just need to expose this to the user.

For mask application, we also need to specify the desired mask dtype (currently, the dataset determines the mask dtype, which is not always what we want). Together with the data dtype, we can then determine the result dtype (see also: np.result_type). I'm not sure about the efficiency of doing np.dot(f4, c8) vs. converting on reading and then doing np.dot(c8, c8) - in my tests they were surprisingly close.

[uellue]

While working on an implementation, I was running into an additional thing to consider: Since our masks are generated with factories, we don't know what dtype they have until they are generated. However, we do need this information on the client to allocate the result buffer, in case the masks contain complex values, for example.

As a solution, I'd propose to allow passing length and dtype as optional parameters with a mask job. If the user sets them properly, she can avoid generating the mask stack on the client just to check what is in it. Most notably, this will be the case for many analyses where we know exactly what's going to be in the mask stack. However, if the mask stack is quick to generate and the user likes convenience, the parameters can be omitted because they are autodetected from an ephemeral mask stack.

As a general design principle, we need a forward and a backward channel for dtype determination: We should allow user overrides, and we should allow autodetection from dataset and mask dtypes.

As far as I am aware, the dtype has to be considered at the following points:

• Allocation of result buffer
• Mask container, in particular the last conditioning step
• Dataset

Channels:

• Pass dtype to dataset (@sk1p 👍 )
• Read dtype from dataset (ds.dtype?)
• Pass dtype to mask container
• Read dtype from mask container

We should probably be able to query independently what the native dtype of a dataset or maks stack is, and what dtype it is configured to produce. A dataset should be able to produce different dtypes depending on what a job requests.

Here are some thoughts on what behavior could be correct:

• Float32 has proven to be a good default for most cases since it is fastest and high numerical precision is usually not a big issue. Integer masks should be converted some time before the mask slices are cached to reuse the converted version, source data probably already within the dataset implementation?
• Converting mask data from float64 to float32 should not throw an error since numpy often switches to float64 "by itself" when generating masks, without the user intending to use it. We should ignore float64 in the masks and use float32 by default.
• Float64 source data from detectors is uncommon. This might only appear as the result of conversions or simulations. Likely the precision and range is not really required.
• If either mask or dataset contain complex numbers, the default result type should be complex64, analogous to float32 for real numbers.
• Converting from complex to float should throw an error since that truncates the imaginary part, which is not what the user intended, usually.
• Users should be able to set float64 or complex128 as result type to work at higher numerical precision to address #227. What is the correct and most efficient method to switch the back-end to higher precision? Converting mask and/or data?
• Supporting integer result types could be interesting for electron counting detectors. This is significantly slower since numpy or torch don't have an optimized back-end for that. For that reason it should only be used if the user explicitly sets it. Otherwise both mask and data should be converted to floats. The dataset should already produce data with the correct dtype. This only makes sense if the mask and the data are integers. Otherwise this should throw an error.