Sparse structures is an exciting studying field that enables you to build, store, modify and retrieve scattered data in a fast and efficient way. Along with algebraic operations, its usage ranges from several Artificial Intelligence areas, Operations Research, simulations and many others places where, e.g., a dense (non-sparse) matrix-vector multiplication (MxV) is too costly or even pointless for very sparse problems.

In this post, I am going to talk about ndsparse, a batteries-included library written in Rust that is intended to provide different types of multidimensional sparse structures where you can choose the format that best suits you.

For the hopeful

Before anything else, ndsparse isn't a multidimensional sparse algebra/arithmetic library (disappointment face) because of its self-contained responsibility and complexity. Furthermore, a really good implementation of such a library would require a titanic amount of research, work and free time that I don't have.

This intended limitation restricts the usage of this project by a lot but it is still useful for store, transform and retrieve use-cases. One can also be a hero and use some of the supported structures as a building foundation for higher level libraries.

Checkout sprs for an awesome "rustic" sparse linear algebra library.

Supported formats

There are a bunch of different 2D sparse structures that determine the space-usage and the asymptotic limit of a given operation, some are generic and others are more problem specific. E.g.: BSR, COO, CSC, CSR, DIA, DOK, ELL, JDS and LIL.

Two (or three, depending on your POV) widely known formats were picked for adaptation, namely, COO and CSC/CSR.

COO (Coordinate)

Probably the most intuitive format, fits gracefully for N-dimensions. Just need a set of indices and its corresponding non-zero elements.

use ndsparse::coo::CooArray;
// As odd as it may seem, this illustration is just a guide to get a grasp of
// a 5D structure.
//
// The order is up to the caller. In this case, the dimensions [a, b, c, d, e] were
// arranged as follows:
//
// a: top to bottom
// b: left to right
// c: front to back
// d: top to bottom
// e: left to right
//
//          ___ ___ ___            ___ ___ ___            ___ ___ ___
//        /   /   /   /\         / 3 /   /   /\         /   /   /   /\
//       /___/___/___/ /\       /_3_/___/___/ /\       /___/___/___/ /\
//      /   /   /   /\/ /\     /   /   /   /\/ /\     /   / 4 /   /\/ /\
//     /___/___/___/ /\/ /    /___/___/___/ /\/ /    /___/_4_/___/ /\/ /
//    /   /   /   /\/ /\/    /   /   /   /\/ /\/    /   /   /   /\/ /\/
//   /___/___/___/ /\/ /    /___/___/___/ /\/ /    /___/___/___/ /\/ /
//  /   /   /   /\/1/\/    /   /   /   /\/ /\/    /   /   /   /\/ /\/
// /___/___/___/ /\/ /    /___/___/___/ /\/ /    /___/___/___/ /\/ /
// \___\___\___\/ /\/     \___\___\___\/ /\/     \___\___\___\/ /\/
//  \___\___\___\/ /       \___\_2_\___\/ /       \___\___\___\/ /
//   \___\___\___\/         \___\___\___\/         \___\___\___\/
//
//          ___ ___ ___            ___ ___ ___            ___ ___ ___
//        /   /   /   /\         /   /   /   /\         /   /   / 6 /\
//       /___/___/___/ /\       /___/___/___/ /\       /___/___/_6_/6/\
//      /   /   /   /\/ /\     /   /   /   /\/ /\     /   /   /   /\/ /\
//     /___/___/___/ /\/ /    /___/___/___/ /\/ /    /___/___/___/ /\/7/
//    /   /   /   /\/ /\/    /   /   /   /\/ /\/    /   /   /   /\/ /\/
//   /___/___/___/ /\/ /    /___/___/___/ /\/ /    /___/___/___/ /\/ /
//  /   /   /   /\/ /\/    /   /   /   /\/ /\/    /   /   /   /\/ /\/
// /___/___/___/ /\/ /    /___/___/___/ /\/ /    /___/___/___/ /\/ /
// \___\___\___\/ /\/     \___\___\___\/ /\/     \___\___\___\/ /\/
//  \___\___\___\/ /       \___\___\___\/ /       \___\___\___\/ /
//   \___\___\___\/         \___\_5_\___\/         \___\___\___\/
let _coo_array_5 = CooArray::new(
    [2, 3, 4, 3, 3],
    [
        ([0, 0, 1, 1, 2].into(), 1),
        ([0, 1, 0, 1, 1].into(), 2),
        ([0, 1, 3, 0, 0].into(), 3),
        ([0, 2, 2, 0, 1].into(), 4),
        ([1, 1, 0, 2, 1].into(), 5),
        ([1, 2, 3, 0, 2].into(), 6),
        ([1, 2, 3, 2, 2].into(), 7),
    ],
);

You might be wondering what these Array and into() things are. Well, we will get there in a minute.

CSL (Compressed Sparse Line)

A generalization of the Compressed Sparse Column (CSC) and Compressed Sparse Row (CSR) formats for N-dimensions. Since all data is compressed line-by-line, this nomenclature came naturally.

Basically, three indexed storage are needed, one for the data itself, one for the line index of each data and one to indicate the number of non-zero elements of each line. Here, each line can also be interpreted as the innermost dimension or the rightmost dimension.

