LU decomposition of banded matrix with partial pivoting

Question

Disclaimer: I'm rusty as can be in this department.

I'm looking into how to implement a banded matrix LU decomposition with partial pivoting ($PA = LU$). So to start with I implemented regular matrix LU decomposition with partial pivot (code based on C example in https://en.wikipedia.org/wiki/LU_decomposition), and I fed it the following initial (banded) matrix:

$$ A^{(0)} = \begin{bmatrix} 5 & 8 & 0 & 0 & 0 \\\\ 35 & 64 & 4 & 0 & 0 \\\\ 0 & 40 & 22 & 9 & 0 \\\\ 0 & 0 & 18 & 85 & 7 \\\\ 0 & 0 & 0 & 8 & 19 \end{bmatrix} $$

So the half band width is 1, and I was hoping that the pivoting would keep the non-zeroes within a half band width of 2. But alas, this does not seem to be the case. See below what unfolds in working through the complete LU decomposition (note that matrix A iterations below contains L + U - the identity matrix, so the intermediate stages in the algorithm implementation):

$$ A^{(1)} = \begin{bmatrix} 35 & 64 & 4 & 0 & 0\\\\ 0.142857 & -1.142857 & -0.571429 & 0 & 0\\\\ 0 & 40 & 22 & 9 & 0\\\\ 0 & 0 & 18 & 85 & 7\\\\ 0 & 0 & 0 & 8 & 19\\\\ \end{bmatrix} $$ Above: row 0 got swapped with row 1. $$ A^{(2)} = \begin{bmatrix} 35 & 64 & 4 & 0 & 0\\\\ 0 & 40 & 22 & 9 & 0\\\\ 0.142857 & -0.028571 & 0.057143 & 0.257143 & 0\\\\ 0 & 0 & 18 & 85 & 7\\\\ 0 & 0 & 0 & 8 & 19\\\\ \end{bmatrix} $$ Above: row 1 got swapped with row 2. $$ A^{(3)} = \begin{bmatrix} 35 & 64 & 4 & 0 & 0\\\\ 0 & 40 & 22 & 9 & 0\\\\ 0 & 0 & 18 & 85 & 7\\\\ 0.142857 & -0.028571 & 0.003175 & -0.012698 & -0.022222\\\\ 0 & 0 & 0 & 8 & 19\\\\ \end{bmatrix} $$ Above: row 2 got swapped with row 3. $$ A^{(4)} = \begin{bmatrix} 35 & 64 & 4 & 0 & 0\\\\ 0 & 40 & 22 & 9 & 0\\\\ 0 & 0 & 18 & 85 & 7\\\\ 0 & 0 & 0 & 8 & 19\\\\ 0.142857 & -0.028571 & 0.003175 & -0.001587 & 0.007937\\\\ \end{bmatrix} $$ Above: final result: row 3 got swapped with row 4.

So the original row 0 got all the way pushed down to row 4, and you see the matrix is now not neatly banded anymore unfortunately.

My question is: am I overlooking something here? The result seems to be correct, if I extract L and U matrices from this, and multiply them, I get back the permutation of $A^{(0)}$.

I've seen a few articles hint that the propagation of non-zeroes is expected to be limited to a half band width that is double of the original half band width, see e.g. https://sep.stanford.edu/data/media/public/oldreports/sep20/20_11.pdf. But the example LU decomposition above seems to contradict this, and a row can propagate as far as it likes. Am I overlooking something?

Answer 1

I'm afraid you're not. I was also struggling around this question, as I was trying to implement a banded LU factorization with partial pivoting by myself. As a matter of fact, partial pivoting expands the bandwidth, but it's possible to make statements only on the bandwidth of the upper triangular part - this becomes equal to $p + q$, being $p$ the lower and $q$ the upper bandwidth, respectively.

Citing Golub and Van Loan (Chapter 4.3.3, emphasis mine):

Thus, pivoting destroys band structure in the sense that U becomes "wider" than A's upper triangle, while nothing at all can be said about the bandwidth of L. However, since the jth column of L is a permutation of the jth Gauss vector $\alpha_j$, it follows that L has at most p + 1 nonzero elements per column.

Golub, Gene H.; Van Loan, Charles F., Matrix computations., Johns Hopkins Studies in the Mathematical Sciences. Baltimore, MD: The Johns Hopkins University Press (ISBN 978-1-4214-0794-4/hbk; 978-1-4214-0859-0/ebook). xxi, 756 p. (2013). ZBL1268.65037.