GitHub

5. How to Define Custom Gates

PauliPropagation.jl is extensible and allows you to define your own gates. Depending on how much you can or want to code, you can definte a gate that works or one that is as fast as it gets. Here will see what you need to define.

Just before, we saw the TransferMapGate and how it can be defined from a matrix representation. Such a definition works and may already be enough for your purposes. But it is not guaranteed to give you the best performance. In this notebook, we will discover more of the inner workings and how to go towards optimized gates.

using PauliPropagation

A static gate that maps one Pauli string to one Pauli string

Let us start by defining a SWAP gate. It is sub-typing from StaticGate, which denotes that it does not take any variable parameters at propagation time. It always acts the same. Types in Julia are typically just name tags carrying data. In this case, you need to give the gate a name that doesn't already exist in your name space, and think about what data it needs to carry: the qubit indices that it swaps. Note that everything about this can be customized as long as you define the action accordingly.

struct CustomSWAPGate <: StaticGate
    qinds::Tuple{Int, Int}  # The two sites to be swapped
end

The action of a SWAP gate on a Pauli string is that it swaps the Paulis on two sites. We can now define a function apply which receives these 3 arguments in this order, as well as potential kwargs: apply(gate::YourGate, pstr, coeff; kwargs...). We can ignore kwargs for now, but you can use them to pass arguments from the top level down to your function. If your custom gate is parametrized, then you have you sub-type ParametrizedGate instead of StaticGate, and the apply function will receive a fourth argument, which is the parameter: apply(gate::YourGate, pstr, coeff, theta; kwargs...). This is how you can define SWAP:

function PauliPropagation.apply(gate::CustomSWAPGate, pstr, coeff; kwargs...)
    # get the Pauli on the first site
    pauli1 = getpauli(pstr, gate.qinds[1])
    # get the Pauli on the second site
    pauli2 = getpauli(pstr, gate.qinds[2])
    
    # set the Pauli on the first site to the second Pauli
    pstr = setpauli(pstr, pauli2, gate.qinds[1])
    # set the Pauli on the second site to the first Pauli
    pstr = setpauli(pstr, pauli1, gate.qinds[2])

    # apply() is always expected to return a tuple of (pstr, coeff) tuples
    return tuple((pstr, coeff))
end

This is it, really.

apply() is expected to always return a tuple of (pstr, coeff) tuples. Here, only one Pauli string and its coefficient are returned, so we need to create a tuple of length one via tuple((pstr, coeff)). Alternatively ((pstr, coeff),) also works. We will see later what to do with gates that create new Pauli strings. Try the SWAP gate on a few very simple examples:

# define the gate
g = CustomSWAPGate((1, 2))

# define a Pauli string
pstr = PauliString(2, :X, 1)
# turn into a PauliSum so it is easier to compare
psum = PauliSum(pstr)
println("Paulis before: ", psum)

swapped_psum = propagate(g, psum)
println("Paulis after: ", swapped_psum)
Paulis before: PauliSum(nqubits: 2, 1 Pauli term: 
 1.0 * XI
)


Paulis after: 

PauliSum(nqubits: 2, 1 Pauli term: 
 1.0 * IX
)
nq = 3

# define the gate
g = CustomSWAPGate((2, 3))

# add some initial Pauli strings
psum = PauliSum(nq) 
add!(psum, :X, 1, 0.1)
add!(psum, :X, 2, 0.2)
add!(psum, [:Z, :Z], [1, 2], 0.3)
add!(psum, [:Y, :X], [2, 3], 0.4)
println("Paulis before: ", psum)

swapped_psum = propagate(g, psum)
println("Paulis after: ", swapped_psum)
Paulis before: PauliSum(nqubits: 3, 4 Pauli terms:
 0.2 * IXI
 0.3 * ZZI
 0.4 * IYX
 0.1 * XII
)


Paulis after: PauliSum(nqubits: 3, 4 Pauli terms:
 0.3 * ZIZ
 0.2 * IIX
 0.4 * IXY
 0.1 * XII
)

This works! You can now freely plug such a gate into a circuit and let it swap your qubits.

Consider a simple more wholistic example where we insert a layer of SWAP gates at the beginning of the circuit, which means at the end of the backpropagation.

nq = 6

# some initial observable
pstr = PauliString(nq, [:X, :X], [3, 4])

# a 1D bricklayer topology
topology = bricklayertopology(nq)

