#  Theoretical realisation of a two qubit quantum controlled-not logic gate
and a single qubit Hadamard logic gate in the anti-Jaynes-Cummings model

Christopher Mayero  Department of Physics and Materials Science, Maseno
University, Private Bag-40105, Maseno, Kenya  Joseph Akeyo Omolo  Department
of Physics and Materials Science, Maseno University, Private Bag-40105,
Maseno, Kenya  Onyango Stephen Okeyo  Department of Physics and Materials
Science, Maseno University, Private Bag-40105, Maseno, Kenya

###### Abstract

We provide a theoretical scheme for realizing a Hadamard and a quantum
controlled-NOT logic gates operations in the anti-Jaynes-Cummings interaction
process. Standard Hadamard operation for a specified initial atomic state is
achieved by setting a specific sum frequency and photon number in the anti-
Jaynes-Cummings qubit state transition operation with the interaction
component of the anti- Jaynes-Cummings Hamiltonian generating the state
transitions. The quantum controlled-NOT logic gate is realised when a single
atomic qubit defined in a two-dimensional Hilbert space is the control qubit
and two non-degenerate and orthogonal polarised cavities defined in a two-
dimensional Hilbert space make the target qubit. With precise choice of
interaction time in the anti-Jaynes-Cummings qubit state transition operations
defined in the anti-Jaynes-Cummings sub-space spanned by normalised but non-
orthogonal basic qubit state vectors, we obtain ideal unit probabilities of
success in the quantum controlled-NOT operations.

Keywords: anti-Jaynes-Cummings, Jaynes-Cummings, Hadamard, controlled-NOT

##  I Introduction

In quantum computers, quantum bits (qubits) [1, 2] are the elementary
information carriers. In such a computer, quantum gates [1, 2, 3] can
manipulate arbitrary multi-partite quantum states [4] including arbitrary
superposition of the computational basis states, which are frequently also
entangled. Thus the logic gates of quantum computation are considerably more
varied than the logic gates of classical computation. In addition, a quantum
computer can solve problems exponentially faster than any classical computer
[5] because by exploiting superposition principle and entanglement allows the
computer to manipulate and store more bits of information than a classical
computer.

In this paper we present a theoretical approach of realizing Hadamard and
controlled-NOT (C-NOT) quantum logic gates which form a universal set for
quantum computation [6, 7, 8]. The important discovery and proof of a
conserved excitation number operator of the AJC Hamiltonian [9] now means that
dynamics generated by the AJC Hamiltonian is exactly solvable, as demonstrated
in the polariton and anti-polariton qubit (photospin qubit) models in [10,
11]. The reformulation developed in [9, 10, 11], drastically simplifies exact
solutions of the AJC model which we shall apply in the present work.

We define the quantum C-NOT gate as that which affects the unitary operation
on two qubits which in a chosen orthonormal basis in $\mathbb{C}^{2}$ gives
the C-NOT operation obtained as

$|a\rangle|b\rangle\rightarrow|a\rangle|a\oplus|b\rangle$ (1)

where $|a\rangle$ is the control qubit, $|b\rangle$ is the target qubit and
$\oplus$ indicates addition modulo 2 [12, 2, 1]. The C-NOT gate transforms
superposition into entanglement thus acts as a measurement gate [1, 2, 12]
fundamental in performing algorithms in quantum computers [13]. Transformation
to a separable state (product state) is realized by applying the C-NOT gate
again. In this case, it is used to implement Bell measurement on the two
qubits [14].

We note here that the JC model has been applied extensively in implementing
C-NOT and Hadamard gate operations. Domokos et al (1995) [15] showed that
using induced transitions between dressed states, it is possible to implement
a C-NOT gate in which a cavity containing at most one photon is the control
qubit and the atom is the target qubit. Later, Vitali, D. et al (2001)[16]
proposed a scheme of implementing a C-NOT gate between two distinct but
identical cavities, acting as control and target qubits respectively. By
passing an atom prepared initially in ground state consecutively between the
two cavities a C-NOT ($cavity\rightarrow{atom}$) and a C-NOT
($atom\rightarrow{cavity}$) is realised with the respective classical fields
S. Saif, F. et al (2001) [17] presented a study of quantum computing by
engineering non-local quantum universal gates based on interaction of a two-
level atom with two modes of electromagnetic field in high Q superconducting
cavity. The two-level atom acted as the control qubit and the two-mode
electromagnetic field served as the the target qubit. In this letter, we apply
an approach similar to that in [17] where we implement a quantum C-NOT gate
operation between two cavities defined in a two-dimensional Hilbert space
spanned by the state vectors $|\mu_{1}\rangle=|1_{A},0_{B}\rangle$ and
$|\mu_{2}\rangle=|0_{A},1_{B}\rangle$ as target qubits. Here $|\mu_{1}\rangle$
expresses the presence of one photon in mode A, when there is no photon in
mode B, and $|\mu_{2}\rangle$ indicates that mode A is in vacuum state and one
photon is present in mode B. The control qubit in this respect is a two-level
atom. The important difference with the approach used in [17] is the model,
i.e, while the initial absolute atom-field ground state $|g,0\rangle$ in the
AJC interaction is affected by atom-cavity coupling, the ground state
$|g,0\rangle$ in the JC model [17] is not affected by atom-cavity coupling. A
similar result was determined independently in [11]. Further, with precise
choice of interaction time in the AJC qubit state transition operations
defined in the AJC qubit sub-space spanned by normalised but non-orthogonal
basics qubit state vectors [10, 11], C-NOT gate operations are realized
between the two cavities.

The Hadamard gate also known as the Walsh-Hadamard gate is a single qubit gate
[2, 1]. The Hadamard transformation is defined as

$\hat{H}=\frac{\hat{\sigma}_{\textit{x}}+\hat{\sigma}_{\textit{z}}}{\sqrt{2}}$
(2)

where it transforms atomic computational basis states $|e\rangle(|0\rangle)$,
$|g\rangle(|1\rangle)$ into diagonal basis states according to

$\displaystyle\hat{H}|e\rangle\rightarrow\frac{|e\rangle+|g\rangle}{\sqrt{2}}\quad;\quad\hat{H}|g\rangle\rightarrow\frac{|e\rangle-|g\rangle}{\sqrt{2}}$
$\displaystyle\hat{H}|0\rangle\rightarrow\frac{|0\rangle+|1\rangle}{\sqrt{2}}\quad;\quad\hat{H}|1\rangle\rightarrow\frac{|0\rangle-|1\rangle}{\sqrt{2}}$
(3)

Vitali, D. et al (2001)[16] showed that one qubit operations can be
implemented on qubits represented by two internal atomic states because it
amounts to applying suitable Rabi pulses. He demonstrated that the most
practical solution on implementing one qubit operations on two Fock states is
sending the atoms through the cavity. If the atom inside the cavity undergoes
a $\frac{\pi}{2}$ pulse one realizes a Hadamard-phase gate. Saif, F. et al
(2001) [17] also showed that it is possible to realise Hadamard operation by a
controlled interaction between a two-mode high Q electromagnetic cavity field
and a two-level atom. In his approach, the two-level atom is the control
qubit, whereas the target qubit is made up of two modes of cavity field.
Precision of the gate operations is realised by precise selection of
interaction times of the two-level atom with the cavity mode. In this paper,
we show that Hadamard operation in the AJC interaction is possible for a
specified initial atomic state by setting a specific sum frequency and photon
number in the anti-Jaynes-Cummings qubit state transition operation [10, 11],
noting that the interaction components of the anti-Jaynes-Cummings Hamiltonian
generates state transitions.

The content of this paper is therefore summarised as follows. Section II
presents an overview of the theoretical model. In sections III and IV
respectively, implementation of a quantum C-NOT and Hadamard gates in the AJC
interaction are presented. Finally section V contains the conclusion.

##  II The model

The quantum Rabi model of a quantized electromagnetic field mode interacting
with a two-level atom is generated by the Hamiltonian [9]

$\hat{H}_{R}=\frac{1}{2}\hbar\omega(\hat{a}^{\dagger}\hat{a}+\hat{a}\hat{a}^{\dagger})+\hbar\omega_{0}\hat{s}_{z}+\hbar\lambda(\hat{a}+\hat{a}^{\dagger})(\hat{s}_{+}+\hat{s}_{-})$
(4)

noting that the free field mode Hamiltonian is expressed in normal and anti-
normal order form
$\frac{1}{2}\hbar\omega(\hat{a}^{\dagger}\hat{a}+\hat{a}\hat{a}^{\dagger})$.
Here, $\omega\hskip 2.84526pt,\hskip 2.84526pt\hat{a}\hskip 2.84526pt,\hskip
2.84526pt\hat{a}^{\dagger}$ are quantized field mode angular frequency,
annihilation and creation operators, while $\omega_{0}\hskip 2.84526pt,\hskip
2.84526pt\hat{s}_{z}\hskip 2.84526pt,\hskip 2.84526pt\hat{s}_{+}\hskip
2.84526pt,\hskip 2.84526pt\hat{s}_{-}$ are atomic state transition angular
frequency and operators. The Rabi Hamiltonian in eq. (4) is expressed in a
symmetrized two-component form [9, 10, 11]

$\hat{H}_{R}=\frac{1}{2}(\hat{H}+\hat{\overline{H}})$ (5)

where $\hat{H}$ is the standard JC Hamiltonian interpreted as a polariton
qubit Hamiltonian expressed in the form [9]

$\displaystyle\hat{H}$ $\displaystyle=$
$\displaystyle\hbar\omega\hat{N}+2\hbar\lambda\hat{A}-\frac{1}{2}\hbar\omega\quad;\quad\hat{N}=\hat{a}^{\dagger}\hat{a}+\hat{s}_{+}\hat{s}_{-}$
$\displaystyle\hat{A}$ $\displaystyle=$
$\displaystyle\alpha\hat{s}_{z}+\hat{a}\hat{s}_{+}+\hat{a}^{\dagger}\hat{s}_{-}\quad;\quad\alpha=\frac{\omega_{0}-\omega}{2\lambda}$
(6)

while $\hat{\overline{H}}$ is the AJC Hamiltonian interpreted as an anti-
polariton qubit Hamiltonian in the form [9]

$\displaystyle\hat{\overline{H}}$ $\displaystyle=$
$\displaystyle\hbar\omega\hat{\overline{N}}+2\hbar\lambda\hat{\overline{A}}-\frac{1}{2}\hbar\omega\quad;\quad\hat{\overline{N}}=\hat{a}\hat{a}^{\dagger}+\hat{s}_{-}\hat{s}_{+}$
$\displaystyle\hat{\overline{A}}$ $\displaystyle=$
$\displaystyle\overline{\alpha}\hat{s}_{z}+\hat{a}\hat{s}_{-}+\hat{a}^{\dagger}\hat{s}_{+}\quad;\quad\overline{\alpha}=\frac{\omega_{0}+\omega}{2\lambda}$
(7)

In eqs. (6) and (7), $\hat{N}$, $\hat{\overline{N}}$ and $\hat{A}$,
$\hat{\overline{A}}$ are the respective polariton and anti-polariton qubit
conserved excitation numbers and state transition operators.

