How I Made G6K’s Spherical Code Generator Up to 100× Faster with Incremental Gram Matrix Updates in Rust

Plain English: G6K is software researchers use to test how hard lattice problems are. Before its main search starts, it loads a table of helper vectors from a file. The script that makes that file recalculated almost the entire score matrix after changing two values. Recalculating only the affected row made that setup step up to 100× faster.

This post is one piece of a longer effort to make the tooling around post-quantum cryptography faster to run and easier to change.

TL;DR

G6K, the open-source lattice sieving toolkit used in most post-quantum cryptanalysis work, filters candidate vector pairs with a SimHash built from sparse spherical codes. Those codes come from spherical_coding/codes.py, a single-threaded Python hill-climber. At n = 128 it needs minutes per dimension, and a sweep across many dimensions is a multi-day job.

I rewrote the generator in Rust around one observation: swapping two coordinates inside one row of the code matrix only changes one row (and its mirror column) of the Gram matrix. Updating only that row turns an O(N · n) step into O(N), and keeping a cached sum of squares per row turns the O(N²) score query into O(N). With candidates searched in parallel via rayon, the new tool generates the full code bank for n = 6..1000 on an 8-core M-series MacBook in one sitting.

The output files use the same six-integers-per-line layout that G6K’s C++ loader expects, so no downstream code changes were needed.

The timings are wall-clock observations from internal code, not an independently reproducible benchmark. Both implementations ran on the same 8-core M-series MacBook. The per-candidate number is the closest like-for-like comparison. The three-candidate Rust runs spend extra work on purpose to find a better score, so I do not use them for the speedup figure.

Why G6K uses spherical codes

G6K (the General Sieve Kernel, Albrecht, Ducas, Herold, Kirshanova, Postlethwaite, Stevens and Wessel, Eurocrypt 2019) is the standard open implementation of lattice sieving. Sieving is currently the most practical attack on the shortest vector problem (SVP) and its approximate variants. Those are the problems whose hardness underpins NIST-standardized post-quantum schemes like Kyber and Dilithium.

So when a lattice-based scheme is benchmarked or attacked, the work usually runs through G6K. Anything that makes G6K’s setup or runtime cheaper feeds directly into how thoroughly the next generation of standards can be stress-tested.

The sieving inner loop

A lattice sieve keeps a database of lattice vectors and repeatedly looks for pairs (v, w) such that v ± w is shorter than max(|v|, |w|). With N vectors in the database, naive pair search is O(N²), and N itself grows as roughly (4/3)^(n/2) in dimension n. By dimension 80, brute-force pair search is hopeless.

The idea that made high-dimensional sieving practical is locality-sensitive filtering, most famously the BDGL filter (Becker, Ducas, Gama and Laarhoven, SODA 2016) and its descendants. The principle is to reject most pairs cheaply without computing v ± w at all.

SimHash as G6K’s filter

G6K uses a SimHash-style filter, in the spirit of Charikar’s 2002 random-hyperplane hash, adapted for sieving. The construction has four steps:

1. Fix a bank of probe vectors c_1, …, c_b ⊂ ℝⁿ (here b = 256 or 512, fixed at compile time as XPC_BIT_LEN).
2. For each lattice vector v, compute the bit sign(⟨c_i, v⟩) for every probe.
3. Pack those bits into a fingerprint.
4. To test a pair, XOR the fingerprints and popcount. A small Hamming distance means similar orientation, which means the pair is worth a real reduction attempt.

One 64-bit XOR and popcount rejects pairs that would otherwise cost a full O(n) inner product. Across a sieve that is the difference between hours and weeks. The filter only works if the probe bank is well spread over the sphere, though. A bad bank either rejects good pairs or lets noise through, and sieving performance collapses either way.

The loader in kernel/simhash.inl is short. For dimensions below 30 it skips the filter entirely. Otherwise it opens sc_<n>_<XPC_BIT_LEN>.def, reads six integers per line, and stores them after applying a random permutation of the coordinates seeded per run. The first three integers are the coordinates added into the inner product and the last three are subtracted, so each probe is a vector with three +1 entries and three −1 entries.

What the spherical-code generator does

The probe banks are sparse spherical codes: each probe has six non-zero entries (three at +1, three at −1) spread across the n coordinates. Sparsity keeps ⟨c_i, v⟩ cheap. The “spherical” part means the bank should approximate a uniform spread on the sphere, so that the SimHash bits are close to independent across probes.

