Double floating precision exercise

Question

today I had to deal with this exercise:

If $x \approx y$, we might expect some cancellation in computing $log(x)- log(y)$. On the other hand, $log(x) - log(y) = log(\frac{x}{y})$, and the latter involves no cancellation. Does this mean that computing $log(\frac{x}{y})$ is likely to give a better result? (Hint: For what value is the log function sensitive (highly elastic)?)

I would assume that in double precision the second expression would be more accurate, but replicating the calculations in R did not show any differences between these two expressions.

x1 <- (1/0.999999999999)
y1 <- (1/0.99999999999999)

log(x1)-log(y1)
log((x1/y1))
1/log(1.00000001)

Answer 1

There is some inaccuracy in the result that is independent of how the result is computed. Basically, even if I did the computation to infinite precision, I would still get some error because the inputs were uncertain. The latter method may help when it is impractical to perform the operation with a bit of extra precision before rounding (as is common in multiple-precision libraries).

As far as independent error, say an estimate $\overline{x}$ of $x$ has a relative error of $\epsilon_x$ and an estimate $\overline{y}$ of $y$ has a relative error of $\epsilon_y$. Then we have $\overline{x}=x\,(1+\epsilon_x)$ and $\overline{y}=y\,(1+\epsilon_y)$. Note that these relative errors may be negative, indicating a smaller result than expected. Then the independent relative error $\epsilon$ in the result is found to be $$\begin{align*} \epsilon&=\frac{\log(\overline{x})-\log(\overline{y})}{\log(x)-\log(y)}-1\\ &=\frac{\log(x)-\log(y)+\log(1+\epsilon_x)-\log(1+\epsilon_y)}{\log(x)-\log(y)}-1\\ &=\frac{\log(1+\epsilon_x)-\log(1+\epsilon_y)}{\log(x)-\log(y)}\\ &=\frac{\log\left(\frac{1+\epsilon_x}{1+\epsilon_y}\right)}{\log\left(\frac{x}{y}\right)}=\frac{\log\left(1+\frac{\epsilon_x-\epsilon_y}{1+\epsilon_y}\right)}{\log\left(\frac{x}{y}\right)}\\ &\approx\frac{\epsilon_x-\epsilon_y}{\log\left(\frac{x}{y}\right)}\\ \end{align*}$$ So if $x$ and $y$ are similarly sized, the resulting relative error will be horrible no matter how much extra precision is added in computing the result (except if the two relative errors cancel out, but you can never plan on that).

The following Python code demonstrates the precision gain (double precision) from using division (checked against mpmath) as well as the independent loss of precision from using inexact values:

def relerr(xbar, x):
    return xbar/x-1

x,y = 1+1e-10,1+1e-10-1e-11
xbar,ybar = x*(1+1e-12),y*(1+1e-12-1e-14)
sub,div,l1p = log(xbar)-log(ybar),log(xbar/ybar),log1p((xbar-ybar)/ybar)
tsub,tdiv,tl1p = log(x)-log(y),log(x/y),log1p((x-y)/y)
mp.prec = 1000
for val in ['sub','div','l1p','tsub','tdiv','tl1p']:
    true = mp.log(mpf(x)/mpf(y))
    print(val+': '+mp.nstr(relerr(mpf(eval(val)),true)))

This prints

sub: 0.000999201
div: 0.000999201
l1p: 0.000999201
tsub: -8.9999e-11
tdiv: 0.0
tl1p: -9.0e-11

You will observe that rewriting using division resulted in zero measurable error compared to the 9e-11 introduced error from the subtraction. You can also see the large increase in relative error (from 1e-12 to almost 0.1%) when using inexact values that are close. If the relative errors have opposite signs, the effect is even worse.

Answer 2

Assume that every single operation has a relative error ≤ eps. Starting with this, what if x = 1.001 and y = 0.999, or x = 10.01 and y = 9.99, or x = 1,001,000,000 and y = 999,000,000? How is log (x / y) affected, and how is log (x) - log (y) affected?