The single biggest specialized package in Equana is not the FFT library or the plotting layer. It is the sparse matrix package: 40 functions covering storage conversion, triplet assembly, arithmetic, reorderings, direct factorizations, eigensolvers, and eight iterative linear solvers. It is also the package almost nobody knows is there, because "sparse linear algebra in a browser tab" sounds like something that shouldn't work. It works, and the reason it works is exactly the reason sparse matrices exist at all: most of the numbers are zero, so you refuse to store or multiply them.

This post walks through how the package is built — the storage formats, the assembly model that makes finite-element codes happy, and the iterative solvers that make large systems tractable — and the memory-and-FLOPs argument for why sparse wins in a tab where you don't get 800 MB to throw at a matrix.

The zeros are the point

A dense n×n matrix of Float64 costs 8 · n² bytes no matter what is in it. At n = 10,000 that is 800 MB — enough on its own to evict everything else from a browser tab's heap. But the matrices that come out of real problems — a finite element mesh, a graph Laplacian, a discretized PDE — are almost all zeros. A node in a mesh only touches its handful of neighbors, so each row of the matrix has a constant number of nonzeros regardless of how big the mesh gets.

A sparse matrix stores only the nonzeros. If 0.1% of a 10,000² matrix is nonzero, the sparse form needs on the order of 0.2 MB instead of 800 MB — a 4,000× saving, and the gap widens as the matrix grows because the dense cost is quadratic while the sparse cost is linear in the number of nonzeros (nnz). In a browser tab, where a SharedArrayBuffer is the whole world and there is no swap to fall back on, that difference is not an optimization. It is the line between "runs" and "tab crashed".

Storage: COO in, CSR to compute

Equana keeps two views of a sparse matrix, and they exist for two different jobs.

Coordinate / triplet (COO) is the assembly format. You describe the matrix as three parallel arrays — row indices, column indices, values — in whatever order is convenient. This is what sparse(i, j, v, m, n) takes, using Equana's 1-based indexing:

# A 3x3 matrix from (row, col, value) triplets — 1-based indices rows = [1, 2, 3, 1] cols = [1, 2, 3, 3] vals = [10, 20, 30, 5] S = sparse(rows, cols, vals, 3, 3) issparse(S) # 1 nnz(S) # 4 full(S) # densify to inspect: [10 0 5; 0 20 0; 0 0 30]

Compressed Sparse Row (CSR) is the compute format. Instead of storing a row index per nonzero, CSR stores the values and their column indices contiguously, plus a short "row pointer" array of length m + 1 that says where each row begins. A row's nonzeros are then values[rowptr[i] .. rowptr[i+1]], and a sparse matrix–vector product — the inner loop of every iterative solver — becomes a tight walk over contiguous memory with no branching on structure. The transpose layout, Compressed Sparse Column (CSC), stores column pointers instead and is the natural form for column operations and many factorizations. sparse(A) converts a dense matrix straight into compressed form:

A = [1, 0, 0; 0, 2, 0; 0, 0, 3] S = sparse(A) # stored compressed; the three structural zeros are dropped nnz(S) # 3

You never touch the pointer arrays by hand. You assemble in triplets, and the package compresses.

Duplicate triplets sum — which is why FEM loves this

Here is the design decision that makes the assembly model actually usable for physics. When two triplets land on the same (i, j) position, their values add rather than overwrite.

That sounds like a small semantic footnote until you write a finite-element or finite-volume assembly loop. In those methods you build the global stiffness matrix element by element: each element contributes a small local matrix, and shared nodes between adjacent elements accumulate contributions at the same global (i, j). Summing-on-collision means you emit one triplet per local contribution and let the matrix add them up — no read-modify-write, no pre-scan to find which entries already exist. The classic discrete Laplacian is the smallest example of the pattern:

# 1D Laplacian on n interior points: tridiagonal [-1, 2, -1] n = 8 rows = zeros(3 * n - 2) cols = zeros(3 * n - 2) vals = zeros(3 * n - 2) idx = 1 for i in 1:n { rows[idx] = i; cols[idx] = i; vals[idx] = 2.0; idx = idx + 1 if i > 1 { rows[idx] = i; cols[idx] = i - 1; vals[idx] = -1.0; idx = idx + 1 } if i < n { rows[idx] = i; cols[idx] = i + 1; vals[idx] = -1.0; idx = idx + 1 } } T = sparse(rows, cols, vals, n, n) nnz(T) # 22 == 3n - 2

