@KofClubs is also working on this (or maybe the work has stalled?). Just FYI.

The Add method contains logic for processing the CounterResetHint, but the Sub method does not. Is this intentional, or could it be an oversight?

Yes, I know. (See the TODO here). I'm not sure
yet how to deal with negative counters in general (which wreaks havoc with the whole "only goes down if there is a counter reset" heuristics).

@KofClubs is also working on this (or maybe the work has stalled?). Just FYI.

I know, I guess the work has stalled as there hasn’t been activity for a long time. If that's not the case, then I'm ok to step aside and not interfere. Otherwise I'm ready to keep on

Thank you very much, I'll have a detailed look ASAP (which is probably not very soon, as I'm about to disappear into my holiday break). If anyone wants to review this in the meantime, be my guest.

I haven't heard anything from @KofClubs , so let's hope we don't have redundant work happening here.

Back from my holiday break, but swamped by the backlog. I hope I'll get to this next week.

@crush-on-anechka @beorn7

I just received an email about @crush-on-anechka 's good work. I'm sorry that I put this work on hold due to the long office hours every day:(

I shall unassign the corresponding issue.

Sorry, I caught a cold and got less done than planned. I need to ask for more patience. (And of course, if somebody else beats me to the review, be my guest.)

The Add method contains logic for processing the CounterResetHint, but the Sub method does not. Is this intentional, or could it be an oversight?

That's because I still haven't made up my mind how to treat counter resets for subtraction. I'll add that once I have an idea.

The reduceResolution function currently performs floating-point additions when merging buckets. Incorporating the Kahan summation algorithm here would add complexity but Is it worth it?

Let's not do it for now. I don't expect the counts in the merged buckets to be hugely different, and usually we'll only add a few buckets anyway.

A few more thoughts:

We currently don't do Kahan summation for the float binary + or - operator. That would require more plumbing in the PromQL engine, and I would also doubt it's terribly useful because it's unlikely to chain a lot of those operations in practical use.

Kahan summation is currently only used in certain aggregations, the corresponding "over time" functions and some specific functions. Many of those only act on float samples. The only aggregations and functions we will need Kahan summation for histograms for will be:

• sum aggregation
• avg aggregation
• sum_over_time function
• avg_over_time function

As far as I can see, we do not need KahanSub for any of those. So we actually do not need those parts.

@beorn7 thank you for the comments! I will be able to proceed in a week, if that's fine

The idea is to only switch to incremental calculation once we would everflow float64. An imprecise average is still better than the Inf or NaN we get when overflowing. Check how it is done for averaging the simple float samples.

After a series of experiments, I found that both incremental and non-incremental average calculation involving large-magnitude values do not behave as expected for simple floats also.

Summing numbers of drastically different magnitudes leads to precision loss. Kahan helps reduce errors, but it still relies on basic + and - operations, which themselves are subject to rounding when the values differ too much in scale.

For example we try to add 2.08169134891958e107 and 1.13729473572312e93.
Converted to the same exponent:
1.13729473572312e93 ≈ 0.0000000000000113729473572312e107
Since float64 holds ~16 significant digits, the smaller number is rounded to 0.00000000000001e107.

While reviewing the current summing algorithm, I noticed some potential issues and experimented with an updated version, nevertheless both of them don't work as expected in certain cases:

input: 1e+220 45 -1e+220
expected avg = 15

actual avg (incremental - current version) = 0
45 is lost - absorbed into 1e+220 and wasn't captured by the compensation term.

actual avg (incremental - updated version) = -2.0901902125471323e+203
on the final step (calculating the delta):
-1e+220 - 5e+219 = -1.4999999999999999e+220 (should be -1.5e+220)
This triggers a huge compensation term: -1.2541141275282797e+204, which corrupts the result.

input: 4.671249346736492e30 1e17 54 -1e17 -4.671249346736492e30 0
expected avg = 9

actual avg (incremental - current version) = 5.8640620148053336e+13
actual avg (incremental - updated version) = -7.237358535746467e+13
Although the values balance out by the end, a massive compensation term (13 orders of magnitude) remains and is applied to the final result, producing incorrect average.

input: 6.6666666666666666e107 1e53 45 -6.6666666666666666e107 -1e53
expected avg = 9

actual avg (incremental - current version) = -2e+52
1e53 and 45 are lost, so basically only the numbers with a magnitude around 1e107 are effectively taken into account, resulting in a final mean of 0.

actual avg (incremental - updated version) = 0
actual avg (non-incremental) = 0
45 is lost - absorbed into a compensation term ~52 orders of magnitude higher.

As you can see, non-incremental calculation also fails when the input values span three different orders of magnitude. If we only had values around ~100 orders of magnitude and others below ~14 orders of magnitude, the calculation would work: the main sum would handle the large numbers, while the compensation would account for the small ones. However, introducing values in the middle range (around ~50 orders of magnitude) corrupts the calculation.

TL;DR:
Are there any constraints or assumptions around the input data? Mixing values with such a wide range of magnitudes (e.g., 1e100, 1e50, and 45) is impossible to handle accurately using float64, whether in incremental or non-incremental calculation. Most likely, my test cases are much more extreme than what happens with real data. So far, I have more questions than answers — sorry about that!

Thanks for your thoughts.

There aren't really any assumptions about the input values. Users can do whatever they want. However, there is no expectation that we can solve all floating point precision issues. Even Kahan summation can take us only so far.

I'm not an expert in floating point arithmetics, but your line of argument that providing values from many different orders of magnitude will eventually break even Kahan summation, sounds plausible.

I usually try to imagine Kahan summation as "getting something like float128 precision even on a float64 FPU" – at somewhat higher cost, which doesn't really make a dent in the PromQL case because fetching all those numbers from the TSDB is much more expensive that the additional summation we have to do for Kahan. But even float128 will eventually run into precision issues.

About your examples:

1e+220 45 -1e+220 should just work fine with the current implementation. In fact, I tried it out, and it worked:

# Test summing and averaging extreme values.
clear

load 10s
	data{test="ten",point="a"} 1e+220
	data{test="ten",point="b"} 45
	data{test="ten",point="c"} -1e220

eval instant at 1m sum(data{test="ten"})
	{} 45

eval instant at 1m avg(data{test="ten"})
	{} 15

Where did you get your numbers from?

If you add other values with an order of magnitude in between, like the following, I can also see the problem:

load 10s
	data{test="ten",point="a"} 1e+220
	data{test="ten",point="b"} 45
	data{test="ten",point="c"} -1e220
	data{test="ten",point="d"} 1e22
	data{test="ten",point="e"} -1e22

But as said, I think that is expected. However, I expect the outcome to be close to zero here, not any gigantic number. Again, I don't know how you got the numbers you quoted above.

Apologies for the lack of clarity in my previous message. I’ve been experimenting with the avg_over_time function, and here’s how I set up the test:

clear
load 10s
  metric 1e+220 45 -1e+220

eval instant at 55s avg_over_time(metric[1m])
  {} 15

To force the incremental calculation path, I explicitly set incrementalMean = true.

When I referred to the “updated version” of the algorithm, I meant this one:

var delta float64
delta, kahanC = kahansum.Inc(-mean, f.F, kahanC)
delta /= count
kahanC /= count / (count - 1)
mean += delta

This is in contrast to the current version:

correctedMean := mean + kahanC
mean, kahanC = kahansum.Inc(f.F/count-correctedMean/count, mean, kahanC)

(Note for the reader to not lose track: We are currently discussion the average calculation for simple floats, not for histograms. We just stumbled upon problems with floats while trying to figure out how to do histograms.)

The current non-incremental version works:

load 10s
  metric14 1e+220 45 -1e+220

eval instant at 55s avg_over_time(metric14[1m])
  {} 15

The problem is if the test is performed at 1m because then the first value is not included (anymore, since v3). This is actually a problem in existing tests. I've created #16566 to fix that, please have a look. (We already had tests further down testing the above scenario.)

Even the incremental versions should give "correct" results within floating point accuracy, of course. I'm not sure if our current "hybrid incremental/Kahav avg calculation" is ideal (or maybe even correct). I found a bit of discussion on Stackoverflow, but it fell short to provide the actual algorithm that combines Kahan and incremental avg calculation.

I think your "updated version" is not quite correct. I have a slightly different idea. Maybe with that idea in place, we can get rid of the whole "switching to incremental calculation if needed" thing. I'm working on a PR to present to you for further discussion.

I have created #16569 for the (hopefully) correct combination of Kahan and incremental mean calulation.

#16569 looks tidy indeed! My thoughts:

• You’ve omitted the non-incremental calculation path, which actually performs better in cases where large positive and negative values should cancel each other out. Unfortunately, even with the Kahan summation improvement, the incremental approach can't handle these scenarios accurately due to float64 rounding limitations.
• I’m a bit puzzled because I haven’t found a test case where Kahan summation clearly outperforms direct summation in terms of accuracy. In fact, there is one example:

load 5s
  metric1 1e10 1e-6 1e-6 1e-6 1e-6 1e-6

eval instant at 55s avg_over_time(metric1[1m])
  {} 1.666666666666675e+09

This yields 1.6666666666666675e+09 with Kahan and 1.666666666666668e+09 with direct summation, but the difference is minimal and the test still passes as it's within the accepted tolerance.

I’ve written additional tests for the approach in #16569 and grouped them by behavior (see comments). metric9, metric10, and metric11 fail due to the issue of large positive and negative values that should cancel out. While metric9 and metric10 produce results that are at least close to expected, metric11 is significantly off.

clear
# All positive values with varying magnitudes
load 5s
  metric1 1e10 1e-6 1e-6 1e-6 1e-6 1e-6
  metric2 5.30921651659898 0.961118537914768 1.62091361305318 0.865089463758091 0.323055185914577 0.951811357687154
  metric3 1.78264e50 0.76342771 1.9592258 7.69805412 458.90154
  metric4 0.76342771 1.9592258 7.69805412 1.78264e50 458.90154
  metric5 1.78264E+50 0.76342771 1.9592258 2.258E+220 7.69805412 458.90154

# Contains negative values; result is dominated by a large-magnitude value
load 5s
  metric6 -1.78264E+50 0.76342771 1.9592258 2.258E+220 7.69805412 -458.90154
  metric7 -1.78264E+50 0.76342771 1.9592258 -2.258E+220 7.69805412 -458.90154
  metric8 -1.78264E+215 0.76342771 1.9592258 2.258E+220 7.69805412 -458.90154

# Contains negative values; large value dominates, but smaller values should cancel each other out
load 5s
  metric9 -1.78264E+215 0.76342771 1.9592258 2.258E+220 7.69805412 1.78264E+215 -458.90154
  metric10 -1.78264E+219 0.76342771 1.9592258 2.3757689E+217 -2.3757689E+217 2.258E+220 7.69805412 1.78264E+219 -458.90154

# High-magnitude values should cancel each other out, leaving a small result
load 5s
  metric11 -2.258E+220 -1.78264E+219 0.76342771 1.9592258 2.3757689E+217 -2.3757689E+217 2.258E+220 7.69805412 1.78264E+219 -458.90154

eval instant at 55s avg_over_time(metric1[1m])
  {} 1.6666666666666675e+09

eval instant at 55s avg_over_time(metric2[1m])
  {} 1.67186744582113
  
eval instant at 55s avg_over_time(metric3[1m])
  {} 3.56528E+49
  
eval instant at 55s avg_over_time(metric4[1m])
  {} 3.56528E+49
  
eval instant at 55s avg_over_time(metric5[1m])
  {} 3.76333333333333E+219
  
eval instant at 55s avg_over_time(metric6[1m])
  {} 3.76333333333333E+219
  
eval instant at 55s avg_over_time(metric7[1m])
  {} -3.76333333333333E+219
  
eval instant at 55s avg_over_time(metric8[1m])
  {} 3.76330362266667E+219
  
eval instant at 55s avg_over_time(metric9[1m])
  {} 3.228571428571429e+219
  
eval instant at 55s avg_over_time(metric10[1m])
  {} 2.5111e+219
  
eval instant at 55s avg_over_time(metric11[1m])
  {} -44.848083237

Thank you. I still have to process your latest response with the depth it deserves. However, I'm still caught between several events at work and in my usual life, but I will hopefully get to it next week.

Thanks for the additional test cases.

I have tried them out, and for all that fail with #16569, they also fail before #16569, just with more or less differing values. Thus, I'm still not convinced that the incremental approach with "correct" Kahan summation as in #16569 is really performing worse than direct Kahan summation with the one division at the end (what I call "simple mean calculation").

I will incorporate your test cases in #16569 and discuss them over there. (Hopefully I'll get to it tomorrow.) We can continue the detailed discussion over there.

About your other comment:

I’m a bit puzzled because I haven’t found a test case where Kahan summation clearly outperforms direct summation in terms of accuracy.

Not quite sure what you are referring to here. We have two dimensions here:

1. Kahan summation vs. regular summation.
2. Incremental mean calculation vs. simple mean calculation.

What do you mean by "direct summation"? There are plenty of test cases already in the test data that would fail without Kahan summation.

One more question about metric9 in your test cases: I assume your expected value has a typo. This series is completely dominated by 2.258E+220, so the mean should essentially be 2.258E+220 / 7 = 3.225714285714286e+219, which is what we see (and what I also get with Python's statistics package). I think you just typo'd in an 8 as the 3rd place after the dot.

For metric10, you said the expected value is 2.5111e+219. But again, this series is completely dominated by 2.258e+220, so the correct average should just be 2.258e+220 / 9 = 2.5088888888888888e+219. Interestingly, this is what happens after #16569 (and also with the Python statistics package), but before #16569, we get something fairly wrong (3.225714285714286e+219).

You're right — I suspect my earlier results were off due to inaccuracies in the averaging method I used. Using Python's statistics is a good idea, I will stick with it going forward.
Regarding your comment on metric10: how did you get 3.225714285714286e+219? It's not arbitrary but aligns with the value from metric9 (2.258e+220 / 7). It looks like the metric10 test case may have been trimmed from 9 to 7 entries. I actually get 2.5088888888888888e+219 in my version before #16569.

Regarding my comment on "Kahan summation outperforming direct summation" — what I actually meant was "Kahan summation vs. regular summation."

Since I was working on the funcAvgOverTime function and we encountered the incremental vs. non-incremental issue, my tests primarily focused on that. However I also ran the same algorithm without the Kahan compensation tracking and the results turned out to be identical.

I then tried to come up with scenarios where Kahan compensation would make a noticeable difference, but I couldn't find any clear-cut cases. This was just a side experiment driven by curiosity, and if there are already tests confirming the effectiveness of the Kahan implementation, I'm happy to disregard my earlier comment.

Regarding your comment on metric10: how did you get 3.225714285714286e+219? It's not arbitrary but aligns with the value from metric9 (2.258e+220 / 7). It looks like the metric10 test case may have been trimmed from 9 to 7 entries. I actually get 2.5088888888888888e+219 in my version before #16569.

You are right. I must have left something out of metric10 when I tried it with the pre-#16569 state. The test actually works fine. I'll remove the comment.

Oven in #16714, we also got some more evidence that the direct mean calculation is doing better than the incremental mean calculation. My plan is now to combine what we had before #16569 with the improved incremental mean calculation.

See #16773 .

Assuming that PR is now really what we want, we have the blueprint for how to do Kahan summation in native histograms. (Sadly, with the adaptive approach still being better than always just doing the incremental mean calculation, things are not as simple as I was hoping for.)

Sounds great! I’ll admit I got a bit lost while tweaking the Kahan algorithm, so I mostly ended up taking a backseat and observing. I’ll jump back into working on the histograms now

I've implemented #16773 avg calculation behavior for histograms. Here are a few notes from my side:

1. I’ve added several TODO(crush-on-anechka) markers where errors are currently not handled. The only error that may occur is ErrHistogramsIncompatibleSchema, and I’m not sure what the correct behavior should be in these cases.
2. While testing the funcAvgOverTime function, I encountered a confusing issue: histograms were being processed multiple times, which led to incorrect results. I noticed that this was due to in-place modification of the input histograms: when I introduced toAdd variable that copies and operates on the input histogram, the issue was gone. Although I didn’t observe this issue in the aggregation function, I applied the same toAdd pattern there for consistency.
3. I wrote new PromQL tests using values of different magnitudes and also tested switching to incremental calculation.

Sounds great. I'll have a look ASAP. (My review queue is quite manageable these days, but I will have limited availability over the next two weeks.)

And now I'm on vacation… sorry for all the delays. I'll be back in action for good starting August 12.

And now I'm on vacation… sorry for all the delays. I'll be back in action for good starting August 12.

That's totally fine. Do you think it's a good time to rebase onto the latest main locally to resolve conflicts?