Implement Sparse Left-Looking LU factorization #1289
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
This allows a sparse matrix to be used for efficient solving with a dense LU decomposition.
It generally follows lectures from Tim Davis
https://www.youtube.com/watch?v=wiL_jIuENkk&ab_channel=TimDavis
(Lectures 1-3), and implements a basic left looking sparse LU factorization.
I primarily focus on CSC matrices, since they are what is focused on in the lecture, identical to what is done in Matlab.
I believe for most use cases this should be efficient and correct, and thus it serves as a useful to building even more customized solvers down the line.
Summary of changes:
Sparsity Pattern Builder: Allowing for building up sparsity patterns, rather than constructing them all in one go. This seemed to be a drawback of the previous design, and I've added a builder to construct Sparsity Patterns. This was necessary when computing the sparsity of the decomposed LU matrix.
Sparse lower triangular solving. For most users, the result will probably want
LUx = b
, withb
a dense vector. While constructing the LU matrix though, it calls a sparseLy = s
, withs
sparse. This allows for building up a sparse LU pattern iteratively. The LU matrix is packed into a single sparse matrix, with the diagonal belonging to U, andL
has implicit 1s along the diagonal.Dense lower/upper triangular solving. This is what the end user will interact with. I've added it to take slices, but I should probably convert it to DMatrices. Slices feel more natural to me though, since it decouples the requirement of the input being a matrix, as it's just solving a vector anyway.
r? @sebcrozet @Andlon
Not sure if y'all have time to review.
Motivation:
Currently I'm trying to reimplement a paper that uses libigl within Rust, but there's a significant amount of library functionality that is missing. Specifically, quadratic programming doesn't exist, so this is a step to making it available.
Closes #1197