# Solving Fractional Differential Equations on a Quantum Computer:
A Variational Approach

Fong Yew Leong  Institute of High Performance Computing (IHPC), Agency for
Science, Technology and Research (A*STAR), 1 Fusionopolis Way, #16-16
Connexis, Singapore 138632, Singapore 0000-0002-0064-0118 Dax Enshan Koh
Institute of High Performance Computing (IHPC), Agency for Science, Technology
and Research (A*STAR), 1 Fusionopolis Way, #16-16 Connexis, Singapore 138632,
Singapore Science, Mathematics and Technology Cluster, Singapore University of
Technology and Design, 8 Somapah Road, Singapore 487372, Singapore
0000-0002-8968-591X Jian Feng Kong  Institute of High Performance Computing
(IHPC), Agency for Science, Technology and Research (A*STAR), 1 Fusionopolis
Way, #16-16 Connexis, Singapore 138632, Singapore 0000-0001-5980-4140 Siong
Thye Goh  Institute of High Performance Computing (IHPC), Agency for Science,
Technology and Research (A*STAR), 1 Fusionopolis Way, #16-16 Connexis,
Singapore 138632, Singapore Singapore Management University, 81 Victoria
Street, 188065, Singapore  0000-0001-7563-0961 Jun Yong Khoo  Institute of
High Performance Computing (IHPC), Agency for Science, Technology and Research
(A*STAR), 1 Fusionopolis Way, #16-16 Connexis, Singapore 138632, Singapore
0000-0003-0908-3343 Wei-Bin Ewe  Institute of High Performance Computing
(IHPC), Agency for Science, Technology and Research (A*STAR), 1 Fusionopolis
Way, #16-16 Connexis, Singapore 138632, Singapore 0000-0002-4600-0634 Hongying
Li  Institute of High Performance Computing (IHPC), Agency for Science,
Technology and Research (A*STAR), 1 Fusionopolis Way, #16-16 Connexis,
Singapore 138632, Singapore 0000-0003-4736-295X Jayne Thompson  Institute of
High Performance Computing (IHPC), Agency for Science, Technology and Research
(A*STAR), 1 Fusionopolis Way, #16-16 Connexis, Singapore 138632, Singapore
Dario Poletti  Science, Mathematics and Technology Cluster and Engineering
Product Development Pillar, Singapore University of Technology and Design, 8
Somapah Road, 487372 Singapore  Centre for Quantum Technologies, National
University of Singapore 117543, Singapore  MajuLab, CNRS-UNS-NUS-NTU
International Joint Research Unit, UMI 3654, Singapore

###### Abstract

We introduce an efficient variational hybrid quantum-classical algorithm
designed for solving Caputo time-fractional partial differential equations.
Our method employs an iterable cost function incorporating a linear
combination of overlap history states. The proposed algorithm is not only
efficient in time complexity, but has lower memory costs compared to classical
methods. Our results indicate that solution fidelity is insensitive to the
fractional index and that gradient evaluation cost scales economically with
the number of time steps. As a proof of concept, we apply our algorithm to
solve a range of fractional partial differential equations commonly
encountered in engineering applications, such as the sub-diffusion equation,
the non-linear Burgers’ equation and a coupled diffusive epidemic model. We
assess quantum hardware performance under realistic noise conditions, further
validating the practical utility of our algorithm.

##  1 Introduction

Differential equations play an important role in the modeling of many
practical engineering problems. Fractional differential equations, which are
generalizations of differential equations to arbitrary non-integer order,
involve fractional derivatives such as $\frac{d^{\alpha}}{dx^{\alpha}}$, where
$\alpha$ is a non-integer real number. These models have been shown to be
capable of modeling more complex processes that exhibit memory effects or non-
local behavior. For example, fractional differential equations have found
applications in fractional advection-dispersion equations [1] and viscoelastic
flow problems [2]. Fractional-order derivatives in time have also been used to
model cancer tumor growth, by numerically fitting the fractional order [3].

Here, we consider numerical solutions to time-fractional diffusive
differential equations of order $\alpha$, which derive from non-Markovian
continuous-time random walks [4]. Finite difference methods are effective and
simple to implement for time-fractional problems [5], where implicit methods
are preferred to explicit methods for numerical stability [6]. However, the
time-fractional derivative creates a _global dependence_ problem [7], where
previous solutions are accessed from memory, rendering storage expensive for
spatially large problems.

Quantum computers can offer new tools for efficiently solving differential
equations. In particular, variational quantum algorithms (VQAs) have emerged
as a leading strategy in the noisy intermediate-scale quantum (NISQ) era [8,
9], featuring low-depth parameterized quantum circuits, alongside classical
optimizers used to minimize a cost function. These VQAs yield quantum states
that encode approximate solutions to partial differential equations. They have
been successfully applied to both linear [10, 11, 12] and nonlinear
differential equations [13, 14] in a diverse range of fields, including
electromagnetics [15], collodial transport [16], fluid dynamics [17, 18, 19],
and finance [20].

In this paper, we propose and implement a variational quantum optimization
scheme for solving time-fractional diffusive partial differential equations in
space and time, either implicitly or semi-implicitly, using simple quantum
overlap measurements for the time integral term. We show that our approach is
not only efficient in time complexity, but can also yield memory cost
advantages in handling the global dependence problem. The rest of the paper is
organized as follows. In Section 2, we introduce the time-fractional diffusion
equation and present our variational quantum optimization scheme for numerical
integration. In Section 3, we discuss the details of the implementation and
evaluate time and memory complexities of our scheme. In Section 4, we present
numerical experiments, including time-fractional Burgers’ equation for
hydrodynamics (Section 4.2), and fractional diffusive epidemic model (Section
4.3). Lastly, we trial our algorithm on noise model simulators and hardware in
Section 5, and conclude with discussions in Section 6.

##  2 Fractional diffusion equations

We consider the fractional diffusion equation [21, 6] in space $x$ and time
$t$,

$\displaystyle D_{t}^{\alpha}u(t,x)=\partial_{xx}u(t,x),$ (1)

where $D_{t}^{\alpha}$ is the Caputo fractional derivative of the $\alpha$-th
order with respect to $t$, defined as

$\displaystyle\frac{D^{\alpha}u(t,x)}{\partial{t^{\alpha}}}=\frac{1}{\Gamma(n-\alpha)}\int^{t}_{0}\frac{\partial^{n}{u(s,x)}}{\partial
s^{n}}\frac{ds}{(t-s)^{\alpha-n+1}},\quad n-1<\alpha<n,$ (2)

where $\Gamma$ denotes the Gamma function and
$n=\lceil\alpha\rceil\in\mathbb{N}$ is the ceiling of $\alpha$.

###  2.1 Difference scheme

We set up a space-time rectangular domain with spatial extent $L$ and temporal
duration $T$, dividing it into a regularly spaced grid with integer size
$N\times M$. The spatial grid points are defined as $x_{i}=ih$ where $h=L/N$
and $i=1,2,\dots,N$, and the temporal grid points as $t_{k}=k\tau$ where
$\tau=T/M$ and $k=1,2,\dots,M$.

The first-order finite difference of the Caputo derivative of order
$(0<\alpha\leq 1)$ is approximated by [5, 22]

$\displaystyle D_{t}^{\alpha}u(t_{k},x_{n})\simeq
g_{\alpha,\tau}\sum_{j=1}^{k}w_{j}^{(\alpha)}(u_{n}^{k-j+1}-u_{n}^{k-j}),$ (3)

where $u_{n}^{k}\equiv u(t_{k},x_{n})$,
$g_{\alpha,\tau}:=\tau^{-\alpha}/\Gamma(2-\alpha)$ and
$w_{j}^{(\alpha)}:=j^{1-\alpha}-(j-1)^{1-\alpha}$. Rearranging Eq. 3 yields a
discrete sum of terms $u_{n}^{j}$ as

$\displaystyle D_{t}^{\alpha}u(t_{k},x_{n})\simeq
g_{\alpha,\tau}\left[u_{n}^{k}-w_{k}u_{n}^{0}+\sum_{j=1}^{k-1}\left(w_{j+1}-w_{j}\right)u_{n}^{k-j}\right],$
(4)

where $w_{1}=1$ by definition.

The finite difference of the spatial derivative is approximated by

$\displaystyle\partial_{xx}u(t_{k},x_{n})\simeq\frac{u_{n+1}^{k}-2u_{n}^{k}+u_{n-1}^{k}}{h^{2}}.$
(5)

Together, the difference scheme for the fractional diffusion equation can be
iterated in time as

$\displaystyle\left(1-a\delta\right)u_{n}^{k}=w_{k}u_{n}^{0}-\sum_{j=1}^{k-1}\left(w_{j+1}-w_{j}\right)u_{n}^{k-j},$
(6)

where $a=g_{\alpha,\tau}^{-1}h^{-2}$ and $\delta
u_{n}:=u_{n+1}-2u_{n}+u_{n-1}$. In a more compact matrix form, this reads

$\displaystyle A\mathbf{u}^{k}=w_{k}\mathbf{u}^{0}-\sum_{j=1}^{k-1}\Delta
w_{j}\mathbf{u}^{k-j},$ (7)

where $\mathbf{u}^{k}:=u(t_{k},x)$, $\Delta w_{j}:=w_{j+1}-w_{j}$, and $A$ is
an $N\times N$ symmetric tridiagonal matrix with $b=1+2a$ on the main
diagonal, $-a$ on the adjacent off-diagonals, and terms in the corners of the
matrix that are specified by the boundary conditions. More precisely, the
matrix $A$ is given by

$\displaystyle A=\begin{bmatrix}c&-a&0&&&\cdots&&d\\\ -a&b&-a&0&&\cdots&&0\\\
0&-a&b&-a&&\cdots&&0\\\ \vdots&&&\ddots&&&&\vdots\\\ 0&&\cdots&&0&-a&b&-a\\\
d&&\cdots&&&0&-a&c\end{bmatrix}\in\mathbb{R}^{N\times N},$ (8)

where $(c,d)=(b,-a)$ for periodic boundary condition, $(c,d)=(b,0)$ for
Dirichlet boundary condition and $(c,d)=(1+a,0)$ for Neumann boundary
condition. The right-hand side of Eq. 7 carries the history of solutions up to
$k-1$.

###  2.2 Variational quantum optimization

Using Dirac notation, we rewrite the iterative fractional diffusion equation
as

$\displaystyle A|\tilde{u}^{k}\rangle=|\tilde{f}^{k-1}\rangle,$ (9)

where we approximate the unnormalized quantum state $|\tilde{u}^{k}\rangle$ at
time $k$ by an ansatz or trial state formed by a set of parameterized
unitaries $R$ and unparameterized entangling unitaries $W$ as,

$\displaystyle|\tilde{u}^{k}\rangle:=r^{k}|u^{k}\rangle=r^{k}\prod_{i=1}^{l}W_{l-i+1}R(\boldsymbol{\theta}_{l-i+1}^{k})|0\rangle,$
(10)

where $r^{k}$ is the norm of the unnormalized state $|\tilde{u}^{k}\rangle$,
$\boldsymbol{\theta}$ is a set of gate parameters and $l$ is the number of
layers in the ansatz. The unnormalized source state $|\tilde{f}^{k-1}\rangle$
is a linear combination of history states in $[0,k-1]$.

Following [11], we define a cost function based on potential energy at time
$k$,

$\displaystyle\mathcal{C}^{k}(r^{k},\boldsymbol{\theta}^{k})=\frac{1}{2}(r^{k})^{2}\langle
u^{k}|A|u^{k}\rangle-\frac{1}{2}r^{k}\left(\langle
u^{k}|\tilde{f}^{k-1}\rangle+\langle\tilde{f}^{k-1}|u^{k}\rangle\right),$ (11)

where

$\displaystyle|\tilde{f}^{k-1}\rangle=w_{k}r^{0}|u^{0}\rangle-\sum_{j=1}^{k-1}\Delta
w_{j}r^{k-j}|u^{k-j}\rangle.$ (12)

With that, $\langle u^{k}|\tilde{f}^{k-1}\rangle$ expands into a list of $k$
measurable overlap terms $\langle u^{k}|u^{[0,k-1]}\rangle$ [23].

Introducing $|f,u\rangle:=2^{-1/2}(|0\rangle|f\rangle+|1\rangle|u\rangle)$, we
have

$\displaystyle\mathcal{C}^{k}(r^{k},\boldsymbol{\theta}^{k})=\frac{1}{2}(r^{k})^{2}\langle
u^{k}|A|u^{k}\rangle-r^{k}\left(w_{k}r^{0}\langle
u^{k},u^{0}\rangle-\sum_{j=1}^{k-1}\Delta w_{j}r^{k-j}\langle
u^{k},u^{k-j}\rangle\right).$ (13)

Here, $\langle f,u\rangle:=\langle f,u|X\otimes\mathbb{I}^{\otimes
n}|f,u\rangle$ is used as a notation shorthand, where $X=|1\rangle\langle
0|+|0\rangle\langle 1|$ is the Pauli-X matrix and
$\mathbb{I}:=\mathbb{I}_{0}+\mathbb{I}_{1}=|0\rangle\langle
0|+|1\rangle\langle 1|$ is the $2\times 2$ identity matrix. Also,
$n=\log_{2}N$ denotes the number of qubits representing $N$ spatial grid
points.

Taking the derivative of $\mathcal{C}^{k}$ with respect to $r^{k}$, the
optimal norm

$\displaystyle r^{k}=\frac{\langle u^{k},\tilde{f}^{k-1}\rangle}{\langle
u^{k}|A|u^{k}\rangle}$ (14)

can be eliminated from the cost function

$\displaystyle\mathcal{C}^{k}(r^{k},\boldsymbol{\theta}^{k})=-\frac{1}{2}\frac{\left({\langle
u^{k},\tilde{f}^{k-1}\rangle}\right)^{2}}{\langle u^{k}|A|u^{k}\rangle},$ (15)

which can be measured using quantum circuits (Fig. 12a,b) and optimized
classically as

$\displaystyle\underset{\boldsymbol{\theta}^{k}}{\arg\min}\
\mathcal{C}^{k}\left(r^{k},\boldsymbol{\theta}^{k}\right),$ (16)

via either gradient-based (evaluate
$\partial_{\boldsymbol{\theta}^{k}}\mathcal{C}^{k}$) or gradient-free methods.
A schematic of the proposed variational quantum algorithm is shown in Fig. 1.

