Faster random circuits

[J-C-Q]

Motivation

Motivated by random circuits, generating a random Clifford operation can be a time-critical function. When trying to implement such circuits using QuantumClifford.jl I noticed many allocations from the random_clifford function. So I tried to implement a version of the random_clifford function that reuses memory, so that repeated generation of random gates doesn't allocate that much.

Implementation

This PR implements a RandDestabMemory type which stores all the memory needed in the already existing random_destabilizer algorithm. An instance of this type can be passed to the random_destabilizer function to reuse the same memory. The algorithm is should be unchanged; however, this PR changes the code in some places to avoid allocations. Some of the functionality gets separated into helper functions, but some loops persist in the main function, which arguably hurts readability. To avoid as many allocations as possible, this PR also rewrites a few other functions like quantum_mallows to reuse memory. Most notably, the calculation of the inverse of delta' and delta_p' doesn't rely on external packages anymore and makes use of the respective structure of the matrices, which increases performance.
This PR also implements this new version of random_destabilizer in the random_brickwork_clifford_circuit function.

TODO

There is still some cleanup to be done. As well as documentation and testing. However, I thought I'd share this idea already to get feedback if you are even interested in integrating this. I'm very open to any criticism and ideas to improve my code, as well as happy to work on this more.

Benchmark

Attached you can find some rudimentary benchmarks of the random_brickwork_clifford_circuit function, first without and then with the new version implemented.

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 83.42%. Comparing base (0d13791) to head (3ee7299).

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #464      +/-   ##
==========================================
+ Coverage   82.73%   83.42%   +0.68%     
==========================================
  Files          70       72       +2     
  Lines        4656     4801     +145     
==========================================
+ Hits         3852     4005     +153     
+ Misses        804      796       -8     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[J-C-Q]

I think this is ready for review.

• I haven't put a changelog entry yet, but can do so if you want me to do it.
• Also, I'm not sure if you want this to be documented in the random_destabilizer doc string, or somewhere else.
• The tests are very basic, just comparing the result to the already existing implementation. Do you want me to add more or is that enough? The existing tests seem to be passing.

Thank you for the positive reaction to this idea! :)

[Krastanov]

This is fantastic, thank you!

I have left a few comments in, mostly focusing on the custom functions you have implemented and checking whether they are necessary on latest julia.

If you can also just add a changelog and bump the version number, I should be able to merge and release this quite quickly.

[Krastanov]

May we add another testset that uses something like @allocated or @allocations to keep track of potential regressions in performance.

[Krastanov]

now that this is very thoroughly tested, let's drop this down to for _ in 1:2 to not make the tests take too long.

[Krastanov]

could you confirm that x .= mod.(x, n) is worse than a new function. On newer julias I get this:

julia> a = rand(Int8, 20,20);

julia> @time _inplace_mod!(a, 5);
  0.000001 seconds

julia> a = rand(Int8, 200,200);

julia> @time _inplace_mod!(a, 5);
  0.000054 seconds

julia> function f(a,n)
       a.=mod.(a,n)
       end;

julia> @time f(a, 5);
  0.017539 seconds (13.95 k allocations: 738.812 KiB, 99.65% compilation time)

julia> @time f(a, 5);
  0.000056 seconds

Both seem equally fast and non-allocating. I prefer we just use broadcast instead of a new _inplace_mod! if there is indeed no difference in performance

[Krastanov]

same question about using broadcast here

julia> a = trues(10,10);

julia> b = rand(Int8, 10, 10);

julia> @time a .= b.==5;
  0.047673 seconds (213.96 k allocations: 10.826 MiB, 99.95% compilation time)

julia> @time a .= b.==5;
  0.000004 seconds (1 allocation: 32 bytes)

julia> function _inplace_equals!(x::BitMatrix, y::Matrix{Int8}, n::Integer)
           @inbounds for i in eachindex(y)
               x[i] = y[i] == n
           end
           return x
       end;

julia> @time _inplace_equals!(a,b,5);
  0.000001 seconds

