scipy.sparse — sparse matrix formats (CSR, CSC), sparse solvers: Exercises

Level 1 — Recognition

These test whether you can read the three CSR/CSC arrays and turn them back into a matrix.

L1.1 — Read a CSR triple back into a matrix

Given a $3 \times 3$ matrix in CSR:

data    = [5, 8, 3, 7]
indices = [0, 2, 1, 2]
indptr  = [0, 2, 3, 4]

Reconstruct the full $3 \times 3$ matrix.

Recall Solution

WHAT we do: slice data/indices per row using indptr. Recall row $i$ lives in data[indptr[i]:indptr[i+1]].

Row 0: indptr[0]:indptr[1] $=$ [0:2] → data[0:2]=[5,8] at columns indices[0:2]=[0,2]. So $A_{0, 0} = 5, A_{0, 2} = 8$ .
Row 1: [2:3] → data=[3] at column indices=[1]. So $A_{1, 1} = 3$ .
Row 2: [3:4] → data=[7] at column indices=[2]. So $A_{2, 2} = 7$ .

$A = 500030807 .$ Sanity check: indptr has length $n + 1 = 4$ ✓, and its last value $= 4 = nnz$ ✓ (there are 4 stored numbers).

L1.2 — Spot the illegal CSR

One of these indptr arrays cannot be valid CSR for a matrix with $nnz = 6$ and $n = 4$ rows. Which, and why?

(a) [0, 1, 2, 5, 6]
(b) [0, 2, 2, 4, 6]
(c) [0, 3, 1, 4, 6]

Recall Solution

