sparsemat / sprs

sparse linear algebra library for rust

sprs, sparse matrices for Rust

sprs implements some sparse matrix data structures and linear algebra algorithms in pure Rust.

The API is a work in progress, and feedback on its rough edges is highly appreciated :)

Features

Structures

• CSR/CSC matrix
• triplet matrix
• Sparse vector

Operations

• sparse matrix / sparse vector product
• sparse matrix / sparse matrix product
• sparse matrix / sparse matrix addition, subtraction
• sparse vector / sparse vector addition, subtraction, dot product
• sparse/dense matrix operations

Algorithms

• Outer iterator on compressed sparse matrices
• sparse vector iteration
• sparse vectors joint non zero iterations
• simple sparse Cholesky decomposition (requires opting into an LGPL license)
• sparse triangular solves with dense right-hand side

Examples

Matrix construction

  use sprs::{CsMat, CsMatOwned, CsVec};
  let eye : CsMatOwned<f64> = CsMat::eye(3);
  let a = CsMat::new_csc((3, 3),
                         vec![0, 2, 4, 5],
                         vec![0, 1, 0, 2, 2],
                         vec![1., 2., 3., 4., 5.]);

Matrix vector multiplication

  use sprs::{CsMat, CsVec};
  let eye = CsMat::eye(5);
  let x = CsVec::new(5, vec![0, 2, 4], vec![1., 2., 3.]);
  let y = &eye * &x;
  assert_eq!(x, y);

Matrix matrix multiplication, addition

  use sprs::{CsMat, CsVec};
  let eye = CsMat::eye(3);
  let a = CsMat::new_csc((3, 3),
                         vec![0, 2, 4, 5],
                         vec![0, 1, 0, 2, 2],
                         vec![1., 2., 3., 4., 5.]);
  let b = &eye * &a;
  assert_eq!(a, b.to_csr());

For a more complete example, be sure to check out the heat diffusion example.

Documentation

Documentation is available at docs.rs.

Changelog

See the changelog.

Minimum Supported Rust Version

The minimum supported Rust version currently is 1.64. Prior to a 1.0 version, bumping the MSRV will not be considered a breaking change, but breakage will be avoided on a best effort basis.

License

Licensed under either of

