Content of Quantum Monte Carlo page

[Den]

The Pennylane page on Quantum Monte Carlo shows an integral without dx

.

In addition, using Pennylane’s provided code, mu is obtained as ~ 0.43, as is on the webpage, but using scipy.integrate.simpson to integrate the product of numpy arrays f(x)p(x) on the same interval (-pi, pi) gives mu ~ 0.087, not ~0.43.
Increasing the x interval to (-10pi, 10*pi), mu becomes ~1.43 using scipy’s simpson and numpy’s trapz functions.

Am I missing something?

The code is like this:

from scipy.stats import norm
from scipy.integrate import simpson
from pennylane import numpy as np
import numpy as np

m    = 5
M    = 2 ** m

xmax = np.pi
xs   = linspace(-xmax, xmax, M)     # Domain

def p(xs):
    "Probability distribution"
    ps1   = np.array([norm().pdf(x) for x in xs]) # Distribution
    ps1   /= np.sum(ps1)    # Normalization
    return ps1

def f(xs):
    "Function to be integrated with distribution"
    return sin(xs) ** 2

simpson(  f(xs)* p(xs), xs)
# Result is ~ 0.087

np.trapz( f(xs)* p(xs), xs)
# Result is ~ 0.087

I’m using python 3.9.13.

[kevinkawchak]

Hello @Den , Welcome to the community. I was able to get the 0.43 from documentation, but receive a linspace error for the code you provided.

[Den]

My bad, I’m used to call the functions straight from the package.
Fixed. Now this code should compile properly.

from scipy.stats import norm
from scipy.integrate import simpson
from pennylane import numpy as np

m    = 5
M    = 2 ** m

xmax = np.pi
xs   = np.linspace(-xmax, xmax, M)     # Domain

def p(xs):
    "Probability distribution"
    ps1   = np.array([norm().pdf(x) for x in xs]) # Distribution
    ps1   /= np.sum(ps1)    # Normalization
    return ps1

def f(xs):
    "Function to be integrated with distribution"
    return np.sin(xs) ** 2

simpson(  f(xs)* p(xs), xs)
# Result is ~ 0.087

np.trapz( f(xs)* p(xs), xs)
# Result is ~ 0.087

[DanielNino27]

Hi @Den,
Welcome to the forum! The cause for the discrepancy between the classical integration and the quantum Monte Carlo example is a little subtle. In your classical integration example, you are re-normalizing an already normalized distribution. The function norm.pdf(x) samples from a normalized standard distribution. So when you call

ps1   /= np.sum(ps1) 

You are normalizing an already normalized distribution. Try removing that line from your classical integration and you should get the correct value.

This begs the question: why does the code in the quantum Monte Carlo function do that? In the quantum Monte Carlo algorithm, you need to embed the probability distribution on a circuit. So it needs to be a discretized version of whatever continuous function you’re looking at and, in this case, it needs to be renormalized over the interval you are sampling.

Hope that helps!

[Den]

You’re right, thanks.
Now the results match.

[DanielNino27]

That’s great!
We are interested in collecting feedback from users on how we can keep improving PennyLane so that we can keep implementing great features. We have a new PennyLane survey that you may be interested in if you would like to provide some feedback on your experience using PennyLane.

[Den]

Two doubts appeared as extension of this conversation:

• Is it allowed for the probability p(x) to be an arbitrary one with arbitrary integration limits as long as p(x) is normalized in the integration interval?
  I tested p(x) as a normalized polynomial ( one that satisfies \int_a^b \mathrm{d} x \, p(x) = 1 between a and b), but the QMC result is different from the classical integration (code below).

from scipy.stats import norm
from scipy.integrate import simpson
from pennylane import numpy as np
import pennylane as qml


def qmci(f, p, xs, m, n):
    N = 2 ** n
    
    tgt_wires = range(m + 1)
    est_wires = range(m + 1, n + m + 1)
    dev       = qml.device("default.qubit", wires=(n + m + 1))

    ps = p(xs)
    ps /= np.sum(p(xs))
    
    @qml.qnode(dev)
    def circuit(ps, f, xs):
        qml.templates.QuantumMonteCarlo(
            ps,
            f,
            target_wires     = tgt_wires,
            estimation_wires = est_wires,
        )
        return qml.probs(est_wires)

    ph_est = np.argmax(circuit(ps, f, xs)[:int(N / 2)]) / N
    res    = (1 - np.cos(np.pi * ph_est)) / 2
    return res

n    = 10
m    = 5
M    = 2 ** m
xmax = np.pi
xs   = np.linspace(-xmax, xmax, M)     # Domain


