Dense, Sparse, and Matrix-Free Operators

Dense storage retains every matrix entry; sparse storage records selected entries and their positions; a matrix-free operator computes its action without assembling the full matrix. The best choice depends on what a solver needs to do. A sparse Hamiltonian can still lead to dense factors, dense eigenvectors, or an unaffordable Krylov basis.

Required background. Array representations fixes the basis, dtype, and operator action.

Storage and arithmetic costs

For an $N\times N$ matrix with $z$ stored entries, dense storage costs $O(N^2)$, while a compressed sparse row representation stores $z$ values, $z$ column indices, and $N+1$ row pointers. Matrix-vector multiplication costs $O(N^2)$ or $O(z)$ respectively. Constants matter: sparse multiplication must load indices and often accesses vector elements irregularly, whereas dense kernels reuse data efficiently.

For complex binary64 values, the value array alone uses $16N^2$ bytes in a dense matrix. At $N=10^5$, that is $160$ GB in decimal units. A tridiagonal matrix with approximately $3N$ entries instead needs about $4.8$ MB of values, plus indices and row pointers. Storing its three diagonals directly is more economical still.

A sparse format is not an instruction to discard small entries. Structural zeros come from the model or discretization; thresholding nonzero entries introduces a new approximation. If $\Delta H$ is the discarded Hermitian part, eigenvalue shifts are bounded in absolute value by $\|\Delta H\|_2$, but a componentwise threshold alone does not supply this norm bound independently of dimension.

A finite-difference operator without a matrix

For a grid function with interior values $u_j$ and zero Dirichlet endpoint values, define

$$(Hu)_j=-t u_{j-1}+(2t+V_j)u_j-tu_{j+1}, \qquad t=\frac{\hbar^2}{2mh^2},$$

with $u_0=u_{N+1}=0$ and $V_j=V(x_j)$. Each output uses only the input at that grid point and its two neighbors. Three shifted vector operations apply $H$ in $O(N)$ work without storing indices. This is a matrix-free stencil. Assemble the same matrix for a small $N$ and compare its action on a complex vector before using the stencil in a large sparse eigensolver.

For $N=3$, $V=0$, and $u=(1,i,-1)^T$, the result is

$$Hu=t(2-i,\,2i,\,-2-i)^T.$$

This hand calculation catches periodic wraparound introduced accidentally by a cyclic array shift. Boundary entries must be handled explicitly.

Solver costs beyond the Hamiltonian

If only a few extremal eigenpairs are needed, Krylov methods may avoid dense diagonalization. Yet storing $m$ trial vectors costs $O(Nm)$ and orthogonalization costs can become substantial. Asking for all eigenvectors ultimately requires $O(N^2)$ output storage even when $H$ is sparse.

Interior eigenvalues often motivate shift-invert actions $(H-\sigma I)^{-1}x$. A sparse factorization can create many additional nonzero entries, called fill-in. Its memory demand depends on ordering and graph structure, not just the original nonzero count. Alternatively, iterative inner solves add another tolerance and preconditioning problem.

The exponential $e^{-iHt/\hbar}$ is generally dense even if $H$ is sparse. To evolve one vector, use an exponential action rather than constructing the full propagator. Similarly, a sparse pure-state vector can quickly become dense under time evolution; sparsity is a property of both the representation and the state at the relevant time.

Choosing an implementation

Use dense storage for small problems, broad spectral information, and reference calculations where simplicity aids verification. Use structured diagonals for fixed-bandwidth stencils, general sparse formats for irregular connectivity, and matrix-free actions when structure makes assembly wasteful. Benchmark the whole calculation, including assembly, format conversion, factorization, and output.

Coordinate-format assembly permits repeated contributions to an entry. Ensure duplicates are summed before reasoning about nonzero counts or comparing a matrix with its adjoint. A matrix-free implementation must expose the correct adjoint action if the solver requests it; for non-Hermitian problems, reusing the forward action is incorrect.

Changing storage should preserve the mathematical problem. Compare actions and residuals on small cases before comparing timings. Faster execution of a differently truncated Hamiltonian is not a storage-only improvement.

Exercises

A calculation stores 80 complex vectors of length $10^6$. How much memory do those vectors require, ignoring allocator overhead?

Solution

$80\times10^6\times16=1.28\times10^9$ bytes, or $1.28$ GB in decimal units. This can exceed the storage of a sparse local Hamiltonian, so the Krylov dimension matters to memory planning.

Why can setting each entry smaller than $\epsilon$ to zero fail to limit the operator error to $\epsilon$?

Solution

Many discarded entries can act coherently. For example, a matrix whose every entry is $\epsilon/2$ has spectral norm $N\epsilon/2$. An entrywise threshold controls individual entries, not their collective operator action.

References

• Bai, Z., J. Demmel, J. Dongarra, A. Ruhe, and H. van der Vorst, eds., Templates for the Solution of Algebraic Eigenvalue Problems, SIAM, 2000, Sparse Matrix Storage Formats — storage schemes and matrix-vector operations.
• Saad, Y., Numerical Methods for Large Eigenvalue Problems, revised ed., SIAM, 2011 — sparse and matrix-free spectral algorithms and their memory requirements.