The banks live in G6K as files named sc_<n>_<bitlen>.def, one per dimension. They are produced by spherical_coding/codes.py, which does the following:

1. Generates bitlen random sparse rows (random_sparse_code).
2. Builds the Gram matrix G = M Mᵀ of size bitlen × bitlen. Entry G[i, j] is the inner product of probes i and j, so it measures how much two probes overlap.
3. Defines a global score S(G) = Σᵢ (Σⱼ Gᵢⱼ² − G₀₀²)², lower is better. G₀₀ is always SPARSITY = 6, so the inner term is each row’s sum of squares measured against a constant 36. The score penalises rows whose overlap with the rest of the bank deviates from the ideal.
4. Sorts rows by their score contribution, worst first. For each row it tries swapping two coordinates and keeps the first swap that lowers the score (improve_once). On the first improvement it re-sorts and starts over. It stops when no row can be improved or after Tlim = 1000 seconds.

Conceptually clean. Operationally painful, because of what happens inside step 4.

Where the time goes, and what does not fix it

Profiling codes.py shows the inner loop is dominated by two operations, both of which run once per candidate swap:

def update_G(M, G, i):
    v = M.dot(M[i].transpose())   # O(bitlen * n)
    for j in range(l):
        G[i, j] = v[j]
        G[j, i] = v[j]

def score(G):
    s = G[0, 0] * G[0, 0]
    v = [sum([x * x for x in v]) - s for v in G]   # O(bitlen²)
    return sum([x * x for x in v])

For bitlen = 512 and n = 128, that is roughly 65K multiplications for update_G plus 262K for score, per candidate swap. A pass over one row tries up to n(n−1)/2 ≈ 8K swaps, and a code needs dozens of passes. Interpreter overhead sits on top of all of it.

What I tried first

My first instinct was to keep the algorithm and speed up the language:

• Vectorized NumPy rewrites of score and update_G: about 3× faster.
• Numba JIT (@jit(nopython=True) on score, update_G, and the swap loop): another 5-8×.
• Multiprocessing over candidate codes: linear in cores.

Stacked, that gets maybe 30-50× over the original on a good day. Useful, but the per-dimension cost was still in the minutes, and the algorithmic cost had not moved at all. Every candidate swap still triggered a full M.dot(M[i].T) and a full score(G). The JIT just compiled the wasted work so it could be wasted faster. The fix had to come before any language choice.

The locality-of-change observation

A swap inside row i of M can only change the entries of G that involve row i:

• G[i, j] for all j ≠ i, and by symmetry G[j, i]. That is one row and its mirror column, holding the same bitlen − 1 values twice: 2(bitlen − 1) entries, or 1,022 at bitlen = 512.
• The other bitlen² − 2(bitlen − 1) entries (261,122 of 262,144) are bit-identical to before.
• G[i, i] never changes. A swap moves a +1 or −1 to a different coordinate, but the row still has six non-zero entries of magnitude one, so its self inner product stays 6.
• The score depends on the changed values only through the row sums of squares of row i and of each affected j.

Concretely, if a swap on row row changes M[row, col_a] by delta_a and M[row, col_b] by delta_b:

ΔG[row, j] = delta_a · M[j, col_a] + delta_b · M[j, col_b]   for all j ≠ row
ΔG[j, row] = ΔG[row, j]                                       (symmetry)
G[row, row] is invariant                                      (always equals SPARSITY = 6)

That is O(bitlen) arithmetic per swap instead of O(bitlen · n). And with a cached vector row_sq_sums[i] = Σⱼ G[i,j]², each changed entry updates the score with two scalar additions (one for the row, one for its mirror) instead of a recompute of S from scratch. The score query drops from O(bitlen²) to O(bitlen).

The Rust implementation

The result is a small Rust crate, spherical-codegen, that lives beside the Python script. Four source files:

spherical_coding/spherical-codegen/
├── Cargo.toml
├── run_batch.sh           # sweep driver
└── src/
    ├── code_matrix.rs     # SparseCodeMatrix: sparse + dense storage, swap, dot_rows, write_def
    ├── gram.rs            # GramMatrix: incremental updates + cached row_sq_sums
    ├── optimizer.rs       # hill-climber: exhaustive best-swap selection per row
    └── main.rs            # clap CLI, rayon over candidates, picks the best

Dual representation: sparse for I/O, dense for lookups