Figure 1: Schematic diagram of variational quantum algorithm for time-
fractional differential equations. The initial solution, encoded by the vector
$\boldsymbol{\theta}^{0}$ containing $nl$ optimized parameters and the norm
$r^{0}$, is stored in classical memory. Subsequent iterations add
$\boldsymbol{\theta}^{k}$ and $r^{k}$ to the growing pool of classically
stored parameters $\boldsymbol{\theta}^{[0,k-1]}$ and norms $r^{[0,k-1]}$,
until the final time-step $M$ is reached. Solutions can be read out
selectively via tomography using up to $\mathcal{O}(n2^{n})$ measurements on
specific $\boldsymbol{\theta}^{k}$ parameterized ansatz corresponding to time
$k$. Sampling complexity can be reduced using partial tomography or shadow
tomography techniques.

##  3 Implementation

###  3.1 Hamiltonian decomposition

The Hamiltonian matrix $A$ (Eq. 8) can be decomposed into

$\displaystyle A=b\mathbb{I}^{\otimes
n}-a\left\\{\underbrace{\mathbb{I}^{\otimes n-1}\otimes
X}_{H_{1}}+\mathcal{S}^{\dagger}\big{[}\underbrace{\mathbb{I}^{\otimes
n-1}\otimes X}_{H_{2}}-\underbrace{\mathbb{I}_{0}^{\otimes n-1}\otimes
X}_{H_{3}}+\underbrace{\mathbb{I}_{0}^{\otimes
n-1}\otimes\mathbb{I}}_{H_{4}}\big{]}\mathcal{S}\right\\}.$ (17)

Since expectation values of the identity operator $\mathbb{I}$ are equal to 1,
i.e. $\langle\phi|\mathbb{I}^{\otimes
n}|\phi\rangle=\langle\phi^{\prime}|\mathbb{I}^{\otimes
n}|\phi^{\prime}\rangle=1$ [11], evaluating the expectation value of $A$
requires only the evaluation of expectation values of the simple terms
($H_{1-2}$ for periodic boundary condition, $H_{1-3}$ Dirichlet boundary
condition and $H_{1-4}$ for Neumann boundary condition). The operator
$\mathcal{S}$ denotes the n-qubit cyclic shift operator,

$\displaystyle\mathcal{S}=\sum^{2^{n}-1}_{i=0}|(i+1)\ \text{mod}\
2^{n}\rangle\langle i|,$ (18)

which can be implemented using relative-phase Toffoli gates, Toffoli gates, a
cnot gate and an X gate [24].

###  3.2 Ansatz initialization

For problems that admit only real solutions, it is preferable to use a real-
amplitude ansatz, represented by Eq. 10 formed by $l$ repeating blocks, each
consisting of a parameterized $R$ layer with one $R_{Y}(\theta)$ gate on each
qubit, followed by unparameterized $W$ layer with a cnot gate between
consecutive qubits [25]:

$\displaystyle|u^{k}\rangle=\prod_{i=l}^{1}\left[\prod_{a=n-1}^{1}C_{a}X_{a+1}\bigotimes_{j=1}^{n}R_{Y}(\theta_{i,j})\right]|0\rangle,$
(19)

where we have employed the product notation with an initial index greater than
the stopping index to signify that the terms in the product are arranged in
decreasing index order [26]. In addition, we consider a real-amplitude ansatz
with circular entanglement, such that

$\displaystyle|u^{k}\rangle=\prod_{i=l}^{1}\left[C_{n}X_{0}\prod_{a=n-1}^{1}C_{a}X_{a+1}\bigotimes_{j=1}^{n}R_{Y}(\theta_{i,j})\right]|0\rangle,\quad
n>2,$ (20)

which has been shown to be efficient for solving partial differential
equations [16]. See Fig. 11 for circuit diagrams of the ansaetze.

To solve Eq. 9, we prepare the initial ansatz state
$|\tilde{u}^{0}\rangle=r^{0}|u(\boldsymbol{\theta}^{0})\rangle$ using
classical optimization on parameters
$\boldsymbol{\theta}=\left(\theta_{i,j}\right)_{(i,j)\in\\{1,\ldots,l\\}\times\\{1,\ldots,n\\}}$
such that

$\displaystyle\boldsymbol{\theta}^{0}\in\underset{\boldsymbol{\theta}^{0}}{\arg\min}\left\|\langle\hat{u}^{0}|u(\boldsymbol{\theta}^{0})\rangle\right\|,$
(21)

where $|\hat{u}^{0}\rangle\equiv\mathbf{u}^{0}/\left\|\mathbf{u}^{0}\right\|$
is the normalized initial classical solution vector and
$r^{0}=\left\|\mathbf{u}^{0}\right\|$ is the initial norm.

###  3.3 Gradient estimation

The gradient of the cost function Eq. 15 can be estimated on quantum computers
using the parameter shift rule. Following [11], here we write down the partial
derivative of the cost function with respect to parameter $\theta^{k}_{i,j}$
indexed by $(i,j)\in[1,nl]$ at $k$, as

$\displaystyle\begin{split}\frac{\partial\mathcal{C}^{k}}{\partial{\theta}_{i}^{k}}=-\frac{1}{2}\left(\frac{\langle
u^{k},\tilde{f}^{k-1}\rangle}{\langle
u^{k}|A|u^{k}\rangle}\right)\left[w_{k}r^{0}\left\langle\frac{\partial
u^{k}}{\partial\theta_{i}^{k}},u^{0}\right\rangle-\sum_{j=1}^{k-1}\Delta
w_{j}r^{k-j}\left\langle\frac{\partial
u^{k}}{\partial\theta_{i}^{k}},u^{k-j}\right\rangle\right]\\\
+\frac{1}{2}\left(\frac{\langle u^{k},\tilde{f}^{k-1}\rangle}{\langle
u^{k}|A|u^{k}\rangle}\right)^{2}\left\langle\frac{\partial
u^{k}}{\partial\theta_{i}^{k}},u^{k}\right\rangle_{A},\end{split}$ (22)

where $\langle f,u\rangle_{A}:=\langle f,u|X\otimes A|f,u\rangle$ is used as a
notational shorthand. The parametric shift is implemented by a $\pi$ rotation
on the $i$-th $R_{Y}$ gate as

$\displaystyle\frac{\partial
u^{k}}{\partial\theta_{i}^{k}}=u^{k}\left(\theta^{\lfloor
i/n\rfloor}_{i\%n}+\pi\right),$ (23)

where $i\%n$ is the remainder $r\in\\{0,1,\ldots,n-1\\}$ obtained when $i$ is
divided by $n$ and $\lfloor i/n\rfloor$ is the integer part of $i/n$. Thus,
for each index $i$ at time $k$, gradient estimation requires at least $k$
overlap measurements.

###  3.4 Time complexity and memory scaling

_Time complexity_. — Here we briefly estimate the time complexity of the
algorithm based on the number of time steps, quantum circuits, gates and shots
required, neglecting classical overheads such as parameter update, iteration,
functional evaluation and initial encoding.

For a number of time steps $M$, the overall time complexity is

$\displaystyle\mathcal{T}\approx
MN_{\mathrm{it}}N_{g}\left[\frac{Ml+(l+n^{2})}{\epsilon^{2}}\right],$ (24)

where $n$ is the number of qubits, $l$ is the number of ansatz layers,
$N_{\mathrm{it}}$ is the mean number of iterations required per time step and
$N_{g}$ is the mean number of evaluations required for gradient estimation per
iteration. Specifically, the terms in square brackets are the gate
complexities for the fractional derivative $Ml$ based on $M/2$ overlap
circuits containing $2l$ parametric gates in depth and for the Hamiltonian
$(l+n^{2})$ based on $\mathcal{O}(1)$ circuits ($2-4$ Eq. 17) containing $l$
parametric gates in depth for the ansatz and $\mathcal{O}(n^{2})$ non-
parametric gates for the shift operator. Each circuit is sampled repeatedly
for $\mathcal{O}(\epsilon^{-2})$ shots per measurement, where $\epsilon$ is
the desired precision. We remark that with quantum amplitude estimation [27],
the number of shots required may be quadratically reduced from
$\mathcal{O}(\epsilon^{-2})$ to $\mathcal{O}(\epsilon^{-1})$. This decrease,
however, entails increased circuit depth, as exemplified by the
$\Omega(\epsilon^{-1})$ circuit depths utilized in [27, 28, 29]. Given these
challenges in reducing circuit depth in quantum amplitude estimation and its
extensions [30, 31, 32, 33], we do not consider this approach in the present
study, deferring it for future exploration.

For $n\sim l$, the gradient-based optimizer $N_{g}\sim\mathcal{O}(nl)$ yields
an overall time complexity

$\displaystyle\mathcal{T}\sim\mathcal{O}\left(\max\left\\{\frac{N_{\mathrm{it}}M^{2}nl^{2}}{\epsilon^{2}},\frac{N_{\mathrm{it}}Mn^{3}l}{\epsilon^{2}}\right\\}\right),$
(25)

scaling as $\mathcal{O}(M^{2})$ for overlap circuits, or $\mathcal{O}(M)$ for
shift operators, the latter scaling applicable to the non-fractional
diffusion/heat equation ($\alpha=1$, [11]). Using simultaneous perturbation
stochastic approximation (SPSA) [34], only two circuit executions are required
per gradient evaluation $N_{g}\sim\mathcal{O}(1)$, independent of the number
of parameters, similar to gradient-free optimization.

The present algorithm is efficient in time complexity with respect to number
of spatial grid points $N(=2^{n})$ [11, 35], neglecting initial encoding
costs. Quantum advantage with respect to time $T$ may be achievable in higher
dimensions (Remark 2 [35], Remark 3 [36]).

_Memory scaling_. — To solve the Caputo derivative Eq. 4, a classical
algorithm assesses $\mathcal{O}(kN)$ parameters from memory at time $k$,
compared to $\mathcal{O}(knl)$ parameters for variational quantum algorithm
Eq. 12. Assuming $n\sim l$, the memory cost of the Caputo sums scales sub-
exponentially as $\mathcal{O}(k\log N)$ (Fig. 1), demonstrating efficient
algorithmic space complexity.

##  4 Numerical experiments

The circuits are implemented in the quantum software framework _Pennylane_
[37] using noiseless statevector emulation on lightning qubit backend.
Parametric updates are performed using the limited-memory Broyden-Fletcher-
Goldfarb-Shanno boxed (L-BFGS-B) optimizer with relative tolerance of
$10^{-6}$ and gradient tolerance of $10^{-6}$. The L-BFGS-B optimizer
evaluates gradients using finite difference without explicit gradient inputs
from Eq. 22. Initial encoding Eq. 21 is performed using SciPy’s minimize
function with an initial null parameter set as $(0,\dots,0)$.

###  4.1 Sub-diffusion equation

We solve for time-fractional $(0<\alpha\leq 1)$ sub-diffusion equation Eq. 1
as an initial value problem with initial condition $u(0,x)=x(1-x)$ and
Dirichlet boundary condition $u(t,0)=u(t,1)=0$ [6].

Figure 2 shows that the variational quantum solutions agree with classical
sub-diffusion solutions in $32\times 32$ space-time domain with qubit-layer
count $(n,l)=(5,4)$, for (a) $\alpha=1$ and (b) $\alpha=0.5$. Figure 3 shows
that low fractional-order solutions tend to diffuse faster for short times,
but slow down for long times, both characteristics of sub-diffusive behavior.

Figure 2: Variational quantum solutions to time-fractional diffusion equation
in $N\times M=32\times 32$ space-time domain with qubit-layer count
$(n,l)=(5,4)$, for (a) $\alpha=1$ and (b) $\alpha=0.5$. Initial and boundary
conditions are $u(0,x)=x(1-x)$ and $u(t,0)=u(t,1)=0$. Colorbar denotes
relative deviation from classical finite difference solution on the same
domain. Figure 3: Variational quantum solutions (as Fig. 2) plotted on (a)
$\alpha-x$ at $t=0.5$, and (b) $t-\alpha$ at $x=0.5$. Colorbar denotes
relative deviation from classical finite difference solution on the same
domain.

The trace error at time $k$

$\displaystyle\epsilon_{\mathrm{tr}}^{k}=\sqrt{1-\left\|\langle\hat{u}^{k}|u(\boldsymbol{\theta}^{k})\rangle\right\|^{2}},$
(26)

can be averaged over the number of time steps $M$ for $0<\alpha\leq 1$, where
$|\hat{u}^{k}\rangle\equiv\mathbf{u}^{k}/\left\|\mathbf{u}^{k}\right\|$ is the
normalized classical solution vector at time $k$.

####  4.1.1 Trace error and optimization

Figure 4a compares the results between the linear real-amplitude hardware-
efficient ansatz Eq. 19 and the circular ansatz Eq. 20, suggesting an improved
solution fidelity for the circular ansatz over the linear ansatz ($n>4$).
According to unitary dependence theory [38], the effect of an additional
circular entangling cnot gate ($C_{n}X_{0}$) on a single layer ansatz ($l=1$)
is to transform the unitary dependence of the first qubit from
$q_{1}:\\{R_{Y}(\theta_{1})\\}$ to $q_{1}:\\{R_{Y}(\theta_{2})\sim
R_{Y}(\theta_{n})\\}$, which effectively increases the connectivity of $q_{1}$
by $n-2$ qubits. This may improve the performance of an ansatz for certain
problems [38].

Figure 4b shows that the number of function evaluations
$\sum_{k=1}^{M}N_{eval}^{k}$, required by the L-BFGS-B gradient-free optimizer
scales sub-linearly with the number of time steps $\mathcal{O}(M^{0.6})$ in
time $t=[0,0.5]$. The overhead costs for classical gradient estimation remains
economical up to $M=2^{6}$.

Figure 4:  (a) Time-averaged trace error $\epsilon_{\mathrm{tr}}$ plots for
linear positive-amplitude hardware-efficient ansatz Eq. 19 (dotted lines) vs.
circular ansatz Eq. 20 (continuous lines), for $0<\alpha\leq 1$. Error bars
denote standard deviation. (b) Number of function evaluations
$\sum_{k=1}^{M}N_{eval}^{k}$ required by L-BFGS-B gradient-free optimizer for
$\alpha=0.5$ (dashed lines) and $\alpha=1$ (continuous lines), scales with
number of time steps $M$ to a fractional power of $0.6$ in time $t=[0,0.5]$.
Indicated $(n,l)$ pairs are selected based on the minimum number of layers $l$
required to encode $|u(\boldsymbol{\theta}^{0})\rangle$ for $n\in\\{3,4,5\\}$.

Figure 5 shows that increasing the number of parameters $nl$ does not improve
the time-averaged trace error $\epsilon_{\mathrm{tr}}$, but instead increases
the number of function evaluations which scales as $\mathcal{O}(nl)$. The
linear scaling suggests that the problem of vanishing gradients due to
excessive parameterization, commonly known as barren plateaus [39], may be
mitigated by short time-steps, such that the initial solution at each time
step is closer to the optimal solution.