Following the physical property established in [11], that for the field mode
in an initial vacuum state only an atom entering the cavity in an initial
excited state $|e\rangle$ couples to the rotating positive frequency field
component in the JC interaction mechanism, while only an atom entering the
cavity in an initial ground state $|g\rangle$ couples to the anti-rotating
negative frequency field component in an AJC interaction mechanism, we
generally take the atom to be in an initial excited state $|e\rangle$ in the
JC model and in an initial ground state $|g\rangle$ in the AJC model.

Considering the AJC dynamics, applying the state transition operator
$\hat{\overline{A}}$ from eq. (7) to the initial atom-field n-photon ground
state vector $|g,n\rangle$, the basic qubit state vectors $|\psi_{gn}\rangle$
and $|\overline{\phi}_{gn}\rangle$ are determined in the form (n=0,1,2,….)
[11]

$|\psi_{gn}\rangle=|g,n\rangle\quad;\quad|\overline{\phi}_{gn}\rangle=-\overline{c}_{gn}|g,n\rangle+\overline{s}_{gn}|e,n+1\rangle$
(8)

with dimensionless interaction parameters $\overline{c}_{gn}$,
$\overline{s}_{gn}$ and Rabi frequency $\overline{R}_{gn}$ defined as

$\displaystyle\overline{c}_{gn}$ $\displaystyle=$
$\displaystyle\frac{\overline{\delta}}{2\overline{R}_{gn}}\quad;\quad\overline{s}_{gn}=\frac{2\lambda\sqrt{n+1}}{\overline{R}_{gn}}\quad;\quad\overline{R}_{gn}=2\lambda{\overline{A}_{gn}}$
$\displaystyle\overline{A}_{gn}$ $\displaystyle=$
$\displaystyle\sqrt{(n+1)+\frac{\overline{\delta}^{2}}{16\lambda^{2}}}\quad;\quad\overline{\delta}=\omega_{0}+\omega$
(9)

where we have introduced sum frequency $\overline{\delta}=\omega_{0}+\omega$
to redefine $\overline{\alpha}$ in eq. (7).

The qubit state vectors in eq. (8) satisfy the qubit state transition
algebraic operations

$\hat{\overline{A}}|\psi_{gn}\rangle=\overline{A}_{gn}|\overline{\phi}_{gn}\rangle\quad;\quad\hat{\overline{A}}|\overline{\phi}_{gn}\rangle=\overline{A}_{gn}|\psi_{gn}\rangle$
(10)

In the AJC qubit subspace spanned by normalized but non-orthogonal basic qubit
state vectors $|\psi_{gn}\rangle$ , $|\overline{\phi}_{gn}\rangle$ the basic
qubit state transition operator $\hat{\overline{\varepsilon}}_{g}$ and
identity operator $\hat{\overline{I}}_{g}$ are introduced according to the
definitions [11]

$\hat{\overline{\varepsilon}}_{g}=\frac{\hat{\overline{A}}}{\overline{A}_{gn}}\quad;\quad\hat{\overline{I}}_{g}=\frac{\hat{\overline{A}}^{2}}{\overline{A}_{gn}^{2}}\quad\Rightarrow\quad\hat{\overline{I}}_{g}=\hat{\overline{\varepsilon}}_{g}^{2}$
(11)

which on substituting into eq. (10) generates the basic qubit state transition
algebraic operations

$\displaystyle\hat{\overline{\varepsilon}}_{g}|\psi_{gn}\rangle$
$\displaystyle=$
$\displaystyle|\overline{\phi}_{gn}\rangle\quad;\quad\hat{\overline{\varepsilon}}_{g}|\overline{\phi}_{gn}\rangle=|\psi_{gn}\rangle$
$\displaystyle\hat{\overline{I}}_{g}|\psi_{gn}\rangle$ $\displaystyle=$
$\displaystyle|\psi_{gn}\rangle\quad;\quad\hat{\overline{I}}_{g}|\overline{\phi}_{g}\rangle=|\overline{\phi}_{g}\rangle$
(12)

The algebraic properties
$\hat{\overline{\varepsilon}}^{2k}=\hat{\overline{I}}_{g}$ and
$\hat{\overline{\varepsilon}}^{2k+1}=\hat{\overline{\varepsilon}}_{g}$ easily
gives the final property [11]

$e^{-i\theta\hat{\overline{\varepsilon}}_{g}}=\cos(\theta)\hat{\overline{I}}_{g}-i\sin(\theta)\hat{\overline{\varepsilon}}_{g}$
(13)

which is useful in evaluating time-evolution operators.

The AJC qubit Hamiltonian defined within the qubit subspace spanned by the
basic qubit state vectors $|\psi_{gn}\rangle$ , $|\overline{\phi}_{gn}\rangle$
is then expressed in terms of the basic qubit states transition operators
$\hat{\overline{\varepsilon}}_{g}$, $\hat{\overline{I}}_{g}$ in the form [11]

$\hat{\overline{H}}_{g}=\hbar\omega(n+\frac{3}{2})\hat{\overline{I}}_{g}+\hbar\overline{R}_{gn}\hat{\overline{\varepsilon}}_{g}$
(14)

##  III Quantum c-not gate operations

In order to realise a C-NOT quantum gate operation in this context, we take a
two-level atom as the control qubit, which is defined in a two dimensional
Hilbert space with $|e\rangle$ and $|g\rangle$ as basis vectors, where
$|e\rangle$ expresses the excited state of the two-level atom and $|g\rangle$
indicates the ground state. Two non-degenerate and orthogonal polarized cavity
modes $C_{A}$ and $C_{B}$ make the target qubit. The target qubit is defined
in two-dimensional Hilbert space spanned by the state vector
$|\mu_{1}\rangle=|1_{A},0_{B}\rangle$, which expresses the presence of one
photon in mode A, when there is no photon in mode B, and the state vector
$|\mu_{2}\rangle=|0_{A},1_{B}\rangle$, which indicates that mode A is in the
vacuum state and one photon is present in mode B.

With reference to the AJC qubit state transition operator in eq.(13), lets
first consider when an atom in ground state $|g\rangle$ enters an
electromagnetic cavity with mode A in vacuum state and a single photon in mode
B. The atom couples to the anti-rotating negative frequency component of the
field mode undergoing an AJC qubit state transition. After the atom interacts
with mode A for a time $t=\frac{\pi}{\overline{R}_{g0}}$, equal to half Rabi
oscillation time, the driving field is modulated such that
$\theta=\overline{R}_{g0}t=2\lambda\overline{A}_{g0}t=\pi$. Redefining [11]

$\overline{\alpha}=\frac{\overline{\delta}}{2\lambda}=\frac{\omega_{0}-\omega+2\omega}{2\lambda}=\frac{\delta}{2\lambda}+\frac{\omega}{\lambda}=\alpha+\frac{\omega}{\lambda}$
(15)

and considering a resonance case where $\delta=\omega_{0}-\omega=0$ with
$\lambda\gg\omega$, $\overline{\alpha}$ becomes very small thus
$\theta=\lambda{t}=\frac{\pi}{2}$, since $\overline{A}_{g0}=1$ determined from
eq. (9). The evolution of this interaction determined by applying the AJC
qubit state transition operation in eq. (13) noting the definition of
$\hat{\overline{I}}_{g}$ and $\hat{\overline{\varepsilon}}_{g}$ [11] in eq.
(11) is of the form

$e^{-i\theta\hat{\overline{\varepsilon}}_{g}}|g,0_{A}\rangle=\cos(\theta)|g,0_{A}\rangle-i\sin(\theta)|e,1_{A}\rangle$
(16)

which reduces to

$|g,0_{A}\rangle\rightarrow-i|e,1_{A}\rangle$ (17)

We observe that the atom interacted with mode A and completed half of the Rabi
oscillation, as a result, it contributed a photon to mode A and evolved to
excited state $|e\rangle$. Now, after the interaction time, it enters mode B
containing a single photon, interacting with the cavity mode as follows

$-ie^{i\theta\hat{\overline{\varepsilon}}_{e}}|e,1_{B}\rangle=-i\cos(\theta)|e,1_{B}\rangle+\sin(\theta)|g,0_{B}\rangle$
(18)

After an interaction with mode B for a time $t_{1}=2t$ such that
$t_{1}=\frac{\pi(\overline{R}_{g0}+\overline{R}_{e1})}{\overline{R}_{g0}\overline{R}_{e1}}$,
the driving field is modulated such that
$\theta=\left(\frac{\overline{R}_{g0}\overline{R}_{e1}}{\overline{R}_{g0}+\overline{R}_{e1}}\right)t=\frac{\pi}{2}$
with $\overline{R}_{g0}=2\lambda{\overline{A}_{g0}}=2\lambda$ since
$\overline{A}_{g0}=1$ and
$\overline{R}_{e1}=2\lambda{\overline{A}_{e1}}=2\lambda$ since
$\overline{A}_{e1}=1$. Therefore, $\lambda{t}=\frac{\pi}{2}$. The form of eq.
(18) results into the evolution

$-i|e,1_{B}\rangle\rightarrow|g,0_{B}\rangle$ (19)

The results in eq. (19) shows that the atom evolves to ground state and
absorbs a photon initially in mode B. Therefore the atom clearly performs a
swapping of the electromagnetic field between the two field modes by
controlled interaction.

When the atom in ground state $|g\rangle$, enters the electromagnetic cavity
containing a single photon in mode A and mode B in vacuum state, the atom and
the field interacts as follows

$e^{-i\theta\hat{\overline{\varepsilon}}_{g}}|g,0_{B}\rangle=\cos(\theta)|g,0_{B}\rangle-\sin(\theta)|e,1_{B}\rangle$
(20)

After an interaction with field mode B for a time
$t=\frac{\pi}{\overline{R}_{g0}}$ equal to half Rabi oscillation time, the
driving field is modulated such that $\theta=\overline{R}_{g0}t=\pi$, with
$\overline{R}_{g0}=2\lambda{\overline{A}}_{g0}=2\lambda$ since
$\overline{A}_{g0}=1$. Therefore $\theta=\lambda{t}=\frac{\pi}{2}$. The form
of eq. (20) results into the evolution

$|g,0_{B}\rangle\rightarrow-|e,1_{B}\rangle$ (21)

The atom then enters mode A containing one photon and interacts as follows

$-e^{i\theta\hat{\overline{\varepsilon}}_{e}}|e,1_{A}\rangle=-\cos(\theta)|e,1_{A}\rangle-i\sin(\theta)|g,0_{A}\rangle$
(22)

After an interaction with the cavity mode for a time $t_{1}=2t$ such that
$t_{1}=\frac{\pi(\overline{R}_{e1}+\overline{R}_{g0})}{\overline{R}_{e1}\overline{R}_{g0}}$
we obtain a driving field modulation
$\theta=\left(\frac{\overline{R}_{e1}\overline{R}_{g0}}{\overline{R}_{e1}+\overline{R}_{g0}}\right)t=\frac{\pi}{2}$,
with $\overline{R}_{e1}=2\lambda{\overline{A}_{e1}}=2\lambda$ since
$\overline{A}_{e1}=1$ and
$\overline{R}_{g0}=2\lambda{\overline{A}_{g0}}=2\lambda$ since
$\overline{A}_{g0}=1$. Therefore $\theta=\lambda{t}=\frac{\pi}{2}$. The form
of eq. (22) results into the evolution

