# A quantum algorithm for simulating non-sparse Hamiltonians

Chunhao Wang Institute for Quantum Computing and School of Computer Science,
University of Waterloo Leonard Wossnig Department of Computer Science,
University College London

###### Abstract

We present a quantum algorithm for simulating the dynamics of Hamiltonians
that are not necessarily sparse. Our algorithm is based on the assumption that
the entries of the Hamiltonian are stored in a data structure that allows for
the efficient preparation of states that encode the rows of the Hamiltonian.
We use a linear combination of quantum walks to achieve a poly-logarithmic
dependence on the precision. The time complexity measured in terms of circuit
depth of our algorithm is $O(t\sqrt{N}\lVert H\rVert\polylog(N,t\lVert
H\rVert,1/\epsilon))$, where $t$ is the evolution time, $N$ is the dimension
of the system, and $\epsilon$ is the error in the final state, which we call
precision. Our algorithm can directly be applied as a subroutine for unitary
Hamiltonians and solving linear systems, achieving a $\widetilde{O}(\sqrt{N})$
dependence for both applications.

##  1 Introduction

Hamiltonian simulation is the problem of simulating the dynamics of quantum
systems. Being the original motivation for quantum computers [16, 17], it has
been shown to be BQP-hard and is hence conjectured not to be classically
solvable in polynomial time, since such an algorithm would imply an efficient
classical algorithm for any problem with an efficient quantum algorithm,
including integer factorization [30]. Following the first proposal by Lloyd
[23] for local Hamiltonians, Aharonov and Ta-Shma gave an efficient algorithm
for sparse Hamiltonians [4]. Subsequently many algorithms have been proposed
which improved the runtime [8, 6, 9, 11, 14, 28, 32, 24, 25], mostly in the so
called sparse Hamiltonian model, and have recently culminated in a few works
with optimal dependence on all (or nearly all respectively) parameters for
sparse Hamiltonians [7, 26].

In this work, we consider a different model for obtaining information about
the Hamiltonian, which is based on the assumption that the entries of the
Hamiltonian are stored in a novel memory architecture inspired by [22]. A
quantum access to this memory model allows us to efficiently prepare states
that encode the rows of the Hamiltonian. Using the ability to prepare these
states in combination with a quantum walk [7], we give a Hamiltonian
simulation algorithm that runs in $\widetilde{O}(\sqrt{N})$ time for non-
sparse Hamiltonians of dimensionality $N$, where $\widetilde{O}(\cdot)$ hides
poly-logarithmic factors.

In [22], a quantum algorithm for recommendation systems was introduced based
on an explicit description of a memory structure, which resulted in a fast
quantum algorithm for estimating singular values for any real matrix. This
fast singular value estimation algorithm is used in a quantum algorithm for
solving dense linear systems [33]. The memory structure in [22, 33] allows us
to prepare states that corresponds to row vectors. In the problem considered
here, i.e., Hamiltonian simulation, the main hurdle to use this memory
structure is that we need to efficiently prepare states (that are different
from the states in [22, 6, 7]) from row vectors; these states in particular
encode more complicated information for each row. The precise definition of
the memory structure, i.e., the model in which we obtain our result, is
presented in Definition 1 (in Sec. 2). The memory structure appears to be
reasonable since it is a quantum extension of the classical random access
memory. A quantum computer with access to this model can efficiently prepare
quantum states that encode information about the rows of the $N\times N$
dimensional Hamiltonian, i.e., with time complexity (circuit depth) of
$\widetilde{O}(\log(N))$. This is proven in Lemma 1.

Using the efficient state preparation, we implement the linear combination of
quantum walks in order to simulate the time evolution for non-sparse
Hamiltonians with only poly-logarithmic dependency on the precision. The main
result of this work is summarized in the following theorem, which we prove in
Sec. 3.

######  Theorem 1 (Non-sparse Hamiltonian Simulation).

Let $H\in\mathbb{C}^{N\times N}$ (with $N=2^{n}$) be a Hermitian matrix stored
in the memory as specified in Definition 1. There exists a quantum algorithm
to simulate the evolution of $H$ for time $t$ and error $\epsilon$ with time
complexity (circuit depth)

$\displaystyle O\left(t\Lambda
n\log^{5/2}(t\Lambda/\epsilon)\frac{\log(t\lVert
H\rVert/\epsilon)}{\log\log(t\lVert H\rVert/\epsilon)}\right),$ (1)

where $\Lambda=\max\\{\left\|H\right\|_{1},\lVert H\rVert\\}$.

In this work, we use the notation $\left\|\cdot\right\|_{1}$ to denote the
induced 1-norm (i.e., maximum absolute column-sum norm), defined as
$\left\|H\right\|_{1}=\max_{j}\sum_{k=0}^{N-1}|H_{jk}|$; we use
$\lVert\cdot\rVert$ to denote the spectral norm, and use
$\left\|\cdot\right\|_{\mathrm{max}}$ to denote the max norm, defined as
$\left\|H\right\|_{\mathrm{max}}=\max_{i,j}|H_{jk}|$.

By the fact that $\left\|H\right\|_{1}\leq\sqrt{N}\lVert H\rVert$ (see [12],
alternatively, a more generalized version of this relation is shown in the
Appendix A), we immediately have the following corollary.

######  Corollary 1.

Let $H\in\mathbb{C}^{N\times N}$ (where $N=2^{n}$) be a Hermitian matrix
stored in the memory as specified in Definition. 1. There exists a quantum
algorithm to simulate the evolution of $H$ for time $t$ and error $\epsilon$
with time complexity (circuit depth)

$\displaystyle O\left(t\sqrt{N}\lVert H\rVert\,n\log^{5/2}(t\sqrt{N}\lVert
H\rVert/\epsilon)\frac{\log(t\lVert H\rVert/\epsilon)}{\log\log(t\lVert
H\rVert/\epsilon)}\right).$ (2)

Remarks.

  1. 1.

In Theorem 1, the circuit depth scales as $\widetilde{O}(\sqrt{N})$. However,
the gate complexity could be as large as $O(N^{2}\log^{2}(N))$ in general
because of the addressing scheme which allows for quantum access to classical
data stored in the data structure as specified in Definition 1. If some
structure of $H$ is promised (e.g., the entries of $H$ repeat in some
pattern), the addressing scheme can be implemented efficiently.

  2. 2.

If $H$ is $d$-sparse (i.e., $H$ has at most $d$ non-zero entries in any row),
the time complexity (circuit depth) of our algorithm is

$\displaystyle O\left(t\sqrt{d}\lVert H\rVert\,n\log^{5/2}(t\sqrt{d}\lVert
H\rVert/\epsilon)\frac{\log(t\lVert H\rVert/\epsilon)}{\log\log(t\lVert
H\rVert/\epsilon)}\right).$ (3)

This follows from Theorem 1 and the fact that
$\left\|H\right\|_{1}\leq\sqrt{d}\lVert H\rVert$ (as shown in Appendix A).

###  1.1 Related Work

It is important to clearly distinguish our result from other works on
Hamiltonian simulation since we use a different framework. Previous quantum
algorithms are based on access models that give constant-size local terms of a
Hamiltonian (i.e., local Hamiltonian specification), a linear combination of
tensor products of Paulis that form a Hamiltonian (i.e., linear combination of
unitaries specification), or an oracle which allows for data-access via a
black-box (i.e., sparse specification, or black-box specification). For non-
sparse Hamiltonians, a suitable model is the black-box specification: querying
the oracle with an index-pair $\ket{i,j}$ returns the corresponding entry of
the Hamiltonian $H$, i.e.,

$O_{H}\ket{i,j}\ket{z}\rightarrow\ket{i,j}\ket{z\oplus H_{i,j}},$ (4)

where $\oplus$ denotes the bit-wise XOR. With access to the black-box
Hamiltonian, simulation with error $\epsilon$ can then be provably performed
with query complexity $O(N^{3/2}\polylog(N)\log(1/\epsilon))$ for dense
Hamiltonians [6], and it was also empirically observed in [6] that for several
classes of Hamiltonians even a complexity of $O(N^{1/2}\log(N))$ queries can
be achieved. Using an efficient implementation of the query oracle, e.g. a
qRAM as presented in [18], then it is possible to achieve the time scaling
given by $O(N^{1/2}\polylog(N))$ for certain classes of Hamiltonians. Whether
this is generally possible for all Hamiltonians was left as an open problem.

We note that a direct comparison between our algorithm and previous works,
such as [6, 7], is not appropriate due to the difference in the memory models
used. The results in [6, 7] are based on the black-box oracle model, or the
sparse Hamiltonian model, while in our work, we assume that the entries of the
Hamiltonian are stored in the data structure which we explicitly give in
Definition 1. In order to perform Hamiltonian simulation we prepare here a
state of the form of Eq. (6) with time complexity (circuit depth)
$O(\polylog(N))$ which allows us to implement a quantum walk operator with
time-dependence on $\Lambda$, which is upper bound by $\sqrt{N}\lVert
H\rVert$. The mechanism is similar to the quantum walk for Hamiltonian
simulation presented in [6], where the state preparation is done using
amplitude amplification. However, in [6] the cost for preparing that state
with black-box access to the Hamiltonian incurs a cost of $O(\sqrt{N})$ in
time.

Given access to a $d$-sparse Hamiltonian oracle, the circuit depth of the
black-box Hamiltonian simulation is given by
$O(td\left\|H\right\|_{\mathrm{max}})$ as shown in [6, 7]. When $H$ is non-
sparse, their results imply the scaling
$O(tN\left\|H\right\|_{\mathrm{max}})$. In applications where
$\left\|H\right\|_{\mathrm{max}}$ is a measure of cost, our result has no
advantage against theirs, as the relation $\lVert
H\rVert\leq\sqrt{N}\left\|H\right\|_{\mathrm{max}}$ implies $O(t\sqrt{H}\lVert
H\rVert)=O(tN\left\|H\right\|_{\mathrm{max}})$. However, in the case where
$\lVert H\rVert$ is a measure of cost, as assumed e.g. in [19, 13], our result
has a quadratic improvement in the dimensionality dependence, as the relation
$\left\|H\right\|_{\mathrm{max}}\leq\lVert H\rVert$ implies
$O(tN\left\|H\right\|_{\mathrm{max}})=O(tN\lVert H\rVert)$. Here we explicitly
show that our memory model can perform the black-box Hamiltonian simulation
and hence is at least as strong as the black-box model. Whether the reverse
direction is possible or not is left open for future work.

###  1.2 Applications

One immediate application of simulating non-sparse Hamiltonians is unitary
implementation: given access to the entries of a unitary $U$, the objective is
to construct a quantum circuit to approximate this $U$ with precision
$\epsilon$. As proposed in [6, 21], unitary implementation can be reduced to
Hamiltonian simulation by considering the Hamiltonian

$\displaystyle H=\begin{pmatrix}0&U\\\ U^{{\dagger}}&0\end{pmatrix},$ (5)

and the fact that $e^{-iH\pi/2}\ket{1}\ket{\psi}=-i\ket{0}U\ket{\psi}$. If the
entries of a unitary matrix can be accessed by a black-box query oracle, it is
shown in [6] that this unitary operator can be implemented with
$O(N^{2/3}\polylog(N)\poly(1/\epsilon))$ queries to the black box. Assume the
entries of $U$ are stored in a data structure as in Definition 1. Using the
same reduction to Hamiltonian simulation, our Hamiltonian simulating algorithm
implies an implementation of $U$ with time complexity (circuit depth)
$O(\sqrt{N}\,\polylog(N,1/\epsilon))$.

Simulating non-sparse Hamiltonians can also be used as subroutines for solving
linear systems of equations for non-sparse matrices. The essence for solving a
linear system $A\vec{x}=\vec{b}$ is to apply $A^{-1}$ on $\vec{b}$, assuming
here for simplicity that $\vec{b}$ is entirely in the column-space of $A$.
When $A$ is $d$-sparse, it is shown in [13] that $A^{-1}$ can be approximated
as a linear combination of unitaries of the form $e^{-iAt}$. An efficient
quantum algorithm for Hamiltonian simulation such as [7] can then be used as a
subroutine, so that this linear system can be solved with gate complexity
$O(d\kappa^{2}\polylog(N,\kappa/\epsilon))$. If $A$ is non-sparse, their
algorithm scales as $O(N)$. Based on the memory model in [22], a quantum
algorithm for solving linear systems for non-sparse matrices was described in
[33], with time complexity (circuit depth)
$O(\kappa^{2}\sqrt{N}\polylog(N)/\epsilon)$. If we assume a similar memory
structure as in [22, 33], using our Hamiltonian simulation algorithm, together
with the _linear combinations of unitaries_ (LCU) decompositions in [13], we
have a quantum algorithm for solving linear systems for non-sparse matrices
with time complexity (circuit depth)
$O(\kappa^{2}\sqrt{N}\polylog(\kappa/\epsilon))$, which is an exponential
improvement in error dependence compared to [33].

In the following we first define our memory model, then describe the algorithm
in detail and finally prove the main results. We then finish with a summary of
this work and a discussion of the benefits and possible drawbacks of our
algorithm.

##  2 Data structure and quantum walk

We first give a precise definition of the memory structure that holds the
entries of the Hamiltonian, and then show how this memory structure can be
used to prepare states that will allow us to perform fast Hamiltonian
simulation even for dense matrices.

######  Definition 1 (Data Structure).