Figure 5:  (a) Time-averaged trace error $\epsilon_{\mathrm{tr}}$ is
insensitive to over-parameterization, shown here for $\alpha=0.5$ (dashed
lines) and $\alpha=1$ (continuous lines). Error bars denote standard
deviation. (b) Number of function evaluations $\sum_{k=1}^{M}N_{eval}^{k}$
scales linearly with number of parameters $nl$.

###  4.2 Time-fractional Burgers’ equation

The fractional Burgers’ equation [40, 41] has found applications in shallow
water waves and waves in bubbly liquids [42]. We consider the following 1D
time-fractional Burgers’ equation [43] as an extension to the time-fractional
diffusion equation Eq. 1,

$\displaystyle
D_{t}^{\alpha}u(t,x)=\nu\partial_{xx}u(t,x)-u(t,x)\partial_{x}u(t,x),\quad
0<\alpha\leq 1,$ (27)

where $D_{t}^{\alpha}$ is the Caputo fractional derivative of the $\alpha$-th
order, as previously defined, and $\nu$ is the fluid kinematic viscosity.
Readers interested in quantum computing algorithms for non-fractional Burgers’
equation ($\alpha=1$) are referred to [44, 45, 13].

####  4.2.1 Semi-implicit scheme

Using explicit source terms, we apply central difference on the non-linear
advection term, such that

$\displaystyle\left(1-a\delta\right)u_{n}^{k}=w_{k}u_{n}^{0}-\sum_{j=1}^{k-1}\Delta
w_{j}u_{n}^{k-j}-b\left(u_{n+1}^{k-1}-u_{n-1}^{k-1}\right)u_{n}^{k-1},$ (28)

where $a=\nu/(g_{\alpha,\tau}h^{2})$ and $b=1/(2g_{\alpha,\tau}h)$. In Dirac
notation, the time-fractional Burgers’ equation can be iterated as Eq. 9

$\displaystyle
A|\tilde{u}^{k}\rangle=|\tilde{f}^{k-1}\rangle-b\left(\tilde{\Lambda}_{+}^{k-1}-\tilde{\Lambda}_{-}^{k-1}\right)|\tilde{u}^{k-1}\rangle,$
(29)

where $\tilde{\Lambda}:=\mathrm{diag}(u)$ is a diagonal matrix with $u$ along
its diagonal, and the subscript $\pm$ denotes either incremental or
decremental shift in each eigenstate under periodic boundary conditions, i.e.
$|2^{n}\rangle=|0\rangle$.

With that, the cost function for Burgers’ equation is

$\displaystyle\mathcal{C}^{k}(r^{k},\boldsymbol{\theta}^{k})=-\frac{1}{2}\frac{\left[\langle
u^{k},\tilde{f}^{k-1}\rangle-b\left(\langle
u^{k},\tilde{\Lambda}_{+}^{k-1}\tilde{u}^{k-1}\rangle-\langle
u^{k},\tilde{\Lambda}_{-}^{k-1}\tilde{u}^{k-1}\rangle\right)\right]^{2}}{\langle
u^{k}|A|u^{k}\rangle},$ (30)

which includes two non-linear terms $\langle
u^{k},\tilde{\Lambda}_{\pm}^{k-1}\tilde{u}^{k-1}\rangle$ evaluated using
redundant copies of $k-1$ quantum states [13, 14] on quantum circuit (Fig.
12c).

####  4.2.2 Test case

We test a 1D bi-directional flow using time-fractional Burgers’ equation Eq.
27 with initial condition $u(0,x)=\sin(2\pi x)$ and Dirichlet boundary
condition $u(0,x)=x(1-x)$ and $u(t,0)=u(t,1)=0$.

Fig. 6(a,c,e) show time series solutions in $32\times 32$ space-time domain of
size $(L,T)=(1,1)$ with qubit-layer count $(n,l)=(5,5)$ for fractional order
$\alpha=\\{1,0.8,0.6\\}$ and kinematic viscosity $\nu=0.02$. The Reynolds
number is $Re\equiv 2u_{c}L/\nu=100$, where the characteristic velocity
$u_{c}=1$. In the inviscid case, an initial sine-wave profile tends towards an
N-wave shock solution. As with Fig. 3, reducing fractional order leads to fast
decay at short times and slow convergence at long times [46], shown here for
$\alpha\in\\{1,0.8,0.6\\}$. Insets show the relative deviations from classical
finite difference solution are less than 2%, without any clear dependence on
$\alpha$.

For semi-implicit time-fractional Burgers’ equation, the Courant number $C\sim
u_{c}\tau^{\alpha}/h$ depends on the fractional index $\alpha$. To verify the
dependence of numerical stability on $t^{\alpha}$, the number of time steps
$M$ is quadrupled from 32 to 128 in the time series solutions Fig. 6(b,d,f).
By inspection, time-fractional solutions are sensitive to time-step $\tau$.
Specifically, for $\alpha=0.6$, convergence in solution accuracy requires
$M\gg N$ (see Fig. 6f inset).

Figure 6:  Time-series solutions of 1D Burgers’ equation plotted in time
intervals of $1/32$ in (a,c,e) $N\times M=32\times 32$ and (b,d,f) $N\times
M=32\times 128$ space-time domain of size $(L,T)=(1,1)$ for qubit-layer count
$(n,l)=(5,5)$, kinematic viscosity $\nu=0.02$ ($Re=100$), and fractional order
$\alpha$ (a,b) 1, (c,d) 0.8, and (e,f) 0.6. Insets (a,c,e) show time-averaged
relative deviation from classical finite difference solution on the same
domain and insets (b,d,f) show relative norm convergence with $M$ at times
$t\in\\{1/32,1/16\\}$. Initial and boundary conditions are $u(0,x)=\sin{(2\pi
x)}$ and $u(t,0)=u(t,1)=0$.

####  4.2.3 Trace error and optimization

Figure 7a shows for $\nu\gtrsim 0.04$, the maximum trace error
$\max\epsilon_{\mathrm{tr}}$ for normal Burgers’ equation ($\alpha=1$) is
significantly greater than for the time-fractional Burgers’ equation
($\alpha<1$), the latter which converges to steady state solutions faster
(Fig. 6). Figure 7b shows the required number of function evaluations $\sum
N_{\mathrm{eval}}$ scales with both $\nu$ and $\alpha$, due to the difficulty
of convergence towards unstable solutions driven by the explicit advection
term.

Together, these results suggest a trade-off between solution fidelity and
gradient evaluation costs, depending on fractional index $\alpha$ and
diffusion parameter $\nu$.

Figure 7: $\alpha$-$\nu$ contour plots of (a) maximum trace error
$\max\epsilon_{\mathrm{tr}}$ and (b) number of function evaluations $\sum
N_{\mathrm{eval}}$ for variational quantum solutions of 1D Burgers’ equation
(same as Fig. 6).

###  4.3 Fractional diffusive epidemic model

A natural extension to the nonlinear differential equations involves exploring
coupled systems of differential equations, as exemplified by the susceptible-
exposed-infectious-recovered (SEIR) epidemic model [44]. Fractional-order
epidemic models have garnered significant attention in the wake of the
COVID-19 pandemic [47, 48, 49]. These models include spatio-temporal models
that account for geographical spread of the disease [50, 51, 52]. A minimal
representation of a fractional diffusive SEIR model is given by [53, Eq.
(24)–(29)]:

$\displaystyle D_{t}^{\alpha_{1}}S$ $\displaystyle=\nu_{1}\partial_{xx}S-\beta
IS-\mu S+\Pi,$ (31) $\displaystyle D_{t}^{\alpha_{2}}E$
$\displaystyle=\nu_{2}\partial_{xx}E+\beta IS-(\sigma+\mu)E,$ (32)
$\displaystyle D_{t}^{\alpha_{3}}I$
$\displaystyle=\nu_{3}\partial_{xx}I+\sigma E-(\rho+\mu)I,$ (33)
$\displaystyle D_{t}^{\alpha_{4}}R$ $\displaystyle=\nu_{4}\partial_{xx}R+\rho
I-\mu R,$ (34)

which partitions a population into four distinct cohorts: susceptible
individuals $S$, exposed individuals $E$, infectious individuals $I$, and
those who have recovered $R$, without accounting for vaccination and
quarantine. The symbols $\alpha_{i}$ and $\nu_{i}$, where $i\in\\{1,2,3,4\\}$,
are the fractional order and the diffusion coefficient corresponding to each
respective $\\{S,E,I,R\\}$ cohort in that order, $\Pi$ is the population
influx rate, $\beta$ is the infective rate, $\mu$ is the death rate, $\sigma$
is the progression rate and $\rho$ is the recovery rate.

####  4.3.1 Coupled semi-implicit scheme

The coupled system of Eqs. 31, 32, 33 and 34 is discretized via finite
difference and expressed in the following iterative form with common time-step
$\tau$,

$\displaystyle\left(\beta
I^{k-1}+\mu+g_{\alpha_{1},\tau}-\frac{\nu_{1}}{h^{2}}\delta\right)S^{k}$
$\displaystyle=g_{\alpha_{1},\tau}f_{\alpha_{1}}\left(S^{[0,k-1]}\right)+\Pi,$
(35)
$\displaystyle\left(\sigma+\mu+g_{\alpha_{2},\tau}-\frac{\nu_{2}}{h^{2}}\delta\right)E^{k}$
$\displaystyle=g_{\alpha_{2},\tau}f_{\alpha_{2}}\left(E^{[0,k-1]}\right)+\beta
I^{k-1}S^{k-1},$ (36)
$\displaystyle\left(\rho+\mu+g_{\alpha_{3},\tau}-\frac{\nu_{3}}{h^{2}}\delta\right)I^{k}$
$\displaystyle=g_{\alpha_{3},\tau}f_{\alpha_{3}}\left(I^{[0,k-1]}\right)+\sigma
E^{k-1},$ (37)
$\displaystyle\left(\mu+g_{\alpha_{4},\tau}-\frac{\nu_{4}}{h^{2}}\delta\right)R^{k}$
$\displaystyle=g_{\alpha_{4},\tau}f_{\alpha_{4}}\left(R^{[0,k-1]}\right)+\rho
I^{k-1},$ (38)

where the operator
$f_{\alpha}\left(u^{[0,k-1]}\right)=w_{k}^{[\alpha]}u^{0}-\sum_{j=1}^{k-1}\Delta
w_{j}^{[\alpha]}u^{k-j}$. Since the $\beta IS$ term is strictly positive, it
can be iterated as $\beta I^{k-1}S^{k}$ (Eq. 35). With only an
$\mathcal{O}(1)$ number of cohorts, the coupled system of equations can be
solved using a separate cost function for each cohort at time $k$, i.e.

$\displaystyle\mathcal{C}^{k}_{S}$
$\displaystyle=\frac{(r^{k}_{S})^{2}}{2}\left(\langle
S^{k}|A_{S}|S^{k}\rangle+\beta\langle
S^{k}|\tilde{\Lambda}^{k-1}_{I}|S^{k}\rangle\right)-r^{k}_{S}\left(g_{\alpha_{1},\tau}\langle
S^{k},\tilde{f}_{S,\alpha_{1}}^{k-1}\rangle+\Pi\langle
S^{k},\tilde{+}^{n}\rangle\right),$ (39) $\displaystyle\mathcal{C}^{k}_{E}$
$\displaystyle=\frac{(r^{k}_{E})^{2}}{2}\langle
E^{k}|A_{E}|E^{k}\rangle-r^{k}_{E}\left(g_{\alpha_{2},\tau}\langle
E^{k},\tilde{f}_{E,\alpha_{2}}^{k-1}\rangle+\beta\langle
E^{k},\tilde{\Lambda}^{k-1}_{I}\tilde{S}^{k-1}\rangle\right),$ (40)
$\displaystyle\mathcal{C}^{k}_{I}$
$\displaystyle=\frac{(r^{k}_{I})^{2}}{2}\langle
I^{k}|A_{I}|I^{k}\rangle-r^{k}_{I}\left(g_{\alpha_{3},\tau}\langle
I^{k},\tilde{f}_{I,\alpha_{3}}^{k-1}\rangle+\sigma\langle
I^{k},\tilde{E}^{k-1}\rangle\right),$ (41) $\displaystyle\mathcal{C}^{k}_{R}$
$\displaystyle=\frac{(r^{k}_{R})^{2}}{2}\langle
R^{k}|A_{R}|R^{k}\rangle-r^{k}_{R}\left(g_{\alpha_{4},\tau}\langle
R^{k},\tilde{f}_{R,\alpha_{4}}^{k-1}\rangle+\rho\langle
R^{k},\tilde{I}^{k-1}\rangle\right),$ (42)

where $|\tilde{+}^{n}\rangle\equiv\sqrt{2^{n}}|+\rangle^{\otimes n}$ is the
(unnormalized) equal superposition over all computational basis states, and
each subscript label $\mathbb{S}\in\\{S,E,I,R\\}$ corresponds to the indicated
cohort. The Hamiltonian matrices $\in\mathbb{R}^{N\times N}$ are
$A_{S}=(\mu+g_{\alpha_{1},\tau})\mathcal{I}-\nu_{1}/h^{2}\mathcal{L}$,
$A_{E}=(\sigma+\mu+g_{\alpha_{2},\tau})\mathcal{I}-\nu_{2}/h^{2}\mathcal{L}$,
$A_{I}=(\rho+\mu+g_{\alpha_{3},\tau})\mathcal{I}-\nu_{3}/h^{2}\mathcal{L}$ and
$A_{R}=(\mu+g_{\alpha_{4},\tau})\mathcal{I}-\nu_{4}/h^{2}\mathcal{L}$, where
$\mathcal{L}$ is a symmetric tridiagonal matrix with $-2$ along its main
diagonal and $1$ along the adjacent off-diagonals, and
$\mathcal{I}=\mathbb{I}^{\otimes n}$ is an $n$-qubit identity matrix. Terms
such as $\langle
S(\boldsymbol{\theta}^{k})|\Lambda^{k-1}_{I}|S(\boldsymbol{\theta}^{k})\rangle$
can be measured using the circuit shown in Fig. 12d.

####  4.3.2 Basic reproduction number

The spread of a transmissible disease is governed by an important quantity
known as the basic reproduction number [54],

$\displaystyle\mathcal{R}=\frac{\beta\sigma\Pi}{\mu(\mu+\sigma)(\mu+\rho)},$
(43)

