Why does lu function yield different lower triangle matrix if I return [L,U] rather than [L, U, P]?

% square matrix A

A=[10,-7,0;-3,2,6;5,-1,5]

A = 3×3

10 -7 0 -3 2 6 5 -1 5

Return only L and U

[L1,U1] = lu(A);

Return L, U and P

[L2,U2,P2] = lu(A);

Compare L1 and L2

L1

L1 = 3×3

1.0000 0 0 -0.3000 -0.0400 1.0000 0.5000 1.0000 0

L2

L2 = 3×3

1.0000 0 0 0.5000 1.0000 0 -0.3000 -0.0400 1.0000

Answer 1

[Christine Tobler]

The LU decomposition really involves three new matrices: An upper-triangular matrix U, a lower-triangular matrix L, and a permutation matrix P. This is what the three-output syntax returns:

A = [10,-7,0;-3,2,6;5,-1,5]

A = 3×3

10 -7 0 -3 2 6 5 -1 5

[L, U, P] = lu(A)

L = 3×3

1.0000 0 0 0.5000 1.0000 0 -0.3000 -0.0400 1.0000

U = 3×3

10.0000 -7.0000 0 0 2.5000 5.0000 0 0 6.2000

P = 3×3

1 0 0 0 0 1 0 1 0

P*A

ans = 3×3

10 -7 0 5 -1 5 -3 2 6

L*U

ans = 3×3

10.0000 -7.0000 0 5.0000 -1.0000 5.0000 -3.0000 2.0000 6.0000

Unfortunately, lu also has a 2-output syntax. However, since it wouldn't be numerically safe to just compute L and U without a permutation matrix, internally we still compute all three matrices, and then return the first output as P'*L

[L2, U2] = lu(A);
L2

L2 = 3×3

1.0000 0 0 -0.3000 -0.0400 1.0000 0.5000 1.0000 0

P'*L

ans = 3×3

1.0000 0 0 -0.3000 -0.0400 1.0000 0.5000 1.0000 0

L2*U

ans = 3×3

10.0000 -7.0000 0 -3.0000 2.0000 6.0000 5.0000 -1.0000 5.0000

A = 3×3

10 -7 0 -3 2 6 5 -1 5

So the result of two-output LU satisfies A == L*U, but the output L isn't a lower-triangular matrix as one might expect.

You can argue that it would be better if the LU function didn't have a two-output syntax at all, but that decision was made a long time ago, and was probably based on the point that most people think of LU as a two-matrix decomposition, without thinking of the necessary permutation vector.

[John Taranto]

Thank you for the details.

Answer 2

[Steven Lord]

From the documentation page for the lu function:

"[L,U] = lu(A) returns an upper triangular matrix U and a matrix L, such that A = L*U. Here, L is a product of the inverse of the permutation matrix and a lower triangular matrix.

[L,U,P] = lu(A) returns an upper triangular matrix U, a lower triangular matrix L, and a permutation matrix P, such that P*A = L*U. The syntax lu(A,'matrix') is identical."

Emphasis added.

[John Taranto]

Thank you... I should have read the documentation page more closely before submitting the question.