$|e,1_{A}\rangle\rightarrow{i}|g,0_{A}\rangle$ (23)

This shows that the atom evolves to ground state and performs a field swapping
by absorbing a photon in mode A.

When the atom in excited state $|e\rangle$ enters mode A in vacuum state,that
is target qubit $|\mu_{2}\rangle$, the atom propagates as a free wave without
coupling to the field mode in vacuum state $|0\rangle$ [11], leaving the
cavity without altering the state of the cavity-field mode. A similar
observation is made when the atom enters cavity B in vacuum state for the case
of target qubit $|\mu_{1}\rangle$.

From the results obtained, we conclude that the target qubit made up of the
electromagnetic field remains unchanged if the control qubit, that is, the
two-level atom, is initially in the excited state $|e\rangle$, while when the
atom is in ground state $|g\rangle$, the cavity states $|0\rangle$ and
$|1\rangle$ flip. We shall refer to this gate as the AJC C-NOT
$(atom\rightarrow{cavity})$

###  III.1 Probability of success of the C-NOT gate

Success probability for the C-NOT gate is given by

$P_{s}=1-(\sin^{2}(\theta_{A})+\cos^{2}(\theta_{A})\sin^{2}(\theta_{B}))$ (24)

In terms of the Rabi frequency we write eq. (24) as

$P_{s}=1-(\sin^{2}(\overline{R}_{A}\Delta{t_{A}})+\cos^{2}(\overline{R}_{A}\Delta{t_{A}})\sin^{2}(\overline{R}_{B}\Delta{t_{B}}))$
(25)

For the case in which the atom is in the ground state $|g\rangle$ and enters
an electromagnetic cavity with mode A in vacuum state and a single photon of
the field mode B

$\overline{R}_{A}=\overline{R}_{{g}{0}}=2\lambda{\overline{A}_{{g}{0}}}=2\lambda$

$\Delta{t_{A}}=\frac{\pi}{\overline{R}_{A}}=\frac{\pi}{\overline{R}_{{g}{0}}}=\frac{\pi}{2\lambda}$

$\overline{R}_{B}=\overline{R}_{{e}{1}}=2\lambda{\overline{A}_{{e}{1}}}=2\lambda$

$\Delta{t_{B}}=\frac{\pi}{2}\frac{(\overline{R}_{A}+\overline{R}_{B})}{\overline{R}_{A}\overline{R}_{B}}=\frac{\pi}{2}\frac{(\overline{R}_{{g}{0}}+\overline{R}_{{e}{1}})}{\overline{R}_{{g}{0}}\overline{R}_{{e}{1}}}=\frac{\pi}{2\lambda}$
(26)

substituting eq. (26) into eq. (25) we obtain

$P_{s}=1-(\sin^{2}(\pi)+\cos^{2}(\pi)\sin^{2}(\pi))=1$ (27)

a unit probability of success.

When the atom in ground state $|g\rangle$ enters an electromagnetic cavity
containing a single photon in mode A, and mode B in the vacuum state

$\overline{R}_{A}=\overline{R}_{{e}{1}}=2\lambda{\overline{A}_{{e}{1}}}=2\lambda$

$\Delta{t_{A}}=\frac{\pi}{2}\frac{(\overline{R}_{A}+\overline{R}_{B})}{\overline{R}_{A}\overline{R}_{B}}=\frac{\pi}{2}\frac{(\overline{R}_{{e}{1}}+\overline{R}_{{g}{0}})}{\overline{R}_{{e}{1}}\overline{R}_{{g}{0}}}=\frac{\pi}{2\lambda}$

$\overline{R}_{B}=\overline{R}_{{g}{0}}=2\lambda{\overline{A}_{{g}{0}}}=2\lambda$

$\Delta{t_{B}}=\frac{\pi}{\overline{R}_{B}}=\frac{\pi}{\overline{R}_{{g}{0}}}=\frac{\pi}{2\lambda}$
(28)

substituting eq. (28) into eq. (25) we obtain

$P_{s}=1-(\sin^{2}(\pi)+\cos^{2}(\pi)\sin^{2}(\pi))=1$ (29)

a unit probability of success.

We observe that success probabilities depend mainly upon the precise selection
of the interaction times of the two level atom with the successive cavity
modes.

##  IV Hadamard logic gate

To realise Hadamard operation in the AJC interaction, we apply the qubit state
transition operation in eq. (11) and the general form in [10, 11]. In this
respect, we define the Hadamard operation at sum frequency
$\overline{\delta}=4\lambda$ and $n=0$ specified for an initial atomic state
$|g\rangle$ as

$\hat{\overline{\varepsilon}}_{g}=\frac{1}{\overline{A}_{g0}}\left(2\hat{s}_{z}+\hat{a}\hat{s}_{-}+\hat{a}^{\dagger}\hat{s}_{+}\right)\hskip
14.22636pt;\hskip 14.22636pt\overline{A}_{g0}=\sqrt{2}$ (30)

The initial atomic state $|g\rangle$ is rotated to

$|g\rangle\rightarrow\frac{1}{\sqrt{2}}(|e\rangle-|g\rangle)$ (31)

Similarly, the Hadamard operation at sum frequency
$\overline{\delta}=4\lambda$ and $n=1$ specified for an initial atomic state
$|e\rangle$ is defined as [11]

$\hat{\overline{\varepsilon}}_{e}=\frac{1}{\overline{A}_{e1}}\left(2\hat{s}_{z}+\hat{a}\hat{s}_{-}+\hat{a}^{\dagger}\hat{s}_{+}\right)\hskip
14.22636pt;\hskip 14.22636pt\overline{A}_{e1}=\sqrt{2}$ (32)

The initial atomic state $|e\rangle$ is rotated to

$|e\rangle\rightarrow\frac{1}{\sqrt{2}}(|e\rangle+|g\rangle)$ (33)

The Hadamard transformations in Eqs. (31) and (33) realised in the AJC
interaction process (AJC model) agree precisely with the standard definition
in Eq. (3).

##  V  Conclusion

In this paper we have shown how to implement quantum C-NOT and Hadamard gates
in the anti-Jaynes-Cummings interaction mechanism. We obtained ideal unit
probabilities of success due to precise selection of interaction times during
the C-NOT gate operations. We also realised efficient Hadamard operations
through application of respective AJC qubit state transitions.

## Acknowledgment

We thank Maseno University, Department of Physics and Materials Science for
providing a conducive environment to do this work.

## References

  * Nielsen and Chuang [2011] M. A. Nielsen and I. L. Chuang, _Quantum Computation and Quantum Information_ (Cambridge University Press, 2011).
  * Scherer [2019] W. Scherer, _Mathematics of Quantum Computing: An Introduction_ (Springer Nature, 2019).
  * Deutsch [1989] D. E. Deutsch, Proc R. Soc. London A 425, 73 (1989).
  * Raussendorf and Briegel [2001] R. Raussendorf and H. J. Briegel, Phys. Rev. Lett. 86, 5188 (2001).
  * Deutsch and Jozsa [1992] D. Deutsch and R. Jozsa, Proc. R. Soc. London A 439, 553 (1992).
  * Barenco et al. [1995a] A. Barenco, C. H. Bennett, R. Cleve, D. P. DiVincenzo, N. Margolus, P. Shor, T. Sleator, J. A. Smolin, and H. Weinfurter, Phys. Rev. A 52, 3457 (1995a).
  * Shi [2002] Y. Shi, arXiv preprint quant-ph/0205115 (2002).
  * Boykin et al. [2000] P. O. Boykin, T. Mor, M. Pulver, V. Roychowdhury, and F. Vatan, Inf. Process. Lett. 75, 101 (2000).
  * Omolo [2017a] J. A. Omolo, Preprint Research Gate, DOI:10.13140/RG.2.2.30936.80647 (2017a).
  * Omolo [2017b] J. A. Omolo, preprint Research Gate, DOI:10.13140/RG.2.2.11833.67683 (2017b).
  * Omolo [2019] J. A. Omolo, preprint Research Gate, DOI: 10.13140/RG.2.2.27331.96807 (2019).
  * Barenco et al. [1995b] A. Barenco, D. Deutsch, A. Ekert, and R. Jozsa, Phys. Rev. Lett. 74, 4083 (1995b).
  * Knill et al. [2002] E. Knill, R. Laflamme, H. Barnum, D. Dalvit, J. Dziarmaga, J. Gubernatis, L. Gurvits, G. Ortiz, L. Viola, and W. Zurek, arXiv preprint quant-ph/0207171 (2002).
  * Braunstein et al. [1992] S. L. Braunstein, A. Mann, and M. Revzen, Phys. Rev. Lett. 68, 3259 (1992).
  * Domokos et al. [1995] P. Domokos, J.-M. Raimond, M. Brune, and S. Haroche, Phys. Rev. A 52, 3554 (1995).
  * Vitali et al. [2001] D. Vitali, V. Giovannetti, and P. Tombesi, in _Macroscopic Quantum Coherence and Quantum Computing_ (Springer, 2001), pp. 235–244.
  * Saif et al. [2007] F. Saif, R. Ul Islam, and M. Javed, J. Russ. Laser Res. 28, 529 (2007).


††thanks: Corresponding author<EMAIL_ADDRESS>
# Uncertainty Quantification for Virtual Diagnostic of Particle Accelerators

Owen Convery  SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA
Lewis Smith  Department of Computer Science, University of Oxford, UK  Yarin
Gal  Department of Computer Science, University of Oxford, UK  Adi Hanuka
SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA

###### Abstract

Virtual Diagnostic (VD) is a computational tool based on deep learning that
can be used to predict a diagnostic output. VDs are especially useful in
systems where measuring the output is invasive, limited, costly or runs the
risk of altering the output. Given a prediction, it is necessary to relay how
reliable that prediction is, i.e. quantify the uncertainty of the prediction.
In this paper, we use ensemble methods and quantile regression neural networks
to explore different ways of creating and analyzing prediction’s uncertainty
on experimental data from the Linac Coherent Light Source at SLAC National
Lab. We aim to accurately and confidently predict the current profile or
longitudinal phase space images of the electron beam. The ability to make
informed decisions under uncertainty is crucial for reliable deployment of
deep learning tools on safety-critical systems as particle accelerators.

##  I Introduction

Particle accelerators serve a wide variety of applications ranging from
chemistry, physics to biology experiments. Those experiments require increased
accuracy of diagnostics tools to measure the beam properties during its
acceleration, transport and delivery to users. Diagnostics must keep pace with
the advance of extreme beam conditions and the increased experiments’
complexity, which presents challenges to the current state-of-the-art
Craievich  _et al._ (2020); Marchetti  _et al._ (2021). Machine learning (ML)
applications for accelerator diagnostics have been recently shown to improve
beam stability Leemann  _et al._ (2019), identify glitchy accelerator
components Tennant  _et al._ (2020), perform optics corrections and detect
faulty beam position monitors Fol  _et al._ (2020); Arpaia  _et al._ (2021),
and control in real-time by combining adaptive feedback and ML Scheinker  _et
al._ (2018).