The rule: indptr must be non-decreasing (each row's start is $\geq$ the previous row's start — you can't un-count nonzeros), start at 0, have length $n + 1 = 5$ , and end at $nnz = 6$ .

(a) [0,1,2,5,6]: increasing, ends at 6 ✓ valid (row counts $1, 1, 3, 1$ ).
(b) [0,2,2,4,6]: non-decreasing (a 2,2 means row 1 is empty — allowed), ends at 6 ✓ valid (counts $2, 0, 2, 2$ ).
(c) [0,3,1,4,6]: goes $3 \to 1$ , which decreases. That would say row 1 has $1 - 3 = - 2$ nonzeros — impossible. ❌ This is the invalid one.

Level 2 — Application

Now you produce the arrays yourself and use them.

L2.1 — Derive CSR by hand

Build the CSR arrays for $B = 010540000206 .$

Recall Solution

Step 1 — scan each row left→right, list (column, value):

Row 0: $(1, 4)$
Row 1: $(0, 1), (2, 2)$
Row 2: — empty —
Row 3: $(0, 5), (2, 6)$

So data = [4, 1, 2, 5, 6], indices = [1, 0, 2, 0, 2].

Step 2 — count nonzeros per row: $[1, 2, 0, 2]$ .

Step 3 — cumulative sum starting at 0: $indptr = [0, 0 + 1, 1 + 2, 3 + 0, 3 + 2] = [0, 1, 3, 3, 5] .$ Note the repeated 3,3 — that empty row 2 (see L1 mistake). Last value $= 5 = nnz$ ✓, length $= n + 1 = 5$ ✓.

See the scan happen. The figure below shows the row-by-row sweep over $B$ : a magenta arrow slides left→right through each row, dropping each nonzero into data/indices; the empty row 2 produces the tell-tale flat step in indptr.

Figure — scipy.sparse — sparse matrix formats (CSR, CSC), sparse solvers

L2.2 — CSC of the same matrix

Now build CSC arrays for the same $B$ above. (CSC scans column by column; indices holds row indices; indptr has length $m + 1$ .)

Recall Solution

Scan each column top→bottom, list (row, value):

Col 0: $(1, 1), (3, 5)$
Col 1: $(0, 4)$
Col 2: $(1, 2), (3, 6)$

So data = [1, 5, 4, 2, 6], indices = [1, 3, 0, 1, 3].

Per-column counts: $[2, 1, 2]$ → cumulative from 0: indptr = [0, 2, 3, 5]. Length $= m + 1 = 4$ ✓, ends at $5 = nnz$ ✓. Same 5 numbers as CSR, reordered by column.

See the difference. In the recap figure, CSR swept rows; here the sweep runs down each column. Same coloured cells, a different traversal order — that reordering is the entire distinction between CSR and CSC.

L2.3 — Compute $y = B x$ using the CSR SpMV formula

Using CSR of $B$ from L2.1 and $x = (1, 10, 100)^{⊤}$ , compute $y = B x$ with the row-slice formula $y_{i} = \sum_{k = indptr [i]}^{indptr [i + 1] - 1} data [k] \cdot x [indices [k]] .$

Recall Solution

Recall data=[4,1,2,5,6], indices=[1,0,2,0,2], indptr=[0,1,3,3,5].

$y_{0}$ : $k \in [0, 1)$ → $4 \cdot x [1] = 4 \cdot 10 = 40$ .
$y_{1}$ : $k \in [1, 3)$ → $1 \cdot x [0] + 2 \cdot x [2] = 1 \cdot 1 + 2 \cdot 100 = 201$ .
$y_{2}$ : $k \in [3, 3)$ → empty sum $= 0$ .
$y_{3}$ : $k \in [3, 5)$ → $5 \cdot x [0] + 6 \cdot x [2] = 5 + 600 = 605$ .

$y = (40, 201, 0, 605)^{⊤} .$ Only 5 multiply–adds were needed (one per nonzero) — cost $O (nnz)$ , not $O (nm) = 12$ .

See the routing. The figure below draws the SpMV as a wiring diagram: each nonzero of $B$ picks the entry of $x$ named by its indices value, multiplies, and drops the product into the $y$ slot named by its row. Row 2 has no wires → $y_{2} = 0$ .

Level 3 — Analysis

Reason about memory, cost, and structure.

L3.1 — Memory crossover

For an $n \times n$ CSR matrix, a dense array stores $n^{2}$ floats. CSR stores $2 nnz + (n + 1)$ numbers (nonzeros data + indices + indptr). Take $n = 1000$ . Below what density $ρ$ does CSR use fewer numbers than dense? Ignore the small $(n + 1)$ term for the estimate, then refine.

Recall Solution

Estimate (drop $n + 1$ ): CSR cheaper when $2 nnz < n^{2} ⟺ nnz < \frac{n ^{2}}{2} ⟺ ρ = \frac{nnz}{n ^{2}} < \frac{1}{2} = 0.5.$ Refine (keep $n + 1$ ): solve $2 nnz + (n + 1) < n^{2}$ : $nnz < \frac{n ^{2} - ( n + 1 )}{2} = \frac{1 0 ^{6} - 1001}{2} = 499499.5,$ so $nnz \leq 499499$ , i.e. $ρ < 0.499499$ . Interpretation: by a raw count of stored numbers, CSR wins below ~50% density. But indices/indptr are integers and each nonzero also needs its integer column — in bytes, and because sparse SpMV has worse cache behaviour per nonzero, the practical rule of thumb from the parent note ( $ρ ≲ 0.1$ ) is much stricter. The crossover in numbers is generous; the crossover in wall-clock speed is far tighter. See Big-O complexity.

See the crossover. The plot below draws both storage costs against density for $n = 1000$ : a flat dense line at $n^{2}$ and a rising CSR line $2 ρ n^{2} + (n + 1)$ . Where they cross ( $ρ \approx 0.4995$ ) is the count-crossover; the shaded band near $ρ ≲ 0.1$ marks where sparse actually pays off in speed.

L3.2 — Row access vs column access cost

In CSR, getting all of row 7 is $O (1)$ to locate plus $O (nnz in that row)$ to read. What is the cost of getting all of column 7 in CSR, and why? What format fixes it?

Recall Solution

Column access in CSR is slow: the nonzeros of column 7 are scattered arbitrarily across data — CSR groups by row, so to find every entry with indices[k]==7 you must scan all $nnz$ entries. Cost $= O (nnz)$ , and it also costs a search inside each row. Fix: use CSC, which groups by column, making column 7 a single contiguous slice data[indptr[7]:indptr[8]] — $O (1)$ to locate. This is exactly why direct solvers (LU decomposition, Cholesky decomposition) prefer CSC: they eliminate column by column.

L3.3 — Fill-in

You factor a sparse LU of an arrow matrix: nonzeros only on the diagonal, the entire first row, and the entire first column ( $n = 5$ ). Sketch (in words) why eliminating the first column before reordering creates a dense trailing block (fill-in), and how reversing the arrow avoids it.

Recall Solution

The arrow (point at top-left): $A = * * * * * * * * * * * * * .$ Read the figure left→right. The left panel shows the original arrow: violet cells are the only stored nonzeros (*) — the full first row, full first column, and the diagonal. The right panel shows what happens after eliminating column 0 first: to zero each $A_{i, 0}$ (rows 1–4) we subtract a multiple of the full row 0, which injects a new nonzero into every column 1–4 of those rows. Those brand-new entries are drawn in orange and marked + — the previously-empty $4 \times 4$ trailing block has become completely dense. The magenta arrow points straight at that catastrophic block.

Reverse the arrow (permute so the dense row/column is last): $P A P^{⊤} = * * * * * * * * * * * * * .$ Now eliminating columns 0–3 each touches only its own diagonal and the last row — no new fill until the final pivot. A good reordering (AMD, RCM) chases exactly this. Takeaway: fill-in depends on ordering, not just on the matrix.

Level 4 — Synthesis

Combine format choice, solver choice, and code.

L4.1 — Choose the solver

For each system, pick the best solver from {spsolve, cg, gmres, lsqr} and justify:

(a) $A$ is $5000 \times 5000$ , symmetric, all eigenvalues positive, sparse.
(b) $A$ is $5000 \times 5000$ , sparse, non-symmetric, you need an exact answer and it fits in memory.
(c) $A$ is $8000 \times 3000$ (tall, non-square), solve in the least-squares sense.
(d) $A$ is huge, non-symmetric, memory-limited.

Recall Solution

(a) cg — symmetric positive-definite (SPD) is exactly CG's home turf; it only needs matrix–vector products $A x$ and no factorization, so it's cheap and memory-light.
(b) spsolve — non-symmetric but you want an exact, robust answer and it fits: direct sparse LU (feed it CSC).
(c) lsqr — non-square, so "solve" means minimise $∥ A x - b ∥_{2}$ ; lsqr/lsmr are built for that.
(d) gmres (or bicgstab) — general non-symmetric, low memory: a Krylov iterative method, optionally with a preconditioner to speed convergence.

L4.2 — Build a 1-D Poisson matrix and count its nonzeros

The 1-D Poisson discretization gives a tridiagonal matrix with stencil $[- 1, 2, - 1]$ : value $2$ on the diagonal, $- 1$ on the sub- and super-diagonal, for size $n \times n$ . (i) How many nonzeros does it have? (ii) What is its density for $n = 1000$ ? (iii) Write the scipy.sparse one-liner to build it in CSC.

Recall Solution

(i) Count the nonzeros. The diagonal has $n$ entries. The sub-diagonal (just below) has $n - 1$ entries, and the super-diagonal (just above) also has $n - 1$ . There is no overlap between these three bands, so $nnz = n + (n - 1) + (n - 1) = 3 n - 2.$ (ii) Density for $n = 1000$ . Plug in: $nnz = 3 (1000) - 2 = 2998$ . The full matrix has $n^{2} = 1 0^{6}$ cells, so $ρ = \frac{nnz}{n ^{2}} = \frac{2998}{1000000} = 2.998 \times 1 0^{- 3} \approx 0.003 = 0.3%.$ Extremely sparse — a dense store would waste 99.7% of its memory on zeros. (iii) One-liner in CSC.

from scipy.sparse import diags
A = diags([-1, 2, -1], [-1, 0, 1], shape=(1000, 1000), format='csc')

The list [-1, 0, 1] gives the three diagonal offsets (below main, main, above main); format='csc' builds it directly in CSC so spsolve won't need an internal conversion.

L4.3 — Predict `indptr` for the tridiagonal (small case)

For the $[- 1, 2, - 1]$ tridiagonal matrix with $n = 4$ , write out data, indices, indptr in CSR by hand.

Recall Solution

$T = 2 - 1 00 - 1 2 - 1 0 0 - 1 2 - 1 00 - 1 2 .$ Scan rows left→right:

Row 0: $(0, 2), (1, - 1)$
Row 1: $(0, - 1), (1, 2), (2, - 1)$
Row 2: $(1, - 1), (2, 2), (3, - 1)$
Row 3: $(2, - 1), (3, 2)$

data = [2, -1, -1, 2, -1, -1, 2, -1, -1, 2] indices = [0, 1, 0, 1, 2, 1, 2, 3, 2, 3] Per-row counts $[2, 3, 3, 2]$ → indptr = [0, 2, 5, 8, 10]. Check: $nnz = 3 n - 2 = 10$ = last indptr value ✓.

See the mapping. The figure below colour-codes the tridiagonal $T$ and draws, for each stored nonzero, an arrow into its exact slot of data/indices, with the indptr bookmarks marking each row boundary. This is the whole "matrix → three arrays" translation made literal.

Level 5 — Mastery

Derive and design.

L5.1 — Derive the SpMV cost from the CSR structure

Prove that computing $y = A x$ via CSR costs exactly $nnz$ multiplications and $nnz$ additions (ignoring loop overhead), and explain why this beats dense $O (n^{2})$ .

Recall Solution

Setup. For each output entry, $y_{i} = \sum_{k = indptr [i]}^{indptr [i + 1] - 1} data [k] \cdot x [indices [k]] .$ The inner loop for row $i$ runs over $r_{i} = indptr [i + 1] - indptr [i]$ terms — that's the number of nonzeros in row $i$ . Each term is one multiply ( $data [k] \cdot x [\dots]$ ) and one add into the running sum.

Sum over all rows. Total multiplies $\sum_{i = 0}^{n - 1} r_{i} = \sum_{i = 0}^{n - 1} (indptr [i + 1] - indptr [i]) = indptr [n] - indptr [0] = nnz - 0 = nnz,$ a telescoping sum — the internal terms cancel, leaving last minus first. Adds are the same count. ∎

See the cancellation. The figure below draws the telescoping sum as a chain of indptr values: each term subtracts the previous bookmark, so every interior value is added once and subtracted once (they annihilate in pairs), leaving only $indptr [n] - indptr [0] = nnz$ .

Why it wins: dense multiply touches every one of the $n^{2}$ entries, doing $n^{2}$ multiplies. When $nnz = O (n)$ (e.g. 5 per row from a stencil), CSR does $O (n)$ work — a factor $\sim n$ faster. For $n = 1 0^{6}$ that's a million-fold saving.

L5.2 — Design: symmetric storage question

You have a symmetric sparse matrix ( $A = A^{⊤}$ ). A colleague says "store only the upper triangle to halve memory." What breaks in a generic CSR SpMV if you do this naively, and how do you fix the SpMV to still compute the full $y = A x$ ?

Recall Solution

What breaks: if you store only entries with $j \geq i$ , then the plain row-slice SpMV computes only the contribution of the upper triangle. Each off-diagonal nonzero $A_{ij}$ ( $i < j$ ) has a mirror $A_{j i} = A_{ij}$ that is no longer stored, so its contribution $A_{j i} x_{i}$ to $y_{j}$ is silently dropped. Result: wrong $y$ for every row that had a mirror entry below the diagonal.

Fix (symmetric SpMV): for each stored nonzero $A_{ij}$ apply it twice by hand: $y_{i} + = A_{ij} x_{j}, and if i \neq = j : y_{j} + = A_{ij} x_{i} .$ The diagonal ( $i = j$ ) is applied once. This restores the full product using half the storage — the standard trick in FEM and Cholesky codes, where matrices are symmetric by construction.

L5.3 — Design a preconditioner rationale

You solve $A x = b$ with CG on an SPD matrix whose condition number is $κ = 1 0^{6}$ ; convergence crawls. Explain, in terms of the CG convergence bound, why a preconditioner $M \approx A$ helps, and what $M = I$ (no preconditioning) vs $M = A$ (perfect) represent as extremes.

Recall Solution

The bound. CG's error after $t$ iterations shrinks roughly like $\frac{∥ e _{t} ∥ _{A}}{∥ e _{0} ∥ _{A}} \leq 2 (\frac{κ - 1}{κ + 1})^{t} .$ Large $κ$ pushes that ratio near 1, so each step barely reduces the error → slow.

Preconditioning replaces $A$ with $M^{- 1} A$ (via a symmetric split), whose condition number $κ (M^{- 1} A)$ we want small. Since $κ$ appears in the bound, cutting $κ$ from $1 0^{6}$ to, say, $1 0^{2}$ changes $κ$ from $1000$ to $10$ — dramatically faster convergence.

Extremes:

$M = I$ : $M^{- 1} A = A$ , no change — plain CG, $κ = 1 0^{6}$ , crawls.
$M = A$ : $M^{- 1} A = I$ , condition number $1$ — CG converges in one step, but forming $M^{- 1} = A^{- 1}$ is the original solve, so it's useless as work.

The art is choosing $M$ cheap to invert yet close to $A$ (e.g. incomplete Cholesky) — a compromise between these two extremes.

Final active recall

Recall One-line answers (cover them)

Invalid indptr symptom? ::: It decreases somewhere (or wrong length / wrong final value).
Empty row in CSR looks like? ::: Two equal consecutive indptr values.
CSC indices holds? ::: row indices.
Column access in CSR costs? ::: $O (nnz)$ — use CSC instead.
SpMV cost, derived how? ::: Telescoping sum of per-row counts $= nnz$ .
$[- 1, 2, - 1]$ tridiagonal nnz? ::: $3 n - 2$ .
CG requires? ::: symmetric positive-definite.
Preconditioner goal? ::: shrink $κ (M^{- 1} A)$ while keeping $M$ cheap to invert.