Quantum monte carlo for different functions except sin^2x

I am working with the example given in qml.quantum_monte_carlo — PennyLane 0.40.0 documentation . This integrates sin^2x.

But if I replace, sin^2x with cos^2x that is anouther function bounded between [0,1] then the example does not give correct result

from scipy.stats import norm, uniform
import numpy as np
m = 5
M = 2 ** m

xmax = 1  # bound to region [-pi, pi]
xs = np.linspace(0, xmax, M)

probs = np.array([uniform().pdf(x) for x in xs])
probs /= np.sum(probs)

func = lambda i: np.cos(xs[i])**2
r_rotations = np.array([2 * np.arcsin(np.sqrt(func(i))) for i in range(M)])
from pennylane.templates.state_preparations.mottonen import (
    _apply_uniform_rotation_dagger as r_unitary,
)
import pennylane as qml
n = 6
N = 2 ** n

a_wires = range(m)
wires = range(m + 1)
target_wire = m
estimation_wires = range(m + 1, n + m + 1)

dev = qml.device("default.qubit", wires=(n + m + 1))

def fn():
    qml.templates.MottonenStatePreparation(np.sqrt(probs), wires=a_wires)
    r_unitary(qml.RY, r_rotations, control_wires=a_wires[::-1], target_wire=target_wire)

@qml.qnode(dev)
def qmc():
    qml.quantum_monte_carlo(fn, wires, target_wire, estimation_wires)()
    return qml.probs(estimation_wires)

phase_estimated = np.argmax(qmc()[:int(N / 2)]) / N
(1 - np.cos(np.pi * phase_estimated)) / 2

What changes should I do in the code to make this work?

Hi @Rochisha_Agarwal , welcome to the Forum!

I’m not sure but I’m assuming you’ll need to change the r_rotations since at the moment you’re using the arcsine, where you might need the arccosine instead.

On the other hand your xs should probably go from -pi/2 to pi/2 so that you can complete one full period of the function.

For the estimated value I don’t know if the formula for \mu holds. You may need to check the paper for this.

Finally, I don’t think there’s any guarantee that this will work for every function. You’ll have to check the assumptions and conditions where this is expected to work in case my suggestions above are not enough.

I hope this helps though, and let me know if the suggestions helped you solve the issue!

Hello, changing arcsin to arccosine doesn’t solve the problem. This code also doesn’t work if I change the function np.sin^2(x) → np.sin^2(x-np.pi/2). This looks like a bigger issue.

The condition for the function is that it should be between [0, 1] for the entire domain. This satisfies for cos^2x.

Hi @Rochisha_Agarwal ,

Have you checked whether the error you get is between the expected margin of error? If not feel free to check the paper. There are some equations regarding this in several parts of the paper. There may be other conditions or limitations mentioned in the paper. For example it’s possible that the model isn’t expressive enough to capture the changes you’re trying to test. I’m not really sure since I’m not an author of the paper. I noticed that the paper has the email of the first author so you can try contacting Patrick Rebentrost.

Just to rule out any code issues I tested a different function and the result was correct. If you change the function to np.sin^2(0.1x) you’ll see that the output closely resembles the correct output of 0.0099. So my guess is that the limitation is in the model itself.