SciPy Cheatsheet

Sparse Matrices

Use this SciPy reference while you build software engineering projects, review code, or refresh the syntax you reach for most.

Core Import

import scipy.sparse as sp
import numpy as np

# Modern API (SciPy 1.8+): prefer *array* over *matrix*
from scipy.sparse import (
    csr_array, csc_array, coo_array, lil_array, dok_array, bsr_array, dia_array,
    eye, identity, diags, block_diag, block_array,
    spdiags, triu, tril, hstack, vstack, kron,
    issparse, isspmatrix
)
from scipy.sparse.linalg import (
    spsolve, spsolve_triangular,
    eigs, eigsh, svds,
    cg, bicg, bicgstab, cgs, gmres, lgmres, minres, qmr, gcrotmk,
    splu, spilu, factorized,
    LinearOperator, aslinearoperator, norm
)

csr_array vs csr_matrix: From SciPy 1.8, array classes are preferred. They follow ndarray semantics (no implicit broadcasting of * as matrix multiply — use @ instead). Old csr_matrix still works but is effectively in maintenance mode.

Sparse Formats

Format	Full Name	Best For	Notes
`csr_array`	Compressed Sparse Row	Row slicing, matrix-vector multiply	Most common; fast `A @ v`
`csc_array`	Compressed Sparse Column	Column slicing, direct solvers	Fast column ops; preferred by UMFPACK
`coo_array`	Coordinate	Building / converting	Duplicate entries allowed (summed)
`lil_array`	List of Lists	Incremental construction	Slow arithmetic; use only to build
`dok_array`	Dictionary of Keys	Random element access/modification	Like a sparse dict
`bsr_array`	Block Sparse Row	Block-structured matrices (FEM)
`dia_array`	Diagonal	Diagonal/banded matrices	Very memory efficient

Construction

From Dense Array

import numpy as np
from scipy.sparse import csr_array

A = np.array([[1, 0, 2], [0, 3, 0], [4, 0, 5]])
sparse_A = csr_array(A)

print(sparse_A.shape)    # (3, 3)
print(sparse_A.nnz)      # 5 (number of stored elements)
print(sparse_A.dtype)    # dtype of stored values

From COO (Most Flexible)

from scipy.sparse import coo_array

row = np.array([0, 0, 1, 2, 2])
col = np.array([0, 2, 1, 0, 2])
data = np.array([1, 2, 3, 4, 5], dtype=float)

A = coo_array((data, (row, col)), shape=(3, 3))
A_csr = A.tocsr()

From CSR Directly

from scipy.sparse import csr_array

# CSR format: data, (row indices sorted, column indices), pointers into arrays
# Easier: build from COO then convert

# Or using indptr, indices, data
indptr  = np.array([0, 2, 3, 5])     # row i has cols in indices[indptr[i]:indptr[i+1]]
indices = np.array([0, 2, 1, 0, 2])  # column indices
data    = np.array([1, 2, 3, 4, 5], dtype=float)
A = csr_array((data, indices, indptr), shape=(3, 3))

Incremental Construction

from scipy.sparse import lil_array, dok_array

# lil_array: fastest for row-wise construction
A = lil_array((1000, 1000))
A[0, 5] = 1.0
A[10, 200] = 3.0
A[500:510, 0] = np.arange(10)     # slice assignment
A_csr = A.tocsr()                  # convert before arithmetic

# dok_array: fast for random access
A = dok_array((1000, 1000))
A[5, 10] = 2.0
A[(5, 15)] = 3.0

Special Matrices

from scipy.sparse import eye_array, diags_array, block_diag, kron

# Identity (sparse ARRAY constructors, SciPy 1.12+)
I = eye_array(5, dtype=float, format='csr')

# Diagonal
d = diags_array([1.0, 2.0, 3.0, 4.0, 5.0])       # main diagonal
d = diags_array([[4., 5.], [1., 2., 3.], [6., 7.]],
                offsets=[-1, 0, 1])               # tridiagonal 3x3

