# Difficulty using LUDecomposition

I'm having trouble using LUDecomposition with pivoting. I read the Mathematica help on this particular command, but I'm still lost. Take a matrix like:

a = {{1, 2, 3}, {2, 4, 1}, {2, 5, 7}}


I evaluated

{lu, p, c} = LUDecomposition[{{1, 2, 3}, {2, 4, 1}, {2, 5, 7}}]


and got

{{{1, 2, 3}, {2, 1, 1}, {2, 0, -5}}, {1, 3, 2}, 0}


I don't really understand how to interpret this result, assuming that it's correct.

• You may also take a look at the last comment in this link – yshk Aug 1 '15 at 23:33

m = {{1, 2, 3}, {2, 4, 1}, {2, 5, 7}};
{lu, p, c} = LUDecomposition[m];
l = lu SparseArray[{i_, j_} /; j < i -> 1, {3, 3}] + IdentityMatrix[3];
u = lu SparseArray[{i_, j_} /; j >= i -> 1, {3, 3}];
l.u == m[[p]]
(* True *)


l.u is equal to a permutation of the rows of m

MatrixForm /@ {l, u}


• It'd be more idiomatic to use UpperTriangularize[] and LowerTriangularize[]. ;) – J. M. will be back soon Aug 1 '15 at 4:34
• @J.M. Not in my mindset, but feel free to edit or post another answer if you think there is something meaningful that I'm leaving aside – Dr. belisarius Aug 1 '15 at 4:48
• Alright, I might do so later… – J. M. will be back soon Aug 1 '15 at 5:00

As noted in the docs for LUDecomposition[], the two triangles are by default returned together as a single array; this is customary for LU decomposition routines, as in the original LINPACK and MATLAB's lu(). In fact, exactly this same format is stored internally by the LinearSolveFunction[] returned by LinearSolve[]:

a = {{1, 2, 3}, {2, 4, 1}, {2, 5, 7}};
{lu, piv, cond} = LUDecomposition[{{1, 2, 3}, {2, 4, 1}, {2, 5, 7}}]
{{{1, 2, 3}, {2, 1, 1}, {2, 0, -5}}, {1, 3, 2}, 1}

ls = LinearSolve[{{1, 2, 3}, {2, 4, 1}, {2, 5, 7}}];
ls[[2, 3]]
{{{1, 2, 3}, {2, 1, 1}, {2, 0, -5}}, {1, 3, 2}, 1}


You will need the first two results to reconstitute your original matrix. (The third result is a valuable diagnostic quantity called the matrix condition number; I should hope that this will eventually be introduced to you in your classes.)

In MATLAB, one could use the utility functions tril() and triu() to extract the triangular factors from the compressed array representing the LU decomposition. In Mathematica, on the other hand, belisarius has shown the way it used to be done. However, since version 7, Mathematica has been able to catch up to MATLAB in this regard, since the functions LowerTriangularize[] and UpperTriangularize[] became built-in:

l = LowerTriangularize[lu, -1] + IdentityMatrix[Length[lu]]
{{1, 0, 0}, {2, 1, 0}, {2, 0, 1}}

u = UpperTriangularize[lu]
{{1, 2, 3}, {0, 1, 1}, {0, 0, -5}}


To reconstitute the permutation matrix, do this:

p = IdentityMatrix[Length[lu]][[piv]]
{{1, 0, 0}, {0, 0, 1}, {0, 1, 0}}


Now, we can easily check the decomposition:

p.a - l.u
{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}}


As an aside, here is how one computes Det[a] using the results of LUDecomposition[a]:

{Signature[piv] Tr[lu, Times], Det[a]}
{5, 5}

• Ups, I understood that you meant to use Upper/LowerTriangularize instead of LUDecomposition. Ha! – Dr. belisarius Aug 1 '15 at 13:43

I don't know why people are complicating this more than it should be.

You did the LUDecomposition function and got this:

 {{{1, 2, 3}, {2, 1, 1}, {2, 0, -5}}, {1, 3, 2}, 1}


All you need to do is grab the first part of this answer, so the matrix:

{{1, 2, 3}, {2, 1, 1}, {2, 0, -5}}


To make it easier input it as an actual matrix form

From here it is super easy to see the answer once you know how to interpret it

For the Lower triangular matrix all you do is:

1. Keep the values in red in the lower part of the matrix
2. Erase the diagonal entries and everything above
3. Replace the diagonal entries by 1's
4. Everything else above by 0's

Like so:

As you can guess the Upper triangular matrix is similar to this:

1. You grab the same matrix with all the entries
2. This time you erase everything below the diagonal and replace it by 0's

Like so:

That is literally all you need to do, it seems long when I explain but this is such a quick way to do this.

This works everytime with all matrices.

• This is what both J.M. and belisarius did though. How would you write it in Mathematica code? – C. E. Jun 23 '17 at 12:22
• That is the whole point though. I am giving a new perspective on the same problem for some people who are just starting. It took me a while to get the hang of Mathematica codes and syntax which is unintuitive if you are learning alone. If you can simplify doing things on Mathematica it's very time saving. I only use the codes when I really have a restriction otherwise why would you? Someone just starting on Mathematica will have no clue what these codes are and most people explain poorly the sequence and relation of things. Look how beautifully this gives you the answer, why do you need code? – TheProteanGirl Jun 24 '17 at 5:32
• ok, I see where you are coming from now. There are many reasons why we assume people want to use the LU decomposition programmatically. Usually, the LU decomposition is just one step in a long chain of computations. Mathematica can be used to make all the computations for a given input. This makes it possible to experiment with lots of different inputs and get answers within milliseconds. With your approach, there would be copy and paste involved for getting each answer. Second, the values in the matrix may not be integers and you may lose precision if you read values off of the matrix form. – C. E. Jun 24 '17 at 9:37
• Third, the problems dealt with in Mathematica can be very large. Imagine trying to use your approach on a problem where the matrices as $10^4$ by $10^4$ elements. This is not unrealistic. So the reason is simply: people don't typically use Mathematica for simple problems where your approach is feasible. The problems are either much larger, require higher precision, or are not the one-off kind of problem. – C. E. Jun 24 '17 at 9:41
• I do remember also thinking that many answers poorly explain what they are doing when I was starting out with Mathematica. This becomes less and less of an issue the more you know about the Mathematica language. In this case, I find J.M.'s answer very readable, but then again I'm not a beginner anymore. – C. E. Jun 24 '17 at 9:51

# Alternative Implementation

Hidden deep in the documentation code is supplied for explicit LU-factorization.

DocumentationHelpLookup@"Compatibility/tutorial/LinearAlgebra/MatrixManipulation"


LUmatrices[mm_?MatrixQ] := Module[
{m},
m=First@LUDecomposition[mm];
{m - # + IdentityMatrix[Length@m], #} &[m*Array[Boole[# <= #2] &, Dimensions@m]]
];


A = {{1, 2, 3}, {2, 4, 1}, {2, 5, 7}};
MatrixForm /@ LUmatrices[a]
`

Here is the output for your matrix.

I do not know how efficient or general this implementation is.