Solving Dense Matrices Fast in the Browser — Adaptive Frequency Sampling and Mixed Precision on WebGPU

Almost all the work in the electromagnetic simulator goes into one task: solving a complex dense matrix with N unknowns at every frequency. The effort grows as N³, so doubling N costs eight times as much and multiplying it by ten costs a thousand times as much.

This article is about what I did to make that fast entirely inside the browser. The result: the N = 2157 structure used to reproduce a paper (21 frequency points) went from about 29 s to 2.1 s (Apple M1).

Measure first: the ceiling of JavaScript

First I measured how far plain JavaScript could go.

The matrix Z is complex symmetric (equal to its transpose; not Hermitian), so an LDLᵀ factorisation that uses the symmetry needs half the work of LU. If a small pivot appears, it falls back automatically to LU with partial pivoting.

Blocking failed because the cause was not the one I expected. I had assumed memory bandwidth was the limit, but the measured 3.7 GFLOP/s corresponds to about 15 GB/s of memory traffic, well short of the M1's bandwidth (about 68 GB/s). The limit was the arithmetic itself in scalar JavaScript, which cannot use vector instructions. That was the end of micro-optimisation inside JavaScript.

I also confirmed that Python with numpy (on Apple's Accelerate) runs at 103–123 GFLOP/s, 12–15 times faster, but passed on it because it would break the "just open it in a browser" value of the tool.

Memory was a problem too. At N = 8000 one matrix is 512 MB, and with working copies the peak reached about 2.9 GB, enough to risk crashing the tab. Using the symmetry of L, Ψ and Z to store only the packed lower triangle, and factorising in place, brought the peak down to about 1 GB.

Fewer frequency points: adaptive frequency sampling

Next I reduced the number of solves itself.

As the previous article explained, the matrix has the form

Z(ω) = jωL − jΨ/ω + Zs(ω)·G

where L, Ψ and G do not depend on frequency. So the idea is to solve exactly at only a few frequencies, and shrink the problem onto the space spanned by those solution vectors (a reduced-order model — the same idea as Momentum's Adaptive Frequency Sampling).

1. Solve exactly at a few frequencies, and orthonormalise the real and imaginary parts of the solutions into a real basis V
2. Form the small matrix Z_r = jω·VᵀLV − jVᵀΨV/ω + Zs·VᵀGV (q×q, with q around 20–60). It can be solved at any frequency in an instant
3. Estimate the residual ‖Z(ω)·V·y − b‖ at every candidate frequency, and solve exactly next where it is largest
4. Stop when the residual falls below 1e-4

Because V is real, the reduced matrix is symmetric just like the original, so reciprocity and passivity are preserved. At the frequencies solved exactly, it matches the exact solution perfectly. Estimating the residual also needs only small q×q computations.

For the default presets, a 41-point sweep needs 8 exact solves, and differs from solving every point exactly by 3e-13. Since adding output points barely changes the run time, the point limit went up to 1001. The small triangles along the frequency axis of the S-parameter plot mark the frequencies solved exactly.

Double-precision answers from a GPU without double precision

Once the exact solves are down to around ten, the single factorisation (about 3 GFLOP/s) becomes the limit. Trying WebGPU, which gives the browser access to the GPU, a bare matrix multiplication ran at 550–630 GFLOP/s. But WebGPU has no double precision (64 bit). A single-precision (32-bit) factorisation does not carry enough significant digits.

So I used mixed-precision iterative refinement:

1. Factorise the matrix in single precision on the GPU (LU in blocks of 32, with matrix products in 64×64 tiles)
2. Solve the system with that factorisation to get an approximate solution x
3. Compute the residual r = b − Z·x in double precision
4. Solve for a correction with the single-precision factors and add it to x; go back to step 2

After 4–7 iterations the answer agrees with a double-precision direct solve on the CPU to about 1e-11 relative. Factorising N = 2157 takes 74 ms on the GPU (about 4 s on the CPU), and one frequency's solve including refinement takes 0.16–0.25 s. The factorisation does not pivot, so if refinement fails to converge at some frequency, that frequency alone switches automatically to a direct solve on the CPU.

The shader compiler was erasing the error

It is tempting to compute the double-precision residual in step 3 on the GPU as well. There is a technique (double-float) that combines two single-precision numbers to get about 48 bits of precision, and at its heart is an "error-free addition" like this:

s = a + b
e = b − (s − a)    // the part lost to rounding

But on Apple's GPU the error term e came out always 0. It turned out that the shader compiler, to go faster, was simplifying (a + b) − a to b mathematically — a transformation that does not hold in floating point, and that breaks the very premise of the technique.

The fix was to make the compiler treat each intermediate result as an integer once. Taking an exclusive OR (XOR) with a 0 passed in from outside and converting back to floating point cuts the value out of the expression tree as far as the compiler is concerned, so it can no longer simplify it. The value itself does not change. With this in place, the results agree with double precision on the CPU to 1e-14 on random vectors.

There is no guarantee that the same will not happen on some other GPU, so the first computation after start-up checks the GPU against the CPU with random vectors, and if they do not agree, the residual is left to the CPU.

Moving the precomputation to the GPU too

With solving fast, assembling the frequency-independent matrices L and Ψ started to stand out (about 7 s at N = 6614). That moved to the GPU as well:

• For cells far apart (about 98–99 % of all pairs), the cell integrals reduce to a point approximation and the computation becomes a simple series. The GPU computes it for every pair
• Pairs of nearby cells, and the diagonal, are computed by the CPU with the original exact formulas and written over the GPU's values. Nearby pairs are found by dividing the space into buckets, and pairs with the same geometry reuse each other's results

The GPU path differs from the original double-precision computation by about 1e-15 relative. The precomputation became about 15 times faster.

Results

The N = 2157 structure used to reproduce the paper, 21 frequency points (Apple M1):

Once time is no longer the constraint, what remains is memory. The GPU side needs about 20·N² bytes and the browser side about 8·N². Running out can crash the tab rather than raise an error, so the limit on unknowns is set from the environment: up to 12000 when a GPU is available, 8000 otherwise. A solid conductor with N = 11402 finished a 41-point sweep in 123 s (14 exact solves).

Lessons

• Measure rather than guess. I tried blocking on the assumption that memory was the limit; the cause was the arithmetic itself
• Doing fewer solves through the algorithm often pays more than doing one solve faster. Adaptive frequency sampling pays more the more points you ask for
• Always pair reduced-precision computation with a way to check it. Mixed-precision iterative refinement has a built-in check in the residual, and a self-test catches unexpected traps such as a compiler optimisation

The next article uses this computation to try to reproduce structures from a paper. → Solving pixelated circuits