Let $H\in\mathbb{C}^{N\times N}$ be a Hermitian matrix (where $N=2^{n}$), and
$\Lambda=\max\\{\left\|H\right\|_{1},\lVert H\rVert\\}$. Each entry $H_{jk}$
is represented with $b$ bits of precision. Define $D$ as an array of $N$
binary trees $D_{j}$ for $j\in\\{0,\ldots,N-1\\}$. Each $D_{j}$ corresponds to
the row $H_{j}$, and its organization is specified by the following rules.

  1. 1.

The leaf node $k$ of the tree $D_{j}$ stores the tuple
$\bigl{(}{H_{jk}^{*}},\frac{\Lambda-\sigma_{j}}{N}\bigr{)}$ corresponding to
the index-entry pair $(j,k,H_{jk})$.

  2. 2.

For the level immediately above the bottom level, i.e., the leaves, the data
stored in any node is determined as follows: suppose the node has two children
storing data $(a,(\Lambda-\sigma_{j})/N)$, and $(b,(\Lambda-\sigma_{j})/N)$
respectively. (Note that $a$ and $b$ are complex numbers.) Then the tuple that
is stored in this node is given by $((|a|+|b|),2(\Lambda-\sigma_{j})/N)$.

  3. 3.

For an internal node above the parents-of-the-leaves level, we calculate the
data stored in it by the following rule: suppose the node has two children
that store tuples $(t_{00},t_{01})$, and $(t_{10},t_{11})$, respectively. Then
this internal node stores the tuple $(t_{00}+t_{10},t_{01}+t_{11})$.

Note that for each binary tree $D_{j}$ in the above data structure, the tuple
stored in an internal (non-leaf) node consists of two real numbers, while for
a leaf node, the tuple stored in it consists of one complex number and one
real number. The root node of $D_{j}$ stores the tuple
$\bigl{(}\sum_{k=0}^{N-1}|H_{jk}^{*}|,\Lambda-\sigma_{j}\bigr{)}$.

Using this memory model, we can efficiently perform the mapping described in
the following technical lemma for efficient state preparation.

######  Lemma 1 (State Preparation).

Let $H\in\mathbb{C}^{N\times N}$ be a Hermitian matrix (where $N=2^{n}$)
stored in the memory as specified in Definition 1. Each entry $H_{jk}$ is
represented with $b$ bits of precision. Then the following holds

  1. 1.

Let $\Lambda=\max\\{\left\|H\right\|_{1},\lVert H\rVert\\}$. A quantum
computer that has access to the data structure can perform the mapping

$\displaystyle\ket{j}\ket{0^{\log
N}}\ket{0}\mapsto\frac{1}{\sqrt{\Lambda}}\ket{j}\sum_{k=0}^{N-1}\ket{k}\left(\sqrt{H_{jk}^{*}}\ket{0}+\sqrt{\frac{\Lambda-\sigma_{j}}{N}}\ket{1}\right),\mbox{
for }j\in\\{0,\ldots,N-1\\},$ (6)

with time complexity (circuit depth) $O(nb^{5/2})$, where
$\sigma_{j}=\sum_{k}|H_{jk}|$, and the square-root satisfies
$\sqrt{H_{jk}}\bigl{(}\sqrt{H_{jk}^{*}}\bigr{)}^{*}=H_{jk}$.

  2. 2.

The size of the memory containing all $N^{2}$ complex entries is
$O(N^{2}\log^{2}(N))$.

In order to perform the mapping we will need the following Lemma, which
enables us to efficiently implement the conditional rotations with complex
numbers.

######  Lemma 2.

Let $\theta,\phi_{0},\phi_{1}\in\mathbb{R}$ and let
$\widetilde{\theta},\widetilde{\phi_{0}},\widetilde{\phi_{1}}$ be the $b$-bit
finite precision representation of $\theta,\phi_{0}$, and $\phi_{1}$,
respectively. Then there exists a unitary $U$ that performs the following
mapping:

$\displaystyle
U:\ket{\widetilde{\phi_{0}}}\ket{\widetilde{\phi_{1}}}\ket{\widetilde{\theta}}\ket{0}\mapsto\ket{\widetilde{\phi_{0}}}\ket{\widetilde{\phi_{1}}}\ket{\widetilde{\theta}}\left(e^{i\widetilde{\phi_{0}}}\cos(\widetilde{\theta})\ket{0}+e^{i\widetilde{\phi_{1}}}\sin(\widetilde{\theta})\ket{1}\right).$
(7)

Moreover, $U$ can be implemented with $O(b)$ 1- or 2-qubit gates.

###### Proof.

Define $U$ as

$\displaystyle
U=\left(\sum_{\widetilde{\phi_{0}}\in\\{0,1\\}^{b}}\ket{\widetilde{\phi_{0}}}\\!\bra{\widetilde{\phi_{0}}}\otimes
e^{i\ket{0}\\!\bra{0}\widetilde{\phi_{0}}}\right)\left(\sum_{\widetilde{\phi_{1}}\in\\{0,1\\}^{b}}\ket{\widetilde{\phi_{1}}}\\!\bra{\widetilde{\phi_{1}}}\otimes
e^{i\ket{1}\\!\bra{1}\widetilde{\phi_{1}}}\right)\left(\sum_{\widetilde{\theta}\in\\{0,1\\}^{b}}\ket{\widetilde{\theta}}\\!\bra{\widetilde{\theta}}\otimes
e^{-iY\widetilde{\theta}}\right),$ (8)

where $Y=\bigl{(}\begin{smallmatrix}0&-i\\\ i&0\end{smallmatrix}\bigr{)}$ is
the Pauli $Y$ matrix.

To implement
$\bigl{(}\sum_{\widetilde{\theta}\in\\{0,1\\}^{b}}\ket{\widetilde{\theta}}\\!\bra{\widetilde{\theta}}\otimes
e^{-iY\widetilde{\theta}}\bigr{)}$, we use one rotation controlled on each
qubit of the first register, with the rotation angles halved for each
successive bit. The other two factors of $U$ can be implemented in a similar
way. Therefore, $U$ can be implemented with $O(b)$ 1- or 2-qubit gates. ∎

Now we are ready to prove Lemma 1.

###### Proof of Lemma 1.

We first describe the construction and size of the memory structure. Our
structure is similar to the model presented in [22] but requires a different
connectivity.

  * •

The qRAM is built from $N$ binary trees $D_{i},i\in[N]$ and we start with an
empty tree.

  * •

When a new entry arrives $(i,j,H_{ij})$ we create or update the leaf node $j$
in the tree $D_{i}$, where the adding of the entry costs $O(\log(N))$ time,
since the depth of the tree for $H\in\mathbb{C}^{N\times N}$ is at most
$\log(N)$. Since the path from the root to the leaf is of length at most
$\log(N)$ (under the assumption that $N=2^{n}$), we have furthermore to update
at most $\log(N)$ nodes, which can be done in time $O(\log(N))$ if we store an
ordered list of the levels in the tree.

  * •

The total time of updating the tree with a new entry is given by
$\log(N)\times\log(N)=\log^{2}(N)$.

  * •

The memory requirements for $k$ entries are given by $O(k\log^{2}(N))$ as for
every entry $(j,k,H_{jk})$ at least $\log(N)$ nodes are added and each node
requires at most $O(\log(N))$ bits.

With this memory structure, we can then perform the mapping specified in Eq.
(6), with the following steps. For each $j$, we start from the root of
$D_{j}$. Starting with the initial state $\ket{j}\ket{0^{\log N}}\ket{0}$,
first apply the rotation (according to the tuple stored in the root node) on
the last register to obtain the state

$\displaystyle\frac{1}{\sqrt{\Lambda}}\ket{0}\ket{0^{\log
N}}\left(\sqrt{\sum_{k=0}^{N-1}|H_{jk}^{*}|}\ket{0}+\sqrt{\Lambda-\sigma_{j}}\ket{1}\right).$
(9)

Then a sequence of conditional rotations is applied on each qubit of the
second register to obtain the state as in Eq. (6). At each level of the binary
tree $D_{j}$, a query to the data structure is made to load the data (stored
in the tuple of each node) into a register in superpositions. Then the
rotation angles will be determined by calculating the square root and
trigonometric functions on the output of the query: this can be implemented
with $O(b^{5/2})$ 1- and 2-qubit gates based on methods using Taylor series
and long multiplication. Then the conditional rotation is applied by the
circuit described in Lemma 2, and the gate complexity for the conditional
rotation is $O(b)$. There are $n=\log(N)$ levels, so the gate complexity
excluding the implementation of the oracle is $O(nb^{5/2})$. To obtain the
quantum access to the classical data structure, a quantum addressing scheme is
required. One addressing scheme described in [18] can be used. Although the
circuit size of this addressing scheme is $\widetilde{O}(N)$ for each $D_{j}$,
its circuit depth is $O(n)$. Therefore, the time complexity (circuit depth)
for preparing the state in Eq. (6) is $O(nb^{5/2})$.

We use the following rules to determine the sign of the square-root of a
complex number: if $H_{jk}$ is not a real negative number, we write
$H_{jk}=re^{i\varphi}$ (for $r\geq 0$ and $-\pi\leq\varphi\leq\pi$) and take
$\sqrt{H_{jk}^{*}}=\sqrt{r}e^{-i\varphi/2}$; when $H_{jk}$ is a real negative
number, we write $H_{jk}=re^{i\pi}$ and take
$\sqrt{H_{jk}^{*}}=\sqrt{r}e^{-i\pi/2}$. With this recipe, we have
$\sqrt{H_{jk}}(\sqrt{H_{jk}^{*}})^{*}=H_{jk}$. ∎

An example of the state preparation procedure can be demonstrated in Fig. 1.
In this example, the initial state (omitting the first register) is
$\ket{00}\ket{0}$. Let $\sigma_{j}=|c_{0}|+|c_{1}|+|c_{2}|+|c_{3}|$. Apply the
first rotation, we obtain the state

$\displaystyle\frac{1}{\sqrt{\Lambda}}\ket{00}\left(\sqrt{|c_{0}|+|c_{1}|+|c_{2}|+|c_{3}|}\ket{0}+\sqrt{\Lambda-\sigma_{j}}\ket{1}\right)=$
(10)
$\displaystyle\frac{1}{\sqrt{\Lambda}}\left(\sqrt{|c_{0}|+|c_{1}|+|c_{2}|+|c_{3}|}\ket{00}\ket{0}+\sqrt{\Lambda-\sigma_{j}}\ket{00}\ket{1}\right).$
(11)

Then apply a rotation on the first qubit of the first register conditioned on
the last register, we obtain the state

$\displaystyle\frac{1}{\sqrt{\Lambda}}\left(\left(\sqrt{|c_{0}|+|c_{1}|}\ket{00}+\sqrt{|c_{2}|+|c_{3}|}\ket{10}\right)\ket{0}+\left(\sqrt{\frac{\Lambda-\sigma_{j}}{2}}\ket{00}+\sqrt{\frac{\Lambda-\sigma_{j}}{2}}\ket{10}\right)\ket{1}\right).$
(12)

Then apply a rotation on the second qubit of the first register conditioned on
the first qubit of the first register and last register, we obtain the state

$\displaystyle\frac{1}{\sqrt{\Lambda}}\left(\sqrt{c_{0}}\ket{00}\ket{0}+\sqrt{c_{1}}\ket{01}\ket{0}+\sqrt{c_{2}}\ket{10}\ket{0}+\sqrt{c_{3}}\ket{11}\ket{0}+\right.$
(13)
$\displaystyle\left.\sqrt{\frac{\Lambda-\sigma_{j}}{4}}\ket{00}\ket{1}+\sqrt{\frac{\Lambda-\sigma_{j}}{4}}\ket{01}\ket{1}+\sqrt{\frac{\Lambda-\sigma_{j}}{4}}\ket{10}\ket{1}+\sqrt{\frac{\Lambda-\sigma_{j}}{4}}\ket{11}\ket{1}\right)=$
(14)
$\displaystyle\frac{1}{\sqrt{\Lambda}}\left(\ket{00}\left(\sqrt{c_{0}}\ket{0}+\sqrt{\frac{\Lambda-\sigma_{j}}{4}}\ket{1}\right)+\ket{01}\left(\sqrt{c_{1}}\ket{0}+\sqrt{\frac{\Lambda-\sigma_{j}}{4}}\ket{1}\right)+\right.$
(15)
$\displaystyle\left.\ket{10}\left(\sqrt{c_{2}}\ket{0}+\sqrt{\frac{\Lambda-\sigma_{j}}{4}}\ket{1}\right)+\ket{11}\left(\sqrt{c_{3}}\ket{0}+\sqrt{\frac{\Lambda-\sigma_{j}}{4}}\ket{1}\right)\right),$
(16)

as required.

Figure 1: A example of the memory structure that allows for efficient state
preparation using a logarithm number of conditional rotations. In this
example, $\sigma_{j}=|c_{0}|+|c_{1}|+|c_{2}|+|c_{3}|$.

Based on the memory structure specified in Definition 1 and the efficient
state preparation in Lemma 1, we construct a quantum walk operator for $H$ as
follows. First define the isometry $T$ as

$\displaystyle
T=\sum_{j=0}^{N-1}\sum_{b\in\\{0,1\\}}(\ket{j}\\!\bra{j}\otimes\ket{b}\\!\bra{b})\otimes\ket{\varphi_{jb}},$
(17)

with $\ket{\varphi_{j1}}=\ket{0}\ket{1}$ and