# Block diagonal
from scipy.sparse import block_diag as sp_block_diag
A = sp_block_diag([A1, A2, A3], format='csr')

# Kronecker product
K = kron(A1, A2, format='csr')

Legacy: eye, diags, and spdiags return sparse matrix objects. On SciPy 1.12+ use eye_array / diags_array, which return sparse arrays and match the rest of the array API.

Format Conversion

csr = A.tocsr()
csc = A.tocsc()
coo = A.tocoo()
lil = A.tolil()
dok = A.todok()
bsr = A.tobsr()
dia = A.todia()

# To/from dense
dense = A.toarray()                  # → numpy ndarray
sparse = csr_array(dense_array)

# Check format
print(A.format)   # 'csr', 'csc', 'coo', etc.

Array Properties & Access

A.shape       # (m, n)
A.nnz         # number of stored elements
A.ndim        # 2
A.dtype
A.data        # stored values
A.indices     # column indices (CSR) or row indices (CSC)
A.indptr      # row pointers (CSR) or column pointers (CSC)

# Element access (slow for CSR/CSC on random access)
val = A[1, 2]
row_slice = A[0, :]    # returns sparse array
col_slice = A[:, 1]

# Nonzero indices
rows, cols = A.nonzero()

# Density
density = A.nnz / (A.shape[0] * A.shape[1])

Arithmetic Operations

# Matrix-vector multiply (key operation)
x = np.random.rand(A.shape[1])
y = A @ x              # fast sparse matvec
y = A.dot(x)           # equivalent

# Matrix-matrix multiply
C = A @ B
C = A.dot(B)

# Element-wise
C = A + B
C = A - B
C = A.multiply(B)      # element-wise multiply (Hadamard)
C = A.power(2)         # element-wise square

# Scalar
C = 2.0 * A
C = A / 3.0

# Transpose
At = A.T
At = A.transpose()

# Conjugate transpose
Ah = A.conj().T

# Diagonal
d = A.diagonal()
d = A.diagonal(k=1)     # superdiagonal

Stacking & Combining

from scipy.sparse import hstack, vstack, block_array

C = hstack([A, B])                    # horizontal concatenation
C = hstack([A, B], format='csr')
C = vstack([A, B])                    # vertical concatenation
C = vstack([A, B, C], format='csc')

# block_array (SciPy 1.11+, ndarray-style)
C = block_array([[A, B], [None, D]], format='csr')
# None → zero block of appropriate size

Sparse Solvers

Direct Solver

from scipy.sparse.linalg import spsolve, spsolve_triangular, splu, factorized

A = csr_array(...)    # must be square
b = np.random.rand(n)

# One-shot solve
x = spsolve(A, b)              # uses UMFPACK/SuperLU
x = spsolve(A, B)              # B is dense matrix → multiple RHS

# Triangular solve
x = spsolve_triangular(A, b, lower=True)

# Factorize once, solve many times (LU)
lu = splu(A.tocsc())           # SuperLU factorization (needs CSC)
x1 = lu.solve(b1)
x2 = lu.solve(b2)
x_T = lu.solve(b, trans='T')  # solve A.T @ x = b

# factorized: returns a callable
solve = factorized(A.tocsc())
x = solve(b)

# ILU preconditioner
ilu = spilu(A.tocsc(), drop_tol=1e-4, fill_factor=10)
M = sp.linalg.LinearOperator(A.shape, ilu.solve)   # use as preconditioner

Iterative Solvers

from scipy.sparse.linalg import cg, bicgstab, gmres, minres

x0 = np.zeros(n)

# Conjugate Gradient (symmetric positive definite)
x, info = cg(A, b, x0=x0, tol=1e-8, maxiter=1000, M=M_precond)

# BiCGSTAB (general, non-symmetric)
x, info = bicgstab(A, b, tol=1e-8, maxiter=1000)