which is defined, in epidemiological terms, as the average number of secondary
cases produced by one infected individual in a susceptible population.
Depending on the value of $\mathcal{R}$, the number of infectious individuals
may increase resulting in an epidemic ($\mathcal{R}>1$), or decrease resulting
in the disease dying out ($\mathcal{R}<1$).

####  4.3.3 Test case

Consider a set of $2^{n}$ spatially connected populations, each with a small
infected cohort and a much larger susceptible cohort. Following [55], we
assume the following set of SEIR parameters: $\Pi=750$ [53], $\mu=0.03325$,
$\beta=5.1\times 10^{-6}$, $\sigma=0.17$, and $\rho=0.1109$, leading to
$\mathcal{R}\approx 0.66$. Individuals are assumed to move around within a
region $[0,1]$ according to Fick’s law, with diffusion coefficents
$\nu_{1}=\nu_{2}=10^{-3}$ and $\nu_{3}=\nu_{4}=5\times 10^{-4}$. Initial
cohort sizes are $S(0,x)=22,500$ ($\sim\Pi/\mu$ for near-steady populations
[55]), $I(0,x)=20e^{-2x}$ (for spatial dependence [53]) and $E(0,x)=R(0,x)=1$.
Neumann boundary conditions are applied, i.e.
$\partial_{x}\mathbb{S}(t,0)=\partial_{x}\mathbb{S}(t,L)=0$, where
$\mathbb{S}\in\\{S,E,I,R\\}$.

In a $16\times 32$ domain of size $(L,T)=(1,100)$ with $(n,l)=(4,5)$, Fig. 8
verifies that $\mathcal{R}<1$ results in a monotonic decrease (a,c,e) in the
infectious cohort $I$, whereas $\mathcal{R}>1$ results in an increase (b,d,f)
in $I$ before eventual decrease at long times (not shown). A decrease in the
fractional index $\alpha$ flattens the spatiotemporal profile of $I$, in such
a way that $I$ decreases more slowly for $\mathcal{R}<1$, and increases more
slowly for $\mathcal{R}>1$.

Figure 8: Infectious $I$ cohort (actual values) in $16\times 32$ space-time
domain of size $(L,T)=(1,100)$ with qubit-layer count $(n,l)=(4,5)$ for each
cohort, depending on basic reproduction number (a,c,e) $\mathcal{R}=0.66$ and
(b,d,f) $\mathcal{R}=1.33$, and fractional index (a,b) $\alpha=1$, (c,d)
$\alpha=0.8$ and (e,f) $\alpha=0.6$.

Similar time-fractional memory effects also apply to other cohorts, namely
susceptible $S$, exposed $E$ and recovered $R$, as shown in Fig. 9 at the
start of an epidemic ($\mathcal{R}>1$). The apparent smoothness of the
solutions, as presented, affirms the potential of the iterative solver in
other applications involving coupled fractional-order derivatives, such as
predator-prey population models [56] and Oldroyd-B viscoelastic flows [2].

Figure 9: (a,b) Susceptible $S$, (c,d) exposed $E$ and (e,f) recovered $R$
cohorts (actual values) in $16\times 32$ space-time domain of size
$(L,T)=(1,100)$ with qubit-layer count $(n,l)=(4,5)$ for each cohort at the
start of an epidemic ($\mathcal{R}=1.33>1$). Fractional index (a,c,e)
$\alpha=1$ and (b,d,f) $\alpha=0.8$.

##  5 Hardware noise

Near-term quantum devices are subjected to hardware noise which significantly
impairs the performance of variational quantum algorithms [57, 58], despite
their inherent resilience to coherent errors [59]. To explore the effect of
noise, we employ a realistic noise model sampled from the IBM Nairobi quantum
hardware device (Falcon r5.11H, version 1.3.3) applied to the _Aer Simulator_
backend provided by the _Pennylane-Qiskit_ plugin. Circuits are transpiled
using default settings with level 3 optimization and SABRE routing. The
optimizer of choice is the Simultaneous Perturbation Stochastic Approximation
(SPSA), which is known for its exceptional scalability and performance under
noisy conditions [57, 58]. We apply the SPSAOptimizer class on _Pennylane_
using default hyperparameters for 200 iterations, equivalent to 400 circuit
evaluations. Data readout is assumed to be noiseless via statevector
emulation.

Here, we solve minimally for the sub-diffusion problem (Section 4.1) in
$N\times M=4\times 2$ using qubit-layer count $(n,l)=(2,1)$ for $\alpha=1$.
Initial conditions are $u(0,x)=x(1-x)$ and boundary conditions are periodic
instead of Dirichlet. Figure 10 shows the effects of hardware noise model on
_Aer Simulator_ based on 40 instances with fixed 10,000 shots per circuit
evaluation. Of note here is the effect of noise on the evaluation of the norm,
which directly affects the weights of the cost function $\mathcal{C}$ for the
subsequent time step. For fair evaluation of each time-step, we restore the
norm to the classical error-free value before each time-step.

Figure 10:  Hardware noise model tests across 40 instances on _Aer Simulator_
using SPSA optimizer set at 200 iterations per time step. (a) Effect of noise
model on overlap and Hamiltonian measurements results in shifting of (b) the
optimized cost function $\mathcal{C}$. (c) Box plots of norm fidelity
$\varphi_{\mathrm{norm}}$ and trace error $\epsilon_{\mathrm{tr}}$ with noise
model in time steps $k=\\{1,2\\}$. Reference values without noise model are
shown as $\varphi_{\mathrm{norm}}\approx 1$ and $\epsilon_{\mathrm{tr}}\approx
0$. (d) Comparing $|u^{(1)}\rangle$ and $|u^{(2)}\rangle$ with noise model
(box plots) to noiseless simulation (line circles). Inset shows hardware test
results across 20 instances obtained from the IBMQ Mumbai 27 qubit device
using Qiskit Runtime set at 200 SPSA iterations for a single time step.

Figure 10a shows that averaged overlap measurements are shifted by the noise
model to $\langle u^{1},u^{0}\rangle=0.686\pm 0.049$ and $\langle
u^{2},u^{1}\rangle=0.682\pm 0.049$, from approximately $1$ for measurements
without noise model. The observed relative noise error is approximately
$30\%$, which is significantly greater than the error from Hamiltonian
measurements $\langle u^{k}|A|u^{k}\rangle$ with noise model (less than $10\%$
by inspection, Fig. 10b). This can be attributed to greater gate errors from
the Toffoli (controlled cnot) gates found in controlled entangling layers of
the overlap circuit. In addition, the fractional noise error of the cost
function Eq. 15

$\displaystyle\left(\eta_{C}\right)^{2}=\left(2\eta_{O}\right)^{2}+\left(\eta_{H}\right)^{2}$
(44)

depends on $\eta_{O}$, the fractional noise error of the overlap measurements
$\langle u^{k},u^{k-1}\rangle$, and $\eta_{H}$, the fractional noise error of
the Hamiltonian measurements $\langle u^{k}|A|u^{k}\rangle$. Given
$\eta_{O}/\eta_{H}\approx 3$, the proportion of noise error of the cost
function due to the noise of overlap measurement is $\sqrt{36/37}\approx
0.986$. This means that the overall performance of the algorithm depends
almost entirely on the noise error of overlap measurements, and very little on
that of Hamiltonian measurements.

Figure 10c shows the norm fidelity

$\displaystyle\varphi_{\mathrm{norm}}^{k}=\left|\frac{r^{k}}{\|u^{k}\|}\right|,$
(45)

where $r^{k}$ is the measured norm Eq. 14 and $\|u^{k}\|$ is the classical
norm, and trace error $\epsilon_{\mathrm{tr}}^{k}$ Eq. 26 at time step
$k\in\\{1,2\\}$. Here, the mean norm fidelities
$\overline{\varphi}_{\mathrm{norm}}^{k}\approx 0.65$ are nearly identical for
$k\in\\{1,2\\}$ due to the classical norm reset, whereas the mean trace errors
are $\overline{\epsilon}_{\mathrm{tr}}^{1}\approx 0.039$ and
$\overline{\epsilon}_{\mathrm{tr}}^{2}\approx 0.052$.

Figure 10d shows the effects of hardware noise on the quantum solution states
$|u^{k}\rangle$ over 40 instances. Although the degree of data scatter at each
node (width of box plot) is of a similar magnitude to the difference between
consecutive classical solutions, the mean readout states
$|\overline{u}^{k}\rangle$ are in general agreement with the respective
classical solutions. Further tests were conducted on actual hardware, in this
case, the IBMQ Mumbai 27 qubit device across 20 instances (see inset) using
Qiskit Runtime sampler and estimator primitives for overlap and Hamiltonian
measurements respectively, set at 200 SPSA iterations for a single time step.

##  6 Discussion

In this study, we proposed and implemented a hybrid variational quantum
algorithm for solving time-fractional diffusive differential equations in
engineering applications, including non-linear equations such as Burgers’
equation (Section 4.2) and coupled systems of equations such as those in
epidemic models (Section 4.3). Our approach is not only efficient in time
complexity, but also offers advantages in terms of memory utilization. While
the classical memory cost of numerically solving these equations scales
linearly with the number of spatial grid points $N$, the quantum memory cost
scales only logarithmically with $N$. This helps to mitigate the global
dependence problem in fractional calculus, where previous solutions must be
repeatedly accessed from memory causing the classical solver to have
prohibitive memory costs when the spatial range is large. In terms of time
complexity, the present iterative scheme is limited by the quadratic scaling
with the number of time steps; this has some implications in the numerical
stability of non-linear fractional equations (see Fig. 6).

The choice of the ansatz in the form of a parameterized circuit is reflected
in its expressibility and entangling capability performance [60] of a
variational quantum algorithm [9]. Hence, as a matter of heuristics, the real-
amplitude ansatze used in this study (Eqs. 19 and 20) are by no means the most
efficient designs for such applications. Designing and optimizing ansatz
architecture, including utilizing adaptive methods [61] or leveraging
geometric tools [62, 63, 64], remains an active area of research [65, 66, 67].

While the hybrid approach offers notable advantages, it is not exempt from the
challenges that are commonly encountered in variational quantum algorithms.
First, unlike a classical computer, the cost of loading initial conditions as
classical data as an _n_ -qubit quantum state could be prohibitive [68].
Encoding schemes, such as parametric optimization (Eq. 21), may incur an
exponential cost in primitive operations since the encoded _n_ -qubit state
occupies a space of dimension $O(2^{n})$. More efficient encoding techniques
are currently being developed; for example, solutions may be represented in
the form of finite Fourier series, using the quantum Fourier transform as part
of the Fourier series loader introduced by [69] and employed by [14]. The use
of the Walsh series instead of the Fourier series may yield an even more
efficient sub-exponential scaling in encoding [70]. In addition, extracting
the complete solution at each of the $M$ time steps would require full readout
of the statevector via quantum state tomography. This experimental procedure
is notorious for its exponential cost [71], which is a bottleneck to quantum
advantage as pointed out by [68]. Therefore, where possible, applications
should consider a partial readout of the (spatial) solutions only at few
selected time points of interest to minimize the expensive overhead due to
tomography, or measuring only integrated quantities like mean, variance or
other moments. In these cases, robust classical shadow tomography techniques
[72, 73, 74, 75, 76] can be employed to further reduce the required sample
complexity.

In closing, the hybrid algorithm investigated in this work presents a nascent
but promising pathway to solving fractional partial differential equations of
the integro-differential form, subjected to potentially prohibitive space or
memory requirements. More work needs to be done to address issues related to
the trainability of variational quantum algorithms, such as barren plateaus
and narrow gorges [39, 77], as well as practical hardware implementation,
including error mitigation [78] and correction.

##  Appendix A Quantum circuit diagrams

In this appendix, we give explicit circuit diagrams for some of the circuits
we used in this paper. Real-amplitude hardware-efficient ansaetze are shown in
Fig. 11. Circuits used for measuring expectations and overlaps are shown in
Fig. 12.

$\Qcircuit<EMAIL_ADDRESS>1}\mbox{layer 2}\mbox{layer $l$}\\\
~{}\\\
\lstick{|0\rangle}\qw\gate{R_{Y}(\theta_{1,1})}\ctrl{1}\qw\qw\qw\qw\qw\qw\gate{R_{Y}(\theta_{2,1})}\ctrl{1}\qw\qw\qw\qw\qw\qw\ldots\gate{R_{Y}(\theta_{l,1})}\ctrl{1}\qw\qw\qw\qw\qw\qw\\\
\lstick{|0\rangle}\qw\gate{R_{Y}(\theta_{1,2})}\targ\ctrl{1}\qw\qw\qw\qw\qw\gate{R_{Y}(\theta_{2,2})}\targ\ctrl{1}\qw\qw\qw\qw\qw\ldots\gate{R_{Y}(\theta_{l,2})}\targ\ctrl{1}\qw\qw\qw\qw\qw\\\
\lstick{|0\rangle}\qw\gate{R_{Y}(\theta_{1,3})}\qw\targ\ctrl{1}\qw\qw\qw\qw\gate{R_{Y}(\theta_{2,3})}\qw\targ\ctrl{1}\qw\qw\qw\qw\ldots\gate{R_{Y}(\theta_{l,3})}\qw\targ\ctrl{1}\qw\qw\qw\qw\\\
\targ\qw\targ\qw\targ\qw\\\ ~{}\\\ ~{}\\\
\lstick{\vdots}\vdots\ddots\ddots\vdots\ddots\\\ ~{}\\\ ~{}\\\ ~{}\\\
\ctrl{1}\qw\ctrl{1}\qw\ctrl{1}\qw\\\
\lstick{|0\rangle}\qw\gate{R_{Y}(\theta_{1,n})}\qw\qw\qw\qw\targ\qw\qw\gate{R_{Y}(\theta_{2,n})}\qw\qw\qw\qw\targ\qw\qw\ldots\gate{R_{Y}(\theta_{l,n})}\qw\qw\qw\qw\targ\qw\qw\gategroup{3}{3}{14}{8}{.6em}{.}\gategroup{3}{11}{14}{16}{.6em}{.}\gategroup{3}{25}{14}{30}{.6em}{.}}$
(a)