$\displaystyle\ket{\varphi_{j0}}=\frac{1}{\sqrt{\Lambda}}\sum_{k=0}^{N-1}\ket{k}\left(\sqrt{H_{jk}^{*}}\ket{0}+\sqrt{\frac{\Lambda-\sigma_{j}}{N}}\ket{1}\right),$
(18)

where $\sigma_{j}=\sum_{k=0}^{N-1}|H_{jk}|$. Let $S$ be the swap operator that
maps $\ket{j_{0}}\ket{b_{0}}\ket{j_{1}}\ket{b_{1}}$ to
$\ket{j_{1}}\ket{b_{1}}\ket{j_{0}}\ket{b_{0}}$, for all
$j_{0},j_{1}\in\\{0,\ldots,N-1\\}$ and $b_{0},b_{1}\in\\{0,1\\}$. We observe
that

$\displaystyle\bra{j}\bra{0}T^{{\dagger}}ST\ket{k}\ket{0}=\frac{\sqrt{H_{jk}}\left(\sqrt{H_{jk}^{*}}\right)^{*}}{\Lambda}=\frac{H_{jk}}{\Lambda},$
(19)

where the second equality is ensured by the choice of the square-root as in
the proof of Lemma 1. This implies that

$\displaystyle(I\otimes\bra{0})T^{{\dagger}}ST(I\otimes\ket{0})=\frac{H}{\Lambda}.$
(20)

The quantum walk operator $U$ is defined as

$\displaystyle U=iS(2TT^{{\dagger}}-I).$ (21)

A more general characterization of the eigenvalues of quantum walks is
presented in [31]. Here we give a specific proof on the relationship between
the eigenvalues of $U$ and $H$ as follows.

######  Lemma 3.

Let the unitary operator $U$ be defined as in Eq. (21), and let $\lambda$ be
an eigenvalue of $H$ with eigenstate $\ket{\lambda}$. It holds that

$\displaystyle U\ket{\mu_{\pm}}=\mu_{\pm}\ket{\mu_{\pm}},$ (22)

where

$\displaystyle\ket{\mu_{\pm}}=$
$\displaystyle(T+i\mu_{\pm}ST)\ket{\lambda}\ket{0},\mbox{ and}$ (23)
$\displaystyle\mu_{\pm}=$ $\displaystyle\pm e^{\pm
i\arcsin(\lambda/\Lambda)}.$ (24)

###### Proof.

By the fact that $T^{{\dagger}}T=I$ and
$(I\otimes\bra{0})T^{{\dagger}}ST(I\otimes\ket{0})=H/\Lambda$, we have

$\displaystyle
U\ket{\mu_{\pm}}=\mu_{\pm}T\ket{\lambda}\ket{0}+i\left(1+\frac{2\lambda
i}{\Lambda}\mu_{\pm}\right)ST\ket{\lambda}.$ (25)

In order for this state being an eigenstate, it must hold that

$\displaystyle 1+\frac{2\lambda i}{\Lambda}\mu_{\pm}=\mu_{\pm}^{2},$ (26)

and the solution is

$\displaystyle\mu_{\pm}=\frac{\lambda
i}{\Lambda}\pm\sqrt{1-\frac{\lambda}{\Lambda}}=\pm e^{\pm
i\arcsin(\lambda/\Lambda)}.$ (27)

∎

##  3 Linear combination of unitaries and Hamiltonian simulation

To see how to convert the quantum walk operator $U$ to the Hamiltonian
simulation, we consider the linear combination of unitaries

$\displaystyle
V_{\infty}=\sum_{m=-\infty}^{\infty}\frac{J_{m}(z)}{\sum_{j={-\infty}}^{\infty}J_{j}(z)}U^{m}=\sum_{m=-\infty}^{\infty}J_{m}(z)U^{m},$
(28)

where $J_{m}(z)$ is the generating function for the Bessel function, and the
second equality follows from the fact that
$\sum_{j={\infty}}^{\infty}J_{j}(z)=1$. The eigenvalues of $V_{\infty}$ is

$\displaystyle\sum_{m=-\infty}^{\infty}J_{m}(z)\mu_{\pm}^{m}=\exp\left(\frac{z}{2}\left(\mu_{\pm}-\frac{1}{\mu_{\pm}}\right)\right)=e^{iz\lambda/\Lambda},$
(29)

where the first equality follows from (9.1.41) in [2], and the second equality
follows from Eq. (24) and the fact that $\sin(x)=(e^{ix}-e^{-ix})/2i$. Now we
consider an approximation of Eq. (28):

$\displaystyle
V_{k}=\sum_{m=-k}^{k}\frac{J_{m}(z)}{\sum_{j=-k}^{k}J_{j}(z)}U^{m}.$ (30)

Note that the eigenvalues of $U_{k}$ are

$\displaystyle\sum_{m=-k}^{k}\frac{J_{m}(z)}{\sum_{j=-k}^{k}J_{j}(z)}\mu_{\pm}^{m}.$
(31)

Note that each eigenvalue of $V_{k}$ does not depend on $\pm$ as
$J_{-m}(Z)=(-1)^{m}J_{m}(Z)$.

To bound the error in this approximation, we have the technical lemma, and the
proof is shown in Appendix B.

######  Lemma 4.

Let $V_{k}$ and $V_{\infty}$ be defined as above. There exists a positive
integer $k$ satisfying $k\geq|z|$ and

$\displaystyle k=O\left(\frac{\log(\lVert
H\rVert/\Lambda\epsilon)}{\log\log(\lVert H\rVert/\Lambda\epsilon)}\right),$
(32)

such that

$\displaystyle\lVert V_{k}-V_{\infty}\rVert\leq\epsilon,$ (33)

In the following, we provide technical lemmas for implementing linear
combination of unitaries. Suppose we are given the implementations of
unitaries $U_{0}$, $U_{1}$, …, $U_{m-1}$, and coefficients
$\alpha_{0},\alpha_{1},\ldots,\alpha_{m}$. Then the unitary

$\displaystyle V=\sum_{j=0}^{m-1}\alpha_{j}U_{j}$ (34)

can be implemented probabilistically by the technique called _linear
combination of unitaries_ (LCU) [29]. Provided
$\sum_{j=0}^{m-1}|\alpha_{j}|\leq 2$, $V$ can be implemented with success
probability $1/4$. To achieve this, we define the $\mbox{multiplexed-}U$
operation, which is denoted by $\mbox{multi-}U$, as

$\displaystyle\mbox{multi-}U\ket{j}\ket{\psi}=\ket{j}U_{j}\ket{\psi}.$ (35)

We summarize the probabilistic implementation of $V$ in the following lemma,
whose proof is shown in Appendix B.

######  Lemma 5.

Let $\mbox{multi-}U$ be defined as above. If $\sum_{j=0}^{m-1}|\alpha_{j}|\leq
2$, then there exists a quantum circuit that maps $\ket{0}\ket{0}\ket{\psi}$
to the state

$\displaystyle\frac{1}{2}\ket{0}\ket{0}\left(\sum_{j=0}^{m-1}U_{j}\ket{\psi}\right)+\frac{\sqrt{3}}{2}\ket{\Phi^{\bot}},$
(36)

where $(\ket{0}\\!\bra{0}\otimes\ket{0}\\!\bra{0}\otimes
I)\ket{\Phi^{\bot}}=0$. Moreover, this quantum circuit uses $O(1)$
applications of $\mbox{multi-}U$ and $O(m)$ 1- and 2-qubit gates.

Let $W$ be the quantum circuit in Lemma 5, and let $P$ be the projector
defined as $P=\ket{0}\\!\bra{0}\otimes\ket{0}\\!\bra{0}\otimes I$. We have

$\displaystyle
PW\ket{0}\ket{0}\ket{\psi}=\frac{1}{2}\ket{0}\ket{0}\sum_{j=0}^{m-1}U_{j}\ket{\psi}.$
(37)

If $\sum_{j=0}^{m-1}U_{j}$ is a unitary operator, one application of the
oblivious amplitude amplification operator $-W(1-2P)W^{{\dagger}}(1-2P)W$
implements $\sum_{j=0}^{m-1}U_{j}$ with certainty. However, in our
application, the unitary operator $\widetilde{W}$ implements an approximation
of $V_{\infty}$ in the sense that

$\displaystyle
P\widetilde{W}\ket{0}\ket{0}\ket{\psi}=\frac{1}{2}\ket{0}\ket{0}V_{k}\ket{\psi},$
(38)

with $\lVert V_{k}-V_{\infty}\rVert\leq\epsilon$. The following lemma shows
that the error caused by the oblivious amplitude amplification is bounded by
$O(\epsilon)$, and the proof is given in Appendix B.

######  Lemma 6.

Let the projector $P$ be defined as above. If a unitary operator
$\widetilde{W}$ satisfies
$P\widetilde{W}\ket{0}\ket{0}\ket{\psi}=\frac{1}{2}\ket{0}\ket{0}\widetilde{V}\ket{\psi}$
where $\lVert\widetilde{V}-V\rVert\leq\epsilon$. Then
$\lVert-\widetilde{W}(I-2P)\widetilde{W}^{{\dagger}}(I-2P)\widetilde{W}\ket{0}\ket{0}\ket{\psi}-\ket{0}\ket{0}V\ket{\psi}\rVert=O(\epsilon)$.

Now we are ready to prove Theorem 1

###### Proof of Theorem 1.

The proof we outline here follows closely the proof given in [7]. The
intuition of this algorithm is to divide the simulation into $O(t\Lambda)$
segments, with each segment simulating $e^{-iH/2}$. To implement each segment,
we use the LCU technique to implement $V_{k}$ defined in Eq. (30), with
coefficients $\alpha_{m}=J_{m}(z)/\sum_{j=-k}^{k}J_{j}(z)$. When $z=-1/2$,
$\sum_{j=-k}^{k}|\alpha_{j}|<2$ and Lemmas 5 and 6 can be applied. By Eq. (29)
and using Lemma 4, set

$\displaystyle k=O\left(\frac{\log(\lVert
H\rVert/\Lambda\epsilon^{\prime})}{\log\log(\lVert
H\rVert/\Lambda\epsilon^{\prime})}\right),$ (39)

and we obtain a segment that simulates $e^{-iH/2\Lambda}$ with error bounded
by $O(\epsilon^{\prime})$. Repeat the segment $O(t\Lambda)$ times with error
$\epsilon^{\prime}=\epsilon/t\Lambda$, and we obtain a simulation of
$e^{-iHt}$ with error bounded by $\epsilon$. We have

$\displaystyle k=O\left(\frac{\log(t\lVert H\rVert/\epsilon)}{\log\log(t\lVert
H\rVert/\epsilon)}\right).$ (40)

