Stop Wasting Memory: The Ultimate Guide to Sparse Matrix Optimization

Introduction

Sparse matrices are a critical concept in modern computer science and data science. They are used to represent large datasets characterized by a high proportion of zero or missing values. In fields ranging from social network analysis(where most users aren’t connected to each other) to finite element analysis and machine learning, storing every single zero is not just inefficient — it is often computationally impossible. By leveraging the “sparsity” of these datasets, we can design specialized data structures that drastically reduce memory usage and accelerate processing speeds.

> Read whole story for free <

What is sparse matrix ?

Let’s consider a matrix M n x m, where n is number of rows and m is number of columns. Each element mij in M (where 0 ≤ i ≤ n-1 and 0 ≤ j ≤ m-1) is defined as a real number from R. We call matrix M sparse if and only if the density of non-zero elements is significantly lower then the number of zero or missing values.

How to represent it in memory ?

The simplest approach to representing a matrix in computer memory is to flatten it into a one-dimensional array. Since memory is linear, we map the 2D coordinates (i,j) to a 1D index using a formula. The most common convention is Row-Major Order, where the index is calculated as:

index = i * #columns + j

In C-style pseudo-code, a basic dense matrix structure looks like this:

struct Matrix
{
    int columns;
    int rows;
    float* elements;
}

While this is straightforward, it is highly inefficient for sparse matrices. If a 10,000 x 10,000 matrix contains only 100 non-zero values, this “dense” representation still allocates memory for all 100 million elements, 99.99% of which are zeros. This waste of space is what motivates the specialized sparse formats discussed below.

Optimize memory layout

Since a sparse matrix M contains only a small percentage of non-zero elements, representing it in a dense format is often wasteful. Instead, we use specialized formats designed to minimize memory footprint and optimize the performance of linear algebra operations like matrix-vector multiplication.

The most popular sparse matrix representations are:

• COO: Coordinate Format
• CSR: Compressed Sparse Row
• CSC: Compressed Sparse Column
• LIL: List of Lists
• SELL: Sliced Ellpack
• *Band Matrix
• *Diagonal Matrix

• This formats are variants of the SELL representation

To illustrate how these formats work, let’s consider a sample sparse matrix A 4 x 4 containing 7 non-zero elements and 9 zero elements (43.75% density):

Coordinate (COO)

Coordinate (COO) is the most intuitive format. It is essentially a list of triplets (row,column,value).

• Best for: Constructing matrices. Since you can append elements in any order, it’s ideal for the initial assembly of a dataset.
• Drawback: It is inefficient for mathematical operations because it requires storing two integers for every single value and offers no fast way to access specific rows or columns.

The Storage Structure

The COO format is defined by three primary arrays of length nnz (number of non-zero elements):

1. row_indices: The row position of each non-zero element.
2. column_indices: The column position of each non-zero element.
3. values: The actual numerical data.

The Mapping Algorithm

To convert a dense matrix M to COO format, follow these steps:

1. Initialize three empty lists: row_indices, column_indices, and values.
2. Iterate through the matrix row by row (i) and column by column (j).
3. For every element mij​:

• If mij =! 0:
• Append i to row_indices.
• Append j to column_indices.
• Append the value min​ to values.

4. Store the total dimensions (n,m) and the count of non-zero elements (nnz).

Example Matrix A in COO format:

The resulting arrays:

Compressed Sparse Row (CSR)

Compressed Sparse Row (CSR) is the “default” for most linear algebra libraries (like SciPy or Intel MKL). It is a compression of the COO format that specifically optimizes for row-wise traversal.

• Best for: Row-based operations, such as Matrix-Vector Multiplication (Ax=b). The row_offsets allow for fast, contiguous access to all elements in a single row.
• Hardware: Excellent for CPUs due to good cache locality when iterating through rows.

The Storage Structure

The CSR format is defined by three primary arrays of length nnz (number of non-zero elements):

1. row_offsets: An array of length n+1 (number of rows + 1). Each entry i indicates the index in the values array where row i begins. The last entry is always nnz.
2. column_indices: The column position of each non-zero element.
3. values: The actual numerical data.