Given readily available input data, virtual diagnostic (VD) tools provide a
shot-to-shot non-invasive measurement of the beam in cases where the
diagnostic has limited resolution or limited availability Emma  _et al._
(2018); Hanuka  _et al._ (2021); Sanchez-Gonzalez  _et al._ (2017). VDs have
the potential to be useful in experiment’s design, setup and optimization
while saving valuable operation time. They could also aid in interpreting
experimental results, especially in cases in which current diagnostics cannot
provide necessary information.

Current VD provides predictive models based on training a neural network
mapping between non-invasive diagnostic input to invasive output measurements
(Emma  _et al._ , 2018; Hanuka  _et al._ , 2021; Emma  _et al._ , 2019). This
type of mapping is known as supervised regression. Previous work has
demonstrated VD to predict the electron beam current profile and Longitudinal
Phase Space (LPS) distribution Marx  _et al._ (2018a) along the accelerator
using either scalar controls Emma  _et al._ (2018) or spectral information
Hanuka  _et al._ (2021) as the non-invasive input to the VD. However, while
current VD methods provide mean prediction only, an essential component in
deploying the virtual diagnostic tool is to quantify the confidence in the
prediction, i.e. estimate an interval presenting the uncertainty in
prediction. A metric of uncertainty and its reliable presentation provide a
way of making informed decisions, that becomes crucial in safety-critical
systems such as particle accelerators.

In general, there are two classes of uncertainty: epistemic and aleatoric
uncertainty Kendall and Gal (2017). Aleatoric, or statistical irreducible
uncertainty, is the uncertainty in the data set that arises from experimental
error or inherent measurement noise. Given the same set of inputs, we may
observe slightly different results. Epistemic uncertainty, also known as
systematic reducible uncertainty, is produced when the model’s knowledge is
limited or hindered, e.g. from missing data. We aim to capture both types of
uncertainty in the VD tool and to incorporate them into the final prediction.
While the neural network can only make point predictions about beam
properties, we can use tools from deep learning to better understand the
uncertainty of the predictions.

There are various ways to quantify uncertainty of machine learning models for
supervised regression problems. These include Bayesian model averaging
approaches such as Gaussian processes Rasmussen (2003) and Bayesian neural
networks (BNNs) (MacKay, 1992), which learn distributions over the parameters
of the network. Compared to Gaussian processes, BNNs have the advantage that
they can easily scale to high dimensional inputs. However, BNNs are
computationally expensive and require modifications to the training procedure.
Other non-Bayesian approaches widely used are bootstrapping, ensembling Bickel
and Freedman (1981); Clarke (2003) and quantile regression Torossian  _et al._
(2020); Sriram  _et al._ (2013). Those non-Bayesian approaches are simple to
implement, could be easily parallelized to scale with large amount of data,
and yields high quality predictive uncertainty estimates Lakshminarayanan  _et
al._ (2017); Dietterich (2000). Ensembles are closely related to Bayesian
model averaging, and ensembles of neural networks can be seen as an ad hoc
approximation of Bayesian neural networks or Gaussian processes for specific
architectures and training conditions Gal (2016); Smith and Gal (2018); He
_et al._ (2020). When each ensemble component is trained with random
initialization _and_ a random data subset, the procedure is known as bootstrap
aggregation, or bagging. Bootstrapping is well studied in the Frequentist
statistics literature and is closely related to Bayesian statistics with
certain priors (Efron, 1982; Efron and Tibshirani, 1994; Bickel and Freedman,
1981).

In this paper, we apply deep learning tools to provide confidence intervals
for virtual diagnostic predictions. We compare the predicted uncertainty
learned by multiple ensemble methods as well as quantile regression using
experimental data from the Linac Coherent Light Source (LCLS) at SLAC National
Lab Emma  _et al._ (2010). We evaluate the robustness of the ensemble methods
to provide accurate mean predictions as well as reliable uncertainty
estimation.

The paper is organized as follows: in section II we first present the VD
architecture. Then we discuss metrics to evaluate the reliability of the
prediction to be used online on the particle accelerator. In section III we
demonstrate the methods and compare them on two experimental datasets from
LCLS; 1D current profiles and 2D longitudinal phase space images. The 1D
current profile is a projection of the LPS image. Finally, in section IV we
discuss re-weighing the data sets, and other neural network architectures to
improve the results.

##  II Methods

In this section we first present the Virtual Diagnostic (VD) neural network
architecture and the data sets. We then present various ensemble and quantile
methods to quantify the uncertainty. Last, we present the metrics used to
evaluate the accuracy of the VD’s mean prediction and uncertainty.

###  II.1 VD Architecture and Data Set

High brightness beam linacs typically operate in single-pass, multi-stage
configurations where a high-density electron beam formed in the RF gun is
accelerated and manipulated prior to delivery to users in an experimental
station. An example of such a facility is the LCLS XFEL at SLAC where the
electron beam traverses through an undulator, and emits coherent X-ray pulses.
Typically, longitudinal phase space (LPS) is destructively measured by an
X-band transverse deflecting cavity (XTCAV) Marx  _et al._ (2018b); Craievich
_et al._ (2020) \- as shown in Figure 1.

In this paper we used two experimental data sets: 2D LPS images or 1D current
profiles, measured at the XTCAV, as the outputs, and the corresponding
spectral information, as can be collected by IR spectrometer, as an input. The
1D current profile is a projection of an LPS image. The outputs have been
centered and cropped about their region of interest as in Ref Emma  _et al._
(2018). Examples of the inputs and outputs are shown in Figure 2. Both data
sets contain 4,046 shots and were randomly shuffled and split 80$\%$ for
training and validation, and 20$\%$ for testing. We then normalized all output
profiles and images. The outputs shown in the paper are all normalized.

Figure 1: Schematic of a generic, linac-driven light source.

Reproduced from Hanuka  _et al._ (2021).

The neural network (NN) architecture we used is a fully connected feed-forward
NN composed of three hidden layers (200, 100, 50) with rectified linear unit
activation function. Other network architectures are discussed in section IV.
For training we used a batch size of 32, 500 epochs and the Adam optimizer
with a fixed learning rate of 0.001 Hanuka  _et al._ (2021). Training the NNs
with a Gaussian likelihood, i.e. to minimize the standard Mean Squared Error
(MSE) loss function on a training set, yields symmetric uncertainty intervals.
To introduce asymmetric uncertainty intervals relative to the mean value, we
trained on quantiles of data as described in subsection II.2. The open source
Keras and TensorFlow libraries were used to build and train the NN module
Chollet  _et al._ (2015); Abadi  _et al._ (2015).

(a) Spectra

(b) Current Profiles

(c) LPS for spectrum $\\#$1

(d) LPS for spectrum $\\#$2

Figure 2: Example of VD input of spectrum in (a) and corresponding outputs;
either 1D current profiles in (b), or images of Longitudinal Phase Space (LPS)
— examples of which are shown in (c) and (d) for two different spectra input.

###  II.2 Uncertainty Quantification Methods

Various ways exist to estimate an interval presenting the uncertainty in the
neural network prediction. In what follows, we describe a few ensemble
methods, wherein a collection of various neural networks is combined together.
In this case, the variance of the ensemble is an estimate of the uncertainty.
We then describe quantile regression to predict multiple quantiles of the
data. We also applied the popular MC dropout technique Gal and Ghahramani
(2016). However, we observed a degraded mean prediction and inaccurate
uncertainty bounds compared to ensembling across the entire dataset. Further
improvement could combine ensembles and MC dropout but this is beyond the
scope of this paper.

####  II.2.1 Ensemble methods

Deep ensembles are ensembles of neural networks, each initialized differently
and optimized to converge independently of the others. Ref. Lakshminarayanan
_et al._ (2017) has shown empirically that explicit ensembling can result in
improved uncertainty estimates when they used large neural networks with non-
convex loss surfaces. Deep ensembles are closely related to Bayesian neural
networks and Gaussian processes. For example, Ref. He  _et al._ (2020) has
recently shown that with a small adaptation to their training procedure, deep
ensemble components can be interpreted as draws from a Gaussian process
posterior with a certain kernel.

We investigated three ways of generating this ensemble: using random
initialization of the NN parameters, using a random subset of the data (e.g.
k-fold cross validation), and bagging (e.g. using both). The predicted current
profile for a test shot $\vec{I}_{{\rm
predicted}}=M^{-1}\sum_{m=1}^{M}\vec{i}_{{\rm predicted},m}$ is the mean
prediction of $M$ neural network predictions $\vec{i}_{{\rm predicted},m}$.
The uncertainty is manifested as the standard deviation of the neural network
predictions $\vec{\sigma}=\sqrt{M^{-1}\sum_{m=1}^{M}(\vec{i}_{{\rm
predicted},m}-\vec{I}_{\rm predicted})^{2}}$.

Random parameters initialization. Random initialization of the neural network
parameters (e.g. weights) is a promising approach for improving prediction
accuracy Lee  _et al._ (2015) and uncertainty Lakshminarayanan  _et al._
(2017). It has been recently shown that random initializations explore
entirely diverse solutions in function space Fort  _et al._ (2020) rather than
collapse to the same solution Gal (2016). Popular initializations used in many
applications are the He and Glorot uniform distributions He  _et al._ (2015);
Glorot and Bengio (2010). Those initializations often lead to quick and
reliable convergence during training. However, to create a more diverse
ensemble of neural networks, a random normal initialization has the potential
to reduce the amount of near identical solutions. Here, we initialized all NN
parameters as independent and identically distributed random variables from
$\mathcal{N}(0,0.05)$ when training on the same predetermined
training/validation data split. Glorot Uniform is used for initialization in
all other parts of this paper.

K-Fold cross validation Rather than using the same training/validation split
for each NN in the ensemble, we randomly select a subset of the data. We fixed
the seed for a Glorot uniform distribution for all ensemble components. We
split the training data set into $K$ partitions for $K$ neural networks. Using
a different partition for validating each model would yield a global model
that can generalize better since it has been trained/validated using many
subsets of the same data. This is a more responsible way to validate models
for ensembles Raschka (2018).

Bagging. Bootstrap AGGregatING (bagging) is an ensemble method that trains
each NN on a different ‘bag’ of data. Each bag randomly (with replacement)
contains $n^{\prime}$ out of $n$ possible data points where $n^{\prime}/n$ is
typically $\sim 60\%$. All ensemble components are initialized with a
different Glorot uniform seed, resulting in a NN model with both random sub-
sampling of the data and random initialization. Bagging is a method commonly
used to reduce variance and avoid over-fitting Opitz and Maclin (1999). 3(a)
shows an example of the measured current profile (normalized) $\vec{I}_{\rm
measured}$, its mean prediction $\vec{I}_{\rm predict}$ and the uncertainty
($\pm 2\vec{\sigma}$) of bagging with 16 ensemble components.

(a) Bagging ensemble

(b) Quantile regression

Figure 3: Normalized measured current profile (red) of shot $\\#729$ and its
prediction (blue) using (a) bagging with 16 ensembles and with $\pm 2\sigma$
uncertainty interval, and (b) quantile regression with 90$\%$ uncertainty
interval.

####  II.2.2 Quantile regression

Neural networks can be trained to predict multiple quantiles of the data
Torossian  _et al._ (2020) by using a tilted loss function

