Sparse Matrices

A sparse matrix is a matrix whose useful structure is carried by relatively few nonzero entries. Instead of storing every one of its $n^2$ entries, a sparse algorithm stores only the nonzeros and the indexing information needed to find them.

Sparse matrices matter in quantum mechanics because many numerical Hamiltonians are local. A finite-difference kinetic energy couples nearby grid points. A local spin-chain Hamiltonian changes only a few sites at a time. A tight-binding model connects neighboring orbitals. These systems may have enormous Hilbert spaces, while each basis vector couples to only a small fraction of all other basis vectors.

The word sparse is not only a visual description. It is a computational contract: algorithms should use the number of nonzero entries, not the number of possible entries, as the main cost scale.

Definition

For an $n\times n$ matrix $A$ , define

\operatorname{nnz}(A) = \#\{(i,j):A_{ij}\ne0\}.

A dense algorithm typically stores and manipulates all $n^2$ entries. A sparse algorithm is useful when

\operatorname{nnz}(A) \ll n^2.

The precise meaning depends on the family of problems. A tridiagonal $n\times n$ matrix has only $3n-2$ nonzero entries, so it remains sparse as $n$ grows. A matrix with $0.1n^2$ nonzeros may look mostly empty on a small plot, but it is still quadratic in storage and usually behaves more like a dense matrix at large $n$ .

Why Quantum Hamiltonians Are Often Sparse

Sparsity usually comes from locality or selection rules:

a one-dimensional finite-difference Hamiltonian couples each grid point only to nearby grid points;
a local potential is diagonal in a position grid;
nearest-neighbor tight-binding models connect only adjacent sites;
few-body spin interactions flip or compare a small number of spins;
angular-momentum selection rules can forbid many matrix elements;
symmetry block decomposition splits one large matrix into smaller independent blocks.

Sparsity is representation-dependent. The same operator can be sparse in one basis and dense in another. A kinetic-energy operator is diagonal in momentum space but a local potential is diagonal in position space. A Fourier method may avoid a sparse matrix entirely and use transforms instead.

For the continuum-to-finite step, see Discretization. For the standard grid derivative construction, see Finite Difference Methods.

Dense Versus Sparse Storage

A dense matrix stores every entry:

A = \begin{pmatrix} A_{11} & A_{12} & \cdots & A_{1n}\\ A_{21} & A_{22} & \cdots & A_{2n}\\ \vdots & \vdots & \ddots & \vdots\\ A_{n1} & A_{n2} & \cdots & A_{nn} \end{pmatrix}.

This is simple and efficient for small matrices, but memory grows like $n^2$ . For $n=10^6$ , dense storage is not a practical representation.

Sparse storage records the nonzero values and enough integer indices to reconstruct their positions. The cost is roughly proportional to $\operatorname{nnz}(A)$ , with an additional indexing overhead. That overhead is worthwhile only when the matrix is genuinely sparse and the algorithm respects the sparse structure.

Common Sparse Formats

Several storage formats are standard. The best one depends on how the matrix is assembled and used.

Format	Idea	Typical use
COO	store triples $(i,j,A_{ij})$	easy assembly from local contributions
CSR	compressed sparse row storage	fast row-wise matrix-vector products
CSC	compressed sparse column storage	column access and some factorizations
DIA	store diagonals by offset	banded finite-difference matrices
BSR	store dense blocks sparsely	spin, orbital, or multi-component systems

COO is convenient while building a matrix because one can append local terms. CSR is usually better once the matrix is fixed, because matrix-vector multiplication reads one row at a time. CSC is the column-oriented analogue and is useful in some direct solvers.

Block sparse formats are common when each grid point or lattice site carries internal degrees of freedom. For example, a spinor wavefunction with two components per grid point often gives a sparse matrix whose nonzero entries are small dense blocks.

Finite-Difference Example

For a one-dimensional grid with Dirichlet boundaries, the standard second-difference kinetic-energy operator has a tridiagonal matrix. Up to constants and the potential term, the Hamiltonian has the form