use ndsparse::csl::CslArray;
// Two cuboids illustrating a [2, 3, 4, 5] 4D in a [w, y, z, x] order, i.e., each "line"
// or 1D representation is a left to right row and each "matrix" or 2D representation
// is filled in a top-down manner.
//
//  w: left to right
//  y: top to bottom
//  z: front to back
//  x: left to right
// 
//          ___ ___ ___ ___ ___            ___ ___ ___ ___ ___
//        /   /   /   / 4 / 5 /\         /   /   /   /   /   /\
//       /___/___/___/_4_/_5_/5/\       /___/___/___/___/___/ /\
//      /   /   /   /   /   /\/ /\     /   /   / 9 /   /   /\/ /\
//     /___/___/___/___/___/ /\/ /    /___/___/_9_/___/___/ /\/ /
//    /   / 3 /   /   /   /\/ /\/    /   /   /   /   /   /\/ /\/
//   /___/_3_/___/___/___/ /\/ /    /___/_ _/___/___/___/ /\/ /
//  / 1 /   /   / 2 /   /\/ /\/    /   /   /   /   /   /\/ /\/
// /_1_/___/___/_2_/___/ /\/8/    /___/___/___/___/___/ /\/ /
// \_1_\___\___\_2_\___\/ /\/     \___\___\___\___\___\/ /\/
//  \___\___\_6_\___\___\/ /       \___\___\___\___\___\/ /
//   \___\___\_7_\___\___\/         \___\___\___\___\___\/
let _csl_array_4_ = CslArray::new(
  [2, 3, 4, 5],
  [1, 2, 3, 4, 5, 6, 7, 8, 9],
  [0, 3, 1, 3, 4, 2, 2, 4, 2],
  [0, 2, 3, 3, 5, 6, 6, 6, 6, 7, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9],
);

Yes, each additional dimension significantly raises the number of line offsets. The good thing is the possibility of slicing sub-dimensions views and the bad thing is obviously the increased space usage.

If you are still confused or don't want to manually create instances, it is easier to construct a valid CSL by using CslLineConstructor.

Rust-ish features

Have you ever heard of Rust? If not, you should probably start getting used to it. Currently, several companies are using it [1] and several programs are being written or rewritten in Rust like Firefox, ripgrep and librsvg. The reasons for such success are numerous: It is very fast, includes high-level facilities, the ownership rules prevent many memory management pitfalls (as well as thread safety), has an incredible community and many other cool things.

Unfortunately, for truly N-dimensional structures, the nightly constant generics feature is a hard requirement and even with it, there are no std implementations for arrays greater than 32 elements [2] [3], which leads to the creation of the ArrayWrapper alternative that is used heavily internally, thus, the Into::into() method conversion from [T; N] to ArrayWrapper<T, N>.

Putting all that aside, ndsparse has a lot of optional features. Initially, considering a #[no_std] and "no alloc" environment, none of them are used by default, giving the user the freedom to choose whatever is required.

[dependencies]
ndsparse = { features = ["alloc", "with_arrayvec", "with_rand", "with_rayon", "with_serde", "with_smallvec", "with_staticvec"], version = "0.2" }

Different backends for storage

Owned structures can use static arrays, heap-allocated vectors or even third-party dependencies like ArrayVec to store the underlying data.

use ndsparse::{
  coo::{CooRef, CooStaticVec},
  csl::{CslRef, CslSmallVec, CslStaticVec}
};

// CSL

let mut csl_small_vec = CslSmallVec::<i32, 2, 5, 26>::default();
let mut csl_static_vec = CslStaticVec::<i32, 2, 5, 26>::default();
csl_small_vec.constructor().next_outermost_dim(5).push_line(&[1, 2], &[0, 3]);
csl_static_vec.constructor().next_outermost_dim(5).push_line(&[1, 2], &[0, 3]);
assert!(
  csl_small_vec.line([0, 0]) == Some(CslRef::new([5], &[1, 2][..], &[0, 3][..], &[0, 2][..]))
);
assert!(csl_small_vec.line([0, 0]) == csl_static_vec.line([0, 0]));

// COO

let coo_ref_data = [([0, 0, 0].into(), 1)];
let coo_ref = CooRef::new([9, 9, 9], &coo_ref_data[..]);
let coo_static_vec = CooStaticVec::new([9, 9, 9], [([0, 0, 0].into(), 1)]);
assert!(coo_ref.value([0, 0, 0]) == Some(&1));
assert!(coo_ref.value([0, 0, 0]) == coo_static_vec.value([0, 0, 0]));

Iterators and parallel iterators

Outermost iterators (first or left most dimension) for CSL can be retrieved in parallel using rayon. Pretty useful for very large structures.

use rayon::prelude::*;
let are_equal = some_csl
  .outermost_rayon_iter()
  .enumerate()
  .all(|(idx, csl_ref)| csl_ref == some_csl.outermost_iter().nth(idx).unwrap());
assert!(are_equal, true);

COO is more straightforward. For example, one can use some_coo.data().par_iter().for_each(|_| {}).

Future

Slicing isn't great for CSL, COO should receive more love, there isn't any agnostic transformation like resize or transpose and more formats could be added. All of these TODO's might be added at some point in the future with enough free time and willingness.

Last but not the least, I think that I invented a new sparse structure that is space-efficient and enables a fine-grained control over sparsity. More on this in a later post.