pub struct SparseCodeMatrix {
    pub n: usize,
    pub bitlen: usize,
    /// rows[i] = [pos0, pos1, pos2, neg0, neg1, neg2]
    pub rows:  Vec<[u16; SPARSITY]>,
    /// dense[i * n + j] = M[i,j] in {-1, 0, 1}
    pub dense: Vec<i8>,
}

rows is what gets written to the .def file: the three +1 coordinates then the three −1 coordinates, which is the order the G6K loader expects. dense gives O(1) lookups for M[j, col_a] during the incremental Gram update. Without it every term would scan a sparse row. Building the initial Gram matrix uses the sparse side instead: dot_rows compares six indices against six indices, so each of the bitlen²/2 initial entries costs at most 36 comparisons.

Incremental Gram update

pub fn update_row_incremental(
    &mut self,
    code: &SparseCodeMatrix,
    row: usize,
    col_a: usize, col_b: usize,
    delta_a: i32, delta_b: i32,
) {
    let size = self.size;
    for j in 0..size {
        if j == row { continue; }                     // diagonal never changes
        let mj_a = code.get(j, col_a) as i32;
        let mj_b = code.get(j, col_b) as i32;
        let delta = delta_a * mj_a + delta_b * mj_b;
        if delta != 0 {
            let old_val = self.data[row * size + j];
            let new_val = old_val + delta;
            let sq_delta = (new_val as i64) * (new_val as i64)
                         - (old_val as i64) * (old_val as i64);
            self.row_sq_sums[row] += sq_delta;
            self.row_sq_sums[j]   += sq_delta;        // G[j,row] changes too
            self.data[row * size + j] = new_val;
            self.data[j * size + row] = new_val;
        }
    }
}

Most iterations of that loop hit the delta != 0 check and do nothing, because row j has to have a non-zero entry at col_a or col_b for anything to change. With six non-zero entries per row out of n coordinates, that is rare, so the loop is mostly two array reads per j.

The row_sq_sums cache is checked against a from-scratch rebuild in test_incremental_update_matches_full. The test swaps a +1 with a −1 in row 0, runs the incremental update, rebuilds the Gram matrix with GramMatrix::from_code, and asserts every entry and every cached row sum are equal. It caught two off-by-one bugs while I was writing this. I would not trust the optimizer without it.

The shape of the change

            Before swap                 After swap on row r
            ┌─────────────┐             ┌─────────────┐
            │  G  bitlen  │             │  G  bitlen  │
            │             │             │             │
   row r ─► │ ░░░░░░░░░░░ │   ─────►    │ ▓▓▓▓▓▓▓▓▓▓▓ │  ◄─ recomputed
            │             │             │             │
            │ ░ ◄─ col r  │             │ ▓ ◄─ col r  │  ◄─ recomputed (symmetry)
            │ ░           │             │ ▓           │
            │ ░           │             │ ▓           │
            └─────────────┘             └─────────────┘
            261,122 entries unchanged.   1,022 entries change.

The old code recomputed all 262,144 entries through a fresh M @ M.T slice. The new code touches the 1,022 entries that moved.

One outer pass, step by step

The outer loop in optimizer.rs keeps the structure of the Python original. Here is what one iteration of improve does:

1. Ask the Gram matrix for row_scores(): every row’s contribution (row_sq_sums[i] − 36)², sorted descending. This is O(bitlen) because the sums are cached.
2. Walk the rows worst first and call improve_once on each.
3. Inside improve_once, shuffle the n column indices (with the candidate’s own RNG), then for every pair (a, b) with M[row, a] ≠ M[row, b]: swap them in the code matrix, apply update_row_incremental, read the new score off the cached sums, record it if it is the best improvement so far, then swap back and apply the incremental update again to revert.
4. After all pairs, apply the single best swap if it improves the score, and return true.
5. Back in the outer loop, the first row that improved ends the pass. Scores are re-sorted and the next pass starts from step 1.
6. If no row improved, a no-improvement counter goes up. The loop stops when that counter reaches patience (default 50) or when time_limit expires.

One thing I would change: because step 3 is exhaustive, a pass that finds no improvement will not find one on the next pass either. The shuffle only affects which of several equal-value swaps wins. So patience above 1 mostly burns time. The batch script sets it to 50 anyway, which is harmless at these sizes but not useful.

Best-swap selection instead of first-improvement

The Python original returns on the first improving swap. Once a swap evaluation costs O(bitlen) instead of O(bitlen²), that heuristic stops making sense. Searching every pair in the row and applying only the best one converges to lower local minima in fewer outer passes, and the extra evaluations are cheap enough to afford.