H = \begin{pmatrix} d_1 & t & 0 & \cdots & 0\\ t & d_2 & t & \cdots & 0\\ 0 & t & d_3 & \ddots & \vdots\\ \vdots & \vdots & \ddots & \ddots & t\\ 0 & 0 & \cdots & t & d_N \end{pmatrix}.

The diagonal entries include the potential and the diagonal kinetic contribution. The off-diagonal entries come from nearest-neighbor coupling. The number of nonzero entries is

\operatorname{nnz}(H) = N+2(N-1) = 3N-2.

This is why a grid with millions of points can still be meaningful for certain one-dimensional or structured problems, even though a dense matrix of the same dimension would be impossible to store.

Tensor-Product Grids

In more than one dimension, separable grid operators often have Kronecker-sum structure. For a two-dimensional Cartesian grid, a schematic Hamiltonian is

H = T_x\otimes I_y + I_x\otimes T_y + V.

Here $T_x$ and $T_y$ are one-dimensional second-difference matrices, and $V$ is diagonal if the potential is evaluated pointwise. If there are $n$ points in each coordinate and $M=n^2$ total grid points, each row still couples to only a small stencil of nearby points. The storage grows like a constant times $M$ , not like $M^2$ .

For $d$ dimensions with a nearest-neighbor stencil, the number of nonzero entries is typically of order $dM$ , where $M$ is the total number of grid points. This scaling is still large because $M=n^d$ grows rapidly with dimension, but sparsity prevents an additional square in the storage.

The tensor-product notation itself is reviewed in Tensor Products.

Matrix-Vector Products

The central sparse operation is not usually forming all eigenvectors. It is applying the matrix to a vector:

y=Ax.

In components, a sparse matrix-vector product reads

y_i = \sum_{j:A_{ij}\ne0} A_{ij}x_j.

The cost is proportional to $\operatorname{nnz}(A)$ rather than $n^2$ . This operation is the engine behind Sparse Eigensolvers, iterative linear solvers, many Time-Stepping Methods, and matrix-free operator implementations.

For a finite-difference Hamiltonian, the same multiplication can often be written directly as a stencil action. In one dimension,

(H\psi)_j = t\psi_{j-1} + d_j\psi_j + t\psi_{j+1}.

This formula may be faster and clearer than explicitly storing the matrix, especially when boundary conditions are simple.

Matrix-Free Operators

Some large quantum calculations never store the sparse matrix at all. They provide a routine that computes $Ax$ for any input vector $x$ . This is called a matrix-free representation.

Matrix-free methods are natural when:

the operator is a simple stencil;
the Hamiltonian is assembled from local terms on demand;
the basis is too large for even sparse storage;
the action uses fast transforms rather than matrix entries;
only a few eigenvalues or time steps are needed.

The tradeoff is that matrix-free operators are harder to inspect. One should still test Hermiticity, norms, boundary behavior, and benchmark cases.

Hermiticity and Weighted Inner Products

For closed-system quantum mechanics, the finite Hamiltonian should represent a self-adjoint operator. In an ordinary Euclidean finite-dimensional inner product, this means

H^\dagger=H.

Sparse assembly can accidentally break Hermiticity. Common causes include adding an off-diagonal entry $H_{ij}$ without adding the conjugate entry $H_{ji}^\ast$ , applying boundary rows asymmetrically, or mixing indexing conventions.

On nonuniform grids or quadrature-based discretizations, the correct inner product may be weighted:

\langle u,v\rangle_W = u^\dagger Wv.

Then self-adjointness is the condition

H^\dagger W = WH.

This is not the same as ordinary matrix symmetry unless $W$ is proportional to the identity. The weighted-inner-product issue is introduced in Discretization, and quadrature weights are treated in Numerical Quadrature.

Sparse Does Not Mean Easy

Sparse storage solves the memory problem for matrix entries, but it does not remove all numerical difficulty.

Sparse direct factorizations can create fill-in: entries that were zero in $A$ become nonzero in triangular factors. The inverse of a sparse matrix is usually dense. Therefore, computing $A^{-1}$ explicitly is almost never the right way to solve a large quantum problem.