# an empty circuit
circuit = Gate[]
# an RX Pauli rotation layer
rxlayer!(circuit, nq)
# an RZ Pauli rotation layer
rzlayer!(circuit, nq)
# an RZZ Pauli rotation layer on the topology 
rzzlayer!(circuit, topology);


# define some parameters
using Random
Random.seed!(420)
nparams = countparameters(circuit)
thetas = randn(nparams)

Now define our custom layer of SWAP gates:

ourSWAPs = Gate[]
# insert swaps flipping the order of qubits 4 to 6
# 4 5 6 -> 5 4 6
push!(ourSWAPs, CustomSWAPGate((4, 5)))
# 5 4 6 -> 5 6 4
push!(ourSWAPs, CustomSWAPGate((5, 6)))
# 5 6 4 -> 6 5 4
push!(ourSWAPs, CustomSWAPGate((4, 5)))
# this flattens the two circuits into one
ourSWAP_circuit = [ourSWAPs..., circuit...]
20-element Vector{Gate}:
 CustomSWAPGate((4, 5))
 CustomSWAPGate((5, 6))
 CustomSWAPGate((4, 5))
 PauliRotation([:X], [1])
 PauliRotation([:X], [2])
 PauliRotation([:X], [3])
 PauliRotation([:X], [4])
 PauliRotation([:X], [5])
 PauliRotation([:X], [6])
 PauliRotation([:Z], [1])
 PauliRotation([:Z], [2])
 PauliRotation([:Z], [3])
 PauliRotation([:Z], [4])
 PauliRotation([:Z], [5])
 PauliRotation([:Z], [6])
 PauliRotation([:Z, :Z], [1, 2])
 PauliRotation([:Z, :Z], [3, 4])
 PauliRotation([:Z, :Z], [5, 6])
 PauliRotation([:Z, :Z], [2, 3])
 PauliRotation([:Z, :Z], [4, 5])

Propagate through the circuit and overlap with the zero state:

our_psum = propagate(ourSWAP_circuit, pstr, thetas)
PauliSum(nqubits: 6, 81 Pauli terms:
 0.0056441 * IYYIZX
 -0.18621 * IZZIIY
 -0.093178 * IIXIYY
 0.021738 * IIZIZY
 0.012669 * IYXIZX
 -0.077652 * IIXIZY
 -0.022844 * IZXIZZ
 0.095416 * IYXIIX
 -0.056601 * IIZIIY
 -0.060448 * IYZIYY
 0.016621 * IYYIIY
 0.010874 * IZYIYY
 0.51708 * IIXIIX
 0.087206 * IIXIZZ
 0.0090624 * IZYIZY
 -0.1473 * IYZIIZ
 -0.027412 * IZXIYZ
 0.071515 * IZZIZY
 -0.050376 * IYZIZY
 -0.0071724 * IIYIIY
  ⋮)
overlapwithzero(our_psum)
0.16795162488187168

This looks okay, but is it correct? One thing you may have noticed is that SWAP is a Clifford operation, i.e., one that takes one Pauli to exactly one other Pauli. We actually have that in our package so we can easily compare.

cliffordSWAPs = Gate[]
# insert swaps flipping the order of qubits 4 to 6
# 4 5 6 -> 5 4 6
push!(cliffordSWAPs, CliffordGate(:SWAP, (4, 5)))
# 5 4 6 -> 5 6 4
push!(cliffordSWAPs, CliffordGate(:SWAP, (5, 6)))
# 5 6 4 -> 6 5 4
push!(cliffordSWAPs, CliffordGate(:SWAP, (4, 5)))
# this flattens the two circuits into one
cliffordSWAP_circuit = [cliffordSWAPs..., circuit...]
20-element Vector{Gate}:
 CliffordGate(:SWAP, [4, 5])
 CliffordGate(:SWAP, [5, 6])
 CliffordGate(:SWAP, [4, 5])
 PauliRotation([:X], [1])
 PauliRotation([:X], [2])
 PauliRotation([:X], [3])
 PauliRotation([:X], [4])
 PauliRotation([:X], [5])
 PauliRotation([:X], [6])
 PauliRotation([:Z], [1])
 PauliRotation([:Z], [2])
 PauliRotation([:Z], [3])
 PauliRotation([:Z], [4])
 PauliRotation([:Z], [5])
 PauliRotation([:Z], [6])
 PauliRotation([:Z, :Z], [1, 2])
 PauliRotation([:Z, :Z], [3, 4])
 PauliRotation([:Z, :Z], [5, 6])
 PauliRotation([:Z, :Z], [2, 3])
 PauliRotation([:Z, :Z], [4, 5])
