Encountering difficulties when trying to use the fidelity as cost function

PennyLane Help
1

Hello! I’m trying to train a PQC using initial states and target states as training statistics, and the fidelity between the output states of the PQC and the target states as my cost function. However, I find it difficult to implement a cost function returning the value 1-fidelity. My original codes of my circuit is long, so let’s take a simple circuit with a 3-bit special unitary as example:

def SU8(params, wires):
# circuit implementing a 3-bit arbitrary operation
    if len(params)<63:
        raise ValueError('The number of parameters should be 63 to construct a SU(8).')
    qml.RZ(params[0], wires[0])
    uniform_CZ(params[1:3], wires[0], wires[1]) # The definition of the uniform controlled rotations will be provided at the end of the topic so that we can concentrate on the problem first
    uniform_CZ(params[3:7], [wires[0], wires[1]], wires[2])
    uniform_CY(params[7:11], [wires[0], wires[1]], wires[2])
    uniform_CZ(params[11:15], [wires[0], wires[1]], wires[2])
    uniform_CY(params[15:19], [wires[0], wires[2]], wires[1])
    uniform_CZ(params[19:23], [wires[0], wires[2]], wires[1])
    uniform_CY(params[23:27], [wires[0], wires[1]], wires[2])
    uniform_CZ(params[27:31], [wires[0], wires[1]], wires[2])
    uniform_CY(params[31:35], [wires[1], wires[2]], wires[0])
    uniform_CZ(params[35:39], [wires[1], wires[2]], wires[0])
    uniform_CY(params[39:43], [wires[0], wires[1]], wires[2])
    uniform_CZ(params[43:47], [wires[0], wires[1]], wires[2])
    uniform_CY(params[47:51], [wires[0], wires[2]], wires[1])
    uniform_CZ(params[51:55], [wires[0], wires[2]], wires[1])
    uniform_CY(params[55:59], [wires[0], wires[1]], wires[2])
    uniform_CZ(params[59:63], [wires[0], wires[1]], wires[2])

N=3 # Here I use an SU(8) circuit to show the problem, so I set up a 3-bit device. The problem will appear no matter what the bit number is.
dev = qml.device('lightning.gpu', wires=N)

I have tried multiple ways to implement the cost function, but I always got errors. For example, the following cost function:

@qml.qnode(dev)
def cost(params, ini_state, tar_state):
    qml.StatePrep(ini_state, range(N))
    SU8(params, range(N))
    return 1-qml.math.fidelity(qml.density_matrix(range(N)), qml.math.dm_from_state_vector(tar_state))

will throw the error TypeError: Input should be a NumPy array or scalar. It is a <class 'pennylane.measurements.state.DensityMatrixMP'> instead. at the return statement. Another cost function below:

@qml.qnode(dev)
def cost(params, ini_state, tar_state):
    qml.StatePrep(ini_state, range(N))
    SU8(params, range(N))
    return 1-np.abs(np.vdot(tar_state, qml.state()))

It will throw the error ValueError: cannot reshape array of size 1 into shape (8,) at the same line. The ini_state is a 8-dimention numpy array generated by numpy.random.rand() and the target_state is generated by target_matrix @ ini_state. I tried to print(type(tar_state)) and get <class 'pennylane.numpy.tensor.tensor'>, which is expected. It seems that the problem is not in my tar_state.

The only method I succeeded until now is using the Swap test:

N=3
dev2 = qml.device('lightning.gpu', wires=2*N+1)

@qml.qnode(dev2)
def cost(params, ini_state, tar_state):
    qml.StatePrep(ini_state, range(N))
    qml.StatePrep(tar_state, range(N+1, 2*N+1))
    SU8(params, range(N))
    return 1-swaptest(range(N), range(N+1, 2*N+1), N) # test the fidelity between range(N) and range(N+1, 2*N+1) with bit N as the ancilla bit. The code will also be provided at the end of the topic.

