julia and sparse matrix computation: is there plan?

[vav...@uwaterloo.ca]

Sparse matrix computation is a key aspect of scientific computing.  The creators of Julia made the smart decision to rely on SuiteSparse, high-quality free C++ software that carries out most of the common sparse matrix factorizations and related operations very efficiently.

The issue that concerns me is that the wrappers in Julia, namely, cholmod.jl, spqr.jl and umfpack.jl, expose only a small fraction of SuiteSparse's capabilities.  Furthermore, the exposed APIs give the Julia programmer significantly fewer capabilities than the sparse matrix suite of Matlab.

And the reason for this limitation has become clear to me over the past few days as I tried to develop a routine to find a right null vector of a sparse matrix using SuiteSparse: it is difficult and tedious to write Julia wrappers for the routines in SuiteSparse.  Even figuring out whether a particular variable is an Int32 or Int64 requires poring over SuiteSparse's nested header files full of conditional platform-dependent compilation directives.  And my resulting wrapper works only on one platform for one version of Julia and one version of SuiteSparse.  It is clear that whoever wrote cholmod.jl, spqr.jl and umfpack.jl put a lot of time and painstaking effort into getting these interfaces correct.

And this raises a question: will Julia ever fully incorporate SuiteSparse?  Will its sparse matrix capabilities ever approach Matlab's?  Who will take the time to write, test and maintain cross-platform Julia wrappers for SuiteSparse?  This is a tedious job that carries little glory - I certainly would not want to do it!  And yet sparse matrix computation is so important.

[Tony Kelman]

> Furthermore, the exposed APIs give the Julia programmer significantly fewer capabilities than the sparse matrix suite of Matlab.

Like what? There is a SuiteSparse.jl package under the JuliaSparse organization that currently has a handful of routines (Julia ports) but will be the home as more of the bindings gradually get migrated there, and out of base Julia. The long-term plan is to be less reliant on SuiteSparse as the only option, and instead to come up with a more general interface that multiple different solver packages could plug into. Outside of base Julia, a SuiteSparse.jl package would be more flexible to wrap as much of the SuiteSparse API as its contributors and maintainers choose.

[vav...@uwaterloo.ca]

Tony,

Here are a few Julia routines sparse matrix routines that I looked for and didn't find.  Maybe I looked in the wrong place.

- Getting the P'*L factor from SuiteSparse's sparse Cholesky back into Julia.  In other words, if F is the cholmod representation of the Cholesky factor, then sparse(F[:L]) works to retrieve the L factor, but sparse(F[:PtL]) does not work to retrieve P'*L (so I wrote my own wrapper).

- A general routine for computing a null vector.  So I wrote my own, based on the squeezed SuiteSparse QR factorization.  (See the next item.)

- A routine to extract the squeezed R factor from the SuiteSparse QR factorization.  I wrote a wrapper for this.

All of this functionality is built into Matlab.

-- Steve

[Tony Kelman]

Those all sound useful and would be worth contributing to the SuiteSparse.jl package if you'd like to do that.

[vav...@uwaterloo.ca]

Tony,

I would be happy to share my prototypes for these wrappers with you or anyone else, but I would not be interested in implementing and testing library versions of them.  This was the point of my first posting.  It's not clear to me who would be interested in that thankless and time-consuming task.  Just like you, I have a "day job"!

-- Steve

[Tony Kelman]

Given that SuiteSparse.jl is a package and not part of base Julia, it doesn't have to meet the quality standards that code contributed to Base usually would. Publicly posted code (with a license or agreement to contribute under an existing license) is better than nothing, and the maintainers of that package can decide what quality bar it would have to meet - in terms of tests and documentation, which is what I assume you mean - before merging it. Even just an open pull request would serve as a starting point for discussion to see who else would be interested.