Add derivatives for PQK

[RobertoFlorez]

Hi!

As discussed in the meeting, here is the implementation of the single-variable RBF PQK derivatives.

It includes 3 relevant additions:

1. evaluate_string_derivatives: given an evaluation string (i.e "dKdx", ...) returns the corresponding matrix
2. evaluate_derivatives: calls evaluate_string_derivative, checks caching and returns a dictionary (loosely based on LowLevelPennylane.evaluate )
3. GaussianOuterKernel.dKdx, GaussianOuterKernel.dKdy, GaussianOuterKernel.dKdxdx, GaussianOuterKernel.dKdxdy : the analytical derivatives of the RBF

About notation:

In the implementation of evaluate_string_derivatives, O correspond to the array of expectation values of the observables that go into the PQK, in the squlearn documentation of the PQK this is refered to as QNN(x). Perhaps, a better name instead of O would be good :) (pylint also did not like this name)

About the correctness:

I have numerically benchmarked with an analytical example evaluate_derivatives for dKdx, dKdy, dKdxdx, dKdxdy, so they should be correctly implemented. I have not benchmarked yet dKdp (as I have not used it for any ODE). In summary, the status of the implemented derivatives is:

1D variable:

• Implemented and numerically benchmarked: dKdx, dKdy, dKdxdx
• Implemented, not yet numerically benchmarked: dKdp
• Not Implemented: dKdop, (more ?)

Multidimensional variable:

• Implemented and numerically benchmarked: dKdx, dKdy
• Not Implemented: dKdxdx, dKdxdy, dKdop, dKdp (more ?)

I am very happy and grateful to receive your feedback to improve the implementation. 🙏 🚀

Thank you very much for your help!

---

---

In case someone wants to double check the derivation, see screenshot.

[David-Kreplin]

Thank you very much, most of the stuff looks alright! I commented on a few code parts that are a little redundant. Especially, the dKdy part is basically the same as the dKdx part except for the sign, so it can be combined into a single code block. Also, you could merge the two functions since the evaluate_string_derivatives function is not really needed. Could you also test the dKdp derivative?

[RobertoFlorez]

Hi!

Again thank you very much for the feedback!

I implemented the requested changes:

• Removed itertools: now nested for-loops
• Merged evaluate_derivatives and evaluate_derivatives_string;
• Merged dKdx and dKdy
• Included py_tests that compare analytical derivatives with numerical derivatives for single and multi variable differentiation (dKdx, dKdy, dKdp)

---

The tests use an analytical expression for the PQK kernel. If there is interest in deriving this equations, see here:

import sympy as sp
from sympy.physics.quantum import TensorProduct

x = sp.symbols("x0", real=True)
y = sp.symbols("y0", real=True)

x1 = sp.symbols("x1", real=True)
y1 = sp.symbols("y1", real=True)

p = sp.symbols("p0", real=True)

p1 = sp.symbols("p1", real=True)


gamma = sp.symbols("gamma", real=True)


def RX(x):
    return sp.Matrix([[sp.cos(x/2), -1j*sp.sin(x/2)], [-1j*sp.sin(x/2), sp.cos(x/2)]])
def RY(y):
    return sp.Matrix([[sp.cos(y/2), -sp.sin(y/2)], [sp.sin(y/2), sp.cos(y/2)]])
def RZ(z):
    return sp.Matrix([[sp.exp(-1j*z/2), 0], [0, sp.exp(1j*z/2)]])


X = sp.Matrix([[0, 1], [1, 0]])
Y = sp.Matrix([[0, -1j], [1j, 0]])
Z = sp.Matrix([[1, 0], [0, -1]])

ket0 = sp.Matrix([[1], [0]])
bra0 = ket0.T


def U_gate(x):
    return RX(x)@RY(p)


def U_gate_2d(x, x1):
    return TensorProduct(RX(x), RX(x1)) @ TensorProduct(RY(p), RY(p1))

rhox_2d = U_gate_2d(x, x1) @ TensorProduct(ket0, ket0) @ TensorProduct(bra0, bra0) @ U_gate_2d(x, x1).H
rhoy_2d = U_gate_2d(y, y1) @ TensorProduct(ket0, ket0) @ TensorProduct(bra0, bra0) @ U_gate_2d(y, y1).H

rhox = U_gate(x) @ ket0 @ bra0 @ U_gate(x).H
rhoy = U_gate(y) @ ket0 @ bra0 @ U_gate(y).H

#Fidelity QK
def FQK(rho1, rho2):
    return sp.trace(rho1 @ rho2)



def PQK(rho1, rho2, gamma):
    def create_observable(n_qubits, operator):
        identity_list = [sp.eye(2) for _ in range(n_qubits)]
        observable_list = []
        for i in range(n_qubits):
            observable_list.append(TensorProduct(*identity_list[:i], operator, *identity_list[i+1:]))
        return observable_list
    
    n_qubits = int(sp.log(rho1.shape[0], 2))
    observable_list = create_observable(n_qubits, X) + create_observable(n_qubits, Y) + create_observable(n_qubits, Z)
    O_sum_x = sum([sp.trace(observable @ rho1) for observable in observable_list])
    O_sum_y = sum([sp.trace(observable @ rho2) for observable in observable_list])
    O_sum = 0
    for observable in observable_list:
        O_sum += ((sp.trace(observable @ rho1).simplify() - sp.trace(observable @ rho2).simplify())**2)

    return sp.exp(-gamma*O_sum)

print(PQK(rhox, rhoy, gamma).simplify())
print(PQK(rhox_2d, rhoy_2d, gamma).simplify())
print(FQK(rhox, rhoy).simplify())

[David-Kreplin]

I did a little bit of simplifying the code. Is ready now.