clifford_psum = propagate(cliffordSWAP_circuit, pstr, thetas)
PauliSum(nqubits: 6, 81 Pauli terms:
 0.0056441 * IYYIZX
 -0.18621 * IZZIIY
 -0.093178 * IIXIYY
 0.021738 * IIZIZY
 0.012669 * IYXIZX
 -0.077652 * IIXIZY
 -0.022844 * IZXIZZ
 0.095416 * IYXIIX
 -0.056601 * IIZIIY
 -0.060448 * IYZIYY
 0.016621 * IYYIIY
 0.010874 * IZYIYY
 0.51708 * IIXIIX
 0.087206 * IIXIZZ
 0.0090624 * IZYIZY
 -0.1473 * IYZIIZ
 -0.027412 * IZXIYZ
 0.071515 * IZZIZY
 -0.050376 * IYZIZY
 -0.0071724 * IIYIIY
  ⋮)

And these produce the same Pauli sum and consequently the same expectation!

our_psum == clifford_psum
true
overlapwithzero(clifford_psum)
0.16795162488187168

A parametrized custom noise channel

The above example was not only a Clifford gate mapping one Pauli string to one Pauli string, but it was also not parametrized. We will now implement a funky custom and parametrized noise gate. Note that it may technically not be a valid quantum channel, but we can still define it and have it be part of the simulation.

struct FunkyNoise <: ParametrizedGate
    qinds::Tuple{Int, Int}  # say it acts on 2 qubits
end

Say we want to define a noise gate with strength p that behaves in the following way:

XX -> (1 - 2p) * XX + p * IX + p * XI
YY -> (1 - 2p) * YY + p * IY + p * YI
ZZ -> (1 - 2p) * ZZ + p * IZ + p * ZI
P -> P else

Let's now implement that.

# if your gate is a subtype of `ParametrizedGate`, `apply()` receives a fourth argument which is the parameter
function PauliPropagation.apply(gate::FunkyNoise, pstr, coeff, p; kwargs...)
    
    # the integer representation of the Paulis sitting on the active sites
    paulis = getpauli(pstr, gate.qinds)

    # check whether the Paulis in `paulis` are the integer representations of our targets
    # "|" is the OR operator
    # if you replaced the (:X, :X), (:Y, :Y), or (:Z, :Z) by their integer representation,
    # i.e., 5, 10, or 15, respectively, it would be faster.
    # Use `symboltoint()` to find out those integers.
    # if ispauli(paulis, 5) | ispauli(paulis, 10) | ispauli(paulis, 15)
    if ispauli(paulis, (:X, :X)) | ispauli(paulis, (:Y, :Y)) | ispauli(paulis, (:Z, :Z))
        
        # calculate the coefficients
        coeff_PP = coeff * (1 - 2p)
        coeff_IP = coeff_PI = coeff * p        

        # create new Paulis that have I in the right position
        # if you replace :I with the integer 0, it will be faster
        pstr_IP = setpauli(pstr, :I, gate.qinds[1])
        pstr_PI = setpauli(pstr, :I, gate.qinds[2])

        return tuple((pstr, coeff_PP), (pstr_IP, coeff_IP), (pstr_PI, coeff_PI))
    end

    # else do nothing and return the original pstr and coefficient pair
    return tuple((pstr, coeff))
end

It was a bit more challenging than above, but we are done!

nq = 2
g = FunkyNoise((1, 2))
p = 0.1
0.1
# this branches
pstr = PauliString(nq, [:X, :X], [1, 2])

# turn into a PauliSum just for consistent display
psum = PauliSum(pstr)
println("Pauli sum before: ", psum)

funky_noise_psum = propagate(g, pstr, p)
println("Pauli sum after: ", funky_noise_psum)
Pauli sum before: PauliSum(nqubits: 2, 1 Pauli term: 
 1.0 * XX
)
Pauli sum after: 

PauliSum(nqubits: 2, 3 Pauli terms:
 0.8 * XX
 0.1 * IX
 0.1 * XI
)
# and this branches
pstr = PauliString(nq, [:Y, :Y], [1, 2])

# turn into a PauliSum for consistent display
psum = PauliSum(pstr)
println("Pauli sum before: ", psum)

