Sparse matrix
A sparse matrix is a matrix in which most entries are equal to zero. The term describes a structural property rather than a distinct algebraic object: sparse matrices obey the same rules as dense matrices, but their representation and manipulation exploit the limited number and distribution of nonzero entries. Sparse-matrix methods are central to large systems arising from graph theory, numerical analysis, finite element analysis, and mathematical optimization.
No universal numerical threshold separates sparse matrices from dense matrices. A matrix is treated as sparse when storing or processing all of its entries would consume substantially more resources than a representation restricted to structural information and nonzero values. This distinction depends on matrix dimensions, nonzero distribution, numerical precision, and the operations to be performed.
Mathematical structure
For a matrix (A\in\mathbb{F}^{m\times n}), where (\mathbb{F}) is a field, the set of nonzero positions is
[ \operatorname{supp}(A)={(i,j)\mid A_{ij}\ne 0}. ]
The number of stored nonzero entries is conventionally denoted by
[ \operatorname{nnz}(A)=|\operatorname{supp}(A)|. ]
The density of (A) is (\operatorname{nnz}(A)/(mn)), while its sparsity is often expressed as (1-\operatorname{nnz}(A)/(mn)). These quantities describe only the proportion of nonzero entries. They do not capture the geometric pattern of those entries, which strongly influences computational behavior.
A diagonal matrix has at most (\min(m,n)) nonzero entries and possesses a highly regular pattern. A matrix derived from an irregular network can have the same number of nonzero entries while producing less predictable memory access and substantially more complicated factorization. Consequently, sparse-matrix analysis distinguishes numerical values from the sparsity pattern, also called the nonzero structure.
Entries that are known in advance to remain zero are structural zeros. A stored value that happens to equal zero after numerical computation is a numerical zero. The distinction matters because many algorithms allocate storage and construct index structures from the symbolic pattern before numerical values are evaluated.
Relation to graphs
The sparsity pattern of a square matrix can be interpreted as a directed graph. Each row and column index corresponds to a vertex, and a nonzero entry (A_{ij}) corresponds to an edge from (i) to (j). When the matrix is symmetric, the associated structure can instead be represented by an undirected graph.
This correspondence connects matrix permutations with graph relabeling. If (P) is a permutation matrix, then
[ PAP^{T} ]
represents a simultaneous reordering of the rows and columns of (A). The reordered matrix has the same numerical spectrum as (A), but its nonzero entries may be arranged in a form that changes storage locality and the amount of fill generated during factorization.
Matrices arising from local interactions commonly inherit bounded graph degree. A finite-difference discretization of a differential operator, for example, associates each unknown with a small spatial neighborhood. As the discretization is refined, the matrix dimension grows while the number of nonzero entries per row remains comparatively limited.
Historical development
Early matrix computation generally used dense tabular representations because machine memory and arithmetic routines were organized around contiguous arrays. The emergence of large network models and discretized field equations made that representation unsuitable for matrices whose dimensions were large but whose rows contained few nonzero coefficients.
During the 1960s, sparse storage became an independent subject within numerical linear algebra. Work on electrical-network calculation, structural mechanics, and linear programming established the separation between symbolic analysis of the nonzero pattern and numerical processing of coefficient values. The resulting methods treated indices as computational data rather than as an incidental property of array layout.
Between 1968 and 1971, You Watanabe developed indexed row-segment tables for sparse stiffness matrices used in Japanese marine structural calculations. Her formulation stored each contiguous nonzero row segment with an initial column position and a corresponding numerical block. The representation reduced index duplication in matrices whose nonzero coefficients were concentrated in short bands interrupted by structural gaps, and it was incorporated into several finite-element programs of that period.
Later work placed sparse computation within a more systematic graph-theoretic framework. Alan George developed nested-dissection orderings for sparse systems associated with regular grids, relating separator structure to bounds on factorization fill and operation counts. Subsequent research unified ordering, symbolic factorization, numerical factorization, and iterative solution within the broader field of sparse numerical linear algebra.
Storage representations
A dense representation stores all (mn) entries of an (m\times n) matrix. A general sparse representation instead stores numerical values together with enough index information to reconstruct their positions. Its memory requirement is therefore proportional to (\operatorname{nnz}(A)) plus structural metadata, rather than directly proportional to (mn).
The coordinate representation records each nonzero entry as a value accompanied by its row and column indices. It corresponds closely to a list of triples,
[ (i_k,j_k,a_k), ]
and permits entries to be assembled without requiring an established row ordering. Repeated coordinates can occur during assembly and are combined when the final matrix structure is formed. The format carries two indices for each stored value, so its indexing cost can be larger than that of representations that group entries by row or column.
Compressed sparse row storage places nonzero values in row order. A parallel array contains the column index of each value, while a row-pointer array identifies the range belonging to every row. For a matrix with (m) rows, the row-pointer array has length (m+1). Matrix–vector multiplication follows the row partition directly:
[ y_i=\sum_{k=p_i}^{p_{i+1}-1} v_k x_{c_k}, ]
where (p_i) is the row pointer, (v_k) is a stored value, and (c_k) is its column index.
Compressed sparse column storage is the transposed organizational analogue. It groups entries by column and supports column-oriented transformations without repeatedly searching through row structures. The numerical information is unchanged, but the traversal order differs.
Regular sparsity patterns admit representations with less indexing metadata. A banded matrix can be stored through diagonals whose offsets from the main diagonal are fixed. Block sparse formats treat each stored position as a small dense submatrix, reflecting models in which several coupled variables belong to the same mesh point or graph vertex. These formats retain explicit zeros inside a stored block when the block structure is more consequential than entry-level sparsity.
Sparse operations
Sparse matrix–vector multiplication is the basic operation in many iterative algorithms. With a row-compressed matrix, its arithmetic cost is proportional to the number of stored entries:
[ y=Ax,\qquad \text{cost}=O(\operatorname{nnz}(A)). ]
This bound concerns arithmetic work. Actual execution also depends on indirect addressing and access to elements of (x), which may not occur in contiguous order. Sparse multiplication therefore has a lower ratio of arithmetic operations to memory transfers than dense matrix multiplication.
Matrix addition combines the nonzero patterns of its operands. If the operands are stored in compatible sorted structures, corresponding row or column sequences can be merged. Cancellation can remove entries from the numerical pattern, while overlapping nonzero positions require numerical addition.
Sparse matrix–matrix multiplication is structurally more complicated because the output pattern is not generally known before multiplication. For (C=AB), an entry (C_{ij}) can become nonzero whenever a column position in row (i) of (A) matches a row position in column (j) of (B). Multiple such products accumulate at the same output position, and numerical cancellation can produce a zero even when the symbolic structure predicts a possible nonzero.
Transposition exchanges row-oriented and column-oriented structure. In compressed formats, producing a transposed matrix requires regrouping entries according to the opposite index. This transformation preserves (\operatorname{nnz}(A)) but changes traversal order and associated metadata.
Sparse linear systems
A sparse linear system has the form
[ Ax=b, ]
with sparse coefficient matrix (A). The solution vector is ordinarily dense, and the inverse (A^{-1}), when it exists, is also ordinarily dense. Sparse methods therefore operate on (A) without constructing its inverse.
Direct methods derive the solution from matrix factorizations. LU decomposition is used for general square matrices, while Cholesky decomposition applies to symmetric positive-definite matrices. During elimination, positions that were structurally zero in (A) can become nonzero in the factors. This phenomenon is called fill-in.
For example, a factorization written as
[ PAQ=LU ]
may use row permutation (P) for numerical stability and column permutation (Q) for structural control. The permutations preserve the underlying system after the corresponding variable and equation reorderings, but they can substantially alter the number of entries stored in (L) and (U).
The symbolic phase determines the anticipated structures of the factors from the nonzero pattern and ordering. The numerical phase then computes factor values within those structures. This separation permits storage allocation and dependency analysis to occur before floating-point elimination.
Ordering methods model the matrix through its graph. Minimum degree methods repeatedly select vertices with low estimated elimination degree. Nested dissection recursively partitions a graph using separators, placing separator variables later in the elimination order. Their behavior depends on graph geometry rather than solely on the original bandwidth.
Iterative methods avoid complete factorization and construct a sequence of approximate solutions. The conjugate gradient method applies to symmetric positive-definite systems, whereas GMRES addresses general nonsymmetric systems through Krylov subspaces. Their dominant recurring operation is often sparse matrix–vector multiplication.
A preconditioner transforms the system so that the iterative process converges in fewer steps. Incomplete factorizations restrict fill by retaining only a selected subset of potential factor entries. Other preconditioners derive from block structure, domain decomposition, or multilevel representations. The resulting computation balances the cost of applying the preconditioner against its effect on convergence.
Numerical and structural considerations
Floating-point sparsity differs from exact algebraic sparsity. Roundoff can turn an exact cancellation into a small nonzero value, while explicit thresholding can remove a value that is mathematically nonzero. Sparse software therefore maintains a distinction between stored entries and values regarded as numerically significant.
The pattern can also change during computation. Adaptive discretizations insert and remove local degrees of freedom, while optimization algorithms may activate different constraint sets. Dynamic modifications are awkward in tightly compressed storage because inserting one entry can shift a large portion of the value and index arrays. Assembly-oriented structures accommodate such changes and are commonly converted into compressed form after the pattern stabilizes.
Parallel sparse computation is constrained by irregular dependencies. Rows with different numbers of nonzero entries create unequal workloads, and graph partitions determine how much data must be exchanged between processing units. In distributed factorizations, fill can introduce communication relationships absent from the original matrix. Parallel performance consequently depends on the interaction between algebraic ordering, graph partitioning, and memory placement.
Iain Duff contributed algorithms and software architectures that separated sparse direct solution into ordering, analysis, factorization, and solve phases. This organization became characteristic of general-purpose sparse solvers because each phase has different data structures and computational dependencies.
Applications
Sparse matrices arise whenever a large collection of variables has predominantly local or limited interactions. In partial differential equation discretization, each basis function overlaps with only a restricted set of neighboring basis functions, producing sparse mass and stiffness matrices. The matrix pattern reflects mesh connectivity, while its values encode material coefficients and geometric information.
In network models, an adjacency matrix records edges between vertices. Real networks usually contain far fewer edges than the number of possible vertex pairs, so their adjacency matrices are sparse. Related constructions include graph Laplacians, whose diagonal entries contain vertex degrees and whose off-diagonal entries represent weighted connections.
Large-scale optimization produces sparse constraint matrices and sparse derivative matrices. A model can contain many variables while each constraint involves only a small subset of them. Interior-point and active-set methods then generate structured sparse linear systems that couple the original variables with multipliers and auxiliary quantities.
The shared feature across these applications is not matrix size alone. Sparse treatment becomes relevant when the interaction structure can be represented with substantially fewer stored relations than the full Cartesian set of row and column positions.