I am computing the nullspace of a sparse rectangular $m$ x $n$ matrix $A$, where $m$ << $n$. I do this by computing the QR decomposition of $A^T$ and extract the $n-m$ right-most columns of the resulting $Q$. Symbolically: $$A^T=QR, Q=[Q_1 Q_2] , Q_2 = Nullspace(A)$$
Using this technique, I get the results I want but I would like it to be faster.
Here is the sparsity pattern of $A^T$ as generated by Mathematica (2,337 x 36,090; 42,900 nonzeros):
and the resulting nullspace (36,090 x 33,753; 8,167,054 nonzeros)
The sparsity pattern arises out of the constraints on the contact of two subdomains in a finite element problem. The coefficients on the left are positive whereas the coefficients on the right are negative. This is due to a sign convention. I ordered the node indices of each subdomain using a nodal connectivity matrix from the finite element mesh and fed that into METIS (which made my conjugate gradient solves MUCH faster!).
I primarily use Eigen and Intel MKL for my matrix algebra, but also use SuiteSparse for their much-faster QR decomposition routine. Although not important for this question, I am including my procedure for interfacing between Eigen and SuiteSparse, just in case it does matter. Here is how I compute the nullspace (kernel) of $A$:
using namespace Eigen; typedef SparseMatrix<double,ColMajor,int> Sparse_t; cholmod_common* cc, Common; cc = &Common; cholmod_l_start(cc); Sparse_t AT; //assume A^T is already defined //unfortunately necessary copy with SuiteSparse_long indices -- minimal impact on performance SparseMatrix<double, ColMajor, SuiteSparse_long> tmp(AT); const int rows = AT.rows(); const int cols = AT.cols(); const int ordering = SPQR_ORDERING_BEST; const double tol = SPQR_NO_TOL; const SuiteSparse_long econ = rows; const int getCTX = 1; cholmod_sparse AT_cm = viewAsCholmod(tmp); cholmod_sparse* I_cm = cholmod_l_speye(rows, rows, CHOLMOD_REAL, cc); cholmod_sparse* Q_cm = NULL; //perform QR factorization and request (Q^T * I)^T = Q SuiteSparseQR<double>(ordering, tol, econ, getCTX, &AT_cm, I_cm, NULL, &Q_cm, NULL, NULL, NULL, NULL, NULL, NULL, cc); Sparse_t kernel = viewAsEigen<Sparse_t::Scalar, Sparse_t::Options, Sparse_t::Index>(*Q_cm).rightCols(rows - cols); cholmod_l_free_sparse(&I_cm, &cc); cholmod_l_free_sparse(&Q_cm, &cc);
From my research on this website and other sources, the QR decomposition seems to be the best way to compute the nullspace of a sparse matrix. Furthermore, SuiteSparse's SPQR appears to be the fastest and most versatile routine for sparse matrices. I have compiled everything on Windows and linked SuiteSparse against Intel MKL (probably not relevant in this problem) and Intel TBB, and I do achieve limited multithreading while using SPQR.
I think my code is optimized for my platform of choice (Windows) on a six-core CPU laptop. I didn't have much luck getting the GPU accelerated features of SPQR to work, but that could be my weak GPU. However, my real question is if there is a better way of approaching this problem. Is there a linear algebra trick I am missing? Maybe some matrix decomposition I can perform on $A^T$ before performing the QR?
Thanks! (I wonder how long it will take someone to delete the pleasantries)
Update: I used umfpack to calculate the LU decomposition first, and then calculate the QR decomposition. $$P_rRAP_c=LU , U^T=QR , Q = [Q_1Q_2] , Nullspace(A) = P_cQ_2$$ I have not posted the additional code because it is a confusing c-style mess! The nullspace is computed in a fraction of the time (including the LU decomposition), but the resulting sparsity pattern is quite bizarre which appears to slow down subsequent operations. My numerical results appear to be unaffected. If $A$ is very sparse, there is negligible improvement in my overall finite element calculations. As $A$ becomes less sparse (more constraints), there is a very significant increase in performance. Maybe a different node ordering strategy would fix this.
Compared to the QR decomposition, the LU decomposition is instantaneous. Is there any way to compute the nullspace directly from a decomposition in the form of $P_rRAP_c=LU$?