$\Qcircuit<EMAIL_ADDRESS>1}\mbox{layer 2}\mbox{layer $l$}\\\
~{}\\\
\lstick{|0\rangle}\qw\gate{R_{Y}(\theta_{1,1})}\ctrl{1}\qw\qw\qw\qw\qw\targ\qw\qw\gate{R_{Y}(\theta_{2,1})}\ctrl{1}\qw\qw\qw\qw\qw\targ\qw\qw\ldots\gate{R_{Y}(\theta_{l,1})}\ctrl{1}\qw\qw\qw\qw\qw\targ\qw\qw\\\
\lstick{|0\rangle}\qw\gate{R_{Y}(\theta_{1,2})}\targ\ctrl{1}\qw\qw\qw\qw\qw\qw\qw\gate{R_{Y}(\theta_{2,2})}\targ\ctrl{1}\qw\qw\qw\qw\qw\qw\qw\ldots\gate{R_{Y}(\theta_{l,2})}\targ\ctrl{1}\qw\qw\qw\qw\qw\qw\qw\\\
\lstick{|0\rangle}\qw\gate{R_{Y}(\theta_{1,3})}\qw\targ\ctrl{1}\qw\qw\qw\qw\qw\qw\gate{R_{Y}(\theta_{2,3})}\qw\targ\ctrl{1}\qw\qw\qw\qw\qw\qw\ldots\gate{R_{Y}(\theta_{l,3})}\qw\targ\ctrl{1}\qw\qw\qw\qw\qw\qw\\\
\targ\qw\targ\qw\targ\qw\\\ ~{}\\\ ~{}\\\
\lstick{\vdots}\vdots\ddots\ddots\vdots\ddots\\\ ~{}\\\ ~{}\\\ ~{}\\\
\ctrl{1}\qw\ctrl{1}\qw\ctrl{1}\qw\\\
\lstick{|0\rangle}\qw\gate{R_{Y}(\theta_{1,n})}\qw\qw\qw\qw\targ\qw\ctrl{-11}\qw\qw\gate{R_{Y}(\theta_{2,n})}\qw\qw\qw\qw\targ\qw\ctrl{-11}\qw\qw\ldots\gate{R_{Y}(\theta_{l,n})}\qw\qw\qw\qw\targ\qw\ctrl{-11}\qw\qw\gategroup{3}{3}{14}{10}{.6em}{.}\gategroup{3}{13}{14}{20}{.6em}{.}\gategroup{3}{29}{14}{36}{.6em}{.}}$
(b)

Figure 11:  (a) Real-amplitude hardware-efficient ansatz (Eq. 19) and (b)
circular real-amplitude ansatz (Eq. 20).

$\Qcircuit@C=0.3em@R=0.05em@!{\lstick{|0\rangle}&\gate{H}\ctrlo{1}\ctrl{1}\gate{H}\meter\qw\\\
\lstick{|0\rangle^{\otimes
n}}\qw\gate{U(\boldsymbol{\theta}^{k})}\gate{U^{j}}\qw\qw\qw}$ (a)

$\Qcircuit<EMAIL_ADDRESS>\lstick{|0\rangle^{\otimes
n}}\qw\gate{U(\boldsymbol{\theta}^{k})}\qw\gate{\mathcal{I}/\mathcal{S}}\qw\meter\qw}$
(b)

$\Qcircuit@C=0.3em@R=0.05em@!{\lstick{|0\rangle}&\gate{H}\ctrlo{1}\ctrl{2}\ctrl{1}\ctrl{1}\gate{H}\meter\qw\\\
\lstick{|0\rangle^{\otimes
n}}\qw\gate{U(\boldsymbol{\theta}^{k})}\qw\gate{U^{k-1}}\ctrl{1}\qw\qw\qw\\\
\lstick{|0\rangle^{\otimes
n}}\qw\qw\gate{U^{k-1}}\gate{\mathcal{S}/\mathcal{S}^{\dagger}}\targ\qw\qw\qw\\\
}$ (c)

$\Qcircuit@C=0.3em@R=0.05em@!{\lstick{|0\rangle}&\gate{H}\ctrl{2}\ctrl{1}\gate{H}\meter\qw\\\
\lstick{|0\rangle^{\otimes
n}}\gate{S(\boldsymbol{\theta}^{k})}\qw\ctrl{1}\qw\qw\qw\\\
\lstick{|0\rangle^{\otimes n}}\qw\gate{I^{k-1}}\targ\qw\qw\qw\\\ }$ (d)

$\Qcircuit<EMAIL_ADDRESS>\\\ \mathcal{H}\\\ \lstick{|0\rangle^{\otimes
n}}\qw\gate{U(\boldsymbol{\theta}^{k})}\gate{U^{k-1}}\gate{\mathcal{I}/\mathcal{S}}\meter\qw}$
(e)

Figure 12:  Circuits for measuring (a) overlap $\langle
u(\boldsymbol{\theta}^{k}),u^{j}\rangle$ and (b) expectation $\langle
U(\boldsymbol{\theta}^{k})|A|U(\boldsymbol{\theta}^{k})\rangle$ (Eq. 15), (c)
nonlinear overlap $\langle u^{k},\Lambda_{\pm}^{k-1}u^{k-1}\rangle$ (Eq. 30),
(d) expectation $\langle
S(\boldsymbol{\theta}^{k})|\Lambda^{k-1}_{I}|S(\boldsymbol{\theta}^{k})\rangle$
(Eq. 39) and (e) higher-order overlap $\langle u^{k},Bu^{k-1}\rangle$ (Eq.
49). $U(\boldsymbol{\theta}^{k})$ denotes an ansatz with parameters
$\boldsymbol{\theta}$ optimized at time $k$ and $U^{j}$ denotes an ansatz with
parameters fixed at time $j$. $\mathcal{I}/\mathcal{S}$ denotes an optional
shift operator (Eq. 18) and $\mathcal{S}/\mathcal{S}^{\dagger}$ denotes either
incremental or decremental shift operator. Operators $S$ and $I$ are the
ansaetze representing the respective variables.

##  Appendix B Higher-order Crank-Nicolson scheme

The quadratic scaling of quantum complexity in time steps Eq. (25) motivates
higher-order time-stepping methods in order to limit time discretization and
its errors. The Crank-Nicolson method is a second-order finite-difference
scheme that retains the stability advantage associated with implicit methods
[79]. Using the trapezoidal rule, the finite difference of the spatial
derivative Eq. (5) takes the midpoint between forward and backward Euler
schemes, yielding

$\displaystyle\partial_{xx}u(t_{k},x_{n})\simeq\frac{1}{2h^{2}}\left[\left(u_{n+1}^{k}-2u_{n}^{k}+u_{n-1}^{k}\right)+\left(u_{n+1}^{k-1}-2u_{n}^{k-1}+u_{n-1}^{k-1}\right)\right].$
(46)

Equating to the finite difference of the Caputo derivative, we obtain

$\displaystyle\left(1-\frac{a}{2}\delta\right)u_{n}^{k}=\frac{a}{2}\delta
u_{n}^{k-1}+w_{k}u_{n}^{0}-\sum_{j=1}^{k-1}\Delta w_{j}u_{n}^{k-j},$ (47)

or in matrix form,

$\displaystyle
A\mathbf{u}^{k}=B\mathbf{u}^{k-1}+w_{k}\mathbf{u}^{0}-\sum_{j=1}^{k-1}\Delta
w_{j}\mathbf{u}^{k-j},$ (48)

where $A=\mathbb{I}-B\in\mathbb{R}^{N\times N}$,
$B=(a/2)\mathcal{L}\in\mathbb{R}^{N\times N}$, $\mathcal{L}$ is a symmetric
tridiagonal matrix with $-2$ along its main diagonal and $1$ along the
adjacent off-diagonals, and $\mathbb{I}$ is the identity matrix. This can be
iterated using variational quantum optimization Eq. 9 as before, with an
updated cost function

$\displaystyle\mathcal{C}^{k}=-\frac{1}{2}\frac{\left(r^{k-1}\langle
u^{k},Bu^{k-1}\rangle+\langle u^{k},\tilde{f}^{k-1}\rangle\right)^{2}}{\langle
u^{k}|A|u^{k}\rangle},$ (49)

which includes a new term $\langle u^{k},Bu^{k-1}\rangle$ to be measured as a
sum of observables ($2-4$, Eq. (17)), using the quantum circuit (Fig. 12e).

##  Appendix C Truncated Caputo integral

Since the coefficient $\Delta w_{j}$ of the $(k-j)$-th term decreases
monotonically with $j$ for constant fractional power index $\alpha$, the
$\mathcal{O}(M^{2})$ time complexity scaling Eq. (25) can be limited to
$\mathcal{O}(M)$ by setting an upper bound to the limits of the Caputo
integral. This is done by replacing the Caputo finite difference derivative in
Eq. (4) with

$\displaystyle D_{t}^{\alpha}u(t_{k},x_{n})\simeq
g_{\alpha,\tau}\left[w_{1}u_{n}^{k}+\sum_{j=1}^{\min\\{k,\xi\\}}\left(w_{j+1}-w_{j}\right)u_{n}^{k-j}\right],$
(50)

where $\xi\in[1,M]$ and $w_{k+1}=0$. For $\xi<k$, the truncation error thus
incurred is

$\displaystyle\varepsilon_{\operatorname{tr}}(\alpha,\xi)\sim\mathcal{O}\left[g_{\alpha,\tau}(\Delta
w_{\xi+1})r^{k-(\xi+1)}\right],$ (51)

where $\Delta w_{\xi+1}=w_{\xi+2}-w_{\xi+1}$. A similar trade-off between
solution fidelity and time complexity applies to classical finite difference
schemes.

## Acknowledgements

We acknowledge the use of IBM Quantum services for this work. The views
expressed are those of the authors, and do not reflect the official policy or
position of IBM or the IBM Quantum team. This research is supported by the
National Research Foundation, Singapore and the Agency for Science, Technology
and Research (A*STAR) under the Quantum Engineering Programme
(NRF2021-QEP2-02-P03), and A*STAR C230917003. This project is supported by the
Singapore Ministry of Health’s National Medical Research Council through the
Programme for Research in Epidemic Preparedness and Response (PREPARE), under
Environmental Transmission & Mitigation Co-operative (PREPARE-CS1-2022-004).
DEK acknowledges funding support from the A*STAR Central Research Fund (CRF)
Award for Use-Inspired Basic Research. DP acknowledges support from the
Ministry of Education Singapore, under the grant MOE-T2EP50120-0019.

## Author Declarations

The authors have no conflicts to disclose.

