Replace Sobol sequences with SZ sequences

[wantonsushi]

Hello again!

This implements SZ sequences as discussed in #103, replacing the generator matrices behind SobolSampler but keeping the name.

Sorry this took a while. The prototype I posted in #103 had a problem I only caught later: dimensions 0 and 1 stopped matching Sobol.

I had built each dimension as the base-4 Pascal matrix of a GF(4) symbol, straight from eq (24). That is the set the paper prints in eq (23), and it is a valid (0, 4)-sequence, but its dimension 1 is the base-4 Pascal matrix rather than the canonical binary one, so SZ[1] != P.

The fix is to pre-multiply each 2x2 block by [[1,1],[0,1]], the self similar block of the binary Pascal matrix. Dimension 1 then reproduces that matrix exactly, dimensions 0 and 1 are bit-identical to Sobol again, and SZ[0] = I, SZ[1] = P holds. This is the Sobol nesting from section 4.3. The factor is block diagonal, so it cannot disturb the base-4 stratification.

Changes

• include/oqmc/owen.h: SZ programs replace the Sobol ones. Two paths as agreed: CPU is the shared right-shift form with the reversal hoisted out of the draw and eq (13) for dimension 3, GPU keeps the left-shift form with the reversal on the output and its own dimension 3 program. The SSE, AVX and Neon paths are gone. The build options, oqmc/arch.h and the CI variants stay for now, since I wasn't certain what to do with them.
• src/tools/cli/matrices.cpp: constructs the SZ matrices and prints the programs and the reference table, as it did for Sobol. The Gruenschloss matrices are removed.
• src/tests/owen.cpp: reference table is now SZ. ShirleyRemapping is unchanged and still passes. 02Sequence now covers (2, 3) as well as (0, 1), which pins the eq (13) relation, and all three stratification tests run at the full 16 bit precision of the index rather than the leading half. Adds ~2s to the suite in Debug.
• src/tools/lib/plot.cpp, python/wrapper.py: plot_error takes a dimension pair, defaulting to (0, 1), so the plots below are possible.

Conveniently, dimensions 0 and 1 are bit-identical to before, only 2 and 3 change. Thus, the blue noise tables stay valid. optimise baked them from shuffledScrambledSobol<2>, so no re-bake needed since they were never tuned for dimensions 2 and 3 anyways.

Stratification

256 points on a 16x16 grid, cells not holding exactly one point marked red. Sobol fails on (1, 3) and (2, 3), SZ passes on all six pairs. SZ also passes the 4x64 and 64x4 splits, where Sobol fails (0, 3) and (2, 3).

Convergence

plot_error, 1024 realizations, 4096 samples, rmse ratio sobol/sz, so >1 means SZ is better. Geometric mean over the last five points of each series, split by sample count since the two differ either side of a power of four:

shape	pair	at 4^k	at 2*4^k
quarter disk	(0,1)	1.00	1.00
quarter gaussian	(0,1)	1.00	1.00
bilinear	(0,1)	1.00	1.00
orientated heaviside	(0,1)	1.00	1.00
quarter disk	(2,3)	1.11	1.13
quarter gaussian	(2,3)	1.38	1.72
bilinear	(2,3)	1.53	2.02
orientated heaviside	(2,3)	1.07	1.12

(0, 1) is exactly 1.00 because those curves are bit-identical. (2, 3) is better on every shape, and at both sample counts.

Timings

benchmark <sampler> samples, gcc 16.1 release, median of 25 interleaved runs pinned to one core, microseconds:

sampler	old scalar	old SSE	old AVX	new
sobol	56532	47944	48385	57649
sobolbn	68471	59649	60308	63336

Level with the old scalar path and still behind the SIMD paths we are dropping, so this is a quality change rather than a speed one at 16 bits.

On the eq (13) chaining, I kept it, but testing against giving dimension 3 its own 9 step program, ratio of independent to chained, so >1 means chaining is faster:

compiler	sobol	sobolbn
gcc 16.1	0.93	1.03
gcc 15	1.01	1.06
clang 22.1	0.99	1.03

sobolbn prefers the chain on all three, by 3% to 6%. sobol is within noise, and gcc 16 is the outlier where the independent program came out 6.7% ahead. So I think chain looks like the better default, and it is fewer constants.

GPU

RTX 4070 Laptop, benchmark <sampler> samples, median of 21 interleaved runs, microseconds:

sampler	old	new
sobol	332	312
sobolbn	844	802

Chaining dimension 3 on the device measured 8% and 4% worse, so I gave it its own program.

NB: these timings needed a benchmark fix first. sobol and lattice have cacheSize == 0, so initialiseCache() issues no device work and CUDA context creation lands inside the timed region: sobol samples reads 151538 against 312 with a warm-up launch. I opened an issue here: #109

Verification

