ENH: New sparse format, Multiply-Compressed Sparse Row (MCSR), for matrices with repeated nonzero values

[cmclean]

Hi scipy developers,

I’d like to propose a new sparse matrix format for SciPy: Multiply-Compressed Sparse Row (MCSR). I have a working implementation (Python + C++ sparsetools) and am looking for feedback before opening a PR.

The motivation for this format is based on the observation that many real-world sparse matrices have rows in which the same nonzero value appears many times. For example:

- Single-cell RNA sequencing (scRNA-seq) count matrices. Lowly-expressed genes produce many identical small integer counts in each of many thousands of cells.

- Adjacency matrices of weighted graphs with discrete edge weights.

- Potentially in quantized or rounded scientific data.

CSR stores every nonzero value individually. If a row has 5,000 nonzero entries but only 10 distinct values, CSR wastes space storing the same number thousands of times. MCSR compresses this redundancy via run-length encoding (RLE) on the values within each row.

Format description

MCSR is a lossless compression of CSR that performs run-length encoding of the nonzero data values. For an M × N matrix, it stores five arrays:

Array	Type	Length	Description
indices	int	nnz	Column indices (sorted by value within each row)
indptr	int	M+1	Row pointer into indices (identical semantics to CSR)
data	Any	nrle	Unique nonzero values per row, sorted per row and concatenated
counts	uint	nrle	Run length for each entry in data. sum(counts) == nnz.
counts_indptr	int	M+1	Row pointer into data and counts.

Where nrle is the total number of (value, count) pairs across all rows.

Code to reconstruct row i is then:

rle_start, rle_end = counts_indptr[i], counts_indptr[i + 1]
row_data = np.repeat(data[rle_start:rle_end], counts[rle_start:rle_end])
row_indices = indices[indptr[i]:indptr[i + 1]]

Both CSR and MCSR store arrays of size nnz and M+1, but MCSR compresses the data values from nnz to 2 * nrle. In practice, on a variety of single cell datasets, this yields in-memory compression ratios MCSR:CSR of 0.51:1.0.

Because the total number of data values stored can be much smaller than CSR, performance of arithmetic operations can also be much more efficient. Evaluation of iter, row/column/global sums, log1p and expm1, and extraction of the top 1k highest values per row across a representative set of 457 single cell datasets show average 3.1- to 17.2-fold performance improvement.

A related work ( Value-Compressed Sparse Column (VCSC): Sparse Matrix Storage for Single-cell Omics Data ) that independently developed essentially the same implementation, but column-oriented, performed additional benchmarking on random datasets and show similar gains and the degradation of performance as data redundancy decreases.

I have an initial implementation at GitHub - cmclean/scipy at multiply-compressed-sparse-row · GitHub . It consists of the Python classes mcsr_array and mcsr_matrix that follow the sparray / spmatrix pattern, construction from dense and sparse matrices, full IndexMixin support, a variety of arithmetic and aggregation functions, and matrix multiplication routines, all optimized using C++ sparsetools kernels. The main code is reasonably solid at this point; while the test coverage is high now that I better understand the testing setup for sparse formats it would need an overhaul to match the current conventions of the library.

So I guess my main question is whether there is interest in adding this sparse format to SciPy? It would be directly beneficial for single cell omics data, particularly if integrated for support in AnnData which I have looked into and would be straightforward. I understand the maintenance burden and am happy to discuss scope. Also happy to discuss other things like the naming, whether a column-oriented analog should exist, or additional performance-critical use cases I should benchmark against. Thanks for considering!

[dschult]

Thanks @cmclean for this proposal and the work you’ve put into preparing a PR to demonstrate it.

There hasn’t been a new format for quite a while, but it might be time to consider that. We are just about to deprecate sparse matrices and we’ve put a bit of effort into n-D sparse arrays. From my perspective you have a motivating use-case. Maintenance is always a concern. How likely are you to be around for the next few years to see this through release and a couple release cycles to get rid of most bugs? A related question, though it might not seem like it: how well-known is this sparse format. You link to a similar format for columns. Are there links to articles about this format? Did you create it on your own? Are there publications or implementations elsewhere?

I don’t think you would need to put a lot of work into supporting this for spmatrix. That is scheduled for deprecation next release and removal 2(+?) releases later. But as you have noticed, it isn’t all that hard to support both. Most of the guts are the same.

