makansij
makansij

Reputation: 9865

How to do complete lu factorization in matlab

I am not able to do complete pivoting in matlab. Say I have some matrix that is not sparse:

A = [0,1,1,1;0,1,0,0;1,1,1,1;0,0,0,-1]

When I try to use the lu command in matlab it does not like it because it is not sparse:

>> [L,U,P,Q] = lu(A) 
Error using lu
Too many output arguments.

Even the docs say that it must be sparse:

[L,U,P,Q] = lu(A) for sparse nonempty A, returns a unit lower triangular matrix L, an upper triangular matrix U, a row permutation matrix P, and a column reordering matrix Q, so that PAQ = L*U. If A is empty or not sparse, lu displays an error message. The statement lu(A,'matrix') returns identical output values.

I have two questions about this:

1) Why must it be sparse? In theory LU decomposition works for non-sparse matrices and so does complete pivoting.

2) What is the appropriate MATLAB method to call to do complete pivoting for non-sparse matrices?

Upvotes: 1

Views: 1834

Answers (3)

Eduardo
Eduardo

Reputation: 1

The DATATYPE must be sparse:

>> A = [0,1,1,1;0,1,0,0;1,1,1,1;0,0,0,-1];
>> A = sparse(A);
>> [L,U,P,Q] = lu(A) 

L =

(1,1) 1 (2,2) 1 (3,3) 1 (4,4) 1

U =

(1,1) 1 (1,2) 1 (2,2) 1 (1,3) 1 (2,3) 1 (3,3) 1 (1,4) 1 (2,4) 1 (4,4) -1

P =

(2,1) 1 (3,2) 1 (1,3) 1 (4,4) 1

Q =

(1,1) 1 (3,2) 1 (2,3) 1 (4,4) 1

Upvotes: 0

Riccardo Scorretti
Riccardo Scorretti

Reputation: 13

In theory LU factorization with full pivoting works fine with any matrix (full or sparse). In practice, it depends on the particular library which Matlab uses to compute it. In the book from Quarteroni et al. "Scientific computing with Matlab and Octave" (Springer, 4th edition, pag. 166) it is stated that if the matrix is full and has no particular structure, "a generic triangular factorization is computed by Gaussian elimination with partial pivoting".

My understanding is that in case your matrix is full, Matlab perform LU factorization by using an algorithm which uses only partial pivoting, thus matrix Q is not defined. If the matrix is sparse, a different library is used (UMFpack) which uses full pivoting. That's why matrix Q is returned only for sparse matrix.

Upvotes: 0

OmG
OmG

Reputation: 18838

As mentioned in comments the error is about the number of outputs. Q is for reordering columns which is specific for sparse matrices. Hence, you can get L, U, and P using the following command:

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

Upvotes: 1

Related Questions