Differential equation solving through quantum annealing
Project description
A differential equation solver using quantum annealing. It can be applied to coupled linear ODEs and PDEs with variable coefficients and inhomogeneous terms. The solution is obtained by finding the optimal weights for a linear combination of "basis functions":
The problem can then be viewed as a minimization one, for the loss function
where E is the equation to be solved, B are the initial/boundary conditions, and the x are the sample points at which they must hold. From this perspective, there is no difference between the equation and the boundary conditions, and they are treated on an equal footing. The loss function is then translated into an Ising model Hamiltonian by means of a binary encoding for the weights and solved in a quantum annealer.
Solutions obtained with qade for the wave equation and a system of coupled first order ODEs (see the examples
directory):
Contents
Installation
> pip install qade
This will install the dependencies numpy
and dwave-ocean-sdk
if not present. Python version 3.7 or higher is required. In order for qade
to send the problem for solution in the D-Wave systems, the API token needs to be configured.
Usage example
As an example, we can solve the Laguerre equation in the [0, 1]
interval with boundary conditions at the extremes:
where L_n
is the nth Laguerre polynomial. To do so, we parametrize the solution as a linear combination of the first few powers of x, and look for the optimal values of the weights such that the equation is satisfied at a set of sample points. This can be done through:
import numpy as np
from scipy.special import eval_laguerre as L
import qade
n = 4 # Laguerre equation parameter
x = np.linspace(0, 1, 20) # Sample points at which the equation is to be evaluated
y = qade.function(n_in=1, n_out=1) # Define the function to be solved for
# Define the equation and boundary conditions
eq = qade.equation(x * y[2] + (1 - x) * y[1] + n * y[0], x)
bcs = qade.equation(y[0] - [1, L(n, 1)], [0, 1])
# Solve the equation using as the basis of functions the monomials [1, x, x^2, x^3]
y_sol = qade.solve([eq, bcs], qade.basis("monomial", 4), scales=(n - 2), num_reads=500)
This illustrates the use of the 4 main functions provided by qade:
-
qade.function(n_in, n_out)
defines a symbolic functiony
ofn_in
inputs andn_out
outputs, whose nth derivative is denoted asy[n]
. -
qade.equation(eq, samples)
defines an equation of the formeq == 0
, withsamples
being an array-like object containing sample points at which it holds. -
qade.basis(name, size_per_dim, ...)
defines the basis of functions in terms of which the solution will be written. -
qade.solve(eqs, basis, ...)
uses a binary encoding to represent the problem of solving the given equations in the given basis as an Ising model, sends it to a D-Wave quantum annealer, and collects it into aSolution
objecty_sol
.
We can now print the weights of the basis functions in the solution and the corresponding through:
print(f"loss = {y_sol.loss:.3}, weights = {y_sol.weights}")
We can also evaluate the solution at any set of points, plot the results, and compare to the analyical solution, the which is nth Laguerre polynomial:
import matplotlib.pyplot as plt
plt.plot(x, y_sol(x), linewidth=5)
plt.plot(x, L(n, x), color="black", linestyle="dashed")
plt.show()
For other examples, see the examples
directory. An example PDE, the wave equation, is solved in wave.py
. A system of coupled first-order ODEs is solved in coupled.py
.
Documentation
qade.function
Creates Function
objects, to be used in definition of the equations. The arguments are:
-
n_in: int
. The number of input variables of the function, i.e. the dimension of the domain. -
n_out: int
. The number of output variables of the function, i.e. the dimension of the target space. Forn_out == 1
,qade.function
returns a singleFunction
. Forn_out > 1
, it returns list ofFunction
objects, one for each output component.
A differential equation or boundary condition is then defined as a linear combination of the derivatives of the functions, plus possible an inhomogeneous term, with coefficients and inhomogeneous term being either scalars or flat array-like objects (numpy arrays, lists, ...) having the same length as the set of sample points at which the equation is evaluated. The examples show how this is done in different practical settings. The notation f[k1, ..., kn]
represents the following expression:
qade.equation
Creates an Equation
, representing any linear differential equation or boundary condition. The arguments are:
-
formula
. The expression for the equation, constructed as a linear combination of derivativesf[k1, ..., kn]
ofFunction
objects, plus the optional inhomogeneous term. -
samples
, an array-like object of shape(n_samples, n_in)
(or just(n_samples,)
, ifn_in == 1
). Defines the sample points from the domain in which the equation must hold.
qade.basis
Creates a Basis
, representing the functions to be linearly combined into the solution. The arguments are:
-
name: str
. The allowed names are:-
"fourier"
. A product of the following functions for each componentx
in the input space: -
"monomial"
. A product of the following functions for each componentx
in the input space: -
"trig"
. A product of the following functions for each componentx
in the input space: -
"gaussian"
. The following functions of the distancer
to each pointx_n
in an equally-spaced grid in the input space (which is assumed to be the box[0, 1]^n_in
): -
"multiquadric"
. The following functions of the distancer
to each pointx_n
in an equally-spaced grid in the input space (which is assumed to be the box[0, 1]^n_in
):
-
-
size_per_dim: int
. The number of functions to use per input variable. The total dimension of the basis will besize_per_dim ** n_in
. -
n_in: int
. Only used by the radial basis functions"gaussian"
and"multiquadric"
to create the grid of points in the domain. -
scale: float
. The lambda scale parameter of the radial basis functions"gaussian"
and"multiquadric"
.
Custom bases
Custom bases can be provided by the user. In this case, instead of using qade.basis
, the user should define a class implementing the methods:
-
dimension(self, n_in: int) -> int
. Returning the total number of functions in the basis. -
derivatives(self, k: int, samples: np.ndarray) -> np.ndarray
. Returning the array of kth derivatives of the basis functions at the given sample points, received as an array of shape(n_samples, n_in)
. The output array should have shape:(n_samples, n_out) + (n_in,) * k
.
qade.solve
Solves the given equations and returns a Solution
object. The arguments are:
-
equations: List[Equation]
. The equations to be solved, including initial/boundary conditions. -
basis: Basis
. The basis to be used in their solution. -
n_spins: int = 3
. The number of spins per weight in the binary encoding. -
n_epochs: int = 10
. The number of epochs in which the problem is solved, using the values obtained in the previous epoch as the center values in the binary encoding of the weights. -
scale_factor: float = 0.5
. The relative change in the scale in the binary encoding of the weights from one epoch to the next. -
anneal_schedule: Optional[List[Tuple[int, int]]] = None
. The schedule to be used by the quantum annealer. See D-Wave Ocean Tools' documentation for details. IfNone
, it will default to a linear schedule of the form[(0, 0), (200, 0)]
. -
num_reads: int = 200
. The number of reads to be performed by the quantum annealer. -
qpu_solver: str = "Advantage_system4.1"
. The QPU to be used. -
centers: ArrayLike = 0
. The initial central values in the binary encoding of the weights. Must be a scalar or an array of sizen_out * basis_dim
. -
scales: ArrayLike = 1
. The initial scales in the binary encoding of the weights. Must be a scalar or an array of sizen_out * basis_dim
. -
verbose: bool = False
. Determines if information about intermediate steps should be printed. -
classical_minimizer: Optional[Callable[[Callable[[np.ndarray], float]], np.ndarray]] = None
. If notNone
, a classical function to be used in place of the quantum procedure for minimizing the loss function, for testing purposes.
Solution objects
Return type of qade.solve
. They are callables that receive an array-like object of sample points, with shape (n_samples, n_in)
(or just (n_samples,)
, if n_in == 1
), and return the corresponding values of the solution, as a numpy array of shape (n_samples, n_out)
. Their attributes are:
-
basis: Basis
. The basis in terms of which the solution has been obtained. -
weights: np.ndarray
, with shape(n_out, n_in)
. The optimal weights providing the solution as a linear combination of the basis functions. -
loss: float
. The value of the loss function at the given weights.
Inspecting the loss function
Two functions are provided to obtain the internal parameters of the loss function, both as a quadratic function of the real-valued weights, and as an Ising model Hamiltonian in the annealer spins.
qade.loss
Returns the (J, h)
parameters of the quadratic loss function Q(w) = w^T J w + h^T w
, where w
is the flattened array of the continuous weights. The arguments are
-
equations: List[Equation]
. The equations to be included. -
basis: Basis
. The basis to be used in their solution.
qade.ising
Returns the (J, h)
parameters of the Ising model Hamiltonian H(w) = w^T J w + h^T w
, where w
is the flattened array of spins. The arguments are:
-
equations: List[Equation]
. The equations to be included. -
basis: Basis
. The basis to be used in their solution. -
n_spins: int = 3
. The number of spins per weight in the binary encoding. -
centers: ArrayLike = 0
. The initial central values in the binary encoding of the weights. Must be a scalar or an array of sizen_out * basis_dim
. -
scales: ArrayLike = 1
. The initial scales in the binary encoding of the weights. Must be a scalar or an array of sizen_out * basis_dim
.
Citation
If you use qade, please cite:
Project details
Download files
Download the file for your platform. If you're not sure which to choose, learn more about installing packages.