• Apache License, Version 2.0, (LICENSE-APACHE or https://www.apache.org/licenses/LICENSE-2.0)
• MIT license (LICENSE-MIT or https://opensource.org/licenses/MIT)

at your option.

Contribution

Unless you explicitly state otherwise, any contribution intentionally submitted for inclusion in the work by you, as defined in the Apache-2.0 license, shall be dual licensed as above, without any additional terms or conditions.

Please see the contribution guidelines for additional information about contributing.

sprs's Issues

Shrink matrices produced in operations

Sparse matrix products may allocate more non zero storage than needed. Shrinking the storage
could be a good idea since the resulting matrix is owned.

Indexing using NnzIndex

Indexing a sparse matrix with an NnzIndex can be done in constant time, this can be a nice option for someone looking into a performant way to modify a particular value multiple times.

Expose a trait to generalize over CsMat variants

This should enable easier and more general APIs.

For instance the vstack, bmat constructors could benefit from this API with an API like

pub fn vstack<N>(mats: AsRef<[CsMatTrait<N>]>)

Documentation for algorithms in the `linalg` module

It would be helpful to have (even rudimentary) documentation for these functions, instead of having to figure out the correct usage from the test cases.

Triplet format constructor

Enable constructing a sparse matrix from the (i, j, x) triplet format.

changelog

It would be great to add a changelog file

Per element application of functions

@dylanede's description of the feature:

One set of very useful things that is missing is support for per-element operations between matrices/vectors. For example if I want to multiply one vector by a constant and add it to another. I think the most practical way for the time being of supporting per-element operations would be to add some functions with differing numbers of arguments that take some borrowed matrices/vectors, and a closure to operate on them to produce the result, or maybe write into an existing matrix/vector. For example, something like

fn per_element1<N, F: FnMut(N) -> N>(a: CsVecView<N>, f: F) -> CsVecOwned<N>;
fn per_element2<N, F: FnMut(N, N) -> N>(a: CsVecView<N>, b: CsVecView<N>, f: F) -> CsVecOwned<N>;
fn per_element3<N, F: FnMut(N, N, N) -> N>(a: CsVecView<N>, b: CsVecView<N>, c: CsVecView<N>, f: F) -> CsVecOwned<N>;

Block constructor

In the spirit of scipy's bmat constructor.

Usage of zip in products may be inefficient

As explained here, using zip to iterate on members with the same length is inefficient bacause the compiler won't be able to remove both bounds checkings.

Some tricks exist for slices, but here the situation is somewhat different, maybe an unsafe block (or a dedicated primitive implemented as unsafe) is necessary.

Integration tests for sparse/dense products

Because it seems doc-tests cannot use other crates.

Hosted documentation needs to be updated

I've noticed that there are features available in the crates.io published version that are not shown in the online documentation, such as LdlNumeric and LdlSymbolic, so really a cargo doc should be redone.

Bind to existing solvers

Offer optional bindings to SuiteSparse's solvers.

Type checked invariants

Use phantom data to ensure at compile time that algorithms expecting e.g. sorted indices get matrices satisfying that property .

Isolate cholesky factorization behind cargo feature

As an adapted code it is probably under LGPL.

Mutable view type

A CsMat variant with mutable data vector, where it is possible to override an element's value with an indexing operation (which should probably panic if the location was not a nnz).

Depend on ndarray for dense matrices

Ergonomics could be improved

After trying out the library for a bit, I've found some places in the API that are a bit painful to work with at the moment. Number one is the (lack of) interop with dense matrices and vectors. It would be very nice to be able to construct a sparse array from a Vec<N>, for example, or dot product a sparse vector with a Vec<N>. I won't mention too much more about dense storage interop since you've got some issues opened about this.

One set of very useful things that is missing is support for per-element operations between matrices/vectors. For example if I want to multiply one vector by a constant and add it to another. I think the most practical way for the time being of supporting per-element operations would be to add some functions with differing numbers of arguments that take some borrowed matrices/vectors, and a closure to operate on them to produce the result, or maybe write into an existing matrix/vector. For example, something like

fn per_element1<N, F: FnMut(N) -> N>(a: CsVecView<N>, f: F) -> CsVecOwned<N>;
fn per_element2<N, F: FnMut(N, N) -> N>(a: CsVecView<N>, b: CsVecView<N>, f: F) -> CsVecOwned<N>;
fn per_element3<N, F: FnMut(N, N, N) -> N>(a: CsVecView<N>, b: CsVecView<N>, c: CsVecView<N>, f: F) -> CsVecOwned<N>;

(If this is to support dense vector/matrix interop, it may be worth considering a trait for vector views that dense vectors/matrices can implement as well as sparse ones)
(Edit: There may need to be two sets of functions - ones for operation only on non-zero elements, and ones that operate on all elements, so that it is possible to have sparse usage of the function f)

Another more basic set of things that is missing is more operations on elements of matrices/vectors. E.g. retrieving and setting a value at a particular position just like a normal matrix. The implementation should handle the case of whether a value is zero or not and act appropriately, including insertion of indices into the storage vectors. To better support this kind of usage you may want to add support for mutable views, and think more carefully about the internal representation.

Speaking of views, it would be very nice to have views that are arbitrary sub-regions of the matrix/vector, as well as regions that extend past the boundaries of the original matrix/vector, returning zeroes for out-of-bounds locations.

Another thing that would need some careful thought would be support for reusing allocated storage and preallocating. This would would improve performance for cases where a particular routine must operate several times, and could eliminate repeated heap (de)allocations.

I'll add more thoughts to this issue if I think of anything else.

Unsafe unchecked constructors

CsMat and CsVec constructors perform structure check by default. However, in some cases, we already know the structure is correct (when taking a view into an existing CsMat/CsVec for instance).

To avoid performing checks twice, and unsafe unchecked constructor could be used.

This should be unsafe as algorithms are free to rely on some properties ensured by structure check to perform unchecked array access. For instance, structure checks guarantee that indices are coherent with the matrix dimensions.

Triangular view

Be able to create a new indptr to have a triangular view into an existing matrix.

Expose arithmetic operations using traits

The traits implementations should pick the appropriate method given the arguments' structures.

Dense traits

Add a DenseMatrix and a DenseVec trait, to enable interopeability with other libraries.

Cleaner generlization of nnz_or_zip

Currently nnz_or_zip is both defined as an extension trait and as a method for sparse vector iterator,
ideally it should be the former.

LU factorization requires unsorted matrices

This might require implementing #16.

Add a changelog

Depend on rendarray for dense mats

Now that view types have landed (rust-ndarray/ndarray#8), the rendarray crate should be able to replace the dense_mats crate, while being more feature complete and actively developped.

OuterIterator should be able to operate by reference

This is needed to remove unnecessary Copy bounds (#31). Also this would be a good step towards mutable iteration.

vstack and bmat macros

Also expose vstack and other similar functions via a macro that would call .borrow() on the arguments. This should enable vstacking with different storages.

Store indptr as isize

As sparse matrices contain two arrays of equal shape, it is impossible to have more non-zeros than the maximum value of an isize.

Also, some factorization algorithms require signed integers to be able to mark visited nodes.

Copy bound is present even when unneeded.

This adds some burden for implmenting new functions.

csr-csr product structure is data dependent

It would be advisable to have the resulting structure depend only on the inputs structure, and not on the numerical zeros. This might mean going back to a workspace using Option.

Handle permutations in Cholesky factorization

Probably hard as new iterators will have to be derived.

ArrayView

Have an array view data structure, and an associated trait enabling to take an array view into a matrix.
ArrayView implementation of Mul will perform a scalar multiplication instead of a matrix product.

Check that dimensions agree in triangular solvers

type alias implementations do not show in doc

Looks related to rust-lang/rust#19381

Can be fixed by explicitly specifying the type in the impl block, as a workaround.

More generic CsMat storage

Currently, CsMat requires that its indices and inptr members have the same storage. However this is inconvenient when creating a matrix view over a CsVec (extra useless fields have to be added to CsVec). Even though it would add some burden for CsMat implementations, it would probably be better (and more generic) to have CsMat more generic. Type aliases should ensure a pleasant path for common cases (ie full view or owned).

Make CsMat::at_inner_outer use CsVec::at

For code reuse.

SpVecView

There's SpMatView for abstracting over matrix implementations, but no SpVecView for abstracting over vector implementations. When writing a function that takes vectors and/or matrices, the signature gets very long without using traits like this, since extra type parameters must be introduced specifying the storage details for each vector/matrix argument. SpMatView goes part of the way towards dealing with this, but SpVecView is really required as well.

Consistent naming

For CsMat and CsVec using Vec storage, methods exposing Vec storage should take similar names.
CsMat and CsVec methods should have similar names when possible.

Friendly interface to sparse/dense dot product

Dot products should work with anything that can be interpreted as a dense vector, so the method should probably take an IntoIterator bound for its argument.

vstack and hstack builders

Be able to construct matrices by stacking them. Easy if the stacking direction agrees with the structure, harder otherwise.

CsMat type name

Can you explain the naming? Why not just Mat type? I see that you use two formats CSR and CSC, so Cs prefix in type name doesn't clarify anything.

Am I missing something? If not, I would consider renaming to Mat. It can still be something like cs::Mat, but in that case users could actually use cs::Mat. It gives more choice to the user.

Arbitrary subivews

Be able to take a view that's an arbitrary block inside a matrix.

Fast paths should be provided for easy paths (such as taking a row subview of a CSR matrix).
It might be nice to be able to use ndarray's indexing syntax for the general path.

Rename at_ familiy of method into get_ for consistency with std

These methods have semantics closely related to the slice's get methods.

Insert element into CsMatOwned

Even though this is inefficient, it can sometimes be useful.

Automate doc generation

It's apparently possible to automate documentation generation for th master branch.

It would also be nice to build the documentation for each version.

Storage should be implemented as an unsafe trait

See rust-lang/rust#29701
The implementation in bluss/rust-ndarray looks like a good example of how to do it.

Relicense under dual MIT/Apache-2.0

This issue was automatically generated. Feel free to close without ceremony if
you do not agree with re-licensing or if it is not possible for other reasons.
Respond to @cmr with any questions or concerns, or pop over to
#rust-offtopic on IRC to discuss.

You're receiving this because someone (perhaps the project maintainer)
published a crates.io package with the license as "MIT" xor "Apache-2.0" and
the repository field pointing here.

TL;DR the Rust ecosystem is largely Apache-2.0. Being available under that
license is good for interoperation. The MIT license as an add-on can be nice
for GPLv2 projects to use your code.

Why?

The MIT license requires reproducing countless copies of the same copyright
header with different names in the copyright field, for every MIT library in
use. The Apache license does not have this drawback. However, this is not the
primary motivation for me creating these issues. The Apache license also has
protections from patent trolls and an explicit contribution licensing clause.
However, the Apache license is incompatible with GPLv2. This is why Rust is
dual-licensed as MIT/Apache (the "primary" license being Apache, MIT only for
GPLv2 compat), and doing so would be wise for this project. This also makes
this crate suitable for inclusion and unrestricted sharing in the Rust
standard distribution and other projects using dual MIT/Apache, such as my
personal ulterior motive, the Robigalia project.

Some ask, "Does this really apply to binary redistributions? Does MIT really
require reproducing the whole thing?" I'm not a lawyer, and I can't give legal
advice, but some Google Android apps include open source attributions using
this interpretation. Others also agree with
it.
But, again, the copyright notice redistribution is not the primary motivation
for the dual-licensing. It's stronger protections to licensees and better
interoperation with the wider Rust ecosystem.

How?

To do this, get explicit approval from each contributor of copyrightable work
(as not all contributions qualify for copyright, due to not being a "creative
work", e.g. a typo fix) and then add the following to your README:

## License

Licensed under either of

 * Apache License, Version 2.0 ([LICENSE-APACHE](LICENSE-APACHE) or http://www.apache.org/licenses/LICENSE-2.0)
 * MIT license ([LICENSE-MIT](LICENSE-MIT) or http://opensource.org/licenses/MIT)

at your option.

### Contribution

Unless you explicitly state otherwise, any contribution intentionally submitted
for inclusion in the work by you, as defined in the Apache-2.0 license, shall be dual licensed as above, without any
additional terms or conditions.

and in your license headers, if you have them, use the following boilerplate
(based on that used in Rust):

// Copyright 2016 sprs developers
//
// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
// http://www.apache.org/licenses/LICENSE-2.0> or the MIT license
// <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your
// option. This file may not be copied, modified, or distributed
// except according to those terms.

It's commonly asked whether license headers are required. I'm not comfortable
making an official recommendation either way, but the Apache license
recommends it in their appendix on how to use the license.

Be sure to add the relevant LICENSE-{MIT,APACHE} files. You can copy these
from the Rust repo for a plain-text
version.

And don't forget to update the license metadata in your Cargo.toml to:

license = "MIT/Apache-2.0"

I'll be going through projects which agree to be relicensed and have approval
by the necessary contributors and doing this changes, so feel free to leave
the heavy lifting to me!

Contributor checkoff

To agree to relicensing, comment with :

I license past and future contributions under the dual MIT/Apache-2.0 license, allowing licensees to chose either at their option.

Or, if you're a contributor, you can check the box in this repo next to your
name. My scripts will pick this exact phrase up and check your checkbox, but
I'll come through and manually review this issue later as well.

• @vbarrielle

sparse triangular solve

Enhance documentation

Most methods lack documentation.

Dot product on ArrayView is impractical

Requiring to take a reference on a view is unnatural, it would be nice to have the same system that enables taking CsVecView by value. However this might require some re-borrowing API in ndarray, or better, waiting for Rust features making that case easier.