By Lemma 5, each segment can be implemented by $O(1)$ application of
$\mbox{multi-}U$ and $O(k)$ 1- and 2-qubit gates, as well as the cost for
computing the coefficients $\alpha_{m}$ for $m\in\\{-k,\ldots,k\\}$. The
coefficients can be classically computed using the methods in [10, 27], and
the cost is $O(k)$ times the number of bits of precision, which is
$O(\log(t\lVert H\rVert/\epsilon)$. The cost for each $\mbox{multi-}U$ is $k$
times the cost for implementing the quantum walk $U$. By Lemma 1, the state in
Eq. (18) can be prepared with time complexity (circuit depth) $O(nb^{5/2})$,
where $b$ is the number of bit of precision. To achieve the overall error
bound $\epsilon$, we choose $b=O(\log(t\Lambda/\epsilon))$. Hence the time
complexity for the state preparation is $O(n\log^{5/2}(t\Lambda/\epsilon))$,
which is also the time complexity for applying the quantum walk $U$.
Therefore, the time complexity for one segment is

$\displaystyle O\left(n\log^{5/2}(t\Lambda/\epsilon)\frac{\log(t\lVert
H\rVert/\epsilon)}{\log\log(t\lVert H\rVert/\epsilon)}\right).$ (41)

Considering $O(t\Lambda)$ segments, the time complexity is as claimed. ∎

##  4 Discussion

We presented a quantum Hamiltonian simulation algorithm that provably runs in
$\widetilde{O}(\sqrt{N})$ time for non-sparse Hamiltonians, which gives a
polynomial speedup compared to previous results. In order to do so, we
required access to a seemingly more powerful memory model for state
preparation that can prepare states that encode the columns and rows of the
Hamiltonian in logarithmic time. Our technique for Hamiltonian simulation
combines ideas from linear combination of quantum walks and an efficient
memory model which prepares a special states to provide the improved
performance.

We note that the used memory structure may require a large overhead in
practice, if the data is not already stored in it. Furthermore it might be
hard to implement such a memory structure physically due to the exponential
amount of quantum resources[1, 3, 15]. However, if the Hamiltonian is highly
structured (i.e., the entries repeat in some pattern), the memory model can be
efficiently implemented. Another potential point of criticism is the required
error rate of such a device, as some computations will require an error rate
per gate of $\epsilon_{err}=O(1/\mathrm{poly}(N))$ to retain a feasible error
rate for applications [5]. Whereas, not all computations might need such low
error rates [5] and hence the feasibility of our algorithm as a subroutine in
an explicit algorithm must be further validated experimentally. To summarize,
our algorithm inherits many problems of Grover’s search algorithm for
unordered classical data.

However, if the above mentioned caveats can be overcome, our algorithm still
supplies a polynomial speedup over known quantum algorithms for Hamiltonian
simulation and furthermore allows for a polynomial speedup in comparison with
the best known classical algorithms for several practical problems, as for
example linear regression or linear systems.

##  5 Acknowledgement

We thank Richard Cleve and Simone Severini for the discussion and comments on
this project. CW acknowledges financial support by a David R. Cheriton
Graduate Scholarship. LW acknowledges financial support by the Royal Society
and Cambridge Quantum Computing Ltd.

## References

  * [1] Scott Aaronson  “Read the fine print”  In _Nature Physics_ 11.4 Nature Publishing Group, 2015, pp. 291
  * [2] Milton Abramowitz and Irene A Stegun  “Handbook of mathematical functions: with formulas, graphs, and mathematical tables”  Courier Corporation, 1964
  * [3] Jeremy Adcock et al.  “Advances in quantum machine learning”  In _arXiv preprint arXiv:1512.02900_ , 2015
  * [4] Dorit Aharonov and Amnon Ta-Shma  “Adiabatic quantum state generation and statistical zero knowledge”  In _Proceedings of the thirty-fifth annual ACM symposium on Theory of computing_ , 2003, pp. 20–29  ACM
  * [5] Srinivasan Arunachalam et al.  “On the robustness of bucket brigade quantum RAM”  In _New Journal of Physics_ 17.12 IOP Publishing, 2015, pp. 123010
  * [6] Dominic W Berry and Andrew M Childs  “Black-box Hamiltonian simulation and unitary implementation”  In _arXiv preprint arXiv:0910.4157_ , 2009
  * [7] Dominic W Berry, Andrew M Childs and Robin Kothari  “Hamiltonian simulation with nearly optimal dependence on all parameters”  In _Foundations of Computer Science (FOCS), 2015 IEEE 56th Annual Symposium on_ , 2015, pp. 792–809  IEEE
  * [8] Dominic W Berry, Graeme Ahokas, Richard Cleve and Barry C Sanders  “Efficient quantum algorithms for simulating sparse Hamiltonians”  In _Communications in Mathematical Physics_ 270.2 Springer, 2007, pp. 359–371
  * [9] Dominic W Berry et al.  “Exponential improvement in precision for simulating sparse Hamiltonians”  In _Forum of Mathematics, Sigma_ 5, 2017  Cambridge University Press
  * [10] W.. Bickley et al.  “Bessel Functions: Part II. Functions of Positive Integer Order”  Cambridge University Press, 1960
  * [11] Andrew M Childs  “On the relationship between continuous-and discrete-time quantum walk”  In _Communications in Mathematical Physics_ 294.2 Springer, 2010, pp. 581–603
  * [12] Andrew M. Childs and Robin Kothari  “Limitations on the simulation of non-sparse Hamiltonians”  In _Quantum Information & Computation_ 10.7&8, 2010, pp. 669–684
  * [13] Andrew M Childs, Robin Kothari and Rolando D Somma  “Quantum algorithm for systems of linear equations with exponentially improved dependence on precision”  In _SIAM Journal on Computing_ 46.6 SIAM, 2017, pp. 1920–1950
  * [14] Andrew M Childs and Nathan Wiebe  “Hamiltonian simulation using linear combinations of unitary operations”  In _arXiv preprint arXiv:1202.5822_ , 2012
  * [15] Carlo Ciliberto et al.  “Quantum machine learning: a classical perspective”  In _Proc. R. Soc. A_ 474.2209, 2018, pp. 20170551  The Royal Society
  * [16] Richard P Feynman  “Simulating physics with computers”  In _International journal of theoretical physics_ 21.6-7 Springer, 1982, pp. 467–488
  * [17] Richard P Feynman  “Quantum mechanical computers”  In _Foundations of physics_ 16.6 Springer, 1986, pp. 507–531
  * [18] Vittorio Giovannetti, Seth Lloyd and Lorenzo Maccone  “Quantum random access memory”  In _Physical review letters_ 100.16 APS, 2008, pp. 160501
  * [19] Aram W Harrow, Avinatan Hassidim and Seth Lloyd  “Quantum algorithm for linear systems of equations”  In _Physical review letters_ 103.15 APS, 2009, pp. 150502
  * [20] Roger A Horn and Charles R Johnson  “Matrix analysis”  Cambridge university press, 1990
  * [21] Stephen P Jordan and Pawel Wocjan  “Efficient quantum circuits for arbitrary sparse unitaries”  In _Physical Review A_ 80.6 APS, 2009, pp. 062301
  * [22] Iordanis Kerenidis and Anupam Prakash  “Quantum recommendation systems”  In _arXiv preprint arXiv:1603.08675_ , 2016
  * [23] Seth Lloyd  “Universal quantum simulators”  In _Science_ JSTOR, 1996, pp. 1073–1078
  * [24] Guang Hao Low and Isaac L Chuang  “Hamiltonian simulation by qubitization”  In _arXiv preprint arXiv:1610.06546_ , 2016
  * [25] Guang Hao Low and Isaac L Chuang  “Hamiltonian Simulation by Uniform Spectral Amplification”  In _arXiv preprint arXiv:1707.05391_ , 2017
  * [26] Guang Hao Low and Isaac L Chuang  “Optimal Hamiltonian simulation by quantum signal processing”  In _Physical review letters_ 118.1 APS, 2017, pp. 010501
  * [27] F… Olver  “Error analysis of Miller’s recurrence algorithm”  In _Mathematics of Computation_ 18.85, 1964, pp. 65–74
  * [28] David Poulin, Angie Qarry, Rolando Somma and Frank Verstraete  “Quantum simulation of time-dependent Hamiltonians and the convenient illusion of Hilbert space”  In _Physical review letters_ 106.17 APS, 2011, pp. 170501
  * [29] R.  “Efficient algorithms in quantum query complexity” PhD Thesis, University of Waterloo, 2014
  * [30] Peter W Shor  “Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer”  In _SIAM review_ 41.2 SIAM, 1999, pp. 303–332
  * [31] Mario Szegedy  “Quantum speed-up of Markov chain based algorithms”  In _Foundations of Computer Science, 2004. Proceedings. 45th Annual IEEE Symposium on_ , 2004, pp. 32–41  IEEE
  * [32] Nathan Wiebe, Dominic W Berry, Peter Hï¿½yer and Barry C Sanders  “Simulating quantum dynamics on a quantum computer”  In _Journal of Physics A: Mathematical and Theoretical_ 44.44 IOP Publishing, 2011, pp. 445308
  * [33] Leonard Wossnig, Zhikuan Zhao and Anupam Prakash  “Quantum Linear System Algorithm for Dense Matrices”  In _Physical review letters_ 120.5 APS, 2018, pp. 050502

##  Appendix A The relation between $\left\|H\right\|_{1}$ and $\lVert
H\rVert$

We prove the following proposition.

######  Proposition 2.

If $A\in\mathbb{C}^{N\times N}$ has at most $d$ non-zero entries in any row,
it holds that $\left\|A\right\|_{1}\leq\sqrt{d}\lVert A\rVert$.

###### Proof.

First observe that $\lVert
A\rVert^{2}=\rho(A^{*}A)\leq\sum_{i}\lambda_{i}(A^{*}A)=\mathrm{Tr}\left(A^{*}A\right)=\lVert
A\rVert^{2}_{F}$, and furthermore we have that $\sum_{ij}|a_{ij}|^{2}\leq
d\max_{j\in[N]}\sum_{i}|a_{ij}|^{2}$. From this we have that $\lVert
A\rVert\leq\sqrt{d}\lVert A\rVert_{1}$ for a $d$-sparse $A$. By Theorem
5.6.18, in [20], we have that $\lVert A\rVert_{1}\leq C_{M}(1,*)\lVert
A\rVert$ for $C_{M}(1,*)=\max_{A\neq 0}\frac{\lVert A\rVert}{\lVert
A\rVert_{1}}$ and using the above we have $C_{M}(mc,*)\leq\sqrt{d}$. Therefore
we find that $\lVert A\rVert_{mc}\leq\sqrt{d}\lVert A\rVert$ as desired. ∎

It immediately follows that $\left\|A\right\|_{1}\leq\sqrt{N}\lVert A\rVert$
for dense $A$, i.e., $s=N$.

##  Appendix B Proofs of technical lemmas

###### Proof of Lemma 6.

The proof outlined here follows closely the proof of Lemma 8 in [7]. Recalling
the definition of $V_{k}$ and $V_{\infty}$, we define the weights in $V_{k}$
by

$\alpha_{m}:=\frac{J_{m}(z)}{C_{k}},$

where $C_{k}=\sum_{l=-k}^{k}J_{l}(z)$. The normalization here is chosen such
that $\sum_{m}a_{m}=1$ which will give the best result [7].
Since

$\sum_{m=-\infty}^{\infty}J_{m}(z)=\sum_{m=[-\infty:-k-1;k+1:\infty]}J_{m}(z)+\sum_{m=-k}^{k}J_{m}(z)=1,$
(42)

observe that we have two error sources. The first one coming from the
truncation of the series, and the second one from the different
renormalization of the terms which introduces an error in the first $\lvert
m\rvert\leq k$ terms in the sum. We therefore start by bounding the
normalization factor $C_{k}$. For the Bessel-functions for all $m$ it holds
that $\lvert J_{m}(z)\rvert\leq\frac{1}{\lvert
m\lvert!}\left\lvert\frac{z}{2}\right\rvert^{|m|}$, since
$J_{-m}(z)=(-1)^{m}J_{m}(z)$ [2][9.1.5]. For $\lvert m\rvert\leq k$ we can
hence find the following bound on the truncated part

$\displaystyle\sum_{m=[-\infty:-k-1;k+1:\infty]}\lvert J_{m}(z)\rvert$
$\displaystyle=2\sum_{m=k+1}^{\infty}\lvert J_{m}(z)\rvert\leq
2\sum_{m=k+1}^{\infty}\frac{\lvert z/2\rvert^{m}}{m!}$
$\displaystyle<2\frac{\lvert
z/2\rvert^{k+1}}{(k+1)!}\sum_{m=k+1}^{\infty}\left(\frac{1}{2}\right)^{m-(k+1)}=\frac{4\lvert
z/2\rvert^{k+1}}{(k+1)!}.$ (43)

Since $\sum_{m}J_{m}(z)=1$, based on the normalization, we hence find that

$\sum_{m=-k}^{k}J_{m}(z)\geq\left(1-\frac{4\lvert
z/2\rvert^{k+1}}{(k+1)!}\right),$ (44)

which is an lower bound on the normalization factor $C_{k}$. Since
$a_{m}=\frac{J_{m}(z)}{C_{k}}$, we and since the correction is small, this
implies that

$a_{m}=J_{m}(z)\left(1+O\left(\frac{\lvert
z/2\rvert^{k+1}}{(k+1)!}\right)\right),$ (45)

and we have a multiplicative error based on the renormalization.

Next we want to bound the error in the truncation before we join the two error
sources.
From Eq. (29) we know that

$e^{iz\lambda/\Lambda}-1=\sum_{m=-\infty}^{\infty}J_{m}(z)(\mu_{\pm}^{m}-1),$
(46)

by the normalization of $\sum_{m}J_{m}(z)$. From this we can see that we can
hence obtain a bound on the truncated $J_{m}(z)$ as follows.

$\sum_{m=-k}^{k}J_{m}(z)(\mu_{\pm}^{m}-1)=e^{iz\lambda/\Lambda}-1-\sum_{\begin{subarray}{c}m=[-\infty:-(k+1);\\\
(k+1):\infty]\end{subarray}}J_{m}(z)(\mu_{\pm}^{m}-1).$ (47)

Therefore we can upper bound the left-hand side in terms of the exact value of
$V_{\infty}$, i.e. $e^{iz\lambda/\Lambda}$ if we can bound the right-most term
in Eq. (47). Using furthermore the bound in Eq. (45) we obtain

$\sum_{m=-k}^{k}a_{m}(z)(\mu_{\pm}^{m}-1)=\left(e^{iz\lambda/\Lambda}-1-\sum_{\begin{subarray}{c}m=[-\infty:-(k+1);\\\
(k+1):\infty]\end{subarray}}J_{m}(z)(\mu_{\pm}^{m}-1)\right)\left(1+O\left(\frac{\lvert
z/2\rvert^{k+1}}{(k+1)!}\right)\right),$ (48)

which reduced with $\left\lvert 2^{iz\lambda/\Lambda}-1\right\rvert\leq\lvert
z\lambda/\Lambda\rvert$ and $|z|\leq k$ to

$\sum_{m=-k}^{k}a_{m}(z)(\mu_{\pm}^{m}-1)=e^{iz\lambda/\Lambda}-1-O\left(\sum_{\begin{subarray}{c}m=[-\infty:-(k+1);\\\
(k+1):\infty]\end{subarray}}J_{m}(z)(\mu_{\pm}^{m}-1)\right).$ (49)

We can then obtain the desired bound $\lVert V_{\infty}-V_{k}\rVert$ by
reordering the above equation, and using that $\sum_{m=-k}^{k}a_{m}(z)=1$ such
that we have

$\lVert
V_{\infty}-V_{k}\rVert=\left\lvert\sum_{m=-k}^{k}a_{m}\mu_{\pm}^{m}-e^{iz\lambda/\Lambda}\right\rvert=O\left(\sum_{\begin{subarray}{c}m=[-\infty:-(k+1);\\\
(k+1):\infty]\end{subarray}}J_{m}(z)(\mu_{\pm}^{m}-1)\right).$ (50)

We hence only need to bound the right-hand side.
For $\mu_{+}$ we can use that $\lvert\mu_{+}^{m}-1\rvert\leq 2\lvert
m\lambda/\Lambda\rvert=:2\lvert m\nu\rvert$ and obtain the bound
$2\frac{\lvert\nu\rvert}{k!}\left\lvert\frac{z}{2}\right\rvert^{k+1}$ [7]. For
the $\mu_{-}$ case we need to refine the analysis and will show that the bound
remains the same. Let $\nu:=\lambda/\Lambda$ as above. First observe that
$J_{m}(z)\mu_{-}^{m}+J_{-m}(z)\mu_{i}^{-m}=J_{-m}(z)\mu_{+}^{-m}+J_{m}(z)\mu_{+}^{m}$,
and it follows that

$\displaystyle\sum_{m=-\infty}^{-(k+1)}J_{m}(z)(\mu_{-}^{m}-1)$
$\displaystyle+\sum_{m=k+1}^{\infty}J_{m}(z)(\mu_{-}^{m}-1)$
$\displaystyle=\sum_{m=-\infty}^{-(k+1)}J_{m}(z)(\mu_{+}^{m}-1)+\sum_{m=k+1}^{\infty}J_{m}(z)(\mu_{+}^{m}-1).$
(51)

Therefore we only need to treat the $\mu_{+}$ case.

$\displaystyle\left\lvert\sum_{\begin{subarray}{c}m=[-\infty:-(k+1);\\\
(k+1):\infty]\end{subarray}}J_{m}(z)(\mu_{+}^{m}-1)\right\rvert$
$\displaystyle\leq
2\left\lvert\sum_{m=k+1}^{\infty}J_{m}(z)(\mu_{+}^{m}-1)\right\rvert$
$\displaystyle\leq 2\sum_{m=k+1}^{\infty}\lvert
J_{m}(z)\rvert\lvert\mu_{+}^{m}-1\rvert$
$\displaystyle=2\sum_{m=k+1}^{\infty}\frac{1}{|m|!}\left\lvert\frac{z}{2}\right\rvert^{|m|}\lvert\mu_{+}^{m}-1\rvert$
$\displaystyle\leq
4\sum_{m=k+1}^{\infty}\frac{1}{|m|!}\left\lvert\frac{z}{2}\right\rvert^{|m|}m\lvert\nu\rvert$
$\displaystyle<\frac{8\lvert\nu\rvert}{(k+1)!}\left\lvert\frac{z}{2}\right\rvert^{k+1}(k+2).$
(52)

Using this bound we hence obtain from Eq. (50),

$\lVert
V_{\infty}-V_{k}\rVert=\left\lvert\sum_{m=-k}^{k}a_{m}\mu_{\pm}^{m}-e^{iz\lambda/\Lambda}\right\rvert\leq
O\left(\frac{\lambda}{k!\
\Lambda}\left\lvert\frac{z}{2}\right\rvert^{k+1}\right)=O\left(\frac{\lVert
H\rVert(z/2)^{k+1}}{\Lambda k!}\right).$ (53)

In order for the above equation being upper-bounded by $\epsilon$, it suffices
to some $k$ that is upper bounded as claimed. ∎

###### Proof of Lemma 5.

Let $s=\sum_{j=0}^{m-1}|\alpha_{j}|$. We first define the unitary operator $B$
to prepare the coefficients:

$\displaystyle
B\ket{0}\ket{0}=\left(\sqrt{\frac{s}{2}}\ket{0}+\sqrt{1-\frac{s}{2}}\ket{1}\right)\otimes\frac{1}{\sqrt{s}}\sum_{j=0}^{m-1}\sqrt{\alpha_{j}}\ket{j}.$
(54)

Define the unitary operator $W$ as $W=(B^{{\dagger}}\otimes
I)(I\otimes\mbox{multi-}U)(B\otimes I)$. We claim that $W$ performs the
desired mapping, as

$\displaystyle W\ket{0}\ket{0}\ket{\psi}=$ $\displaystyle(B^{{\dagger}}\otimes
I)(I\otimes\mbox{multi-}U)(B\otimes I)\ket{0}\ket{0}\ket{\psi}$ (55)
$\displaystyle=$ $\displaystyle\frac{1}{\sqrt{2}}(B^{{\dagger}}\otimes
I)\ket{0}\sum_{j=0}^{m-1}\sqrt{\alpha_{j}}\ket{j}U_{j}\ket{\psi}+\sqrt{\frac{2-s}{2s}}(B^{{\dagger}}\otimes
I)\ket{1}\sum_{j=0}^{m-1}\sqrt{\alpha_{j}}\ket{j}U_{j}\ket{\psi}$ (56)
$\displaystyle=$
$\displaystyle\frac{1}{2}\ket{0}\ket{0}\sum_{j=0}^{m-1}\alpha_{j}U_{j}\ket{\psi}+\sqrt{\gamma}\ket{\Phi^{\bot}},$
(57)

where $\ket{\Phi^{\bot}}$ is a state satisfying
$(\ket{0}\\!\bra{0}\otimes\ket{0}\\!\bra{0}\otimes I)\ket{\Phi^{\bot}}=0$, and
$\gamma$ is some normalization factor.

The number of applications of $\mbox{multi-}U$ is constant, as in the
definition of $W$. To implement the unitary operator $B$, $O(m)$ 1- and
2-qubit gates suffice. ∎

###### Proof of Lemma 6.

We have

$\displaystyle-\widetilde{W}(I-2P)\widetilde{W}^{{\dagger}}(I-2P)\widetilde{W}\ket{0}\ket{0}\ket{\psi}=$
$\displaystyle(\widetilde{W}+2P\widetilde{W}-4\widetilde{W}P\widetilde{W}^{{\dagger}}PW)\ket{0}\ket{0}\ket{\psi}$
(58) $\displaystyle=$
$\displaystyle(\widetilde{W}+2P\widetilde{W}-4\widetilde{W}P\widetilde{W}^{{\dagger}}PPWP)\ket{0}\ket{0}\ket{\psi}$
(59) $\displaystyle=$
$\displaystyle\widetilde{W}\ket{0}\ket{0}\ket{\psi}+\ket{0}\ket{0}\widetilde{V}\ket{\psi}-\widetilde{W}\left(\ket{0}\ket{0}\widetilde{V}^{{\dagger}}\widetilde{V}\ket{\psi}\right)$
(60)

Because $\lVert\widetilde{V}-V\rVert\leq\epsilon$ and $V$ is a unitary
operator, we have
$\lVert\widetilde{V}^{{\dagger}}\widetilde{V}-I\rVert=O(\epsilon)$. Therefore,
we have

$\displaystyle\lVert-\widetilde{W}(I-2P)\widetilde{W}^{{\dagger}}(I-2P)\widetilde{W}\ket{0}\ket{0}\ket{\psi}-\ket{0}\ket{0}\widetilde{V}\ket{\psi}\rVert=O(\epsilon).$
(61)

Thus

$\displaystyle\lVert-\widetilde{W}(I-2P)\widetilde{W}^{{\dagger}}(I-2P)\widetilde{W}\ket{0}\ket{0}\ket{\psi}-\ket{0}\ket{0}V\ket{\psi}\rVert=O(\epsilon).$
(62)

∎


Online dynamic flat-field correction for MHz Microscopy data at European XFEL

Sarlota Birnsteinovaa,*, Danilo E. Ferreira de Limaa, Egor Soboleva, Henry J.
Kirkwooda, Valerio Belluccia, Richard J. Beana, Chan Kima, Jayanath C. P.
Koliyadua, Tokushi Satoa, Fabio Dall’Antoniaa, Eleni Myrto Asimakopouloub,
Zisheng Yaob, Khachiwan Buakora,b, Yuhe Zhangb, Alke Meentsc, Henry N.
Chapmanc,d, Adrian P. Mancusoa,e,f, Pablo Villanueva-Perezb, Patrik Vagovičc,a

a European XFEL GmbH, Schenefeld Germany

b Synchrotron Radiation Research and NanoLund, Lund University Sweden

c Center for Free-Electron Laser Science (CFEL), DESY, Hamburg Germany

d University of Hamburg, Hamburg, Germany

e Diamond Light Source Ltd., Harwell Science and Innovation Campus, Didcot,
OX11 0DE United Kingdom

f Department of Chemistry and Physics, La Trobe Institute for Molecular
Science, La Trobe University, Melbourne, Victoria Australia

*<EMAIL_ADDRESS>
## Abstract

The X-ray microscopy technique at the European X-ray free-electron laser
(EuXFEL), operating at a MHz repetition rate, provides superior contrast and
spatial-temporal resolution compared to typical microscopy techniques at other
X-ray sources. In both online visualization and offline data analysis for
microscopy experiments, baseline normalization is essential for further
processing steps such as phase retrieval and modal decomposition. In addition,
access to normalized projections during data acquisition can play an important
role in decision-making and improve the quality of the data. However, the
stochastic nature of XFEL sources hinders the use of existing flat-flied
normalization methods during MHz X-ray microscopy experiments. Here, we
present an online dynamic flat-field correction method based on principal
component analysis of dynamically evolving flat-field images. The method is
used for the normalization of individual X-ray projections and has been
implemented as an online analysis tool at the Single Particles, Clusters, and
Biomolecules and Serial Femtosecond Crystallography (SPB/SFX) instrument of
EuXFEL.

##  1 Introduction

In the past few decades, X-ray microscopy was established as an important tool
for the study of a large range of applications, from biology to material
research, at third and fourth-generation X-ray synchrotron sources, where the
achievable spatial resolution and sensitivity have been pushed to the limits
[1]. Recently, fast megahertz (MHz)-rate X-ray microscopy synchronized to
individual pulses has been pioneered at the Advanced Photon Source (APS) [2]
and at the ESRF - the European Synchrotron [3, 4], which made it possible to
image fast stochastic phenomena. MHz X-ray microscopy experiment with frame
acquisition synchronized to individual X-ray pulses has been demonstrated
recently at European XFEL (EuXFEL) [5]. This development was enabled by the
unique properties of X-ray Free-Electron Lasers (XFELs). In particular, EuXFEL
is currently the only XFEL source providing intense pulses at MHz rates [6],
which opens up the potential for novel methods and applications. EuXFEL
provides two to four orders of magnitude more photons per pulse ($\sim
10^{12}-10^{13}$) than the most brilliant synchrotron, i.e ESRF. EuXFEL also
provides high spacial coherence, with a beam path of 1 km and a maximum of 4
$\mu$rad divergence, with photon energy up to 30 keV. The unique properties of
EuXFEL come with unique issues, which complicate the acquisition of high-
quality MHz microscopic data. XFEL’s X-ray pulses are generated through the
SASE process [7], which leads to strongly stochastic spatial and temporal
fluctuations. Another type of aberration appearing in microscopy data is often
referred to as fixed-pattern noise. It originates from the different responses
of the detector’s pixels and imperfections in the optic systems, such as
scintillator screens, lenses, etc. Ideally, image aberrations should be
corrected, so that only the signal variations originating from the X-rays’
interaction with a sample are observed. Methods aiming to achieve this goal
are often referred to as “flat-field correction” methods.

Flat-field correction methods rely on using reference images averaged over
many realizations of acquired images. Reference images acquired with beam
illumination are referred to as “flat-field images”, while images acquired
without beam illumination are referred to as “dark-field images”. However, as
a stationary method, traditional flat-field correction, which assumes a flat-
filed that is stable with time, is unable to sufficiently normalize stochastic
effects, such as those stemming from SASE processes which cause variations in
beam filed from XFEL pulse to pulse.

There are several methods used for correcting both these effects [8, 9, 10,
11]. The issue of normalization of a dynamically changing flat-field has been
tackled by a method [8], in which the reference flat-fields are represented in
a latent lower-dimensional vector space and an “effective” flat-field is
chosen per input image in the latent space. This method, called dynamic flat-
field correction, was introduced in Ref. [8] for synchrotron data. The
application of a similar method was demonstrated on EuXFEL’s full-field
imaging data [9] of high magnification and in combination with deep-learning
approaches [10].

In this paper, we present an _online_ method used at EuXFEL to perform dynamic
flat-field normalization for MHz microscopy data. EuXFEL users often require
fast feedback when performing an experiment, such that they may adapt their
configurations to improve the quality of the data taken. In such conditions,
obtaining corrected images during the experiment may be the deciding factor
between a successful experiment, or a failure to obtain reliable data. The
online flat-field correction method described in this paper is optimized to
provide output to the user as data arrives, minimising the delay between data-
acquisition and analysis. The method proposed here also provides an indication
of whether the flat-field images collected as references are sufficient, or
whether a new dataset must be collected due to changes in the experimental
conditions. This helps to enhance the method’s performance, as a dataset of
flat-fields that sufficiently describes the illumination at a given moment of
measurement has a crucial impact on the resulting quality of the normalization
when affected by the fluctuations of SASE illumination. Moreover, the
normalized data provided during experiments can be further re-used as input to
other processing tools, for which pre-processed data are essential. Our main
objective is to ensure high quality image data at high speeds. We are using a
software tool provided by EuXFEL, enabling us to use online data and follow
with the normalization procedure in near real-time. The resulting online
normalization tool for MHz microscopy was developed and implemented at the
Single Particles, Clusters, and Biomolecules and Serial Femtosecond
Crystallography (SPB/SFX) instrument of EuXFEL. Our online implementation
based on the dynamic method led to a major improvement in the quality of
corrected images compared to results using conventional flat-field correction
approach.

The manuscript is organised as follows: in Section 2, we describe flat-field
correction methods. The normalization algorithm and implementation are
described in Section 3. We discuss the results and the performance of
implementation using test data in Section 4, followed by conclusion and
discussion that are given in Section 5.

##  2 Flat-field correction

In this section, we briefly review of the conventional flat-field correction
method. We then expand on the dynamic flat-field correction method, which is
implemented in this work. While we closely follow [8, 10], differences in the
methods are pointed out as needed.

A stationary, conventional flat-field correction removes the effects of
spatially uneven illumination functions. Such effects may be caused by
variations in the detector response in different pixels or modules, or effects
arising from the scintillator and optics configurations. One assumes that
there has been a previous data acquisition of $N_{f}$ flat-field images
$f_{j}$ and $N_{d}$ dark-fields $d_{k}$, where $j\in[1,N_{f}]$ and
$k\in[1,N_{d}]$ refer to the reference images’ index. In a stationary flat-
field correction, one takes advantage only of the average flat-field and dark-
field images, $\bar{f}=\frac{1}{N_{f}}\sum_{j}f_{j}$ and
$\bar{d}=\frac{1}{N_{d}}\sum_{k}d_{k}$.

A simple procedure to correct the uneven detector response would be to correct
each raw input image $p_{i}$ by calculating:

$p_{i}^{c}=\frac{p_{i}-\bar{d}}{\bar{f}-\bar{d}}.$ (1)

As mentioned previously, due to the character of SASE processes, noise is
shot-to-shot-dependent, introducing a change in the intensity profile on the
scintillator screen between shots. This behaviour cannot be captured solely by
the average flat-field image, causing the conventional method to fail for such
conditions. To correct for the stochastic effects observed in SASE sources
such as EuXFEL, a dynamic flat-field correction substitutes the average flat-
field with an effective flat-field $f_{i}^{\prime}$, which depends on the
collected sample image $p_{i}$ itself. With this alternative procedure, the
corrected image is

$p_{i}^{c}=\frac{p_{i}-\bar{d}}{f^{\prime}_{i}-\bar{d}}.$ (2)

In principle, the true illumination function may be any image. However, some
simplifying assumptions allow for an approximate estimate of $f_{i}^{\prime}$.
Firstly, it is natural to assume that such flat-field could be expressed as a
function of the previously collected flat-fields $f_{j}$. While such a
function may be extremely complex, a simple first-order approximation would be
that $f_{i}^{\prime}$ is a linear combination of the collected flat fields.
Namely, $f_{i}^{\prime}=\sum_{j}w_{j}f_{j}$. For the purposes of an online
flat-field correction, this initial assumption leads to a fast and reliable
implementation, which may be supplemented by further corrections offline.

While one may approach the issue of the approximation of an effective flat-
field with as many weights as there are collected flat-fields, one may further
simplify such a description by assuming that only a few collected flat-fields
are the main contributors to the sum. A data-reduction technique could,
therefore, be of use to reduce the number of coefficients required. Assuming
the flat-fields collected are samples of a Gaussian distribution, one may use
Principal Components Analysis (PCA) [12, 13] to rewrite the linear flat-field
combination as a linear combination of the principal components and keep only
the components with the largest variance. Discarded components would then
contribute little to the linear combination. In general, it is feasible to use
$M$ components, where $M\ll N_{f}$. In such a setting, we expand the effective
flat-field $f^{\prime}_{i}$ as:

$f^{\prime}_{i}=\bar{f}+\sum_{m}^{M}u_{m}w_{i\,m},$ (3)

where $u_{m}$ is the $m$-th principal component of the mean-subtracted flat-
fields, $\\{f_{j}-\bar{f}\\}$, and $w_{i\,m}$ are free parameters to be
identified online.

While a functional form is available to parameterize the effective flat-field
$f_{i}^{\prime}$, one must determine a procedure to select the weights
$w_{i\,m}$, such that the effective flat-field is uniquely chosen for a given
sample image $p_{i}$. A constraint must be imposed on the corrected image for
a meaningful choice of weights.

A key constraint is to preserve details originating from the interaction of
X-rays with the sample. Such details appear through sharp transitions or edges
in the image, separating the signal from the background content. The free
parameters $w_{i\,m}$ must, therefore, preserve such edge effects, while
removing spurious variations in the images caused by the uneven detector
response. Additionally, one would like to remove spurious variations in the
images added by an uneven flat-field. The regularization technique of Total
Variation (TV) denoising is often used in image processing to achieve exactly
those goals, by minimizing the TV of the corrected image, while aiming to
obtain an image as similar as possible to the uncorrected one. Details
concerning calculation of weights $w_{i}$ can be found in the following
articles [8, 14].

The method described here consists of two separate steps. Initially, reference
flat-fields and dark-fields are acquired and PCA is used to obtain the most
relevant principal components of the flat-field dataset. During data-
acquisition with a sample, for each individual frame, weights appropriate for
the effective flat-field are found. After such calculation, the final
corrected frame is obtained by applying the weights in the linear combination
in Eq. 2.

Such a procedure must be implemented in a way that it can be run online at
adequately fast rates to keep up with the data-acquisition rate of EuXFEL. We
expand on the technical implementation in the next section.

##  3 Implementation

This section provides a schematic overview of the two steps of the main online
normalization algorithm outlined in the previous section. Selected details of
the implementation of the algorithm and specific features for visualization
applied at EuXFEL are detailed here.

The schematic overview of the first part of the algorithm, employing PCA, can
be found in Fig. 1 (a), whereas Fig. 1 (b) summarises the second part which
estimates an individual flat-field for every sample image and visualizes the
resulting normalized data.

Next, we list some of the specific functions and details used in the
implementation of the algorithm. For the total variation minimization of the
objective function a quasi-Newton algorithm is applied [15]. The precision
parameter was tuned to maximise processing rates of the correction algorithm.
Adjustments have been made to preserve the quality of the corrected images
while increasing the processing speed.

The speed of the algorithm is further increased by down-sampling all the
images entering the minimization procedure [8, 10], causing also reduction of
noise in corrected images. Here, in order not to introduce more time-consuming
steps, we omitted any additional algorithms aiming exclusively for noise
reduction, which may be included at later stages for offline processing.

Figure 1: Algorithm overviews of the first part with principal component
analysis of flat-field dataset (a) and the second part performing the online
dynamic flat-field correction and visualization (b).

This work demonstrates a novel implementation of the dynamic flat-field
correction algorithm to be used at EuXFEL’s infrastructure, which can be used
as a normalization technique during experiments. Moreover, considering the
goal of user-friendliness of the online normalization tool, EuXFEL’s software
framework Karabo [16] is used here for online visualization. One of the two
implementations uses the _metropc_ framework and the _extra-metro_ package
[17], which is integrated within the Karabo framework and its Graphical User
Interface (GUI). The _metropc_ framework allows for flexible adjusting of a
script during an experiment. It allows the visualization of the corrected
incoming data, in addition to raw sample images. The second implementation of
the online dynamic flat-field correction algorithm is done using Karabo bridge
[17] to access online data in the Karabo pipeline and perform parallelized
analysis and further visualization using the Qt toolkit [18].

##  4 Results

We have described here the results obtained by our dynamic flat-field
correction implementation at EuXFEL. We demonstrated our algorithm on a
Venturi tube dataset [19]. Data were acquired in March 2021 at the SPB/SFX
instrument [20] of EuXFEL with a photon energy of $9.3$ keV. A fast-frame-rate
Shimadzu HPV-X2 camera, coupled to a $8\mu m$ thick LYSO:Ce scintillator, was
used as an imaging system with a 10$\times$ magnification [5]. Datasets
consisting of 68 dark-field and 67 flat-field trains were taken to test our
method. Each train consisted of 128 images with a size of $250\times 400$
pixels. The effective pixel size was $3.2\mu m$. The exposure time of an image
was ultra-fast and was given by the X-ray pulse duration translated into the
latent image at the scintillator. The duration of the signal emitted by the
scintillator was given by the scintillator emission time $\tau$ (LYSO:Ce
$\tau\sim$ 40 ns). The frequency of acquired frames within a train was $1.128$
MHz and trains were collected at frequencies from $0.08-0.1$ Hz, limited by
the camera idle time. The performance of the dynamic flat-field correction
algorithm on the aforementioned dataset is demonstrated in this section.

The first step toward the normalization of online data is the processing of
flat-field and dark-field data. Our aim was to acquire hundreds to thousands
of flat-field images, which, at the current acquisition speed, requires up to
$10$ minutes of measurements. After the first step of the algorithm, as
described in Fig. 1 (a), a preview of the cumulative sum of explained variance
ratio of principal components can be viewed. This allows monitoring of the
rate at which the selected number of principal components explain variance in
the flat-field dataset. In addition, it enables the increase or decrease of
the number of components before proceeding with the normalization of sample
images.

The flat-field dataset consists of approximately 10 000 flat-field images with
a few instances shown in Fig. 2 (a), part (b) shows the mean flat-field image
calculated over the whole dataset. In Fig. 3, selected principal components
are depicted for the index of principal component $m=\\{1,2,3,16\\}$. The
first principal component, labeled $m=1$, is similar to the mean flat-field
image in Fig. 2 (b) and explains the largest fraction of the variance in the
dataset. The following two components, $m\in\\{2,3\\}$, capture the intensity
variation in horizontal and vertical directions. Components $m\leq 4$, show
mostly horizontal changes in intensity with pronounced stripe features, as
depicted on the component $m=16$. The cumulative sum of explained variance
ratio of the flat-field dataset is shown in Fig. 3 for the first $25$
principal components. Figure 3 (b) shows that the first three components
explain approximately $83\%$ of the variance. Adding more components does not
increase it significantly, and after component $3$ the increases in cumulative
sum of explained variance are minor. However, as will be discussed in more
detail later in this chapter, in order to correct for the majority of vertical
features in sample images, one needs to include also those components with
seemingly low values of explained variance, as they obtain variations needed
for a high-quality correction.

After identifying the most important components describing flat-field
illumination, one can continue with the normalization procedure of sample
images described in Section 3. To better visualize the performance of dynamic
flat-field correction on the test data, in Fig. 4 we plot the sample image
(a), its counterparts corrected by both conventional (b) and dynamic method
(c). The default value for a number of principal components taken into
consideration here is 20. Fig. 4 illustrates that the dynamic method (c) more
successfully removes vertical features than the conventional correction method
(b). We can notice vertical features along the whole width of an image
changing position from frame to frame, which is captured in examples (a)
Fig.2. Moreover, the lack of those features on the average flat-field image
(b) Fig.2 causes the inability of conventional method successfully remove such
artefacts from acquired images. Even fine stripes, appearing on average flat-
field (b) Fig.2 on the left side, are not removed using conventional
normalisation method, in Fig.4 (b), due to the slight change of their position
during measurement. Another challenging detail to eliminate is the circular
spot from the middle of the sample image (a). Similar to the fine stripes from
the previous example, better, however not perfect, reduction of the spot is
reached only with the dynamic method (c) Fig.2.

The quality of flat-field corrected images was mostly assessed visually. To
estimate the method’s performance quantitatively, we calculated pixel value
spread (pvs) as a deviation of an image from its average pixel value

$\text{pvs}=\sqrt{\frac{\sum_{l}^{N^{\prime}}(p_{l}-\mu)^{2}}{N^{\prime}}},\,\,\mu=\sum_{l}^{N^{\prime}}p_{l}$
(4)

where index $l$ iterates over pixels of image $p$. We were using only ”flat”
parts of images, without any sample, in order to estimate the amount of noise
present in an image. Higher pvs values are expected for uncorrected images
with uneven intensity profile and lower values for images, which were
corrected in a way, that reduces most of the uneven structures and noisy
pixels.

In order to assess the quality of correction of images of the same position
within a train we calculated TV and pvs for each of 128 images acquired in one
train and averaged over 15 trains. Average values for both variables and their
standard deviations are shown in Fig.5. While the pvs and TV value varies for
each frame, one may see in Fig. 5 that the dynamic flat-field corrected images
have a consistently lower TV (a) and spread (b), suggesting the lowering of
noise level in the chosen sample-free part of the images. In contrast, the
conventional flat-field correction shows very little benefit, in comparison to
no correction at all, as for TV and pvs values no statistically significant
difference was found for uncorrected image and images corrected using
conventional method. This result was expected based on the visual comparison
in Fig.4, in which no visible removal of vertical features was observed.
Standard deviation expressed as colored area for uncorrected images and both
normalization methods in Fig. 5 describes train-to-train changes in TV and pvs
value for a given frame in dataset. On average, the pvs value for a
dynamically corrected image was found to be 0.069, while the spread for
conventionally corrected images was 0.098. Uncorrected images have a pvs of
0.102.

Additionally, we studied the quality of normalized images estimated by TV
values as a function of a varying number of principal components. Resulting
corrected images were also assessed visually. Fig. 6 describes changes in TV
value due to a varying number of principal components used in the
normalization procedure. One can notice a decrease in average value of TV for
numbers of principal components $15$ and $20$, followed by an increase
starting after the $20$ components mark. For a visual comparison, Fig. 7
provides examples of flat-field corrected images using a varying number of
principal components in the normalization algorithm. Fig. 7 illustrates that
including 2-10 principal components in the normalization procedure leads to
non-optimal flat-field correction results. Fine line structures are still
visible on the left part of images for 2 and 4 component cases. Vertical
features along the whole width are not fully eliminated using 10 principal
components. Increasing number of principal components up to 10 results in
decreased visibility of vertical features, which also supports observation of
decreasing TV values in Fig. 6. The circular spot in the middle is becoming
less noticeable with higher number of principal components, although can be
detected on all corrected images up to the case using 100 components.
Including around twenty principal component seems to remove majority of the
unwanted features and results for 40 and 100 components in Fig. 6 are not
leading to visibly better results. TV values in Fig. 6 are growing after 20
components, which may be caused by overfitting. The results listed in Figs. 6
and 7 indicate that the optimal number of principal components is around $20$,
where most of the features appearing on flat fields are removed.

Considering the idle time of the camera, our goal was to reach a correction
before the next acquired train arrives. The idle time of the camera is around
$15$s, meaning the camera is able to acquire one in 150 of EuXFEL’s trains.
This consideration does not take into account set up time spent for the
analysis of flat- and dark-field datasets, as this process usually takes a few
minutes and needs to be repeated on an hourly basis, whenever new flat- and
dark-fields are taken. The proposed algorithm and both its implementations
mentioned in Sec. 3 were able to be faster and keep up with the volume of
incoming data. A parallelized version reached correction time of approximately
$1$s for a dataset consisting of one train with $128$ images. Tests were
performed at online cluster INTEL, Gold-6140 CPU @ 2.30 GHz, 768G. Moreover,
the processing speeds should be able to sustain the load with anticipated
updates, which aim to shorten the idle time of the camera. The current
performance of the algorithm is shown to be suitable for its use during
experiments as a near real-time analysis tool.

Figure 2: Example instances (a) of the test flat-field dataset consisting of
approximately 10 000 images and the average flat-field image (b) calculated
over the whole test dataset. Figure 3: An example of principal components
$u_{m}$ for $m=\\{1,2,3,16\\}$ (a), obtained from the test flat-field data
shown in Fig. 2. The cumulative sum of explained variance ratio of the first
25 principal components is depicted in subfigure (b). Figure 4: Comparison of
a sample image (a) without any normalization, conventionally flat-field
corrected image (b) and dynamic flat-field corrected data (c). Only the bottom
section of images ($(50,400)$ pixels) is shown here from originally sized
images ($(250,400)$ pixels). Figure 5: Comparison of total variation (a) and
pixel value spread (b) for uncorrected images and images corrected by
conventional and dynamical method. Values of TV and pvs are calculated
depending on their position within a train, where the maximum frame number is
128. Both TV and pvs for each frame number are averaged over 15 trains and
train-to-train changes are captured by their standard deviation expressed as a
colored area. Figure 6: Dependence of TV of images corrected by the dynamical
method using a different number of principal components. Values are calculated
for one train and are averaged over 128 images with gray colored area given by
standard deviation of frame-to-frame changes within a train. Figure 7:
Comparison of uncorrected sample image and flat-field corrected images using
the dynamical method with $[2,4,10,20,40,100]$ principal components.

##  5 Conclusion and Discussion

We have demonstrated a method for dynamic flat-field correction on MHz
repetition rate EuXFEL data. It has been shown that, for X-ray MHz microscopy
data [10], the dynamic method results in an improved correction compared to
the conventional flat-field correction method, taking into account temporal
variations in the flat-field dataset, which is characteristic for XFEL
facilities. To compare the performance of conventional and dynamic methods
applied to test data, we have calculated total variation, and the pixel value
spread for the area of images without any sample. Both quantitative comparison
and visual assessment have shown improved results with the dynamic
normalization method compared to conventional method. The normalisation
algorithm has been implemented online to visualize data and flat-field
corrected data during experiments at the SPB/SFX instrument of European XFEL.
The implementation was able to correct and visualize incoming data in near
real-time, before the next data from the MHz-frame-rate camera arrived. In
future, we plan to extend the use by implementing new data analysis tools such
as data decomposition methods and adding more features to the GUI.

## Acknowledgments

This work was performed within R&D project “MHz microscopy at EuXFEL: From
demonstration to method”, 2020 - 2022.

## References

  * 1. A. Morgan, M. Prasciolu, A. Andrejczuk, J. Krzywinski, A. Meents, D. Pennicard, H. Graafsma, A. Barty, R. Bean, M. Barthelmess, D. Oberthuer, O. Yefanov, A. Aquila, H. Chapman, and S. Bajt, “High numerical aperture multilayer laue lenses,” Scientific Reports, vol. 5, p. 9892, 06 2015.
  * 2. K. Fezzaa and Y. Wang, “Ultrafast x-ray phase-contrast imaging of the initial coalescence phase of two water droplets,” Phys. Rev. Lett., vol. 100, p. 104501, Mar 2008.
  * 3. M. P. Olbinado, X. Just, J.-L. Gelet, P. Lhuissier, M. Scheel, P. Vagovic, T. Sato, R. Graceffa, J. Schulz, A. Mancuso, J. Morse, and A. Rack, “Mhz frame rate hard x-ray phase-contrast imaging using synchrotron radiation,” Opt. Express, vol. 25, pp. 13857–13871, Jun 2017.
  * 4. M. P. Olbinado, V. Cantelli, O. Mathon, S. Pascarelli, J. Grenzer, A. Pelka, M. Roedel, I. Prencipe, A. L. Garcia, U. Helbig, D. Kraus, U. Schramm, T. Cowan, M. Scheel, P. Pradel, T. D. Resseguier, and A. Rack, “Ultra high-speed x-ray imaging of laser-driven shock compression using synchrotron light,” Journal of Physics D: Applied Physics, vol. 51, p. 055601, jan 2018\.
  * 5. P. Vagovič, T. Sato, L. Mikeš, G. Mills, R. Graceffa, F. Mattsson, P. Villanueva-Perez, A. Ershov, T. Faragó, J. Uličný, H. Kirkwood, R. Letrun, R. Mokso, M.-C. Zdora, M. P. Olbinado, A. Rack, T. Baumbach, J. Schulz, A. Meents, H. N. Chapman, and A. P. Mancuso, “Megahertz x-ray microscopy at x-ray free-electron laser and synchrotron sources,” Optica, vol. 6, pp. 1106–1109, Sep 2019.
  * 6. W. Decking et al., “A mhz-repetition-rate hard x-ray free-electron laser driven by a superconducting linear accelerator,” Nature Photonics, vol. 14, no. 6, pp. 391–397, 2020.
  * 7. S. Milton, E. Gluskin, N. Arnold, C. Benson, W. Berg, S. Biedron, M. Borland, Y. Chae, and R. Dejus, “Exponential gain and saturation of a self-amplified spontaneous emission free-electron laser,” Science, vol. 292, 06 2001.
  * 8. V. V. Nieuwenhove, J. D. Beenhouwer, F. D. Carlo, L. Mancini, F. Marone, and J. Sijbers, “Dynamic intensity normalization using eigen flat fields in x-ray imaging,” Opt. Express, vol. 23, pp. 27975–27989, Oct 2015.
  * 9. J. Hagemann, M. Vassholz, H. Hoeppe, M. Osterhoff, J. M. Rosselló, R. Mettin, F. Seiboth, A. Schropp, J. Möller, J. Hallmann, C. Kim, M. Scholz, U. Boesenberg, R. Schaffer, A. Zozulya, W. Lu, R. Shayduk, A. Madsen, C. G. Schroer, and T. Salditt, “Single-pulse phase-contrast imaging at free-electron lasers in the hard X-ray regime,” Journal of Synchrotron Radiation, vol. 28, pp. 52–63, Jan 2021.
  * 10. K. Buakor, Y. Zhang, S. Birnsteinova, V. Bellucci, T. Sato, H. Kirkwood, A. P. Mancuso, P. Vagovic, and P. Villanueva-Perez, “Shot-to-shot flat-field correction at x-ray free-electron lasers,” Opt. Express, vol. 30, pp. 10633–10644, Mar 2022.
  * 11. D. S. Hodge, A. F. T. Leong, S. Pandolfi, K. Kurzer-Ogul, D. S. Montgomery, H. Aluie, C. Bolme, T. Carver, E. Cunningham, C. B. Curry, M. Dayton, F.-J. Decker, E. Galtier, P. Hart, D. Khaghani, H. J. Lee, K. Li, Y. Liu, K. Ramos, J. Shang, S. Vetter, B. Nagler, R. L. Sandberg, and A. E. Gleason, “Multi-frame, ultrafast, x-ray microscope for imaging shockwave dynamics,” Opt. Express, vol. 30, pp. 38405–38422, Oct 2022.
  * 12. K. Pearson F.R.S., “Liii. on lines and planes of closest fit to systems of points in space,” The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, vol. 2, no. 11, pp. 559–572, 1901.
  * 13. H. Hotelling, “Relations between two sets of variates,” Biometrika, vol. 28, no. 3/4, pp. 321–377, 1936.
  * 14. A. Chambolle, “An algorithm for total variation minimization and applications,” J. Math. Imaging Vis., vol. 20, p. 89–97, jan 2004.
  * 15. R. Byrd, P. Lu, J. Nocedal, and C. Zhu, “A limited memory algorithm for bound constrained optimization,” SIAM Journal of Scientific Computing, vol. 16, pp. 1190–1208, Sept. 1995.
  * 16. B. Heisen, D. Boukhelef, S. Esenov, S. Hauf, I. Kozlova, L. Maia, A. Parenti, J. Szuba, K. Weger, K. Wrona, and C. Youngman, “Karabo: An Integrated Software Framework Combining Control, Data Management, and Scientific Computing Tasks,” in Proceedings, 14th International Conference on Accelerator $\&$ Large Experimental Physics Control Systems (ICALEPCS 2013) : San Francisco, CA, USA, October 6-11, 2013, p. FRCOAAB02, 14th International Conference on Accelerator $\&$ Large Experimental Physics Control Systems, San Francisco (USA), 6 Oct 2013 - 11 Oct 2013, Oct 2013.
  * 17. H. Fangohr et al., “Data Analysis Support in Karabo at European XFEL,” in Proc. of International Conference on Accelerator and Large Experimental Control Systems (ICALEPCS’17), Barcelona, Spain, 8-13 October 2017, no. 16 in International Conference on Accelerator and Large Experimental Control Systems, (Geneva, Switzerland), pp. 245–252, JACoW, Jan. 2018.  https://doi.org/10.18429/JACoW-ICALEPCS2017-TUCPA01.
  * 18. PyQT, “Pyqt reference guide,” 2012.
  * 19. H. Soyama and J. Hoshino, “Enhancing the aggressive intensity of hydrodynamic cavitation through a venturi tube by increasing the pressure in the region where the bubbles collapse,” AIP Advances, vol. 6, p. 045113, 04 2016.
  * 20. A. P. Mancuso, A. Aquila, L. Batchelor, R. J. Bean, J. Bielecki, G. Borchers, K. Doerner, K. Giewekemeyer, R. Graceffa, O. D. Kelsey, Y. Kim, H. J. Kirkwood, A. Legrand, R. Letrun, B. Manning, L. Lopez Morillo, M. Messerschmidt, G. Mills, S. Raabe, N. Reimers, A. Round, T. Sato, J. Schulz, C. Signe Takem, M. Sikorski, S. Stern, P. Thute, P. Vagovič, B. Weinhausen, and T. Tschentscher, “The Single Particles, Clusters and Biomolecules and Serial Femtosecond Crystallography instrument of the European XFEL: initial installation,” Journal of Synchrotron Radiation, vol. 26, pp. 660–676, May 2019.


# AI-Assisted Diagnosis for Covid-19 CXR Screening: From Data Collection to
Clinical Validation

###### Abstract

In this paper, we present the major results from the Covid Radiographic
imaging System based on AI (Co.R.S.A.) project, which took place in Italy.
This project aims to develop a state-of-the-art AI-based system for diagnosing
Covid-19 pneumonia from Chest X-ray (CXR) images. The contributions of this
work are manyfold: the release of the public CORDA dataset, a deep learning
pipeline for Covid-19 detection, and the clinical validation of the developed
solution by expert radiologists. The proposed detection model is based on a
two-step approach that, paired with state-of-the-art debiasing, provides
reliable results. Most importantly, our investigation includes the actual
usage of the diagnosis aid tool by radiologists, allowing us to assess the
real benefits in terms of accuracy and time efficiency. Project homepage:
https://corsa.di.unito.it/

Index Terms—  Covid-19, CAD, Clinical Impact, Deep Learning, Medical Image
Analysis

(a)

(b)

Fig. 1: AI-assisted diagnosis help radiologists in making more accurate and
faster dagnosis. In (a) we report the average ROC curve of a pool of
radiologists during external validation, in a blind setting (no AI) and in an
AI-assisted setting. In (b) we compare the average diagnosis time (per image)
in a blind vs. assisted setting. We notice that the average AUC increases in
the AI-assisted setting, while the diagnosis is quicker on average.

##  1 Introduction

During the beginning of 2020, Covid-19 rapidly spread in China and out into
multiple countries worldwide [1]. As of November 2023, there were more than
770 million confirmed cases worldwide, with almost 7 million
deaths111https://covid19.who.int/. It is now widely recognized how
radiological imaging can contribute to identifying Covid-19 even in the early
stages of the disease: Chest X-rays (CXR) and computed tomography (CT) can be
helpful diagnostic means. Clinical evidence shows that CT represents the gold
standard for diagnosing Covid-19 [2]. However, its use in the initial
management of the suspected patient is limited by logistical problems
(sanitization of the room/suitable dressing/undressing areas of the
operators), radioprotection, and, finally, costs. For these reasons, CXR is
often chosen in emergency departments as the first radiological method for the
ease of execution, the low dose of ionizing radiation, and the more contained
costs (if compared to CT).

In 2022 the Piedmont Region funded the Covid Radiographic imaging System based
on AI (Co.R.S.A.) project, a collaborative initiative by the University of
Turin, two hospital radiological units (A.O.U Cittá della Salute e della
Scienza, Azienda Sanitaria ASL TO3) and an enterprise (REGOLA), to develop and
validate in the field a state-of-the-art AI-based system for aiding the
diagnosis of Covid-19 pneumonia from CXR images. In this work, we describe the
main project milestones and contributions:

  1. 1.

the release of the CORDA dataset;

  2. 2.

a robust deep learning pipeline for Covid-19 detection;

  3. 3.

an external clinical validation by expert radiologists.

Among the mentioned contributions, the latter is probably the most impactful:
indeed, papers on Covid-19 identification from CXR proliferated in recent
years, but in most cases a complete validation of the system by radiology
doctors is lacking. Although the emergency setting of the Covid-19 pandemic
peak is fortunately over, the groundwork built by this project can serve as a
solid foundation for quick-starting future responses to epidemics, should the
need arise.

##  2 The CORDA dataset

CORDA [3] contains 3843 images of different modalities, with 1601 CXR and 2242
CT images222In the scope of this work, we only considered CXR images.. The
dataset was made publicly available in January
2023333https://zenodo.org/record/7821611. This dataset aims to provide a
multi-center collection of radiographic images for Covid-19 detection, in
order to build more robust machine learning algorithms and models. The dataset
curation is part of the ongoing project Co.R.S.A. CORDA comprises four
different Italian hospitals:

  1. 1.

A.O.U. Città della Salute e delle Scienza (Molinette), Torino (CDSS);

  2. 2.

A.O.U. San Luigi Gonzata, Orbassano (SLG);

  3. 3.

A.O. Mauriziano, Torino (MRZ);

  4. 4.

Centro Cardiologico Monzino, Milano (MNZ).

Dataset composition In Tab. 1 we show the distributions of Covid-19 positive
and negative CXR cases across the four institutions that contributed to CORDA;
the number of images acquired by Computed Radiography (CR) and Digital
Radiography (DR) is reported as well. It can be observed that the samples of
different institutes are inevitably unbalanced and potentially present
different imaging characteristics. As shown in [4, 5] this can negatively
impact deep learning models, limiting their generalization capability. These
biases will be considered in Sec. 3.1 with proper regularization techniques to
limit their effect in the model training phase.

Institution | Pos. cases | Neg. cases | CR/DR
---|---|---|---
CDSS | 362 | 183 | 401/144
SLG | 250 | 477 | 713/14
MRZ | 138 | 35 | 163/10
MNZ | 156 | 0 | 63/93
Total | 906 | 695 | 1340/261
Table 1: CORDA composition: Covid-19 positive and negative patients, imaging
modality CR/DR.

##  3 Automatic Diagnosis of Covid-19

In this section, we present a method for Covid-19 diagnosis based on deep
learning (DL). The method that we propose consists of two steps: first, we
pretrain a deep neural network (DNN) on a large-scale CXR dataset, with the
aim of detecting objective radiological findings, then we apply transfer
learning to train a Covid-19 classifier on the CORDA dataset. As previously
shown in [6], this approach provides reliable results. Additionally, we also
employ a state-of-the-art debiasing technique, in order to address the
presence of possible biases and spurious information related to the
acquisition site in the data, a phenomenon commonly known as _site-effect_
[4].

###  3.1 Method

| CDSS | SLG | MRZ | MNZ | Avg
---|---|---|---|---|---
Baseline | 0.68 | 0.83 | 0.74 | 0.87 | 0.78
FairKL | 0.70 | 0.85 | 0.77 | 0.88 | 0.80
Table 2: Results of Covid-19 detection on the CORDA dataset, in terms of
balanced accuracy.

Pretraining on objective radiological findings. We leverage a large-scale
dataset, _CheXpert_ [7], which contains annotation for different kinds of
common radiological findings that can be observed in CXR images (like opacity,
pleural effusion, cardiomegaly, etc.). This large dataset is well suited for
multi-label classification tasks; in fact, more than one finding can be
commonly observed simultaneously in ill patients’ lungs. CheXpert provides 14
different types of observations for each image in the dataset. We use the
setup described in [6], consisting of a slightly modified DenseNet-121 [8]
architecture.
Covid-19 Prediction To obtain a Covid-19 classifier, we employ the DNN encoder
pre-trained on CheXpert as a frozen feature extractor on the CORDA dataset. We
train a fully connected binary classifier for the final prediction, using the
standard binary cross-entropy loss (BCE).
Site regularization From a debiasing point of view, we define as bias-aligned
all those samples that share the same acquisition site w.r.t. to a given
_anchor_ sample, and bias-conflicting the ones which do not. Also, we define
as _positive_ the samples which share the same target class as the anchor. We
employ the recently proposed FairKL [9] regularization technique, which aims
at minimizing the Kullback-Leibler divergence of the distance distributions in
the latent space of positive bias-aligned $B_{+,b}$ and positive bias-
conflicting and $B_{+,b^{\prime}}$:

$\mathcal{R}^{FairKL}=D_{KL}(B_{+,b}||B_{+,b^{\prime}})$ (1)

This regularization term aims to make samples of the same class
indistinguishable in the latent space based on the acquisition site. This
choice is motivated by the possible correlation between the target label,
modality, and source institution, as shown in Tab. 1. The final objective
function we optimize is thus
$J=\mathcal{L}^{BCE}+\lambda\mathcal{R}^{FairKL}$, with $\lambda\geq 0$.

###  3.2 Results

We train the final classifier (on the frozen backbone) for 100 epochs, using
SGD as optimizer, with a learning rate of 0.01, decayed with a cosine
schedule, a batch size of 32 and an image size of 448x448. We also employ an
augmentation scheme composed of random crop, rotation, and cutout. For each
experiment, we perform a 4-fold cross-validation on the CORDA dataset, and we
report the mean results in Tab. 2, evaluated separately on each institution
test set. We report as _baseline_ a simple transfer learning employing only
BCE for training. With this approach, we achieve a 78% balanced accuracy in
Covid-19 classification. By employing FairKL (with $\lambda=1$), we achieve a
consistent improvement with respect to the baseline, improving the results for
each institution and achieving an average balance accuracy of 80%.

##  4 Clinical Validation

Fig. 2: Example of AI report in our custom DICOM viewer. This report is shown
to the evaluating radiologist alongside the CXR image. The information
includes the predicted probability of Covid-19 infection and other relevant
lung pathologies. For easiness of readability, probabilities $>0.5$ are marked
in red (except for the ”No Finding” class).

In this section, we focus on assessing the impact of AI-aided diagnosis in the
radiologists’ workflow, one of the key aspects of our project. To this end,
the clinical validation focuses on measuring two KPIs: _i)_ the accuracy of
the radiologists’ diagnosis _ii.)_ the time needed to formulate it.