The package also gives you spdiags to build banded matrices directly from their diagonals, spconvert to turn [i j v] row data into a matrix, and the random generators sprand, sprandn, and sprandsym for test problems. find(S) runs the assembly in reverse, handing back the (i, j, v) triplets so you can transform the nonzeros and rebuild.

Iterative solvers: never form the inverse

Once you have a large sparse A, the question is how to solve Ax = b. A dense LU factorization is O(n³) and, worse, fills in the zeros — the factors of a sparse matrix are usually far denser than the matrix itself, which destroys the whole memory argument. For large systems the answer is iterative solvers, and this is where the package spends most of its 40 functions.

Iterative Krylov methods never form A⁻¹. They only ever multiply A by a vector, and each such product costs O(nnz) — linear in the number of stored nonzeros. Starting from a guess, they build up the solution in the subspace spanned by b, Ab, A²b, … and refine until the residual is small enough. The one to reach for depends on what you know about the matrix:

# Symmetric positive definite (e.g. the Laplacian above): conjugate gradient f = ones(n) (u, flag, relres, iters) = pcg(T, f, 1e-10, 100) flag # 0 = converged iters # iterations actually used

pcg is preconditioned conjugate gradient, the fastest choice when A is symmetric positive definite — the case that diffusion, elasticity, and Poisson problems all produce. For general nonsymmetric systems there is gmres (with an optional restart to cap memory) and bicgstab; for symmetric-indefinite systems minres; for rectangular least-squares problems lsqr; and cgs, bicg, and tfqmr round out the set. Every solver takes optional tol and maxit arguments and, when you ask for multiple return values, hands back convergence diagnostics — a flag, the relative residual, and the iteration count — so you can see whether it actually converged instead of trusting it blindly.

When iteration stalls on an ill-conditioned matrix, the package has the usual escape hatches: incomplete factorizations ilu and ichol as preconditioners, reorderings (amd, colamd, symrcm) that reduce fill-in, and direct solvers — spchol, splu, and spsolve — for moderate systems where you'd rather factorize once and reuse. On the Tauri desktop build those direct solves route to SuiteSparse; in the browser they run the WebAssembly kernels.

The FLOP argument, made concrete

Memory is half the story. The other half is arithmetic. A dense matrix–vector product is 2n² floating-point operations; the sparse version is 2·nnz. For the Laplacian above, nnz ≈ 3n, so one matrix–vector product drops from 2n² to 6n — and since a Krylov solver's whole cost is dominated by those products, the entire solve inherits the same linear-versus-quadratic gap. On a 10,000-unknown system that is the difference between hundreds of millions of operations per iteration and tens of thousands. In a tab that is running the interpreter, the WASM kernels, and your plotting all in the same heap, spending FLOPs only on nonzeros is what keeps an interactive notebook interactive.

Where the limits are

The package computes a small number of extremal eigenvalues and singular values (eigs, svds) rather than full spectra — full decompositions of a large sparse matrix are infeasible and usually not what you want. Element indexing goes through a dense view (full(S)(i, j)), so random single-element access is not the fast path; the fast path is bulk operations and solves. And nzmax equals nnz because the storage grows dynamically rather than reserving a fixed pool. None of these are surprises — they are the standard shape of a sparse library — but they are worth stating plainly so you build the way the package wants to be built: assemble in triplets, work in bulk, solve iteratively.

Forty functions is a lot of surface area for something that ships inside a browser tab. But sparse linear algebra is exactly the domain where "runs in the browser" stops being a novelty and starts being useful: the problems are big enough that the memory savings matter, structured enough that iterative solvers converge, and common enough — every PDE, every mesh, every large graph — that having them a keystroke away changes what you can do without leaving the tab.

Try it at equana.dev — the Sparse Matrices and Linear Solvers tutorials run without an account —. If you want the machinery underneath, the shared-memory story is in zero-copy WASM over SharedArrayBuffer, and the same "make a hard algorithm run anywhere" spirit shows up in radix-2 Cooley–Tukey vs Bluestein.