However, this will expand the bit number from N to 2*N+1, significantly increasing the training time and memory usage. Is there any simple method to implement a fidelity cost function without expanding the bit number?

Here is the code of uniform controlled rotations and Swap Test if needed:

def uniform_CR(params, c_wires, target, op):
    if type(c_wires) == int:
        c_wires=[c_wires]
    n = 2 ** len(c_wires)
    if len(params) < n:
        raise ValueError(f'The number of parameters is not enough. Need {n} parameters, given {len(params)}.')
    counter = 0
    for i in range(n):
        table = bin(i)[2:].zfill(len(c_wires))
        match op:
            case 'X':
                qml.ctrl(qml.RX(params[counter], target), c_wires, [bool(int(x)) for x in table])
            case 'Y':
                qml.ctrl(qml.RY(params[counter], target), c_wires, [bool(int(x)) for x in table])
            case 'Z':
                qml.ctrl(qml.RZ(params[counter], target), c_wires, [bool(int(x)) for x in table])
            case _:
                raise ValueError('Operations should be X, Y, or Z.')
        counter += 1

def uniform_CX(params, c_wires, target):
    uniform_CR(params, c_wires, target, 'X')

def uniform_CY(params, c_wires, target):
    uniform_CR(params, c_wires, target, 'Y')

def uniform_CZ(params, c_wires, target):
    uniform_CR(params, c_wires, target, 'Z')

def swaptest(wiresA, wiresB, ancillary:int):
    if len(wiresA)!=len(wiresB):
        raise ArithmeticError('Cannot evaluate the fidelity of 2 groups of qubits with different numbers because the fidelity of them is undefined.')
    for i in range(len(wiresA)):
        if not i:
            a=qml.SWAP(wires=[wiresA[0], wiresB[0]])
        else:
            a=a @ qml.SWAP(wires=[wiresA[i], wiresB[i]])
    qml.Hadamard(ancillary)
    qml.ctrl(a, ancillary)
    qml.Hadamard(ancillary)
    return qml.expval(qml.Z(ancillary))

Finally, the qml.about():

Name: PennyLane
Version: 0.40.0
Summary: PennyLane is a cross-platform Python library for quantum computing, quantum machine learning, and quantum chemistry. Train a quantum computer the same way as a neural network.
Home-page: https://github.com/PennyLaneAI/pennylane
Author: 
Author-email: 
License: Apache License 2.0
Location: /home/timmy/.local/lib/python3.12/site-packages
Requires: appdirs, autograd, autoray, cachetools, diastatic-malt, networkx, numpy, packaging, pennylane-lightning, requests, rustworkx, scipy, tomlkit, typing-extensions
Required-by: PennyLane_Lightning, PennyLane_Lightning_GPU

Platform info:           Linux-5.15.167.4-microsoft-standard-WSL2-x86_64-with-glibc2.39
Python version:          3.12.3
Numpy version:           2.2.2
Scipy version:           1.14.1
Installed devices:
- lightning.gpu (PennyLane_Lightning_GPU-0.40.0)
- lightning.qubit (PennyLane_Lightning-0.40.0)
- default.clifford (PennyLane-0.40.0)
- default.gaussian (PennyLane-0.40.0)
- default.mixed (PennyLane-0.40.0)
- default.qubit (PennyLane-0.40.0)
- default.qutrit (PennyLane-0.40.0)
- default.qutrit.mixed (PennyLane-0.40.0)
- default.tensor (PennyLane-0.40.0)
- null.qubit (PennyLane-0.40.0)
- reference.qubit (PennyLane-0.40.0)
None
2

Hi @Tz_19,

It looks like your optimizer was having some issues with the output of your cost function. It looks like it expected something of shape (8,), but in the first case your output was a matrix and in the second case it was a scalar.

Something you can do is obtain the fidelity by calculating the expectation value of a Hermitian operator as shown in our Data-reuploading classifier demo . In the demo you can see the following:

We construct an observable corresponding to the output label using the Hermitian operator. The expectation value of the observable gives the overlap or fidelity.

