Search code examples
matlabmatrixlinear-algebrasparse-matrix

Sparse LU with Partial Pivoting

I would like to perform LU decomposition with partial pivoting on a sparse matrix. It seems that full pivoting is very fast and efficient for sparse matrices and the partial is not efficient for sparse matrices. My guess is that it is not supported or optimized for sparse.

A=randn(1e4).*(rand(1e4)<0.0001); 
S=sparse(A);

tic; [l,u,p]=lu(A); toc
Elapsed time is 8.699264 seconds.

tic; [l,u,p,q]=lu(S); toc
Elapsed time is 0.006430 seconds.

The second one, of full pivoting is extremely faster (by a factor of 1400)

My question is, how could it be? Shouldn't the partial pivoting LU be more efficient when the matrix is sparse, and always (or almost always) faster than the full pivoting?

Does anyone have an idea how can I perform fast LU with partial pivoting on sparse matrices?

Thanks, Gil

Solution

  • I feel like I have to clarify a few things:

    1. Are you aware that you are working on the full matrix in the [l u p] version, and not the sparse matrix?

    2. Partial pivoting is used to obtain numerical stability, not increase performance.

    3. Full pivoting is used to reduce the amount of fill that occurs when you factorize the sparse matrix (not possible for full matrices). By ordering both the rows and columns in an optimal way, the performance increases significantly.

    The reason for the increasing speed is because the non-zeros are alligned along the diagonal. However:

    "Sparse has to do with the structure of the entries: a band of zeroes around the diagonal, zeroes outside it." is not correct.

    A sparse matrix is a matrix that contains mostly zeros, and is represented three vectors (rows, cols and values).

    sparse() is a perfectly valid way of converting a full matrix into a sparse matrix. Duffymo's statement: "If it takes a full matrix and maps it into a sparse matrix data structure, then of course it's going to be slower.", is not correct, as long as the full matrix contains mostly zeros.

    Try the following:

    S = sprand(100,100,0.01);
[l, u, p] = lu(S);
spy(l)
figure
spy(u)

    Now, do the following:

    [ll, uu, pp, qq] = lu(S);
spy(ll);
figure
spy(uu);

    Look at the structures of l and u. You might get a clue why things are faster with full pivoting.

    Also, it is not really the factorizing part that is faster, it's the forward substitution part later on, when you try to solve this thing.