2
$\begingroup$

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?

$\endgroup$
2
  • $\begingroup$ Is your last equation line correct in the first equality? I see $(x,y)-(\tilde x,\tilde y)=(-ϵ_4x,-ϵ_5y)$ so that the norm has an upper bound of $3ϵ_m\|(x,y)\|$ and the coefficient gets to be independent of $(x,y)$. $\endgroup$ Nov 3, 2021 at 8:23
  • $\begingroup$ @LutzLehmann Thanks! $\endgroup$
    – user56202
    Jun 30, 2022 at 18:37

1 Answer 1

6
$\begingroup$

(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|$$

$\endgroup$
5
  • 1
    $\begingroup$ This is very clear and easy to follow. Note that the relative errors $\epsilon_4$ and $\epsilon_5$ can be negative, so technically an absolute value sign is needed. $\endgroup$ Nov 3, 2021 at 13:23
  • $\begingroup$ Thank you for the very detailed response. $\endgroup$
    – user56202
    Nov 3, 2021 at 13:47
  • 1
    $\begingroup$ I think using the 2-norm would save us some work, since we don't have to split into cases. $\endgroup$
    – user56202
    Nov 3, 2021 at 13:55
  • $\begingroup$ @CarlChristian Good observation. I corrected it now. Thank you :) $\endgroup$ Nov 4, 2021 at 1:37
  • 1
    $\begingroup$ @user56202 Even the 2-norm proof will have some sort of $\max$ in it. We, mathematicians, have a dangerous but simple trick to avoid checking cases: "w.l.o.g. assume that $\max\{a,b\} = a$" :D $\endgroup$ Nov 4, 2021 at 1:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.