5.4.13Scientific Computing (Python)

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

2,010 words9 min readdifficulty · medium

WHY sparse at all?

WHY it matters: A dense 106×10610^6\times 10^6 matrix needs 101210^{12} floats =8=8 TB — impossible. If each row has ~5 nonzeros, nnz=5×106\text{nnz}=5\times10^6, storable in ~120 MB. Same problem, 105×\sim 10^5\times less memory.


The three formats you must know

1. COO — Coordinate (the "scratchpad")

2. CSR — Compressed Sparse Row

3. CSC — Compressed Sparse Column


Deriving CSR from scratch (Feynman-style)

Take A=(1000002000300405006000).A=\begin{pmatrix}10 & 0 & 0 & 0\\ 0 & 20 & 0 & 0\\ 30 & 0 & 40 & 50\\ 0 & 60 & 0 & 0\end{pmatrix}.

Step 1 — list nonzeros row by row. Why? CSR's ordering convention is "scan each row left→right." Row0: (0,10)(0,10). Row1: (1,20)(1,20). Row2: (0,30),(2,40),(3,50)(0,30),(2,40),(3,50). Row3: (1,60)(1,60).

So data = [10, 20, 30, 40, 50, 60], indices = [0, 1, 0, 2, 3, 1].

Step 2 — count nonzeros per row. Counts =[1,1,3,1]=[1,1,3,1]. Why? indptr is the cumulative sum of these, telling us boundaries.

Step 3 — cumulative sum, starting at 0. indptr=[0,  0+1,  1+1,  2+3,  5+1]=[0,1,2,5,6].\text{indptr}=[0,\;0{+}1,\;1{+}1,\;2{+}3,\;5{+}1]=[0,1,2,5,6]. Why this step? Row ii lives in data[indptr[i]:indptr[i+1]]. Check row 2: data[2:5]=[30,40,50] ✓. The final entry =nnz=6=\text{nnz}=6.

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

Worked example: building & using it in code

import numpy as np
from scipy.sparse import coo_matrix, csr_matrix
from scipy.sparse.linalg import spsolve
 
# Build in COO (easy), convert to CSR (fast)
rows = [0, 1, 2, 2, 2, 3]
cols = [0, 1, 0, 2, 3, 1]
vals = [10,20,30,40,50,60]
A = coo_matrix((vals, (rows, cols)), shape=(4,4)).tocsr()
 
print(A.indptr)    # [0 1 2 5 6]   <- matches our derivation
print(A.indices)   # [0 1 0 2 3 1]
print(A.data)      # [10 20 30 40 50 60]

Why .tocsr()? COO can't do fast arithmetic or solves; CSR can.

Solving Ax=bAx=b

from scipy.sparse import diags
from scipy.sparse.linalg import spsolve, cg
 
n = 1000
# 1D Poisson: tridiagonal [-1, 2, -1]
A = diags([-1, 2, -1], [-1, 0, 1], shape=(n, n), format='csc')
b = np.ones(n)
 
x_direct = spsolve(A, b)          # sparse LU (uses CSC)
x_iter, info = cg(A, b)           # conjugate gradient (SPD only)

Why CSC for spsolve? Sparse LU factorization processes columns; passing CSC avoids an internal conversion. Why cg here? This AA is symmetric positive definite, so Conjugate Gradient converges fast and never forms a dense factorization.


Common mistakes (steel-manned)


Active recall

Recall Quick self-test (cover the answers)
  • What three arrays define CSR? → data, indices, indptr.
  • Length of indptr? → n+1n+1.
  • Where do row ii's nonzeros live? → data[indptr[i]:indptr[i+1]].
  • Which format do direct solvers prefer? → CSC.
  • SpMV cost? → O(nnz)O(\text{nnz}).
Recall Feynman: explain to a 12-year-old

Imagine a giant spreadsheet that's almost all empty boxes. Instead of writing down every empty box, you only write the boxes that have a number — plus a little note saying "row 2 starts here, row 3 starts here." That note is indptr: it's like bookmarks telling you where each row's numbers begin, so you can flip straight to any row. Because you skipped all the blanks, the computer reads way less and finishes way faster.

