Provide support for solving sparse linear systems of equations
Issue #58
new
Description
So far the Blaze library is missing the functionality to solve sparse linear systems of equations (i.e. a linear system with a sparse system matrix). Since this feature is important for many use cases, Blaze should provide an API to compute the solution for sparse linear systems of equations. The API should both provide an automatic selection of the most efficient way to solve the problem as well as enable a user to manually select the appropriate algorithm.
Conceptual Example
blaze::CompressedMatrix<double,blaze::columnMajor> A;
blaze::DynamicMatrix<double::blaze::columnMajor> B, X;
blaze::DynamicVector<double,blaze::columnVector> x, b;
// ... Resizing and initialization
// Compute the solution of the LSE with automatic selection of the decomposition algorithm
solve( A * x = b );
// Compute the solution of the LSE by means of an LU decomposition
solve<byLU>( A * x = b );
// Compute the solutions of the LSE by means of an Cholesky (LLH) decomposition
solve<byLLH>( A * X = B );
Tasks
- design an API for computing the solution of sparse linear systems of equations
- provide a full documentation of the feature
- ensure compatibility with all existing matrix classes
- ensure compatibility with all existing matrix expressions
- guarantee maximum performance for all matrix types
- add the necessary number of test cases for the entire functionality
Comments (15)
-
-
I have put my work so far up on Github. You can find it here: https://github.com/tjolsen/BlazeIterative . If you'd prefer that it not be posted publicly until it is part of Blaze, let me know and I'm happy to move it into a development branch of this repository.
-
reporter
Hi Tyler!
Thanks a lot for your effort and the link to your current source code. Please give us some time to get an impression of your work and to reflect on the questions you asked.
Best regards,
Klaus!
-
reporter
Hi Tyler!
Thanks again for your efforts to create BlazeIterative. As promised we have taken a closer look at BlazeIterative and gained an impression.
From a functionality point of view it is a valuable addition to Blaze, but we believe that in the current form it is not yet ready to be integrated into Blaze. Mostly this is because it doesn't "feel" like Blaze since you are using your own coding guidelines. Additionally, most of the documentation is missing.
However, since the content is valuable and since we currently don't have much time to integrate this amount of code since we are busy with wrapping up Blaze 3.2 anyway, we have a suggestion: We would like to make BlazeIterative an official Blaze project. This means that on the Blaze homepage we would add a link and description of your library in a new section called "Blaze Projects".
From our point of view the notion of Blaze projects would be a big win for Blaze. Your project would be the first. We see a couple of advantages for all of us. Most importantly, Blaze would benefit from your experience and work and would get access to a collection of iterative solvers. Second, we could continue the work on Blaze 3.2 and stay on schedule. Third, you could continue the work on BlazeIterative independently and at your own pace. And last but not least, we would not have to nail down all the details right now.
Let us know what you think of the idea.
Best regards,
Klaus!
-
Hi Klaus,
That sounds like a fantastic arrangement, and I agree that it will make life easiest for everyone. I will continue to host the project on my Github page. Going forward, I'd love to have your input/advice/feedback about the progress of the library and any features/algorithms that you'd especially like to see. I aim to make it as easy to use with Blaze as possible, so if you have a style guide that I could follow while developing the library (including API design), it would be very helpful.
Best,
Tyler
-
reporter
Hi Tyler!
BlazeIterative is now officially a Blaze project (see the link on the overview page). Additionally we announced this in the news section.
Unfortunately we don't have a style guide. The best thing you can do is to take a look at Blaze classes (as for instance
DynamicVector). However, we will try to support you as well as possible.Thanks again for your efforts to create and maintain BlazeIterative,
Best regards,
Klaus!
-
Hi,
Is BlazeIterative the way to go, or is there any sparse matrix solver within Blaze? Maybe an LU solver?
Best regards, Santiago
-
reporter
Hi Santiago!
Unfortunately Blaze currently does not support any sparse matrix solver itself. There is also no sparse matrix decomposition available. Therefore BlazeIterative is currently the best approach. If you are interested in native support of sparse linear solvers, you can vote for this issue and watch it in order to keep track of the progress of this feature.
Best regards,
Klaus!
-
I'll chime in here. BlazeIterative had been put on a temporary pause in order for me to finish up my thesis work. However, I just defended yesterday, so I plan to continue development shortly. I haven't settled on an API design (eg, how to handle/represent preconditioners, etc), so if anyone has suggestions I'm happy to hear them.
-
Hi,
Are there any updates on this issue? I am interested in sparse factorizations such as Chollesky, LU and QR.
-
reporter
Hi René!
Unfortunately the work on BlazeIterative seems to have been abandoned. Therefore there is no recent update on this issue. However, I agree that this would be very valuable and important addition to Blaze. We have put it high up in our priority list. Additionally I encourage you to implement some of these algorithms yourself and to provide a pull request.
Best regards,
Klaus!
-
Hi Klaus,
Thanks for the update. I might just do that!
-
Hi Klaus, as a step toward full support for solving sparse linear equations, would it be feasible to add support for the basic infrastructure so that I (and others) can contribute more specialized algorithms? Specifically, I’m considering to add an Incomplete LU decomposition class, which is only missing the capability to solve upper and lower sparse triangular matrix equations in Blaze:
// My algorithm provides these: CompressedMatrix<double> l; CompressedMatrix<double> u; // TODO: sparse upper and lower triangular matrix solves: DynamicVector<double> solution; solve(solution, declupp(u), b); // b is the r.h.s. in Ax = b solve(solution, declunilow(l), solution);
-
reporter
Hi Nils!
That sounds like a reasonable plan and a real step forward. Thanks for your offer. It will take some time, tough, until I can pick up the work on that. More specifically, I’m currently writing a C++ book, which needs to be finished till mid of May. I’ll not be able to do anything before that. Thanks again,
Best regards,
Klaus!
-
Hi Klaus, thanks for your time and all the best for your book project 👍
- Log in to comment
Hi Klaus,
I'm interested in putting together a collection of iterative solvers for Blaze, which would be a useful component of a larger sparse-linear-system framework. There are a couple of items to consider before I get too deep in, however.
1) API design: I favor the tag-dispatch system used by ViennaCL. It eliminates the need to manually specify template parameters to select a solver, and it provides a uniform (and unobtrusive) way to specify convergence criteria, maximum iterations, convergence history logging, etc. I have some work done in this direction already and it feels reasonably natural to use.
2) Preconditioners: It may be necessary to introduce a "Preconditioner" type, which supports a (ideally cheap) "M^{-1}*x" interface. They tend to be similar in API to other types of matrix decompositions, so it might be a good idea to define a MatrixDecomposition type which encompasses both the direct solvers and preconditioners. They should be able to be precomputed, if desired, in order to make re-use possible (this is true of full LU/Cholesky/QR decompositions as well).
3) Interfacing to 3rd-party solvers: Unfortunately, it seems difficult to interface to 3rd party solvers (at least for the sparse matrices) because you use a different storage format than most libraries. The use of ValueIndexPair prevents you from pulling out the pointers to the standard CSR/CSC members and passing them into 3rd party libraries (eg, SuiteSparse). Instead, it forces the user to make a copy into a more standardized format manually, which is obviously undesirable. It should be noted that the inability to extract pointers also prevents anyone from using existing GPU-accelerated libraries (like ViennaCL) with Blaze.
Let me know if you're interested in this work and if you have any thoughts about the above points. I'll keep working on some of the unpreconditioned solvers (which are relatively easy to convert to their preconditioned forms, once design decisions about that have been made). Hopefully this will prove useful at the very least as a starting point for another round of implementation.
Best, Tyler Olsen