Implement CHOLMOD solver

[artofscience]

According to @aatmdelissen we may expect an increase in speed of about 5-10 x

https://pypi.org/project/scikit-sparse/

CHOLMOD cholesky

            if not hasattr(self, 'inv'):
                self.inv = cholmod.analyze(mat, ordering_method='best')
            self.inv.cholesky_inplace(mat)
            self.u = self.inv(rhs)
            self.adjointsolve = lambda b: self.inv(b)

[aatmdelissen]

So currently we are using cvxopt's implementation of CHOLMOD, which is also pretty fast. The 5-10 speedup is compared to scipy's spsolve. I did some benchmarking of the best cholesky solvers I could find, on a 300x100 mesh over 20 iterations. Average timings per solve:

• Scikit-sparse: 145±11 ms
• Intel Pardiso: 160±10 ms
• CVXOPT: 171±16 ms (only once factororize symbolic + factorrize numeric each iteration)
• CVXOPT: 226±6ms (linsolve each iteration)

It seems that scikit is fastest, but the difference is not very significant. CVXOPT is easiest to install out of the three, so I'd suggest sticking with that one. One easy improvement we could make is to do the symbolic factorization only once, and after that only do numeric factorizations, which gives about a 30% speedup.

[artofscience]

Arnoud, interesting comparison :) Sticking with cvxopt sounds like a good idea, especially if we also want to use it for solving the convex subproblems, then we do not have extra dependency. Could you share how to do the factorize symbolic + numeric?

[aatmdelissen]

This implementation should work

import cvxopt
import cvxopt.cholmod
class SolverSparseCholeskyCVXOPT:
    def update(self, A):
        Kcoo = A.tocoo()
        K = cvxopt.spmatrix(Kcoo.data, Kcoo.row.astype(int), Kcoo.col.astype(int))
        if not hasattr(self, 'inv'):
            self.inv = cvxopt.cholmod.symbolic(K)
        cvxopt.cholmod.numeric(K, self.inv)
        return self
    def solve(self, rhs):
        B = cvxopt.matrix(rhs)
        nrhs = 1 if rhs.ndim == 1 else rhs.shape[1]
        cvxopt.cholmod.solve(self.inv, B, nrhs=nrhs)
        return np.array(B).flatten() if nrhs == 1 else np.array(B)

[artofscience]

sweet :)