Adding Symmetry to Pauli Propagation

using PauliPropagation
using PauliPropagation: _translatetolowestinteger, inttostring

In this notebook, we show how symmetry can be added in Pauli Propagation. The basic idea is to merge Pauli strings if they are related by a symmetry.

The related publication is Leveraging Symmetry Merging in Pauli Propagation.

Lowest integer representation (canonical representation)

For instance, consider the following operator:

\[ P = XII + IXI + IIX. \]

Clearly, the Pauli string $P$ is translational invariant. Suppose we can just work with a representative $XII$, then the effective Pauli string is:

\[ \tilde{P} = 3 XII. \]

# Let us define the initial Pauli string
nq = 3
p = PauliSum(nq)
for i in 1:nq
    add!(p, :X, i)
end
println("Initial Pauli string: ", p)

Initial Pauli string: PauliSum(nqubits: 3, 3 Pauli terms:
 1.0 * IXI
 1.0 * IIX
 1.0 * XII
)

First we discuss how to find the representatives.

Consider the translation symmetru in $1d$, we adopt representatives that are given by the lowest integers.

# To find the representative, we define a function that maps a Pauli string to an integer.
# Consider the second Pauli string:
pstr = PauliString(nq, :X, 2)
# The integer representation is given by:
println("Integer representation of ", pstr, " is ", pstr.term)

# We can then translate the integer to the lowest integer in its symmetry class:
rep_int = _translatetolowestinteger(pstr.term, nq)
println("The representative integer is ", rep_int)

# We can verify that this lowest integer representation corresponds to XII
println("The representative Pauli string is XII: ", PauliString(nq, :X, 1).term)

Integer representation of PauliString(nqubits: 3, 1.0 * IXI)

 is 4
The representative integer is 

1
The representative Pauli string is XII: 1

# To get the effective representation of the original Pauli string, we can do:
translationmerge(p)

PauliSum(nqubits: 3, 1 Pauli term: 
 3.0 * XII
)

As another example, consider the two Pauli strings $XIZI$ and $ZIXI$. Their representative is given by the latter because we want the smaller Pauli to be in the higher position and identities as far as possible.

# Let's verify that is the case:
nq = 4
pstr = PauliString(nq, [:X, :Z], [2, 4])
println("The integer representation of ", pstr, " is ", pstr.term)
println("The lowest integer representation of ", pstr, " is ", _translatetolowestinteger(pstr.term, nq))
pstr = PauliString(nq, [:Z, :X], [1, 3])
println("The integer representation of ", pstr, " is ", pstr.term)  
# println("The lowest integer representation of ", pstr, " is ", _translatetolowestinteger(pstr.term, nq))

The integer representation of PauliString(nqubits: 4, 1.0 * IXIZ)

 is 196
The lowest integer representation of PauliString(nqubits: 4, 1.0 * IXIZ) is 19
The integer representation of PauliString(nqubits: 4, 1.0 * ZIXI) is 19

We can also apply this in $2d$. It's a bit tricky to visualize a $2d$ Pauli string but we can use the helper function

# Imagine we have a 2x3 lattice with periodic boundary conditions.
# We can represent a Pauli string on this lattice using a 6-qubit Pauli string.
# For example, consider the following Pauli string:
nx, ny = 2, 3
nq = nx * ny
pstr = PauliString(nq, [:X, :Y, :Z], [1, 3, 5])
println("The integer representation of ", pstr, " is ", pstr.term)
# To find its representative under 2D translations, we can use the following helper function:
println("The 2d visualization of ", pstr, " is:" )
println(inttostring(pstr.term, nx, ny))

The integer representation of PauliString(nqubits: 6, 1.0 * XIYIZI)

 is 801
The 2d visualization of PauliString(nqubits: 6, 1.0 * XIYIZI) is:
XI
YI
ZI

# Now consider another Pauli string that is a translation of the previous one:
pstr2 = PauliString(nq, [:X, :Y, :Z], [2, 4, 6])
println("The integer representation of ", pstr2, " is ", pstr2.term)
println("The 2d visualization of ", pstr2, " is:" )
println(inttostring(pstr2.term, nx, ny))

pstr3 = PauliString(nq, [:Z, :X, :Y], [4, 6, 2])
println("The integer representation of ", pstr3, " is ", pstr3.term)
println("The 2d visualization of ", pstr3, " is:" )
println(inttostring(pstr3.term, nx, ny))

The integer representation of PauliString(nqubits: 6, 1.0 * IXIYIZ) is 3204
The 2d visualization of PauliString(nqubits: 6, 1.0 * IXIYIZ) is:
IX
IY
IZ

The integer representation of PauliString(nqubits: 6, 1.0 * IYIZIX) is 1224
The 2d visualization of PauliString(nqubits: 6, 1.0 * IYIZIX) is:
IY
IZ
IX

# Turns out for 2d finding the lowest integer representation can be speed up 
# by precomputing "masks" for translations in x and y directions. 
# See `_translatetolowestinteger` in the source code for details.
# We will just call the function `translationmerge` here.

