6. Numerical Certificates

Numerical certificates is what we call protocoles for numerically estimating the truncation error we make during propagation. Currently, we only support a limited number of gates, i.e., PauliRotations and CliffordGates.

using PauliPropagation
using Plots

nq = 25
pstr = PauliString(nq, :Z, 13)

nl = 12
topo = bricklayertopology(nq; periodic=true)

circuit = hardwareefficientcircuit(nq, nl; topology=topo);
nparams = countparameters(circuit)

We defined the usual setup, but now, instead of calling propagate(), we use estimatemse(). This function will take the almost all of our default truncations, as well as a custom truncation function, and the number of Monte Carlo samples to estimate the error with. The one truncation that is not supported here is common min_abs_coeff. It looses its meaning in the average error we compute here. But more about that later.

Testing Weight Truncation

Let us test the mean square error we make via weight truncation, averaging over the full parameter space of the circuit's parametrized gates.

Ws = 0:10
mcsamples = 100_000
@time weight_errors = [estimatemse(circuit, pstr, mcsamples; max_weight=W) for W in Ws];
# where errors are zero, replace with small number to avoid log(0) error
weight_errors[weight_errors .<= 1e-10] .= 1e-10

 14.141109 seconds (4.76 M allocations: 257.939 MiB, 2.85% gc time, 16.89% compilation time)

plot(Ws, weight_errors, yscale=:log10, xlabel="Max Pauli Weight", ylabel="MSE", ylim=(1e-5, 1e0), label="", linewidth=2, marker=:circle)

We see that the average error decause quite quickly with the maximum Pauli weight, in fact roughly exponentially. Keep in mind though, that the simulation time with propagate will also scale exponentially in this truncation.

Testing Frequency Truncation

What we call frequency takes the place of the min_abs_coeff coefficient truncation in the average case. We define frequency as the number of splits/branchings at PauliRotation gates. If a path has split $l$ times, then the average coefficient of that path will be $(\frac12)^l$. Note however that the error we measure here is an over-estimation! Because we merge paths in propagate(), we can effectively propagate higher frequency paths with a lower frequency path at no additional computational cost. We don't do merging in the numerical certificate, which is why the average error will be over-estimated.

freqs = 20:20:200
mcsamples = 100000

@time freq_errors = [estimatemse(circuit, pstr, mcsamples; max_freq=freq) for freq in freqs]
# where errors are zero, replace with small number to avoid log(0) error
freq_errors[freq_errors .<= 1e-10] .= 1e-10

 11.970912 seconds (750.81 k allocations: 82.369 MiB, 2.32% gc time, 3.82% compilation time)

plot(freqs, freq_errors, yscale=:log10, xlabel="Max Frequency", ylabel="MSE", ylim=(1e-5, 1e0), label="", linewidth=2, marker=:circle)

Again, we have strongly decaying mean square error, but note that absolute value of the truncation are significantly higher that the values formaximum Pauli weight. We only reach around 1e-5 error with a frequency truncation around 180. This is a very large number and should be combined with other truncations.

More practically, if we know that the parameter range of the PauliRotation gates is small, we can employ the max_sins truncation.

Test Small-Angle Truncation

Here is a loop that additionally passes the radius r for the parameters, i.e., the range around zero where they are uniformly sampled from. We thus compute the average error in a small angle range. For large radii, we expect to need significantly higher max_sins truncation values.

nsins = 0:5:30
rs = 0.1:0.2:1.1
mcsamples = 100_000

pl = plot(yscale=:log10, xlabel="Max Sins", ylabel="MSE", ylim=(1e-5, 1e0))

for r in rs
    @time sins_errors = [estimatemse(circuit, pstr, mcsamples, r; max_sins=ns) for ns in nsins]
    # where errors are zero, replace with small number to avoid log(0) error
    sins_errors[sins_errors .<= 1e-10] .= 1e-10;

    plot!(nsins, sins_errors, label="r=$(round(r, sigdigits=2))", linewidth=2, marker=:circle)
end
pl

  5.583287 seconds (1.01 M allocations: 80.986 MiB, 13.26% compilation time)
  4.408430 seconds (50.74 k allocations: 33.814 MiB, 0.37% gc time)


  4.857387 seconds (50.74 k allocations: 33.814 MiB, 0.16% gc time)
  5.976303 seconds (50.74 k allocations: 33.814 MiB, 0.07% gc time)


  5.196280 seconds (50.74 k allocations: 33.814 MiB, 0.05% gc time)
  6.440426 seconds (50.74 k allocations: 33.814 MiB, 0.05% gc time)

Combining Truncations

To get the most out of Pauli propagation, you probably want to combine truncations. Let us combine max_sins with max_weight.

nsins = 0:5:30
rs = 0.1:0.2:1.1

pl = plot(yscale=:log10, xlabel="Max Sins", ylabel="MSE", ylim=(1e-5, 1e0))

for r in rs
    @time sins_errors = [estimatemse(circuit, pstr, 100_000, r; max_sins=ns, max_weight=6) for ns in nsins]
    # where errors are zero, replace with small number to avoid log(0) error
    sins_errors[sins_errors .<= 1e-10] .= 1e-10;

    plot!(nsins, sins_errors, label="r=$(round(r, sigdigits=2))", linewidth=2, marker=:circle)
end
pl

  6.664841 seconds (1.11 M allocations: 85.749 MiB, 0.29% gc time, 6.37% compilation time)
  6.175398 seconds (50.73 k allocations: 33.814 MiB, 0.11% gc time)


  6.940708 seconds (50.73 k allocations: 33.814 MiB, 0.06% gc time)
 10.239134 seconds (50.73 k allocations: 33.814 MiB, 0.03% gc time)


  7.725097 seconds (50.73 k allocations: 33.814 MiB, 4.45% gc time)
  9.463250 seconds (50.73 k allocations: 33.814 MiB)

The error is pretty much identical. But we in the following we see that the runtime of propagate() is significantly faster with both truncations:

using Random
Random.seed!(42)
thetas = randn(countparameters(circuit))*0.5

With both truncations:

@time psum_both = propagate(circuit, pstr, thetas; max_sins=15, max_weight=6);
print("Number of Paulis: ", length(psum_both))

  2.749529 seconds (1.54 M allocations: 90.100 MiB, 43.77% compilation time)
Number of Paulis: 

35799

With only max_sins truncation:

@time psum_sins = propagate(circuit, pstr, thetas; max_sins=15);
print("Number of Paulis: ", length(psum_sins))

  4.473521 seconds (322.44 k allocations: 52.208 MiB, 0.46% gc time, 6.53% compilation time)
Number of Paulis: 

78430

Significantly slower for a potentially insignificant accuracy improvement.

overlapwithzero(psum_both), overlapwithzero(psum_sins)

(-0.5001034339076172, -0.4982398358741927)

overlapwithzero(psum_both) - overlapwithzero(psum_sins)

-0.0018635980334245628