###  4.1 Materials and tools

Validation data As validation data, 100 external images were collected at ASL
To3, of which 50 are from Covid-19 patients and 50 are control cases. These
images are not part of the CORDA dataset and were not used for previous
phases. This external validation also allows us to assess whether the proposed
method can generalize to novel data and sites.
Dicom viewer For this validation, we developed a custom DICOM image viewer
which allows DNNs predictions to be shown alongside the analyzed image. Our
software automatically collects the diagnosis made by the radiologist for each
image and the time taken to make it. Primary functionalities of a standard
DICOM viewer are included, such as image manipulation (e.g. zoom, translation)
and adjustable windowing (VOI LUT). Fig. 2 depicts an example of DNNs
predictions as shown to the evaluating radiologist in our viewer.

(a) ROC-AUC

(b) Diagnosis Time

Fig. 3: Break-out of diagnosis performance by radiologist. (a) shows the AUC
score of blind and AI-assisted evaluation. While for some radiologists,
especially those with higher blind scores we notice a slight decrease in AUC,
the improvement is consistent in those who achieve a lower base score. On
average, the AUC increases from 0.85 to 0.88, as shown in Fig.1. (b) compares
the blind diagnosis time with the AI-assisted diagnosis time. The blue line
represents a linear regression line, and the dashed black line represents the
identity. We can see that, on average, the diagnosis time decreases as the
fitted line lies below the identity for almost all radiologists.

