Converting circuit gates to Clifford + T gates

PennyLane Help
1

Dear pennylane team,
I have created a circuit by performing the following operations:

First, I created the Hamiltonian.

symbols = ["H", "O", "H"]
coordinates = np.array([-0.0399, -0.0038, 0.0, 1.5780, 0.8540, 0.0, 2.7909, -0.5159, 0.0])
hamiltonian, num_qubits = qchem.molecular_hamiltonian(symbols, coordinates)

Second, I passed the Hamiltonian down to the ApproxTimeEvolution function.

wires = range(num_qubits)
num_wires = num_qubits

dev = qml.device('default.qubit', wires=wires)


@qml.qnode(dev)
def circuit(time):
    ApproxTimeEvolution(hamiltonian, time, 1)
    return [qml.expval(qml.PauliZ(wires=i)) for i in wires]

Next, I would like to see what the circuit is made of. So I run the spec_func and get the following not-informative output:

{'gate_sizes': defaultdict(int, {14: 1}),
 'gate_types': defaultdict(int, {'ApproxTimeEvolution': 1}),
 'num_operations': 1,
 'num_observables': 14,
 'num_diagonalizing_gates': 0,
 'num_used_wires': 14,
 'depth': 1,
 'num_trainable_params': 1,
 'num_device_wires': 14,
 'device_name': 'default.qubit.autograd',
 'expansion_strategy': 'gradient',
 'gradient_options': {},
 'interface': 'autograd',
 'diff_method': 'best',
 'gradient_fn': 'backprop'}`

The problem with this is that it doesn’t show what gate types are used but just says ApproxTimeEvolution.

So, I have two questions:
1. How do I actually see what kind of gates are used in the circuit?
2. I would like the circuit to use exclusively the Clifford + T gate set. How do I do this?

Thanks!

2

Hi @Einar_Gabbassov, it’s great to see you in the Forum!

  1. You can add the keyword argument expansion_strategy='device'.
    You can also use print(qml.draw(circuit, expansion_strategy='device')(time)) and replace time with the value you want.

  2. If you want your circuit to use exclusively a set of gates then you can find a device which uses such gates or create a custom device with the gates that you want. You can learn how to create your own device in this blog post. Then expansion_strategy='device' should decompose the circuit into the gates in that device.

Please let me know if this works for you!

3

Hi @CatalinaAlbornoz, thanks for tips.

Unfortunattely, the suggested workflow doesn’t yield the desired gates.
It actually produces something else. Here is the code:

import numpy as np
import pennylane as qml
from pennylane import qchem
from pennylane.templates import ApproxTimeEvolution

class MyDevice(qml.QubitDevice):
    name = 'Clifford_T_device'
    short_name = "custom.qubit"
    pennylane_requires = ">=0.23"
    version = "0.0.1"
    author = "Maurice"
    # Here, I define the desired gate set to be used
    operations = {"H", "S", "CNOT", "T"}
    observables = {"PauliX", "PauliZ"}

    def __init__(self, num_wires):
        super().__init__(wires=num_wires, shots=None)

        # create the initial state
        self._state = np.array([1, 0])

        # create a variable for future copies of the state
        self._pre_rotated_state = None

    @property
    def state(self):
        return self._pre_rotated_state

    @classmethod
    def capabilities(cls):
        capabilities = super().capabilities().copy()
        capabilities.update(
            returns_state=True,
            supports_finite_shots=False,
            supports_tensor_observables=False
        )
        return capabilities

    def apply(self, operations, rotations=None, **kwargs):
        for op in operations:
            # We update the state by applying the matrix representation of the gate
            self._state = qml.matrix(op) @ self._state

        # store the pre-rotated state
        self._pre_rotated_state = self._state.copy()

        # apply the circuit rotations
        for rot in rotations or []:
            self._state = qml.matrix(rot) @ self._state

    def analytic_probability(self, wires=None):
        if self._state is None:
            return None

        real = self._real(self._state)
        imag = self._imag(self._state)
        prob = self.marginal_prob(real ** 2 + imag ** 2, wires)
        return prob

    def reset(self):
        """Reset the device"""
        self._state = np.array([1, 0])

hamiltonian, num_qubits = qchem.molecular_hamiltonian(symbols, coordinates)
wires = range(num_qubits)
num_wires = num_qubits

dev = MyDevice(num_wires)

time=1
trotter_steps=1
@qml.qnode(dev)
def circuit(time):
    ApproxTimeEvolution(hamiltonian, time, trotter_steps)
    return [qml.expval(qml.PauliZ(wires=i)) for i in wires]

specs_func = qml.specs(circuit, expansion_strategy='device')
print(specs_func(1))
t = circuit.qtape.to_openqasm()
print(t)

Even though the indicated gate set is universal, I feel like for the suggested approach to work we need some sort of Solovay-Kitaev algorithm, which will assist the conversion process. I’m not sure if the suggested approach anyhow takes this into account. Please correct me if I’m wrong.

For you convenience, I also attach the output

OPENQASM 2.0;
include "qelib1.inc";
qreg q[4];
creg c[4];
rz(0.4740443816356016) q[0];
rz(0.4740443816356016) q[1];
cx q[1],q[0];
rz(0.36906324698635556) q[0];
cx q[1],q[0];
rx(1.5707963267948966) q[0];
rz(1.5707963267948966) q[1];
rx(1.5707963267948966) q[1];
rz(1.5707963267948966) q[1];
rz(1.5707963267948966) q[2];
rx(1.5707963267948966) q[2];
rz(1.5707963267948966) q[2];
rx(1.5707963267948966) q[3];
cx q[3],q[2];
cx q[2],q[1];
cx q[1],q[0];
rz(0.08210294405757057) q[0];
cx q[1],q[0];
cx q[2],q[1];
cx q[3],q[2];
rx(-1.5707963267948966) q[0];
rz(1.5707963267948966) q[1];
rx(1.5707963267948966) q[1];
rz(1.5707963267948966) q[1];
rz(1.5707963267948966) q[2];
rx(1.5707963267948966) q[2];
rz(1.5707963267948966) q[2];
rx(-1.5707963267948966) q[3];
rx(1.5707963267948966) q[0];
rx(1.5707963267948966) q[1];
rz(1.5707963267948966) q[2];
rx(1.5707963267948966) q[2];
rz(1.5707963267948966) q[2];
rz(1.5707963267948966) q[3];
rx(1.5707963267948966) q[3];
rz(1.5707963267948966) q[3];
cx q[3],q[2];
cx q[2],q[1];
cx q[1],q[0];
rz(-0.08210294405757057) q[0];
cx q[1],q[0];
cx q[2],q[1];
cx q[3],q[2];
rx(-1.5707963267948966) q[0];
rx(-1.5707963267948966) q[1];
rz(1.5707963267948966) q[2];
rx(1.5707963267948966) q[2];
rz(1.5707963267948966) q[2];
rz(1.5707963267948966) q[3];
rx(1.5707963267948966) q[3];
rz(1.5707963267948966) q[3];
rz(1.5707963267948966) q[0];
rx(1.5707963267948966) q[0];
rz(1.5707963267948966) q[0];
rz(1.5707963267948966) q[1];
rx(1.5707963267948966) q[1];
rz(1.5707963267948966) q[1];
rx(1.5707963267948966) q[2];
rx(1.5707963267948966) q[3];
cx q[3],q[2];
cx q[2],q[1];
cx q[1],q[0];
rz(-0.08210294405757057) q[0];
cx q[1],q[0];
cx q[2],q[1];
cx q[3],q[2];
rz(1.5707963267948966) q[0];
rx(1.5707963267948966) q[0];
rz(1.5707963267948966) q[0];
rz(1.5707963267948966) q[1];
rx(1.5707963267948966) q[1];
rz(1.5707963267948966) q[1];
rx(-1.5707963267948966) q[2];
rx(-1.5707963267948966) q[3];
rz(1.5707963267948966) q[0];
rx(1.5707963267948966) q[0];
rz(1.5707963267948966) q[0];
rx(1.5707963267948966) q[1];
rx(1.5707963267948966) q[2];
rz(1.5707963267948966) q[3];
rx(1.5707963267948966) q[3];
rz(1.5707963267948966) q[3];
cx q[3],q[2];
cx q[2],q[1];
cx q[1],q[0];
rz(0.08210294405757057) q[0];
cx q[1],q[0];
cx q[2],q[1];
cx q[3],q[2];
rz(1.5707963267948966) q[0];
rx(1.5707963267948966) q[0];
rz(1.5707963267948966) q[0];
rx(-1.5707963267948966) q[1];
rx(-1.5707963267948966) q[2];
rz(1.5707963267948966) q[3];
rx(1.5707963267948966) q[3];
rz(1.5707963267948966) q[3];
rz(-0.9211542548980106) q[2];
cx q[2],q[0];
rz(0.2812358847255707) q[0];
cx q[2],q[0];
rz(-0.9211542548980107) q[3];
cx q[3],q[0];
rz(0.3633388287831413) q[0];
cx q[3],q[0];
cx q[2],q[1];
rz(0.3633388287831413) q[1];
cx q[2],q[1];
cx q[3],q[1];
rz(0.2812358847255707) q[1];
cx q[3],q[1];
cx q[3],q[2];
rz(0.38350121257667436) q[2];
cx q[3],q[2];
measure q[0] -> c[0];
measure q[1] -> c[1];
measure q[2] -> c[2];
measure q[3] -> c[3];
4

Hey @Einar_Gabbassov,

Does the decomposition() method work for you? E.g.,

>>> hamiltonian, num_qubits = qchem.molecular_hamiltonian(["H", "H"], np.array([[0, 0, 0.66], [0, 0, -0.66]]))
>>> qml.ApproxTimeEvolution(hamiltonian, 1, 1).decomposition()
[PauliRot(tensor(0.35590551, requires_grad=False), wires=[0]), PauliRot(tensor(0.35590551, requires_grad=False), wires=[1]), PauliRot(tensor(0.34133649, requires_grad=False), wires=[0, 1]), PauliRot(tensor(0.08945978, requires_grad=False), wires=[0, 1, 2, 3]), PauliRot(tensor(-0.08945978, requires_grad=False), wires=[0, 1, 2, 3]), PauliRot(tensor(-0.08945978, requires_grad=False), wires=[0, 1, 2, 3]), PauliRot(tensor(0.08945978, requires_grad=False), wires=[0, 1, 2, 3]), PauliRot(tensor(-0.48698732, requires_grad=False), wires=[2]), PauliRot(tensor(0.2460383, requires_grad=False), wires=[0, 2]), PauliRot(tensor(-0.48698732, requires_grad=False), wires=[3]), PauliRot(tensor(0.33549807, requires_grad=False), wires=[0, 3]), PauliRot(tensor(0.33549807, requires_grad=False), wires=[1, 2]), PauliRot(tensor(0.2460383, requires_grad=False), wires=[1, 3]), PauliRot(tensor(0.35269298, requires_grad=False), wires=[2, 3])]

If those operations aren’t decomposed enough, you can decompose them further:

>>> for op in qml.ApproxTimeEvolution(hamiltonian, 1, 1).decomposition():
...     print(op.decomposition())
... 
[MultiRZ(tensor(0.35590551, requires_grad=False), wires=[0])]
[MultiRZ(tensor(0.35590551, requires_grad=False), wires=[1])]
[MultiRZ(tensor(0.34133649, requires_grad=False), wires=[0, 1])]
[RX(1.5707963267948966, wires=[0]), Hadamard(wires=[1]), Hadamard(wires=[2]), RX(1.5707963267948966, wires=[3]), MultiRZ(tensor(0.08945978, requires_grad=False), wires=[0, 1, 2, 3]), RX(-1.5707963267948966, wires=[0]), Hadamard(wires=[1]), Hadamard(wires=[2]), RX(-1.5707963267948966, wires=[3])]
[RX(1.5707963267948966, wires=[0]), RX(1.5707963267948966, wires=[1]), Hadamard(wires=[2]), Hadamard(wires=[3]), MultiRZ(tensor(-0.08945978, requires_grad=False), wires=[0, 1, 2, 3]), RX(-1.5707963267948966, wires=[0]), RX(-1.5707963267948966, wires=[1]), Hadamard(wires=[2]), Hadamard(wires=[3])]
[Hadamard(wires=[0]), Hadamard(wires=[1]), RX(1.5707963267948966, wires=[2]), RX(1.5707963267948966, wires=[3]), MultiRZ(tensor(-0.08945978, requires_grad=False), wires=[0, 1, 2, 3]), Hadamard(wires=[0]), Hadamard(wires=[1]), RX(-1.5707963267948966, wires=[2]), RX(-1.5707963267948966, wires=[3])]
[Hadamard(wires=[0]), RX(1.5707963267948966, wires=[1]), RX(1.5707963267948966, wires=[2]), Hadamard(wires=[3]), MultiRZ(tensor(0.08945978, requires_grad=False), wires=[0, 1, 2, 3]), Hadamard(wires=[0]), RX(-1.5707963267948966, wires=[1]), RX(-1.5707963267948966, wires=[2]), Hadamard(wires=[3])]
[MultiRZ(tensor(-0.48698732, requires_grad=False), wires=[2])]
[MultiRZ(tensor(0.2460383, requires_grad=False), wires=[0, 2])]
[MultiRZ(tensor(-0.48698732, requires_grad=False), wires=[3])]
[MultiRZ(tensor(0.33549807, requires_grad=False), wires=[0, 3])]
[MultiRZ(tensor(0.33549807, requires_grad=False), wires=[1, 2])]
[MultiRZ(tensor(0.2460383, requires_grad=False), wires=[1, 3])]
[MultiRZ(tensor(0.35269298, requires_grad=False), wires=[2, 3])]
5

hi @isaacdevlugt and @CatalinaAlbornoz
The call qml.ApproxTimeEvolution(hamiltonian, 1, 1).decomposition() outputs the list of gates. This is great! Thanks. Also, thanks a lot for the tip on how to decompoe even further. This was extremly useful.

However, the second question is still pending. My final goal is not to preview the abstract gates produced by pennylane but actually take the resulting circuit and approximate it in terms of “usefull” gates which are Cliffotd + T gates. Is it possible to do that?

6

Right! Apologies about that. I joined this thread a little late in the game :sweat_smile:.

Currently, PennyLane doesn’t support what you’re after specifically. I think this would be a very nice feature to have, though. As such, I created a Feature Request issue! Feel free to contribute to it if you’d like and track its progress.