Faster random circuits (#464)

Author: J-C-Q

CHANGELOG.md

@@ -5,6 +5,11 @@
 
 # News
 
+## v0.9.17 - 2025-02-18
+
+- New memory structure and matrix inversion function for `random_destabilizer`, to reduce allocations and speed up repeated generation of many random destabilizers.
+- Improvements to allocations in `apply!`
+
 ## v0.9.16 - 2024-12-29
 
 - 100× faster unbiased `random_pauli`.

Project.toml

@@ -1,7 +1,7 @@
 name = "QuantumClifford"
 uuid = "0525e862-1e90-11e9-3e4d-1b39d7109de1"
 authors = ["Stefan Krastanov <[email protected]> and QuantumSavory community members"]
-version = "0.9.16"
+version = "0.9.17"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"

src/randoms.jl

@@ -171,8 +171,199 @@ random_destabilizer(n::Int; phases::Bool=true) = random_destabilizer(GLOBAL_RNG,
 random_destabilizer(rng::AbstractRNG, r::Int, n::Int; phases::Bool=true) = MixedDestabilizer(random_destabilizer(rng,n;phases),r)
 random_destabilizer(r::Int, n::Int; phases::Bool=true) = random_destabilizer(GLOBAL_RNG,r,n; phases)
 
+# Reuse memory for faster generation of many random destabilizers
+struct RandDestabMemory{N,T<:Integer}
+    F1::Matrix{Int8}
+    F2::Matrix{Int8}
+    hadamard::BitVector
+    perm::Vector{T}
+    had_idxs::Vector{T}
+    perm_inds::Vector{T}
+    U::Matrix{Int8}
+    xzs::BitMatrix
+    phasesarray::Vector{UInt8}
+    phase_options::Vector{UInt8}
+    arr::Vector{T}
+    n::Int64
+
+    function RandDestabMemory(n::Integer)
+        T = typeof(n)
+        F1 = zeros(Int8, 2n, 2n)
+        F2 = zeros(Int8, 2n, 2n)
+        hadamard = falses(n)
+        perm = zeros(T, n)
+        had_idxs = zeros(T, n)
+        perm_inds = zeros(T, 2n)
+        U = zeros(Int8, 2n, 2n)
+        xzs = falses(2n, 2n)
+        phasesarray = zeros(UInt8, 2n)
+        phase_options = [0x0, 0x2]
+        arr = collect(1:n)
+        new{n,T}(F1, F2, hadamard, perm, had_idxs, perm_inds, U, xzs, phasesarray, phase_options, arr, n)
+    end
+end
+function _reset!(memory::RandDestabMemory)
+    fill!(memory.F1, 0)
+    fill!(memory.F2, 0)
+    fill!(memory.hadamard, false)
+    fill!(memory.perm, 0)
+    fill!(memory.had_idxs, 0)
+    fill!(memory.perm_inds, 0)
+    fill!(memory.U, 0)
+    fill!(memory.xzs, false)
+    fill!(memory.phasesarray, 0)
+    for i in eachindex(memory.arr)
+        @inbounds memory.arr[i] = i
+    end
+end
+
+# Allocation free inverse of upper trinagular int matrix with 1 on diagonal.
+function _inv!(inverse, A)
+    for i in 2:size(A, 2)
+        for j in 1:i-1
+            @inbounds factor = A[j, i] # diagonal always 1 => only integers in inverse matrix
+            for k in 1:j
+                @inbounds inverse[k, i] -= factor * inverse[k, j]
+            end
+        end
+    end
+    nothing
+end
+
+function random_destabilizer(rng::AbstractRNG, memory::RandDestabMemory; phases::Bool=true)
+    n = memory.n
+    # reset the working memory
+    _reset!(memory)
+
+    hadamard = memory.hadamard
+    perm = memory.perm
+    _quantum_mallows!(rng, hadamard, perm, memory.arr)
+    had_idxs = memory.had_idxs
+    j = 1
+    for i in eachindex(hadamard)
+        @inbounds if hadamard[i]
+            @inbounds had_idxs[j] = i
+            j += 1
+        end
+    end
+
+    # delta, delta', gamma, gamma' appear in the canonical form
+    # of a Clifford operator (Eq. 3/Theorem 1)
+    # delta is unit lower triangular, gamma is symmetric
+    F1 = memory.F1
+    F2 = memory.F2
+    delta = @view F1[1:n, 1:n]
+    delta_p = @view F2[1:n, 1:n]
+    prod = @view F1[n+1:2n, 1:n]
+    prod_p = @view F2[n+1:2n, 1:n]
+    gamma = @view F1[1:n, n+1:2n]
+    gamma_p = @view F2[1:n, n+1:2n]
+    inv_delta = @view F1[n+1:2n, n+1:2n]
+    inv_delta_p = @view F2[n+1:2n, n+1:2n]
+    for i in 1:n
+        delta[i, i] = 1
+        delta_p[i, i] = 1
+        gamma_p[i, i] = rand(rng, 0x0:0x1)::UInt8
+    end
+
+    # gamma_ii is zero if h[i] = 0
+    for idx in had_idxs
+        if idx != 0
+            gamma[idx, idx] = rand(rng, 0x0:0x1)::UInt8
+        end
+    end
+
+    # gamma' and delta' are unconstrained on the lower triangular
+    fill_tril(rng, gamma_p, n, symmetric=true) # I think this should be fill_tril!
+    fill_tril(rng, delta_p, n)
+
+    # off diagonal: gamma, delta must obey conditions C1-C5
+    for row in 1:n, col in 1:row-1
+        if hadamard[row] && hadamard[col]
+            gamma[row, col] = gamma[col, row] = rand(rng, 0x0:0x1)::UInt8
+            # otherwise delta[row,col] must be zero by C4
+            if perm[row] > perm[col]
+                delta[row, col] = rand(rng, 0x0:0x1)::UInt8
+            end
+        elseif hadamard[row] && (!hadamard[col]) && perm[row] < perm[col]
+            # C5 imposes delta[row, col] = 0 for h[row]=1, h[col]=0
+            # if perm[row] > perm[col] then C2 imposes gamma[row,col] = 0
+            gamma[row, col] = gamma[col, row] = rand(rng, 0x0:0x1)::UInt8
+        elseif (!hadamard[row]) && hadamard[col]
+            delta[row, col] = rand(rng, 0x0:0x1)::UInt8
+            # not sure what condition imposes this
+            if perm[row] > perm[col]
+                gamma[row, col] = gamma[col, row] = rand(rng, 0x0:0x1)::UInt8
+            end
+        elseif (!hadamard[row]) && (!hadamard[col]) && perm[row] < perm[col]
+            # C1 imposes gamma[row, col] = 0 for h[row]=h[col] = 0
+            # if perm[row] > perm[col] then C3 imposes delta[row,col] = 0
+            delta[row, col] = rand(rng, 0x0:0x1)::UInt8
+        end
+    end
+
+    # now construct the tableau representation for F(I, Gamma, Delta)
+    mul!(prod, gamma, delta)
+    mul!(prod_p, gamma_p, delta_p)
+    for i in n+1:2n
+        @inbounds F1[i, i] = 1
+        @inbounds F2[i, i] = 1
+    end
+    _inv!(inv_delta, delta')
+    _inv!(inv_delta_p, delta_p')
+
+    # block matrix form
+    F1 .= mod.(F1, 2)
+    F2 .= mod.(F2, 2)
+    gamma .= 0
+    gamma_p .= 0
+
+    # apply qubit permutation S to F2
+    perm_inds = memory.perm_inds
+    U = memory.U
+    for i in 1:n
+        @inbounds perm_inds[i] = perm[i]
+        @inbounds perm_inds[i+n] = perm[i] + n
+    end
+    for (i, e) in enumerate(perm_inds)
+        for j in 1:2n
+            U[i, j] = F2[e, j]
+        end
+    end
+
+    # apply layer of hadamards
+    for i in 1:n
+        if had_idxs[i] != 0
+            for j in 1:2n
+                @inbounds (U[had_idxs[i], j], U[had_idxs[i]+n, j]) = (U[had_idxs[i]+n, j], U[had_idxs[i], j])
+            end
+        end
+    end
+
+    # apply F1
+    xzs = memory.xzs
+    fill!(F2, 0)
+    mul!(F2, F1, U)
+
+    F2 .= mod.(F2, 2)
+    xzs .= F2 .== 1
+
+    # random Pauli matrix just amounts to phases on the stabilizer tableau
+    phasesarray = memory.phasesarray
+    if phases
+        for i in 1:2n
+            phasesarray[i] = rand(rng, memory.phase_options)
+        end
+    end
+
+    return Destabilizer(Tableau(phasesarray, xzs))
+end
+
+
+
 """A random Clifford operator generated by the Bravyi-Maslov Algorithm 2 from [bravyi2020hadamard](@cite)."""
 random_clifford(rng::AbstractRNG, n::Int; phases::Bool=true) = CliffordOperator(random_destabilizer(rng, n; phases))
+random_clifford(rng::AbstractRNG, memory::RandDestabMemory; phases::Bool=true) = CliffordOperator(random_destabilizer(rng, memory; phases))
 random_clifford(n::Int; phases::Bool=true) = random_clifford(GLOBAL_RNG, n::Int; phases)
 
 """A random Stabilizer tableau generated by the Bravyi-Maslov Algorithm 2 from [bravyi2020hadamard](@cite)."""
@@ -196,23 +387,49 @@ function nemo_inv(a, n)::Matrix{UInt8}
     return collect(UInt8.(inverted.==1)) # maybe there is a better way to do the conversion
 end
 
+
 """Sample (h, S) from the distribution P_n(h, S) from Bravyi and Maslov Algorithm 1."""
 function quantum_mallows(rng::AbstractRNG, n::Int) # each one is benchmarked in benchmarks/quantum_mallows.jl
     arr = collect(1:n)
     hadamard = falses(n)
     perm = zeros(Int64, n)
+
+    return _quantum_mallows!(rng, hadamard, perm, arr)
+end
+
+function _quantum_mallows!(
+    rng::AbstractRNG,
+    hadamard::BitVector,
+    perm::Vector{T},
+    arr::Vector{T}) where {T<:Integer} # inplace version
+
+    n = length(perm)
     for idx in 1:n
-        m = length(arr)
+        m = n - idx + 1
         # sample h_i from given prob distribution
         l = sample_geometric_2(rng, 2 * m)
         weight = 2 * m - l
         hadamard[idx] = (weight < m)
-        k = weight < m ? weight : 2*m - weight - 1
-        perm[idx] = popat!(arr, k + 1)
+        k = weight < m ? weight : 2 * m - weight - 1
+        perm[idx] = _popat!(arr, k + 1)
     end
     return hadamard, perm
 end
 
+function _popat!(arr, n)
+    m = length(arr)
+    at = arr[n]
+
+    for i in n:(m-1)
+        arr[i] = arr[i+1]
+    end
+
+    arr[end] = 0
+
+    return at
+end
+
+
 """ This function samples a number from 1 to `n` where `n >= 1`
     probability of outputting `i` is proportional to `2^i`"""
 function sample_geometric_2(rng::AbstractRNG, n::Integer)
@@ -263,6 +480,7 @@ function random_brickwork_clifford_circuit(rng::AbstractRNG, lattice_size::NTupl
     cartesian = CartesianIndices(lattice_size)
     dim = length(lattice_size)
     nqubits = prod(lattice_size)
+    working_memory = RandDestabMemory(2)
     for i in 1:nlayers
         gate_direction = (i - 1) % dim + 1
         l = lattice_size[gate_direction]
@@ -273,7 +491,7 @@ function random_brickwork_clifford_circuit(rng::AbstractRNG, lattice_size::NTupl
                 cardk = cardj
                 cardk[gate_direction] = cardk[gate_direction] + 1
                 k = LinearIndices(cartesian)[cardk...]
-                push!(circ, SparseGate(random_clifford(rng, 2), [j, k]))
+                push!(circ, SparseGate(random_clifford(rng, working_memory), [j, k]))
             end
         end
     end

test/test_allocations.jl

@@ -1,6 +1,6 @@
-@testitem "Allocation checks" tags=[:alloccc] begin
-    using QuantumClifford
-    using QuantumClifford: mul_left!
+@testitem "Allocation checks" begin
+    using QuantumClifford: mul_left!, RandDestabMemory, Tableau
+    using Random
     n = Threads.nthreads()
     allocated(f::F) where {F} = @allocations f()
     @testset "apply! mul_left! canonicalize!" begin
@@ -44,6 +44,20 @@
             @test allocated(f6) <= 2
         end
     end
+    @testset "random_destabilizer" begin
+        N = 100
+        memory = RandDestabMemory(N)
+        f1() = RandDestabMemory(N)
+        f1()
+        f2() = random_destabilizer(Random.GLOBAL_RNG, memory)
+        f2()
+        f3() = Destabilizer(Tableau(memory.phasesarray, memory.xzs))
+        f3()
+        @test allocated(f1) < 12.5 * N^2 + 50 * N + 1000
+        if VERSION >= v"1.11"
+            @test abs(allocated(f2) - allocated(f3)) / allocated(f3) < 0.05
+        end
+    end
     @testset "project!" begin
         N = 100
         d = random_destabilizer(N)

test/test_random.jl

@@ -1,4 +1,5 @@
 @testitem "Random" begin
+    using Random
     using QuantumClifford
     using QuantumClifford: stab_looks_good, destab_looks_good, mixed_stab_looks_good, mixed_destab_looks_good
 
@@ -32,6 +33,21 @@
         end
     end
 
+    @testset "Random sampling of operators memory reuse" begin
+        for n in [1, test_sizes..., 200, 500]
+            workingmemory = QuantumClifford.RandDestabMemory(n)
+            for _ in 1:2
+                seed = rand(1:100000)
+                rng = Random.GLOBAL_RNG
+                Random.seed!(rng, seed)
+                non_reuse_version = random_destabilizer(rng, n)
+                Random.seed!(rng, seed)
+                reuse_version = random_destabilizer(rng, workingmemory)
+                @test non_reuse_version == reuse_version
+            end
+        end
+    end
+
     @testset "Random Paulis" begin
         for n in [1, test_sizes..., 200,500]
             @test all((random_pauli(n).phase[] == 0  for _ in 1:100))