${{\rm\mathcal{L}}(\zeta_{t}|\tau)}=\begin{cases}\tau\zeta_{t}&\text{if
$\zeta_{t}\geq 0$}\\\ (\tau-1)\zeta_{t}&\text{if $\zeta_{t}<0$}\\\
\end{cases}$ (1)

where ${{\zeta_{t}}}={I}_{{\rm measured},t}-{I}_{{\rm predicted},t}$ for each
time index $t$. The average loss over the entire shot of length $T$ is
${T^{-1}{\sum}_{t=1}^{T}{\mathcal{L}}(\zeta_{t}|\tau)}$. In practice, a
separate NN is used to predict each quantile. Each is trained using the tilted
loss function but with a different $\tau$ value corresponding to the desired
quantile, 0.5 being the median prediction.

3(b) shows 19 quantiles between 0.05 and 0.95. When plotting multiple quantile
lines at once, we can gain better insight into where the ground truth may lie.
We used the median quantile $(\tau=0.50)$ for calculating the MSE of the
prediction.

###  II.3 Metrics for Evaluation

In order to compare different approaches to quantify uncertainty, we need to
have a way to evaluate the desirable characteristics of the uncertainty. Those
metrics of the predictive uncertainty quality will be used in real-time during
the machine operations as well to indicate the validation of the VD’s
prediction. We used a custom accuracy metric to evaluate the model’s
performance in capturing the ground truth. At each time index ($t$) in the
shot, we check whether the measured ground truth value ($I_{{\rm
measured},t}$) was within an upper $I_{{\rm upper},t}$ and lower $I_{{\rm
lower},t}$ prediction bounds. We then weight it by the measured value to
penalize the near-zero noise while focusing on important features in the
signal. The prediction accuracy of a shot with a length $T$ is defined as:

${{\rm Accuracy}=\frac{\sum_{t=1}^{T}\alpha_{t}\cdot I_{{\rm
measured},t}^{2}}{\sum_{t=1}^{T}I^{2}_{{\rm measured},t}}}$ (2)

where ${\rm\alpha}_{t}=\mathbbm{1}(I_{{\rm lower}{,t}}<I_{{\rm
measured}{,t}}<I_{{\rm upper}{,t}})$. Unless otherwise stated, for a symmetric
loss function we used bounds of $I_{{\rm upper}{,t}},I_{{\rm
lower}{,t}}=I_{{\rm predicted}{,t}}\pm 2\sigma_{t}$ where $\sigma_{t}$ is
predictive standard deviation at time $t$. For the tilted loss function, the
NN predictions using $\tau=0.05$ and $\tau=0.95$ were used for $I_{\rm lower}$
and $I_{\rm upper}$ respectively. We chose to use $\pm 2\sigma$ of the
ensembles and $0.90$ $(0.95-0.05)$ of the quantiles for these calculations
because we wanted to capture $\geq 90\%$ of the uncertainties correlated with
each method.

Since the measured ground truth will not be available when operating the
machine in real-time, in order to to discern a good prediction from a poor
prediction we correlate between the mean squared error between the measured
and predicted values - ${\rm MSE}={T^{-1}\sum_{t=1}^{T}{(I_{{\rm
measured},t}-I_{{\rm predicted},t})}^{2}}$ \- for a vector length of $T$ and
the maximum predicted uncertainty $\sigma_{\rm max}=\max\\{\vec{\sigma}\\}$. A
positive correlation indicates that the model can accurately inform which
predictions have low vs. high uncertainty before seeing the ground truth. If
choosing to retain/reject predictions in practice, a retention curve can be
used to set a proper threshold for deployment; we retain (and predict) shots
with maximum standard deviation per shot smaller than a given threshold.

##  III Results

In this section, we evaluate the methods described in subsection II.2 with the
metrics discussed in subsection II.3 on the 1D current profiles. Next, we
apply those techniques to the 2D LPS image data set. Finally, we discuss
problems that are unique to the images and how they can be relieved.

###  III.1 1D Current Profiles

Before exploring the quantified uncertainty using the 1D current profile data
set, we first found the ensemble method which yielded the best MSE on the mean
prediction. Figure 4 shows the MSE vs. ensemble size for the three ensemble
methods. Bagging has a better MSE than random initialization and performs
slightly worse than the K-Fold ensemble for small ensemble size. Since bagging
with 16 ensemble components yielded the smallest MSE, we evaluated this
model’s performance using the correlation plot and retention curve described
in subsection II.3.

Figure 4: MSE vs. ensemble size for three ensemble methods: Random
initializations (red), K-fold cross validation (green), bagging (blue) with
standard deviation from five runs.  Figure 5: Retention Curve with the
threshold ($\theta$) set by $\sigma_{\rm max}$ to evaluate the Bagging
ensemble method. When a Retention${}_{\theta}[\%]$ of the predictions are
accepted, MSEθ indicates the VD performance on that specific prediction
subset.

(a) Test Shot

(b) Out-of-distribution

(c) Spectra

Figure 6: Comparison of (a) test shot in the same distribution as the training
set, and (b) out-of-distribution (OOD) test shot. The uncertainty intervals
are predicted with an ensemble model of Bagging with 16 components and 90%
quantiles (dashed). The spectra of both shots are compared in (c).

Figure 5 shows retention curve, which can be used to set a rejection threshold
for deployment of the VD. Here, ${\rm MSE}_{\theta}$ is the average MSE for
all test shots with $\sigma_{\rm max}$ lower than a given threshold $\theta$
displayed on the colorbar. The ${\rm Retention}_{\theta}$ on the x-axis
describes the percentage of test shots we retain vs. reject in the process of
choosing a threshold $\theta$. If a VD shot prediction yields a $\sigma_{\rm
max}$ higher than the set $\theta$ value, we ignore the prediction. The
retention curve is helpful for balancing the rejection/threshold trade-off for
a specific use case and is comparable to a Receiver Operating Characteristic
(ROC) curve used for binary classification tasks Gneiting and Vogel (2018). In
our case, we see a desirable positive correlation between ${\rm MSE}_{\theta}$
and ${\rm Retention}_{\theta}$. In deployment, the retention curve can be
interpreted as follows: if the experimental conditions can tolerate shots with
averaged $\rm MSE_{\theta}\leq 2$e-5, by choosing $\theta=1.97$e-2, we can
retain $80\%$ of the shots and yield a low $\rm MSE_{\theta}$ of 1.85e-05.

Over the entire test set, the average MSE of bagging with 16 ensemble
components was 3.65e-05 with an accuracy of 0.452 using $\pm 1\sigma$ and
0.706 using $\pm 2\sigma$ for computing $\alpha$ in Equation 2. The average
MSE for quantile regression method (using $\tau=0.50$) was 4.05e-05. Accuracy
was 0.973 using quantiles 0.05-0.95 and 0.648 using quantiles 0.25-0.75. 6(a)
shows a test shot and the predicted uncertainty from the bagging ensemble
($\pm 2\sigma$) and the $90\%$ of the quantiles to represent epistemic and
aleatoric uncertainty respectively Tagasovska and Lopez-Paz (2019).

In order to evaluate the robustness of the model we use it to predict an out-
of-distribution (OOD) sample. This means that the test shot has a different
probability distribution than the training data set Ovadia  _et al._ (2019).
Therefore, it is expected that the model would manifest higher prediction
uncertainty. In particle accelerators, shifts in the distribution could come,
for example, from a different operation mode or slow drift over time. In those
cases, it is critical to have reliable predictive uncertainty. As an example,
we compare a prediction of a test shot with the same distribution as the
training data set (6(a)), and a prediction for an out-of-distribution input
from a different operation mode of the accelerator (6(b)). The spectrum inputs
to the VD are shown in 6(c). It is evident that the $\pm 2\sigma$ interval
from the ensembling method captures most of the uncertainty in the measured
OOD shot.

This shows that the VD can identify out-of-distribution shots, wherein the
quantiles and ensemble predict abnormally large uncertainty ranges, and these
shots can be rejected. Previous works suggest that we can treat the predictive
standard deviation from ensembles as epistemic uncertainty Lakshminarayanan
_et al._ (2017), and quantile regression as doing maximum likelihood with an
asymmetric Laplace distribution, hence capturing aleatoric uncertainty Tishby
_et al._ (1989). Further development of this technique would allow users to
know whether the predicted uncertainty comes from lack of useful data, or
variation within that training data.

###  III.2 2D LPS Images

We used similar architecture to train ensembles with 2D longitudinal phase
space images. The average MSE of the entire data set is 6.714e-04, and the
accuracy with $\pm 1\sigma$ is 0.291, and with $\pm 2\sigma$ is 0.538. The
accuracy of the image predictions was calculated using Equation 2 evaluated
across an additional dimension:

${{\rm Accuracy}=\frac{\sum_{t,e=1}^{T,E}\alpha_{t,e}\cdot L_{{\rm
measured},t,e}^{2}}{\sum_{t,e=1}^{T,E}L^{2}_{{\rm measured},t,e}}}$ (3)

where ${\rm\alpha}_{t,e}=\mathbbm{1}(|L_{{\rm measured}{,t,e}}|<L_{{\rm
predicted}{,t,e}}\pm 2\sigma_{t,e})s$.

Figure 7 shows two examples of shots in the test set. The panels from left to
right show the measured LPS image, the difference between the measured and
predicted, the accuracy metric (Eq. 2), and the prediction uncertainty
($\sigma$). The blue and red colors on the difference panel indicates the
positive and negative differences of the predicted LPS respectively. The red
mask on the accuracy panel indicates regions where the measured value falls
within the predicted uncertainty.

(a) Translational error - Shot $\\#762$

(b) Shape error - Shot $\\#789$

Figure 7: 2D longitudinal phase space (LPS) measurements with (a)
translational error and (b) shape error examples.

