Gate synthesis/decomposition in Pennylane

Hi all, I am trying to perform Hamiltonian simulation in pennylane and I am facing some issues although the concept is very simple. The actual code is quite large so here is a part of what I am trying to do: during gate synthesis when I am reducing the L_{2} norm between the actual gate output(e.g. Hadamard) and parameterized gate (here it would be SU(2)(\alpha, \beta, \gamma) for single wire==qubit) via GradientDescentOptimizer the final output i.e. the SU(2) with optimized \alpha_{optimzed}, \beta_{optimzed}, \gamma_{optimzed} is not providing the exact Hadamard. I know that if we break the Hadamard it will simply give R_y(-\pi/2).\sigma_x but still, I am getting

[[ 0.68851887-0.22800782j, -0.68840703+0.00706844j],
[ 0.68840703+0.00706844j, 0.68851887+0.22800782j]]

instead of

[[ 0.70710678, 0.70710678],
[ 0.70710678, -0.70710678]]

so there is always an issue with the -ve sign.

The code :

dev=qml.device("default.qubit",wires=1)
@qml.qnode(dev)
def gate(l):
    #qml.BasisState(np.array([0], requires_grad=False), wires=0)
    qml.Rot(l[0],l[1],l[2],wires=0)
    #qml.Hadamard(wires=0)

    return qml.state()

@qml.qnode(dev)
def had_op():
    qml.Hadamard(wires=0)

    return qml.state()

def f(l):  #f=cost
    
    return np.linalg.norm(np.real(gate(l))-had_op())**2

np.random.seed(0)
init_param=np.random.rand(3)
params=init_param
print(params)

cost_history = [  ]
steps=2000
for it in range(steps):
    params, cost = qml.GradientDescentOptimizer(stepsize=0.01).step_and_cost(f, params)
    if (it+1)%200==0:
        
        print("the value of cost at step {:} is {:15f}".format(it+1,f(params)))

        cost_history.append(cost)
print(cost,params)


qml.Rot(params[0],params[1],params[2],wires=0).matrix >>

[[ 0.68851887-0.22800782j, -0.68840703+0.00706844j],
       [ 0.68840703+0.00706844j,  0.68851887+0.22800782j]]

Kindly have a look if you can help me in this regard! Thanks in advance

Hey @Arun!

I had a go at running your code with stepsize=0.1 and 4000 steps, and the resulting matrix looked much closer:

array([[ 0.70538886-0.06962842j, -0.70538878+0.0022627j ],
       [ 0.70538878+0.0022627j ,  0.70538886+0.06962842j]])

Note that this is not quite the Hadamard matrix, and its decomposition can be achieved using

qml.PhaseShift(np.pi / 2, wires=0)
qml.RX(np.pi / 2, wires=0)
qml.PhaseShift(np.pi / 2, wires=0)

Perhaps introducing trainable parameters into those gates could help you optimize.

Hey @Tom_Bromley,

Thanks for your response. Can you just clarify that whether it is possible to decompose Hadamard in terms of only rotation gates? As we know H=R_y(-\pi/2).(e^{i\pi/2} R_x(\pi)) but I wish to obtain this decomposition via stochastic optimization!

Hey @Arun!

Can you just clarify that whether it is possible to decompose Hadamard in terms of only rotation gates?

Do you mean specifically in terms of qml.Rot? The decomposition above into phase shifts and X rotations is a standard decomposition into single-qubit rotations. The code below shows how you can optimize the variational parameters for this decomposition:

import pennylane as qml
from pennylane import numpy as np

dev = qml.device("default.qubit.autograd", wires=1)

def decomposition(params):
    qml.PhaseShift(params[0], wires=0)
    qml.RX(params[1], wires=0)
    qml.PhaseShift(params[2], wires=0)

@qml.qnode(dev, diff_method="backprop")
def f(params, state, apply_inv=True):
    qml.BasisState(np.array([state]), wires=0)
    decomposition(params)
    if apply_inv:
        qml.Hadamard(wires=0).inv()
    return qml.state()

params = np.random.random(3)
Id = np.eye(2, requires_grad=False)

def cost_state(params, state):
    return np.sum(np.abs(f(params, state) - Id[state]))

def cost(params):
    return sum(cost_state(params, i) for i in range(2))

opt = qml.GradientDescentOptimizer(stepsize=0.01)

for i in range(300):
    params = opt.step(cost, params)

    if i % 20 == 0:
        print(f"Cost at step {i}:", cost(params))
        
def get_unitary(params):
    return np.array([f(params, i, apply_inv=False) for i in range(2)]).T

print("\nOptimized unitary:\n", get_unitary(params))
print("Ideal unitary:\n", qml.Hadamard.matrix)

print("\nOptimized params:", params.numpy())
print("Expected params: ", np.array([np.pi / 2] * 3))

Note that this optimization ensures that both input states |0> and |1> are mapped to their expected outputs.

I believe using Rot alone can only approximate the Hadamard gate up to a global phase, so won’t work in the above example. You could consider changing the cost function from the norm to the fidelity (which is insensitive to global phase) and using just Rot will then probably work.

Thanks!

1 Like

Hi @Tom_Bromley,
Thank you very much for your beautiful response and the above solution is perfect. And as a response to your query, I thought we can decompose any unitary in terms of rotations and a global phase because that’s what we need to produce any general unitary \in SU(2) i.e. U=e^{i\phi} R_{n1}(\theta_1)R_{n2}(\theta_2)R_{n1}(\theta_3). But still the above solution worth it.

Thanks again :slight_smile: ,

No problem!

U=e^{i\phi} R_{n1}(\theta_1)R_{n2}(\theta_2)R_{n1}(\theta_3)

Right! You will probably have some luck with combining a Rot gate with a PhaseShift gate, using 4 parameters overall. Good luck!

1 Like

Thanks, @Tom_Bromley again :slight_smile:,
your suggestion(Phase + Rotation) works properly with a total of 1000 steps.

1 Like