I have some code working to invert a sparse matrix:

Code: Select all

static bool sparseInvert(const Eigen::SparseMatrix<double> &mat, Eigen::SparseMatrix<double> &imat)
{
   Eigen::SparseLU< Eigen::SparseMatrix<double> > solver;
   solver.isSymmetric(true);
   solver.analyzePattern(mat);
   solver.factorize(mat);
   if (solver.info() != Eigen::Success) {
      mxLog("LU factorization failed");
      return true;
   }
   int size = mat.cols();
   Eigen::VectorXd bv;
   bv = Eigen::VectorXd::Zero(size);
   for (int cx=0; cx < size; ++cx) {
      Eigen::VectorXd xv;
      if (cx) bv(cx-1) = 0;
      bv(cx) = 1;
      xv = solver.solve(bv);
      for (int rx=0; rx < size; ++rx) {
         if (xv(rx) == 0) continue; // can produce 10x more coeff than input
         //if (fabs(xv(rx)) < 1e-4) continue; // doesn't work
         imat.coeffRef(rx, cx) = xv(rx);
      }
   }
   //mxLog("nonzero input %d output %d", mat.nonZeros(), imat.nonZeros());
   return false;
}


I am running some benchmarks now, but early results suggest this is much slower than BLAS dsytrf/dsytri for all but the sparsest of the sparse matrices. Am I doing something dumb that would slow things down? Is SparseLU not as fast as I had hoped?

Sparse solver are designed to handle large matrice (usually from 10000^2 to 10000000^2) having only a very few non-zeros per columns/rows (50 non-zero per column/row is already quite a lot). Moreover, you should never explicitly compute the inverse of a matrix: this is numerically unsafe, usually much more expensive, and even for very sparse matrices, the inverse might be completely dense thus breaking the advantage of a sparse storage.