# DISTRIBUTIONS
# Polynomial A normalized under integration limits
rdn = np.arange(0.1,0.4,0.1)     # Arbitrary sequence of polynomial coefficients
A = lambda xs: (rdn[2]*xs**2 + rdn[1]**xs + rdn[0] )/104.49637640566247

# Normal distribution for comparison
p = lambda xs: norm().pdf(xs)

# INTEGRATION FUNCTION
f = lambda i: np.sin(xs[i]) ** 2 # Used for QMC integration
F = lambda xs: np.sin(xs) ** 2   # Used for classical integration


qmci( f, A, xs, m, n) 
# Result is 0.2670117521160169

simpson( F(xs)* A(xs), xs)
# Result is 0.31218885373330574

• Also, I noticed that the classical and QMC integrations diverge when the integration limits are set as (-10 \pi, 10\pi) with the normal distribution as is in the documentation (code below for copy and paste).
• Is there a way to know the appropriate integration limits of convergence?

from scipy.stats import norm
from scipy.integrate import simpson
from pennylane import numpy as np
import pennylane as qml


def qmci(f, p, xs, m, n):
    N = 2 ** n
    
    tgt_wires = range(m + 1)
    est_wires = range(m + 1, n + m + 1)
    dev       = qml.device("default.qubit", wires=(n + m + 1))

    ps = p(xs)
    ps /= np.sum(p(xs))
    
    @qml.qnode(dev)
    def circuit(ps, f, xs):
        qml.templates.QuantumMonteCarlo(
            ps,
            f,
            target_wires     = tgt_wires,
            estimation_wires = est_wires,
        )
        return qml.probs(est_wires)

    ph_est = np.argmax(circuit(ps, f, xs)[:int(N / 2)]) / N
    res    = (1 - np.cos(np.pi * ph_est)) / 2
    return res

n    = 10
m    = 5
M    = 2 ** m
xmax = 10*np.pi
xs   = np.linspace(-xmax, xmax, M)     # Domain


# DISTRIBUTIONS
# Polynomial A normalized under integration limits
#rdn = np.arange(0.1,0.4,0.1)     # Arbitrary sequence of polynomial coefficients
#A = lambda xs: (rdn[2]*xs**2 + rdn[1]**xs + rdn[0] )/104.49637640566247

# Normal distribution for comparison
p = lambda xs: norm().pdf(xs)

# INTEGRATION FUNCTION
f = lambda i: np.sin(xs[i]) ** 2 # Used for QMC integration
F = lambda xs: np.sin(xs) ** 2   # Used for classical integration


qmci( f, p, xs, m, n) 
# Result is 0.29178521995118134

simpson( F(xs)* p(xs), xs)
# Result is 0.6971081255327303

[isaacdevlugt]

Hey @Den, thanks for the follow up questions . I’ll answer each question individually:

Is it allowed for the probability p(x) to be an arbitrary one with arbitrary integration limits as long as p(x) is normalized in the integration interval?

Sort of, but some details are missing. The probability distribution must be normalized under its entire domain (i.e., \int_{\mathcal{D}} p(x) dx = 1) when you use something like Simpsons rule in this context, and for the QMC routine the distribution it must be normalized just in the region you’re integrating over.

I tested p(x) as a normalized polynomial ( one that satisfies ∫badxp(x)=1 between a and b), but the QMC result is different from the classical integration (code below).

There are some problems with this. Firstly, the function A(x) = 0.1 + 0.2x + 0.3x^2 is not a normalized function over it’s entire domain, and therefore it’s not a probability distribution. While the QMC routine will be tricked into thinking that it indeed is a valid probability distribution because you’re only showing it in a smaller interval and you’ve normalized it with ps /= np.sum(p(xs)), comparing it to Simpson’s rule is not apples-to-apples.

Also, I noticed that the classical and QMC integrations diverge when the integration limits are set as (−10π,10π) with the normal distribution as is in the documentation (code below for copy and paste).

I think this is just a result of a typo in your code. You have A(x) = 0.1 + 0.2^x + 0.3x^2

A = lambda xs: (rdn[2]*xs**2 + rdn[1]**xs + rdn[0] )/104.49637640566247

Is there a way to know the appropriate integration limits of convergence?

This really depends on the form of p(x) and the function f(x) . In the case of p being a standard normal distribution, p(x) \rightarrow 0 when x \rightarrow \pm \infty, so you can safely integrate within a finite region

Hope this helps! Monte Carlo simulations can be nuanced and involve tricky numerics.