PYTHON / SCIPY
Sparse matrices and signal processing
Build and multiply SciPy sparse matrices in the right format, and see how an FIR filter is the same thing as a banded matrix.
What you will learn
- Assemble matrices in COO/LIL, convert to CSR or CSC before arithmetic or solves
- Read a CSR matrix's data, indices and indptr arrays to know what is stored
- Express a difference or moving-average filter as a banded sparse operator
- Use scipy.signal.lfilter and find_peaks on 1-D arrays and predict the edge behaviour
Understanding Sparse matrices and signal processing
A SciPy sparse matrix does not store zeros. The CSR format keeps three arrays: data with the nonzero values row by row, indices with the column of each value, and indptr with n+1 offsets saying where each row's slice begins and ends. That is why row slicing and matrix-vector products are fast (you walk contiguous memory once per nonzero) and why inserting a new entry is slow: a value in the middle of row 3 forces every later element of data and indices to shift and every later indptr entry to increase.
The formats exist because no single layout is good at everything. COO stores plain (row, col, value) triples, so appending is free and duplicates are allowed and summed on conversion; LIL keeps a Python list per row, which makes A[i, j] = v cheap. CSR is built for row-oriented arithmetic, CSC for column operations and for the LU factorisations behind spsolve, and DIA for matrices that are pure diagonals. The normal workflow is build in COO or LIL, call .tocsr() or .tocsc() once, then compute.
The link to signal processing is that a finite impulse response filter is a linear operator with the same coefficients repeated on every row, that is, a banded Toeplitz matrix. Applying scipy.signal.lfilter and multiplying by that banded sparse matrix compute the same sums; the filter routine is leaner because it never stores the band, while the matrix form lets you compose the operator with others, transpose it, or put it inside a least-squares or smoothing problem and solve. The one place the two views can differ is the boundary: lfilter assumes zero initial state, and a matrix built to a fixed shape encodes whatever edge convention you gave it.
import numpy as np
from scipy import sparse, signal
x = np.array([0.0, 1.0, 3.0, 6.0, 10.0, 15.0, 21.0, 28.0])
n = x.size
# FIR filter: y[k] = x[k] - x[k-1]
b = np.array([1.0, -1.0])
y_filter = signal.lfilter(b, [1.0], x)
# the same operator as a banded sparse matrix
D = sparse.diags([-1.0, 1.0], [-1, 0], shape=(n, n), format="csr")
y_matrix = D @ x
print(y_filter)
print(y_matrix)
print(np.allclose(y_filter, y_matrix))
print(D.format, D.nnz, D.shape[0] * D.shape[1])A sparse matrix is a compressed description of a mostly-zero linear operator, and most fixed-coefficient filters are exactly that kind of operator.
Worked examples
Inspecting CSR internals
Shows how COO triples with duplicates become canonical CSR arrays.
import numpy as np
from scipy import sparse
rows = np.array([0, 0, 1, 2, 2, 2])
cols = np.array([0, 2, 1, 0, 0, 2])
vals = np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0])
M = sparse.coo_matrix((vals, (rows, cols)), shape=(3, 3))
print(M.nnz)
C = M.tocsr()
print(C.nnz)
print(C.toarray())
print(C.data)
print(C.indices)
print(C.indptr)Example explained
Line 1COO reports nnz 6 because it stores the two triples that both target (2, 0) separately.
Line 2tocsr() sums those duplicates into the single value 9, so nnz drops to 5.
Line 3data holds the values in row-major order and indices holds their columns.
Line 4indptr [0, 2, 3, 5] says row 0 uses data[0:2], row 1 uses data[2:3], row 2 uses data[3:5].
Solving a banded system
Solves a tridiagonal system with spsolve and checks the recovered vector.
import numpy as np
from scipy import sparse
from scipy.sparse.linalg import spsolve
n = 5
A = sparse.diags([-1.0, 2.0, -1.0], [-1, 0, 1], shape=(n, n), format="csc")
x_true = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
b = A @ x_true
print(b)
x = spsolve(A, b)
print(np.round(x, 10))
print(A.nnz, n * n)Example explained
Line 1A is the second-difference operator, so A @ x_true is zero everywhere the sequence is linear.
Line 2format="csc" is chosen because spsolve factorises column-wise and warns if given another format.
Line 3The LU factors contain values like 4/3 that are inexact in binary, so np.round trims errors near 1e-16.
Line 4A.nnz is 13 of 25 possible entries: 5 on the diagonal plus 4 above and 4 below.
Locating peaks in a signal
Uses find_peaks with a height threshold and reads the returned properties.
import numpy as np
from scipy import signal
x = np.array([0.0, 2.0, 1.0, 5.0, 1.0, 4.0, 0.0, 1.0, 0.5])
peaks, props = signal.find_peaks(x, height=1.5)
print(peaks)
print(props["peak_heights"])Example explained
Line 1find_peaks returns sample indices, not values, so x[peaks] gives the amplitudes.
Line 2Index 7 is a genuine local maximum but its height of 1.0 fails the height=1.5 filter.
Line 3Passing height also makes the function return a peak_heights entry in the properties dict.
Line 4The first and last samples can never be peaks because they have only one neighbour.
Important notes
CSR costs roughly one value plus one index per nonzero plus one offset per row, so a matrix with more than about a third of its entries filled can use more memory than the dense version.
scipy.sparse now ships both the older matrix classes and the newer array classes; * means matrix product on the former and elementwise product on the latter, so write @ when you mean a matrix product.
Common mistakes
Filling a csr_matrix with A[i, j] = v inside a loop: SciPy raises SparseEfficiencyWarning and each assignment rebuilds the internal arrays, turning an O(nnz) build into O(nnz squared).
Passing a sparse matrix to np.dot or np.array: NumPy treats it as an opaque object and you get an object-dtype result or a wrong answer instead of the product, so always use A @ x or A.toarray() first.
Comparing lfilter output against a hand-built matrix product and blaming the matrix for a mismatch in the first few samples, when the real cause is that lfilter starts from a zero initial state.
Try it yourself
Change, predict, then run
Build the 10x10 second-difference matrix with sparse.diags([1, -2, 1], [-1, 0, 1]), apply it to np.arange(10)**3, and confirm the interior rows match np.diff(x, n=2).
Open the Python workspaceCheck your understanding
You must build a 10^6 x 10^6 finite-difference operator by inserting entries one at a time, then multiply it by thousands of vectors. What is the right sequence?
- Collect the entries in COO or LIL, convert once to CSR, then multiply
- Create an empty CSR matrix and assign each entry directly into it
- Build the operator dense with NumPy, then call sparse.csr_matrix on it
- Build in CSR for the insertions, then convert to COO for the multiplications
Show answer
COO and LIL make incremental insertion cheap, and a single conversion gives CSR the contiguous layout its matrix-vector kernel needs. Assigning into CSR directly looks natural but each insertion shifts the data, indices and indptr arrays, so the build cost becomes quadratic; a dense 10^6 x 10^6 array would need 8 terabytes, and COO has no row structure to make repeated products fast.