Unexpected Lower Triangular Matrix given by Matlab

Question

I was trying to find the lower triangular matrix L and upper triangular matrix U of an A matrix using the [L, U]= lu(A) command. For some reason, matlab won't give me a proper L matrix. Below are my original A matrix as well as the L matrix and U matrix produced by matlab

A = [4 2 -5 1;-8 0 9 7;-32 -4 43 18;24 4 -22 -8]

L = [-0.125 1 0 0;0.25 0.6 -0.2 1;1 0 0 0;-0.75 0.6 1 0]

U = [-32 -4 43 18;0 1.5 0.375 3.25;0 0 10 3.33;0 0 0 1]

What exactly is going on here?

Answer 1

The $\mathbf{LU}$ decomposition is not a decomposition of a matrix per se, but rather the decomposition of the linear system represented by that matrix. In order for this decomposition to be carried out, the equations will need to be ordered in a certain way. That is to say that the rows of the matrix representation of the system need to be ordered in a certain way. Remember that the matrix representation of a linear system is not unique. For instance,

$$ \left[ \begin{array}{cc} a_1 & a_2 \\ b_1 & b_2 \end{array} \right] \left[ \begin{array}{c} x \\ y \end{array} \right] = \left[ \begin{array}{c} c \\ d \end{array} \right] $$

is the exact same system as

$$ \left[ \begin{array}{cc} b_1 & b_2 \\ a_1 & a_2 \end{array} \right] \left[ \begin{array}{c} x \\ y \end{array} \right] = \left[ \begin{array}{c} d \\ c \end{array} \right]. $$

So we can permute the rows of a matrix and not alter the system as long as the right hand side is permuted in the same manner.

With the above mind, the $\mathbf{LU}$ decomposition of the linear system represented by the matrix $\mathbf{A}$ is given by

$$ \mathbf{PA} =\mathbf{LU}, $$

where $\mathbf{P}$ is the permutation matrix that puts the rows of $\mathbf{A}$ in the necessary order. What Matlab is returning when you only give two output arguments to lu() (i.e., [L,U]=lu(A)), is the combination of $\mathbf{P}^\mathrm{T} \mathbf{L}$ in the variable L. If you were to specify three output arguments instead ([L,U,P] = lu(A)), then you would get the upper and lower triangular matrices as well as the permutation matrix.

The documentation for lu would, in my opinion, be better worded of it used something other than L for the output argument for the [L,U] = lu(A) case. For instance, [T,U] = lu(A), where T = P'*L and P is the permutation matrix and L the lower triangular matrix.

Answer 2

See the first example in the documentation: L is a permuted lower triangular matrix (aka. psychologically lower triangular matrix).

Example 1

Start with

A = [ 1    2    3
      4    5    6
      7    8    0 ];

To see the LU factorization, call lu with two output arguments. [L1,U] = lu(A)

L1 =
    0.1429    1.0000         0
    0.5714    0.5000    1.0000
    1.0000         0         0

U =
    7.0000    8.0000         0
         0    0.8571    3.0000
         0         0    4.5000

Notice that L1 is a permutation of a lower triangular matrix: if you switch rows 2 and 3, and then switch rows 1 and 2, the resulting matrix is lower triangular and has 1s on the diagonal.