SuiteSparse icon indicating copy to clipboard operation
SuiteSparse copied to clipboard

SPQR solve vs UMFPACK normal equations solve

Open mottelet opened this issue 11 months ago • 8 comments

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.

mottelet avatar Mar 08 '24 09:03 mottelet

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.

DrTimothyAldenDavis avatar Mar 11 '24 16:03 DrTimothyAldenDavis

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.

mottelet avatar Mar 11 '24 16:03 mottelet

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.

DrTimothyAldenDavis avatar Mar 11 '24 16:03 DrTimothyAldenDavis

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.

mottelet avatar Mar 11 '24 17:03 mottelet

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.

DrTimothyAldenDavis avatar Mar 11 '24 17:03 DrTimothyAldenDavis

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

mottelet avatar Mar 11 '24 17:03 mottelet

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.

DrTimothyAldenDavis avatar Mar 11 '24 17:03 DrTimothyAldenDavis

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!

mottelet avatar Mar 11 '24 17:03 mottelet