# The representative Pauli string is given by
psum = PauliSum([pstr, pstr2, pstr3])
rep_int_2d = translationmerge(psum, nx, ny)

PauliSum(nqubits: 6, 1 Pauli term: 
 3.0 * YIZIXI
)

You can also define the merging function for your own symmetry. Let's see an example of this using permutation symmetry

We can define the representative of Paulis under permutation symmetry by counting the number of Pauli $X, Y, Z$

function _permutationcanonicalform(pstr::PauliStringType)
    if pstr == 0
        return pstr
    end

    num_x = countx(pstr)
    num_y = county(pstr)
    num_z = countz(pstr)

    # construct the canonical form PauliString
    pstr_canonical = identitylike(pstr)
    for ii in 1:num_x
        pstr_canonical = setpauli(pstr_canonical, :X, ii)
    end
    for ii in (num_x+1):(num_x+num_y)
        pstr_canonical = setpauli(pstr_canonical, :Y, ii)
    end
    for ii in (num_x+num_y+1):(num_x+num_y+num_z)
        pstr_canonical = setpauli(pstr_canonical, :Z, ii)
    end

    return pstr_canonical
end

_permutationcanonicalform (generic function with 1 method)

# Now define a shorthand for merging under permutation symmetry
permutationmerge(psum::PauliSum) = symmetrymerge(
    psum, (pstr -> _permutationcanonicalform(pstr))
)

permutationmerge (generic function with 1 method)

# Let's see some examples of merging!
pstr1 = PauliString(4, [:X, :Y, :Z], [1, 2, 3])
pstr2 = PauliString(4, [:Y, :Z, :X], [2, 3, 4])
pstr3 = PauliString(4, [:Z, :X, :Y], [1, 3, 4])
psum = PauliSum([pstr1, pstr2, pstr3])
rep_perm = permutationmerge(psum)
println("The representative under permutation symmetry is: ", rep_perm)

The representative under permutation symmetry is: PauliSum(nqubits: 4, 1 Pauli term: 
 3.0 * XYZI
)

Let's use translation symmetry in propagation of a Hamiltonian evolution

Consider the dynamics by the tilted-field Ising model implemented in tiltedtfitrottercircuit.

\[ U(\Delta t) = \prod_{a=1}^{l} \prod_{j=1}^n e^{-i \Delta t g_x X_j} \prod_{j=1}^{n-1} e^{-i \Delta t g_z Z_j} e^{-i \Delta t g_z Z_j Z_{j+1}}. \]

# Let's set up the circuit
nq = 11
periodic = true  # periodic boundary condition for translation symmetry
deltat = 0.1
layers = [3, 5]  # number of layers to look at
circ = tiltedtfitrottercircuit(nq, 1; topology=bricklayertopology(nq, periodic=periodic))

33-element Vector{Gate}:
 PauliRotation([:Z, :Z], [1, 2])
 PauliRotation([:Z, :Z], [3, 4])
 PauliRotation([:Z, :Z], [5, 6])
 PauliRotation([:Z, :Z], [7, 8])
 PauliRotation([:Z, :Z], [9, 10])
 PauliRotation([:Z, :Z], [11, 1])
 PauliRotation([:Z, :Z], [2, 3])
 PauliRotation([:Z, :Z], [4, 5])
 PauliRotation([:Z, :Z], [6, 7])
 PauliRotation([:Z, :Z], [8, 9])
 PauliRotation([:Z, :Z], [10, 11])
 PauliRotation([:Z], [1])
 PauliRotation([:Z], [2])
 ⋮
 PauliRotation([:Z], [11])
 PauliRotation([:X], [1])
 PauliRotation([:X], [2])
 PauliRotation([:X], [3])
 PauliRotation([:X], [4])
 PauliRotation([:X], [5])
 PauliRotation([:X], [6])
 PauliRotation([:X], [7])
 PauliRotation([:X], [8])
 PauliRotation([:X], [9])
 PauliRotation([:X], [10])
 PauliRotation([:X], [11])

# Helper functions for the evolution angles
function get_tfim_couplings(circuit, deltat::Float64; J=1, hz=0.9045, hx=1.4)
    thetas = []

    indices_gate_zz = getparameterindices(circuit, PauliRotation, [:Z, :Z])
    indices_gate_x = getparameterindices(circuit, PauliRotation, [:X])
    indices_gate_z = getparameterindices(circuit, PauliRotation, [:Z])

    for (idx, gate) in enumerate(circuit)

        if idx in indices_gate_zz
            coupling_strength = - deltat * J
        elseif idx in indices_gate_x
            coupling_strength = - deltat * hx
        elseif idx in indices_gate_z
            coupling_strength = - deltat * hz
        else
            throw(ArgumentError("Gate at index $idx is neither ZZ, X nor Z rotation"))
        end
        push!(thetas, coupling_strength * 2) # e^i t P / 2 * 2
    end
            
    return thetas
