An 8 × 8 matrix with only 10 entries that are not zero. Watch them leave the grid as a short list and fold into three compact arrays.
Storing only the non-zero entries
Ten numbers of an 8 × 8 matrix, as triples and as compressed rows
void to_csr(const struct triple *t, int k, double *val, int *col, int *start){ for (int r = 0; r <= N; r++) start[r] = 0; for (int i = 0; i < k; i++) start[t[i].row + 1]++; for (int r = 0; r < N; r++) start[r + 1] += start[r]; for (int i = 0; i < k; i++) { val[i] = t[i].val; col[i] = t[i].col; }} void csr_mul(const double *val, const int *col, const int *start, const double *x, double *y){ for (int r = 0; r < N; r++) { double s = 0.0; for (int p = start[r]; p < start[r + 1]; p++) s += val[p] * x[col[p]]; y[r] = s; }}int main(void){ const struct triple a[10] = { { 0, 0, 4 }, { 0, 5, 1 }, { 1, 1, 3 }, { 2, 0, 2 }, { 2, 6, 5 }, { 3, 3, 7 }, { 4, 7, 6 }, { 5, 2, 8 }, { 6, 6, 9 }, { 7, 4, 1 } }; double val[MAXK], x[N], y[N]; int col[MAXK], start[N + 1]; to_csr(a, 10, val, col, start); printf("start:"); for (int r = 0; r <= N; r++) printf(" %d", start[r]); for (int j = 0; j < N; j++) x[j] = j + 1; csr_mul(val, col, start, x, y); printf("\ny ="); for (int r = 0; r < N; r++) printf(" %g", y[r]); const struct triple b[4] = { { 0, 5, -1 }, { 2, 2, 3 }, { 6, 6, 1 }, { 7, 0, 2 } }; struct triple sum[MAXK]; int k = add_triples(a, 10, b, 4, sum); printf("\nA + B has %d entries:", k); for (int i = 0; i < k; i++) printf(" (%d,%d,%g)", sum[i].row, sum[i].col, sum[i].val); printf("\n"); return 0;}An 8 × 8 matrix with 10 non-zero entries, written as triples in row order.
- numbers stored
- 0 for 64 cells
Output
Describe this picture
The picture is an 8 × 8 grid with the row and column numbers beside it, and under it three strips of boxes. In the first view, “Triples”, the strips are row, col and val, 10 boxes each, the triples of the array a: box k of the three strips holds entry k’s row, column and value. Each grid cell with an entry is drawn solid and tied to its boxes; the other 54 cells are blank. Each row of the grid has a faint tint, carried down to its boxes. Under the grid is a line of text with the current formula. Under the picture are the lines of sparse.c the clip runs, with the line about to run marked, one sentence on the current step, a readout “numbers stored”, the control, and the output. The clip plays once, about 17 seconds. The grid fills with its 10 entries, and each one flies down into its triple: rows 0 0 1 2 2 3 4 5 6 7, columns 0 5 1 0 6 3 7 2 6 4, values 4 1 3 2 5 7 6 8 9 1. Then to_csr counts the entries of each row, takes running totals and copies the values and columns. Next csr_mul multiplies by x = 1, 2, …, 8 a row a beat, lighting the row’s entries: y[0] = 4·x[0] + 1·x[5] = 10, and so on to y[7] = 5. Then the four entries of a second matrix B light, at (0, 5), (2, 2), (6, 6) and (7, 0), and the grid becomes A + B: (0, 5) has cancelled and gone, (2, 2) and (7, 0) are new, and (6, 6) is 10. The strips now hold the 11 triples of sum, and the clip rests on the output. The control, “Stored as”, switches to “Compressed rows”: the strips become val and col, 10 boxes each, and start, 9 boxes. They hold unknown values until to_csr fills them; the values fly from the grid into val when it copies them, and start reads 0 2 3 5 6 7 8 9 10. Then an arrow runs from each row’s start up to the position of that row’s first entry; while csr_mul works on a row, only that row’s arrow shows. The addition works on triples in both views. The readout counts the numbers in the strips at the current step, against the grid’s 64 cells: 30 once the triples are written, 29 once the compressed rows are, and 33 for the 11 triples of the sum. The transport plays and pauses, steps one beat back or forward, and scrubs; a paused step is kept in the page address as sparse.step, and the view as sparse.form.
Only the entries that are there
A sparse matrix is one in which almost every entry is zero. The example has 10 non-zero entries out of 64. Storing it in full keeps 54 zeros; storing only the 10 entries needs a way to say where each one is.
Triples
The simplest form, the coordinate form, writes each entry as a triple, its row, its column and its value:
struct triple {
int row, col;
double val;
};The 10 entries are 10 triples, 30 numbers. The program keeps them sorted by row, then by column. row * N + col is that order as one number, so two triples compare by it: (2, 6) is 22, before (3, 3) at 27.
Compressed rows
In sorted triples the row numbers repeat: 0 0 1 2 2 3 4 5 6 7. Compressed sparse row form, CSR, keeps the values and the columns, k of each, and replaces the k row numbers by n + 1 row starts: start[r] is the position of row r’s first entry, so row r’s entries are positions start[r] to start[r + 1] − 1.
val: 4 1 3 2 5 7 6 8 9 1
col: 0 5 1 0 6 3 7 2 6 4
start: 0 2 3 5 6 7 8 9 10Row 2’s entries are positions start[2] = 3 to start[3] − 1 = 4: the 2 in column 0 and the 5 in column 6. CSR does not store the row of each entry, only where each row starts. That is 10 + 10 + 9 = 29 numbers, against 30 for triples and 64 for the full grid.
The code and its cost
Building the row starts
to_csr makes the starts in three loops: clear them, count each row’s entries into the slot after it, and take running totals.
for (int r = 0; r <= N; r++)
start[r] = 0;
for (int i = 0; i < k; i++)
start[t[i].row + 1]++;
for (int r = 0; r < N; r++)
start[r + 1] += start[r];After counting, start holds 0 2 1 2 1 1 1 1 1; after the running totals, 0 2 3 5 6 7 8 9 10, each the number of entries in the rows above. Counting sort, a later lesson, uses the same running totals. The values and columns are then copied across unchanged, since the triples are already in row order. The program prints start: 0 2 3 5 6 7 8 9 10.
Multiplying by a vector
For each row, csr_mul runs over that row’s positions only:
for (int r = 0; r < N; r++) {
double s = 0.0;
for (int p = start[r]; p < start[r + 1]; p++)
s += val[p] * x[col[p]];
y[r] = s;
}col[p] picks the entry of x that the value meets. Each stored entry is used once, so the product takes k multiplications, here 10, where the full grid takes n² = 64. With x = (1, 2, …, 8) the program prints
y = 10 6 37 28 48 24 63 5Adding two matrices
Two sorted triple lists add by the merge of merging arrays, comparing positions row * N + col:
int ka = a[i].row * N + a[i].col;
int kb = b[j].row * N + b[j].col;
if (ka < kb) {
out[k++] = a[i++];
} else if (kb < ka) {
out[k++] = b[j++];
} else {
double s = a[i].val + b[j].val;
if (s != 0.0) {
out[k] = a[i];
out[k].val = s;
k++;
}
i++;
j++;
}When both lists have an entry at the same position, the values add, and a sum of 0 is dropped, so the result stays free of stored zeros. B = (0, 5, −1), (2, 2, 3), (6, 6, 1), (7, 0, 2) cancels A’s (0, 5) entry, adds to (6, 6) and brings two new entries:
A + B has 11 entries: (0,0,4) (1,1,3) (2,0,2) (2,2,3) (2,6,5) (3,3,7) (4,7,6) (5,2,8) (6,6,10) (7,0,2) (7,4,1)The merge takes at most k₁ + k₂ steps, here 14.
Transposing
To transpose a matrix in CSR, do what to_csr does with the columns instead of the rows: count the entries in each column, take running totals to get where each column starts, then place each entry at its column’s next free position. Three loops, each over k or n.
When it pays
Triples take 3k numbers and CSR 2k + n + 1. Sparse storage pays while that is less than n², the full grid. For the example, 30 and 29 against 64. A matrix half full of non-zeros is cheaper to store in full.
Running it yourself
gcc -std=c11 -Wall -Wextra -pedantic sparse.c -o sparse
./sparseWorked example
This matrix. Full: 64 numbers. Triples: 3 × 10 = 30. CSR: 2 × 10 + 8 + 1 = 29.
A web of pages. Take 1,000,000 web pages with 10 links each, as a matrix with a 1 where page i links to page j. It has 10⁷ non-zero entries. CSR needs about 2 × 10⁷ + 10⁶ numbers; the full grid has 10¹² cells.
One entry of y. Row 2 has entries 2 in column 0 and 5 in column 6, so y[2] = 2 × x[0] + 5 × x[6] = 2 × 1 + 5 × 7 = 37.
Two checks. Which positions of val hold row 2? (3 and 4.) How many entries does row 6 of A + B have? (1, the value 10.)
Where it’s used
Circuits and structures. Simulating a circuit or a building frame gives a matrix in which each node touches only a few others.
Web pages. The link matrix of the web, used to rank pages, is far too large to store in full.
Ratings. A table of the films each user has rated is mostly empty: each user rates a handful of films.
Solvers. Iterative solvers need only products of the matrix with a vector, which CSR does in k steps. A graph stored as adjacency lists, a later lesson, is CSR’s rows; see also graphs and networks.
Reference card
| Code or idea | Meaning |
|---|---|
| (row, col, val) | a triple: 3k numbers for k entries |
row * N + col | the sort order of triples as one number |
val, col, start | CSR: 2k + n + 1 numbers |
start[r] to start[r + 1] - 1 | the positions of row r |
s += val[p] * x[col[p]] | CSR times a vector: k multiplications |
merge on row * N + col | adding triple lists; a zero sum is dropped |
Next, sparse matrices as linked lists let entries be inserted and removed without shifting.
Platform: C11, gcc on x86-64 Linux: int 4 bytes, double 8, pointers 8, little-endian. Addresses in the pictures are short stand-ins, by one rule for the whole course; sizes are the real ones.