###  4.2 Validation Protocol

To perform the clinical evaluation, we employ a pool of 6 radiologists from
Ospedale di Rivoli and ASL To3, with different experience levels. They are
summarized in Tab. 3. For each image, a severity score [10] is assigned by
each radiologist upon visual inspection, ranging from 0 (healthy) to 18
(maximum severity). During diagnosis, the radiologists do not know the correct
label for a given patient. The experiment is repeated in two settings, with a
wash-out period in between:

  * •

_blind_ setting: radiologists can only leverage the image to make a diagnosis

  * •

_AI-assisted_ setting: the predicted probabilities of our DNNs are shown
alongside the image.

Initials | Affiliation | Years of Experience
---|---|---
S. N. | ASL To3 Rivoli | 1
V. N. | ASL To3 Rivoli | 1
D. C. | ASL To3 Rivoli | 2
M. B. | ASL To3 Rivoli | 5
F. L. | ASL To3 Rivoli | 10
P. M. | ASL To3 Rivoli | 15
Table 3: List of the radiologists involved in the study.

###  4.3 Validation Results

Fig. 1 shows the overall results of the validation in terms of prediction AUC
and diagnosis time. In Fig. 1(a), the AUC is computed on the radiologists’
average severity score assigned for each patient. We observe an improvement in
the assisted setting compared to the blind evaluation. In Fig. 1(b), we plot
the diagnosis time in the two settings: on average (blue line), the time
required for diagnosis decreases in the assisted setting. In Fig. 3, we show a
breakout of the results for each participating radiologist for prediction AUC
and diagnosis time. In most cases, we observe an improved performance in the
assisted setting, with some significant improvement (e.g. AUC increases from
0.85 to 0.96 for radiologist #1). These results show that by employing AI
systems in a real-world clinical setting, more accurate and faster diagnosis
can be achieved.

##  5 Conclusions

In this work, we have presented the results from the Co.R.S.A. projects, which
aimed at _i.)_ building an open and freely available data collection for
Covid-19 diagnosis from CXR images _ii.)_ developing a robust deep-learning-
based pipeline for automatic classification _iii.)_ performing a clinical
validation in a real-world context, with a pool of expert radiologists. This
project represents the results of our efforts in Covid-19 detection over the
past few years. Furthermore, while the critical phase of the Covid-19 pandemic
has passed, the foundation built by this project can serve as a robust basis
for promptly initiating responses to future epidemics if necessary, thanks to
the collaboration developed with different major hospitals and radiology
units.