• OwenTest.SobolReversedIndex checks the closed forms against the matrix products for every 16-bit index and dimension.
• The matrices tool output matches the owen.h programs and the test table. The tool also prints a host program for dimension 3, which owen.h does not use since it chains.
• generate sobol is byte-identical across scalar, AVX, Debug and GPU builds, and dimensions 0 and 1 are byte-identical to main.
• Tests pass (173/173), clang-format and clang-tidy clean.

Not done here: the README performance plots and pair plots will require regenerating, and need to decide on removing the SIMD build options and CI variants.

Edit: added GPU timings.

[fpsunflower]

Nice! Its cool to see the supporting code for how the matrices are derived.

Would it be possible to adapt the plot script to check the other dimension pairs as well? The pairs (0,2),(0,3),(1,2),(1,3) are not expected to be as good as (0,1) and (2,3) but it would be interesting to see how they behave compared to sobol.

Of course the final test is to check inside a full blown renderer ... In my experience its helpful to check things like motion blurred scenes with objects moving in all directions or volumetric shadows to catch particular interactions between dimension pairs.

[wantonsushi]

@fpsunflower Yes, I already adapted the plot script to take any pair of dimensions. Here are the remaining pairs, same protocol as before: 1024 realizations, 4096 samples, rmse ratio sobol/sz (so >1 means SZ is better). Geometric mean
over the last five points of each series, split by sample count.

shape	pair	at 4^k	at 2*4^k
quarter disk	(0, 2)	1.02	0.98
quarter gaussian	(0, 2)	0.92	0.71
bilinear	(0, 2)	0.93	0.68
orientated heaviside	(0, 2)	1.03	0.93
quarter disk	(0, 3)	1.19	1.13
quarter gaussian	(0, 3)	1.46	1.10
bilinear	(0, 3)	1.56	1.13
orientated heaviside	(0, 3)	1.03	1.00
quarter disk	(1, 2)	0.97	0.97
quarter gaussian	(1, 2)	1.09	0.63
bilinear	(1, 2)	1.08	0.57
orientated heaviside	(1, 2)	0.99	0.97
quarter disk	(1, 3)	1.03	0.90
quarter gaussian	(1, 3)	1.25	0.88
bilinear	(1, 3)	1.32	0.86
orientated heaviside	(1, 3)	1.14	0.92

I tried the same process with a few different seeds and every cell moved by at most 0.02.

I'll test the differences in my hobby engine when I get the chance. I'll test motion blur / volumetric shadows, as advised.

[joshbainbridge]

Those plots look great. This will be a great addition. Thank you for putting in the time here to also validate so thoroughly.

The implementation in oqmc/owen.h is good. I’ve added a few comments, but just on code layout more than anything. Leaving oqmc/arch.h and the build setup in place I think is the right call for now. We might want to leverage that again. If enough time goes by and we don’t, it can always be removed at a later date.

I was expecting we would need to change dimensions 0 and 1. Nice work on keeping that consistent and avoiding the need to update the blue noise tables.

Also great to see the large Sobol matrices gone. Always happy to have less code in the repo, and having the full derivation in the codebase.

Tests are good. I’d agree, running at full precision is the right call, and worth the 2 second cost.

One thing to note, we likely need to regenerate the readme-images. If you are using the Nix flake env, docker container, or have Jupyter Notebook available, then this can be done with just readme-images. If not, then I’m happy to help out and run that once you're happy with the PR.

[wantonsushi]

Pushed changes. Most of it is as suggested in your review:

• owen.h: I went with your suggested names szToReversed() and szReversedBasis(), and both functions are defined everywhere (no #if). Dimension 3 has its own program in both, chaining only in shuffledScrambledSobol.
• tests/owen.cpp: added a test for each function checking it against the matrix table, and a third test for the chaining. stratumIndex gets the assert and the twoPow32 you suggested, though I kept the int parameter since that's what the callers all pass.
• matrices.cpp: switched to std::array, added the comment above matrices, and renamed the printer to printReversedMatrices.
• CHANGELOG.md: mentions the rename.
• Regenerated the readme images (only pair-plot-sobol should be affected since 0,1 is the same?).

Also, one thing I'd like to mention: I tried this SZ implementation in my own engine and in my pbrt fork, and any scene where a 3D draw is used heavily, it looks to be slightly worse than Sobol at most sample counts. But in general, I'm pretty happy with this code.

[joshbainbridge]

Hey @wantonsushi. I'd agree, code changes look great. Thank you for addressing those comments. I think this is almost ready to go in. Just need a rebase and squash with a descriptive commit message.

On the practical tests, that is interesting. We'll run some tests as well. It might be an unlucky pair (xu, xv, yu, yv) does not project well at non-power-of-four. Eyeballing the updated plot in the diff, xv appears to have some notable gaps at 128 samples. But I think this is still a step in the right direction, as I'm more confident in the guaranttes that SZ provides over what we had with traditional Sobol matrices from S. Joe and F. Y. Kuo.