## Data Availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * [1] Taohua Liu and Muzhou Hou.  A Fast Implicit Finite Difference Method for Fractional Advection-Dispersion Equations with Fractional Derivative Boundary Conditions.  Advances in Mathematical Physics, 2017(4):1–8, 2017.
  * [2] Fang Wang and Yu Wang.  A Finite Difference Method for Solving Unsteady Fractional Oldroyd-B Viscoelastic Flow Based on Caputo Derivative.  Advances in Mathematical Physics, 2023:1–22, apr 2023.
  * [3] Olaniyi Samuel Iyiola and FD Zaman.  A fractional diffusion equation model for cancer tumor.  AIP Advances, 4(10), 2014.
  * [4] V.E. Lynch, B.A. Carreras, D. Del-Castillo-Negrete, K.M. Ferreira-Mejias, and H.R. Hicks.  Numerical methods for the solution of partial differential equations of fractional order.  Journal of Computational Physics, 192(2):406–421, dec 2003.
  * [5] Changpin Li and Fanhai Zeng.  The Finite Difference Methods for Fractional Ordinary Differential Equations.  Numerical Functional Analysis and Optimization, 34(2):149–179, feb 2013.
  * [6] Diego A. Murio.  Implicit finite difference approximation for time fractional diffusion equations.  Computers and Mathematics with Applications, 56(4):1138–1145, 2008\.
  * [7] Yumin Lin and Chuanju Xu.  Finite difference/spectral approximations for the time-fractional diffusion equation.  Journal of Computational Physics, 225(2):1533–1552, aug 2007.
  * [8] Kishor Bharti, Alba Cervera-Lierta, Thi Ha Kyaw, Tobias Haug, Sumner Alperin-Lea, Abhinav Anand, Matthias Degroote, Hermanni Heimonen, Jakob S. Kottmann, Tim Menke, Wai-Keong Mok, Sukin Sim, Leong-Chuan Kwek, and Alán Aspuru-Guzik.  Noisy intermediate-scale quantum algorithms.  Reviews of Modern Physics, 94(1):015004, feb 2022.
  * [9] M. Cerezo, Andrew Arrasmith, Ryan Babbush, Simon C. Benjamin, Suguru Endo, Keisuke Fujii, Jarrod R. McClean, Kosuke Mitarai, Xiao Yuan, Lukasz Cincio, and Patrick J. Coles.  Variational quantum algorithms.  Nature Reviews Physics, 3(9):625–644, Sep 2021.
  * [10] Hui-Min Li, Zhi-Xi Wang, and Shao-Ming Fei.  Variational quantum algorithms for Poisson equations based on the decomposition of sparse Hamiltonians.  Physical Review A, 108(3):032418, 2023.
  * [11] Yuki Sato, Ruho Kondo, Satoshi Koide, Hideki Takamatsu, and Nobuyuki Imoto.  Variational quantum algorithm based on the minimum potential energy for solving the Poisson equation.  Physical Review A, 104(5):052409, 2021.
  * [12] Mazen Ali and Matthias Kabel.  Performance study of variational quantum algorithms for solving the Poisson equation on a quantum computer.  Phys. Rev. Appl., 20:014054, Jul 2023.
  * [13] Michael Lubasch, Jaewoo Joo, Pierre Moinier, Martin Kiffner, and Dieter Jaksch.  Variational quantum algorithms for nonlinear problems.  Physical Review A, 101(1):010301, jan 2020.
  * [14] Abhijat Sarma, Thomas W Watts, Mudassir Moosa, Yilian Liu, and Peter L McMahon.  Quantum variational solving of nonlinear and multi-dimensional partial differential equations.  arXiv preprint arXiv:2311.01531, 2023.
  * [15] Wei-Bin Ewe, Dax Enshan Koh, Siong Thye Goh, Hong-Son Chu, and Ching Eng Png.  Variational quantum-based simulation of waveguide modes.  IEEE Transactions on Microwave Theory and Techniques, 70(5):2517–2525, 2022.
  * [16] Fong Yew Leong, Dax Enshan Koh, Wei-Bin Ewe, and Jian Feng Kong.  Variational quantum simulation of partial differential equations: applications in colloidal transport.  International Journal of Numerical Methods for Heat & Fluid Flow, 33(11):3669–3690, July 2023.
  * [17] Dieter Jaksch, Peyman Givi, Andrew J Daley, and Thomas Rung.  Variational quantum algorithms for computational fluid dynamics.  AIAA journal, 61(5):1885–1894, 2023.
  * [18] Pete Rigas.  Variational quantum algorithm for measurement extraction from the Navier-Stokes, Einstein, Maxwell, Boussniesq-type, Lin-Tsien, Camassa-Holm, Drinfeld-Sokolov-Wilson, and Hunter-Saxton equations.  arXiv preprint arXiv:2209.07714, 2022.
  * [19] Fong Yew Leong, Wei-Bin Ewe, and Dax Enshan Koh.  Variational quantum evolution equation solver.  Scientific reports, 12(1):10817, jun 2022.
  * [20] Kenji Kubo, Koichi Miyamoto, Kosuke Mitarai, and Keisuke Fujii.  Pricing multi-asset derivatives by variational quantum algorithms.  arXiv preprint arXiv:2207.01277, 2022.
  * [21] Francesco Mainardi, Gianni Pagnini, and Rudolf Gorenflo.  Some aspects of fractional diffusion equations of single and distributed order.  Applied Mathematics and Computation, 187(1):295–305, 2007.
  * [22] Sadeq Taha Abdulazeez and Mahmut Modanli.  Solutions of fractional order pseudo-hyperbolic telegraph partial differential equations using finite difference method.  Alexandria Engineering Journal, 61(12):12443–12451, 2022.
  * [23] Alternatively, $|\tilde{f}^{k-1}\rangle$ can be prepared separately through post-selection of linear combination of unitaries [80] using a quantum multiplexer, to be measured directly as one single overlap term $\langle u^{k}|\tilde{f}^{k-1}\rangle$.
  * [24] Dmitri Maslov.  Advantages of using relative-phase Toffoli gates with an application to multiple control Toffoli optimization.  Physical Review A, 93(2):022311, feb 2016.
  * [25] Hedayat Alghassi, Amol Deshmukh, Noelle Ibrahim, Nicolas Robles, Stefan Woerner, and Christa Zoufal.  A variational quantum algorithm for the Feynman-Kac formula.  Quantum, 6:730, June 2022.
  * [26] For example, in decreasing index order, we have $\prod_{i=3}^{1}{A_{i}=A_{3}A_{2}A_{1}}$.
  * [27] Gilles Brassard, Peter Hoyer, Michele Mosca, and Alain Tapp.  Quantum amplitude amplification and estimation.  Contemporary Mathematics, 305:53–74, 2002.
  * [28] Yohichi Suzuki, Shumpei Uno, Rudy Raymond, Tomoki Tanaka, Tamiya Onodera, and Naoki Yamamoto.  Amplitude estimation without phase estimation.  Quantum Information Processing, 19(2), January 2020.
  * [29] Scott Aaronson and Patrick Rall.  Quantum approximate counting, simplified.  In Symposium on simplicity in algorithms, pages 24–32. SIAM, 2020\.
  * [30] Tudor Giurgica-Tiron, Iordanis Kerenidis, Farrokh Labib, Anupam Prakash, and William Zeng.  Low depth algorithms for quantum amplitude estimation.  Quantum, 6:745, June 2022.
  * [31] Dax Enshan Koh, Guoming Wang, Peter D. Johnson, and Yudong Cao.  Foundations for Bayesian inference with engineered likelihood functions for robust amplitude estimation.  Journal of Mathematical Physics, 63(5):052202, 05 2022.
  * [32] Guoming Wang, Dax Enshan Koh, Peter D Johnson, and Yudong Cao.  Minimizing estimation runtime on noisy quantum computers.  PRX Quantum, 2(1):010346, 2021.
  * [33] Peter D Johnson, Alexander A Kunitsa, Jérôme F Gonthier, Maxwell D Radin, Corneliu Buda, Eric J Doskocil, Clena M Abuan, and Jhonathan Romero.  Reducing the cost of energy estimation in the variational quantum eigensolver algorithm with robust amplitude estimation.  arXiv preprint arXiv:2203.07275, 2022.
  * [34] James C Spall.  Multivariate stochastic approximation using a simultaneous perturbation gradient approximation.  IEEE transactions on automatic control, 37(3):332–341, 1992.
  * [35] Yuki Sato, Ruho Kondo, Ikko Hamamura, Tamiya Onodera, and Naoki Yamamoto.  Hamiltonian simulation for time-evolving partial differential equation by scalable quantum circuits.  arXiv preprint arXiv:2402.18398, 2024.
  * [36] Junpeng Hu, Shi Jin, Nana Liu, and Lei Zhang.  Quantum Circuits for partial differential equations via Schrödingerisation.  arXiv preprint arXiv:2403.10032, 2024.
  * [37] Ville Bergholm, Josh Izaac, Maria Schuld, Christian Gogolin, M Sohaib Alam, Shahnawaz Ahmed, Juan Miguel Arrazola, Carsten Blank, Alain Delgado, Soran Jahangiri, et al.  Pennylane: Automatic differentiation of hybrid quantum-classical computations.  arXiv preprint arXiv:1811.04968, 2018.
  * [38] Zixuan Hu and Sabre Kais.  The unitary dependence theory for characterizing quantum circuits and states.  Communications Physics, 6(1), 2023.
  * [39] M. Cerezo, Akira Sone, Tyler Volkoff, Lukasz Cincio, and Patrick J. Coles.  Cost function dependent barren plateaus in shallow parametrized quantum circuits.  Nature Communications, 12(1), dec 2021.
  * [40] Yufeng Xu and Om P. Agrawal.  Numerical solutions and analysis of diffusion for new generalized fractional Burgers equation.  Fractional Calculus and Applied Analysis, 16(3):709–736, sep 2013\.
  * [41] N. Sugimoto.  Burgers equation with a fractional derivative; hereditary effects on nonlinear acoustic waves.  Journal of Fluid Mechanics, 225:631–653, apr 1991.
  * [42] Roberto Garra.  Fractional-calculus model for temperature and pressure waves in fluid-saturated porous rocks.  Physical Review E, 84(3):036605, sep 2011.
  * [43] Vijitha Mukundan and Ashish Awasthi.  Linearized Implicit Numerical Method for Burgers’ Equation.  Nonlinear Engineering, 5(4):219–234, jan 2016.
  * [44] Jin Peng Liu, Herman Øie Kolden, Hari K. Krovi, Nuno F. Loureiro, Konstantina Trivisa, and Andrew M. Childs.  Efficient quantum algorithm for dissipative nonlinear differential equations.  Proceedings of the National Academy of Sciences of the United States of America, 118(35), aug 2021.
  * [45] Furkan Oz, Rohit K. S. S. Vuppala, Kursat Kara, and Frank Gaitan.  Solving Burgers’ equation with quantum computing.  Quantum Information Processing, 21(1):30, jan 2022.
  * [46] Tayyaba Akram, Muhammad Abbas, Muhammad Bilal Riaz, Ahmad Izani Ismail, and Norhashidah Mohd Ali.  An efficient numerical technique for solving time fractional Burgers equation.  Alexandria Engineering Journal, 59(4):2201–2220, aug 2020.
  * [47] Zeeshan Ali, Faranak Rabiei, Mohammad M Rashidi, and Touraj Khodadadi.  A fractional-order mathematical model for COVID-19 outbreak with the effect of symptomatic and asymptomatic transmissions.  Eur. Phys. J. Plus, 137:395, 2022.
  * [48] Shahram Rezapour, Hakimeh Mohammadi, and Mohammad Esmael Samei.  SEIR epidemic model for COVID-19 transmission by Caputo derivative of fractional order.  Advances in Difference Equations, 2020:490, 12 2020.
  * [49] Muhammad Altaf Khan and Abdon Atangana.  Modeling the dynamics of novel coronavirus (2019-nCov) with fractional derivative.  Alexandria Engineering Journal, 59(4):2379–2389, aug 2020.
  * [50] Benedict Barnes, Martin Anokye, Mohammed Muniru Iddrisu, Bismark Gawu, and Emmanuel Afrifa.  The analytic solutions of the fractional-order model for the spatial epidemiology of the COVID-19 infection.  Advances in Mathematical Physics, 2023:1–19, 5 2023.
  * [51] Nushrat Nazia, Zahid Ahmad Butt, Melanie Lyn Bedard, Wang-Choi Tang, Hibah Sehar, and Jane Law.  Methods used in the spatial and spatiotemporal analysis of COVID-19 epidemiology: A systematic review.  International journal of environmental research and public health, 19, 7 2022.
  * [52] Bisong Hu, Pan Ning, Jingyu Qiu, Vincent Tao, Adam Thomas Devlin, Haiying Chen, Jinfeng Wang, and Hui Lin.  Modeling the complete spatiotemporal spread of the COVID-19 epidemic in mainland china.  International Journal of Infectious Diseases, 110:247–257, 9 2021\.
  * [53] Adedapo Ismaila Alaje and Morufu Oyedunsi Olayiwola.  A fractional-order mathematical model for examining the spatiotemporal spread of COVID-19 in the presence of vaccine distribution.  Healthcare Analytics, 4:100230, 12 2023.
  * [54] Pauline van den Driessche.  Reproduction numbers of infectious disease models.  Infectious Disease Modelling, 2(3):288–303, aug 2017.
  * [55] Brice Kammegne, Kayode Oshinubi, Oluwatosin Babasola, Olumuyiwa James Peter, Olumide Babatope Longe, Roseline Bosede Ogunrinde, Emmanuel Olurotimi Titiloye, Roseline Toyin Abah, and Jacques Demongeot.  Mathematical Modelling of the Spatial Distribution of a COVID-19 Outbreak with Vaccination Using Diffusion Equation.  Pathogens, 12(1):88, jan 2023.
  * [56] H. Jafari, R.M. Ganji, N.S. Nkomo, and Y.P. Lv.  A numerical study of fractional order population dynamics model.  Results in Physics, 27:104456, 2021.
  * [57] Aidan Pellow-Jarman, Ilya Sinayskiy, Anban Pillay, and Francesco Petruccione.  A comparison of various classical optimizers for a variational quantum linear solver.  Quantum Information Processing, 20(6):202, jun 2021.
  * [58] Xavier Bonet-Monroig, Hao Wang, Diederick Vermetten, Bruno Senjean, Charles Moussa, Thomas Bäck, Vedran Dunjko, and Thomas E. O’Brien.  Performance comparison of optimization methods on variational quantum algorithms.  Physical Review A, 107(3):032407, mar 2023.
  * [59] Enrico Fontana, Nathan Fitzpatrick, David Muñoz Ramo, Ross Duncan, and Ivan Rungger.  Evaluating the noise resilience of variational quantum algorithms.  Physical Review A, 104(2), nov 2020.
  * [60] Sukin Sim, Peter D. Johnson, and Alán Aspuru-Guzik.  Expressibility and Entangling Capability of Parameterized Quantum Circuits for Hybrid Quantum-Classical Algorithms.  Advanced Quantum Technologies, 2(12):1–15, 2019.
  * [61] Harper R. Grimsley, Sophia E. Economou, Edwin Barnes, and Nicholas J. Mayhall.  An adaptive variational algorithm for exact molecular simulations on a quantum computer.  Nature Communications, 10(1):3007, jul 2019.
  * [62] Amara Katabarwa, Sukin Sim, Dax Enshan Koh, and Pierre-Luc Dallaire-Demers.  Connecting geometry and performance of two-qubit parameterized quantum circuits.  Quantum, 6:782, August 2022.
  * [63] Tobias Haug, Kishor Bharti, and M.S. Kim.  Capacity and quantum geometry of parametrized quantum circuits.  PRX Quantum, 2:040309, Oct 2021.
  * [64] Tobias Haug and M. S. Kim.  Natural parametrized quantum circuit.  Phys. Rev. A, 106:052611, Nov 2022.
  * [65] Junhan Qin.  Review of ansatz designing techniques for variational quantum algorithms.  Journal of Physics: Conference Series, 2634(1):012043, nov 2023\.
  * [66] Jia-Bin You, Dax Enshan Koh, Jian Feng Kong, Wen-Jun Ding, Ching Eng Png, and Lin Wu.  Exploring variational quantum eigensolver ansatzes for the long-range XY model.  arXiv preprint arXiv:2109.00288, 2021.
  * [67] Jules Tilly, Hongxiang Chen, Shuxiang Cao, Dario Picozzi, Kanav Setia, Ying Li, Edward Grant, Leonard Wossnig, Ivan Rungger, George H Booth, et al.  The variational quantum eigensolver: a review of methods and best practices.  Physics Reports, 986:1–128, 2022.
  * [68] Scott Aaronson.  Read the fine print.  Nature Physics, 11(4):291–293, apr 2015.
  * [69] Mudassir Moosa, Thomas W Watts, Yiyou Chen, Abhijat Sarma, and Peter L McMahon.  Linear-depth quantum circuits for loading Fourier approximations of arbitrary functions.  Quantum Science and Technology, 9(1):015002, jan 2024.
  * [70] Julien Zylberman and Fabrice Debbasch.  Efficient quantum state preparation with Walsh series.  Phys. Rev. A, 109:042401, Apr 2024.
  * [71] Ryan O’Donnell and John Wright.  Efficient quantum tomography.  In Proceedings of the forty-eighth annual ACM symposium on Theory of Computing, pages 899–912, 2016.
  * [72] Hsin-Yuan Huang, Richard Kueng, and John Preskill.  Predicting many properties of a quantum system from very few measurements.  Nature Physics, 16(10):1050–1057, 2020.
  * [73] Kaifeng Bu, Dax Enshan Koh, Roy J. Garcia, and Arthur Jaffe.  Classical shadows with Pauli-invariant unitary ensembles.  npj Quantum Information, 10:6, 2024.
  * [74] Senrui Chen, Wenjun Yu, Pei Zeng, and Steven T. Flammia.  Robust shadow estimation.  PRX Quantum, 2:030348, Sep 2021.
  * [75] Dax Enshan Koh and Sabee Grewal.  Classical Shadows With Noise.  Quantum, 6:776, August 2022.
  * [76] Daniel Grier, Hakop Pashayan, and Luke Schaeffer.  Sample-optimal classical shadows for pure states.  arXiv preprint arXiv:2211.11810, 2022.
  * [77] Andrew Arrasmith, Zoë Holmes, M Cerezo, and Patrick J Coles.  Equivalence of quantum barren plateaus to cost concentration and narrow gorges.  Quantum Science and Technology, 7(4):045015, aug 2022.
  * [78] Zhenyu Cai, Ryan Babbush, Simon C. Benjamin, Suguru Endo, William J. Huggins, Ying Li, Jarrod R. McClean, and Thomas E. O’Brien.  Quantum error mitigation.  Rev. Mod. Phys., 95:045005, Dec 2023.
  * [79] Umair Ali, Farah Aini Abdullah, and Ahmad Izani Ismail.  Crank-Nicolson finite difference method for two-dimensional fractional sub-diffusion equation.  Journal of Interpolation and Approximation in Scientific Computing, 2017(2):18–29, 2017.
  * [80] Andrew M Childs and Nathan Wiebe.  Hamiltonian simulation using linear combinations of unitary operations.  Quantum Information & Computation, 12(11&12):901–924, 2012.


