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

cmclean · July 14, 2026, 1:15pm

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 · July 15, 2026, 1:05am

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 · July 15, 2026, 8:01am

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 · July 15, 2026, 2:21pm

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.