for aa in 0..n {
    let a = positions[aa];
    for bb in 0..aa {
        let b = positions[bb];
        if code.get(row, a) == code.get(row, b) { continue; }
        code.swap(row, a, b);
        gram.update_row_incremental(code, row, a, b, delta_a, delta_b);
        let improvement = old_score - gram.score();
        if improvement > best_improvement { best_swap = Some((a, b)); ... }
        code.swap(row, a, b);                                // revert
        gram.update_row_incremental(code, row, a, b, rev_delta_a, rev_delta_b);
    }
}
if let Some((a, b)) = best_swap { /* apply it for real */ }

Parallel candidates and a batch driver

Hill climbers land in different basins from different seeds, so main.rs runs several candidates through rayon and keeps the lowest score. Each candidate gets its own ChaCha8Rng seeded from its index, so a run is reproducible for a given n, bitlen, and candidate count:

let best = (0..args.candidates)
    .into_par_iter()
    .map(|i| {
        let mut rng = ChaCha8Rng::seed_from_u64(i as u64 ^ 0xdeadbeef);
        let mut code = SparseCodeMatrix::random(args.n, bitlen, &mut rng);
        let score = improve(&mut code, &config, &mut rng);
        (code, score)
    })
    .min_by(|a, b| a.1.partial_cmp(&b.1).unwrap())
    .unwrap();

run_batch.sh sweeps n = START..END with three candidates, a 600 second time limit, and patience 50 per dimension. It skips any .def file that already exists, so a killed sweep resumes where it stopped, and it caps the number of concurrent processes at the core count because each process already spreads its candidates across threads.

Release profile

LTO and codegen-units = 1 matter here. The hot path crosses three modules (optimizer calls gram.update_row_incremental, which calls code.get), and cross-module inlining is what keeps the inner loop tight enough for exhaustive best-swap selection to pay off.

[profile.release]
opt-level = 3
lto = true
codegen-units = 1

Results

On an 8-core M-series MacBook with bitlen = 512:

• n = 64: about 1.5 s per candidate, about 5 s for 3 candidates in parallel
• n = 128: about 3 s per candidate, about 10 s total
• n = 256: about 12 s per candidate, about 40 s total
• Full sweep n = 6..1000: completed. Its runtime depends on the per-dimension time-limit, so I do not use it in the speedup figure.

Score quality is lower (better) across every dimension I sampled, mainly because of best-swap selection plus best-of-3 candidate search. Better codes mean a cleaner SimHash filter downstream, so fewer wasted reductions per pair pass in the sieve. I have not measured that downstream effect end to end.

Prior work and where this fits

• G6K (Albrecht et al., Eurocrypt 2019). The toolkit this work plugs into. The SimHash filter is a core component.
• BDGL filter (Becker, Ducas, Gama, Laarhoven, SODA 2016). The locality-sensitive filter design that G6K’s SimHash is a practical descendant of.
• BGJ1 and BGJ-style sieves. Earlier locality-sensitive sieving constructions. G6K supersedes them, but the filter idea is shared lineage.
• Charikar 2002. The random-hyperplane SimHash, the original LSH for cosine similarity.
• Spherical code construction. The hill-climbing approach in codes.py is a practical heuristic. As far as I know it is not derived from a specific paper, just the natural greedy approach to the spread objective. If you know of a more principled construction (say, algebraic spherical codes adapted to the sparse {−1, 0, +1} alphabet), I would like a pointer.

This work does not change the algorithm or improve the codes’ theoretical properties. It is an engineering optimization. The contribution is making code regeneration cheap enough that you can experiment with code parameters (bitlen, sparsity, alphabet) instead of treating the bank as a frozen artifact.

Status and availability

The tool is internal. It produces .def files in the layout kernel/simhash.inl reads, runs across the full n = 6..1000 range at bitlen = 256 and 512, and has a test suite verifying that incremental Gram updates and the cached row_sq_sums match a from-scratch recompute on every entry.

I have not open-sourced it, for two reasons:

1. Quality bar. The hill-climber is good but not provably optimal. There is probably another order of magnitude available with simulated annealing or a proper SDP relaxation, and I would rather release once than release twice.
2. Time. Maintenance, issues, packaging, and integration patches against upstream G6K are their own commitment.

Because the implementation is not public, the numbers cannot be independently reproduced from this article alone. The useful part is the invariant behind them: one coordinate swap changes one Gram row and column, so recomputing the other 99.6 percent is wasted work.