DrTimothyAldenDavis / SuiteSparse

The official SuiteSparse library: a suite of sparse matrix algorithms authored or co-authored by Tim Davis, Texas A&M University.
https://people.engr.tamu.edu/davis/suitesparse.html
Other
1.15k stars 259 forks source link

SPQR solve vs UMFPACK normal equations solve #774

Open mottelet opened 7 months ago

mottelet commented 7 months ago

Hello,

I tried to compare the performances of SPQR solve vs UMFPACK to solve a sparse least squares problem. The problem is a Laplacian mesh deformation application where we use 2D finite differences Laplacian on a square, leading to a (n+p) x n matrix A with a pentadiagonal n x n upper block and a p x n lower block of selected lines of the identity matrix. In the test case n=10000 and p=398. The test has been done on a mac mini M1 (first generation) with Scilab (www.scilab.org) which already has gateways to umfpack and I have used a homemade gateway to SPQR (via Eigen API). My results show that solving the normal equations with UMFPACK is roughly 10x faster than directly solving the least squares problem with SPQR. Did I miss something here ?

S.

DrTimothyAldenDavis commented 6 months ago

That is to be expected. QR factorization is intriniscally more costly than LU factorization.

When you are using the SPQR solve, are you getting Q as a regular matrix, or as a Householder representation? Q as a regular matrix could be very costlty.

mottelet commented 6 months ago

Hello

That is to be expected. QR factorization is intriniscally more costly than LU factorization.

When you are using the SPQR solve, are you getting Q as a regular matrix, or as a Householder representation? Q as a regular matrix could be very costlty.

yes, at least twice slower. But in the use case I just use the built-in "solve" method, so Q is not explicitely formed or used in the program. So I was trying to identify where could be the eventual bottlenecks, but Eigen:: methods are using Eigen::Maps to internal cholmod_sparse objects so there is not copy of any kind.

S.

DrTimothyAldenDavis commented 6 months ago

yes, at least twice slower.

That is for the dense case, not the sparse case. There is no simple 2x factor for the sparse case because of how fill-reducing orderings work, and because of how the sparse Householders work as supposed to sparse LU. A 10x difference is possible, depending on the sparsity pattern.

Since you're using the built-in SPQR solve, you won't have the issue for forming Q. It was worth a check though.

mottelet commented 6 months ago

That is for the dense case, not the sparse case. There is no simple 2x factor for the sparse case because of how fill-reducing orderings work, and because of how the sparse Householders work as supposed to sparse LU. A 10x difference is possible, depending on the sparsity pattern.

OK. As I said in my first message the matrix is just a finite differences of 2D Laplacian on a cartesian grid with additional rows of the identity matrix. I though that since A'*A would be less sparse normal equations + LU would be slower.

DrTimothyAldenDavis commented 6 months ago

There's no guarantee. The underlying ordering method, probably COLAMD, is seeing 2 different matrices. Since COLAMD is using a heuristic to solve and NP-hard problem, it can never guarantee anything.

mottelet commented 6 months ago

There's no guarantee. The underlying ordering method, probably COLAMD, is seeing 2 different matrices. Since COLAMD is using a heuristic to solve and NP-hard problem, it can never guarantee anything.

I used SPQR_ORDERING_BEST

DrTimothyAldenDavis commented 6 months ago

It's still just a heuristic. In addition, even if the ordering is optimal, sparse QR has more work to do that a sparse LU, in how it factorizes and assembles its frontal matrices.

mottelet commented 6 months ago

OK. Nevertheless, I notice an expected difference of 2 to 3 orders of magnitude when comparing the residuals of SPQR vs Normal Equations+UmfPack so everything is fine!