Computing $\frac{x - y}{x - z}$ when $x,y,z$ are close to each other

Question

What is the most stable way to compute $$\frac{x - y}{x - z}$$ when $x$, $y$, and $z$ are all close to each other? I would like to compute expressions of this form in low precision on a GPU, but when I do so, I find that rounding errors really mess up the computation.

A straightforward implementation subtracts two numbers that are close to each other, and divides two numbers that are near zero. Neither of these operations is advisable.

I played around with the expression for a bit, trying to find a mathematical transformation that gives the same result without doing questionable numerical operations. However, couldn't figure out any good transformations.

Answer 1

If your inputs are $x,y,z$, this computation is not unstable, but ill-conditioned. That's worse, because it means that a small change in your input (such as a previous approximation as a floating-point number) causes a large change in your output.

In practical terms, this means that this is the wrong subproblem to focus on. If all you have is $x,y,z$ stored as floating point numbers, then you have already lost the battle.

Either you have some more properties on $x,y,z$ (such as them coming from a finite difference, see my question linked in a comment by NNN, or there is no hope of obtaining an accurate result.

Answer 2

Main source of instability here is that when all numbers are very close to each other,- expression approaches indeterminate form $\frac 00$. Add GPU rounding errors and it becomes obvious why you get unstable results. We can't get rid much of floating point rounding errors, because it's built-in into GPU floating point format, but we can try to get rid of indeterminate form mathematical uncertainty. For that try to compute expression $$ \frac {1}{1-z/x} - \frac {y}{x-z} ,$$

which is mathematically equivalent to your given one, but is more stable, because now you will not have indeterminate forms $0/0$, assuming $y \neq 0$ and $z \neq 0$.

As a test shot, I have tried such variables,

x=1.00010000521010  
y=1.00010000521020  
z=1.00010000521030

Then I have changed $x$ accordingly ,

|       x        |
|----------------|
|1.00010000521015|
|1.00010000521018|
|1.00010000521019|

Then I have noticed that relative error between steps of old and new expressions differs radically :

|   Err_Old (%)   |   Err_New (%)   |
|-----------------|-----------------|
|33.4072431633408 | 33.4814814817031|
|49.9999999999999 | 24.9999999997805|
|45.4545454545456 | 9.09090909099427|

Which shows that second expression converges gracefully, while original one- jumps in error and later re-considers back.

Note: do not expect black magic here, because these may be false result, due to not enough points taken (I have played with LibreCalc, so...). Besides results may depend on exact $x,y,z$ values taken and how they converges. In addition , GPU floating point inaccuracy may even spoil this small advantage. However, worth a try.

HTH!