The two examples show predictions with translational error (shot $\\#$762 in
7(a)) and shape error (shot $\\#$789 in 7(b)) with corresponding MSEs of
2.737e-4 and 1.096e-3. We classified a prediction as one with translational
error when the positive and negative differences are shifted with respect to
each other. For the translation error in shot $\\#$762, shifting the mean
prediction up by a single pixel brings a 85.0$\%$ decrease in MSE and a
41.8$\%$ increase in accuracy. In order to determine how to shift the
prediction, the ‘center of mass’ of the prediction and the ground truth are
calculated and then translated to match. In section IV, we further discuss
ways to leverage the spatial connectivity in these images to reduce such
errors.

##  IV Discussion

In this section we discuss two techniques aiming to improve the LPS image
neural network predictions. First we present the concept of a transposed
convolution which can leverage the spatial connectivity in each data sample to
produce better predictions. Then we present a specific “Bottleneck” NN
architecture that was applied and the insights derived from it.

Transposed convolutions. Transposed convolutions are a mathematical operation
that allow for the up-scaling of data. Filters are trained to learn features
which are then projected on to a feature map larger than the input. Transposed
convolutions can be especially useful in regression tasks where the output
involves important features and connectivity Dumoulin and Visin (2018). The 2D
outputs of the LPS image dataset are essentially images rather than vectors of
unrelated measurements. Therefore, it makes sense to model architectures
better fit for image generation which take advantage of the structured nature
of the output. If we train filters to learn specific features in the data set
(tails, curves, bright spots, etc.), the NN model would be better equipped to
handle the data set at hand.

Bottleneck architecture. Bottleneck-shaped architecture is a NN design in
which the number of neurons in the middle layer is smaller than in the other
layers. The number of neurons per layer gradually decrease from the input to
the middle layer, and then gradually increase from the middle layer to the
output. The first part serves as an encoder or compressor which condenses
essential information, and the second part serves as a decoder. The middle
layer allows for the input to be recorded as an $N$ dimensional vector. This
embedding vector can be viewed as a representation of the instance in those
$N$ dimensions, with similar data points having a small distance between them,
and distinct ones being very distant. Bottleneck architectures can reduce
over-fitting by decreasing the system complexity Perin  _et al._ (2020).

Figure 8: t-SNE clustering of the compressed vectors of the bottleneck
architecture for all test shots, with six LPS examples of the outputs shown on
the top. Notably, shots are grouped by shape and maximum current in the latent
space.

The encoder was a dense NN with layers of size 256, 128, and an embedding
dimension of $N=64$. The decoder had transposed convolution layers of 8 2x2
filters, 8 4x4 filters, and 16 16x16 filters. Using this architecture, we
trained a NN which achieved a negligibly improved average MSE of 8.833e-4.
Nevertheless, a strength of the bottleneck architecture is its ability to
condense information about the input and output relation into a single vector.
In a high dimensional space, we expect that similar shots would aggregate
close to each other, and distinctly different shots would be distant from one
another.

In order to visualize this, we used the t-SNE clustering algorithm Arora  _et
al._ (2018) on all the embedding vectors for shots in the test set. After
investigating multiple parameter values, we used a perplexity of $30$, a
learning rate of $200$, and observed that the algorithm converged after $1000$
iterations. Figure 8 shows a visualization of the latent space of these
vectors using two components. The maximum current for each shot is shown in
color bar. The six examples shown illustrate that shots are grouped by shape
and maximum current in the latent space.

##  V Conclusions

In this paper we presented several deep learning approaches to incorporate
uncertainty in the virtual diagnostic (VD) tool. Deep learning models as
neural networks have great predictive power, but they often suffer from over-
fitting and provide over-confident predictions. Here, we compared three
ensemble methods and quantile regression as a way to provide accurate mean
predictions as well as correctly capturing confidence intervals. Specifically,
we considered a VD trained on non-invasive spectral measurements of the
electron beam to predict the 2D longitudinal phase space or a 1D current
profile.

The UQ methods presented in the paper were shown to be robust against out-of-
distribution inputs by providing un-confident predictions on data it was not
trained with. In comparison to a simple dense neural network, it was shown
that a more tailored architecture can be used to exploit information about the
data and offer better explainability during inference. A principled approach
to quantifying uncertainty is crucial for the deployment of the virtual
diagnostic tool, especially for safety-critical systems as particle
accelerators.

## Acknowledgments

This work was supported by the Department of Energy, Laboratory Directed
Research and Development program at SLAC National Accelerator Laboratory,
under contract DE-AC02-76SF00515.

##  Appendix A Re-Weighting the Current Profile Data Set

Upon further analysis of the data set, it became clear that certain current
profile features are under-represented. As a result, the prediction of
profiles with higher peak current presents high MSE. Typically, using models
trained on this unbalanced data would be reasonable assuming nominal usage of
the machine. However, although we have fewer examples of high peak shots, we
would like to train the model to give them equal importance as the lower peak
shots. Therefore, we weighted the output (current profile) such that we gave
instances with lower representation more weight during training, and instances
with over-represented samples less weight. Re-weighting the data not only
greatly reduced the MSE for many shots with high peak shots, but it also
slightly reduced the overall MSE in other cases.

In order to determine the weight for each shot, we made a 100-bin histogram of
$I_{\rm max}$ for all $S=\sum_{i=1}^{100}n_{i}$ shots in the training set.
Denoting $n_{i}$ as the number of shots in the $i$-th bin, the weight of the
shots in that bin is ${w_{i}}=c\cdot(n_{i})^{-0.5}$ with a normalizing
constant ${c}=S\cdot(\sum_{i=1}^{100}n_{i}^{0.5})^{-1}$. We tried other
constants for the weighting power, but found that 0.5 as weighting constant
performed the best. The overall MSE of the test set was 4.601e-5 for the
unweighted NN, and 4.478e-5 for the re-weighted NN (both averaged over 5
runs).

## References

  * Craievich  _et al._ (2020) P. Craievich, M. Bopp, H.-H. Braun, A. Citterio, R. Fortunati, R. Ganter, T. Kleeb, F. Marcellini, M. Pedrozzi, E. Prat, S. Reiche, K. Rolli, R. Sieber, A. Grudiev, W. L. Millar, N. Catalan-Lasheras, G. McMonagle, S. Pitman, V. d. P. Romano, K. T. Szypula, W. Wuensch, B. Marchetti, R. Assmann, F. Christie, B. Conrad, R. D’Arcy, M. Foese, P. G. Caminal, M. Hoffmann, M. Huening, R. Jonas, O. Krebs, S. Lederer, D. Marx, J. Osterhoff, M. Reukauff, H. Schlarb, S. Schreiber, G. Tews, M. Vogt, A. d. Z. Wagner, and S. Wesch, Phys. Rev. Accel. Beams 23, 112001 (2020).
  * Marchetti  _et al._ (2021) B. Marchetti, A. Grudiev, P. Craievich, R. Assmann, H. H. Braun, N. Catalan Lasheras, F. Christie, R. D’Arcy, R. Fortunati, R. Ganter, P. González Caminal, M. Hoffmann, M. Huening, S. M. Jaster-Merz, R. Jonas, F. Marcellini, D. Marx, G. McMonagle, J. Osterhoff, M. Pedrozzi, E. Prat Costa, S. Reiche, M. Reukauff, S. Schreiber, G. Tews, M. Vogt, S. Wesch, and W. Wuensch, Scientific Reports 11, 3560 (2021).
  * Leemann  _et al._ (2019) S. C. Leemann, S. Liu, A. Hexemer, M. A. Marcus, C. N. Melton, H. Nishimura, and C. Sun, Phys. Rev. Lett. 123, 194801 (2019).
  * Tennant  _et al._ (2020) C. Tennant, A. Carpenter, T. Powers, A. Shabalina Solopova, L. Vidyaratne, and K. Iftekharuddin, Phys. Rev. Accel. Beams 23, 114601 (2020).
  * Fol  _et al._ (2020) E. Fol, R. Tomás, J. Coello de Portugal, and G. Franchetti, Phys. Rev. Accel. Beams 23, 102805 (2020).
  * Arpaia  _et al._ (2021) P. Arpaia, G. Azzopardi, F. Blanc, G. Bregliozzi, X. Buffat, L. Coyle, E. Fol, F. Giordano, M. Giovannozzi, T. Pieloni, R. Prevete, S. Redaelli, B. Salvachua, B. Salvant, M. Schenk, M. S. Camillocci, R. Tomás, G. Valentino, F. Van der Veken, and J. Wenninger, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 985, 164652 (2021).
  * Scheinker  _et al._ (2018) A. Scheinker, A. Edelen, D. Bohler, C. Emma, and A. Lutman, Phys. Rev. Lett. 121, 044801 (2018).
  * Emma  _et al._ (2018) C. Emma, A. Edelen, M. J. Hogan, B. O’Shea, G. White, and V. Yakimenko, Phys. Rev. Accel. Beams 21, 112802 (2018).
  * Hanuka  _et al._ (2021) A. Hanuka, C. Emma, T. Maxwell, A. Fisher, B. Jacobson, M. J. Hogan, and Z. Huang, Nature Scientific Reports 11, 2945 (2021).
  * Sanchez-Gonzalez  _et al._ (2017) A. Sanchez-Gonzalez, P. Micaelli, C. Olivier, T. R. Barillot, M. Ilchen, A. A. Lutman, A. Marinelli, T. Maxwell, A. Achner, M. Agåker, N. Berrah, C. Bostedt, J. D. Bozek, J. Buck, P. H. Bucksbaum, S. C. Montero, B. Cooper, J. P. Cryan, M. Dong, R. Feifel, L. J. Frasinski, H. Fukuzawa, A. Galler, G. Hartmann, N. Hartmann, W. Helml, A. S. Johnson, A. Knie, A. O. Lindahl, J. Liu, K. Motomura, M. Mucke, C. O’Grady, J. E. Rubensson, E. R. Simpson, R. J. Squibb, C. Såthe, K. Ueda, M. Vacher, D. J. Walke, V. Zhaunerchyk, R. N. Coffee, and J. P. Marangos, Nature Communications 8, 15461 (2017).
  * Emma  _et al._ (2019) C. Emma, M. D. Alverson, A. Edelen, A. Hanuka, M. Hogan, B. O’Shea, D. Storey, G. White, and V. Yakimenko, Proceedings of the 8th International Beam Instrumentation Conference (2019), 10.18429/JACOW-IBIC2019-THBO01.
  * Marx  _et al._ (2018a) D. Marx, R. Assmann, P. Craievich, U. Dorda, A. Grudiev, and B. Marchetti, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 909, 374 (2018a), 3rd European Advanced Accelerator Concepts workshop (EAAC2017).
  * Kendall and Gal (2017) A. Kendall and Y. Gal, “What uncertainties do we need in bayesian deep learning for computer vision?”  (2017), arXiv:1703.04977 .
  * Rasmussen (2003) C. E. Rasmussen, in  _Summer School on Machine Learning_ (Springer, 2003) pp. 63–71.
  * MacKay (1992) D. J. MacKay,  _Bayesian methods for adaptive models_ , Ph.D. thesis, California Institute of Technology (1992).
  * Bickel and Freedman (1981) P. J. Bickel and D. A. Freedman, The annals of statistics , 1196 (1981).
  * Clarke (2003) B. Clarke, Journal of Machine Learning Research 4, 683 (2003).
  * Torossian  _et al._ (2020) L. Torossian, V. Picheny, R. Faivre, and A. Garivier, “A review on quantile regression for stochastic computer experiments,”  (2020), arXiv:1901.07874 .
  * Sriram  _et al._ (2013) K. Sriram, R. Ramamoorthi, and P. Ghosh, Bayesian Analysis 8, 479 (2013).
  * Lakshminarayanan  _et al._ (2017) B. Lakshminarayanan, A. Pritzel, and C. Blundell, in  _Advances in neural information processing systems_ (2017) pp. 6402–6413.
  * Dietterich (2000) T. G. Dietterich, in  _Multiple Classifier Systems_ (Springer Berlin Heidelberg, Berlin, Heidelberg, 2000) pp. 1–15.
  * Gal (2016) Y. Gal,  _Uncertainty in deep learning_ , Ph.D. thesis, University of Cambridge (2016).
  * Smith and Gal (2018) L. Smith and Y. Gal, UAI (2018).
  * He  _et al._ (2020) B. He, B. Lakshminarayanan, and Y. W. Teh, arXiv:2007.05864 (2020).
  * Efron (1982) B. Efron,  _The jackknife, the bootstrap and other resampling plans_ (SIAM, 1982).
  * Efron and Tibshirani (1994) B. Efron and R. J. Tibshirani,  _An introduction to the bootstrap_ (CRC press, 1994).
  * Emma  _et al._ (2010) P. Emma, R. Akre, J. Arthur, R. Bionta, C. Bostedt, J. Bozek, A. Brachmann, P. Bucksbaum, R. Coffee, F.-J. Decker, Y. Ding, D. Dowell, S. Edstrom, A. Fisher, J. Frisch, S. Gilevich, J. Hastings, G. Hays, P. Hering, Z. Huang, R. Iverson, H. Loos, M. Messerschmidt, A. Miahnahri, S. Moeller, H.-D. Nuhn, G. Pile, D. Ratner, J. Rzepiela, D. Schultz, T. Smith, P. Stefan, H. Tompkins, J. Turner, J. Welch, W. White, J. Wu, G. Yocky, and J. Galayda, Nature Photonics 4, 641 (2010).
  * Marx  _et al._ (2018b) D. Marx, R. Assmann, P. Craievich, U. Dorda, A. Grudiev, and B. Marchetti, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 909, 374 (2018b), 3rd European Advanced Accelerator Concepts workshop (EAAC2017).
  * Chollet  _et al._ (2015) F. Chollet _et al._ , “Keras,” https://keras.io (2015).
  * Abadi  _et al._ (2015) M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng, “TensorFlow: Large-scale machine learning on heterogeneous systems,”  (2015), software available from tensorflow.org.
  * Gal and Ghahramani (2016) Y. Gal and Z. Ghahramani, in  _Proceedings of The 33rd International Conference on Machine Learning_ (2016) pp. 1050–1059.
  * Lee  _et al._ (2015) S. Lee, S. Purushwalkam, M. Cogswell, D. Crandall, and D. Batra, “Why m heads are better than one: Training a diverse ensemble of deep networks,”  (2015), arXiv:1511.06314 .
  * Fort  _et al._ (2020) S. Fort, H. Hu, and B. Lakshminarayanan, “Deep ensembles: A loss landscape perspective,”  (2020), arXiv:1912.02757 .
  * He  _et al._ (2015) K. He, X. Zhang, S. Ren, and J. Sun, in  _Proceedings of the 2015 IEEE International Conference on Computer Vision (ICCV)_, ICCV ’15 (IEEE Computer Society, USA, 2015) p. 1026–1034.
  * Glorot and Bengio (2010) X. Glorot and Y. Bengio, in  _Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics_, Proceedings of Machine Learning Research, Vol. 9, edited by Y. W. Teh and M. Titterington (PMLR, Chia Laguna Resort, Sardinia, Italy, 2010) pp. 249–256.
  * Raschka (2018) S. Raschka, “Model evaluation, model selection, and algorithm selection in machine learning,”  (2018), arXiv:1811.12808 .
  * Opitz and Maclin (1999) D. Opitz and R. Maclin, Journal of Artificial Intelligence Research 11, 169–198 (1999).
  * Gneiting and Vogel (2018) T. Gneiting and P. Vogel, “Receiver operating characteristic (roc) curves,”  (2018), arXiv:1809.04808 .
  * Tagasovska and Lopez-Paz (2019) N. Tagasovska and D. Lopez-Paz, in  _Advances in Neural Information Processing Systems_ (2019).
  * Ovadia  _et al._ (2019) Y. Ovadia, E. Fertig, J. Ren, Z. Nado, D. Sculley, S. Nowozin, J. V. Dillon, B. Lakshminarayanan, and J. Snoek, Advances in neural information processing systems (2019), 1906.02530 .
  * Tishby  _et al._ (1989) Tishby, Levin, and Solla, in  _International 1989 Joint Conference on Neural Networks_ (1989) pp. 403–409 vol.2.
  * Dumoulin and Visin (2018) V. Dumoulin and F. Visin, “A guide to convolution arithmetic for deep learning,”  (2018), arXiv:1603.07285 .
  * Perin  _et al._ (2020) G. Perin, I. Buhan, and S. Picek, IACR Cryptol. ePrint Arch. 2020, 58 (2020).
  * Arora  _et al._ (2018) S. Arora, W. Hu, and P. K. Kothari, in  _Proceedings of the 31st Conference On Learning Theory_, Proceedings of Machine Learning Research, Vol. 75, edited by S. Bubeck, V. Perchet, and P. Rigollet (PMLR, 2018) pp. 1455–1462.


# Weakly Supervised Semantic Segmentation of Remote Sensing Images for Tree
Species Classification based on Explanation Methods

###### Abstract

The collection of a high number of pixel-based labeled training samples for
tree species identification is time consuming and costly in operational
forestry applications. To address this problem, in this paper we investigate
the effectiveness of explanation methods for deep neural networks in
performing weakly supervised semantic segmentation using only image-level
labels. Specifically, we consider four methods: i) class activation maps
(CAM); ii) gradient-based CAM; iii) pixel correlation module; and iv) self-
enhancing maps (SEM). We compare these methods with each other using both
quantitative and qualitative measures of their segmentation accuracy, as well
as their computational requirements. Experimental results obtained on an
aerial image archive show that: i) considered explanation techniques are
highly relevant for the identification of tree species with weak supervision;
and ii) the SEM outperforms the other considered methods. The code for this
paper is publicly available at https://git.tu-berlin.de/rsim/rs_wsss.

