# Differential equation resolution

## Introduction

We present here a Perceval implementation of a Quantum Machine Learning algorithm for solving differential equations. Its aims is to approximate the solution to the differential equation considered in :

$\frac{d f}{d x}+\lambda f(x)(\kappa+\tan (\lambda x))=0$

with boundary condition $$f(0)=f_{0}$$. The analytical solution is $$f(x)=f_0\exp (-\kappa \lambda x) \cos (\lambda x)$$.

### QML Loss Function Definition

In order to use QML to solve this differential equation, we first need to derive from it a loss function whose minimum is associated to its analytical solution.

Let $$F\left[\left\{d^{m} f / d x^{m}\right\}_{m},f, x\right]=0$$ be a general differential equation verified by $$f(x)$$, where $$F[.]$$ is an operator acting on $$f(x)$$, its derivatives and $$x$$. For the solving of a differential equation, the loss function described in  consists of two terms

$\mathcal{L}_{\boldsymbol{\theta}}\left[\left\{d^{m} g / d x^{m}\right\}_{m},g, x\right]:=\mathcal{L}_{\boldsymbol{\theta}}^{(\mathrm{diff})}\left[\left\{d^{m} g / d x^{m}\right\}_{m},g, x\right]+\mathcal{L}_{\boldsymbol{\theta}}^{(\text {boundary})}[g, x].$

The first term $$\mathcal{L}_{\boldsymbol{\theta}}^{(\mathrm{diff})}$$ corresponds to the differential equation which has been discretised over a fixed regular grid of $$M$$ points noted $$x_i$$:

$\mathcal{L}_{\boldsymbol{\theta}}^{(\mathrm{diff})}\left[\left\{d^{m} g / d x^{m}\right\}_{m},g, x\right]:=\frac{1}{M} \sum_{i=1}^{M} L\left(F\left[d_{x}^m g\left(x_{i}\right), g\left(x_{i}\right), x_{i}\right], 0\right),$

where $$L(a,b) := (a - b)^2$$ is the squared distance between two arguments. The second term $$\mathcal{L}_{\boldsymbol{\theta}}^{(\text {boundary })}$$ is associated to the initial conditions of our desired solution. It is defined as:

$\mathcal{L}_{\boldsymbol{\theta}}^{\text {(boundary) }}[g, x]:=\eta L\left(g(x_0), f_{0}\right),$

where $$\eta$$ is the weight granted to the boundary condition and $$f_{0}$$ is given by $$f(x_0) = f_0$$.

Given a function approximator $$f^{(n)}(x, \boldsymbol{\theta}, \boldsymbol{\lambda})$$, the loss function above will be minimised using a classical algorithm, updating the parameters $$\boldsymbol{\theta}$$ based on samples obtained using a quantum device.

### Quantum circuit architecture

The feature map used is presented in [2,3,4]. The quantum circuit architecture from  is expressed as $$\mathcal{U}(x, \boldsymbol{\theta}):=\mathcal{W}^{(2)}\left(\boldsymbol{\theta}_{2}\right) \mathcal{S}(x) \mathcal{W}^{(1)}\left(\boldsymbol{\theta}_{1}\right).$$ The phase-shift operator $$\mathcal{S}(x)$$ incorporates the $$x$$ dependency of the function we wish to approximate. It is sandwiched between two universal interferometers $$\mathcal{W}^{(1)}(\boldsymbol{\theta_1})$$ and $$\mathcal{W}^{(2)}(\boldsymbol{\theta_2})$$, where the beam-splitter parameters $$\boldsymbol{\theta_1}$$ and $$\boldsymbol{\theta_2}$$ of this mesh architecture are tunable to enable training of the circuit. The output measurement operator, noted $$\mathcal{M}(\boldsymbol{\lambda})$$, is the projection on the Fock states obtained using photon-number resolving detectors, multiplied by some coefficients $$\boldsymbol{\lambda}$$ which can also be tunable. Formally, we have:

$\mathcal{M}(\boldsymbol{\lambda}) = \sum_{\mathbf{\left | n^{(f)}\right \rangle}}\lambda_{\mathbf{\left | n^{(f)}\right \rangle}}\mathbf{\left | n^{(f)}\right \rangle}\mathbf{\left \langle n^{(f)}\right |},$

where the sum is taken over all $$\binom{n+m-1}{n}$$ possible Fock states considering $$n$$ photons in $$m$$ modes. Let $$\mathbf{\left | n^{(i)}\right \rangle} = \left |n^{(i)}_1,n^{(i)}_2,\dots,n^{(i)}_m\right \rangle$$ be the input state consisting of $$n$$ photons where $$n^{(i)}_j$$ is the number of photons in input mode $$j$$. Given these elements, the circuit’s output $$f^{(n)}(x, \boldsymbol{\theta}, \boldsymbol{\lambda})$$ is given by the following expectation value:

$f^{(n)}(x, \boldsymbol{\theta}, \boldsymbol{\lambda})=\left\langle\mathbf{n}^{(i)}\left|\mathcal{U}^{\dagger}(x, \boldsymbol{\theta}) \mathcal{M}(\boldsymbol{\lambda}) \mathcal{U}(x, \boldsymbol{\theta})\right| \mathbf{n}^{(i)}\right\rangle.$

This expression can be rewritten as the following Fourier series 

$f^{(n)}(x, \boldsymbol{\theta}, \boldsymbol{\lambda})=\sum_{\omega \in \Omega_{n}} c_{\omega}(\boldsymbol{\theta}, \boldsymbol{\lambda}) e^{i \omega x},$

where $$\Omega_n = \{-n, -n+1, \dots, n-1, n \}$$ is the frequency spectrum one can reach with $$n$$ incoming photons and $$\{c_\omega(\boldsymbol{\theta}, \boldsymbol{\lambda})\}$$ are the Fourier coefficients. The $$\boldsymbol{\lambda}$$ parameters are sampled randomly in the interval $$[-a;a]$$, with $$a$$ a randomly chosen integer. $$f^{(n)}(x, \boldsymbol{\theta}, \boldsymbol{\lambda})$$ will serve as a function approximator for this chosen differential equation. Differentiation in the loss function is discretised as $$\frac{df}{dx} \simeq \frac{f(x+\Delta x) - f(x-\Delta x)}{2\Delta x}$$.

$$n, m,$$ and $$\boldsymbol{\lambda}$$ are variable parameters defined below. $$\Delta x$$ is the mesh size.

## Perceval Simulation

### Initialisation

:

import perceval as pcvl
import numpy as np
from math import comb
from scipy.optimize import minimize
import time
import matplotlib.pyplot as plt
import matplotlib as mpl
import tqdm as tqdm


We will run this notebook with 4 photons. We could use more photons, but the result with 4 photons is already satisfying.

:

nphotons = 4


### Differential equation parameters

We define here the value of the differential equation parameters and boundary condition $$\lambda, \kappa, f_0$$.

:

# Differential equation parameters
lambd = 8
kappa = 0.1

def F(u_prime, u, x):       # DE, works with numpy arrays
return u_prime + lambd * u * (kappa + np.tan(lambd * x))

:

# Boundary condition (f(x_0)=f_0)
x_0 = 0
f_0 = 1

:

# Modeling parameters
n_grid = 50    # number of grid points of the discretized differential equation
range_min = 0  # minimum of the interval on which we wish to approximate our function
range_max = 1  # maximum of the interval on which we wish to approximate our function
X = np.linspace(range_min, range_max-range_min, n_grid)  # Optimisation grid

:

# Differential equation's exact solution - for comparison
def u(x):
return np.exp(- kappa*lambd*x)*np.cos(lambd*x)

:

# Parameters of the quantum machine learning procedure
N = nphotons              # Number of photons
m = nphotons              # Number of modes
eta = 5                   # weight granted to the initial condition
a = 200                   # Approximate boundaries of the interval that the image of the trial function can cover
fock_dim = comb(N + m - 1, N)
# lambda coefficients for all the possible outputs
lambda_random = 2 * a * np.random.rand(fock_dim) - a
# dx serves for the numerical differentiation of f
dx = (range_max-range_min) / (n_grid - 1)

:

# Input state with N photons and m modes
input_state = pcvl.BasicState(*N+*(m-N))
print(input_state)

|1,1,1,1>


## Definition of the circuit

We will generate a Haar-random initial unitary using QR decomposition built in Perceval Matrix.random_unitary, the circuit is defined by the combination of 3 sub-circuits - the intermediate phase is a parameter.

:

"Haar unitary parameters"
# number of parameters used for the two universal interferometers (2*m**2 per interferometer)
parameters = np.random.normal(size=4*m**2)

px = pcvl.P("px")
c = pcvl.Unitary(pcvl.Matrix.random_unitary(m, parameters[:2 * m ** 2]), name="W1")\
// (0, pcvl.PS(px))\
// pcvl.Unitary(pcvl.Matrix.random_unitary(m, parameters[2 * m ** 2:]), name="W2")

simulator_backend = pcvl.BackendFactory().get_backend("SLOS")
s1 = simulator_backend(pcvl.Matrix.random_unitary(m))
s1.compile(input_state)

pcvl.pdisplay(c)

: ### Expectation value and loss function computation

The expectation value of the measurement operator $$\mathcal{M}(\boldsymbol{\lambda})$$ is obtained directly from Fock state probabilities computed by Perceval. Given this expectation value, the code snippet below computes the loss function defined in the Introduction.

Note the use of the all_prob simulator method giving directly access to the probabilities of all possible output states, including null probabilities. This calculation is optimized in SLOS backend.

:

def computation(params):
global current_loss
global computation_count
"compute the loss function of a given differential equation in order for it to be optimized"
computation_count += 1
f_theta_0 = 0  # boundary condition
coefs = lambda_random  # coefficients of the M observable
# initial condition with the two universal interferometers and the phase shift in the middle
U_1 = pcvl.Matrix.random_unitary(m, params[:2 * m ** 2])
U_2 = pcvl.Matrix.random_unitary(m, params[2 * m ** 2:])

px = pcvl.P("x")
c = pcvl.Unitary(U_2) // (0, pcvl.PS(px)) // pcvl.Unitary(U_1)

px.set_value(np.pi * x_0)
U = c.compute_unitary(use_symbolic=False)
s1.U = U
f_theta_0 = np.sum(np.multiply(s1.all_prob(input_state), coefs))

# boundary condition given a weight eta
loss = eta * (f_theta_0 - f_0) ** 2 * len(X)

# Y is before the domain we are interested in (used for differentiation), x_0 is at Y
Y = np.zeros(n_grid + 2)

# x_0 is at the beginning of the domain, already calculated
Y = f_theta_0

px.set_value(np.pi * (range_min - dx))
s1.U = c.compute_unitary(use_symbolic=False)
Y = np.sum(np.multiply(s1.all_prob(input_state), coefs))

for i in range(1, n_grid):
x = X[i]
px.set_value(np.pi * x)
s1.U = c.compute_unitary(use_symbolic=False)
Y[i + 1] = np.sum(np.multiply(s1.all_prob(input_state), coefs))

px.set_value(np.pi * (range_max + dx))
s1.U = c.compute_unitary(use_symbolic=False)
Y[n_grid + 1] = np.sum(np.multiply(s1.all_prob(input_state), coefs))

# Differentiation
Y_prime = (Y[2:] - Y[:-2])/(2*dx)

