1
$\begingroup$

I have come across the following rectangular 5 x 10 matrix and carried out a LU decomposition of it, in the form PA = LU. The following matrices were obtained by function scipy.linalg.lu from module scipy to Python language. At first I thought was a computational issue, but after a contributor obtained the same result in both matlab and R, I decided to post it here (See: https://stackoverflow.com/questions/25186560/is-the-upper-triangular-matrix-in-function-scipy-linalg-lu-always-in-row-echelon)

\begin{equation} A = \begin{bmatrix} 1 & 1 & 1 & 1 & 0 & 1 & 1 & 1 & 1 & 0 \\ 1 & 1 & 0 & 0 & 1 & 0 & 1 & 0 & 1 & 1 \\ 1 & 1 & 0 & 0 & 0 & 1 & 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 1 & 1 & 0 & 1 & 0 & 0 & 1 \\ 1 & 1 & 0 & 0 & 1 & 1 & 1 & 1 & 1 & 1 \end{bmatrix} \end{equation}

\begin{equation} U = \begin{bmatrix} 1 & 1 & 1 & 1 & 0 & 1 & 1 & 1 & 1 & 0 \\ 0 & 1 & 0 & 1 & 1 & 0 & 1 & 0 & 0 & 1 \\ 0 & 0 & -1 & -1 & 0 & 0 & 0 & -1 & -1 & 0 \\ 0 & 0 & 0 & 0 & 1 & -1 & 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & 1 & 1 \end{bmatrix} \end{equation}

Notice that U is upper triangular, though it is not in row echelon form. Element U[5,5] = 1, but it should be zero, since the pivot for row 4 is immediately above it. According to Strang (1988, p.72), the U factor of a LU decomposition of a rectangular matrix should be in row echelon form, with the following characteristics:

"(i) The nonzero rows come first - otherwise there would have been row exchanges - and pivots are the first nonzero entries in those rows. (ii) Below each pivot is a column of zeros, obtained by elimination. (iii) Each pivot lies to the right of the pivot in the row above. This produces the staircase pattern."

STRANG, G. Linear Algebra and its applications. 3rd Ed., Thomson, 1988.

The L factor and permutation matrix P are below. Notice in P that Gaussian elimination swapped rows 2 and 4.

\begin{equation} L = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 1 & 0 & 1 & 0 & 0 \\ 1 & 0 & 1 & 1 & 0 \\ 1 & 0 & 1 & 0 & 1 \end{bmatrix} \end{equation}

\begin{equation} P = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ \end{bmatrix} \end{equation}

Say one tries to cure the problem by cancelling out element U[5,5] = 1, by subtracting row 5 from row 4, and finally producing the row echelon form. L will not be modified, since L[5,5] will remain equals 1, and $PA \neq LU$. I suspect that Gaussian elimination does not always produce the row echelon form.

$\endgroup$
  • $\begingroup$ If we write the decomposition you've already done as $A=PLU$ (i.e. our "L" is lower triangular up to permutation of rows), then clearing below the last pivot amounts to writing $A=PLE^{-1}EU$, where $E$ is the elementary matrix representing the row operation, which is the identity plus a matrix which is zero everywhere except at [5,4] and -1 at [5,4]. Then $E^{-1}$ is the identity except for a 1 at [5,4], which inserts a 1 into (PL)[5,4]. So the next step does change $L$. $\endgroup$ – Ian Aug 8 '14 at 17:42
  • $\begingroup$ Being clearer, we have $A=\begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ 1 & 0 & 1 & 1 & 0 \\ 1 & 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 1 & 0 & 1 & 1 & 1 \end{bmatrix} \begin{bmatrix} 1 & 1 & 1 & 1 & 0 & 1 & 1 & 1 & 1 & 0 \\ 0 & 1 & 0 & 1 & 1 & 0 & 1 & 0 & 0 & 1 \\ 0 & 0 & -1 & -1 & 0 & 0 & 0 & -1 & -1 & 0 \\ 0 & 0 & 0 & 0 & 1 & -1 & 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 & 0 \end{bmatrix}$ $\endgroup$ – Ian Aug 8 '14 at 17:47
1
$\begingroup$

According to Strang (1988, p.72), the U factor of a LU decomposition of a rectangular matrix should be in row echelon form,

What you can conclude from this is that in the 1988 book by Strang, the LU decomposition is defined so that $U$ is in row echelon form. By picking up another linear algebra book, you may find LU described differently. Also, none of those books constitute the documentation of numerical algebra routines.

The LAPACK routine dgetrf (which is probably used by the aforementioned software) promises the following output:

The factorization has the form $A = P L U$ where $P$ is a permutation matrix, $L$ is lower triangular with unit diagonal elements (lower trapezoidal if $m > n$), and $U$ is upper triangular (upper trapezoidal if $m < n$).

No promise is made to produce row-echelon form. A simpler example than yours is $$ A = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 \end{pmatrix} $$ for which the aforementioned software returns $U=A$.

This may happen for square matrices, too: e.g., $U=A$ is also returned for $$ A = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & 1 \end{pmatrix} $$

$\endgroup$
  • $\begingroup$ Thanks, I did not know there was no standard LU decomposition. I have always seen in textbooks the form described by Strang. $\endgroup$ – Anselmo Jan 18 '16 at 17:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.