Index Terms—  Tree species mapping, weakly supervised learning, semantic
segmentation, explanation methods, remote sensing.

##  1 Introduction

Accurate identification of tree species by the analysis of remote sensing (RS)
images is important for various forestry applications. By reducing cost
intensive on-site surveys, it significantly supports public authorities,
conservation agencies and private owners in forest mapping and management.
Tree species mapping can be achieved by using semantic segmentation methods,
which aim to predict pixel-wise classification results on RS images. Deep
learning (DL) based semantic segmentation methods have recently seen a rise in
popularity in the context of tree species classification [1]. Most DL models
require a high amount of labeled samples to optimize all parameters and reach
a high performance of tree species classification. The labeling of samples can
be achieved based on: 1) in situ ground surveys; 2) the expert interpretation
of color composites (image photo-interpretation); or 3) hybrid solutions where
both photo-interpretation and ground surveys are exploited [2]. Collection of
a sufficient number of high quality pixel-level labels associated to tree
species can be time consuming, complex and costly.

To address this problem, in this paper we focus our attention on weakly
supervised semantic segmentation (WSSS) for tree species segmentation, which
rely on weak supervision (i.e., image-level labels). The use of image-level
labels can significantly reduce the annotation cost and effort in forestry
applications. However, such supervision indicates only the existence of
certain tree species assigned to images without their exact pixel based
location information (which is essential for obtaining tree species maps). To
obtain accurate tree species maps at the pixel-level, in this paper we
investigate the effectiveness of explainable neural networks in the context of
WSSS. Explanation methods are capable of generating explanations in the form
of pixel-level heatmaps that can be highly relevant for providing a pixel-
level tree species map from a DL model trained using image level labels. In
this paper, we consider four explanation methods: i) class activation maps
(CAM) [3]; ii) Gradient-based CAM (GradCAM) [4]; iii) pixel correlation module
(PCM) [5]; and iv) self-enhancing maps (SEM) [6]. All considered methods have
been experimentally compared in terms of their: 1) capability to perform
accurate tree species mapping; 2) model complexity; and 3) segmentation time.
This work is the first, to the best of our knowledge, that explores the
potential of pixel-level explanation heatmaps for tree species mapping in a
weak supervision framework.

##  2 Tree Species Mapping from Image Level Labels

Let $\mathcal{X}=\\{\textbf{X}_{1},\cdots,\textbf{X}_{M}\\}$ denote an archive
of $M$ RS images acquired over forestry areas, where the $m$-th image within
the archive is represented by $\textbf{X}_{m}$. We assume that a training set
$\mathcal{X}_{T}\subset\mathcal{X}$ is available where each image
$\textbf{X}_{n}$ within the set is associated with a multi-label vector
$\textbf{y}_{n}=[y_{n,1},\cdots,y_{n,S}]$, where $S$ is the total number of
classes, indicating which tree species are present. Here, $y_{n,s}$ takes the
value of 1 if the class $s$ is present in $\textbf{X}_{n}$, and 0 otherwise.
The set of all multi-label vectors corresponding to $\mathcal{X}_{T}$ is thus
denoted as $\mathcal{Y}_{T}$.

To obtain tree-species maps by using image-level labels (and thus weak
supervision), we consider a convolutional neural network (CNN) which learns a
mapping $\hat{\textbf{y}}_{n}=U(\textbf{X}_{n})$, where $U(\cdot)$ represents
the CNN. However, the predicted output vector $\hat{\textbf{y}}_{n}$ by the
direct application of a CNN does not provide any spatial information regarding
class locations or extents, which are required for tree species mapping. To
obtain the map at pixel-level, we investigate CAM, GradCAM, PCM, and SEM in
the context of WSSS and exploit the class specific heatmaps (which show the
probability of a given class at each pixel) derived by these methods.
Specifically, feature maps $\textbf{F}\in\mathbb{R}^{C\times H\times W}$ ($C$,
$H$ and $W$ represent the channels, height, and width of the feature maps,
respectively) from the final layer of the trained model $U(\cdot)$
corresponding to an image $\textbf{X}_{n}$ are exploited. The importance of
each feature map with respect to each class is derived and a weighted linear
combination of the feature maps form the class specific heatmaps. Each class
specific heatmap has its values normalized between 0-1, thus acting as pixel
probabilities. Heatmaps are thus used to obtain the semantic segmentation by
assigning class labels to pixels where the class has the highest probability.
Specific details on how the heatmaps are derived from the four methods are
given in the following sections.

###  2.1 Class Activation Maps

CAMs aim to highlight the image regions which are used by the CNN to classify
an image to a given class. This is achieved by establishing a relationship
between the set of feature maps F and a given class. To this end, a 1x1
convolution layer can be applied after F, which takes F as input and outputs
$S$ feature maps, generating a set of CAMs ($\textbf{A}\in\mathbb{R}^{S\times
H\times W}$) [3]. Thus, for a given class $s$, we obtain its CAM as such:

$\textbf{A}_{s}=\sum_{c=1}^{C}w^{c}_{s}\textbf{F}^{c}$ (1)

where $w^{c}_{s}$ is the weight of importance for the $c^{th}$ feature map
with respect to the $s^{th}$ class.

The CAM is then averaged and used as the final class prediction scores (image
level prediction):

$\hat{y}_{n,s}=\frac{1}{HW}\sum_{i=1}^{H}\sum_{j=1}^{W}\textbf{A}_{s,i,j}$ (2)

where $i$ and $j$ are the row and column indices of $\textbf{A}_{s}$,
respectively.

Given this formulation, CAMs can be obtained in an end-to-end manner without
any post-processing steps.

###  2.2 Gradient Based CAM (GradCAM)

GradCAM [4] is a reformulation of CAM. Instead of using learned weights from a
1x1 convolution to determine the importance of each feature map, GradCAM uses
the average gradient associated with each feature map in F. Specifically, we
first feed an image $\textbf{X}_{n}$ through the network in order to obtain
the vector of predicted class scores $\hat{\textbf{y}}_{n}$. Then, we perform
back-propagation with respect to a single class prediction $\hat{y}_{n,s}$ in
order to obtain the gradients for the feature maps F with respect to
$\hat{y}_{n,s}$. Finally, we take the mean of the gradients in order to obtain
the importance of feature map $\textbf{F}^{c}$ for class $s$:

$w^{c}_{s}=\frac{1}{HW}\sum_{i=0}^{W}\sum_{j=0}^{H}\frac{\delta\hat{y}_{n,s}}{\delta\textbf{F}^{c}_{i,j}}$
(3)

We can thus obtain the GradCAM for a given class by replacing $w^{c}_{s}$ in
(1) with the definition given in (3). One distinct advantage of GradCAM is
that it can be applied to any network architecture as it gets the importance
weights of each feature map using gradient scores. However, the WSSS cannot be
obtained in an end-to-end manner, as it requires the back-propagation post-
processing.

We would like to emphasize that both CAM and GradCAM have been recently
applied for land cover mapping in RS. For details, we refer the reader to [7,
8].

###  2.3 Pixel Correlation Module (PCM)

The PCM [5] aims to enhance CAM in order to improve the segmentation. PCM
avoids the need for post-processing as the module is integrated into the
network itself, thus allowing for end-to-end processing. The module consists
of three parallel 1x1 convolutional layers which takes feature maps from
different levels of $U(\cdot)$, along with the original image, as input in
order to embed them into latent space. The embedded outputs are then
concatenated and passed through one additional 1x1 convolution which gives us
a matrix of feature maps $\textbf{K}\in\mathbb{R}^{C\times H\times W}$ which
correspond to an image $\textbf{X}_{n}$ and contain a combination of high and
low level features. From this matrix, attention scores ($attn_{q,r}$) are
obtained which provide the similarity between the features at two spatial
locations (e.g. $q$ and $r$) within K through a self-attention mechanism as:

$attn_{q,r}=\frac{\textbf{K}_{q}^{T}\textbf{K}_{r}}{||\textbf{K}_{q}||\cdot||\textbf{K}_{r}||}$
(4)

where $q$ and $r$ each represent a pixel location (e.g. $q=(0,0)$, $r=(0,1)$).

For a given class $s$, the PCM value at pixel location $q$
($\textbf{P}_{s,q}$) is obtained by taking the weighted sum of all values in
$\textbf{A}_{s}$, where the weight values are the cosine similarity estimated
between the feature vector at $q$ and all other feature vectors across the
spatial dimension, as in (4). This is be formulated as:

$\textbf{P}_{s,q}=\sum_{\forall r}attn_{q,r}\times\textbf{A}_{s,r}$ (5)

As opposed to [5] where PCM is trained via supervision from a Siamese network
and equivariant cross regularization loss, our work trains the module as a
separate branch. To this end, we placed PCM after the feature extraction
backbone and trained it using class predictions made from the output of the
module. Thus, when training the model which includes PCM, we use two binary
cross entropy loss terms, one from the classification head, and one from the
PCM module.

###  2.4 Self-Enhancement Maps (SEM)

SEM [6] also aims to enhance the CAM. SEM works on the principle that feature
vectors located within target object regions have higher similarity than those
within different class regions. Here, input image $\textbf{X}_{n}$ is passed
through the network $U(\cdot)$ in order to obtain A and the feature maps F.
Given a class specific $\textbf{A}_{s}$, $E$ seed coordinates corresponding to
locations with the largest values are chosen. The feature vectors
corresponding to the locations of the $E$ seed points are then extracted from
F along the channel dimension. The cosine similarity defined in (4) is then
calculated between the $E$ seed feature vectors and all other feature vectors
in F. This results in $E$ similarity maps, where the value for a given pixel
$r$ within the $e^{th}$ similarity map is the cosine similarity score (4)
between feature at $r$ and the $e^{th}$ feature. The final output for the
given class is then obtained by taking the maximum value at each pixel across
the $E$ similarity maps.

Image

Reference Map

CAM [3]

GradCAM [4]

PCM [5]

SEM [6]

a)

b)

c)

Fig. 1: Examples of images, their reference and segmentation maps obtained by
using different explanation methods in the framework of WSSS.

For each of the methods from the previous four sections, we obtain the final
tree species map by applying the argmax function across the outputs of the
respective method for classes predicted above a threshold $\tau$ at the image
level.

##  3 Experimental Results

Experiments were carried out on a dataset acquired from 2012 - 2020 across the
German federal state of Lower Saxony. The dataset contains a total of 19,995
aerial images, each of which consists of RGB and near infrared bands with a
spatial resolution of 0.2m. The dataset was randomly divided into training
(70%), validation (15%), and test (15%) sets. Dominant tree species within
stand polygons which overlapped an image patch were used to attain the image-
level labels. To evaluate the WSSS results, we rasterized the polygons
corresponding to the test set in order to obtain pixel-wise labels. A total of
five classes were extracted from the polygons: Pine (Pinus spp.), Spruce
(Picea spp.), Beech (Fagus spp.), Oak (Quercus spp.), and Cleared. Cleared
represents any open areas without tree crown cover (e.g. meadows, clear-cuts,
water bodies, etc.). For the experiments we used DeepLabv3+ [9]. Given that
the architecture was designed for semantic segmentation, we replace the
segmentation head with a multi-label classification head in order to obtain
predictions at the image level. For evaluation of the WSSS performance we
utilized $F_{1}$, and also examined the number of model parameters (in
millions) and the segmentation time per image. All results were obtained using
$\tau=0.5$. The number of seeds used in SEM was chosen empirically based on
the performance across the validation image set, resulting in a seed value of
$E=10$.

Table 1: Comparison of the considered explanation methods in terms of $F_{1}$ scores, the number (in millions) of parameters (# Param) and the semantic segmentation time (Seg. Time) per image (in millisecond).  Metric | CAM[3] | GradCAM[4] | PCM[5] | SEM[6]
---|---|---|---|---
$F_{1}$ (%) | 69.88 | 77.18 | 76.70 | 79.17
# Param | 11.65 | 11.65 | 11.69 | 11.65
Seg. Time | 50.4 | 402.9 | 54.9 | 54.5

From Table 1 we can see that SEM provides the highest $F_{1}$ score, whereas
CAM leads to the lowest $F_{1}$ score. Looking at the qualitative results in
Fig. 1, we see that CAM is able to correctly localize class regions, but
oversegments the dominant class within the image. The performance of CAM is
greatly improved when SEM is applied (Fig. 1a, b). Predictions made by SEM
often follow the contours of tree stands (Fig. 1a, b). GradCAM often made
predictions visually similar to SEM (Fig. 1a). However, GradCAM also predicted
incorrect regions (Fig. 1b), leading to a degradation in performance. PCM was
able to correctly localize class regions, but tended to oversegment the
dominant class within the image, leading to misclassified regions (Fig. 1a,
b). This was likely due to how the module was trained, as the original
implementation involved a Siamese network with a cross-regularization loss.
This suggests that PCM may not be a viable addition to standard CNN
architectures. One area where all methods struggled was when two classes of
the same type (e.g. broadleaf or coniferous) were predicted within an image
(Fig. 1c). Here, the classes are not well separated within the learned feature
space of the classifier, leading to misclassifications. In regards to
segmentation time, CAM had the lowest time, while GradCAM had the highest
time. This was due to the need for post-processing step involving a forward
and backward pass through the network for each predicted class within an image
with GradCAM. SEM also required a post-processing step, but managed to
maintain a low segmentation time. This is because SEM only performs self-
attention for a given number of seeds, thus making it less computational than
typical self-attention mechanisms. The main drawback with SEM is that more
computational time must be invested into tuning the seed hyperparameter, which
can become time consuming with larger datasets. With respect to model
parameters, only PCM increased model parameters. However, the additional
parameters did not achieve any benefits in segmentation performance with
respect to GradCAM and SEM.

##  4 Conclusion

In this paper, we have studied four explanation methods in the framework of
WSSS to obtain pixel-level tree species maps using training samples annotated
by image-level labels. In detail, we have investigated: i) CAM; ii) GradCAM;
iii) PCM; and iv) SEM based on their semantic segmentation performance, number
of model parameters, and semantic segmentation time. The theoretical and
experimental analysis show that for tree species segmentation problems, SEM
can be chosen as it: i) yields the highest segmentation accuracy; ii) provides
the lowest model complexity; and iii) requires a low segmentation time. As a
future work, we plan to assess the ability of explanation methods in providing
pixel-level pseudo-labels for training a semantic segmentation model.

##  5 Acknowledgments

This work is funded by the European Research Council (ERC) through the
ERC-2017-STG BigEarth Project under Grant 759764 and by the German Federal
Ministry of Education and Research through TreeSatAI project under Grant
01lS20014A. We also thank the State Forest of Lower Saxony for providing the
aerial images and stand level data.

## References

  * [1] F. Schiefer, T. Kattenborn, A. Frick, J. Frey, P. Schall, B. Koch, and S. Schmidtlein,  “Mapping forest tree species in high resolution uav-based rgb-imagery by means of convolutional neural networks,”  ISPRS Journal of Photogrammetry and Remote Sensing, vol. 170, pp. 205–215, 2020.
  * [2] B. Demir, C. Persello, and L. Bruzzone,  “Batch-mode active-learning methods for the interactive classification of remote sensing images,”  IEEE Transactions on Geoscience and Remote Sensing, vol. 49, no. 3, pp. 1014–1031, 2011.
  * [3] X. Zhang, Y. Wei, J. Feng, Y. Yang, and T. S. Huang,  “Adversarial complementary learning for weakly supervised object localization,”  in IEEE Conference on Computer Vision and Pattern Recognition, 2018, pp. 1325–1334.
  * [4] R. R. Selvaraju, M. Cogswell, A. Das, R. Vedantam, D. Parikh, and D. Batra,  “Grad-cam: Visual explanations from deep networks via gradient-based localization,”  in IEEE International Conference on Computer Vision, 2017, pp. 618–626.
  * [5] Y. Wang, J. Zhang, M. Kan, S. Shan, and X. Chen,  “Self-supervised equivariant attention mechanism for weakly supervised semantic segmentation,”  in IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2020, pp. 12275–12284.
  * [6] X. Zhang, Z. Wei, Y. Yang, and F. Wu,  “Rethinking localization map: Towards accurate object perception with self-enhancement maps,”  arXiv preprint arXiv:2006.05220, 2020.
  * [7] Y. Li, W. Chen, Y. Zhang, C. Tao, R. Xiao, and Y. Tan,  “Accurate cloud detection in high-resolution remote sensing imagery by weakly supervised deep learning,”  Remote Sensing of Environment, vol. 250, pp. 112045, 2020.
  * [8] L. Chan, M. S. Hosseini, and K. N. Plataniotis,  “A comprehensive analysis of weakly-supervised semantic segmentation in different image domains,”  International Journal of Computer Vision, vol. 129, no. 2, pp. 361–384, 2021.
  * [9] L. C. Chen, Y. Zhu, G. Papandreou, F. Schroff, and H. Adam,  “Encoder-decoder with atrous separable convolution for semantic image segmentation,”  in Proceedings of the European Conference on Computer Vision, 2018, pp. 801–818.