Hilbert Schmidt loss does not reproduce the desired unitary

Hello all,

I am trying to compile an arbitrary unitary to a quantum circuit using qml.HilbertSchmidt. The training is successful (converges to 10^(-5)), but when i transform the circuit to a matrix, the trace distance between the target unitary is large (2). Does anyone know what is going on?
And i know it is possible since they do it in a paper https://journals.aps.org/prx/abstract/10.1103/PhysRevX.11.041036

For instance

U = [[1,0,0,0],[0,sqrt(2)/2,sqrt(2)/2,0],[0,-sqrt(2)/2,sqrt(2)/2,0],[0,0,0,1]]
N= 2
depth = 5
u_wires = np.arange(N)
v_wires = np.arange(N,2*N)

with qml.tape.QuantumTape(do_queue=False) as u_tape:
    qml.QubitUnitary(U, wires=u_wires)

def v_function(weights,wires):
    D = np.shape(weights)[-1]
    for d in range(D):
        #rotation
        for i in range(N):
            qml.RY(weights[i,it,d],wires=wires[i])
               
        for i in range(N-1):
                qml.CNOT(wires=[wires[i],wires[(i+1)%N]])


dev = qml.device("lightning.qubit", wires=2*N)
@qml.qnode(dev)
def hilbert_test(v_params, v_function, v_wires, u_tape):
    qml.HilbertSchmidt(v_params, v_function=v_function, v_wires=v_wires, u_tape=u_tape)
    return qml.probs(u_tape.wires + v_wires)

def cost_hst(parameters, v_function, v_wires, u_tape):
    return (1 - hilbert_test(v_params=parameters, v_function=lambda p:v_function(p,v_wires), v_wires=v_wires, u_tape=u_tape)[0])


loss_function = lambda v_params: cost_hst(v_params, v_function, v_wires, u_tape)
opt = qml.GradientDescentOptimizer(stepsize=0.04)
parameters = qml.numpy.array(np.random.normal(0,1,size=(N,1,depth)), requires_grad=True)

cost = []
for _ in range(200):

    parameters = opt.step(loss_function, parameters)
    cost.append(loss_function(parameters))

def f(params,wires):
    v_function(params,wires)

but if i then print the corresponding matrix with the trained parameters, using the technique explained here Printing the matrix of a circuit during optimization, i obtain

 [[-0.38304712+0.28415248j  0.44274339+0.20951309j     0.02119074+0.46715616j
-0.51994317+0.20878219j]
 [-0.20400715-0.52963304j  0.54954541-0.55876896j -0.0565359 +0.10391468j
 0.1190677 -0.18834896j]
[-0.18387383+0.2727543j  -0.03200142+0.07518786j  -0.27096935+0.5267104j
0.72912911-0.0513877j ]
 [-0.42870162+0.39801083j  0.35100433+0.12640534j  0.0427307 -0.64393117j
 0.20939861-0.24145166j]]

which is far away from what i am looking for (U).

Thanks
Cheers
Oriel

Hi @Oriel!
I tried to run your code and I don’t get small losses.
This is the script I used (with some modifications to make it work for me):

import pennylane as qml
from pennylane import numpy as np

U = [[1,0,0,0],[0,np.sqrt(2)/2,np.sqrt(2)/2,0],[0,-np.sqrt(2)/2,np.sqrt(2)/2,0],[0,0,0,1]]
N= 2
depth = 5
u_wires = range(N)
v_wires = range(N,2*N)

with qml.tape.QuantumTape(do_queue=False) as u_tape:
    qml.QubitUnitary(U, wires=u_wires)

def v_function(weights,wires):
    D = np.shape(weights)[-1]
    for d in range(D):
        #rotation
        for i in range(N):
            qml.RY(weights[i,0,d],wires=wires[i])
               
        for i in range(N-1):
                qml.CNOT(wires=[wires[i],wires[(i+1)%N]])


dev = qml.device("default.qubit", wires=2*N)
@qml.qnode(dev)
def hilbert_test(v_params, v_function, v_wires, u_tape):
    qml.HilbertSchmidt(v_params, v_function=v_function, v_wires=v_wires, u_tape=u_tape)
    return qml.probs(u_tape.wires + v_wires)

def cost_hst(parameters, v_function, v_wires, u_tape):
    return (1 - hilbert_test(v_params=parameters, v_function=lambda p:v_function(p,v_wires), v_wires=v_wires, u_tape=u_tape)[0])


loss_function = lambda v_params: cost_hst(v_params, v_function, v_wires, u_tape)
opt = qml.GradientDescentOptimizer(stepsize=0.04)
parameters = qml.numpy.array(np.random.normal(0,1,size=(N,1,depth)), requires_grad=True)
cost = []
for _ in range(200):

    parameters = opt.step(loss_function, parameters)
    cost.append(loss_function(parameters))

def f(params,wires):
    v_function(params,wires)
    
import matplotlib.pyplot as plt

plt.plot(cost)

The error is 0.75.
Make sure again that the loss is the same and let me know :blush:

1 Like

Hi, actually i wanted to make my code easier to understand and so it does not work anymore. If you add a RX and RY rotation with different parameters after the first RY rotation (si you have a general rotation) it should get to arbitrary low values.
Thanks for checking this up

Hi @Oriel, I’m glad Guillermo’s answer solved your problem. Please let us know if you find any other issues going forward!

Hi, i am afraid not. what i meant is that if i use a RY, RX, RY rotation, i get 10^-5 loss, but if i print the matrix of the circuit, this is not the same as the target matrix. Maybe i am doing something wrong by retrieving the matrix but for me it is simply v_function. I can add some details on how i extract the matrix tomorrow. Thanks

You are right, I have been checking the code and there is a bug when doing the inverse of V. I will try to fix it as soon as possible. Thank you very much for the warning! :rocket:

1 Like

Hey again! I have been talking to the team and the code is fine as it is following the format of the original paper. That is, the solution you are going to get is V* instead of V, so you will have to do the conjugate of the circuit if you want to see that the matrix is the same. Here is a small script that I have used to help me.

dev = qml.device("default.qubit", wires = v_wires)
@qml.qnode(dev)
def circuit():
    v_function(parameters)
    return qml.probs(wires = v_wires)

tape = circuit.qtape

with qml.tape.QuantumTape() as tape2:
    for op in tape.operations:
        qml.adjoint(op)
        
tape2.operations
for i in qml.matrix(tape2):
    print([round(float(abs(j)),2) for j in i])

[1.0, 0.04, 0.05, 0.02]
[0.02, 0.7, 0.72, 0.02]
[0.07, 0.71, 0.69, 0.05]
[0.01, 0.04, 0.04, 1.0]

Right now we have to use tapes to solve this but I have been told that they are going to work on a qml.conj() which will make the work easier in the future.