funky_noise_psum = propagate(g, pstr, p)
println("Pauli sum after: ", funky_noise_psum)
Pauli sum before: PauliSum(nqubits: 2, 1 Pauli term: 
 1.0 * YY
)
Pauli sum after: PauliSum(nqubits: 2, 3 Pauli terms:
 0.1 * YI
 0.8 * YY
 0.1 * IY
)
# and this branches
pstr = PauliString(nq, [:Z, :Z], [1, 2])

# turn into a PauliSum for consistent display
psum = PauliSum(pstr)
println("Pauli sum before: ", psum)

funky_noise_psum = propagate(g, pstr, p)
println("Pauli sum after: ", funky_noise_psum)
Pauli sum before: PauliSum(nqubits: 2, 1 Pauli term: 
 1.0 * ZZ
)
Pauli sum after: PauliSum(nqubits: 2, 3 Pauli terms:
 0.8 * ZZ
 0.1 * IZ
 0.1 * ZI
)
# but not this
pstr = PauliString(nq, [:Y, :Z], [1, 2])

# turn into a PauliSum for consistent display
psum = PauliSum(pstr)
println("Pauli sum before: ", psum)

funky_noise_psum = propagate(g, pstr, p)
println("Pauli sum after: ", funky_noise_psum)
Pauli sum before: PauliSum(nqubits: 2, 1 Pauli term: 
 1.0 * YZ
)
Pauli sum after: PauliSum(nqubits: 2, 1 Pauli term: 
 1.0 * YZ
)
# or this
pstr = PauliString(nq, [:I, :X], [1, 2])

# turn into a PauliSum for consistent display
psum = PauliSum(pstr)
println("Pauli sum before: ", psum)

funky_noise_psum = propagate(g, pstr, p)
println("Pauli sum after: ", funky_noise_psum)
Pauli sum before: PauliSum(nqubits: 2, 1 Pauli term: 
 1.0 * IX
)
Pauli sum after: PauliSum(nqubits: 2, 1 Pauli term: 
 1.0 * IX
)
# or this
pstr = PauliString(nq, [:I, :I], [1, 2])

# turn into a PauliSum for consistent display
psum = PauliSum(pstr)
println("Pauli sum before: ", psum)

funky_noise_psum = propagate(g, pstr, p)
println("Pauli sum after: ", funky_noise_psum)
Pauli sum before: PauliSum(nqubits: 2, 1 Pauli term: 
 1.0 * II
)
Pauli sum after: PauliSum(nqubits: 2, 1 Pauli term: 
 1.0 * II
)

A noise channel with several parameters

You can also define gates that take several parameters. The way this can be done is by inserting a vector of tuple of values into the parameter vector. For example, thetas = [0.1, [0.3, -1.0], 2.1] would work and assume that the second gate takes two parameters. However, some functions like countparameters() may not work yet with such a case. Here we define our new custom gate: A depolarizing noise channel with flexible noise strengths per X, Y, or Z Pauli.

struct CustomDepolNoise <: ParametrizedGate
    qind::Int  # say it acts on 1 qubit - i.e. a single qubit gate takes a single integer as input 
end

Here we just specify which site it acts on. Now onto the apply() function:

function PauliPropagation.apply(gate::CustomDepolNoise, pstr, coeff, paramvec; kwargs...)    
    # note that to take multiple parameters as an input you can just stack them in a single vector input - 'paramvec' here   
    
    # the integer representation of the Pauli sitting on the active site
    pauli = getpauli(pstr, gate.qind)
    
    # Here we use that :X is represented by the integer 1, :Y by 2 and :Z by 3
    # update the coefficient if applicable, the Pauli string remains unchanged
    if ispauli(pauli, 1) 
        coeff = coeff * (1 - paramvec[1])  
    elseif ispauli(pauli, 2) 
        coeff = coeff * (1 - paramvec[2])       
    elseif ispauli(pauli, 3) 
        coeff = coeff * (1 - paramvec[3])       
    end

    # to be super fancy, we could have returned early if pauli == 0
    # and otherwise index into paramvec by the pauli itself, like paramvec[pauli]
    
    return tuple((pstr, coeff))
end
# set up a single-qubit example
nq = 1

g =  CustomDepolNoise(1)

px = 0.1
py = 0.2
pz = 0.3

pvec = [px, py, pz]
3-element Vector{Float64}:
 0.1
 0.2
 0.3
