compute: Floating point summation performed in fixed point domain

[lluki]

What version of Materialize are you using?

v0.27.0-alpha.23.post (469bd67f9)

How did you install Materialize?

Built from source

What is the issue?

Floating point summation uses a fixed point accumulator that does not have double precision.

These are mostly follow up observations from #14695

The use of a fixed point accumulator results in the following observable problems:

Constant expression optimizer is not aware of limited summation domain

# Quick dance to avoid optimizations
materialize=> CREATE TABLE t (a int8);
CREATE TABLE
materialize=> INSERT INTO t VALUES (1);
INSERT 0 1

# This query will not be optimized away and go through the reduce operator. Wrong result.
materialize=> SELECT SUM(1e-39::double) FROM t;
 sum 
-----
   0
(1 row)

# This query will be optimized away and produces the correct result
materialize=> SELECT SUM(1e-39::double);
  sum  
-------
 1e-39
(1 row)

# Same again for really big values. Wrong result
materialize=> SELECT SUM(1e38::double) FROM t;
          sum           
------------------------
 1.0141204801825835e+31
(1 row)

# optimized away, correct.
materialize=> SELECT SUM(1e38::double);
  sum  
-------
 1e+38
(1 row)

materialize=>

Unexpected wrapping arithmetic

Normal floating point arithmetic over/underflows to +inf/-inf.

fn main() {                                                                                                                                                                                                                                                           
    println!("float overflow={:?}", f64::MAX + f64::MAX);                                                                                                                                                                                                             
}

=> float overflow=inf

On the other hand we wrap (note that 1e31 is essentially our max float):

# Materialize in release mode, we do wrapping arithmetic
materialize=> select * from d;
   a   
-------
 1e+31
 1e+31
(2 rows)

materialize=> select SUM(a) from d;
           sum           
-------------------------
 -2.8240960365167115e+29
(1 row)

Postgres in the same example produces an error:

# PSQL creates an error on overflow
lukas=> CREATE TABLE d (a double precision);
CREATE TABLE
lukas=> insert into d values (1e308),(1e308);
INSERT 0 2
lukas=> select sum(a) from d;
ERROR:  value out of range: overflow

Inconsistent use of wrapping arithmetic

Users expect that + behaves as SUM, but it does not, as the former uses float arithmetic while the later uses fixed point:

# Materialize in release mode, we do wrapping arithmetic
materialize=> select * from d;
   a   
-------
 1e+31
 1e+31
(2 rows)

materialize=> select SUM(a) from d;
           sum           
-------------------------
 -2.8240960365167115e+29
(1 row)

# normal addition is not affected
materialize=> CREATE TABLE d1 (a double, b double);
CREATE TABLE
materialize=> insert into d1 values (1e31, 1e31);
INSERT 0 1
materialize=> SELECT a + b from d1;
 ?column? 
----------
    2e+31
(1 row)

Effects on integer aggregation

Since bigint (uint8 and int8) perform their aggregation in double domain, the effects can be visible there as well. The following is planned as summation of t^2 in double domain. The plan is correct to account for the big intermediate values that might occur. However, due to our limited precision of double aggregation, this does not work correctly. In this example we calculate 9223372036854775807 ~= 2^63 squared which comfortable overflows our fixed precision aggregator. (math.log10((2**63)**2) ~= 38 and we support only 31)

materialize=> CREATE TABLE t (a int8);
CREATE TABLE
materialize=> INSERT INTO t VALUES (13), (32767), (2147483647), (9223372036854775807);
INSERT 0 4
materialize=> SELECT STDDEV(a) FROM t;
ERROR:  Evaluation error: cannot take square root of a negative number

Possible Fixes

1. Remove floating point from materialize. The aggregation is comparable to NUMERIC aggregation.
2. Use a double and work around the timely expectation of a group
3. Use a deterministic datatype that has higher precision than a machine double 3.1. Fixed point with higher precision + range than double 3.2. Use arbitrary precision 3.3. Use a reproducible float point accumulator

Thoughts for 3.3:

Note that reproducible accumulation does not mean the accumulator behaves like a a group, i.e.: accum_old = accum; accum += big; accum -= big; assert(accum_old == accum) is not guaranteed. But what is guaranteed is that the order additions will yield the bitwise same result (unlike a double). It's not clear to me if these guarantees are enough. It would certainly make each replica produce bitwise the same results, no matter the execution order.

• The tech report from Demmel, lot's of math that justifies the approach: https://www2.eecs.berkeley.edu/Pubs/TechRpts/2015/EECS-2015-229.pdf
• Demmel and H. D. Nguyen. “Parallel Reproducible Summation.” DOI: 10.1109/TC.2014.2345391: Introduces the parallel rsum.
• Here's a MIT licensed C++ implementation https://github.com/r-barnes/reproducible_floating_sums
• https://arxiv.org/pdf/1802.09883.pdf: Uses demmel in a RDBMS context, extends the demmel et al. to work with SIMD and gives achieves a fat speedup in (3x slowdown over normal float arithmetic). Unfortunately no source code available.

I think 3.2. or 3.3. are the way to go, they both offer the correct precision. 3.2 is faster to implement but slower (it might do memory allocation for example when summing). 3.3. is really nice, but the lack of a ready made implementation means some implementation effort.

Relevant log output

No response

[lluki]

See also this slack discussion: https://materializeinc.slack.com/archives/C02PPB50ZHS/p1664977591417439

[teskje]

While rug/gmp is not an option for use due to licensing, there appear to be some alternatives we could use:

• https://crates.io/crates/dashu
• https://crates.io/crates/astro-float

astro-float provides benchmark results that show it is significantly slower than rug, unfortunately: https://github.com/stencillogic/bigfloat-bench. I couldn't find anything on dashu but I wouldn't hold my breath. And we also don't know how much slower those libraries are than imprecise float math. But still, if we value correctness over efficiency, we should give those a try.

[tokenrove]

I would like to consider implementing/porting Ahrens's design; the accumulator is kind of complex but I think not actually that big a deal given that there are a variety of well-defined tests and two MIT-licensed implementations we can study.

[teskje]

If we can resolve the concern Lukas raised above, I agree that sounds like the best solution. If it turns out that implementing a Rust port is too hard, we could alternatively attempt a wrapper around the C++ port instead.

Note that reproducible accumulation does not mean the accumulator behaves like a a group, i.e.: accum_old = accum; accum += big; accum -= big; assert(accum_old == accum) is not guaranteed. But what is guaranteed is that the order additions will yield the bitwise same result (unlike a double). It's not clear to me if these guarantees are enough.

Looking at the code in reduce.rs ISTM that the accumulator only needs to be a semi group, so maybe we are good?

[tokenrove]

After discussing this in timely office hours, I think arbitrary precision floats would probably be fine, and less work than implementing Ahrens and making sure it works for all operations (though I would be happy to do this work). I had previously considered suggesting https://bellard.org/libbf/ in its infinite-precision mode, but then thought we needed a library for arbitrary precision rationals; however I suspect infinite precision floats are fine, too, on thinking about it. And I see that dashu supports rationals and floats anyway, so we should probably just use this.

[tokenrove]

Also worth noting that ReproBLAS itself could be wrapped, although it might be overkill consider how we'd normally use it; it's BSD licensed.