2
$\begingroup$

I am trying to solve problem $15.1$ from Numerical Linear Algebra by Trefethen and Bau, which reads

Determine whether the algorithm is backward stable, stable but not backward stable, or unstable.

Data: None. Solution: $e$, computed by summing $\sum_{j = 0}^{\infty} \frac {1}{j!}$ from left to right, and stopping when a summand is reached of magnitude < $\epsilon_{machine}.$

Now, the definition of backward stability given in the book is that for a problem $f$ with data space $X$, an algorithm $\tilde f$ is stable if for all $x$ there exists an $\tilde x$ such that $$||\tilde f(x) - f(\tilde x)|| = \mathcal O(\epsilon_m)$$ with $$||x - \tilde x|| / ||x|| = \mathcal O(\epsilon_m)$$

Since we have no "data", I am pretty confused on how to apply this definition. I am fairly certain that for a fixed $\epsilon_m$, $\tilde f(x)$ should be $\sum_{j = 0}^{n} \frac {1}{j!}$ where $n$ is the first summand less than $\epsilon_m$ (this has it's own problems, because the algorithm depends on $\epsilon_m$, which is weird). I am tempted to say $f(\tilde x)$ should be $e$ (even though it doesn't really match the definition) but that feels too much like forward-error analysis, which I know is distinct from stability. As far as interpreting $||x - \tilde x|| / ||x|| = \mathcal O(\epsilon_m)$, I am at a total loss.

I also thought of interpreting the problem space $X$ as all the set of all real number sequences, and $$f((a_k)_{k = 0}^{\infty}) = \sum_{j = 0}^{\infty}a_k.$$ since at least it would give a way to interpret $||x - \tilde x|| / ||x||$. But I don't think this is right either, because I don't think we should vary the sequence itself.

EDIT: In case it's relevant, here is the full exercise:

enter image description here

Axioms $(13.5)$ and $(13.5)$ are just standard axioms of floating point arithmetic. Namely, that for real $x$ in between the max and min value of the floating point system, $x$ can be rounded to some $\hat x = x(1 + \epsilon)$ with $|\epsilon| < \epsilon_{machine}$, and if $*$ is any of the $4$ arithmetic operations and $@$ is its floating point counterpart, then $x@y = x*y(1 + \epsilon)$ where $|\epsilon| < \epsilon_{machine}$.

The following are the definitions of stability and backward stability given in the book. I didn't see any definition for unstable, so I assume it just means "not stable".

Stability: enter image description here

Backward Stability: enter image description here

$\endgroup$
7
  • $\begingroup$ Could just be lack of coffee but I find the question confusing as well. I wonder if you could instead assess the stability of calculating $e^x$ for $x$ in some interval by summing $x^j / j!$? And then see what happens near $x = 1$... $\endgroup$ Nov 3 '21 at 17:34
  • $\begingroup$ @DanielShapero I omitted a few details in the problem statement that I really don't think are relevant. But just to make 100% sure, I added a picture of the exercise. The OP is about part d. $\endgroup$
    – user56202
    Nov 3 '21 at 18:27
  • $\begingroup$ Perhaps you mean "about part e"? Part d is a different problem. One might take the tack that since the data space is empty, any criterion defined by "for all $x$ in $X$" becomes vacuously true. $\endgroup$
    – hardmath
    Nov 3 '21 at 20:32
  • $\begingroup$ @hardmath Yes sorry, I meant part e. I think that we probably need to come up with a new definition of stability/backward stability that still captures its spirit. The spirit of backward stability, as famously put, is "an algorithm is backward stable if it gives the correct solution to almost the same question" $\endgroup$
    – user56202
    Nov 3 '21 at 21:04
  • 1
    $\begingroup$ I suspect the spirit of the book exercise is "Be sure to follow the definitions as given in the text." That is an example of the context available to you but not to Readers without additional effort in posting. $\endgroup$
    – hardmath
    Nov 3 '21 at 21:30
5
$\begingroup$

As hardmath mentioned in the comments above, if data space $D$ is empty, then the first half of the condition "for all $X\in D$ there exist $Y\in D$ such that $||X-Y||/||X||\approx\epsilon_m$ and $A(X)=P(Y)$" becomes vacuously true. So for backward stability analysis you only need to show the second part. In this case, it is simply proving that the algorithm gives you the best approximation to $e$ in the floating-point system. We can equivalently write this as $A_{l\in\{1,2\}} = \hat{e}$ where $\hat{e} = fl(e)$. (1 will be the left-to-right algorithm, 2 will be the right-to-left algorithm.)

I will partially solve the left-to-right case, because it demonstrates something. After it is demonstrated, right-to-left case is trivial.

First, we need to define the term "unit roundoff". Given a floating-point system $F$, the smallest $0<\mu \in F$ such $fl(1+\mu) > 1$ (where $fl$ is the rounding function) is called the unit roundoff. (Depending on the rounding function, unit roundoff and machine epsilon may have the same value, but usually unit roundoff is smaller than machine epsilon. The nitpicky difference is that $\epsilon_m = next(1)- 1$ where $next(1)$ is the smallest number in the floating-point system which is larger than 1).