for pauli in [:I, :X, :Y, :Z]
    pstr = PauliString(nq, pauli, 1)
    
    # turn into a PauliSum just for consistent display
    psum = PauliSum(pstr)
    println("Paulis before: ", psum)
    
    depolnoise_psum = propagate(g, psum, pvec)
    println("Paulis after: ", depolnoise_psum)
    # create some space
    println("\n")
end
Paulis before: PauliSum(nqubits: 1, 1 Pauli term: 
 1.0 * I
)
Paulis after: 

PauliSum(nqubits: 1, 1 Pauli term: 
 1.0 * I
)


Paulis before: PauliSum(nqubits: 1, 1 Pauli term: 
 1.0 * X
)
Paulis after: PauliSum(nqubits: 1, 1 Pauli term: 
 0.9 * X
)


Paulis before: PauliSum(nqubits: 1, 1 Pauli term: 
 1.0 * Y
)
Paulis after: PauliSum(nqubits: 1, 1 Pauli term: 
 0.8 * Y
)


Paulis before: PauliSum(nqubits: 1, 1 Pauli term: 
 1.0 * Z
)
Paulis after: PauliSum(nqubits: 1, 1 Pauli term: 
 0.7 * Z
)

This is exactly what we expect, and we could even verify those results with some damping operations we have. Note that they are not valid quantum channels!

gx =  PauliPropagation.PauliXDamping(1)
gy =  PauliPropagation.PauliYDamping(1)
gz =  PauliPropagation.PauliZDamping(1)

circuit = Gate[gx, gy, gz]

for pauli in [:I, :X, :Y, :Z]
    pstr = PauliString(nq, pauli, 1)
    
    # turn into a PauliSum just for consistent display
    psum = PauliSum(pstr)
    println("Paulis before: ", psum)
    
    depolnoise_psum = propagate(circuit, psum, pvec)
    println("Paulis after: ", depolnoise_psum)
    # create some space
    println("\n")
end
Paulis before: PauliSum(nqubits: 1, 1 Pauli term: 
 1.0 * I
)
Paulis after: 

PauliSum(nqubits: 1, 1 Pauli term: 
 1.0 * I
)


Paulis before: PauliSum(nqubits: 1, 1 Pauli term: 
 1.0 * X
)
Paulis after: PauliSum(nqubits: 1, 1 Pauli term: 
 0.9 * X
)


Paulis before: PauliSum(nqubits: 1, 1 Pauli term: 
 1.0 * Y
)
Paulis after: PauliSum(nqubits: 1, 1 Pauli term: 
 0.8 * Y
)


Paulis before: PauliSum(nqubits: 1, 1 Pauli term: 
 1.0 * Z
)
Paulis after: PauliSum(nqubits: 1, 1 Pauli term: 
 0.7 * Z
)

One important thing to keep in mind when using multi-parameter gates:

If you use them within circuits, some parameters will be numbers, others will be vectors of numbers. This means you cannot easily work with functions like countparameters(circuit) or randn(nparams). countparameters(circuit) will return the correct length of the theta vector but not know whether gates use several parameters. And randn(nparams) will create a vector with Float64 entries. Once defined and typed, you cannot change the entries to Vector{Float64}. See the following cautionary examples:

circuit = [PauliRotation(:X, 1), CustomDepolNoise(1)]
# returns 2 even though CustomDepolNoise takes 3 parameters
# but 2 is the correct length of the full parameter vector
nparams = countparameters(circuit)
2
thetas = randn(nparams)
@show typeof(thetas)

# this doesn't work because thetas is already strictly typed to Float64 entries
# thetas[2] = pvec
typeof(thetas) = Vector{Float64}





Vector{Float64} (alias for Array{Float64, 1})
# this works
# for type-stability consider typing it a bit more concretely like `Union{Float64, Vector{Float64}}[randn(), pvec]`
thetas = [randn(), pvec]
2-element Vector{Any}:
 1.128502633232417
  [0.1, 0.2, 0.3]

This is how you can implement custom gates - by defining their action on Pauli strings (in the Heisenberg picture).

The branching noise gate, for example, could be implemented more efficiently via some lower level details, but this certainly works. This way of defining apply() for a gate that either branches into three or remains as one Pauli string is also usually not optimal because the compiler cannot infer from the input types whether the output will be a length 3 tuple or a length 1 tuple. This is related to type-stability (some Julia docs) and a topic of a more advanced example notebook.