6292f12A benchmark chapter, and offline benchmarks to compare it againstJeremy Magland 1Every chapter so far has posed a problem shaped to show one method at its best. This closing chapter does the opposite. It fixes one problem that is unfair to everybody and runs the tutorial's four families of methods on it, so that the accuracy-against-cost question the introduction promised can finally be answered on level ground: same domain, same data, same error points, same clock.
3The domain is the *Pac-Man*,
5$$\Omega = \{(r,\theta) : 0 < r < 1,\ 0 < \theta < 3\pi/2\},$$
7the unit disk with the quadrant $x > 0$, $y < 0$ removed. Its boundary is a smooth arc and two straight edges; the edges meet the arc in two right-angle corners, and each other in a re-entrant corner of interior angle $3\pi/2$ at the origin. Both difficulties of Part "Beyond the square" are therefore present at once: the boundary cuts through grid cells, as in Chapter 10, and the re-entrant corner carries the $r^{2/3}$ singularity of Chapter 9. The corner has the same angle as the L-shape's, at the same origin, so the exact solutions of Chapter 9 (`lp_lexact`) work here unchanged: the singular function $r^{2/3}\sin(2\theta/3)$, the smooth $\log|z - z_0|$, and their sum, which is the generic case. As always, the data is the exact solution's boundary trace, and the error is known everywhere.
9<!--fig:pacman-->
11### Four methods, four handicaps
13Each method is the one you have already met, pushed off its home ground.
15- **Shortley-Weller + SOR.** The five-point scheme, with the stencil legs shortened wherever a neighbour falls outside the arc, so the boundary condition is imposed on the arc itself rather than on a staircase. The price is symmetry: the shortened stencil is not a symmetric matrix, conjugate gradients no longer applies, and the solver here is SOR. On smooth data the scheme is cleanly second order, which is Shortley and Weller's achievement; on generic data the corner takes that away, exactly as in Chapter 9.
16- **The harmonic basis.** The method of particular solutions, with the Chapter 9 basis: harmonic polynomials plus the corner functions $r^{2k/3}\sin(2k\theta/3)$, fitted by least squares at boundary points clustered towards the corner. Nothing about a curved arc troubles it, since the basis never sees the boundary's shape, only samples of the data on it.
17- **The double layer on graded panels.** The Nyström method of Chapter 10 loses its trapezoid rule here, because the rule's spectral accuracy needed a smooth periodic integrand and the corners break both the geometry and the density: at each corner the density itself is singular, near the re-entrant one like $r^{-1/3}$. The replacement is composite 12-point Gauss-Legendre on panels refined dyadically towards the three corners, to a depth $d$ that is the method's resolution parameter.
18- **Chebyshev on the polar rectangle.** In polar coordinates the Pac-Man is the rectangle $(r, \theta) \in (0,1) \times (0, 3\pi/2)$, Laplace's equation becomes $r^2 u_{rr} + r u_r + u_{\theta\theta} = 0$, and the matrix decomposition of Chapter 7 applies with two different one-dimensional operators. The whole side $r = 0$ of the rectangle is the corner, a single boundary point with a known value, so the coordinate degeneracy costs nothing to impose. What it costs is convergence: the $r^{2/3}$ singularity now sits on an edge of the expansion domain, and the geometric convergence of Chapters 4 and 7 drops to an algebraic rate.
20<!--lab:field-->
22### The benchmark
24The sweep below runs all four methods over their own resolution ladders and measures the maximum error at a fixed lattice of interior points, kept a little away from the boundary (the reason is the close-evaluation problem of Chapter 10, which would otherwise charge the double layer for a difficulty every method of this family has and this chapter is not about) and reaching within $0.125$ of the corner. Timings are for the solve only, warmed up, as everywhere in this tutorial.
26<!--lab:bench-->
28On the generic data the four curves tell one story with four morals. The difference scheme converges at the corner-limited order $4/3$ and buys its accuracy with by far the most unknowns: at $N = 128$ it holds about $38\,000$ of them and stands at $4 \times 10^{-4}$. The polar Chebyshev method does somewhat better than second order in this range, but it has lost its exponential rate to the corner; on the `smooth` data, where nothing is singular, it reaches $10^{-10}$ with $40$ points per direction and is the second-best method on the board, which is a useful reminder that the ranking is a property of the problem, not of the methods. The graded double layer converges geometrically in the depth $d$, but with the modest ratio $2^{2/3} \approx 1.6$ per unit of depth, which is the corner exponent speaking: the leftover error is what the unresolved $r^{-1/3}$ tail of the density costs, about $2^{-2d/3}$. And the harmonic basis, whose every element satisfies the equation and whose two singular families are both in the span, is exact to roundoff with 97 unknowns in a few milliseconds, ten to twelve orders of magnitude ahead of everything else at equal cost. The conditions attached to that verdict are the ones stated at the end of Chapter 9, and they have not weakened: the equation must be homogeneous and linear with known solutions, and the singularities must be known in advance.
30Note that the first plot cannot rank the methods, because each refines a different parameter; only the plot against time can. That is the plot this tutorial has been building towards since the introduction.
32### Benchmarks beyond the browser
34Everything above ran in your browser through numbl, and the browser is a particular place to do linear algebra: the dense solves run on a TypeScript fallback rather than a native LAPACK, and the loops run through numbl's JIT. How much of the picture is the method, and how much the platform? To answer that, the repository carries offline ports of the two best methods, the harmonic basis and the graded double layer, for real MATLAB, Python with NumPy, and Julia:
36```
37python3 bench/python/bench_pacman.py
38julia bench/julia/bench_pacman.jl
39matlab -batch "run('bench/matlab/bench_pacman.m')"
40```
1a8bb6aMake the cross-language benchmark measure the languagesJeremy Magland 42Each writes a small JSON file of resolutions, errors and times; the "load offline results" button under the work-precision plot draws them into it as dashed curves. The MATLAB harness does not port anything: it runs this chapter's driver and script, unchanged, in real MATLAB, which is the most direct numbl comparison available. The Python and Julia scripts are ports, and the check that they are honest ports is in the plot itself: their error curves must land on the browser's, since error depends only on the mathematics, while their time curves may not.
1a8bb6aMake the cross-language benchmark measure the languagesJeremy Magland 44The reason for the upload, rather than a table of times measured on the authors' machine, is comparability. A timing from someone else's hardware shares no clock with the numbers your browser just produced; a timing from the scripts above, run on the machine you are reading this on, does. Loading a file recorded elsewhere is still allowed, and the legend will say what it was, but the one comparison the plot can defend is browser and native on the same machine.
46### What it takes for such a comparison to mean anything
48A cross-language benchmark measures the author's fluency in each language unless care is taken, and the first version of these scripts measured exactly that. Three disciplines are what make the curves comparable, and they are worth stating because they are easy to get wrong.
50**Each port is written in the fastest shape its own language offers, and those shapes are not the same shape.** The Nyström matrix can be built by a scalar double loop over its entries, or by broadcasting the coordinate differences into $n \times n$ arrays. The loop is the faster shape in three of the four environments here, by about a factor of two against implicit expansion in both MATLAB and numbl and by more than ten against it in Julia, and it is what the chapter's script and both of those ports use. In NumPy the same loop would be hopeless, so the Python port assembles with whole-array expressions, and does so in place: written as the obvious sequence of broadcast statements it allocates five $n \times n$ temporaries and spends longer moving memory than the dense solve spends on arithmetic. Holding the *code* fixed across languages would have been the unfair choice, not the fair one.
52**The same linear algebra must be asked for.** NumPy's `lstsq` is an SVD, while MATLAB's and Julia's backslash are both pivoted QR on a rectangular system. Leaving that alone had Python solving a harder problem than the others under the same label; the port now calls a QR driver. The errors then agree exactly on the double layer, and to within a factor of two at the $10^{-15}$ roundoff floor of the basis method, where the pivoting order is all that separates them.
54**Timing needs a protocol, applied everywhere.** Every environment here, including the page, keeps the fastest of several runs, with at least three runs whenever a solve is cheap enough to allow it. A single measurement of a fifty-millisecond solve varies by a factor of two with whatever else the machine is doing, which is enough to reorder neighbouring points and to make a language look slow that is not.
56With those three in place, on the author's machine (a 10-core i7-1355U), the four environments land within a factor of two of each other, and the remaining spread has causes one can name. MATLAB and Julia are within about ten per cent up to a $1300 \times 1300$ system; beyond that MATLAB pulls ahead, and its BLAS is using ten threads where Julia's default is six, which is why the legend records the thread count. Python trails by 1.3 to 1.7 times, essentially all of it in the assembly, which is the one stage where NumPy cannot express the fast shape. numbl runs about twice the time of native MATLAB on the same script, and much more than that on the smallest problems, where what is being measured is the fixed cost of a call rather than the arithmetic in it.
58Before the three disciplines were in place, the same scripts made Julia look about ten times faster than the other two. Nothing about that number was a property of Julia.