Adding the Sparse Cholesky solvers

[adamConnerSax]

I'm new to Eigen and the haskell wrapper--which looks amazing!--so forgive me if this is an already answered and/or obvious question.

I see that the wrapper for Sparse solvers provides the LU and QR but not the simplical Cholesky solvers. Is that a licensing issue (since those are LGPL instead of MPL2), because they depend on external libraries (CHOLMOD ?) or something else?

I'd like to use them. I'm trying to replicate some of the functionality of Rs lme4 package and the simplical cholesky stuff is at the heart of it, at least for problems of any significant size.

I'd be happy to work on a PR--though I have no idea how hard it is!--or a separate library that has eigen as a dependency if that somehow helps with the licensing, dependency issues?

Anyway, any info would be much appreciated! I'm looking forward to using eigen!

[chessai]

Hey, thanks for your interest in this library!

Yes - this (core) library makes it a point to avoid anything LGPL. However, I have thought about a separate library (named eigen-lgpl?) that implements a lot of the LGPL stuff. If you'd like to take a stab at it, please go ahead!

[adamConnerSax]

Great! I'll probably start by just trying to get it working as a piece of the thing I'm working on but if it's straightforward enough, I'll roll the eigen parts out as a separate lib. And I'll try to adhere to your style, etc. when possible.

[chessai]

Thanks!

[adamConnerSax]

Okay. I'm trying to understand things and am deeply confused but that will probably get better.
My rough sense is that I need to have

• cbits/eigen-lgpl-sparse-la.cpp (with the typedefs and switches and c++ stubs for the FFI calls? More on this below)
• src/Eigen/Solver/SimplicalCholesky, with

1. the newtype for SimplicalLLT and SimplicalLDLT and
2. class instances for Solver, DirectSolver, and Code with an implementation assigning ints to the solver/ordering combos. See below.
3. Perhaps an additional function to match the vectorD class method for returning the diagonal matrix D in the LDLt decompositions.

• Perhaps a foreign import call... for vectorD in something like Eigen/LGPL/Internal.hsc

Some questions:

A. I need to add instances to those switch staements in eigen-sparse-la.cpp. Is there some best way to just add my cases (in cbits/eigen-lgpl-sparse-la.cpp) and refer the rest back to the function in the original library? I don't quite see how to do this since all the solver infrastructure, including that switch statement, is in eigen and I'm not looking to duplicate anything, just add the simplical solver cases.
B. I'm not quite sure what name or args vectorD ought to have in the in the Internal module. My best guess was
#api3 sparse_la_vectorD, "CSolverPtr -> Ptr (C a) -> CInt -> CInt -> IO String"
but only because I am wildly guessing that you handle vectors as Matrices with only one column or row and other functions returning a matrix seem to end in Ptr (C a) -> CInt -> Cint
C. Related to "B", I notice that the other functions returning decomposition components, like matrixL, etc. are unimplemented. Can that be fixed? I need to get the underlying L in the LLt decomposition. Or at least the log of its determinant.

Any help would be much appreciated here!

[adamConnerSax]

I feel like this is going to need to be more of a partial fork than an add-on. Because the solvers work via that switch statement, I need to reroute them all to a new switch to pick up the new cases and that's true of many functions. Hopefully I can fork only SparseLA and Internal.
Is there some way to re-design how eigen dispatches to the underlying c++ which would make it easier to add bits like this? Maybe it's just as simple as adding hooks into the switch statements which an add-on lib could use to add cases? I'm not clear enough yet...

[adamConnerSax]

And I think this might be a little beyond me, at least without more help. This fork idea seems bad, since any API change at the eigen level will mean re-writing the whole thing. But I don't see how else to do it.
And I'm going to need to get the actual Cholesky factor out as well as the permutation matrix that analyze computes. But there aren't any examples of that I can follow since the factor out functions matrixL, etc. are unimplemented and I see no examples which produce a dense vector (which is what the permutation matrix is).
Also, I need to be able to stack sparse matrices horizontally and vertically, something I don't see how to do in the c++ library or the haskell.
I'm sure all those things can be done. But having so much unknown makes it harder for me to guess how much work it might be to get to where it's useful for my particular problem.
If you could talk me through a bit of that, so I can see what needs to be done and what the right way to do it might be, that would help me decide if it makes sense to jump in.

[chessai]

I deeply apologise for taking so long to respond. This issue became active when I was predisposed, and then I never circled back. The current interface to the C++ is a pretty hacky, inherited blob of sinister lurking bugs. These bugs have been smashed out AFAIK, but like you say, the addition of more switch/case makes this marshalling's brittleness rear its ugly head. If you have a way to improve this situation, I would gladly take suggestions.

Furthermore, and this may affect your work on any of this, I've recently stumbled upon/played around with an interface to the statically verified/"safe" wrapper that this library provides that will significantly increase ergonomics. However, this will be a breaking change (This package doesn't follow PVP anyway, but still users might have a tough time migrating. The introduction of the first static interface was tough enough.)

Please reach out to me more on this issue and I will try hard to be more prompt. Sorry again.

[adamConnerSax]

No problem!

I've switched, for now, to interfacing directly with CHOLMOD--via hcholmod and some additions for sparse-linear-algebra.
But I think I'd rather work with a more complete sparse library. So maybe I'll circle back once you've settled on an interface via the "safe" wrapper.

Thanks for the work!