loss += np.sum((F(Y_prime, Y[1:-1], X))**2)

current_loss = loss / len(X)
return current_loss


### Classical optimisation

Finally the code below performs the optimisation procedure using the loss function defined in the previous section. To this end, we use a Broyden–Fletcher–Goldfarb–Shanno (BFGS) optimiser  from the SciPy library.

:

def callbackF(parameters):
"""callback function called by scipy.optimize.minimize allowing to monitor progress"""
global current_loss
global computation_count
global loss_evolution
global start_time
now = time.time()
pbar.set_description("M= %d Loss: %0.5f #computations: %d elapsed: %0.5f" %
(m, current_loss, computation_count, now-start_time))
pbar.update(1)
loss_evolution.append((current_loss, now-start_time))
computation_count = 0
start_time = now

:

computation_count = 0
current_loss = 0
start_time = time.time()
loss_evolution = []

pbar = tqdm.tqdm()
res = minimize(computation, parameters, callback=callbackF, method='BFGS', options={'gtol': 1E-2})

M= 4 Loss: 0.00259 #computations: 195 elapsed: 1.53568: : 117it [01:43,  1.19s/it]


After the optimisation procedure has been completed, the optimal unitary parameters (in res.x) can be used to determine the quantum circuit beam-splitter and phase-shifter angles for an experimental realisation.

:

print("Unitary parameters", res.x)

Unitary parameters [ 0.06288015 -0.46766255 -1.35750088 -0.69245981  1.91915527 -0.3664978
2.33459189 -0.81260598  0.00469419 -0.24181104  0.96650738  1.87575626
-0.42499621 -0.46514721 -0.40457874 -2.42729821  0.1573016  -0.73709542
0.61004258  1.84471732  1.73285762 -0.20159712  0.25432091  1.03383559
-0.30688057  2.90045031  0.27206763  0.82699894 -0.27334301 -1.01200369
-1.58195825  1.76439029 -2.74254204  0.69107417  1.48317302 -0.43958469
-0.90807896  1.53240125  0.86739936  0.33513812 -0.82602653 -0.88284405
0.44669912  0.37906319 -0.03029149  0.14336839 -0.17139438 -2.05323973
-0.11085795 -0.09207871 -0.08488444 -0.36721392 -0.66922089  0.35282081
0.53744991 -0.13106399 -0.68348908  0.65520211 -1.04250027  1.02240978
0.55371764 -0.85351974 -0.10360431  0.64527886]


### Plotting the approximation

We now plot the result of our optimisation in order to compare the QML algorithm’s output and the analytical solution.

:

def plot_solution(m, N, X, optim_params, lambda_random):
Y = []
U_1 = pcvl.Matrix.random_unitary(m, optim_params[:2 * m ** 2])
U_2 = pcvl.Matrix.random_unitary(m, optim_params[2 * m ** 2:])
px = pcvl.P("x")
c = pcvl.Unitary(U_2) // (0, pcvl.PS(px)) // pcvl.Unitary(U_1)

for x in X:
px.set_value(np.pi * x)
U = c.compute_unitary(use_symbolic=False)
s1.U = U
f_theta = np.sum(np.multiply(s1.all_prob(input_state), lambda_random))
Y.append(f_theta)
exact = u(X)
plt.plot(X, Y, label="Approximation with {} photons".format(N))

:

X = np.linspace(range_min, range_max, 200)

# Change the plot size
default_figsize = mpl.rcParamsDefault['figure.figsize']
mpl.rcParams['figure.figsize'] = [2 * value for value in default_figsize]

plot_solution(m, N, X, res.x, lambda_random)

plt.plot(X, u(X), 'r', label='Analytical solution')
plt.legend()
plt.show() :

plt.plot([v for v in loss_evolution])
plt.yscale("log")
plt.xlabel("Number of epochs")
plt.ylabel("Loss function value")

:

Text(0, 0.5, 'Loss function value') 