julia> f(a,b,n) = (a.=b.==n);

julia> @time f(a,b,5);
  0.032928 seconds (26.55 k allocations: 1.273 MiB, 99.97% compilation time)

julia> @time f(a,b,5);
  0.000003 seconds

[Krastanov]

what about LinearAlgebra.mul! -- is it measurably worse?

[Krastanov]

to avoid code repetition (really, more to avoid the need to fix a bug in two separate places in the future), could you delete most of the content of quantum_mallows and have it internally just call your new quantum_mallows! together with whatever objects needs to be allocated?

[J-C-Q]

Thank you for the thoughtful feedback!

I've made the following changes based on your suggestions:

• Removed the unnecessary _inplace_mod!, _inplace_equals, and _mul! functions, as they were not faster than native broadcasting and LinearAlgebra implementations.
• Reduced the number of random tests and added a test to track allocations.
• Updated quantum_mallows to call the in-place version.
• Added a changelog entry and bumped the version number (hopefully, I did this correctly!).

Thanks again for your support—I really appreciate it!

[J-C-Q]

Well, apparently there is a big difference in allocations between 1.10 and 1.11. Somehow LinearAlgebra.mul! allocates in 1.10 but doesn't for 1.11 (at least for this use case).

[Fe-r-oz]

Well, apparently there is a big difference in allocations between 1.10 and 1.11. Somehow LinearAlgebra.mul! allocates in 1.10 but doesn't for 1.11 (at least for this use case).

This Julia issue (@allocated false positive on v1.12?) is due to false positives by @allocated. This diagnostic PR might be helpful: #467 (comment). There are some errors in CI documentation as well that are also reproduced in diagnostic PR #467. I think this issue started happening by the time Julia 1.11.3 was released on CI here.

Edit:

The CI documentation errors are all due to errors to run the @example blocks possibly stemming from QCMakieExt

Error: failed to run `@example` block in src/canonicalization.md:16-21
```@example
using QuantumClifford, CairoMakie
f=Figure()
stabilizerplot_axis(f[1,1], canonicalize!(random_stabilizer(20,30)))
f

Error:

ArgumentError: component type FixedPointNumbers.N0f8 is an 8-bit type representing 256 values from 0.0 to 1.0,
      but the values (0x001b9e77, 0x001b9e77, 0x001b9e77) do not lie within this range.
      See the READMEs for FixedPointNumbers and ColorTypes for more information.
    Stacktrace:

The CI doc errors were due to this line:

QuantumClifford.jl/ext/QuantumCliffordMakieExt/QuantumCliffordMakieExt.jl

Line 20 in 0d13791

colormap = Makie.cgrad([:lightgray,Makie.RGB(0x1b9e77),Makie.RGB(0xd95f02),Makie.RGB(0x7570b3)], 4, categorical = true),

The CI documentation errors are fixed in #473 :)

[Krastanov]

As long as it is working fine on 1.11, I am not too worried about 1.10. In this library I feel comfortable relying on improvements in base julia to get simpler code (admittedly a freedom that bigger foundational libraries can not afford to themselves).

Could you summarize whether there are any outstanding issues on 1.11?

[Krastanov]

The allocation test fails are unrelated to this PR, right?

[J-C-Q]

Alright, since LinearAlgebra.mul! does not allocate and is faster in 1.11, we should use it.
To test allocations as expected, I changed the new test to only run on 1.11 and upwards. It seems to pass there. However, there are some other allocation tests that do not pass right now, but they are not related to this PR!

[Krastanov]

oh, one last thing: in your new allocation tests, make sure that you have called f1,f2,f3 once before you run the allocation checks on them, to get them compiled (otherwise depending on what the compiler does, you might end up measuring its performance instead of just f1,f2,f3)

[Krastanov]

Wonderful! (well, not great about the new failures, but you are right, they are unrelated)

I will plan to merge this shortly and look into the new failures separately.

Thank you for the contributions, it is much appreciated!