3
$\begingroup$

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):

enter image description here

and the resulting nullspace (36,090 x 33,753; 8,167,054 nonzeros)

enter image description here

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.

Nullspace of $U^T

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$?

Improve this question
$\endgroup$
10
  • 1
    $\begingroup$ Are you sure that $A$ has full column rank? Note that the QR factorization without column permutation isn't reliable in detecting the rank of a matrix or giving you the null space if the matrix is rank deficient. $\endgroup$
    – Brian Borchers
    Commented Jun 23, 2020 at 18:04
  • $\begingroup$ There are also rank revealing LU factorization algorithms that you could apply to your matrix. $\endgroup$
    – Brian Borchers
    Commented Jun 23, 2020 at 18:05
  • 1
    $\begingroup$ The SuiteSparse QR factorization has column pivoting and is rank revealing, so that's OK. It sounds as though your constraints are unlikely to allow for a rank deficient matrix, but it's still a good idea to check that the rank is n-m. $\endgroup$
    – Brian Borchers
    Commented Jun 23, 2020 at 19:36
  • 1
    $\begingroup$ You don't need an orthogonal basis for this. Actually, there's little point in finding a basis for $N(A)$. Rather you could just write down the KKT conditions for the QP with $Ax=0$ constraints to get a linear system of $m+n$ equations in $m+n$ unknowns. Assuming that $M$ is positive definite, you can reduce this KKT system to a system of $m$ equations in $m$ unknowns and solve it by sparse Cholesky factorization. $\endgroup$
    – Brian Borchers
    Commented Jun 24, 2020 at 3:38
  • 1
    $\begingroup$ So at this point, your original problem was $\min x^{T}Mx$ subject to $Ax=0$, right? Why not revise this question to ask for quick ways of solving that problem. $\endgroup$
    – Brian Borchers
    Commented Jun 24, 2020 at 3:55

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.