In the section on fidelity loss you can see how they construct the circuit and the cost function.

I don’t know how optimized this will be but you can try it out and see if this works for your specific case.

If you’re still having issues please make sure to share a minimal reproducible example of your code, so that we can copy-paste it and try to replicate your issue.

I hope this helps!

3

Hello @CatalinaAlbornoz , thank you for your reply! I have succeed in writting a fidelity cost function by the help of the demo you mentioned in your reply. Here is the feasible code below:

@qml.qnode(dev)
def circuit(params, ini_state, tar_state):
    qml.StatePrep(ini_state, range(N))
    SU8(params, range(N))
    return qml.expval(qml.Hermitian(qml.math.dm_from_state_vector(tar_state), wires=range(N)))

def costf(params, ini_state_batch, tar_state_batch):
    cost=0
    for i in range(len(ini_state_batch)):
        cost+=1-circuit(params, ini_state_batch[i], tar_state_batch[i])
    return cost/len(ini_state_batch)

I find some rules to ensure the codes run correctly:
· The qml.qnode(dev) statement should only appears before the function defining the whole circuit (circuit() in my case). It should not appear before the component circuits (for example, SU8() in my case) and the cost function.
· The component circuits should not return anything.
· The fidelity should be measured out instead of calculating from the final state vector or density matrix.

Let me explain that. Multiple qml.qnode(dev) appearing in the same code will cause errors, so it should appear only once. If qml.qnode(dev) is not put before the whole circuit (for example, put before the cost function):

def circuit(params, ini_state, tar_state):
    qml.StatePrep(ini_state, range(N))
    SU8(params, range(N))
    return qml.expval(qml.Hermitian(qml.math.dm_from_state_vector(tar_state), wires=range(N)))

@qml.qnode(dev)
def costf(params, ini_state_batch, tar_state_batch):
    cost=0
    for i in range(len(ini_state_batch)):
        cost+=1-circuit(params, ini_state_batch[i], tar_state_batch[i])
    return cost/len(ini_state_batch)

The type of the returned object from circuit() will be ExpectationMP rather than a value. It is not supported by any operation including the most basic +, - and *, causing the errors like unsupported operand type(s) for -: 'int' and 'ExpectationMP'. If we return qml.state() or qml.density_matrix(range(N)) in circuit(), The type of the returned object will be StateMP and DensityMP instead of the expected numpy array, causing errors too.

The reason why the component circuits should not return anything is the same. It will prevent us from getting unwanted StateMP, DensityMP and ExpectationMP.

If we try to get the final state vector or the density matrix from the simulation and calculate the fidelity manually, for example:

@qml.qnode(dev)
def circuit(params, ini_state, tar_state):
    qml.StatePrep(ini_state, range(N))
    SU8(params, range(N))
    return qml.state()

def costf(params, ini_state_batch, tar_state_batch):
    cost=0
    for i in range(len(ini_state_batch)):
        final_state=circuit(params, ini_state_batch[i], tar_state_batch[i])
        cost+=1-np.abs(np.dot(final_state, tar_state_batch[i]))**2
    return cost/len(ini_state_batch)

It will cause the following error: ValueError: Computing the gradient of circuits that return the state with the parameter-shift rule gradient transform is not supported, as it is a hardware-compatible method.. If we try to return the density matrix in circuit() and use the qml.math.fidelity() to calculate the fidelty, the same error will occur.

4

Hi @Tz_19 , thank you for sharing your learnings.

You’re right in the use of @qml.qnode. It shouldn’t be used everywhere, in fact it should only be used when you have a quantum circuit that returns a measurement! You can learn more about qnodes in the Quantum circuits section of the documentation or in the PennyLane Fundamentals module of the Codebook. In fact I strongly recommend that you explore this Codebook module since it might teach you a lot on how to use PennyLane and save you a lot of time debugging!

Finally, if you’re interested in differentiation and optimization you might like the Gradients and training section of the documentation.

I hope you find these resources useful. Enjoy using PennyLane!