Prove Solving a Lower Triangular Matrix By Forward Substitution is Backwards Stable I'm taking a class in scientific computing and we are working on proving stability of certain algorithms. Unfortunately, at this stage, everything is proof-based, and I have little to no experience in proofs. Any suggestions are greatly appreciated.
Let L be a nonsingular lower triangular matrix. How can I prove solving Lx = b by forward substitution is backwards stable?
I know forward substitution is defined by the algo:
$$
\ y_i = {l_{ii}}^{-1}* ( b_i-\sum^{i-1}_{j=1} l_{ij}*y_{i})
$$
and that backwards stability implies 
$$
\frac{||x - \tilde{x} ||}{||x||}
$$
where $$
\tilde{x} = $$ 
the computer approximation of some real x 
But I'm lost as to how to even get this started. Please advise!
P.S. this is my first interaction with LaTxT-- and I'm in love
 A: Although some parts of the question seem to be missing, here we go.
Algorithms in numerical linear algebra (worth adding the proper tag) are usually classified from the point of view of the numerical stability as forward and backward stable (excluding the exotic definitions like mixed stability etc.). Say, that for some data $x$ you want to compute the solution $y$ given by some function $f$: $y=f(x)$. Because our world (and computers as well) is not ideal, we get some $\tilde{y}$ instead of $y$. We distinguish here between the forward error, that is, some norm of $y-\tilde{y}$, and backward error. The backward error answers the question of how much we have to perturb the input data $x$ for which we actually solved the problem, that is, we look for $\tilde{x}$ such that $\tilde{y}=f(\tilde{x})$. This problem is not usually uniquely solvable (there might be many $\tilde{x}$ satisfying this condition), so one usually seeks $\tilde{x}$ which is close to the original data $x$.
Forward and backward errors are usually related using the perturbation theory associated with the problem you want to solve and the condition numbers comes into play here. In general, one has something like
$$
\text{forward error} \leq \text{condition number} \times \text{backward error}.
$$
Note, that the condition number does not necessarily mean the condition number of some matrix. This number depends on: what norms do you use and relatively to what do you measure them, on the data (and OK, possibly on some matrix, and the problem you solve (square nonsingular linear systems, least squares problems, the problem answered by 42 etc.).
Now to show that the forward substitution $Ly=b$ is backward stable, one can show that your computed solution $\tilde{y}$ satisfies $\tilde{L}\tilde{y}=b$ for some $\tilde{L}$ which is close to the original $L$; hopefully close up to something reasonably related to the machine precision used. Of course, since $b$ is also given (it is the data of the problem), one can also want to perturb that guy, but for the forward substitution there is actually really no gain in the final result.
In order to proceed, one needs a model of the finite precision arithmetic: e.g., assume that the result of $a~\mathrm{op}~b$ in the FPA satisfies:
$$\mathrm{fl}(a~\mathrm{op}~b)=(a~\mathrm{op}~b)(1+\delta), \qquad |\delta|\leq u,
\qquad\mathrm{op}\in\{+,-,*,/\}$$
where $u$ is usually called the unit roundoff and closely related to the so called machine precision (e.g., $u\approx 10^{-16}$ for the standard IEEE double precision floating point arithmetic).
In exact arithmetic, we have
$$
l_{ii}y_i = b_i-\sum_{j=1}^{i-1}l_{ij}y_j.
$$
Note that usually the step $i$ of the forward substitution is implemented as follows:

*

*set $w^{(1)}:=b_i$


*do $w^{(j+1)}:=w^{(j)}-l_{ij}y_j$ for $j=1,\ldots,i-1$,


*compute $y_i:=w^{(i)}/l_{ii}$.
Using the model of the FPA above, we have for $j=1,\ldots,i-1$:
$$
\tilde{w}^{(j+1)} = (\tilde{w}^{(j)}-l_{ij}\tilde{y}_j(1+\delta_j))(1+\delta_j'),
\qquad |\delta_j|\leq u, \quad |\delta_j'|\leq u.
$$
and
$$
l_{ii}\tilde{y}_i=(1+\delta_i)\tilde{w}^{(i)}, \qquad |\delta_i|\leq u.
$$
Putting this together gives
$$
\frac{l_{ii}\tilde{y}_i}{1+\delta_i}
=
b_i(1+\delta_1')\cdots(1+\delta_{i-1}')
-
\sum_{j=1}^{i-1} l_{ij}\tilde{y}_j(1+\delta_j)(1+\delta_j')\cdots(1+\delta_{i-1}')
$$
and dividing by the factor in the term containing $b_i$ leads to
$$
\frac{l_{ii}\tilde{y}_i}{(1+\delta_i)(1+\delta_1')\cdots(1+\delta_{i-1}')}
=
b_i
-
\sum_{j=1}^{i-1} l_{ij}\tilde{y}_j\frac{1+\delta_j}{(1+\delta_1')\cdots(1+\delta_{j-1}')}
$$
or by the proper definition of the components of a lower triangular matrix $\tilde{L}$:
$$
\tilde{l}_{ii}\tilde{y}_i
=
b_i
-
\sum_{j=1}^{i-1} \tilde{l}_{ij}\tilde{y}_j.
$$
Now it remains to bound the terms containing the $\delta's$. Although it is possible to make it more precise, one is usually happy with some first order expansion in terms of $u$ (to have a bound "$\text{something}\times u+O(u^2)$". In particular,
$$
\tilde{l}_{ii}=\frac{l_{ii}}{(1+\delta_i)(1+\delta_1')\cdots(1+\delta_{i-1}')}
=(1+\epsilon_{ii})l_{ii}, \qquad |\epsilon_{ii}|\leq iu+O(u^2)
$$
and
$$
\tilde{l}_{ij}=\frac{l_{ij}(1+\delta_j)}{(1+\delta_1')\cdots(1+\delta_{j-1}')}
=l_{ij}(1+\epsilon_{ij}), \qquad |\epsilon_{ij}|\leq ju+O(u^2).
$$
Since the factors by $u$ are all not greater than $n$, where $n$ is the dimension of the matrix $L$, this can be written in a more concise form
$$
(L+\Delta L)\tilde{y}=b, \qquad |\Delta L|\leq nu|L|+O(u^2),
$$
where both the absolute value and the inequality are to be understood component-wise.