What problem do sparse matrices solve?
Matrices that are mostly zeros waste O(n²) memory/compute; sparse storage keeps only nonzeros + coordinates, giving O(nnz) cost.
What three arrays make up CSR format?
data (nonzeros row-by-row), indices (column index of each nonzero), indptr (length n+1, row-start offsets).
How do you get the nonzeros of row i in CSR?
data[indptr[i]:indptr[i+1]] at columns indices[indptr[i]:indptr[i+1]].
What is the length and last value of CSR indptr?
Length n+1; last value equals nnz (total nonzeros).
How does CSC differ from CSR?
CSC stores column-by-column: indices hold row indices, indptr marks column starts; fast column slicing and used by direct solvers.
Why build a sparse matrix in COO (or LIL/DOK) then convert?
COO/LIL allow cheap incremental insertion; CSR/CSC make insertion O(nnz). Build once, convert once.
What is the cost of sparse matrix-vector product y=Ax in CSR?
O(nnz), summing data[k]*x[indices[k]] over each row's slice.
When should you NOT use a sparse matrix?
When density is high (≳10%): index overhead (2·nnz ints) can make sparse bigger and slower than dense.
Which scipy function does a direct sparse solve and what format does it like?
spsolve (sparse LU); prefers CSC input.
When is cg (conjugate gradient) valid, and what to use otherwise?
cg needs symmetric positive-definite A; use gmres/bicgstab for general matrices.
What is fill-in?
Zeros that become nonzero during LU/Cholesky factorization, increasing memory; iterative solvers avoid it.
In modern scipy.sparse, how do you do matrix multiply vs elementwise multiply?
A @ x for matrix product; A.multiply(B) for elementwise.

Connections

  • Dense matrices (NumPy ndarray) — the baseline sparse improves on.
  • LU decomposition and Cholesky decomposition — what spsolve does under the hood.
  • Conjugate Gradient method / Krylov subspace methods — iterative cg, gmres.
  • Preconditioning — accelerates iterative sparse solvers.
  • Finite element method / PDE discretization — main source of sparse systems.
  • Graph adjacency matrices — graphs are naturally sparse.
  • Big-O complexity — why O(nnz) beats O(n²).

Concept Map

quantified by

small rho pays off

build easily in

convert to

convert to

compresses row array via

fast right-multiply

transpose of

preferred by

Sparse matrix mostly zeros

Density rho = nnz over n*m

COO coordinate triples

CSR compressed row

CSC compressed column

indptr cumulative row starts

Matrix-vector Ax

Direct solvers LU Cholesky

Huge memory savings

Hinglish (regional understanding)

Intuition Hinglish mein samjho

Dekho, bahut badi matrices — jaise PDE solve karte waqt ya graphs me — mostly zeros hoti hain. Agar tum saare zeros bhi store karoge to memory aur time dono waste honge (n2n^2 cost). Sparse matrix ka funda simple hai: sirf nonzero values aur unke coordinates store karo. Isse 101210^{12} numbers ki jagah sirf nnz\text{nnz} numbers bachte hain — kabhi-kabhi 8 TB ki jagah 100 MB!

CSR (Compressed Sparse Row) format me teen arrays hote hain: data (nonzeros, row-by-row left to right), indices (har value ka column number), aur indptr (length n+1n+1, jo batata hai ki har row data me kahan se start hoti hai). indptr actually bookmark jaisa hai — row ii ke nonzeros milte hain data[indptr[i]:indptr[i+1]] me. Isi liye row access O(1)O(1) hota hai aur AxAx multiply super fast ho jata hai. CSC bilkul ulta hai — column-by-column — aur direct solvers (LU factorization) isko prefer karte hain.

Practical workflow yaad rakhna: matrix COO me build karo (entries dalna easy), phir .tocsr() ya .tocsc() se convert karo. CSR me directly A[i,j]=v karne ki galti mat karna — wo har baar arrays shift karta hai, slow aur warning deta hai. Solve karne ke liye: general matrix ke liye spsolve (CSC do), symmetric positive-definite ke liye cg, aur baaki iterative ke liye gmres/bicgstab. Bas yaad rakho: "COO to collect, CSR to compute, CSC to crack." Ye samajh loge to scientific computing me bade-bade systems easily handle kar paoge.

Go deeper — visual, from zero

Test yourself — Scientific Computing (Python)

Connections