Why is subtraction a stable operation?

Question

In Numerical Linear Algebra by Trefethen & Bau, it is claimed that subtraction is backward stable. Here is the proof:

Let $f(x, y) = x-y$ and let $\tilde f(x,y)$ be the answer you get when doing $x-y$ on a computer. We have $$\begin{align*}\tilde f(x, y) &= (x(1+\epsilon_1) - y(1 + \epsilon_2))(1 + \epsilon_3)\\ \end{align*}$$ where $|\epsilon_1|, |\epsilon_2|, |\epsilon_3| < \epsilon_m = \epsilon_{machine}$. Now distributing, we have $$\begin{align*}(1+\epsilon_2)(1+\epsilon_3) &= 1 + \epsilon_2 + \epsilon_3 + \epsilon_2 \epsilon_3\\ &=1 + \epsilon_4 \end{align*}$$ where $|\epsilon_4| \le 2 \epsilon_m + \epsilon_m^2 \le 3 \epsilon_m$, if $\epsilon_m < 1.$ Similarly, $(1+\epsilon_2)(1+\epsilon_3) \le 1 + \epsilon_5$ where $\epsilon_5 \le 3 \epsilon_m$ if $\epsilon_m \le 1.$

Thus \begin{align*} \tilde f(x, y) = x(1 + \epsilon_4) - y(1 + \epsilon_5) = f(\tilde x, \tilde y). \end{align*}

Now here's they key part. Using the 2-norm, the relative error is $$\frac {||(x, y) - (\tilde x, \tilde y)||}{||(x, y)||} = \frac {||(\epsilon_4, \epsilon_5)||}{||(x, y)||} \le \frac {3 \sqrt2}{||(x-y)||} \epsilon_m = c(x, y) \epsilon_m$$ The problem is that the coefficient of $\epsilon_m$ above depends on the input $(x, y)$. There doesn't exist one $x$ that works for all $(x, y)$. But this is what is needed for backward stability. So why is subtraction backward stable?

Answer 1

(My apologies for the terrible typesetting)

Let me start with the definition of backward stability of an algorithm:

An algorithm $A$ to solve a problem $P$ is called backward stable, if for all $x$ in the input space there exists $y$ such that $||y-x||/||x||\approx \epsilon_m$ and $A(x) = P(y)$.

Trefethen & Bau summarizes this as "Backward stable algorithms give exact answers to nearly right questions" (or something like that, I don't have the book. I am looking at my lecture notes, and I am not a good note-taker).

Now, the problem of subtraction can be written as $P(X)=x_1-x_2$ where $X=(x_1,x_2)$. So you need to show that given the algorithm you described, i.e. $A(X) = \tilde{f}(x_1,x_2)$, for all $X$ in the input space there exists a $Y$ such that $||X-Y||/||X||\approx \epsilon_m$ and $A(X)=P(Y)$.

If we choose $||\cdot||$ as the $\infty-$norm, then the proof you provided above (except for the last step) gives us the answer:

$y_1 = x_1(1+\epsilon_4)$, which implies $y_1-x_1 = x_1\epsilon_4$ which in turn implies $ |y_1-x_1|/|x_1| = |\epsilon_4|$ (similar calculation for $y_2$, $x_2$). From this observation, we can conclude $$||X-Y||/||X||=\frac{\max\{|y_1-x_1|,|y_2-x_2|\}}{\max\{|x_1|,|x_2|\}}\approx \epsilon_m$$

To see why this holds, let's consider all 4 cases:

1. $\max\{|y_1-x_1|,|y_2-x_2|\} = |y_1-x_1|$ and $\max\{|x_1|,|x_2|\} = |x_1|$: $$||X-Y||/||X||=\frac{\max\{|y_1-x_1|,|y_2-x_2|\}}{\max\{|x_1|,|x_2|\}} = \frac{|y_1-x_1|}{|x_1|}=|\epsilon_4|$$

2. $\max\{|y_1-x_1|,|y_2-x_2|\} = |y_2-x_2|$ and $\max\{|x_1|,|x_2|\} = |x_2|$: $$||X-Y||/||X||=\frac{\max\{|y_1-x_1|,|y_2-x_2|\}}{\max\{|x_1|,|x_2|\}} = \frac{|y_2-x_2|}{|x_2|}=|\epsilon_5|$$

3. $\max\{|y_1-x_1|,|y_2-x_2|\} = |y_1-x_1|$ and $\max\{|x_1|,|x_2|\} = |x_2|$: $$||X-Y||/||X||=\frac{\max\{|y_1-x_1|,|y_2-x_2|\}}{\max\{|x_1|,|x_2|\}} = \frac{|y_1-x_1|}{|x_2|}\leq \frac{|y_1-x_1|}{|x_1|} =|\epsilon_4|$$

4. $\max\{|y_1-x_1|,|y_2-x_2|\} = |y_2-x_2|$ and $\max\{|x_1|,|x_2|\} = |x_1|$: $$||X-Y||/||X||=\frac{\max\{|y_1-x_1|,|y_2-x_2|\}}{\max\{|x_1|,|x_2|\}} = \frac{|y_2-x_2|}{|x_1|}\leq \frac{|y_2-x_2|}{|x_2|} =|\epsilon_5|$$