Backward Error Analysis

The result of the last theorem has been made use of by J.H. Wilkinson in a technique known as `backward error analysis'. In this, instead of bounding the rounding errors at each stage of a calculation, we find a perturbed equation whose exact solution is the same as the approximate solution we compute for the exact equation. Because of its importance, we give a brief outline of its application to Gaussian elimination.

Because of rounding errors, we do not obtain exact L,U factors; suppose we actually get L_c,U_c, which are such that

 

L_cU_c=A+E. (1.2)

Returning for a moment to the exact solution, we can think of solving $A\mbox{\boldmath$\space x $ } =LU\mbox{\boldmath$\space x $ } =\mbox{\boldmath$\space f $ }$ in two parts, the first of which is solving

$$L\mbox{\boldmath$ f $ } ^{(n)}=\mbox{\boldmath$ f $ } ,$$

corresponding to the reduction of the right hand side, and the second of which is solving

$$U\mbox{\boldmath$ x $ } =\mbox{\boldmath$ f $ } ^{(n)},$$

where L=[m_ij],  U=[a_ij⁽ⁿ⁾]. Because of rounding error, suppose we actually compute $\mbox{\boldmath$\space f $ } _c^{(n)}$ which satisfies

$$(L_c+\delta L)\mbox{\boldmath$ f $ } _c^{(n)}=\mbox{\boldmath$ f $ },$$ (1.3)

where we have used the already computed value of L. Suppose we also compute $\mbox{\boldmath$\space x $ } _c$, which is the exact solution of

$$(U_c+\delta U)\mbox{\boldmath$ x $ } _c=\mbox{\boldmath$ f $ } _c^{(n)}.$$ (1.4)

From (4.3) and (4.4), the computed solution satisfies

$$(L_c+\delta L)(U_c+\delta U)\mbox{\boldmath$ x $ } _c=\mbox{\boldmath$ f $ },$$

or

$$L_cU_c+\delta LU_c+L_c\delta U_c+\delta L\delta U)\mbox{\boldmath$ x $ } _c=\mbox{\boldmath$ f $ } .$$

We can write this in the form

$$(A+\delta A)(\mbox{\boldmath$ x $ } +\delta \mbox{\boldmath$ x $ } )=\mbox{\boldmath$ f $ }$$ (1.5)

by using (4.2) and then setting

 

$$\delta A E+\delta LU_c+L_c\delta U+\delta L\delta U$$ = (1.6)

$$\delta \mbox{\boldmath$ x $ } \mbox{\boldmath$ x $ } _c-\mbox{\boldmath$ x $ } .$$ =

We can now use Theorem 4.1.3 with $\delta \mbox{\boldmath$\space f $ } =\mbox{\bf0}$ to bound the error in the solution, provided we can bound $\left\Vert \kern.05em \delta A \kern.05em \right\Vert$. Taking norms in (4.6) and using norm inequalities, we find

$$\left\Vert \kern.05em \delta A \kern.05em \right\Vert\leq \le... ...ht\Vert.\left\Vert \kern.05em \delta U \kern.05em \right\Vert.$$

Wilkinson derives (at some length!) bounds for all these items in terms of computer arithmetic error bounds and with assumptions about the size of the matrix elements. We shall not go into that here, and in practice numerical bounds are rarely computed. However, in the process of finding a bound for $\left\Vert \kern.05em E \kern.05em \right\Vert$ it turns out that the errors are minimised if the multipliers are bounded by 1 in modulus. We thus use a pivoting strategy that ensures this, the commonest being to use partial pivoting, with the pivot of maximum modulus being chosen at each stage. It is rather disappointing that this clever analysis does not produce a rather more positive outcome.