# JSR-Net: A Deep Network for Joint Spatial-Radon Domain CT Reconstruction
from incomplete data

###### Abstract

CT image reconstruction from incomplete data, such as sparse views and limited
angle reconstruction, is an important and challenging problem in medical
imaging. This work proposes a new deep convolutional neural network (CNN),
called JSR-Net, that jointly reconstructs CT images and their associated Radon
domain projections. JSR-Net combines the traditional model-based approach with
deep architecture design of deep learning. A hybrid loss function is adapted
to improve the performance of the JSR-Net making it more effective in
protecting important image structures. Numerical experiments demonstrate that
JSR-Net outperforms some latest model-based reconstruction methods, as well as
a recently proposed deep model.

Index Terms—  Joint spatial-Radon domain reconstruction, Sparse-view CT,
Limited angle CT, Convolutional neural networks, Deep learning

##  1 Introduction

Deep learning is widely used in Natural language processing, speech
recognition, computer vision and many other fields in recent years [1]. Deep
models have surpassed traditional handcrafted models and human experts in many
tasks in imaging science such as image classification [2, 3]. Recently, deep
learning has attracted much attention from the medical imaging community [4].

In this work, we focus on the problem of CT image reconstruction from
incomplete data [5], i.e., sparse views and limited angle reconstruction
problems. Traditional CT image reconstruction algorithms include filtered
back-projection(FBP) [5], algebraic reconstruction technique(ART) [6] and
model-based iterative reconstruction[7, 8]. However, these traditional methods
cannot effectively utilize large image data sets, which limits their
performance in various image reconstruction tasks.

The CNN based models are capable of learning multi-scale image features from
large data sets with a cascade of simple modules. Although training a deep CNN
can be time-consuming, it is very efficient to use in the validation stage. In
[9], the authors parameterized the proximal operator in the alternating
direction method of multipliers (ADMM) algorithm using CNN for magnetic
resonance imaging (MRI). Their deep model benefits from both the handcrafted
image reconstruction model and CNN and is superior to some of the
regularization based MR image reconstruction models. In [10], the authors
unrolled the primal-dual hybrid gradient(PDHG) algorithm to form a feed-
forward deep network for CT image reconstruction. Then, the proximal operator
and hyper-parameters in the PDHG algorithm were both approximated by CNN
learned from the training data set. This new approach provides significant
improvements in imaging quality when compared to the TV based variational
model. A similar idea was also adopted to learn an iterative scheme for the
nonlinear inverse problem of CT imaging [11].

To suppress the artifacts induced by the incomplete data and noise, [12]
proposed a joint spatial-Radon domain reconstruction(JSR) model for sparse
view CT imaging. This model can better guarantee the data consistency in the
Radon domain which leads to better image quality. Similar idea to the JSR
model is used later for positron emission tomography [13]. To further protect
image features, [14] adopted a data-driven regularization in the JSR model.
More recently, a re-weighting strategy was introduced to the JSR model to
reduce the metal artifacts in multi-chromatic energy CT [15].

In view of the advantages of CNN for unrolled iterative scheme learning and
the JSR model in CT image reconstruction, we propose a new deep architecture
combining deep CNN with JSR model. This newly proposed model will be referred
to as the JSR-Net. In order to fully exploit the advantage of the JSR-Net, a
hybrid loss function is introduced to protect important image structures in
the reconstructed images further. The JSR-Net is then validated and compared
on large CT image data set with simulated projection data for both sparse view
and limited angle problem.

##  2  Method

In this section, we first review the JSR model and its associated iterative
optimization algorithm. Then, we present the proposed JSR-Net that is inspired
by the iterative optimization algorithm.

###  2.1 The JSR Model and Algorithm

The JSR model proposed in [12] reads as follows

$\displaystyle\min_{\bm{u},{\bm{f}}}\mathcal{F}({\bm{u}},{\bm{f}},{\bm{Y}})+\|\bm{\lambda}_{1}\bm{W}_{1}{\bm{u}}\|_{1,2}+\|\bm{\lambda}_{2}\bm{W}_{2}{\bm{f}}\|_{1,2},$
(1)

where the data fidelity term is defined by
$\mathcal{F}({\bm{u}},{\bm{f}},{\bm{Y}})=\frac{1}{2}\|R_{\bm{\Gamma}^{c}}({\bm{f}}-{\bm{Y}})\|^{2}+\frac{\alpha}{2}\|R_{\bm{\Gamma}}({\bm{\mathcal{P}}}{\bm{u}}-{\bm{f}})\|^{2}+\frac{\gamma}{2}\|R_{\bm{\Gamma}^{c}}({\bm{\mathcal{P}}}{\bm{u}}-{\bm{Y}})\|^{2}$
, and $R_{\bm{\Gamma}}$ is the restriction operator with respect to the domain
with missing data $\bm{\Gamma}$. Here, $\bm{\Gamma}^{c}$ denotes the
complement of $\bm{\Gamma}$ and it indicates the region with measured data.
${\bm{\mathcal{P}}}$ is the Radon transform, ${\bm{Y}}$ is the measured
projection data. The parameters $\alpha$ and $\gamma$ are properly chosen to
balance the data consistency in Radon and spatial domain. $\bm{W}_{i},i=1,2,$
is the sparsity promoting transform such as tight wavelet frame transform
[16], and $\|\cdot\|$ is the $\ell_{2}$ norm. $\bm{\lambda}_{i},i=1,2,$ is the
multi-indexed regularization parameter to balance the sparsity prior and the
data fidelity. The special norm $\|\cdot\|_{1,2}$ was firstly introduced in
[17].

The optimization problem in (1) can be solved by the ADMM algorithm with
multi-blocks by introducing the augmented Lagrangian function :

$\displaystyle\mathcal{L}(\bm{X},\bm{d},\bm{b})=\mathcal{F}(\bm{X},{\bm{Y}})+\mathcal{R}(\bm{d})-\langle\bm{b},\bm{W}\bm{X}-\bm{d}\rangle+\frac{1}{2}\|\bm{WX}-\bm{d}\|_{\bm{\mu}}^{2}$
(2)

where  $\bm{d}=\left(\begin{matrix}\bm{d}_{1}\\\
\bm{d}_{2}\end{matrix}\right),$ $\bm{X}=\left(\begin{matrix}\bm{u}\\\
\bm{f}\end{matrix}\right),$ $\bm{\mu}=\left(\begin{matrix}\mu_{1}\\\
\mu_{2}\end{matrix}\right)>\bm{0},$ and
$\bm{W}=\left(\begin{matrix}\bm{W}_{1}\\\ \bm{W}_{2}\end{matrix}\right)$ .
Here,  $\bm{b}=\left(\begin{matrix}\bm{b}_{1}\\\
\bm{b}_{2}\end{matrix}\right)$ is the dual variable. Furthermore,
$\mathcal{R}(\bm{d})=\|\bm{\lambda}_{1}\cdot\bm{d}_{1}\|_{1,2}+\|\bm{\lambda}_{2}\cdot\bm{d}_{2}\|_{1,2}$
and
$\|\bm{WX}-\bm{d}\|_{\bm{\mu}}^{2}=\mu_{1}\|{\bm{W}}_{1}{\bm{u}}-{\bm{d}}_{1}\|^{2}+\mu_{2}\|{\bm{W}}_{2}{\bm{f}}-{\bm{d}}_{2}\|^{2}$
. Then the ADMM scheme with primal and dual variables updated alternatively
can be computed with closed-form formula for each subproblem. More details on
the transition from (1) to (2) can refer to [12]. We omit the detailed
derivation and present the algorithm in Algorithm 1.

Algorithm 1 JSR algorithm for (2)

1:Initialization: ${\bm{b}}_{1}^{0}={\bm{b}}_{2}^{0}=0$,

2:while stop criterion is not met do

3:update ${\bm{u}}$:  $\displaystyle{\bm{u}}^{k+1}$ $\displaystyle=$
$\displaystyle\mathcal{A}^{-1}\left[\alpha{\bm{\mathcal{P}}}^{\top}R_{\bm{\Gamma}}{\bm{f}}^{k}+\mathcal{B}+\mu_{1}{\bm{W}}_{1}^{\top}({\bm{d}}_{1}^{k}-{\bm{b}}_{1}^{k})\right]$
$\displaystyle{\bm{d}}_{1}^{k+1}$ $\displaystyle=$
$\displaystyle\mathcal{T}_{\bm{\lambda}_{1}/\mu_{1}}({\bm{W}}_{1}{\bm{u}}^{k+1}+{\bm{b}}_{1}^{k})$
$\displaystyle{\bm{b}}_{1}^{k+1}$ $\displaystyle=$
$\displaystyle{\bm{b}}_{1}^{k}+({\bm{W}}_{1}{\bm{u}}^{k+1}-{\bm{d}}_{1}^{k+1})$
(3) where $\mathcal{A}=P^{\top}(\alpha R_{\bm{\Gamma}}+\gamma
R_{\bm{\Gamma}^{c}}){\bm{\mathcal{P}}}+\mu_{1}$ and
$\mathcal{B}=\gamma{\bm{\mathcal{P}}}^{\top}R_{\bm{\Gamma}^{c}}{\bm{Y}}$

4:update ${\bm{f}}$:  $\displaystyle{\bm{f}}^{k+1}$ $\displaystyle=$
$\displaystyle\mathcal{C}^{-1}\left[\alpha
R_{\bm{\Gamma}}{\bm{\mathcal{P}}}{\bm{u}}^{k+1}+\mathcal{D}+\mu_{2}{\bm{W}}_{2}^{\top}({\bm{d}}_{2}^{k}-{\bm{b}}_{2}^{k})\right]$
$\displaystyle{\bm{d}}_{2}^{k+1}$ $\displaystyle=$
$\displaystyle\mathcal{T}_{\bm{\lambda}_{2}/\mu_{2}}({\bm{W}}_{2}{\bm{f}}^{k+1}+{\bm{b}}_{2}^{k})$
$\displaystyle{\bm{b}}_{2}^{k+1}$ $\displaystyle=$
$\displaystyle{\bm{b}}_{2}^{k}+({\bm{W}}_{2}f^{k+1}-{\bm{d}}_{2}^{k+1})$ (4)
where $\mathcal{C}=\alpha R_{\bm{\Gamma}}+R_{\bm{\Gamma}^{c}}+\mu_{2}$ and
$\mathcal{D}=R_{\bm{\Gamma}^{c}}{\bm{Y}}.$

5:end while

6:Output: ${\bm{u}}^{\ast}$

###  2.2 JSR-Net

The idea of the design of JSR-Net is to unroll Algorithm 1 and approximate
some of the operators using neural networks. In this way, we form a deep feed-
forward network which is the proposed JSR-Net. To be more precise, we use CNNs
to approximate the inverse operators $\mathcal{A}^{-1}$ and
$\mathcal{C}^{-1}$, and the thresholding operator $\mathcal{T}$ in (3) and
(4). For simplicity, we denote these CNNs as $\mathcal{N}(\cdot;\Theta)$.
Then, we obtain the JSR-Net as shown in Algorithm 2, where all the parameters
of the CNNs and $\alpha,\gamma,\mu_{1},\mu_{2}$ are trainable parameters.

We adopt the same CNN architecture (with a different set of learnable
parameters that varies for each iteration $k$) to approximate the matrix
inversions in both of the ${\bm{u}}^{k}$\- and ${\bm{f}}^{k}$-subproblem. The
thresholding operator that appears in the ${\bm{d}}_{1}^{k}$\- and
${\bm{d}}_{2}^{k}$-subproblem are approximated by a CNN with the same
architecture but different from that of the ${\bm{u}}^{k}$\- and
${\bm{f}}^{k}$-subproblem. To be more precise, the matrix inversion in the
${\bm{u}}^{k}$-subproblem (same for the ${\bm{f}}^{k}$-subproblem) is
approximated by a CNN with a 3 level DenseNet architecture [18] followed by an
LM-ResNet structure [19]. At each level of DenseNet, the input is processed by
the composition of a Convolution(Conv) layer followed by PReLU activation [3].
While for the thresholding operators of the ${\bm{d}}_{1}^{k}$\- and
${\bm{d}}_{2}^{k}$-subproblem, the CNN is composed of 3 convolution layers
that comprises consecutive operations “Conv $\to$ PReLU $\to$ Conv ”. Note
that, it is possible to approximate ${\bm{d}}_{1}^{k}$ and ${\bm{d}}_{2}^{k}$
by choosing a single convolution layer with ReLU activation.

Algorithm 2 JSR-Net

1:Initialization:
${\bm{b}}_{1},{\bm{b}}_{2},{\bm{u}},{\bm{f}},{\bm{W}}_{1},{\bm{W}}_{2},\mathcal{N}(\cdot)$

2:for k=0:N do

3:update ${\bm{u}}$:  $\displaystyle{\bm{u}}^{k+1}$ $\displaystyle=$
$\displaystyle\mathcal{N}_{u}(\left[{\bm{\mathcal{P}}}^{\top}R_{\bm{\Gamma}}{\bm{f}}^{k},\mathcal{B},{\bm{W}}_{1}^{\top}({\bm{d}}_{1}^{k}-{\bm{b}}_{1}^{k})\right];\Theta^{k}_{u})$
$\displaystyle{\bm{d}}_{1}^{k+1}$ $\displaystyle=$
$\displaystyle\mathcal{N}_{d_{1}}({\bm{W}}_{1}{\bm{u}}^{k+1}+{\bm{b}}_{1}^{k};\Theta^{k}_{d_{1}})$
$\displaystyle{\bm{b}}_{1}^{k+1}$ $\displaystyle=$
$\displaystyle{\bm{b}}_{1}^{k}+({\bm{W}}_{1}{\bm{u}}^{k+1}-{\bm{d}}_{1}^{k+1})$
(5) where
$\mathcal{B}=\gamma{\bm{\mathcal{P}}}^{\top}R_{\bm{\Gamma}^{c}}{\bm{Y}}$

4:update ${\bm{f}}$:  $\displaystyle{\bm{f}}^{k+1}$ $\displaystyle=$
$\displaystyle\mathcal{N}_{f}(\left[R_{\bm{\Gamma}}{\bm{\mathcal{P}}}{\bm{u}}^{k+1},\mathcal{D},{\bm{W}}_{2}^{\top}({\bm{d}}_{2}^{k}-{\bm{b}}_{2}^{k})\right];\Theta_{f}^{k})$
$\displaystyle{\bm{d}}_{2}^{k+1}$ $\displaystyle=$
$\displaystyle\mathcal{N}_{d_{2}}({\bm{W}}_{2}{\bm{f}}^{k+1}+{\bm{b}}_{2}^{k};\Theta^{k}_{d_{2}})$
$\displaystyle{\bm{b}}_{2}^{k+1}$ $\displaystyle=$
$\displaystyle{\bm{b}}_{2}^{k}+({\bm{W}}_{2}{\bm{f}}^{k+1}-{\bm{d}}_{2}^{k+1})$
(6) where $\mathcal{D}=R_{\bm{\Gamma}^{c}}{\bm{Y}}$

