Skip to content

Pairwise summation for dx & dy: minor improvement in computational accuracy#7364

Open
halmaia wants to merge 1 commit intoOSGeo:mainfrom
halmaia:patch-1
Open

Pairwise summation for dx & dy: minor improvement in computational accuracy#7364
halmaia wants to merge 1 commit intoOSGeo:mainfrom
halmaia:patch-1

Conversation

@halmaia
Copy link
Copy Markdown

@halmaia halmaia commented May 2, 2026

The original
dx = ((c1 + c4 + c4 + c7) - (c3 + c6 + c6 + c9)) / H;
dy = ((c7 + c8 + c8 + c9) - (c1 + c2 + c2 + c3)) / V;
were replaced to
dx = ((c1 + c7) + (c4 + c4) - ((c3 + c9) + (c6 + c6))) / H;
dy = ((c7 + c9) + (c8 + c8) - ((c1 + c3) + (c2 + c2))) / V;
Pairwise summation reduces cancellation errors, so this form handles edge cases more reliably.
Minor cosmetics at variable "key".

…curacy.

The original
dx = ((c1 + c4 + c4 + c7) - (c3 + c6 + c6 + c9)) / H;
dy = ((c7 + c8 + c8 + c9) - (c1 + c2 + c2 + c3)) / V;
were replaced to
dx = ((c1 + c7) + (c4 + c4) - ((c3 + c9) + (c6 + c6))) / H;
dy = ((c7 + c9) + (c8 + c8) - ((c1 + c3) + (c2 + c2))) / V;
Pairwise summation reduces cancellation errors, so this form handles edge cases more reliably.
Minor cosmetics at variable "key".
@github-actions github-actions Bot added raster Related to raster data processing C Related code is in C module labels May 2, 2026
@echoix
Copy link
Copy Markdown
Member

echoix commented May 2, 2026

I have no opinion on if it works as good or not, so I would approve if no one has opinions in a couple days.

However, how does this compare to Kahan sums, like was used somewhere else in this project last year?
https://en.wikipedia.org/wiki/Kahan_summation_algorithm

@halmaia
Copy link
Copy Markdown
Author

halmaia commented May 2, 2026 via email


dx = ((c1 + c4 + c4 + c7) - (c3 + c6 + c6 + c9)) / H;
dy = ((c7 + c8 + c8 + c9) - (c1 + c2 + c2 + c3)) / V;
/* Pairwise summation reduces cancellation errors, this form handles edge cases more reliably. */
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[pre-commit] reported by reviewdog 🐶

Suggested change
/* Pairwise summation reduces cancellation errors, this form handles edge cases more reliably. */
/* Pairwise summation reduces cancellation errors, this form
* handles edge cases more reliably. */

@wenzeslaus
Copy link
Copy Markdown
Member

I'm generally in favor, but we want to be careful not to change results needlessly. Can you please comment on how this changes result and in which cases? CI with tests is all over the place now, but from what you are saying many computations might be just the same, or not? We do have an errata category where we put things like change of results due to better formulas or fixes. Does this fit? Do you have a minimal reproducible example?

@echoix
Copy link
Copy Markdown
Member

echoix commented May 2, 2026

The ubuntugis ppa is down. Is it related to the Ubuntu DDoS two days ago?

@metzm
Copy link
Copy Markdown
Contributor

metzm commented May 3, 2026

I am all for computational accuracy and numerical stability. In this case, input is elevation and values in a 3x3 window are summed up. With elevation, we usually do not have extreme values that could cause problems with double precision floating point limits. Moreover, the values in a 3x3 window should be similar to each other, not causing any cancellation errors when summing them up. Cancellation errors typically occur when adding a (absolute) very small number to a (absolute) very large number. This is highly unlikely in this case. Moreover, this PR replaces the summations of four values with the summations of 2x2 values. Pairwise summation is apparently meant for a much larger number of values.

OTOH, this PR does not do any harm and could indeed improve computational accuracy for edge cases. An example would be appreciated.

There are other occurrences of of summations of four or more values in r.slope.aspect/main.c. Shouldn't these also be addressed with added parentheses to use pairwise summations?

@metzm
Copy link
Copy Markdown
Contributor

metzm commented May 3, 2026

This PR is not only about summation, but summation and subtraction.

dx = ((c1 + c4 + c4 + c7) - (c3 + c6 + c6 + c9)) / H;

is without inner parentheses

dx = (c1 + c4 + c4 + c7 - c3 - c6 - c6 - c9) / H;

and re-ordered to only sum up differences which should all be similarly small

 dx = ((c1 - c3) + 2 * (c4 - c6) + (c7 - c9)) / H;

Just an example for various options to improve computational accuracy.

@metzm
Copy link
Copy Markdown
Contributor

metzm commented May 3, 2026

Theory says that the sum of differences is safer than the difference of sums. Thus

dx = ((c1 - c3) + 2 * (c4 - c6) + (c7 - c9)) / H;

is the proper solution, similar for dy. Please adjust accordingly.

Comment on lines +814 to +815
dx = ((c1 + c7) + (c4 + c4) - ((c3 + c9) + (c6 + c6))) / H;
dy = ((c7 + c9) + (c8 + c8) - ((c1 + c3) + (c2 + c2))) / V;
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Replace the difference of sums with the sum of differences.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

C Related code is in C module raster Related to raster data processing

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants