7. Record custom path properties

A path is what we call the entire past of a propagating Pauli through a quantum circuit. By default we don't record anything about the past, but only store the current Pauli string and its current coefficient. The PathProperties types gives us a wrapper around the coefficients that you can use to record some of the history of your Pauli string. This allows one to learn more about the propagation through the circuit or even truncate on different properties. For example, we use the PauliFreqTracker under the hood to track how often we branch on Pauli rotations, allowing you to truncate via max_freq or max_sins.

using PauliPropagation

One implemented example is the PauliFreqTracker, which tracks how often a Pauli string split at a PauliRotation gate, we call that frequency, and whether it received a cos or sin coefficient.

The PauliFreqTracker

coeff = 1.3
PauliFreqTracker(coeff)

PauliFreqTracker{Float64}(coeff=1.3, nsins=0, ncos=0, freq=0)

nq = 32
PauliString(nq, :X, 16, PauliFreqTracker(coeff))

PauliString(nqubits: 32, PauliFreqTracker(1.3) * IIIIIIIIIIIIIIIXIIII...)

We can either wrap all of our coefficients in this type, or we use the wrapcoefficients() function.

pstr = PauliString(nq, :X, 16, coeff)
wrapped_pstr = wrapcoefficients(pstr, PauliFreqTracker)

PauliString(nqubits: 32, PauliFreqTracker(1.3) * IIIIIIIIIIIIIIIXIIII...)

This wrapping also works for a PauliSum:

psum = PauliSum(nq)
add!(psum, :X, 16, coeff)

PauliSum(nqubits: 32, 1 Pauli term: 
 1.3 * IIIIIIIIIIIIIIIXIIII...
)

wrapped_psum = wrapcoefficients(psum, PauliFreqTracker)

PauliSum(nqubits: 32, 1 Pauli term: 
 PauliFreqTracker(1.3) * IIIIIIIIIIIIIIIXIIII...
)

Let's apply a PauliRotation to it.

pauli_rot = PauliRotation(:Y, 16)

propagate!(pauli_rot, wrapped_psum, 0.4)

PauliSum(nqubits: 32, 2 Pauli terms:
 PauliFreqTracker(1.1974) * IIIIIIIIIIIIIIIXIIII...
 PauliFreqTracker(0.50624) * IIIIIIIIIIIIIIIZIIII...
)

coefficients(wrapped_psum)

ValueIterator for a Dict{UInt64, PauliFreqTracker{Float64}} with 2 entries. Values:
  PauliFreqTracker{Float64}(coeff=1.1973792922037507, nsins=0, ncos=1, freq=1)
  PauliFreqTracker{Float64}(coeff=0.5062438450012458, nsins=1, ncos=0, freq=1)

We see that the term belonging to the IIIIIIIIIIIIIIIXIIII... Pauli string received ncos=1 cosine prefactors, and split a total of freq=1 times. Similarly the other time with nsins=1 sine prefactors. Note that the coeff field of PauliFreqTracker carries the coefficient as in the normal propagation. In fact, this coeff field is very important and we recommend that you have it on your PathProperties type. If you just want to use max_freq and max_sin truncation, then you don't have to manually wrap the coefficients. We will do that automatically and then unwrap the coefficients again:

# this truncates any Pauli with any sine coefficient
propagate(pauli_rot, psum, 0.4; max_sins=0)


## but the in-place propagate!() function will not automatically convert
# propagate!(pauli_rot, psum, 0.4; max_sins=0)

PauliSum(nqubits: 32, 1 Pauli term: 
 1.1974 * IIIIIIIIIIIIIIIXIIII...
)

A custom PathProperties type

Now let us define a custom PathProperties type. For example, we can define a type that records how often a Clifford gate reduced the weight (number of non-Identity Paulis) of a propagating Pauli String, how often it left it unchanged, and how often it increased it.

struct CliffordWeightTracker <: PathProperties
    coeff::Float64  # for simplicity use Float64, but in general you might want a parametrized number type
    reduced::Int
    unchanged::Int
    increased::Int
end

Now let's try to wrap the coefficient of a PauliString in this type.

## this will throw an error if you run it
# wrapcoefficients(pstr, CliffordWeightTracker)

This did not work, because if you want wrapcoefficients() to do the job, you need to define a constructor for your type expecting only the coefficient. You can still manually construct the wrapped Pauli string by filling in the fields that aren't coeff with zeros, or you defined the default constructor.

CliffordWeightTracker(coeff) = CliffordWeightTracker(coeff, 0, 0, 0)

CliffordWeightTracker

cwt_pstr = wrapcoefficients(pstr, CliffordWeightTracker)

PauliString(nqubits: 32, CliffordWeightTracker(1.3) * IIIIIIIIIIIIIIIXIIII...)

And propagation already works!

cliff_gate = CliffordGate(:H, 16)

result_psum = propagate(cliff_gate, cwt_pstr)

PauliSum(nqubits: 32, 1 Pauli term: 
 CliffordWeightTracker(1.3) * IIIIIIIIIIIIIIIZIIII...
)

coefficients(result_psum)

ValueIterator for a Dict{UInt64, CliffordWeightTracker} with 1 entry. Values:
  CliffordWeightTracker(coeff=1.3, reduced=0, unchanged=0, increased=0)

The propagation natively worked because we have a coeff field on our custom type (you can restart the notebook and change the name of the coeff field to see where it fails). But now we need to define how CliffordWeightTracker increments its fields, i.e., how it knows what a Clifford gate does to it.

The way that CliffordGates are defined is via the applytoall!() and apply() functions. applytoall!() loads the corresponding look-up map for the gate and passes it to the specific apply(). See our example notebook on custom gates for more information. We now need to define a version of apply() that is slightly more specifically typed so that when propagating with the CliffordWeightTracker, the propagation will dispatch to that code. For reference, the following is how the default code looks like:

function PauliPropagation.apply(gate::CliffordGate, pstr, coeff, lookup_map; kwargs...)
    # the lookup array carries the new Paulis + sign for every occuring old Pauli combination

    qinds = gate.qinds

    # this integer carries the active Paulis on its bits
    lookup_int = getpauli(pstr, qinds)

    # this integer can be used to index into the array returning the new Paulis
    # +1 because Julia is 1-indexed and lookup_int is 0-indexed
    partial_pstr, sign = lookup_map[lookup_int+1]

    # insert the bits of the new Pauli into the old Pauli
    pstr = setpauli(pstr, partial_pstr, qinds)

    coeff *= sign

    # always a length-1 tuple, which will be compiled away
    return ((pstr, coeff),)
end

We will still use this apply() function, but wrap it with specific logic for CliffordWeightTracker.

@inline function PauliPropagation.apply(gate::CliffordGate, pstr, wrappedcoeff::CliffordWeightTracker, lookup_map; kwargs...)
    
    # normal Clifford logic here
    new_pstr, new_coeff = only(apply(gate, pstr, wrappedcoeff.coeff, lookup_map; kwargs...))
    
    old_weight = countweight(pstr)
    new_weight = countweight(new_pstr)
    
    if new_weight < old_weight
        new_wrappedcoeff = CliffordWeightTracker(new_coeff, wrappedcoeff.reduced+1, wrappedcoeff.unchanged, wrappedcoeff.increased)
    elseif new_weight == old_weight
        new_wrappedcoeff = CliffordWeightTracker(new_coeff, wrappedcoeff.reduced, wrappedcoeff.unchanged+1, wrappedcoeff.increased)
    else
        new_wrappedcoeff = CliffordWeightTracker(new_coeff, wrappedcoeff.reduced, wrappedcoeff.unchanged, wrappedcoeff.increased+1)
    end
    
    # same return format as the Clifford apply() function
    return ((new_pstr, new_wrappedcoeff),)
end

cliff_gate = CliffordGate(:H, 16)

print("Before: ", cwt_pstr)

result_psum = propagate(cliff_gate, cwt_pstr)

Before: PauliString(nqubits: 32, CliffordWeightTracker(1.3) * IIIIIIIIIIIIIIIXIIII...)




PauliSum(nqubits: 32, 1 Pauli term: 
 CliffordWeightTracker(1.3) * IIIIIIIIIIIIIIIZIIII...
)

coefficients(result_psum)

ValueIterator for a Dict{UInt64, CliffordWeightTracker} with 1 entry. Values:
  CliffordWeightTracker(coeff=1.3, reduced=0, unchanged=1, increased=0)

And now increasing weight,

cliff_gate = CliffordGate(:CNOT, [16, 17])

print("Before: ", cwt_pstr)

result_psum2 = propagate(cliff_gate, cwt_pstr)

Before: PauliString(nqubits: 32, CliffordWeightTracker(1.3) * IIIIIIIIIIIIIIIXIIII...)




PauliSum(nqubits: 32, 1 Pauli term: 
 CliffordWeightTracker(1.3) * IIIIIIIIIIIIIIIXXIII...
)

coefficients(result_psum2)

ValueIterator for a Dict{UInt64, CliffordWeightTracker} with 1 entry. Values:
  CliffordWeightTracker(coeff=1.3, reduced=0, unchanged=0, increased=1)

It seems to work!

Let's run an actual circuit.

# 3 layers of the `efficientsu2` circuit, which contains CNOT gates
circuit = efficientsu2circuit(nq, 3);
thetas = randn(countparameters(circuit))

@time result_psum = propagate(circuit, cwt_pstr, thetas)

  1.068854 seconds (385.16 k allocations: 53.057 MiB, 2.56% gc time, 69.02% compilation time)





PauliSum(nqubits: 32, 130708 Pauli terms:
 CliffordWeightTracker(0.00022377) * IIIIIIIIIIIIIIYXIYIY...
 CliffordWeightTracker(0.0034411) * IIIIIIIIIIIIXZIXXYIZ...
 CliffordWeightTracker(-0.00068928) * IIIIIIIIIIIIIIZYIIXY...
 CliffordWeightTracker(-2.2387e-6) * IIIIIIIIIIIIZXXYXZZX...
 CliffordWeightTracker(-0.00021566) * IIIIIIIIIIIIZZXZXXXX...
 CliffordWeightTracker(2.9664e-7) * IIIIIIIIIIIIZXZYIXXY...
 CliffordWeightTracker(0.00082356) * IIIIIIIIIIIIZZXIYXYX...
 CliffordWeightTracker(-3.7748e-6) * IIIIIIIIIIIIXZXYIXXZ...
 CliffordWeightTracker(0.0005538) * IIIIIIIIIIIIXZZZYXZX...
 CliffordWeightTracker(0.0020224) * IIIIIIIIIIIIZZIZIYXY...
 CliffordWeightTracker(2.3089e-6) * IIIIIIIIIIIIXXXZYZZZ...
 CliffordWeightTracker(-2.3351e-6) * IIIIIIIIIIIIZXZIYYZI...
 CliffordWeightTracker(2.187e-5) * IIIIIIIIIIIIIIZZXZXX...
 CliffordWeightTracker(2.9883e-5) * IIIIIIIIIIIIZXXXXIXX...
 CliffordWeightTracker(-1.6147e-6) * IIIIIIIIIIIIZZXIXZXI...
 CliffordWeightTracker(-2.6473e-5) * IIIIIIIIIIIIZZZIIZXY...
 CliffordWeightTracker(0.00015148) * IIIIIIIIIIIIXXZXXIYZ...
 CliffordWeightTracker(0.00011416) * IIIIIIIIIIIIIIXYYIXX...
 CliffordWeightTracker(-2.2066e-5) * IIIIIIIIIIIIIIXIXYXZ...
 CliffordWeightTracker(-8.5667e-6) * IIIIIIIIIIIIZXZZIXYI...
  ⋮)

Make a nice plot from it:

reduced = 0
unchanged = 0
increased = 0
for cliff_weight_tracker in coefficients(result_psum)
    reduced += cliff_weight_tracker.reduced
    unchanged += cliff_weight_tracker.unchanged
    increased += cliff_weight_tracker.increased
end

# calculate the average
n_terms = length(result_psum)
reduced /= n_terms
unchanged /= n_terms
increased /= n_terms

using Plots

bar([reduced, unchanged, increased], label="", xticks=([1, 2, 3], ["reduced", "unchanged", "increased"]), ylabel="Count")

I worked! But a few notes on the result:

First, why is "unchanged" dominating? Because we start with a local Z Pauli with Identity on all but one site. This means that the behavior is dominated by commutations, i.e., CNOT gates that don't contribute at all. We can quickly fix this:

@inline function PauliPropagation.apply(gate::CliffordGate, pstr, path_prop::CliffordWeightTracker, lookup_map; kwargs...)

    # normal Clifford logic here
    new_pstr, new_coeff = only(apply(gate, pstr, path_prop.coeff, lookup_map; kwargs...))

    if !(new_pstr==pstr && new_coeff==path_prop.coeff)  # only track when it didn't fully commute
        old_weight = countweight(pstr)
        new_weight = countweight(new_pstr)

        if new_weight < old_weight
            path_prop = CliffordWeightTracker(new_coeff, path_prop.reduced+1, path_prop.unchanged, path_prop.increased)
        elseif new_weight == old_weight
            path_prop = CliffordWeightTracker(new_coeff, path_prop.reduced, path_prop.unchanged+1, path_prop.increased)
        else
            path_prop = CliffordWeightTracker(new_coeff, path_prop.reduced, path_prop.unchanged, path_prop.increased+1)
        end
    end

    # same return format as the Clifford apply() function
    return ((new_pstr, path_prop),)
end

@time result_psum = propagate(circuit, cwt_pstr, thetas)

  0.921216 seconds (52.77 k allocations: 36.716 MiB, 48.65% gc time, 22.31% compilation time: 100% of which was recompilation)

Plot it up:

reduced = 0
unchanged = 0
increased = 0
for cliff_weight_tracker in coefficients(result_psum)
    reduced += cliff_weight_tracker.reduced
    unchanged += cliff_weight_tracker.unchanged
    increased += cliff_weight_tracker.increased
end

# calculate the average
n_terms = length(result_psum)
reduced /= n_terms
unchanged /= n_terms
increased /= n_terms

bar([reduced, unchanged, increased], label="", xticks=([1, 2, 3], ["reduced", "unchanged", "increased"]), ylabel="Count", title="Removing commutation")

Great, this is fixed. But now, why are the counts so low? This brings us to the second point:

-> By default, merging paths calls + between two PathProperties objects. We have defined + so that it adds the coeff fields together, but takes the minimum of the other fields over the two summands. This current behavior makes sense if you want to truncate based on your recorded properties, and higher values are associated with higher computational resources. This is the case, for example, for the PauliFreqTracker with the max_freq and max_sins truncations. After merging, both objects propagate at the cost of one, so we choose the one that is furthest from being truncated. The huge caveat is that this assumes that incrementing the fields brings them closer to being truncated.

Let's for example define a + operation for CliffordWeightTracker that takes the maximum, just to show how it would work.

function PauliPropagation.:+(path_prop1::CliffordWeightTracker, path_prop2::CliffordWeightTracker)
    return CliffordWeightTracker(
        path_prop1.coeff + path_prop2.coeff,
        max(path_prop1.reduced, path_prop2.reduced),
        max(path_prop1.unchanged, path_prop2.unchanged),
        max(path_prop1.increased, path_prop2.increased)
    )
end

@time result_psum = propagate(circuit, cwt_pstr, thetas)

  0.603841 seconds (194.25 k allocations: 43.348 MiB, 71.07% compilation time: 100% of which was recompilation)

Plot it up:

reduced = 0
unchanged = 0
increased = 0
for cliff_weight_tracker in coefficients(result_psum)
    reduced += cliff_weight_tracker.reduced
    unchanged += cliff_weight_tracker.unchanged
    increased += cliff_weight_tracker.increased
end

# calculate the average
n_terms = length(result_psum)
reduced /= n_terms
unchanged /= n_terms
increased /= n_terms

bar([reduced, unchanged, increased], label="", xticks=([1, 2, 3], ["reduced", "unchanged", "increased"]), ylabel="Count", title="Custom `+` function")

This now tells you that there exist Pauli strings that had their weight increased 10 times, but none that had their weight reduced 10 times (or even 4 times). Enjoy recording custom path properties!