How can I implement the function lu(A) in MATLAB so that L*U is directly A and I also get the real L matrix?
When I use [L,U] = lu(A), MATLAB doesn't give me the right L matrix. When I use [L,U,P] = lu(A), I need to implement P*A = L*U, but I only want to multiply L*U to receive A.
MATLAB's lu always performs pivoting by default. If you had for example a diagonal coefficient that was equal to 0 when you tried to do the conventional LU decomposition algorithm, it will not work as the diagonal coefficients are required when performing the Gaussian elimination to create the upper triangular matrix U so you would get a divide by zero error. Pivoting is required to ensure that the decomposition is stable.
However, if you can guarantee that the diagonal coefficients of your matrix are non-zero, it is very simple but you will have to write this on your own. All you have to do is perform Gaussian elimination on the matrix and reduce the matrix into reduced echelon form. The result reduced echelon form matrix is U while the coefficients required to remove the lower triangular part of L in Gaussian elimination would be placed in the lower triangular half to make U.
Something like this could work, assuming your matrix is stored in A. Remember that I'm assuming a square matrix here. The implementation of the non-pivoting LU decomposition algorithm is placed in a MATLAB function file called lu_nopivot:
function [L, U] = lu_nopivot(A)
n = size(A, 1); % Obtain number of rows (should equal number of columns)
L = eye(n); % Start L off as identity and populate the lower triangular half slowly
for k = 1 : n
% For each row k, access columns from k+1 to the end and divide by
% the diagonal coefficient at A(k ,k)
L(k + 1 : n, k) = A(k + 1 : n, k) / A(k, k);
% For each row k+1 to the end, perform Gaussian elimination
% In the end, A will contain U
for l = k + 1 : n
A(l, :) = A(l, :) - L(l, k) * A(k, :);
end
end
U = A;
end
As a running example, suppose we have the following 3 x 3 matrix:
>> rng(123)
>> A = randi(10, 3, 3)
A =
7 6 10
3 8 7
3 5 5
Running the algorithm gives us:
>> [L,U] = lu_nopivot(A)
L =
1.0000 0 0
0.4286 1.0000 0
0.4286 0.4474 1.0000
U =
7.0000 6.0000 10.0000
0 5.4286 2.7143
0 0 -0.5000
Multiplying L and U together gives:
>> L*U
ans =
7 6 10
3 8 7
3 5 5
... which is the original matrix A.
You could use this hack (though as already mentioned, you might lose numerical stability):
[L, U] = lu(sparse(A), 0)
You might want to consider doing LDU decomposition instead of unpivoted LU. See, LU without pivoting is numerically unstable - even for matrices that are full rank and invertible. The simple algorithm provided above shows why - there is division by each diagonal element of the matrix involved. Thus, if there is a zero anywhere on the diagonal, decomposition fails, even though the matrix could still be non-singular.
Wikipedia talks a little about LDU decomposition here:
https://en.wikipedia.org/wiki/LU_decomposition#LDU_decomposition
without citing an algorithm. It cites the following textbook for proof of existence:
Horn, Roger A.; Johnson, Charles R. (1985), Matrix Analysis, Cambridge University Press, ISBN 978-0-521-38632-6. See Section 3.5.
LDU is guaranteed to exist (at least for an invertible matrix), it is numerically stable, and it is also unique (provided that both L and U are constrained to have unit elements on the diagonal).
Then, if for any reason "D" gets in your way, you can absorb the diagonal matrix D into either L (L:=LD) or U (U:=DU), or split it symmetrically between L and U (such as L:=L*sqrt(D) and U:=sqrt(D)*U), or however you want to do it. There is an infinite number of ways to split LDU into LU, and this is why LU decomposition is not unique.