Conditioning is a separate issue. A sparse matrix can be ill conditioned, and a dense matrix can be well conditioned. For example, refining a finite-difference grid often makes derivative matrices larger and more ill conditioned, even though they remain sparse. The sensitivity language is developed in Conditioning and Stability.

Sparse Eigenvalue Problems

Dense diagonalization is appropriate when all eigenvalues and eigenvectors of a modest matrix are needed. Sparse Hamiltonians are different. One usually asks for a few extremal or interior eigenvalues, such as the ground state and low-lying excitations.

Those computations are built from repeated matrix-vector products, orthogonalization, residual checks, and convergence tests. The important point on this page is that sparse storage makes the matrix-vector products possible. The eigensolver algorithms themselves are treated in Sparse Eigensolvers.

Do not expect sparse methods to compute the entire spectrum of a huge Hamiltonian cheaply. Sparsity helps selected computations; it does not defeat Hilbert-space dimension.

Practical Assembly Checklist

When assembling a sparse Hamiltonian, record:

the basis or grid ordering;
the indexing convention;
the boundary rows or boundary terms;
the storage format used for assembly and for computation;
whether duplicate COO entries are summed;
the inner product in which Hermiticity is tested;
the number of degrees of freedom and $\operatorname{nnz}(H)$ ;
the cost of one matrix-vector product;
the benchmark problem used to test the implementation.

These details are not low-level clutter. They determine whether a numerical result is reproducible.

Common Mistakes

Building a dense matrix first and converting it to sparse afterward.
Assuming a matrix is sparse because a small plot looks mostly empty.
Confusing sparse storage with good conditioning.
Forgetting the conjugate partner of an off-diagonal complex entry.
Treating $H^\dagger=H$ as the right test when the discretized inner product is weighted.
Computing an explicit inverse instead of using solves or matrix-vector products.
Ignoring fill-in during sparse factorization.
Expecting a sparse eigensolver to return all eigenpairs efficiently.
Changing basis and assuming sparsity survives.
Reporting a large sparse calculation without $\operatorname{nnz}$ , boundary conventions, or residual checks.

Cross-Links

References

T. A. Davis, Direct Methods for Sparse Linear Systems, SIAM, 2006.
Y. Saad, Iterative Methods for Sparse Linear Systems, 2nd ed., SIAM, 2003.
Y. Saad, Numerical Methods for Large Eigenvalue Problems, revised ed., SIAM, 2011.
G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed., Johns Hopkins University Press, 2013.
L. N. Trefethen and D. Bau, Numerical Linear Algebra, SIAM, 1997.
J. M. Thijssen, Computational Physics, 2nd ed., Cambridge University Press, 2007.

Exercises

A tridiagonal $N\times N$ matrix has nonzero entries only on the main diagonal and the two neighboring diagonals. Count its nonzero entries.

Solution

The main diagonal contributes $N$ entries. The upper diagonal has $N-1$ entries, and the lower diagonal has $N-1$ entries. Therefore

\operatorname{nnz} = N+(N-1)+(N-1) = 3N-2.

On a square $n\times n$ grid, suppose each interior point couples to itself and its four nearest neighbors. How does the number of nonzero matrix entries scale with the total number $M=n^2$ of grid points?

Solution

Each interior row has about five nonzero entries, while boundary rows have fewer. Thus

\operatorname{nnz}(H) \sim 5M

up to boundary corrections. The scaling is linear in the number of grid points, not quadratic in $M$ .

Why can a sparse finite-difference Hamiltonian fail to be Hermitian even when the intended differential operator is self-adjoint?

Solution

The discrete assembly may treat boundary rows asymmetrically, omit conjugate off-diagonal entries, use inconsistent indexing, or test symmetry in the wrong inner product. Self-adjointness is a property of the discretized operator with its chosen domain and inner product, not only of the formal differential expression.

Explain why explicitly forming $A^{-1}$ is usually a bad strategy for a large sparse linear problem.

Solution

Even when $A$ is sparse, its inverse is usually dense. Forming $A^{-1}$ can destroy the memory advantage of sparsity and amplify numerical errors. Large sparse problems are usually solved by sparse factorizations, iterative methods, or matrix-vector-product based algorithms, depending on the structure and conditioning.