end

get_tfim_couplings (generic function with 1 method)

# Build evolution angles for TFIM model
nl = maximum(layers)
thetas = get_tfim_couplings(circ, deltat) # angle per layer

# initialie a list to store results
results = []

for merging in (:full, :none)  # :full means translation merging

    # Consider the inital pauli string
    p_evolved = PauliString(nq, :Z, div(nq, 2) + 1)
    println("Initial Pauli string: ", p_evolved)  

    for l in collect(1:maximum(layers))
      
        # Evolution 1 full Trotter layer
        p_evolved = propagate(circ, p_evolved, thetas; layer=1)

        if merging == :none
            # no merging
        elseif merging == :full
            p_evolved = translationmerge(p_evolved)
        end

        if l in layers
            println("Layer $l with merging type: $merging")
            overlap = overlapwithzero(p_evolved)
            push!(results, (
                merging=merging, layer=l, nterms=length(p_evolved),
                overlap=overlap
                )
            )
        end
    end
end    

Initial Pauli string: PauliString(nqubits: 11, 1.0 * IIIIIZIIIII)


Layer 3 with merging type: full


Layer 5 with merging type: full


Initial Pauli string: PauliString(nqubits: 11, 1.0 * IIIIIZIIIII)
Layer 3 with merging type: none
Layer 5 with merging type: none

# Check the results:
results

4-element Vector{Any}:
 (merging = :full, layer = 3, nterms = 1422, overlap = 0.7309034078384146)
 (merging = :full, layer = 5, nterms = 12172, overlap = 0.5258230815192094)
 (merging = :none, layer = 3, nterms = 1972, overlap = 0.7309034078384142)
 (merging = :none, layer = 5, nterms = 34617, overlap = 0.5258230802947047)

# verify that the expectatin values of symmetry merged Pauli strings are the same as those without merging
for l in layers
    overlap_none = filter(r -> r.merging == :none && r.layer == l, results)[1].overlap
    overlap_full = filter(r -> r.merging == :full && r.layer == l, results)[1].overlap
    @assert isapprox(overlap_none, overlap_full)
end

Notice that the minimal symmetric circuit layer is actually not necessary the full layer! For instance, instead of working with the trotter layer implemented in tiltedtfitrottercircuit. let's merge after each sub-layers as the follows: apply

\[ U(\Delta t)_Z = \prod_{j=1}^{n-1} e^{-i \Delta t g_z Z_j} e^{-i \Delta t g_z Z_j Z_{j+1}}, \]

merge and apply

\[ U(\Delta t)_X = \prod_{j=1}^n e^{-i \Delta t g_x X_j}, \]

merge. The two sublayers are each translation invariant but form the total Trotter layer.

# let's do the above merging for the sublayers
circ_z = circ[1:2*nq]  # get the ZZ and Z layers
circ_x = circ[2*nq+1:end]  # get the X layers
thetas_z = thetas[1:2*nq]
thetas_x = thetas[2*nq+1:end]

p_evolved = PauliString(nq, :Z, div(nq, 2) + 1)
println("Initial Pauli string: ", p_evolved)  

for l in collect(1:maximum(layers))

    # Evolution under Z layer
    p_evolved = propagate(circ_z, p_evolved, thetas_z; layer=1)
    p_evolved = translationmerge(p_evolved)

    # Evolution under X layer
    p_evolved = propagate(circ_x, p_evolved, thetas_x; layer=1)
    p_evolved = translationmerge(p_evolved)

    if l in layers

        overlap = overlapwithzero(p_evolved)
        push!(results, (
            merging=:sublayer, layer=l, nterms=length(p_evolved),
            overlap=overlap
            )
        )
    
    end
end            

Initial Pauli string: PauliString(nqubits: 11, 1.0 * IIIIIZIIIII)

results

6-element Vector{Any}:
 (merging = :full, layer = 3, nterms = 1422, overlap = 0.7309034078384146)
 (merging = :full, layer = 5, nterms = 12172, overlap = 0.5258230815192094)
 (merging = :none, layer = 3, nterms = 1972, overlap = 0.7309034078384142)
 (merging = :none, layer = 5, nterms = 34617, overlap = 0.5258230802947047)
 (merging = :sublayer, layer = 3, nterms = 300, overlap = 0.7309034078384145)
 (merging = :sublayer, layer = 5, nterms = 6753, overlap = 0.5258230815192095)

! Interestingly merging every sublayers reduce the number of Pauli terms even further!

# verify that the expectatin values of symmetry merged Pauli strings are the same as those without merging
for l in layers
    overlap_none = filter(r -> r.merging == :none && r.layer == l, results)[1].overlap
    overlap_full = filter(r -> r.merging == :full && r.layer == l, results)[1].overlap
    overlap_sublayer = filter(r -> r.merging == :sublayer && r.layer == l, results)[1].overlap
    @assert isapprox(overlap_none, overlap_full)
    @assert isapprox(overlap_none, overlap_sublayer)
end