# GMRES (general, non-symmetric, more memory)
x, info = gmres(A, b, tol=1e-8, restart=30, maxiter=100)

# MINRES (symmetric, possibly indefinite)
x, info = minres(A, b, tol=1e-8)

# info: 0 = success, >0 = convergence not achieved, <0 = illegal input
if info == 0:
    print("Converged")
else:
    print(f"Did not converge (info={info})")

# Callback to track convergence
residuals = []
x, _ = gmres(A, b, callback=lambda r: residuals.append(r))

Solver	Matrix Type	Notes
`cg`	SPD	Fastest for well-conditioned SPD
`minres`	Symmetric (indefinite ok)
`bicg`	General	Less stable than bicgstab
`bicgstab`	General	More stable than bicg
`cgs`	General	Squared CG
`gmres`	General	Robust; stores Krylov basis (memory)
`lgmres`	General	Less memory than gmres
`gcrotmk`	General	Flexible restart, memory efficient
`qmr`	General	Quasi-minimal residual

Eigenvalue Solvers (ARPACK)

from scipy.sparse.linalg import eigs, eigsh

# eigs: general (non-symmetric) eigenvalues
vals, vecs = eigs(A, k=6, which='LM')   # 6 largest magnitude eigenvalues
vals, vecs = eigs(A, k=6, which='SM')   # 6 smallest magnitude
vals, vecs = eigs(A, k=6, which='LR')   # largest real part

# eigsh: symmetric (real symmetric / Hermitian)
vals, vecs = eigsh(A, k=6, which='LM')
vals, vecs = eigsh(A, k=6, which='SM', sigma=0.0)  # near sigma=0 (shift-invert)

# Shift-invert (for interior eigenvalues)
vals, vecs = eigsh(A, k=6, sigma=1.0, which='LM')

# which options: 'LM', 'SM', 'LR', 'SR', 'LI', 'SI' (complex only)

SVD

from scipy.sparse.linalg import svds

# k largest singular values/vectors
U, s, Vh = svds(A, k=5)
# Note: results are in *ascending* order — reverse if needed
U, s, Vh = svds(A, k=5)
idx = np.argsort(s)[::-1]
U, s, Vh = U[:, idx], s[idx], Vh[idx, :]

Sparse Norm

from scipy.sparse.linalg import norm

n1 = norm(A, ord=1)      # max column sum
ni = norm(A, ord=np.inf) # max row sum
nf = norm(A, ord='fro')  # Frobenius

Linear Operators

from scipy.sparse.linalg import LinearOperator

# Define a matrix-free operator (e.g., FFT-based)
def matvec(x):
    return np.fft.ifft(np.fft.fft(x) * freq_domain_multiplier)

def rmatvec(x):
    return np.conj(matvec(np.conj(x)))

op = LinearOperator(shape=(n, n), matvec=matvec, rmatvec=rmatvec, dtype=complex)

# Use with iterative solvers
x, info = gmres(op, b)
vals, vecs = eigs(op, k=5)

Common Gotchas

*`csr_array csr_array is element-wise** in new array API. Use @ for matrix multiplication. Old csr_matrix * csr_matrix` was matrix multiply — a major semantic difference.

Convert to CSC for direct solvers. spsolve and splu internally convert to CSC if needed, which silently copies the data. Pass CSC from the start to avoid the copy.

lil_array is slow for arithmetic. Use it only to build the matrix incrementally, then convert with .tocsr() before any computation.

eigs/eigsh cannot return all eigenvalues. Use k < n - 1. For all eigenvalues of a sparse matrix, convert to dense and use scipy.linalg.eigh.

ARPACK which='SM' is numerically unstable without shift-invert. Use sigma=0.0 (shift-invert mode) to find eigenvalues near zero reliably.

Iterative solver tol is relative to the initial residual ||b||. If the initial guess is bad and b is large, you may need a lower tol than expected.