My searches show that MCSR is already used in the sparse world for Modified Compressed Spasre Row. And there are many ways to do multiple compressed formats. I know of one where the indptr itself gets compressed a second time (useful for matrices with many empty rows). Your proposed format leaves the indices array unchanged. It really focuses on compressing the data array. So maybe Compressed Row Data with an abbreviation “crd" would be less confusing (could be cdr and cdc but from a typo perspective that seems very close to csr and csc). Anyway, I haven’t searched on that name or abbreviation or anything. Its’ probably something you should think about. All the current formats have a 3-char format code. I’m not sure how important that is, but it’d be nice style in some sense.

Other considerations always seems to include the extra build time and .so filesize. But really the main focus is about whether this will help users in a way that supports the extra maintenance and code complexity. I didn’t look for what is missing from the branch you’ve put together. How much of the sparse function/method space does this implement so far? Does it even make sense to implement everything, or would it make more sense to have it specialize the core functions and then convert to other formats when other operations are needed.

I suspect it will help with review and understanding the format if you can write snippets which convert this to csr and/or coo – probably csr. You’ve got a great 3-liner to show getting one row of data and indices from data, indices, indptr, counts, counts_indptr. I notice that that 3-liner does not sort the resulting indices. We could use A.sort_indices Also, how about the other direction? Is it easy to build the new format from “csr” format using numpy tools? (if not, that’s ok. but it helps understand the format better)

[dschult]

Another question I have: this sparse format would only be effective for discrete valued dtypes. Are you looking to limit the dtype of this format? How do you handle floating point values?

[cmclean]

Hi Dan,

Thank you for taking the time to write these thoughtful responses, I appreciate it. To answer your questions:

I would be happy to do maintenance to support the format.

Re how well known it is, that arXiv post I linked to is the only other description I am aware of. It appears to have an independent implementation (GitHub - Seth-Wolfgang/PyVSparse · GitHub) that is a thin wrapper over their C++ implementation. I was not aware of this when writing the version at scipy/scipy/sparse/_mcsr.py at multiply-compressed-sparse-row · cmclean/scipy · GitHub but the latter is built specifically targeting the scipy.sparse codebase so its conventions are much more similar/easier to integrate.

Acknowledged re spmatrix / sparray support, it seems mostly straightforward to support and then handle deprecation.

Re naming, acknowledged that we could find a better name, I didn’t spend a huge amount of time thinking about it yet… "crd" is kind of nice as it could mean “compressed run-length-encoded data”, or "rer" for “run-length encoded row”.

Regarding coverage so far, it natively supports arbitrary slicing, all sum aggregations, and a few arithmetic functions of interest (log1p and expm1, though adding these are trivial), and conversion to/from csr. Other operations like reshape are supported via a conversion to CSR first, using its implementation, and then converting back. No mutating ops are currently supported but could easily have the same (albeit inefficient) path to support. I agree that specializing these core functions that are certainly useful, and delegating to CSR implementations for others until there is a demonstrated need for native support, is a good path forward.

Conversion to/from CSR is already supported. To CSR is pure numpy (scipy/scipy/sparse/_mcsr.py at multiply-compressed-sparse-row · cmclean/scipy · GitHub), reproduced below for convenience. As you noted on the row conversion, this does rely on CSR’s already implemented sort_indices to bring it to canonical format:

  def tocsr(self, copy: bool = True) -> csr_matrix:
    """Convert to a CSR sparse matrix."""
    del copy  # Unused -- it is always true.
    csr_indptr = np.array(self.indptr)
    csr_indices = np.array(self.indices)
    # Expand the RLE data into a CSR data array.
    csr_data = np.repeat(self.data, self.counts)

    csr_cls = csr_array if isinstance(self, sparray) else csr_matrix
    csr = csr_cls(
        (csr_data, csr_indices, csr_indptr), shape=self._shape, copy=False
    )
    # This is a valid CSR matrix, but the data is sorted by value. Use the
    # CSR optimized sort_indices() to restore column-sorted order.
    csr.sort_indices()
    return csr

The constructor of an mscr_{matrix,array} can take in a CSR matrix as argument and delegates to the _from_csr() function (scipy/scipy/sparse/_mcsr.py at multiply-compressed-sparse-row · cmclean/scipy · GitHub). As currently written, this iterates through each row of the CSR matrix and orders and compresses each row via the following numpy operations:

def _csr_row_to_mcsr(
    indices: np.ndarray,
    data: np.ndarray,
) -> tuple[np.ndarray, int, np.ndarray, int, np.ndarray]:
  """Convert a CSR row to MCSR format.

  A single CSR row is represented by two 1-D arrays of identical length,
  corresponding to the non-zero column indices and their associated data values.

  The analogous MCSR representation uses three 1-D arrays:
    - The non-zero column indices, which must be sorted by data value.
    - The run lengths of each data value.
    - The unique data values, of the same length as the run lengths.

  This function returns these three arrays, along with the two unique lengths.

  Args: indices : ndarray Column indices of non-zero values. data : ndarray Data
  values corresponding to the column indices.

  Returns:
    sorted_indices : ndarray
        Column indices of nonzero values, sorted by value.
    len(sorted_indices) : int
        Number of nonzero values (nnz).
    counts : ndarray
        Run lengths corresponding to each data entry.
    len(counts) : int
        Number of RLE entries (unique nonzero value runs).
    unique_vals : ndarray
        Sorted unique nonzero values.
  """
  sort_order = np.argsort(data)
  sorted_indices = indices[sort_order]
  sorted_data = data[sort_order]
  # Compress the sorted data using RLE.
  unique_vals, counts = np.unique(sorted_data, return_counts=True)
  return sorted_indices, len(sorted_indices), counts, len(counts), unique_vals

The rest of the _from_csr function is just stitching each of these row representations into arrays.

Finally, regarding the question about float support, the PR does support it now. You’re right that for data that is generated natively in float this may be a bad format choice – in worst case this will increase storage usage by nnz since in addition to storing data we store the array of associated counts (all with entry 1). In practice with single-cell data we do still see performance improvement though because the floats arise from a row-constant scaling factor applied to integer values, so the float representations are truly bit identical.

[cmclean]

Any other thoughts on this proposed enhancement? Or how do we move from discussion to a decision on whether I should push forward with cleaning up the implementation for a PR or drop it?

thanks,
Cory

[dschult]

The feeling I’ve gleaned from short discussions and history of formats is that this is probably too narrow an application to be added to SciPy proper. I’d be interested to see how hard it would be to add a “plug-in” scipy sparse format in a separate package. Basically I’m thinking of a package that include the python and C++ code from your branch – but compiles and packages it separately from scipy in a way that makes it directly compatible with the scipy API. I’m not very experienced with packaging so this is probably a fantasy that could get bogged down in a dependency swamp.

So, do I have your permission to use your commit to play with that approach? Would you want to do that? I think there are other formats that people would be willing to create that SciPy would not choose to support. I just don’t know how hard that would be to get working well.

[lucascolley]

@dschult

I’d be interested to see how hard it would be to add a “plug-in” scipy sparse format in a separate package. Basically I’m thinking of a package that include the python and C++ code from your branch – but compiles and packages it separately from scipy in a way that makes it directly compatible with the scipy API. I’m not very experienced with packaging so this is probably a fantasy that could get bogged down in a dependency swamp.

The packaging problem is not infeasible (and I’d be happy to help), but it matters how exactly you envision it being ‘directly compatible’. Working as input/output with other SciPy functions that accept existing sparse formats? Constructible via asformat — SciPy v1.18.0 Manual? Present in the scipy.sparse namespace as a class?

[dschult]

I think it is definitely worth thinking about what would be needed to allow outside formats to work with other scipy functions that accept sparse formats.

For “asformat”, it is straightforward to create tocsr methods on the new format class, but not so easy to go the other direction. I suppose monkeypatching in some way might work – and maybe there would be a way to register a plugin format that asformat could use to try something like a new_format.fromcsr. The outside package would need to define to/from methods.

And I was not expecting at all that a class would be present in the scipy.sparse namespace as a class. But again, if registering a plugin format works, perhaps even that could be done.

Perhaps it is best to take small steps. Can we get the branch to work for scipy functions that accept existing sparse formats? But if you have insight into how a registering system could work it’d be cool to think about that too. :}

[jorenham]

@dschult

I suppose monkeypatching in some way might work – and maybe there would be a way to register a plugin format that asformat could use to try something like a new_format.fromcsr. The outside package would need to define to/from methods.

That’s not something that can be expressed in Python’s type system, unfortunately.

[dschult]

@cmclean do I have your permission to use your commit to play with creating a plugin system for externally supported formats?

[cmclean]

Apologies for the delay, I was away without email for the entire last week and missed this message. Yes, please feel free to use this code however you’d like!