Skip to main content

Subsection 6.4.4 Numerical stability of linear solve via LU factorization with partial pivoting

The analysis of LU factorization without partial pivoting is related to that of LU factorization with partial pivoting as follows:

  • We have shown that LU factorization with partial pivoting is equivalent to the LU factorization without partial pivoting on a pre-permuted matrix: \(PA = LU \text{,}\) where \(P\) is a permutation matrix.

  • The permutation (exchanging of rows) doesn’t involve any floating point operations and therefore does not generate error.

It can therefore be argued that, as a result, the error that is accumulated is equivalent with or without partial pivoting

More slowly, what if we took the following approach to LU factorization with partial pivoting:

  • Compute the LU factorization with partial pivoting yielding the pivot matrix \(P \text{,}\) the unit lower triangular matrix \(L \text{,}\) and the upper triangular matrix \(U \text{.}\) In exact arithmetic this would mean these matrices are related by \(P A = L U.\)

  • In practice, no error exists in \(P \) (except that a wrong index of a row with which to pivot may result from roundoff error in the intermediate results in matrix \(A\)) and approximate factors \(\check L \) and \(\check U \) are computed.

  • If we now took the pivot matrix \(P \) and formed \(B = P A \) (without incurring error since rows are merely swapped) and then computed the LU factorization of \(B \text{,}\) then the computed \(L \) and \(U \) would equal exactly the \(\check L \) and \(\check U\) that resulted from computing the LU factorization with row pivoting with \(A \) in floating point arithmetic. Why? Because the exact same computations are performed although possibly with data that is temporarily in a different place in the matrix at the time of that computation.

  • We know that therefore \(\check L \) and \(\check U \) satisfy

    \begin{equation*} B + \Delta\!B = \check L \check U, \mbox{ where } \vert \Delta\!B \vert \leq \gamma_n \vert \check L \vert \vert \check U \vert . \text{.} \end{equation*}

We conclude that

\begin{equation*} P A + \Delta\!B = \check L \check U, \mbox{ where } \vert \Delta\!B \vert \leq \gamma_n \vert \check L \vert \vert \check U \vert \end{equation*}

or, equivalently,

\begin{equation*} P ( A + \Delta\!A) = \check L \check U, \mbox{ where } P \vert \Delta\!A \vert \leq \gamma_n \vert \check L \vert \vert \check U \vert \end{equation*}

where \(\Delta\!B = P \Delta\!A \) and we note that \(P \vert \Delta\!A \vert = \vert P \Delta\!A \vert \) (taking the absolute value of a matrix and then swapping rows yields the same matrix as when one first swaps the rows and then takes the absolute value).