Gauss-Jordan elimation unstable?

Question

Well for finding the inverse of any matrix (by hand) I learned to use the "Guass-Jordan elimination".

However today I was looking it up again on wolfram|alpha. And what struck me is the line:

The procedure is numerically unstable unless pivoting (exchanging rows and columns as appropriate) is used.

I do not fully understand that line, or the consequences of it. Does this mean that without row swapping the method may not find a solution? How so, why would it not find one, as the steps look pretty linear to me? - without going back to previous lines?

And what are the consequences of using this method in an automated system then?

Answer 1

You can find the solution, of course, and you will find the correct one if all the calculations are done exactly, ie symbolically.

But if you have something like $\sqrt 3 / 10$, the computer software will return an approximation precise to, say, 16 decimal digits. But an approximation means that there will be errors in the results.

If you do not pivot, these errors tend to "pile up", so that it is unstable and your end result may contain a sizable error and thus cannot be considered correct.

If you pivot, the error stay small, so while the solution won't be exactly right you can treat it as an excellent approximation of the "true" solution.

Why is it important to pivot? Because otherwise you will have catastrophic cancellation; take a look at http://en.wikipedia.org/wiki/Loss_of_significance

Answer 2

Consider the two linear systems: \begin{align*} A\mathbf{x}&=\mathbf{b}. \\ (A+\varepsilon_{ij})\mathbf{x}&=\mathbf{b}. \end{align*} where $\varepsilon_{ij}$ is a matrix of small perturbations/errors which changes the coefficients only very slightly. It turns out that the solutions to these two systems, whose parameters are very similar, can be very different indeed.

This comes from the numerical analysis principle which says that

If any calculation if undefined when a certain quantity is zero, the same calculation is likely to suffer numerical problems when the quantity is 'small.

This applies to Gaussian elimination when a pivot is non-zero but 'small'.

As an example, solve the following linear system using Gaussian elimination using four significant figures: $$\left(\begin{array}{cc}.0001& .5 \\ .4 & -.3\end{array}\right)\left(\begin{array}{cc}x_1\\ x_2\end{array}\right)=\left(\begin{array}{cc}.5\\ .1\end{array}\right).$$ You will get $y=1$ and $x=0$

The correct, exact solution to this is $$x=\frac{20000}{20003}\approx0.99975,\text{ and }y=\frac{19999}{20003}\approx 0.9998.$$

You might say, well just do exact arithmetic. In reality if you have a complicated linear system you are going to use a computer programme to solve the system and can computers handle exact arithmetic morryah infinite precision... Hence we need a way of avoiding these massive errors. Note that in the above example we had a very small pivot $0.0001$. Note how you would go from $0.4$ and $-0.3$ to $0$ and $-2000$: you would multiply by the small pivot $a_{11}=0.0001$ by the large multiple $m_1=10,000$ and then take $0.4r_1$ from $r_2$: \begin{align*} 0.4=a_{21}&\longrightarrow a_{21}-m_1a_{11} \\ -0.3=a_{22}&\longrightarrow a_{22}-m_1a_{12} \end{align*} A small pivot such as $0.0001$ can be small with respect to other elements of the row, such as $a_{12}=0.5$ and $b_1=0.5$, and this causes the elements $a_{22}$ and $b_2$ to `disappear' so we appear to be dealing with the system $$\left(\begin{array}{cc}.0001& .5 \\ .4 & 0\end{array}\right)\left(\begin{array}{cc}x_1\\ x_2\end{array}\right)=\left(\begin{array}{cc}.5\\ 0\end{array}\right).$$ which is clearly very different to the original one. Note that according to the principle above $a_{21}-m_1a_{11}\neq 0$ when $a_{11}=0$ so this should result in a problem when $a_{11}\approx0$, as it is here.

Hence the problem is with small pivots. Based on this an obvious idea would be to avoid unnecessarily small pivots. The way we do this is via a method called partial pivoting. Briefly, we make sure that the largest pivot (in magnitude: $-10$ is larger in magnitude than $0.4$ say) in any column or sub-column is the primary pivot. If this is not the case we swap rows to achieve this.