<jupyter_start><jupyter_text>Differential equation resolution IntroductionWe 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 \[1\]:$$\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)=\exp (-\kappa \lambda x) \cos (\lambda x)+ f_0 - 1$. QML Loss Function DefinitionIn 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 \[1\] 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 architectureThe feature map used is presented in \[2,3,4\]. The quantum circuit architecture from \[4\] 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 \[4\]$$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{\Delta f}{\Delta x}$.$n, m, \boldsymbol{\lambda}$ and $\Delta x $ are variable parameters defined below.  Perceval Simulation Initialisation<jupyter_code>from IPython import display
import perceval as pcvl
import perceval.lib.phys as phys
import numpy as np
from math import comb
from scipy.optimize import minimize
import time
import matplotlib.pyplot as plt
import tqdm as tqdm<jupyter_output><empty_output><jupyter_text>We will run this notebook with 6 photons<jupyter_code>nphotons = 6<jupyter_output><empty_output><jupyter_text>Differential equation parametersWe define here the value of the differential equation parameters and boundary condition $\lambda, \kappa, f_0$.<jupyter_code># Differential equation parameters
lambd = 8
kappa = 0.1

def F(u_prime, u, x):       # DE
    return (u_prime(x) + lambd * u(x) * (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) / 1000
# Input state with N photons and m modes
input_state = pcvl.BasicState([1]*N+[0]*(m-N))
print(input_state)<jupyter_output>|1,1,1,1,1,1><jupyter_text>Definition of the circuitWe 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.<jupyter_code>"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)

simulator_backend = pcvl.BackendFactory().get_backend("SLOS")

px = pcvl.P("px")
c = (pcvl.Circuit(m, pcvl.Matrix.random_unitary(m, parameters[:2 * m ** 2]), name="W1")
     // (0, phys.PS(px))
     // pcvl.Circuit(m, pcvl.Matrix.random_unitary(m, parameters[2 * m ** 2:]), name="W2"))
pcvl.pdisplay(c)
s1 = simulator_backend(pcvl.Matrix.random_unitary(m))
s1.compile(input_state)<jupyter_output><empty_output><jupyter_text>Expectation value and loss function computationThe 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. This calculation is optimized in SLOS backend.<jupyter_code>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.Circuit(m, U_2) // (0, phys.PS(px)) // pcvl.Circuit(m, U_1)

    px.set_value(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)

    def ldiff(x):
        i, j = 0, 0
        px.set_value(x)
        s1.U = c.compute_unitary()
        f_theta = np.sum(np.multiply(s1.all_prob(input_state), coefs))
        px.set_value(x+dx)
        s1.U = c.compute_unitary()
        f_theta_dx = np.sum(np.multiply(s1.all_prob(input_state), coefs))
        diff_f_theta = (f_theta_dx - f_theta) / dx
        return (diff_f_theta + lambd * (kappa + np.tan(lambd * x)) * f_theta) ** 2

    for x in X:
        loss += ldiff(x)

    current_loss = loss / len(X)
    return current_loss<jupyter_output><empty_output><jupyter_text>Classical optimisationFinally 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 \[5\] from the SciPy library.<jupyter_code>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})<jupyter_output>M= 6 Loss: 0.00969 #computations: 145 elapsed: 0.71338: : 243it [04:33,  1.07s/it]<jupyter_text>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.<jupyter_code>print("Unitary parameters", res.x)<jupyter_output>Unitary parameters [-2.76520236e+00 -1.91575408e+00 -1.37475078e+00  8.86247487e-01
 -6.39263138e-01 -3.69530446e+00 -4.55679771e-01  3.86299771e-01
 -2.07979490e-01 -4.33443151e-01  4.81646341e-01 -1.43837756e+00
 -3.44653481e+00  4.99906775e-01 -1.46893144e+00  2.15232837e+00
  1.23956750e+00 -5.98775125e-02  9.08920960e-01  1.19851003e+00
  7.69252675e-01 -5.51722277e-01  6.99316491e-01 -2.55988798e-01
 -2.81322628e+00 -3.83386441e+00 -4.84602436e-01 -2.78318984e+00
 -1.28650431e+00 -2.14988072e+00  1.33089167e+00  2.77576102e+00
  9.22228650e-01  2.93134122e+00  1.87176516e+00 -3.25563533e-01
 -1.34065979e+00 -1.78080863e+00 -5.04636399e-01  6.02804874e-01
 -7.43099604e-01 -5.21557531e-01  2.51661835e-01  1.32066880e+00
  9.54530217e-01  1.81284473e+00 -2.21049084e-01 -9.70741306e-01
 -1.49871020e-01  5.26171396e-01 -2.19475642e+00  1.06043978e+00
  1.31218541e+00  4.97280837e-03  2.54422862e+00  2.09262941e-01
 -1.72950655e+00 -9.26785314e-01 -9.08425078e-01  8.76610172e-01
  1.71[...]<jupyter_text>Plotting the approximationWe now plot the result of our optimisation in order to compare the QML algorithm's output and the analytical solution.<jupyter_code>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:])
    simulator_backend = pcvl.BackendFactory().get_backend("SLOS")
    px = pcvl.P("x")
    c = pcvl.Circuit(m, U_2) // (0, phys.PS(px)) // pcvl.Circuit(m, U_1)

    for x in X:
        px.set_value(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)
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[0] for v in loss_evolution])
plt.yscale("log")
plt.xlabel("Number of epochs")
plt.ylabel("Loss function value")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14111485222596334, 'fraction_numerical': 0.061803217358772915, 'mean_word_length': 4.451060358890701, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 31, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '32855', 'license': 'MIT', 'n_tokens_mistral': 4881, 'n_tokens_neox': 4357, 'n_words': 1293, 'path': 'docs/source/notebooks/Differential equation solving.ipynb', 'repo_name': 'Quandela/Perceval', 'size': '13365'}