PHYSICAL REVIEW A 101, 010301(R) (2020)
Rapid Communications
Variational quantum algorithms for nonlinear problems
Michael Lubasch ,
1
Jaewoo Joo,
1
Pierre Moinier,
2
Martin Kiffner ,
3,1
and Dieter Jaksch
1,3
1
Clarendon Laboratory, University of Oxford, Parks Road, Oxford OX1 3PU, United Kingdom
2
BAE Systems, Computational Engineering, Buckingham House, FPC 267, PO Box 5, Filton, Bristol BS34 7QW, United Kingdom
3
Centre for Quantum Technologies, National University of Singapore, 3 Science Drive 2, Singapore 117543
(Received 29 July 2019; revised manuscript received 11 December 2019; published 6 January 2020)
We show that nonlinear problems including nonlinear partial differential equations can be efficiently solved by
variational quantum computing. We achieve this by utilizing multiple copies of variational quantum states to treat
nonlinearities efficiently and by introducing tensor networks as a programming paradigm. The key concepts of
the algorithm are demonstrated for the nonlinear Schrödinger equation as a canonical example. We numerically
show that the variational quantum ansatz can be exponentially more efficient than matrix product states and
present experimental proof-of-principle results obtained on an IBM Q device.
DOI: 10.1103/PhysRevA.101.010301
Nonlinear problems are ubiquitous in all fields of science
and engineering and often appear in the form of nonlinear
partial differential equations (PDEs). Standard numerical ap-
proaches seek solutions to PDEs on discrete grids. However,
many problems of interest require extremely large grid sizes to
achieve accurate results, in particular in the presence of unsta-
ble or chaotic behavior that is typical for nonlinear problems
[1–3]. Examples include large-scale simulations for reliable
weather forecasts [4–6] and computational fluid dynamics
[7–9].
Quantum computers promise to solve problems that are
intractable on conventional, i.e., standard classical, computers
through their quantum-enhanced capabilities. In the context of
PDEs, it has been realized that quantum computers can solve
the Schrödinger equation faster than conventional computers
[10–12] and these ideas have been generalized recently to
other linear PDEs [13–18]. However, nonlinear problems are
intrinsically difficult to solve on a quantum computer due to
the linear nature of the underlying framework of quantum
mechanics.
Recently, the concept of variational quantum computing
(VQC) has attracted considerable interest [19–32] for solving
optimization problems. Variational quantum computing is a
quantum-classical hybrid approach where the evaluation of
the cost function C(λ) is delegated to a quantum computer,
while the optimization of variational parameters λ is per-
formed on a conventional classical computer. The concept of
VQC has been applied, e.g., to simulating the dynamics of
strongly correlated electrons through nonequilibrium dynam-
ical mean field theory [22,23,33,34], and quantum chemistry
calculations were successfully carried out on existing noisy
superconducting [21,25,26] and ion quantum computers [28].
We extend and adapt the concept of VQC to solving
nonlinear problems efficiently on a quantum computer by
virtue of two key concepts. First, we introduce a quan-
tum nonlinear processing unit (QNPU) that efficiently calcu-
lates nonlinear functions of the form F = f
(1)
∗
r
j=1
(O
j
f
( j)
)
for VQC. Measuring the ancilla qubit connected to the
QNPU as shown in Fig. 1(a) directly yields the sum of all
function values
k
Re{F
k
}, where Re{·} denotes the real part.
The functions f
( j)
are encoded in variational n-qubit states
|ψ (λ)
j
=
ˆ
U
j
(λ)|0 created by networks of the form shown in
Fig. 1(b). The same function f
(i)
= f
( j)
may appear multiple
times by choosing
ˆ
U
i
(λ) =
ˆ
U
j
(λ). Second, we use tensor
networks as a programming paradigm for QNPUs to create
optimized circuits that efficiently calculate linear operators O
j
acting on functions f
( j)
. In this way all quantum resources
of nonlinear VQC scale polynomially with the number of
qubits, which represents an exponential reduction compared
with some conventional algorithms.
The variational states |ψ (λ )
j
represent N = 2
n
values of
the functions f
( j)
which form a trial solution to the problem
of interest. The cost function C(λ) for nonlinear VQC is built
up from outputs of different QNPUs that are then processed
classically to iteratively determine the optimal set λ. Large
grid sizes that are intractable on a conventional computer
require only n 20 qubits, which is within the reach of noisy
intermediate-scale quantum (NISQ) devices. In addition, the
scheme is applicable to other types of nonlinear problems that
can be solved via the minimization of a cost function C(λ)
[36].
We demonstrate the concept and performance of nonlinear
VQC by emulating it classically for the canonical example of
the time-independent one-dimensional nonlinear Schrödinger
equation
−
1
2
d
2
dx
2
+V (x) + g|f (x)|
2
f (x) = Ef( x), (1)
where V is an external potential and g denotes the s trength
of the nonlinearity. We also implement nonlinear VQC for
Eq. (1) on IBM quantum computers to establish its feasibility
on current NISQ devices. Proof-of-principle results are shown
in Figs. 1(c) and 1(d), demonstrating excellent agreement
with numerically exact solutions. The nonlinear Schrödinger
equation and its generalization to higher dimensions describes
various physical phenomena ranging from Bose-Einstein con-
densation to light propagation in nonlinear media [37–42].
In particular, below we consider Eq. (1) with quasiperiodic
2469-9926/2020/101(1)/010301(7) 010301-1 ©2020 American Physical Society
MICHAEL LUBASCH et al. PHYSICAL REVIEW A 101, 010301(R) (2020)
036
300
6000
0
0
2
4
FIG. 1. (a) Quantum network for summing a nonlinear function
of the form F = f
(1)
∗
r
j=1
(O
j
f
( j )
). The ancilla qubit on the top
line undergoes Hadamard gates
ˆ
H and controls operations of the
rest of the network via the control port (CP) and is measured
in the computational basis. Starting from product states |0 of n
qubits shown as thick lines, variational states |ψ(λ)
j
=
ˆ
U
j
(λ)|0
representing the functions f
( j )
are created and fed to the QNPU
through input ports (IPs). The QNPU contains problem-specific
quantum networks defining the linear operators O
j
and has the output
ports (OPs). (b) Network
ˆ
U (λ)ofdepthd = 5 with n = 6. The
values λ ={λ
1
, λ
2
,...} determine the form of the two-qubit gates.
(c) Cost function C(λ) for the nonlinear Schrödinger equation (1)for
a harmonic potential V , a single variational parameter λ,andg = 10
4
(solid line) and g = 10 (dashed line). The arrows indicate optimal
values λ = λ
min
. (d) Solutions |f (x)|
2
for λ = λ
min
with N = 4grid
points and periodic boundary conditions f (b) = f (a)forg = 10
4
(solid line) and g = 10 (dashed line). The circles and squares are
numerically exact results. Experimental data in (c) and (d) were
obtained, from a modified circuit, on IBM Q devices (see [35]for
details).
potentials V that are in the focus of current cold-atom ex-
periments [43] and that make Eq. (1) challenging to solve
numerically. The methods used for this equation here are
straightforwardly modified to handle other nonlinear terms
and time-dependent problems as illustrated in [35]forthe
Burgers equation appearing in fluid dynamics.
The ground state of Eq. (1) can be found by minimizing the
cost function
C =K
c
+P
c
+I
c
, (2)
where K
c
, P
c
, and I
c
are the mean kinetic, potential,
and interaction energies, respectively. In Eq. (2)the ·
c
denote averages with respect to a single real-valued function
f
(1)
≡ f on the interval [a, b] satisfying the normalization
condition
b
a
|f ( x)|
2
dx = 1.
In line with standard numerical approaches [44–46], we
apply the finite-difference method (FDM) to Eq. (1) and
discretize the interval [a, b]intoN equidistant grid points
x
k
= a + h
N
k, where h
N
= /N is the grid spacing, = b − a
is the length of the interval, and k ∈{0,...,N − 1}. Each grid
point is associated with a variational parameter f
k
that approx-
imates the continuous solution f (x
k
)atx
k
. Furthermore, we
impose periodic boundary conditions, i.e., f
N
= f
0
, and the
normalization condition imposed on the continuous functions
f translates to
1 = h
N
N−1
k=0
|f
k
|
2
=
N−1
k=0
|ψ
k
|
2
, (3)
where ψ
k
=
√
h
N
f
k
. Note that the condition on the set of
parameters {ψ
k
} is independent of the grid spacing, and in
the following we consider optimizing the cost function with
respect to them.
All averages ·
c
in Eq. (2) can be approximated by
corresponding expressions of the discrete problem · .We
find [44–46] ·
c
=·+E
grid
, where E
grid
∝ 1/N
2
is the
error associated with the trapezoidal rule when transforming
integrals into sums and
K = −
1
2
1
h
2
N
N−1
k=0
ψ
∗
k
(ψ
k+1
− 2ψ
k
+ ψ
k−1
), (4a)
P =
N−1
k=0
[ψ
∗
k
V (x
k
)ψ
k
], (4b)
I =
1
2
g
h
N
N−1
k=0
|ψ
k
|
4
. (4c)
Note that K in Eq. (4a) uses an FDM representation of
the second-order derivative in Eq. (1).
To evaluate the terms in Eq. (4) on a quantum com-
puter we consider quantum registers with n qubits and basis
states |q=|q
1
···q
n
=|q
1
⊗···⊗|q
n
, where q
j
∈{0, 1}
denotes the computational states of qubit j. Regarding the
sequence q
1
...q
n
= binary(k) as the binary representation of
the integer k =
n
j=1
q
j
2
n−j
, we encode all N = 2
n
ampli-
tudes ψ
k
in the normalized state
|ψ=
N−1
k=0
ψ
k
|binary(k). (5)
We prepare the quantum register in a variational state |ψ(λ)
via the quantum circuit
ˆ
U (λ) of depth d shown in Fig. 1(b).
We consider depths d ∝ poly(n) such that the quantum ansatz
requires exponentially fewer parameters N than standard
classical schemes with N parameters. Note that the number
of variational parameters scales like N ∝ nd for our quantum
ansatz [35]. The power of this ansatz is rooted in the fact
that it encompasses all matrix product states (MPSs) [47–50]
with bond dimension χ ∼ poly(n)[35]. Since polynomials
and Fourier series [51,52] can be efficiently represented by
MPSs, the quantum ansatz simultaneously contains universal
basis functions that are capable of approximating a large class
of solutions to nonlinear problems efficiently. Furthermore,
we show below that the quantum ansatz is capable of
storing solutions with exponentially fewer resources than the
classically optimized MPS ansatz. Note that the number of
variational parameters scales like N ∝ nχ
2
for MPSs of bond
dimension χ [35].
Figure 2(a) demonstrates the basic working principle of
the QNPU for the nonlinear term I . The effect of the
010301-2
VARIATIONAL QUANTUM ALGORITHMS FOR NONLINEAR … PHYSICAL REVIEW A 101, 010301(R) (2020)
FIG. 2. (a) QNPU circuit calculating the nonlinear term |ψ|
4
.
The networks are shown for n = 3 and all input ports are fed the
same variational quantum states created by
ˆ
U (λ). IP2 is fed
ˆ
U (λ)
and IP3 is fed
ˆ
U
∗
(λ). (b) QNPU circuit for working out the potential
energy term
˜
V |ψ|
2
.
controlled-NOT operations between pairs of qubits is to pro-
vide a pointwise multiplication with the ancilla thus measur-
ing
k
|ψ
k
|
4
.InFig.2(b) we show the circuit for measuring
P . The unitary
ˆ
V ˆ=O
1
encodes function values of the
external potential V . A copy of ψ is effectively multiplied
pointwise with the external potential by controlled-
NOT gates
to give
k
˜
V
k
|ψ
k
|
2
. Similarly, multiplying ψ with their shifted
versions using adder circuits (see [35] for details) allows
evaluating the kinetic energy term.
The measured expectation value of the ancilla qubit is di-
rectly related to the desired quantities as I = gˆσ
z
I
anc
/2h
N
for the nonlinear term, P = αˆσ
z
P
anc
for the potential
energy, and K = (1 −ˆσ
z
K
anc
)/h
2
N
for the kinetic energy
[53,54]. Furthermore, derivatives of the cost function, as
required by some minimization algorithms [36,55,56], can
be evaluated by combining the ideas presented here with the
quantum circuits discussed in [24,30,57,58].
The unitary network
ˆ
V represents scaled function values
˜
V
k
of the external potential where
N−1
k=0
|
˜
V
k
|
2
= 1, and α>0
is a scaling parameter such that V
k
= α
˜
V
k
. Efficient quantum
circuits
ˆ
V for measuring P can be systematically obtained
by establishing tensor networks as a programming paradigm.
To this end we expand the external potential in polynomials
or Fourier series
˜
V (x) ≈
J
j
c
j
b
j
(x), where b
j
(x) are basis
functions and c
j
are expansion coefficients [35]. In the case
of Fourier series of order J, the approximate potential is rep-
resented by an MPS of bond dimension χ = J [51,52]. Next
we write the MPSs in terms of n −log χ unitaries [59–61],
where · is the ceiling function. Each of these unitaries acts
on 2χ qubits and can be decomposed in terms of elementary
two-qubit gates [35,62–65]. An upper bound for the depth of
the resulting quantum circuit is d
>
9n[
23
48
(2χ )
2
+
4
3
][64].
The depth thus scales polynomially with the number of qubits
n and with χ , and many problems of interest show an even
more advantageous scaling. For example, in the following we
consider the potential
V (x) = s
1
sin(κ
1
x) + s
2
sin(κ
2
x)(6)
and set κ
2
= 2κ
1
/(1 +
√
5). This potential realizes an incom-
mensurate bichromatic lattice where the ratio s
1
/s
2
determines
the amount of disorder in the lattice [66]. The trap potential
V (x)inEq.(6) is exactly represented by an MPS of bond
dimension χ = 4. The depth of the corresponding quantum
circuit d = 5(n − 2) + 1 d
>
is much smaller than the up-
per bound [35].
Next we analyze the Monte Carlo sampling error [44]
associated with the measurement of the ancilla qubit. We
denote the absolute sampling error associated with quantity
X by E
X
MC
and the corresponding relative error is [35]
P
MC
=
E
P
MC
P
= C
P
1
√
M
, (7a)
K
MC
=
E
K
MC
K
≈ C
K
N
N
min
1
√
M
, (7b)
I
MC
=
E
I
MC
I
≈ C
I
N
N
min
1
√
M
. (7c)
In these equations, we assume N N
min
and N
min
=
/
min
is the minimal number of grid points for resolving the
smallest length scale
min
of the problem. The parameters C
X
in Eqs. (7) are of the order of unity [35] and all sampling
errors decrease with the number of samples M as 1/
√
M.
While the relative error associated with the potential term in
Eq. (7a) is independent of the number of grid points,
K
MC
and
I
MC
increase linearly with N/N
min
. It follows that increasing
the grid size requires larger values of M in order to keep
the sampling error small. However, the grid error scales like
E
grid
∝ 1/N
2
= 2
−2n
for N N
min
[35]. We thus conclude
that only moderate ratios N/N
min
> 1 and therefore relatively
small values of M are needed in order to achieve accurate
solutions with small grid errors.
The quantum ansatz in Fig. 1(b) is inspired from ten-
sor network theory and can be regarded, for example, as
the Trotter decomposition of the time t evolution operator
exp[−iH(t )t/¯h] of a time-dependent spin Hamiltonian H (t )
with arbitrary short-range interactions acting on the initial
state |0 [67]. Similarly to the coupled cluster ansatz in
quantum chemistry VQC calculations [19], there is currently
no known efficient classical ansatz for this state [68,69]. From
the VQC perspective, our quantum ansatz in Fig. 1(b) is
composed of generic two-qubit gates that in an experiment
are accurately approximated by short sequences of gates if
a sufficiently tunable or universal gate set is experimentally
available [70]. We envisage that this quantum ansatz is more
efficient than methods based on an MPS ansatz on a classical
computer like the multigrid renormalization (MGR) method
in [56]. This MGR method is the most efficient and accurate
classical algorithm known to us for the problem considered
here, for which it can already be exponentially faster than
standard classical algorithms. A comparison with this pow-
erful classical method, which is based on variational classical
MPSs, will allow us to validate the superior variational power
of the quantum ansatz in Fig. 1(b) on a quantum computer.
Note that the difficulty of solving a classical optimization
problem with many variables does not go away by using the
quantum ansatz, as the actual optimization is classical and
there is no quantum advantage there. The quantum advantage
stems solely from the faster evaluation of the cost function for
our quantum ansatz in Fig. 1(b) on a quantum computer.
010301-3
MICHAEL LUBASCH et al. PHYSICAL REVIEW A 101, 010301(R) (2020)
10
-2
10
-1
10
0
0.02 0.04 0.06
32 64 128
0
500
0.01 0.03 0.05
32 64 128
0
2
10
-2
10
0
0.0 0.5 1.0
10
-12
10
-6
10
0
FIG. 3. (a) Numerically exact solution |f (x)|
2
of Eq. (1) on a log-
arithmic scale, V (x)inEq.(6) with s
1
= 2 × 10
4
and κ
1
= 2π × 32.
The green thin line (blue thick line) is for s
1
/s
2
= 200 (s
1
/s
2
= 2).
(b) A log-log plot of the IPR (top panel) and lin-log plot of S
max
of the exact solution |ψ
exact
(bottom panel) as a function of κ
1
.
Green crosses (blue circles) correspond to s
1
/s
2
= 200 (s
1
/s
2
= 2).
(c) Representation error
R
of the exact solution for the quantum
ansatz (thick lines) and the MPS ansatz (thin lines) as a function of
N /N. All curves are for s
1
/s
2
= 2 and correspond to s
1
= 2 × 10
4
and κ
1
= 2π × 32 (green dash-dotted line), s
1
= 8 × 10
4
and κ
1
=
2π × 64 (blue dashed line), and s
1
= 3.2 × 10
5
and κ
1
= 2π × 128
(purple solid line). The inset shows N as a function of κ
1
for
R
= 0.05 and s
1
/s
2
= 2. Blue squares (green triangles) correspond
to the quantum ansatz (MPS ansatz). All curves in (a)–(c) are for
N = 2
13
= 8192 grid points and g = 50.
To provide numerical evidence for this we first obtain the
numerically exact solution of Eq. (1) on the interval [0, 1] via
the MGR algorithm [56] and by allowing for the maximal
bond dimension χ of the MPS ansatz [ 35]. In this case the
numerically exact solution is described by N = 2
n
parameters
like in other conventional algorithms. The results are shown
in Fig. 3(a) for two different values of s
1
/s
2
. In the weakly
disordered regime s
1
/s
2
1, |f (x)|
2
varies on the length
scale set by 1/κ
1
. On the contrary, the strongly disordered
regime s
1
/s
2
≈ 1 is characterized by strongly localized so-
lutions in space. The localization of the wave function can
be quantified using the inverse participation r atio (IPR) [71]
IPR = (N
N−1
k=0
|ψ
k
|
4
)
−1
. We show the IPR in Fig. 3(b) as a
function of κ
1
(top panel) and find that it stays constant for
s
1
/s
2
1. On the other hand, the IPR decreases according
to a power law with κ
1
for s
1
/s
2
≈ 1, showing that the lo-
calized character of the wave function increases dramatically
with κ
1
.
Next we encode the function values of the numerically
exact solution in the state |ψ
exact
via Eq. (5) and calculate the
maximum bipartite entanglement entropy S
max
of all possible
bipartitions of the n-qubit wave function |ψ
exact
. The quantity
S
max
is a measure of the entanglement of |ψ
exact
and is shown
in the bottom panel of Fig. 3(b).ThevalueofS
max
is small
and stays constant with κ
1
for s
1
/s
2
1 i n the weakly disor-
dered regime. Contrary to this, for s
1
/s
2
≈ 1 in the strongly
disordered regime we observe that S
max
∝ log(κ
1
).
The entanglement measure S
max
provides a useful nec-
essary criterion for efficient MPS approximations [72]. A
MPS of bond dimension χ can at most contain an amount
of entanglement S
max
[MPS] log
2
(χ )[47,50]. For MPSs to
be efficient, we require that χ scales at most polynomially
with n, i.e., χ = poly(n), so that S
max
[MPS] log[poly(n)]
and therefore MPS can only capture small amounts of entan-
glement efficiently. The small values of S
max
for s
1
/s
2
1
suggest that MPSs work well in the weakly disordered regime
and indeed we have confirmed numerically that this is true.
Therefore, in the following we focus on the strongly disor-
dered regime s
1
/s
2
≈ 1. In this regime MPS cannot be an ef-
ficient approximation for large values of κ
1
: The total number
of variational parameters of MPSs N [MPS] depends quadrat-
ically on χ [35], i.e., N [MPS] ∝ χ
2
, and χ needs to grow
polynomially with κ
1
, i.e., χ ∝ poly(κ
1
), to satisfy the ob-
served entanglement scaling, such that N [MPS] ∝ poly(κ
1
).
Our quantum ansatz (QA) of Fig. 1(b) can capture much
larger amounts of entanglement S
max
[QA] d efficiently [69]
and therefore this ansatz can be an efficient approximation
for large values of κ
1
: Because the total number of varia-
tional parameters N depends linearly on d for the quantum
ansatz [35], i.e., N [QA] ∝ d, and d just needs to grow
logarithmically, i.e., d ∝ log(κ
1
), for the observed entan-
glement requirements, we conclude that N [QA] ∝ log(κ
1
).
These entanglement considerations show that the quantum
ansatz has the potential to be exponentially more efficient than
MPSs in the strongly disordered regime for increasing values
of κ
1
.
To quantitatively analyze and demonstrate the efficiency
of the quantum ansatz in this regime, we obtain the set of
parameters λ that maximize the fidelity F = |ψ
exact
|ψ (λ)|
for different depths d [35]. The infidelity
R
= 1 − F is thus
a measure of the error when approximating the exact solution
by this ansatz, and in the following we refer to
R
as the
representation error. As shown in Fig. 3(c), the representation
error decreases exponentially as a function of N for all
values of κ
1
and therefore we obtain accurate solutions for
N /N 1. Even for the largest value of κ
1
= 2π × 128 and
R
≈ 10
−2
, we find N /N ≈ 0.04, so we require only 4% of
the full number of parameters needed in conventional algo-
rithms. Most importantly, the inset of Fig. 3(c) shows that, to
obtain a fixed representation error of
R
= 0.05, the number of
parameters of our quantum ansatz needs to grow as N [QA] ∝
log (κ
1
). This numerical analysis therefore confirms our ex-
pectation, from the entanglement arguments in the preceding
paragraph, that the quantum ansatz efficiently approximates
solutions in the strongly disordered regime even for large
values of κ
1
. This is possible because the quantum ansatz
captures the required entanglement S
max
∝ log(κ
1
) by means
of just the small number of parameters N [QA] ∝ log(κ
1
).
The entanglement capabilities of MPSs imply that N [MPS] ∝
poly(κ
1
) has to be fulfilled. Therefore, we conclude that our
quantum ansatz is exponentially more efficient than MPSs in
the strongly disordered regime for growing values of κ
1
.This
010301-4
VARIATIONAL QUANTUM ALGORITHMS FOR NONLINEAR … PHYSICAL REVIEW A 101, 010301(R) (2020)
key finding shows that nonlinear VQC can be exponentially
more efficient than optimized classical variational schemes
that are based on the MPS ansatz.
We now propose a calculation that illustrates the
exponential advantage of our quantum ansatz on a quantum
computer particularly clearly. We consider the ground-state
problem of Eq. (1) with the quasirandom external potential of
Eq. (6) on the interval [0, 1] discretized by N = 2
n
equidistant
grid points and use the FDM from above that leads to
Eqs. (4a)–(4c). We choose κ
1
= 2π × 2
n/2
and assume that
n = 4, 6, 8,... qubits store the variational quantum ansatz.
The smallest wavelength present in our external potential
of Eq. (6) is determined by κ
1
and reads 2π/κ
1
= 2
−n/2
.
The interval [0, 1] accommodates 1/2
−n/2
= 2
n/2
of such
wavelengths and our grid of N = 2
n
equidistant points re-
solves each such wavelength using N/2
n/2
= 2
n/2
grid points.
Therefore, the randomness of our quasirandom external
potential of Eq. (6) as well as its FDM resolution grows
exponentially with the number of qubits n.Atthesametime,
all FDM errors (using an FDM representation of the Laplace
operator, the trapezoidal rule for integration, and so on)
decrease polynomially with the grid spacing 1/N = 2
−n
and
thus exponentially with n. Therefore, with growing values of
n, our FDM representation of Eq. (1) converges exponentially
fast to the continuous problem and the accuracy of random-
ness in Eq. (6) grows exponentially. Our proposal is to solve
this problem for the strongly disordered regime s
1
/s
2
≈ 1in
Eq. (6), i.e., compute the corresponding ground states, using
an increasing number n = 4 , 6, 8,.... Based on the numerical
analysis of Fig. 3, we conjecture that MPSs require resources
N [MPS] ∝ poly(κ
1
) = poly(2π × 2
n/2
) = exp(n)growing
exponentially with n, whereas the quantum ansatz of Fig. 1(b)
just needs resources N [QA] ∝ log(κ
1
) = log(2π × 2
n/2
) =
poly(n) increasing polynomially with n. This calculation
can be used to demonstrate quantum supremacy as, by
successively increasing n, the limit of what is computationally
possible on a classical computer (using MPSs or exact classi-
cal methods) is reached quickly and our quantum ansatz on a
quantum computer remains efficient far beyond this classical
limit.
The quantum hardware requirements for this quantum
supremacy calculation go beyond the current capabilities
of available NISQ devices [73]. Nevertheless, the superior
performance of our quantum ansatz is relevant for current
NISQ devices as only the most efficient variational states can
succeed in the presence of the current experimental errors.
We test the feasibility of nonlinear VQC on NISQ devices by
calculating the ground state of Eq. (1) for a simple harmonic
potential and a single variational parameter on an IBM Q
device utilizing further network optimizations (see [35]for
details). The experimental implementation of the nonlinear
VQC algorithm is able to identify the optimal variational
parameter with an error of less than 10%, leading to excellent
agreement of the ground-state solutions with exact numerical
solutions [cf. Figs. 1(c) and 1(d)].
The methods presented here are readily modified to two-
and three-dimensional problems, with an overhead scaling
linearly in the number of dimensions, and can be applied to
a broad range of nonlinear terms and differential operators.
An exciting prospect for future work would be to utilize
intermediate-scale quantum computers to solve nonlinear
problems on grid sizes beyond the scope of conventional
computers.
M.L. and D.J. are grateful for funding from the Networked
Quantum Information Technologies Hub of the UK National
Quantum Technology Programme as well as from the EP-
SRC grant “Tensor network theory for strongly correlated
quantum systems” (No. EP/K038311/1). We acknowledge
support from the EPSRC National Quantum Technology Hub
in Networked Quantum Information Technology (Grant No.
EP/M013243/1). M.K. and D.J. acknowledge financial sup-
port from the National Research Foundation and the Ministry
of Education, Singapore.
[1] S. Wiggins, Introduction to Applied Nonlinear Dynamical Sys-
tems and Chaos (Springer, New York, 2003).
[2] L. A. Aguirre and C. Letellier, Modeling nonlinear dynamics
and chaos: A review, Math. Probl. Eng. 2009, 238960 (2009).
[3] S. H. Strogatz, Nonlinear Dynamics and Chaos (Westview,
Boulder, 2015).
[4] P. Lynch, The Emergence of Numerical Weather Prediction
(Cambridge University Press, Cambridge, 2006).
[5] T.T.Warner,Numerical Weather and Climate Prediction (Cam-
bridge University Press, Cambridge, 2011).
[6] R. A. Pielke, Sr., Mesoscale Meteorological Modeling (Aca-
demic, Waltham, 2013).
[7] L. M. Milne-Thomson, Theoretical Aerodynamics (Dover, New
York, 1973).
[8] J. H. Ferziger and M. Peri
´
c, Computational Methods for Fluid
Dynamics (Springer, Berlin, 2002).
[9] H. K. Versteeg and W. Malalasekera, An Introduction to Com-
putational Fluid Dynamics (Pearson Education, Harlow, 2007).
[10] S. Lloyd, Universal quantum simulators, Science 273, 1073
(1996).
[11] A. Aspuru-Guzik, A. D. Dutoi, P. J. Love, and M. Head-Gordon,
Simulated quantum computation of molecular energies, Science
309, 1704 (2005).
[12] I. Kassal, S. P. Jordan, P. J. Love, M. Mohseni, and A. Aspuru-
Guzik, Polynomial-time quantum algorithm for the simulation
of chemical dynamics, Proc. Natl. Acad. Sci. USA 105, 18681
(2008).
[13] A. W. Harrow, A. Hassidim, and S. Lloyd, Quantum Algorithm
for Linear Systems of Equations, Phys. Rev. Lett. 103, 150502
(2009).
[14] D. W. Berry, High-order quantum algorithm for solving linear
differential equations, J. Phys. A: Math. Theor. 47, 105301
(2014).
[15] D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D.
Somma, Simulating Hamiltonian Dynamics with a Truncated
Taylor Series, Phys.Rev.Lett.114, 090502 (2015).
010301-5
![FIG. 1. (a) Quantum network for summing a nonlinear function of the form F = f (1)∗ ∏rj=1(Oj f ( j) ). The ancilla qubit on the top line undergoes Hadamard gates Ĥ and controls operations of the rest of the network via the control port (CP) and is measured in the computational basis. Starting from product states |0〉 of n qubits shown as thick lines, variational states |ψ (λ) j〉 = Ûj (λ)|0〉 representing the functions f ( j) are created and fed to the QNPU through input ports (IPs). The QNPU contains problem-specific quantum networks defining the linear operators Oj and has the output ports (OPs). (b) Network Û (λ) of depth d = 5 with n = 6. The values λ = {λ1,λ2, . . .} determine the form of the two-qubit gates. (c) Cost function C(λ) for the nonlinear Schrödinger equation (1) for a harmonic potential V , a single variational parameter λ, and g = 104 (solid line) and g = 10 (dashed line). The arrows indicate optimal values λ = λmin. (d) Solutions | f (x)|2 for λ = λmin with N = 4 grid points and periodic boundary conditions f (b) = f (a) for g = 104 (solid line) and g = 10 (dashed line). The circles and squares are numerically exact results. Experimental data in (c) and (d) were obtained, from a modified circuit, on IBM Q devices (see [35] for details).](/figures/fig-1-a-quantum-network-for-summing-a-nonlinear-function-of-17nrz4b8.webp)