The Mapping Algorithm

1. Initialize values and column_indices exactly as in COO (sorted by row).
2. Initialize row_offsets with a size of n+1.
3. Set row_offsets[0] = 0.
4. For each row i, set row_offsets[i+1] equal to the total number of non-zero elements found in all rows up to i.

Example Matrix A in CSR format:

Compressed Sparse Column (CSC)

Compressed Sparse Column (CSC) is the transposed version of CSR.

• Best for: Column-based operations or when you need to slice data by features (common in statistics and machine learning).
• Hardware: Similar performance to CSR, but optimized for algorithms that iterate through columns first.

The Storage Structure

CSC uses three primary arrays to store the data:

1. **values**: An array of length nnz containing the non-zero elements, traversed column-by-column.
2. **row_indices**: An array of length nnz containing the row position for each value.
3. **column_offsets**: An array of length m+1 (number of columns + 1). Entry j indicates where column j starts in the values array.

The Mapping Algorithm

1. Iterate through the matrix column-by-column (j=0 \dots m−1).
2. For each column, identify the non-zero elements and store their values in the values array and their row positions in row_indices.
3. Construct the column_offsets array where column_offsets[j] stores the cumulative count of non-zero elements found in all columns prior to j.
4. The final entry, column_offsets[m], is always nnz.

Example Matrix A in CSC format:

List of Lists (LIL)

The List of Lists (LIL) format is a row-based sparse matrix representation. It is essentially a collection of linked lists or dynamic arrays, where each entry in the main list represents a row of the matrix.

• Best for: Incremental construction and modifications. Because it uses dynamic structures for each row, adding a new non-zero element is relatively fast compared to rigid formats like CSR. It also supports efficient row-slicing.
• Drawback: Like COO, it is inefficient for arithmetic operations (addition, multiplication). It also has a higher memory overhead than CSR because it stores many small list objects rather than a few large, contiguous arrays.

The Storage Structure

The LIL format is defined by two primary lists of lists, where the outer list has a length equal to the number of rows (n):

• rows: A list where each element i is another list containing the column indices of the non-zero elements in row i.
• data: A list where each element i is another list containing the actual numerical values corresponding to the column indices in the rows list.

The Mapping Algorithm

To convert a dense matrix M to LIL format, follow these steps:

1. Initialize two empty lists of lists, rows and data, each with a length equal to the number of matrix rows.
2. Iterate through the matrix row by row (i):

• For every element mij​ in the current row:
• If mij ​=! 0:
• Append the column index j to the i-th list in rows.
• Append the value mij​ to the i-th list in data.

3. Store the total dimensions (n,m) of the matrix.

Example Matrix A in LIL format:

Sliced Ellpack (SELL)

Sliced Ellpack (SELL) and its variants are designed to overcome the limitations of the “irregular” data found in CSR when running on SIMD (Single Instruction, Multiple Data) architectures like GPUs.

• Best for: High-performance computing (HPC) on GPUs. By padding rows to a uniform length within a “slice,” it ensures that threads in a GPU warp don’t idle while waiting for longer rows to finish.
• Drawback: If the matrix has very inconsistent row lengths, the padding can waste a significant amount of memory.

The Storage Structure

SELL describes the matrix using the following parameters:

1. Slice Height (C): The number of rows grouped into a single slice (often matched to the GPU warp size, e.g., 32 or 64).
2. **values**: An array containing the non-zero elements and the necessary zero-padding for each slice.
3. **column_indices**: The column positions for each entry in the values array (padded with −1 or a sentinel value).
4. **slice_offsets**: Indices indicating where each slice starts in the values and column_indices arrays.

Example Matrix A in SELL format:

Conclusion

If you are just starting, start with COO to build your matrix and then convert to CSR for computation. If you find your bottleneck is GPU memory bandwidth and your rows are of similar length, investigate SELL.

Ultimately, the goal of sparse representation is to turn an O(N²) problem into an O(nnz) problem — where nnz is the number of non-zero elements. By choosing the right format, you ensure that your code spends more time calculating and less time "hunting" for data in memory.