Now, we write down the algorithm $\oplus_{j=0}^n \frac{1}{\otimes_{i=1}^j i}$ where $n$ is the smallest positive integer such that $1/n! < \epsilon_m$ ($n=18$ in double precision). For the sake of simplicity I will assume that the floating-point multiplication is exact (it will not be, and you should take that into account; but that is not the most important point in this exercise), so the algorithm becomes $\oplus_{j=0}^n \frac{1}{j!}$. Define $S_k = \oplus_{j=0}^k \frac{1}{j!}$, and compute a few terms:

$S_0 = 1$, $S_1 = 2$, $S_2 = 2.5$, $S_3 = 2.\bar{6}$, ... As you notice, $S_k$ is growing as $k$ grows, which is as expected. Also, $s_{k+1}$, the next term is getting smaller quite fast. This is where unit roundoff becomes important, I am adding an ever-decreasing number to an ever-increasing number. When is it going to stop making any difference? The answer turns out to be when $s_{k+1} < \mu S_k$ (left as an exercise to you, or you can just check the earlier chapters of the book). It is intuitively clear that $k$ where this condition occurs is less than or equal to $n$. More careful analysis would tell you that $k$ is strictly less than $n$ (in double precision, for example, $k=17$). This has a consequence, you cannot make that last necessary correction, which is necessary to obtain $\hat{e}$. You will be limited to $\hat{\hat{e}}$ which is different from $\hat{e}$. (The difference in double precision is $|\hat{\hat{e}}-\hat{e}|=4.4\times 10^{-16}$, while $\hat{e} = 2.718281828459045\color{red}09$ and $e=2.718281828459045\color{red}23536$, with difference ~$2.2\times 10^{-16}$).

As summary, under the assumption that the floating-point multiplication is exact, the first algorithm $A_1$ gives you a result where the error in the result is not optimal. Assuming $\hat{e}$ is the best approximation of $e$ in the floating-point system F, $A_1-\hat{e}\neq 0$ which means that it doesn't satisfy the definition. (You can also prove this rigorously, my explanation above is hand-wavy and experimental. But I was trying to demonstrate the point-of-failure to you -which is kind of obscured in rigorous arguments.)

$\endgroup$
12
  • $\begingroup$ Thank you so much, this is an amazing answer. I appreciate how you made everything concrete. I understand the technical points you made. But I am confused when, in the first paragraph, you say "you only need to show the second condition...". The definition of stability (at least the one given in my book, of which I uploaded a picture) is that an algorithm $\tilde f$ is stable if $$\forall x \in X \exists \tilde x \in X: ||\tilde f(x) - f(\tilde x)|| = \mathcal{O}(\epsilon_{machine}) \text{ and } ||x - \tilde x|| < \epsilon_{machine}.$$ This definition only as one condition, so (continued) $\endgroup$
    – user56202
    Nov 4 '21 at 6:15
  • $\begingroup$ (continued) shouldn't we just say the definition is vacuously true, and we are done? You say the second condition is "the algorithm gives you an answer close enough to $e$". Would you mind writing that condition more formally? Do you mean that a stable, or backward stable algorithm should give us $\text{fl}(e)$, the closest number in the floating point system $F$ to $e$? Lastly, you declare two algorithms $A_{l \in \{1, 2\}}$, but we only use one algorithm $A_1$. What did you mean for the second algorithm to be? $\endgroup$
    – user56202
    Nov 4 '21 at 6:20
  • $\begingroup$ Sorry, forgot to make the errors relative in my definition of stability. $\endgroup$
    – user56202
    Nov 4 '21 at 6:23
  • $\begingroup$ What I meant by the second condition was the non-vacuous part of the conjuction. I edited my answer to reflect this. $A_2$ is the algorithm when you sum from right-to-left, which I don't analyze in my answer. But using the observation, it should be easy to prove that $A_2=fl(e)$, hence backward stable. At this point, I should emphasize that my analysis of $A_1$ is not complete. I am barely arguing (in no way rigorous, wouldn't get any points in an exam) that $A_1$ is not backward stable; however, I am not commenting on whether it is stable or unstable. You should still do that work. $\endgroup$ Nov 4 '21 at 8:26
  • $\begingroup$ Thank you for your reply. The $\forall$ quantifier applies to the whole conjunction; that is, The definition is of the form $\forall x \in X [\exists \tilde x [P(x, \tilde x) \wedge Q(x, \tilde x)]]$, or more simply $\forall x \in X [\exists \tilde x R(x, \tilde x)].$ So as soon as $X$ is empty, it seems like the entire definition becomes true; I do not see a non-vacuous part. $\endgroup$
    – user56202
    Nov 4 '21 at 13:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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