5:end for

6:Output: ${\bm{u}}^{\ast}$

##  3 Network training

###  3.1 Loss function

To prevent over-smoothing, we adopt a hybrid loss function that blends SSIM
loss, MSE and semantic loss to guide the training process. Our structure-
semantic-$\ell_{2}$ (SS2) hybrid loss function is defined as

$\mathcal{L}_{SS2}=\theta_{1}\mathcal{L}_{SSIM}+\mathcal{L}_{MSE}+\theta_{3}\mathcal{L}_{sem},$
(7)

where  $\mathcal{L}_{SSIM}=\sum(1-SSIM({\bm{u}}_{rec},{\bm{u}}_{truth}))$ is
the error summation of SSIM over mini batch with respect to the reconstructed
image ${\bm{u}}_{rec}$ and ground truth ${\bm{u}}_{truth}$.
$\mathcal{L}_{MSE}=\mathcal{L}_{S}+\theta_{2}\mathcal{L}_{R}$ is the
$\ell_{2}$ loss in spatial domain and Radon domain defined by
$\mathcal{L}_{S}=\sum\|{\bm{u}}_{rec}-{\bm{u}}_{truth}\|^{2}$, and
$\mathcal{L}_{R}=\sum\|R_{\bm{\Gamma}^{c}}{\bm{\mathcal{P}}}({\bm{u}}_{rec}-{\bm{u}}_{truth})\|^{2}$,
respectively. $\mathcal{L}_{sem}$ is the $\ell_{2}$ norm of the difference of
the level sets (i.e. image contents) of ${\bm{u}}_{rec}$ and
${\bm{u}}_{truth}$ as
$\mathcal{L}_{sem}=\|sem({\bm{u}}_{rec})-sem({\bm{u}}_{truth})\|^{2}$.

In the training stage, the loss function (7) is minimized by the stochastic
gradient descent algorithm with adaptive moment estimation(ADAM)[20]. The ADAM
optimizer is adopted with a cosine annealing learning rate, at training step
$t$, defined by $\eta=\frac{0.001}{2}(1+cos(\pi\frac{t}{t_{max}}))$; the
second raw moment parameter $\beta=0.99$ and the gradient norm is clipped by
$1.0$. The maximum training step $t_{max}$ is set to $10000$ for all the
experiments. The training is conducted on Tensorflow 1.3.0 with a Titan Xp GPU
with memory 10.75G.

###  3.2 Training database

The real clinical data “the 2016 NIH-AAPM-Mayo Clinic Low Dose CT Grand
Challenge” [21] from Mayo Clinic is used for the training and validation. The
data set contains abdomen CT images from five patients. It contains 1684
slices of images with $512\times 512$ pixels. We randomly select 50 slices for
validation and the rest for training.

##  4 Experiments

The simulated CT imaging system has 1024 detectors with fan beam geometry. For
the sparse view CT reconstruction, 90 available views uniformly distribute
over $360^{\circ}$. For the limited angle CT reconstruction, 150 views are
collected which are measured by $1^{\circ}$ per scanning step. The measured
projection data, denoted by ${\bm{Y}}$, is contaminated by additive white
Gaussian noise with variance $\sigma=0.05$ (for simplicity). A more general
model with Poisson-Gaussian noise data [22] will be considered in the future.

We compare the reconstruction results of JSR-Net to FBP, PD-Net [10], JSR
model(Algorithm 1). To have a fair comparison, we choose $N=5$ in JSR-Net
(Algorithm 2) so that it has the similar number of trainable parameters as PD-
Net. The parameters in Algorithm 1 are empirically tuned for optimal
performance. All the hyper-parameters in JSR-Net are fixed for the simulations
of both sparse view and limited angle.

The source code of the PD-Net is publicly available. We only modify the
forward imaging operator ${\bm{\mathcal{P}}}$ by composing it with the data
restriction operator $R_{\bm{\Gamma}^{c}}$ to obtain the final imaging system.
All the other hyper-parameters are set as default. To demonstrate the
advantage of the JSR-Net, we also train the PD-Net using the SS2 hybrid loss
function (7).

###  4.1 Sparse view CT image reconstruction

We evaluate the performance of the JSR-Net for sparse view CT image
reconstruction. Fig. 1 shows an example of the ground truth and reconstructed
images from FBP, PD-Net, JSR model(Algorithm 1) and JSR-Net(Algorithm 2). The
gray scale window for FBP is set to [0, 0.5], and [0, 1] for the rest. The
error maps are shown in Fig. 1(h) and Fig. 1(i).

The reconstruction results in Fig. 1 shows that both the JSR-Net and PD-Net
outperform FBP and can better preserve sharp features than the JSR model. The
PD-Net with the proposed SS2 hybrid loss function (7) has better performance
than using merely the $\ell_{2}$ loss which was originally adopted by [10].
The JSR-Net with the SS2 hybrid loss function has the best overall
performance, and it performs better than using only the $\ell_{2}$ loss in
term of suppressing streak artifacts. Note that, since the semantic term
${\mathcal{L}}_{sem}$ is not crucial and empirically has little performance
improvement for sparse view problem, we set $\theta_{3}=0$ in this experiment.
Error maps (with grayscale window [-0.2, 0.2]) in Fig. 1(h) and (i) show that
JSR-Net can better reconstruct image details than PD-Net.

###  4.2 Limited angle CT image reconstruction

We evaluate the performance of JSR-Net for limited angle CT image
reconstruction. Fig. 2 shows one of the example from test set with the ground
truth and the reconstructed images from FBP, PD-Net [10], JSR model(Algorithm
1) and JSR-Net(Algorithm 2).

The gray scale window is set the same as before except that the error map
window is set to [-0.4,0.4]. Fig. 2 shows that both the PD-Net and JSR-Net
outperform FBP and the JSR model, especially at the top region of the image.
The use of the SS2 hybrid loss can improve the results from PD-Net over the
$\ell_{2}$ loss but only mildly. This is because the PD-Net does not have a
joint spatial Radon recovery mechanism. In contrast, the JSR-Net can
significantly benefit from the SS2 hybrid loss and it generates the best
results overall compared models. Streak artifacts are observable in Fig. 2 (c)
and (f).

(a) |  (b) |  (c)
---|---|---
(d) |  (e) |  (f)
(g) |  (h) |  (i)
Fig. 1: Sparse view CT image reconstruction. (a)Ground truth; (b)FBP; (c)PD-Net, $\ell_{2}$; (d)PD-Net, SS2; (e)JSR model; (f) JSR-Net, $\ell_{2}$; (g)JSR-Net, SS2; (h)Error map of PD-Net, SS2; (i)Error map of JSR-Net, SS2. (a) |  (b) |  (c)
---|---|---
(d) |  (e) |  (f)
(g) |  (h) |  (i)
Fig. 2: Limited angle CT image reconstruction. (a)Ground truth; (b)FBP; (c)PD-Net, $\ell_{2}$; (d)PD-Net, SS2; (e)JSR model; (f) JSR-Net, $\ell_{2}$; (g)JSR-Net, SS2; (h)Error map of PD-Net; SS2, (i)Error map of JSR-Net, SS2. Table 1:  Incomplete data CT image reconstruction.  Tasks | Models | Qual. Meas.
---|---|---
SSIM | PSNR | NRMSE | MSE
Sparse view CT | FBP | 0.6173 | 17.25 | 1.078 | 0.0189
PD-Net, $\ell_{2}$ | 0.8709 | 28.54 | 0.1453 | 0.0014
PD-Net, SS2 | 0.8844 | 30.68 | 0.1134 | 0.0009
JSR model | 0.8088 | 26.64 | 0.1866 | 0.0022
JSR-Net,$\ell_{2}$ | 0.8271 | 27.68 | 0.1604 | 0.0017
JSR-Net,SS2 | 0.9081 | 31.59 | 0.1022 | 0.0007
Limited angle CT | FBP | 0.4826 | 15.91 | 1.5143 | 0.0257
PD-Net, $\ell_{2}$ | 0.8778 | 26.43 | 0.1852 | 0.0023
PD-Net, SS2 | 0.88 | 27.44 | 0.1648 | 0.0018
JSR model | 0.8317 | 25.38 | 0.2174 | 0.0029
JSR-Net,$\ell_{2}$ | 0.7337 | 23.72 | 0.253 | 0.0042
| JSR-Net, SS2 | 0.9076 | 27.31 | 0.1674 | 0.0019

###  4.3 Quantitative comparison

To quantitatively compare the performance of FBP, PD-Net, JSR model and JSR-
Net, Table 1 reports four different quality metrics of the reconstruction
results: i.e., SSIM, PSNR, MSE, and NMSE. We observe that JSR-Net trained with
SS2 loss outperforms PD-Net for sparse view CT and has comparable results with
PD-Net for limited angle CT.

##  5 Conclusion and Future work

In this work, we propose a new CNN, named JSR-Net, by unfolding the joint
spatial-Radon domain image reconstruction algorithm for incomplete data CT
imaging. A hybrid loss, containing MSE, SSIM and semantic segmentation, is
designed to train the proposed JSR-Net. Numerical results show improved
performance of JSR-Net with hybrid loss than state-of-the-art approaches.

In the future, we will try to design more reliable loss function that is more
effective in preserving details in the reconstructed image. Furthermore, the
network architecture design and applications such as interior/exterior CT are
also worth to explore with the proposed JSR-Net.

## References

  * [1] Yann LeCun, Yoshua Bengio, and Geoffrey Hinton,  “Deep learning,”  nature, vol. 521, no. 7553, pp. 436, 2015.
  * [2] Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton,  “Imagenet classification with deep convolutional neural networks,”  in NIPS, 2012, pp. 1097–1105.
  * [3] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun,  “Delving deep into rectifiers: Surpassing human-level performance on imagenet classification,”  in ICCV, 2015, pp. 1026–1034.
  * [4] Ge Wang, Jong Chu Ye, Klaus Mueller, and Jeffrey A Fessler,  “Image reconstruction is a new frontier of machine learning.,”  IEEE Trans. Med. Imaging, vol. 37, no. 6, pp. 1289–1296, 2018.
  * [5] Frank Natterer,  The mathematics of computerized tomography,  SIAM, 2001.
  * [6] Anders H Andersen,  “Algebraic reconstruction in ct from limited views,”  IEEE Trans. Med. Imaging, vol. 8, no. 1, pp. 50–55, 1989.
  * [7] Xun Jia, Bin Dong, Yifei Lou, and Steve B. Jiang,  “Gpu-based iterative cone-beam ct reconstruction using tight frame regularization.,”  Phys. Med. Biol., vol. 56, no. 13, pp. 3787–807, 2011.
  * [8] Zhiqiang Chen, Xin Jin, Liang Li, and Ge Wang,  “A limited-angle ct reconstruction method based on anisotropic tv minimization,”  Phys. Med. Biol., vol. 58, no. 7, pp. 2119, 2013.
  * [9] Yan Yang, Jian Sun, Huibin Li, and Zongben Xu,  “Deep admm-net for compressive sensing mri,”  in NIPS, 2016, pp. 10–18.
  * [10] Jonas Adler and Ozan Öktem,  “Learned primal-dual reconstruction,”  IEEE Trans. Med. Imaging, vol. 37, no. 6, pp. 1322–1332, 2018.
  * [11] Jonas Adler and Ozan Öktem,  “Solving ill-posed inverse problems using iterative deep neural networks,”  Inverse Probl., vol. 33, no. 124007 (24pp), 2017.
  * [12] Bin Dong, Jia Li, and Zuowei Shen,  “X-ray ct image reconstruction via wavelet frame based regularization and radon domain inpainting,”  J. Sci. Comput., vol. 54, no. 2, pp. 333–349, 2013.
  * [13] Martin Burger, Jahn Müller, Evangelos Papoutsellis, and Carola-Bibiane Schönlieb,  “Total variation regularization in measurement and image space for pet reconstruction,”  Inverse Probl., vol. 30, no. 10, pp. 105003, 2014.
  * [14] Ruohan Zhan and Bin Dong,  “Ct image reconstruction by spatial-radon domain data-driven tight frame regularization,”  SIAM J. Imaging Sci., vol. 9, no. 3, pp. 1063–1083, 2016.
  * [15] Haimiao Zhang, Bin Dong, and Baodong Liu,  “A reweighted joint spatial-radon domain ct image reconstruction model for metal artifact reduction,”  SIAM J. Imaging Sci., vol. 11, no. 1, pp. 707–733, 2018.
  * [16] Bin Dong and Zuowei Shen,  “MRA-Based Wavelet Frames and Applications,”  IAS Lecture Notes Series, Summer Program on “The Mathematics of Image Processing”, Park City Mathematics Institute, 2010.
  * [17] Jian-Feng Cai, Bin Dong, Stanley Osher, and Zuowei Shen,  “Image restorations: Total variation, wavelet frames, and beyond,”  J. Amer. Math. Soc., vol. 25, no. 4, pp. 1033–1089, 2012.
  * [18] Gao Huang, Zhuang Liu, Laurens Van Der Maaten, and Kilian Q Weinberger,  “Densely connected convolutional networks.,”  in CVPR, 2017, vol. 1, pp. 2261–2269.
  * [19] Yiping Lu, Aoxiao Zhong, Quanzheng Li, and Bin Dong,  “Beyond finite layer neural networks: Bridging deep architectures and numerical differential equations,”  in ICML, 2018, pp. 3276–3285.
  * [20] Diederik P. Kingma and Jimmy Lei Ba,  “Adam: A method for stochastic optimization,”  in ICLR, 2015, pp. 1–13.
  * [21] Cynthia McCollough and et al.,  “Low-dose ct for the detection and classification of metastatic liver lesions: Results of the 2016 low dose ct grand challenge,”  Med. Phys., vol. 44, no. 10, 2017.
  * [22] Xuehang Zheng, Saiprasad Ravishankar, Yong Long, and Jeffrey A Fessler,  “Pwls-ultra: An efficient clustering and learning-based approach for low-dose 3d ct image reconstruction,”  IEEE Trans. Med. Imaging, vol. 37, no. 6, pp. 1498–1510, 2018.