Compliance with ethical standards This study was performed in line with the
principles of the Declaration of Helsinki. Approval was granted by Comitato
Etico Interaziendale A.O.U. San Luigi di Orbassano, AA.SS.LL. TO3,TO4,TO5
(10/24/2022/ No. 153/2022).

Acknowledgments The project was funded through the “INFRA-P2 - Potenziamento
di laboratori di prova ed infrastrutture di ricerca già esistenti nella
disponibilità di organismi di ricerca pubblici e progetti di ricerca e
sviluppo finalizzati al contrasto della pandemia Covid-19” call for proposals.
A special acknowledgement goes to the other partners of the Co.R.S.A. project:
Città della Scienza e della Salute and REGOLA s.r.l.. We wish to tank the
valuable support of all the radiology units that contributed to the CORDA
dataset.

## References

  * [1] Z. Zu et al.,  “Coronavirus disease 2019 (covid-19): A perspective from china,”  Radiology, 2020.
  * [2] H. Shi et al.,  “Radiological findings from 81 patients with covid-19 pneumonia in wuhan, china: a descriptive study,”  The Lancet Infectious Diseases, 2020.
  * [3] M. Alesina et al.,  “Corda dataset,” https://doi.org/10.5281/zenodo.7821611, Apr. 2023.
  * [4] B. Glocker et al.,  “Machine learning with multi-site imaging data: An empirical study on the impact of scanner effects,”  2019\.
  * [5] E. Tartaglione et al.,  “Unveiling covid-19 from chest x-ray with deep learning: a hurdles race with small data,”  IJERPH, 2020.
  * [6] C. A. Barbano et al.,  “A two-step radiologist-like approach for covid-19 computer-aided diagnosis from chest x-ray images,”  in ICIAP, 2022.
  * [7] J. Irvin et al.,  “Chexpert: A large chest radiograph dataset with uncertainty labels and expert comparison,”  in AAAI, 2019.
  * [8] G. Huang et al.,  “Densely connected convolutional networks,”  in CVPR, 2017.
  * [9] C. A. Barbano et al.,  “Unbiased supervised contrastive learning,”  in ICLR, 2023.
  * [10] A. Borghesi and R. Maroldi,  “Covid-19 outbreak in italy: experimental chest x-ray scoring system for quantifying and monitoring disease progression,”  La Radiologia Medica, 2020.