# Estimating Seroprevalence of SARS-CoV-2 in Ohio: A Bayesian Multilevel
Poststratification Approach with Multiple Diagnostic Tests

David Kline  Center for Biostatistics, Department of Biomedical Informatics,
The Ohio State University, Columbus, Ohio  Zehang Li  Department of
Statistics, University of California Santa Cruz, Santa Cruz, California  Yue
Chu  Department of Sociology, The Ohio State University, Columbus, Ohio  Jon
Wakefield  Department of Statistics and Department of Biostatistics,
University of Washington, Seattle, Washington  William C. Miller  Division of
Epidemiology, College of Public Health, The Ohio State University, Columbus,
Ohio  Abigail Norris Turner  Division of Infectious Diseases, College of
Medicine, The Ohio State University, Columbus, Ohio  Samuel J. Clark
Department of Sociology, The Ohio State University, Columbus, Ohio

###### Abstract

Globally the SARS-CoV-2 coronavirus has infected more than 59 million people
and killed more than 1.39 million. Designing and monitoring interventions to
slow and stop the spread of the virus require knowledge of how many people
have been and are currently infected, where they live, and how they interact.
The first step is an accurate assessment of the population prevalence of past
infections. There are very few population-representative prevalence studies of
the SARS-CoV-2 coronavirus, and only two American states – Indiana and
Connecticut – have reported probability-based sample surveys that characterize
state-wide prevalence of the SARS-CoV-2 coronavirus. One of the difficulties
is the fact that the tests to detect and characterize SARS-CoV-2 coronavirus
antibodies are new, not well characterized, and generally function poorly.
During July, 2020, a survey representing all adults in the State of Ohio in
the United States collected biomarkers and information on protective behavior
related to the SARS-CoV-2 coronavirus. Several features of the survey make it
difficult to estimate past prevalence: 1) a low response rate, 2) very low
number of positive cases, and 3) the fact that multiple, poor quality
serological tests were used to detect SARS-CoV-2 antibodies. We describe a new
Bayesian approach for analyzing the biomarker data that simultaneously
addresses these challenges and characterizes the potential effect of selective
response. The model does not require survey sample weights, accounts for
multiple, imperfect antibody test results, and characterizes uncertainty
related to the sample survey and the multiple, imperfect, potentially
correlated tests.

##  1 Introduction

Slowing or stopping the spread of a new virus for which a vaccine does not
exist starts with two key pieces of information. One, what fraction of the
population has been infected and is thereby potentially less or not
susceptible to future infection, and two, what fraction of the population is
currently infected and potentially able to infect others. Together with a
basic understanding of the infection process, this information roughly
characterizes the potential for the epidemic to grow. This is invaluable to
public health officials and policy makers who have the responsibility to
manage the epidemic and protect the public.

As of this writing (late November, 2020), the global SARS-CoV-2 coronavirus
(hereafter “CV19”) pandemic has infected more than 59 million people and
killed more than 1.39 million [1]. Basic epidemiological information to
describe the pandemic is scarce because the virus is new and the pandemic
exploded rapidly. In its place are a wide variety of indicators based on
convenient, mostly non-representative, or indirectly related data – counts of
all-cause deaths (e.g. [2, 3, 4]), facility-based testing results for
symptomatic patients, non-representative samples, and in many situations,
results from inadequately characterized tests that perform poorly [5].
Franceschi and colleagues [6] identify 37 CV19 prevalence studies from
nineteen countries. Most present results from non-representative, otherwise
special, or very small study populations. Just fourteen represent large enough
populations to be of policy interest – national or state-level – and use a
probability-based sample from a credible sampling frame to produce results
that could represent the population of interest – Asia: one [7], Europe: seven
[8, 9, 10, 11, 12, 13, 14], North America: two [15, 16], and South America:
four [17, 18, 19, 20]. The two in North America are the State of Indiana [15]
and the State of Connecticut [16] in the United States.

Conducting population-representative biomarker surveys is difficult –
particularly in the United States. Good sampling frames exist in a variety of
forms (tax rolls, telephone numbers, etc.), but recruiting willing respondents
is exceptionally difficult and likely affected by selection relative to the
outcome of interest. Both of the studies in the United States had low response
rates for the full interview with valid CV19 test results – Indiana: 23.4% and
Connecticut: 7.8%. Further complicating analysis, there were few positive
tests among those who did respond – Indiana: $47/3,658=1.3\%$ for the PCR test
of current infection, and $38/3,658=1.0\%$ for antibody test of ever infected;
and Connecticut: $23/567=4.1\%$ for the antibody test of ever infected. Both
studies described concern that the non-responding participants were likely to
be at higher risk of infection with CV19. Finally, like all CV19 immunology
investigations to date, both studies struggled with poor quality antibody
tests whose unfavorable performance characteristics were not well understood;
see [7] for an overview of these issues.

Statistical analysis of data like these is difficult. First, the low response
rate requires extensive recalibration of the sampling weights, and in the
worst case, there may be sampling units with no respondents at all. Second,
the very small number of positive cases pushes the asymptotic (large sample)
assumptions of frequentist methods to their limits and can break them. Third,
the imperfect and poorly characterized antibody tests potentially add a lot of
uncertainty that must be reflected in the results, particularly in low
prevalence settings [21]. Fourth, when results from multiple tests with
different performance characteristics are combined, the joint result must be
accurately described and its uncertainty propagated to the final estimate of
prevalence – importantly, including the possibility that results from
individual tests are correlated. Finally, if there is selection on the
outcome, then the effect of this must be understood. In our review of the
literature we did not find an existing method that addresses all of these
challenges in a unified way.

Here, we describe an analytical approach developed to produce estimates of
past infection with CV19 using data from a probability-based household survey
representing adults in the State of Ohio in the United States. Like the CV19
prevalence studies in Indiana and Connecticut, the Ohio survey had a low
response rate, few positive cases, and the possibility of selective response.
Additionally, the Ohio survey used multiple imperfect antibody tests for the
same antibodies, resulting in the need to quantify uncertainty in the joint
result and account for possible dependence among results.

To overcome these challenges, we weave together two well-established modeling
frameworks into a single coherent approach. We utilize the literature on
modeling multiple imperfect diagnostic tests through the use of a Bayesian
latent class model (e.g. [22, 23]). This enables us to combine information
across tests to infer the true latent infection status of a participant while
incorporating uncertainty about the characteristics of the tests. We use the
latent infection status to generate model-based estimates of the population
prevalence using multilevel regression and poststratification [24, 21]. These
approaches are integrated into a single Bayesian model that allows for the
full propagation of uncertainty, exact inferences, and the ability to specify
informative priors using external information. By doing so we produce
estimates that reflect all available information and uncertainty.

##  2 Methods

The purpose of this study is to estimate the prevalence of past CV19
infections in the State of Ohio using three separate antibody tests given to
randomly selected adult participants. We know that each antibody test is
imperfect, and there is no gold standard for detecting prior CV19 infection.
Prevalence estimates based on a single imperfect test are always biased, but
particularly in the case of CV19 infection rates, which are low [21]. To
mitigate that bias and incorporate variability due to error in the testing
results, we will take a Bayesian latent class approach for modeling multiple
diagnostic tests. Our approach will be based on combining a fixed effects
framework for modeling conditional dependence across multiple diagnostic tests
[22, 23] with a model-based analysis using multilevel regression and
poststratification [24] to acknowledge the complex design aspects of the
survey.

###  2.1 Survey Design

The survey was designed to provide policymakers with a “quick” overall
snapshot of the prevalence of prior infection at the state level. The survey
sampling scheme was designed as a stratified two-stage cluster sample. Strata
were defined by 8 administrative regions used by the state. Within each
region, 30 census tracts were randomly selected with probability proportional
to size (PPS) based on total population size. Then within a selected tract, 5
households were randomly selected and one adult (at least 18 years of age)
within each household was randomly selected to participate in the study. Thus,
the planned target sample size was 1200 participants. The study was conducted
from July 9-28, 2020.

Each participant in the study provided biological samples which would be put
through a series of diagnostic tests. In this paper, we will focus on prior
infection as determined by the presence of SARS-CoV-2 antibodies. Since
antibody tests for CV19 are new to the market and of varying quality, we
processed participant samples using 3 different antibody tests. Specifically,
we used the Abbott IgG, Liaison IgG, and Epitope IgM tests.

###  2.2 Model

For each participant in the study, indexed by $i=1,\ldots,n$, let
$\textbf{T}_{i}=(T_{i1},T_{i2},T_{i3})$ be indicators of a positive test
result of the Abbott IgG, Liaison IgG, and Epitope IgM, respectively. Let
$D_{i}$ be the unobserved indicator of whether participant $i$ had a prior
infection of CV19. This latent indicator of prior infection is our primary
outcome of interest. Analysis methods for multiple diagnostic tests without a
gold standard hinge on assumptions related to conditional independence [23].
We will assume that $(T_{i1},T_{i2})$ and $T_{i3}$ are independent given the
true infection status. This implies that conditional on the true presence of
prior infection, we assume tests for the same antibody are dependent and tests
for different antibodies are independent. Based on the underlying design of
the tests and what they target, we believe these assumptions are reasonable
for this specific application.

Table 1: Mapping from the three binary test results to the multinomial response vector and corresponding probabilities. $Y_{i}$ | $T_{i1}$ | $T_{i2}$ | $T_{i3}$ | $\textbf{p}_{i}$
---|---|---|---|---
1 | 1 | 1 | 1 | $p_{i1}$
2 | 1 | 1 | 0 | $p_{i2}$
3 | 1 | 0 | 1 | $p_{i3}$
4 | 0 | 1 | 1 | $p_{i4}$
5 | 1 | 0 | 0 | $p_{i5}$
6 | 0 | 1 | 0 | $p_{i6}$
7 | 0 | 0 | 1 | $p_{i7}$
8 | 0 | 0 | 0 | $p_{i8}$

Given the assumptions stated above, we can consider this a problem with two
conditionally independent sets of tests: the two IgG tests and the IgM test.
Thus, we can decompose the joint probability as

$\displaystyle\Pr(T_{i1},T_{i2},T_{i3}|D_{i})=\Pr(T_{i1},T_{i2}|D_{i})\Pr(T_{i3}|D_{i}).$

By doing so, each conditional probability on the right hand side can be
estimated following a fixed effects approach [22]. Since each test result is
binary, this leads to 8 potential combinations of results, shown in Table 1,
and suggests the following distribution:

$\displaystyle Y_{i}|D_{i}\sim\text{Multinomial}(1,\textbf{p}_{i})$

where $Y_{i}$ is an indicator of participant $i$’s result pattern and
$\textbf{p}_{i}=\textbf{p}_{i}(D_{i})$ is a vector of length 8 where each
element is the probability of a result pattern.

To construct the probability vector $\textbf{p}_{i}$, we first need to define
the conditional probabilities within each antibody. Let $S_{j}$ be the
sensitivity and $C_{j}$ be the specificity of test $j$. Since we have two IgG
tests ($j=1,2$), we allow for their results to be correlated. Let $R^{1}_{12}$
be the covariance between the results of test $1$ and $2$ given the infection
status is positive and $R^{0}_{12}$ be the covariance when the infection
status is negative. Then for the joint probabilities of the IgG test results
[22], we have

$\displaystyle\Pr(T_{i1}=t_{i1},T_{i2}=t_{i2}|D_{i}=1)=S_{1}^{t_{i1}}S_{2}^{t_{i2}}(1-S_{1})^{1-t_{i1}}(1-S_{2})^{1-t_{i2}}+(-1)^{t_{i1}+t_{i2}}R^{1}_{12}$
(1)
$\displaystyle\Pr(T_{i1}=t_{i1},T_{i2}=t_{i2}|D_{i}=0)=C_{1}^{1-t_{i1}}C_{2}^{1-t_{i2}}(1-C_{1})^{t_{i1}}(1-C_{2})^{t_{i2}}+(-1)^{t_{i1}+t_{i2}}R^{0}_{12}.$

Since we only have one IgM test ($j=3$), we have

$\displaystyle\Pr(T_{i3}=t_{i3}|D_{i}=d_{i})=d_{i}S_{3}^{t_{i3}}(1-S_{3})^{1-t_{i3}}+(1-d_{i})C_{3}^{1-t_{i3}}(1-C_{3})^{t_{i3}}.$
(2)

We then use these conditional probabilities to construct the probabilities in
the vector $\textbf{p}_{i}$ in the multinomial distribution above. Those eight
probabilities can calculated with the following general equation for each test
result pattern

$\displaystyle\Pr(T_{i1}=t_{i1},T_{i2}=t_{i2},T_{i3}=t_{i3}|D_{i}=d_{i})=\Pr(t_{i1},t_{i2}|d_{i})\Pr(t_{i3}|d_{i}).$

These probabilities are the individual elements of the vector $\textbf{p}_{i}$
which is a function of $D_{i},\textbf{S},\textbf{C},\textbf{R}$ where
$\textbf{S}=(S_{1},S_{2},S_{3})$, $\textbf{C}=(C_{1},C_{2},C_{3})$ and
$\textbf{R}=(R^{1}_{12},R^{0}_{12})$.

The latent true infection status is the primary process of interest as we are
interested in estimating the prevalence of prior CV19 infections. We assume

$\displaystyle D_{i}|\pi_{i}\sim\text{Bernoulli}(\pi_{i})$

where $\pi_{i}$ is the probability of prior infection for participant $i$,
which will be determined by region, strata, and census tract. Specifically, we
let $r[i]$, $s[i]$, and $t[i]$ refer to the region, strata, and census tract
for participant $i$ and assume

$\displaystyle\text{logit}(\pi_{i})=\alpha_{r[i]}+\textbf{X}_{s[i]}\boldsymbol{\beta}^{s}+\textbf{X}_{t[i]}\boldsymbol{\beta}^{t}+b_{t[i]}$
(3)

where $\alpha_{r}$ is a region-specific random intercept, $\textbf{X}_{s}$ is
a vector of stratum indicators with fixed effects vector
$\boldsymbol{\beta}^{s}$, $\textbf{X}_{t}$ is a vector of census tract
covariates with fixed effects vector $\boldsymbol{\beta}^{t}$, and $b_{t}$ is
a random effect for census tract. In this analysis, $\textbf{X}_{s}$ includes
indicators of participant age group and sex and $\textbf{X}_{t}$ is a scalar
and corresponds to the log of the total census tract population. We code the
groups using sum to 0 contrasts, and the log population was standardized to
have 0 mean and standard deviation of 1 across the entire state. We assume
$\alpha_{r}\stackrel{{\scriptstyle iid}}{{\sim}}N(\alpha,\sigma^{2})$. Since
the fixed effects are coded to sum to 0, $\alpha$ reflects the overall mean on
the logit scale and $\alpha_{r}$ can be interpreted as regional means. We also
assume $b_{t}\stackrel{{\scriptstyle iid}}{{\sim}}N(0,\tau^{2})$, which
accounts for correlation between individuals in the same census tract.

At this point, we link the diagnostic test model to a multilevel regression
and poststratification approach [24, 21] for estimating population prevalence.
In Equation [3], we specified a multilevel logistic regression model for the
probability of prior infection. Since true infection was rare in our sample,
we were unable to fit a saturated model. Instead, we chose to use a
hierarchical model with fixed strata effects. By doing so, we assume that
while regions may have different probabilities of infection, the relative
ordering of the strata will be the same across regions (i.e., there is no
interaction between region and strata). Effects for region and census tract
population are included to account for the characteristics of the underlying
survey design (stratification and PPS sampling), making the selection process
ignorable [25].

To obtain the population prevalence, we calculate

$\displaystyle\pi=\frac{\sum_{r}\sum_{s}\sum_{t}\pi_{rst}P_{rst}}{\sum_{r}\sum_{s}\sum_{t}P_{rst}}$

where $\pi_{rst}$ is the prevalence and $P_{rst}$ is the adult population in
stratum $s$ in Census tract $t$ in region $r$. The prevalence contribution for
region $r$, strata $s$ and tract $t$ is

$\displaystyle\pi_{rst}=\mbox{expit}(\alpha_{r}+\textbf{X}_{s}\boldsymbol{\beta}^{s}+\textbf{X}_{t}\boldsymbol{\beta}^{t}+b_{t})$

where $b_{t}\stackrel{{\scriptstyle iid}}{{\sim}}N(0,\tau^{2})$.

Since we are fitting the model in the Bayesian paradigm, we must specify prior
distributions on all unknown parameters. The main reason we chose a fixed
effects model for the test results was because it is directly determined by
test sensitivity and specificity. This allows us to transparently incorporate
prior information using the validation data on the package insert (as of
October 2020) for each test. Based on this information, we let

$\displaystyle S_{1}\sim\text{Beta}(109,13)\hskip
36.135ptC_{1}\sim\text{Beta}(1066,4)$ $\displaystyle
S_{2}\sim\text{Beta}(96,39)\hskip 36.135ptC_{2}\sim\text{Beta}(1074,16)$
$\displaystyle S_{3}\sim\text{Beta}(9,11)\hskip
36.135ptC_{3}\sim\text{Beta}(54,0.1).$

The remaining parameters for the diagnostic testing part of the model are the
covariance parameters. We enforce necessary constraints on these parameters by
specifying independent uniform prior distributions for each parameter
restricted to its allowable range [22]. Assuming only positive dependence
between tests, we have

$\displaystyle
R^{1}_{12}\sim\text{Uniform}(0,\text{min}(S_{1},S_{2})-S_{1}S_{2})$
$\displaystyle
R^{0}_{12}\sim\text{Uniform}(0,\text{min}(C_{1},C_{2})-C_{1}C_{2}).$

For the multilevel regression part of the model, we assume each element of
$\boldsymbol{\beta}$ is independently normally distributed with 0 mean and
variance of 9. We assume $\tau$ and $\sigma$ have independent uniform
distributions on $(0,5)$. We assume $\alpha\sim N(\text{logit}(0.03),1)$ to
reflect prior belief that the prevalence of prior infection is around 3% and
puts 95% of the probability between 0.4% and 18.0%.

By putting everything together, we have the full model

$\displaystyle
p(\mathbf{S},\mathbf{C},\mathbf{R},\mathbf{D},\boldsymbol{\alpha},\boldsymbol{\beta},\mathbf{b},\tau^{2},\alpha,\sigma^{2}|\mathbf{T})$
$\displaystyle\propto\prod_{i=1}^{n}\Pr(\mathbf{T}_{i}|\mathbf{S},\mathbf{C},\mathbf{R},D_{i})\times
p(\mathbf{R}|\mathbf{S},\mathbf{C})\times p(\mathbf{S},\mathbf{C})$
$\displaystyle\times\prod_{i=1}^{n}\Pr(D_{i}|\boldsymbol{\alpha},\boldsymbol{\beta},\mathbf{b})\times
p(\mathbf{b}|\tau^{2})\times p(\boldsymbol{\alpha}|\alpha,\sigma^{2})\times
p({\boldsymbol{\beta}})p(\alpha)p(\tau^{2})p(\sigma^{2})$

where $\textbf{T}=(\textbf{T}_{1},\ldots,\textbf{T}_{n})$,
$\textbf{S}=(S_{1},S_{2},S_{3})$, $\textbf{C}=(C_{1},C_{2},C_{3})$,
$\textbf{R}=(R^{1}_{12},R^{0}_{12})$, $\textbf{D}=(D_{1},\ldots,D_{n})$,
$\boldsymbol{\beta}=(\boldsymbol{\beta}^{s},\boldsymbol{\beta}^{t})$,
$\boldsymbol{\alpha}$ is the vector of region-specific intercepts, and b is
the vector of census tract random effects. To compute the posterior
distribution, the model was fit using a Markov Chain Monte Carlo (MCMC)
algorithm implemented in R [26] using NIMBLE [27]. The algorithm was run for
500,000 iterations, discarding the first 250,000 as burn-in, and thinning the
remaining iterations by keeping every 20th draw. Convergence was assessed by
visually inspecting trace plots. Posterior distributions are summarized by the
posterior mean and 95% highest posterior density (HPD) credible interval. Code
is available online at: https://github.com/sinafala/bayes-prevalence.

###  2.3 Missing Test Results

Some participants did not have results for all 3 tests considered. However,
since our primary interest is in the latent infection status of the
participant, we can still gain some information about this from the test
results that are available. Test result pattern probabilities would still be
calculated as above corresponding to the test or pair of tests available for
the participant. For example, if a participant only had the IgG test results,
there would be 4 patterns of results that would be defined by the
probabilities in Equation [1]. If an IgG and IgM test are observed, the 4
probabilities would be computed from products of

$\displaystyle\Pr(T_{ij}=t_{ij}|D_{i}=d_{i})=d_{i}S_{j}^{t_{ij}}(1-S_{j})^{1-t_{ij}}+(1-d_{i})C_{j}^{1-t_{ij}}(1-C_{j})^{t_{ij}}.$
(4)

A single available test result would simply lead to a Bernoulli distribution
of whether the result was positive ($t_{ij}=1$) with the probability
calculated as in Equation [4] with the appropriate parameters for the observed
test.

###  2.4 Sensitivity Analyses

For the analysis of any survey, it is important to account for the potential
of non-response bias. We note that our primary analysis is valid assuming a
non-response mechanism that is ignorable given stratum. This is likely a
strong assumption and so we will consider sensitivity analyses under several
assumptions of potential non-ignorable scenarios. To do so, assume

$\displaystyle\pi_{rst}=\pi^{R}_{rst}(1-p^{N}_{rst})+\pi^{N}_{rst}p^{N}_{rst}$
(5)

where $\pi^{R}_{rst}$ and $\pi^{N}_{rst}$ are the prevalence among the
responders and non-responders and $p^{N}_{rst}$ is the probability of non-
response. Note that when ignorability is satisfied, we have
$\pi_{rst}=\pi^{R}_{rst}=\pi^{N}_{rst}$, which results in the estimates from
the primary analysis. Due to the survey design, we can only obtain household
non-response rates by region, and so, we assume $p^{N}_{rst}=p^{N}_{r}$ for
all strata and tracts. In order to carry out a sensitivity analysis, we will
assume the prevalence among non-responders is:

$\displaystyle\pi^{N}_{rst}=\lambda\pi^{R}_{rst}$ (6)

where $\lambda$ is the prevalence ratio in non-responders compared to
responders and $\pi^{R}_{rst}$ is the prevalence estimated from the logistic
regression described above. We will vary $\lambda$ to explore different
scenarios and assess the sensitivity of our estimates to changes in the
prevalence in non-responders.

##  3 Results

Table 2: Descriptive statistics of participants with at least one antibody test result who were included in the analysis of prior infection prevalence. Variable | Count | Proportion
---|---|---
Age 18-44 | 175 | 0.262
Age 45-64 | 239 | 0.358
Age 65 and over | 253 | 0.379
Male | 275 | 0.412
Female | 392 | 0.588
Central | 82 | 0.123
East Central | 90 | 0.135
Northeast | 75 | 0.112
Northwest | 156 | 0.234
Southeast | 77 | 0.115
Southeast Central | 56 | 0.084
Southwest | 65 | 0.097
West Central | 66 | 0.099

A total of 727 adults participated in the survey. To be included in the
analysis, participants had to have at least one antibody test result and have
age and sex recorded. Of the 727 participants, 667 (92%) were included in the
analysis. Characteristics of those included in the analysis are shown in Table
2. We observe that our included participants tend to be older and female. We
also observe that 23.4% of included participants were from the Northwest
region of Ohio.

Our primary goal is to estimate the statewide prevalence of prior infection
with CV19. Based on our model, the posterior mean prevalence is 1.3% with a
95% credible interval of (0.2%, 2.7%). This corresponds to approximately
117,000 adults. As noted in Section 2, the prevalence estimates are based on
the latent infection status estimates inferred from the antibody test results.
In Figure 1, we show the posterior mean probability of prior infection for all
667 participants included in the analysis. This illustrates that there were
very few participants in the study that were estimated to actually have prior
infection. In fact, only 17 (2.5%) participants had a posterior probability of
past infection greater than 1%. As expected, those with the highest estimated
probabilities were those for whom there was agreement across the tests. In
general, there was little observed agreement across the 3 diagnostic tests. Of
the 39 participants with at least one positive result, only 3 (7.7%) had a
positive result on more than one test. For parameters that had informative
prior distributions, Figures S1 and S2 show prior and posterior densities to
illustrate learning based on the observed data.

(a) Posterior Probability of Infection

(b) Posterior Probability of Infection >1%

Figure 1: Posterior probability of past infection of CV19 for all participants
and limited to those participants with a probability of greater than 1%.

We now turn to the sensitivity analysis that accounted for potential non-
ignorable non-response. The household-level non-response rates by region are
shown in Table 3. We believe that non-responders were most likely to have a
higher prevalence of past infection than responders. In Figure 2, we show that
if we assume the prevalence of past infection in non-responders was three
times that of responders, the upper bound of the 95% credible interval is at a
prevalence of 7%. Thus, across the range of reasonable scenarios that we
considered, we observe rates of past infection that do not dramatically differ
from the estimates in our primary analysis. In addition, all of the estimates
in the sensitivity analysis would be aligned with a similar public health and
policy response as they are not compatible with having reached herd immunity.

Table 3: Household-level non-response rates by region. Region | Rate
---|---
Central | 0.783
East Central | 0.841
Northeast | 0.840
Northwest | 0.832
Southeast | 0.752
Southeast Central | 0.790
Southwest | 0.810
West Central | 0.818
Figure 2: Posterior mean and 95% highest posterior density credible intervals
for the prevalence of past infection under several scenarios of non-ignorable
non-response. The prevalence ratio is $\lambda$ in Equation [6].

##  4 Discussion

In this paper, we present an approach to coherently integrate multiple
imperfect diagnostic tests and a model-based analysis. By using this approach,
we are able to estimate the past prevalence of CV19 in the State of Ohio in
the United States while appropriately accounting for uncertainty in multiple
antibody tests and leveraging the strengths of a designed survey. Through the
Bayesian paradigm, we are also able to incorporate external information
through informative prior distributions and provide exact inferences for the
prevalence estimates. Our approach provides policy makers with the best
possible estimates of prior infection given the survey results and our current
knowledge of the quality of the antibody tests.

Through our study, we estimate the magnitude of prior CV19 infection to be
low, roughly 117,000 adults. However, this still reflects a higher burden of
infection than the reported number of cases – approximately 95,300 after
removing deaths through the end of July 2020. Our sensitivity analysis shows
that even with a prevalence ratio of 3 between non-responders and responders,
the upper bound of the credible interval is a prevalence of about 7%. Thus,
even in hypothetical scenarios with large selective response, we conclude that
the vast majority of Ohio adults still remain susceptible to CV19 infection.
This implies that the state must remain vigilant and continue to deploy non-
pharmaceutical interventions like masks and social distancing to limit the
spread of CV19.

Methodologically, we developed a coherent statistical framework for analyzing
seroprevalence surveys with multiple imperfect diagnostic tests. Typically,
one might rely on a single test or the creation of deterministic rules that
combine test results to assess positive cases. In contrast, we utilized
existing methodology for combining the results of multiple diagnostic tests to
fully incorporate information on infection status from each test. We then
connected the estimated infection status to a multilevel regression and
poststratification approach for generating population-based estimates. Through
our fully Bayesian, model-based analysis, we are able to appropriately
propagate uncertainty and provide exact inferences while synthesizing the
information from each type of test administered. Our approach allows us to
make full use of the data collected by the survey while also accounting for
its quality.

Our study has significant strengths: a representative random sample; employing
multiple diagnostic tests; and conducting a fully Bayesian analysis. There are
also limitations. The survey was designed to estimate CV19 prevalence at the
level of the whole state and was not designed for estimation within any
subgroups or smaller geographies. A future study should consider oversampling
high risk subgroups or regions to enable inference specific to those important
populations. Non-response was a major issue that was addressed in the field
using all available resources. We do not have detailed information on non-
responders that could have been incorporated into a more sophisticated model
for missing data. However, based on our sensitivity analysis, we do not
believe the substantive policy implications of the study would change based on
a more complicated missing data model. Finally, the informative prior
distributions for antibody test characteristics are based on validation data
presented in the test inserts – we are unable to verify the rigor and quality
of those validation studies.

In conclusion, we have developed a new statistical analysis framework for
analyzing seroprevalence survey data with few positive cases and results from
multiple imperfect diagnostic tests. Since estimates of prior and current
prevalence are relevant to policy, seroprevalence studies of CV19 are becoming
increasingly important and more frequently conducted. Our methodological
approach is a critical component to ensuring that serology data are analyzed
in a way that is consistent both with the design of the survey and with the
inherent limitations in the accuracy of antibody tests. This enables policy
makers to have access to the best available estimates that also fully and
honestly account for all of the sources of uncertainty that contribute to the
quantification of CV19 infection.

## References

  * [1] Ensheng Dong, Hongru Du, and Lauren Gardner.  An interactive web-based dashboard to track COVID-19 in real time.  The Lancet Infectious Diseases, 20(5):533–534, 2020.
  * [2] National Center for Health Statistics and others.  Excess deaths associated with COVID-19.  Center for Disease Control and Prevention, Washington, DC, 2020\.
  * [3] Daniel M Weinberger, Jenny Chen, Ted Cohen, Forrest W Crawford, Farzad Mostashari, Don Olson, Virginia E Pitzer, Nicholas G Reich, Marcus Russi, Lone Simonsen, et al.  Estimation of excess deaths associated with the COVID-19 pandemic in the United States, March to May 2020.  Journal of the American Medical Association Internal Medicine, 180(10):1336–1344, 2020.
  * [4] Lauren M Rossen, Amy M Branum, Farida B Ahmad, Paul Sutton, and Robert N Anderson.  Excess deaths associated with COVID-19, by age and race and ethnicity—United States, January 26–October 3, 2020.  Morbidity and Mortality Weekly Report, 69(42):1522, 2020.
  * [5] Mayara Lisboa Bastos, Gamuchirai Tavaziva, Syed Kunal Abidi, Jonathon R Campbell, Louis-Patrick Haraoui, James C Johnston, Zhiyi Lan, Stephanie Law, Emily MacLean, Anete Trajman, et al.  Diagnostic accuracy of serological tests for COVID-19: systematic review and meta-analysis.  British Medical Journal, 370, 2020.
  * [6] Vinicius Bonetti Franceschi, Andressa Schneiders Santos, Andressa Barreto Glaeser, Janini Cristina Paiz, Gabriel Dickin Caldana, Carem Luana Machado Lessa, Amanda de Menezes Mayer, Julia Goncalves Kuchle, Paulo Ricardo Gazzola Zen, Alvaro Vigo, et al.  Population-based prevalence surveys during the COVID-19 pandemic: a systematic review.  medRxiv, 2020.
  * [7] Maryam Shakiba, Seyed Saeed Hashemi Nazari, Fardin Mehrabian, Seyed Mahmoud Rezvani, Zahra Ghasempour, and Abtin Heidarzadeh.  Seroprevalence of COVID-19 virus infection in Guilan province, Iran.  medRxiv, 2020.
  * [8] Helen Ward, Christina J Atchison, Matthew Whitaker, Kylie EC Ainslie, Joshua Elliot, Lucy C Okell, Rozlyn Redd, Deborah Ashby, Christl A Donnelly, Wendy Barclay, et al.  Antibody prevalence for SARS-CoV-2 in England following first peak of the pandemic: REACT2 study in 100,000 adults.  medRxiv, 2020.
  * [9] Steven Riley, Kylie EC Ainslie, Oliver Eales, Benjamin Jeffrey, Caroline E Walters, Christina J Atchison, Peter J Diggle, Deborah Ashby, Christl A Donnelly, Graham Cooke, et al.  Community prevalence of SARS-CoV-2 virus in England during May 2020: REACT study.  medRxiv, 2020.
  * [10] P Maver Vodičar, A Oštrbenk Valenčak, Blaž Zupan, T Avšič Županc, Slavko Kurdija, Miša Korva, Miroslav Petrovec, Janez Demšar, Nataša Knap, Erik Štrumbelj, et al.  Low prevalence of active COVID-19 in Slovenia: a nationwide population study of a probability-based sample.  Clinical Microbiology and Infection, 2020.
  * [11] Béla Merkely, Attila J Szabó, Annamária Kosztin, Ervin Berényi, Andor Sebestyén, Csaba Lengyel, Gergő Merkely, Júlia Karády, István Várkonyi, Csaba Papp, et al.  Novel coronavirus epidemic in the Hungarian population, a cross-sectional nationwide survey to support the exit policy in Hungary.  GeroScience, 42(4):1063–1074, 2020.
  * [12] Marina Pollán, Beatriz Pérez-Gómez, Roberto Pastor-Barriuso, Jesús Oteo, Miguel A Hernán, Mayte Pérez-Olmeda, Jose L Sanmartín, Aurora Fernández-García, Israel Cruz, Nerea Fernández de Larrea, et al.  Prevalence of SARS-CoV-2 in Spain (ENE-COVID): a nationwide, population-based seroepidemiological study.  The Lancet, 396(10250):535–544, 2020.
  * [13] Chantal J Snoeck, Michel Vaillant, Tamir Abdelrahman, Venkata P Satagopam, Jonathan D Turner, Katy Beaumont, Clarissa PC Gomes, Joelle Veronique Fritz, Valerie E Schröder, Anne Kaysen, et al.  Prevalence of SARS-CoV-2 infection in the Luxembourgish population: the CON-VINCE study.  medRxiv, 2020.
  * [14] Daniel F Gudbjartsson, Agnar Helgason, Hakon Jonsson, Olafur T Magnusson, Pall Melsted, Gudmundur L Norddahl, Jona Saemundsdottir, Asgeir Sigurdsson, Patrick Sulem, Arna B Agustsdottir, et al.  Spread of SARS-CoV-2 in the Icelandic population.  New England Journal of Medicine, 2020.
  * [15] Nir Menachemi, Constantin T Yiannoutsos, Brian E Dixon, Thomas J Duszynski, William F Fadel, Kara K Wools-Kaloustian, Nadia Unruh Needleman, Kristina Box, Virginia Caine, Connor Norwood, et al.  Population point prevalence of SARS-CoV-2 infection based on a statewide random sample—Indiana, April 25–29, 2020.  Morbidity and Mortality Weekly Report, 69(29):960, 2020.
  * [16] Shiwani Mahajan, Rajesh Srinivasan, Carrie A Redlich, Sara K Huston, Kelly M Anastasio, Lisa Cashman, Dorothy S Massey, Andrew Dugan, Dan Witters, Jenny Marlar, et al.  Seroprevalence of SARS-CoV-2-specific IgG antibodies among adults living in Connecticut: Post-infection prevalence (PIP) study.  The American Journal of Medicine, 2020.
  * [17] Cristiana Costa Gomes, Crispim Cerutti, Eliana Zandonade, Ethel Leonor Noia Maciel, Filomena Euridice Carvalho de Alencar, Gilton Luiz Almada, Orlei Amaral Cardoso, Pablo Medeiros Jabor, Raphael Lubiana Zanotti, Tania Queiroz Reuter, et al.  A population-based study of the prevalence of COVID-19 infection in Espírito Santo, Brazil: methodology and results of the first stage.  medRxiv, 2020.
  * [18] Antônio Augusto Moura da Silva, Lídio Gonçalves Lima Neto, Conceição de Maria Pedrozo, Silva de Azevedo, Léa Márcia Melo da Costa, Maylla Luana Barbosa Martins Bragança, Allan Kardec Duailibe Barros Filho, Bernardo Bastos Wittlin, Bruno Feres de Souza, Bruno Luciano Carneiro Alves de Oliveira, et al.  Population-based seroprevalence of SARS-CoV-2 is more than halfway through the herd immunity threshold in the State of Maranhao, Brazil.  medRxiv, 2020.
  * [19] Mariângela F Silveira, Aluísio JD Barros, Bernardo L Horta, Lúcia C Pellanda, Gabriel D Victora, Odir A Dellagostin, Claudio J Struchiner, Marcelo N Burattini, Andréia RM Valim, Evelise M Berlezi, et al.  Population-based surveys of antibodies against SARS-CoV-2 in Southern Brazil.  Nature Medicine, 26(8):1196–1199, 2020.
  * [20] Pedro Hallal, Fernando Hartwig, Bernardo Horta, Gabriel D Victora, Mariangela Silveira, Claudio Struchiner, Luis Paulo Vidaletti, Nelson Neumann, Lucia C Pellanda, Odir A Dellagostin, et al.  Remarkable variability in SARS-CoV-2 antibodies across Brazilian regions: nationwide serological household survey in 27 states.  medRxiv, 2020.
  * [21] Andrew Gelman and Bob Carpenter.  Bayesian analysis of tests with unknown specificity and sensitivity.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 69(5):1269–1283, 2020.
  * [22] Nandini Dendukuri and Lawrence Joseph.  Bayesian approaches to modeling the conditional dependence between multiple diagnostic tests.  Biometrics, 57(1):158–167, 2001.
  * [23] Zhuoyu Wang, Nandini Dendukuri, Heather J. Zar, and Lawrence Joseph.  Modeling conditional dependence among multiple diagnostic tests.  Statistics in Medicine, 36(30):4843–4859, 2017.
  * [24] Andrew Gelman.  Struggles with survey weighting and regression modeling.  Statist. Sci., 22(2):153–164, 05 2007.
  * [25] Susanna Makela, Yajuan Si, and Andrew Gelman.  Bayesian inference under cluster sampling with probability proportional to size.  Statistics in Medicine, 37(26):3849–3868, 2018.
  * [26] R Core Team.  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria, 2019.
  * [27] P. de Valpine, D. Turek, C.J. Paciorek, C. Anderson-Bergman, D. Temple Lang, and R. Bodik.  Programming with models: writing statistical algorithms for general model structures with NIMBLE.  Journal of Computational and Graphical Statistics, 26:403–417, 2017\.


# Hybrid Quantum-Classical Photonic Neural Networks

Tristan Austin<EMAIL_ADDRESS>Centre for Nanophotonics, Department of
Physics, Engineering Physics, and Astronomy, Queen’s University.  Simon
Bilodeau  Department of Electrical Engineering, Princeton University  Andrew
Hayman  Centre for Nanophotonics, Department of Physics, Engineering Physics,
and Astronomy, Queen’s University.  Nir Rotenberg  Centre for Nanophotonics,
Department of Physics, Engineering Physics, and Astronomy, Queen’s University.
Bhavin J. Shastri<EMAIL_ADDRESS>Centre for Nanophotonics, Department of
Physics, Engineering Physics, and Astronomy, Queen’s University.  Smith
Engineering, Department of Electrical Engineering, Queen’s University

###### Abstract

Neuromorphic (brain-inspired) photonics leverages photonic chips to accelerate
artificial intelligence [1], offering high-speed and energy efficient
solutions for use in RF communication [2], tensor processing [3], and data
classification [4, 5]. However, the limited physical size of integrated
photonic hardware constrains network complexity and computational capacity. In
light of recent advances in photonic quantum technology [6], it is natural to
utilize quantum exponential speedup to scale photonic neural network
capabilities. Here we show a combination of classical network layers with
trainable continuous variable quantum circuits yields hybrid networks with
improved trainability and accuracy. On a classification task, hybrid networks
achieve the same performance when benchmarked against fully classical networks
that are twice the size. When the bit precision of the optimized networks is
reduced through added noise, the hybrid networks still achieve greater
accuracy when evaluated at state of the art bit precision. These hybrid
quantum classical networks demonstrate a unique route to improve computational
capacity of integrated photonic neural networks without increasing the
physical network size.

††preprint: APS/123-QED

##  I Introduction

Neuromorphic photonics is a promising accelerator for artificial intelligence
and brain-inspired computing [1], capable of quickly diagonalizing matrices
[7], separating mixed RF signals [2], and classifying spoken vowels and
handwritten digits [5, 8]. Neuromorphic networks consist of layers of
interconnected neurons, as sketched in Fig. 1, and, like the brain, the
complexity of these networks is largely determined by the number of
connections between these neurons. On neuromorphic photonic platforms, such as
silicon [9], silicon nitride [10], and lithium niobate [11], these neurons are
implemented by circuitry comprised of photonic resonators, waveguides,
modulators, and detectors [12, 13], resulting in high bandwidth, low loss, and
ultralow latency networks [14, 1]. Yet the relatively large size of these
optical components, particularly in comparison with their electronic
counterparts, limits the size of neuromorphic photonic networks and hence
their complexity and power.

Here, we explore a new route to increasing the network complexity, namely by
replacing layers of classical neurons with a quantum neural network. More
specifically, we envision using a continuous variable (CV) photonic network
[15, 16, 17], of the kind in which quantum advantage was recently demonstrated
[6, 18], which is built on the same platforms, using the same devices, as
typical neuromorphic photonic devices. Using an exemplary classification task,
we show that such hybrid quantum-classical neuromorphic networks can
outperform fully classical networks both in terms of ease of training and
overall accuracy. This advantage is more significant for smaller networks with
$\lesssim 350$ weights (for reference, a state-of-the-art 3-layer $6\times
6\times 6\times 6$ integrated photonic network has 98 weights [19]). Finally,
we show that within this range, hybrid neuromorphic networks are sufficiently
robust to noise, reinforcing that hybrid quantum-classical photonics may play
an important role in unlocking the full potential of brain-inspired hardware.

Figure 1: Summary of hybrid and classical neural network architectures. (a) A
one layer two qumode hybrid neural network, network A, with input, encoding,
CV neural network [15], and output layers. Data is inputted into the classical
input layer, outputs from the inputs layer are then inputted into the encoding
gates [20, 21]. A single CV neural net layer is labelled as
$\mathcal{L}^{(h)}$. This specific network has 120 parameters. The gates in
the quantum neural network are described in Section II and Supplementary Fig.
LABEL:fig:CV_gate_progression and Table LABEL:fig:gates. (b) An equivalent all
classical network to the hybrid network, network B. Here, the hidden layer,
$W^{(h)}$ is a classical layer with two neurons. In the lower right is a
mathematical description of operations happening at each layer of the network.
This network has 124 parameters. (c) Feature 1 and Feature 6 from the
validation set of the classification dataset. The shaded region in each
subplot surrounds all the samples of the corresponding class. (d) The number
of parameters in both the hybrid and classical networks as a function of the
number modes. The inset shows the scaling of the input, hidden, and output
layers for the two qumode hybrid network.

##  II Hybrid Continuous Variable Neural Networks

Figure 2: Demonstration of hybrid and classical networks learning to classify
the synthetically generated dataset. (a) The maximum output activations for an
untrained 120 parameter hybrid network, network A. The x-axis is the class the
network predicted while the colour corresponds to the true class. The inset
shows the maximum output activations for all samples predicted as class three.
(b) The maximum output activation for an untrained 124 parameter network B.
(c) The maximum output activations for network A after it has been trained (d)
The maximum output activations for a trained network B. (e) The training and
validation accuracy curves for network A and B.

We begin by constructing models of both hybrid and fully classical neural
networks, as sketched out in Fig. 1a and b, respectively. The difference
between the two is that the hidden layer of the classical network (shaded
region of Fig. 1b) is replaced by a CV quantum neural network (CVQNN) [15] as
shown by the shaded region of the hybrid network (Fig. 1a). Both cases have
identical input $\left(W^{(i)}\right)$ and output $\left(W^{(o)}\right)$
layers, typically fully-connected classical feed-forward networks, whose
dimensions depend on the task to be performed. Here, we select a
classification task because, while challenging, it is well-suited for feed-
forward neural networks. More specifically, we synthetically generate 1000
samples [22] (700 for training and 300 for validation) equally distributed
between 4 classes, each with 8 features (see Methods VII.1). An example of
this distribution, for 2 of the 8 features is shown in Fig. 1c, where the area
encompassed by each feature in each of the classes is shown by a shaded
region. The overlap of these regions demonstrates that it is impossible to
classify each sample simply based on 2 features, and explains why all 8
features must be considered. As a consequence, our input layer has 8 neurons,
while the output layer contains 4 neurons.

To create a hybrid neural network, the classical hidden layers
$\left(W^{(h)}\right)$ are replaced with an encoding layer [20, 21] followed
by a CV quantum neural network ($\mathcal{L}^{(h)}$, CVQNN) [15], as shown in
Fig. 1a. Such CVQNNs can be trained via backpropagation [23] (see
Supplementary Note LABEL:sec:parameter_shift) and they, along with the
encoding layer, can be realized with photonic elements on-chip [24, 25, 6,
26], meaning that both can be seamlessly integrated with classical photonic
neural networks. First, a series of displacement $\left(\mathcal{D}\right)$,
squeezing $\left(\mathcal{S}\right)$ and non-Gaussian Kerr $\left(\Phi\right)$
gates encodes the classical information into quantum modes (qumodes) (see
Supplementary Note LABEL:sec:gates for more information on the quantum gates.)
Note that the Kerr gates provide a non-linearity that, while difficult to
simulate classically [27], results in highly trainable and more performant
networks [28].

Once encoded into qumodes, the information flows through the CVQNN as a
quantum state, before exiting through the classical output layer. As shown in
Fig. 1a, the CVQNN is comprised of $\mathcal{D}$, $\mathcal{S}$ and $\Phi$
gates, parameterized by $s$, $\alpha$ and $\kappa$, respectively, interspersed
with linear interferometers $\mathcal{U}$ that are controlled by a set of
phase shifters, $\vec{\theta}_{i}$ and $\vec{\phi}_{i}$, that can be used to
perform any unitary operation on the qumodes. Altogether, $s$, $\alpha$,
$\kappa$, $\vec{\theta}$ and $\vec{\phi}$ are the parameters of the CVQNN
that, together, with the weights of the classical input and output layers, are
trained using conventional gradient descent.

In sum, each network is characterized by $P$ parameters, which depend on its
type and geometry. For a classical network, this is simply the number of
weights in weight matrix $W$ and bias vector $b$ summed over all layers. In
contrast, for a hybrid network,

$P=5M(I+1)+L\left(M(M-1)+7M\right)+O(M+1),$ (1)

where $I$ is the input dimension, $M$ is the number of qumodes, $L$ is the
number of layers, and $O$ is the output dimension. Here, the first and third
terms represent the number of parameters in the classical input and output
layers, while the middle term is the number of parameters in the quantum
layer. Fig. 1d shows the scaling of $P$ as a function of the number of qumodes
in the CVQNN (hidden classical layer), for 1-, 3- and 5-layer networks. For
this range, it is possible to find both classical and hybrid networks with $P$
ranging between about 100 and 600. We are therefore able to identify both
types of networks with similar number of parameters (and hence, expected
complexity) and, in what follows, we set out to compare their performance.

##  III Training and validation

Figure 3: Statistical violin plots of all hybrid and classical networks
trained. (a) Annotated violins of the 120 parameter hybrid network and 124
parameter classical network. Shaded blue and orange region shows the region of
well trained networks that scored better than 72% (accuracy of linear support
vector machine model).The point on each violin is the average of the well
trained networks with the error bars showing one standard deviation on this
mean. The red shaded region are considered poorly trained networks while
networks below the solid red failed to optimize at all. The white line is
networks A and B seen in Fig. 2 and Fig.1ab. (c) Violin plot of all networks
trained in this study. Each violin is a statistical fit (kernel density
estimation) of all networks with a given number of parameters where each grey
horizontal line is the maximum accuracy of one network during training. The
line is the mean of the upper of the well trained networks networks and the
shaded region is the standard deviation on that mean. (b) Bar plot showing the
percentage of networks which failed to beat the non-probabilistic linear model
($<72\%$ accuracy, poorly trained).

We begin by constructing models of the neural networks, implementing classical
layers in Tensorflow [29] and quantum layers using Pennylane and
Strawberryfields [30, 31] (see Methods VII.3 and VII.2). We assess each
network accuracy, before, during and after training, by comparing the class
that it assigns (i.e., output neuron $o_{i}$ with the highest activation) to
each of the 700 training (or 300 validation) samples. The overall network
accuracy is the number of correct classifications over the size of the
training (validation) set.

As an example, we consider the accuracy of a 120 and 124 parameter hybrid and
classical neural network, corresponding to the 1-layer, 2-qumode and 2-hidden
neuron networks shown in Fig. 1a and b, respectively. We first show how well
the hybrid network can classify the validation set at the start of the
training, plotting the predicted class for each sample (color coded by the
true class) along with the success probability assigned (i.e., value of the
highest output neuron). We observe an almost-random distribution of true class
vs. predicted class for the unoptimized hybrid network, with an overall
accuracy of only $0.58$. Similarly, the maximum output activation for this
unoptimized hybrid network is poor, peaking at $0.61$ for class 3 (see inset),
but generally well-below 0.4. That is, before training, the hybrid network is
both ‘hesitant’ and inaccurate. The results of the untrained classical
network, shown in Fig. 2b, are qualitatively similar but poorer, with an
overall accuracy of only $0.47$ and peak output probability of $0.38$.

We train the network parameters over 200 epochs, updating the weights 22 times
per epoch. We calculate the network accuracy after each epoch, using both the
training (solid curves) and validation (dashed curves) sets, showing the
results for both exemplary networks in Fig. 2e. Note that the overlay of the
curves indicates that there is no overfitting. As expected, the network
accuracy increase, reaching $0.87$ and $0.79$ for the hybrid and classical
systems, demonstrating the advantage of incorporating quantum layers.

The advantage of incorporating quantum layers can be seen more clearly when
considering the classification performance for both hybrid and classical
networks after training (at epoch 195), shown in Fig. 2c and d, respectively.
The hybrid network consistently classifies all samples correctly, with a
maximal probability output over $0.80$ across all 4 classes. In contrast,
while the fully classical network identifies class 1 and 4 samples, with
output probabilities reaching $0.91$, it struggles with class 2 and 3 samples.
For class 3 samples, in particular, the trained fully-classical network is
already outperformed by the untrained hybrid network.

##  IV Larger Networks

To more fully explore the difference between hybrid and classical photonic
neural networks, we repeat our analysis of the previous section with networks
of varying qumodes/neurons and layers (c.f., Fig. 1d). In total we train 342
hybrid networks and 518 classical networks, split in groups of 10 to 30
(hybrid) and 20 (classical) in sizes that range from 69 to 590 parameters. For
the classical devices, this corresponds to roughly single-layer $8\times 8$
and $24\times 24$ classical networks. To determine a threshold for well
trained networks a non-probabilistic linear support vector machine was fit to
the data achieving an accuracy of 0.72 (see Supplementary Note
LABEL:sec:linear_fitting). Neural networks above this threshold are considered
well-trained, while any network with an accuracy around $0.25$ is considered
failed (as a random selection should result in an accuracy of 0.25). The
aggregate results for the 120 (hybrid) and 124 (classical) parameter networks
are shown in Fig. 3a. This result further reinforces the power of the quantum
layers, as it makes clear that many more highly accurate hybrid networks are
trained than classical ones. More specifically, the average well-trained
hybrid network has an accuracy of $0.85\pm 0.01$ while the fully classical
networks have an average accuracy of only $0.79\pm 0.03$, where the larger
variance directly follows from the difficulty of training. In fact, for this
network size, only $6.9\%$ of hybrid networks had an average accuracy below
$72\%$ and none failed in training, in contrast with the $50\%$ of classical
networks with an accuracy below the threshold and the $4.2\%$ that failed.
That is, as expected, the network performance and trainability significantly
improves due to an increased informational capacity and improved optimization
landscape [28].

We repeat this analysis for all network sizes and present the results in Fig.
3c. Here, we observe that for all networks with 316 parameters or less, the
average accuracy of the well-trained hybrid networks exceeds that of the
classical networks, while the accuracy distribution for the hybrid networks
remains much smaller than that of the classical devices. Hybrid networks were
also much more likely to optimize beyond the $0.72$ accuracy threshold, with
95.9% of all trained hybrid networks reaching this threshold compared to the
76.3% of classical networks (the rate of poorly trained networks can be seen
in Fig.3b). This further indicates that hybrid networks are much more likely
to optimize successfully. More so, the average accuracy of $0.85\pm 0.03$ for
the 120 parameter hybrid network is only matched by classical networks at 235
parameters ($0.86\pm 0.04$). That is, the inclusion of quantum layers
consistently allows smaller hybrid networks to outperform similarly sized
classical networks and even those that are significantly larger.

Although Fig. 3c suggests that there is no benefit to including quantum layers
for networks with more than 316 parameters, this is likely not true. Rather,
this may be a direct consequence of the difficulty of simulating large quantum
networks which limits the cutoff dimension for our large CVQNNs to 5 (while
smaller networks were simulated up to a cutoff of 11). Since the amount of
information that can be encoded into each mode of the network is dependent on
this, so to will the accuracy of the network (see Supplementary Note
LABEL:sec:simulation_cutoff for more details), while also requiring longer
optimization times. We note, however, that in contrast to in simulations,
experimentally increasing the cutoff dimensions is straightforward.

##  V Robustness to Noise

Figure 4: Analysis of the impact of noise on hybrid and classical network
performance. (a) The validation accuracy of a 120 parameter hybrid network and
124 parameter classical network as a function of the ENOB of the parameters.
(b) The validation accuracy of a 120 parameter hybrid network with different
ENOB on the different gate parameters. Kerr and interferometer gates are the
most sensitive to bit precision.

Photonic neural networks are inherently analog systems of which precision is
ultimately limited by noise [32]. This noise comes from a combination of
thermal fluctuations, noise in the weight control currents, vibrational
changes in input coupling, or detector noise. Regardless of the source, noise
will ultimately determine how well a given weight in the neural network can be
actuated. For convenient comparison to digital bit precision, noise is often
converted to an effective number of bits (ENOB) which is determined according
to the Shannon-Hartley theorem as,

$ENOB=\log_{2}\left(1+\frac{w_{max}-w_{min}}{\sigma}\right).$ (2)

Here, $w_{max}$ and $w_{min}$ are the respective allowed maximum and minimum
values of the given weight $w$, and $\sigma$ is the standard deviation of the
noise on $w$. Thus, the ENOB represents the signal to noise ratio for an
analog signal.

For a classical layer of a photonic neural network, all parameters have the
range $[w_{min},w_{max}]=[-1,1]$ due to optical transmission in a balanced
photo detection scheme being fully in one detector, $+1$, or fully in the
other detector, $-1$ [33, 9, 34]. In contrast, a quantum layer has both
amplitude and phase parameters, where for phase $[w_{min},w_{max}]=[0,2\pi]$
and for amplitude parameters $[w_{min},w_{max}]=[0,a_{max}]$, where $a_{max}$
is the maximum amplitude parameter in the whole network. Thus, we can chose a
$\sigma$ for the noise and from that calculate the ENOB for the network and,
by rerunning our simulations calculate the network accuracy in the presence of
noise.

Fig. 4a shows exemplary curves for the network accuracy as a function of the
ENOB for both the 120 parameter hybrid network (blue) and the 124 parameter
classical network. For each ENOB value, the validation accuracy was calculated
ten times for different randomized noise values, with the mean shown by the
dark curve and the shaded region representing one standard deviation. The
classical network achieves near ideal accuracy, (10% worse than ideal
performance) at 5.5 bits of precision while the hybrid network requires 6.3
bits of precision. Both of these values are below the state-of-the-art 9 bits
of precision that are achievable with photonic neural networks [35].

Focusing on just the noise in the quantum neural network, we show how the
hybrid network accuracy depends on the ENOB of the different gate parameters
in Fig. 4b. To achieve near ideal accuracy, in this case 0.78, the
displacement, squeezing, kerr, and interferometer gates must have 5.0, 3.5,
6.0, and 5.0 ($\pm 0.5$) bits of precision, respectively. This demonstrates
that the network performance is most sensitive to the operation of
$\mathcal{U}$, which provides the same functionality as the weights of a
classical layer, and $\Phi$, which provides the non-linearity. The relatively
low required bit precision on the squeeze gate indicates that the precise
magnitude of the squeezing is not important, only that squeezing is present
(agreeing with current opinions in the field [36]). Regardless, we observe
that the overall network performance is robust in the presence of noise,
simply needing to exceed an experimentally feasible value.

##  VI Discussion and Conclusion

Here, we have shown that hybrid neural networks based on the continuous
variable quantum formalism could provide real performance and trainability
benefits over fully classical networks of the same size, providing a route to
increased complexity with fewer photonic resources. As an example, we showed
that for a complex classification task, a small hybrid network with 120
trainable parameters achieved an average well trained accuracy of $0.85\pm
0.03$, similar in performance to a much larger 235 parameter classical
network. Across a range of network sizes, the hybrid networks were always more
trainable and in most cases achieved greater performance. Encouragingly,
hybrid network achieved peak performance with only 6.3 effective bits of
precision, well within the capabilities of current photonic platforms [35, 8].
We expect that such networks could also be trained with noise present to
explore the performance benefits of noise aware training.

Hybrid networks, such as we envision, are implementable on a integrated
photonic platform and would also be directly compatible with existing photonic
communication platforms [6, 37]. While we focused on a classification task, it
is clear that they could be applied to realize other functionalities. For
example, in conjunction with direct integration to existing communication
platforms, hybrid networks have potential applications in quantum repeaters
[38] and photonic edge computing [39]. Scaling to larger and more complex
neural networks is a requirement to solve difficult computational tasks often
seen in these applications. Typically, several multiplexing approaches such as
wavelength, time, and spatial division multiplexing are used to provide denser
processing, thereby increasing the number of parameters, and adding more
computational power [40]. However, multiplexing approaches have drawbacks such
as slower speeds, larger footprints, and more complicated electrical control.
Here we have shown that hybrid networks provide an additional approach to
scale network complexity and solve more difficult computational tasks with
fewer parameters.

##  VII Methods

###  VII.1 Synthetic Dataset

Data was generated using the Scikit-Learn [22] `make_classification()`
function. As mentioned in the text, this dataset featured 8 inputs, 4 classes,
and 1000 samples. All features were informative, meaning they included non-
redundant information, and each class featured 3 clusters of datapoints. On
average, the clusters of datapoints were seperated by an eight dimensional
vector of length 3.0 (the default is 1). 2% of samples were assigned a random
class to increase difficulty. The dataset was normalized to the domain $[0,1]$
prior to training to simulate input transmission seen in photonic circuits
[34]. The same random state was used each time the fata was generated.

###  VII.2 Classical Model

The classical models (and layers) were implemented in the Tensorflow
framework. ReLU non-linearities were used between each of the classical
layers. All classical layers were implemented using Keras Dense layers with
weight clipping beyond the domain $[-1,1]$ [34]. This simulated the range of
transmission values available on photonic hardware. The size of the classical
hidden layers was determined using an algorithm which, for each classical
layer, matched the number of parameters in the corresponding layer of the
hybrid network. In this way, both the total number of layers (circuit depth)
and total number of parameters were equivalent. The final output of the
network used a softmax non-linearity to convert the output activations to a
probability distribution for classification.

###  VII.3 Hybrid Model

The hybrid models were fully implemented in the Tensorflow framework [29] with
the Pennylane quantum machine learning Tensorflow interface [30] and the
Strawberryfields Tensorflow based simulator [31]. This avoided array
conversions between different frameworks and enabled automated gradient
calculations, GPU support, and access to machine learning functions. A Fock
basis simulator was used as it can simulate the non-gaussian gate
($\hat{O}\sim e^{\hat{n}^{2}}$) used in the non-linear activation function
(See Supplementary Note LABEL:sec:simulation_cutoff). The Fock simulator,
fully implemented in Tensorflow, records the amplitudes for each of the photon
number states in the Fock basis as gates are applied.

Data begins at the input layer which, unlike the classical network, has no
activation function. After the input layer, is the encoding layer into the
quantum circuit. The encoding layer accepts 5 inputs for each qumode in the
circuit. An amplitude and phase for the squeeze gate, an amplitude and phase
for the displacement gate, and an amplitude and phase for the Kerr gate. The
inputs into the encoding layer are scaled to the domain $[0,2\pi]$, for phase
parameters, and to $[0,a_{max}]$, for amplitude parameters, using the
following functions:

$\tilde{y}_{i}^{(amplitude)}=a_{max}\mathrm{Sig}(y_{i})$
$\tilde{y}_{j}^{(phase)}=2\pi\mathrm{Sig}(y_{j})$

Where $\tilde{y}_{\\{i,j\\}}^{(amplitude,phase)}$ is the scaled value,
$\mathrm{Sig}(y)$ is the sigmoid function, $a_{max}$ is the maximum
displacement/squeezing amplitude, and $y_{\\{i,j\\}}$ is the output from the
previous layer. $a_{max}$ was determined using a brute force method outlined
in Supplementary Note LABEL:sec:trace_normalization. Acting on this encoded
state is the CV quantum neural network layer made up of N-port interferometers
($\mathcal{U}_{1},\mathcal{U}_{2}$), squeezing gates ($\mathcal{S}$),
displacement gates ($\mathcal{D})$, and Kerr non-linearities ($\Phi$). The
first sequence of gates, $\mathcal{U}_{2}\mathcal{S}\mathcal{U}_{1}$,
represents the matrix multiplication step, the displacement gate $\mathcal{D}$
represents the bias, and the kerr non-linearity, $\Phi$, acts as the
activation function [15]. The amplitude parameters were initialized uniformly
on the domain $[0,a_{max}]$ and phase parameters on the domain $[0,2\pi]$. All
amplitude values were L1 regularized to improve state normalization during
training. For mathematical and graphical descriptions of these gates refer to
Supplementary Fig. LABEL:fig:CV_gate_progression and Table LABEL:tab:gates.
Models with 2-4 qumodes and 1-5 layers were trained. Homodyne detection in the
position basis ($\hat{x}$) is used to measure the output value of each of the
qumodes. These measurements are then inputted into the final classical output
layer which is used to make a classification.

###  VII.4 Training

All networks were trained using conventional gradient descent and
backpropagation. The Adam optimizer with a learning rate of $\eta=0.001$ was
used with a batch size of 32. A categorical cross entropy loss function was
used to calculate the loss as this is a multi class classification task. 10-30
hybrid networks (three different cutoff dimension values, 10 for each cutoff)
and 20 classical were trained for each network size. Training was conducted on
the Digital Research Alliance of Canada’s Graham cluster.

##  VIII Acknowledgments

All authors acknowledge the support of the Natural Sciences and Engineering
Research Council of Canada (NSERC), Digital Research Alliance of Canada
(DRAC), and Queen’s University. TA and BJS acknowledge the support of the
Vector Institute. BJS is supported by the Canada Research Chairs program. NR
and BJS acknowledge the support of the Canadian Foundation for Innovation
(CFI).

##  IX Author Contributions

The project was originally conceived by TA and BJS. AH provided initial code
and, with SB, useful discussion. TA performed all simulations to generated
data and figured. TA prepared the original manuscript which was edited by BJS
and NR. Research led by BJS and NR.

## References

  * [1] Bhavin J. Shastri, Alexander N. Tait, T. Ferreira de Lima, Wolfram H. P. Pernice, Harish Bhaskaran, C. D. Wright, and Paul R. Prucnal.  Photonics for artificial intelligence and neuromorphic computing.  15(2):102–114.  ISSN 1749-4885, 1749-4893.  doi: 10.1038/s41566-020-00754-y.  URL http://www.nature.com/articles/s41566-020-00754-y.
  * Zhang et al. [a] Weipeng Zhang, Joshua C. Lederman, Thomas Ferreira de Lima, Jiawei Zhang, Simon Bilodeau, Leila Hudson, Alexander Tait, Bhavin J. Shastri, and Paul R. Prucnal.  A system-on-chip microwave photonic processor solves dynamic RF interference in real time with picosecond latency.  13(1):14, a.  ISSN 2047-7538.  doi: 10.1038/s41377-023-01362-5.  URL https://www.nature.com/articles/s41377-023-01362-5.  Publisher: Nature Publishing Group.
  * [3] Hailong Zhou, Jianji Dong, Junwei Cheng, Wenchan Dong, Chaoran Huang, Yichen Shen, Qiming Zhang, Min Gu, Chao Qian, Hongsheng Chen, Zhichao Ruan, and Xinliang Zhang.  Photonic matrix multiplication lights up photonic accelerator and beyond.  11(1):30.  ISSN 2047-7538.  doi: 10.1038/s41377-022-00717-8.  URL https://www.nature.com/articles/s41377-022-00717-8.  Publisher: Nature Publishing Group.
  * [4] Muhammad Kashif and Saif Al-Kuwari.  Design space exploration of hybrid quantum–classical neural networks.  10(23).  ISSN 2079-9292.  doi: 10.3390/electronics10232980.  URL https://www.mdpi.com/2079-9292/10/23/2980.
  * [5] Yichen Shen, Nicholas C. Harris, Scott Skirlo, Mihika Prabhu, Tom Baehr-Jones, Michael Hochberg, Xin Sun, Shijie Zhao, Hugo Larochelle, Dirk Englund, and Marin Soljačić.  Deep learning with coherent nanophotonic circuits.  11(7):441–446.  ISSN 1749-4893.  doi: 10.1038/nphoton.2017.93.  URL https://www.nature.com/articles/nphoton.2017.93.  Publisher: Nature Publishing Group.
  * [6] J. M. Arrazola, V. Bergholm, K. Brádler, T. R. Bromley, M. J. Collins, I. Dhand, A. Fumagalli, T. Gerrits, A. Goussev, L. G. Helt, J. Hundal, T. Isacsson, R. B. Israel, J. Izaac, S. Jahangiri, R. Janik, N. Killoran, S. P. Kumar, J. Lavoie, A. E. Lita, D. H. Mahler, M. Menotti, B. Morrison, S. W. Nam, L. Neuhaus, H. Y. Qi, N. Quesada, A. Repingon, K. K. Sabapathy, M. Schuld, D. Su, J. Swinarton, A. Száva, K. Tan, P. Tan, V. D. Vaidya, Z. Vernon, Z. Zabaneh, and Y. Zhang.  Quantum circuits with many photons on a programmable nanophotonic chip.  591(7848):54–60.  ISSN 1476-4687.  doi: 10.1038/s41586-021-03202-1.  URL https://www.nature.com/articles/s41586-021-03202-1.  Publisher: Nature Publishing Group.
  * [7] Kun Liao, Chentong Li, Tianxiang Dai, Chuyu Zhong, Hongtao Lin, Xiaoyong Hu, and Qihuang Gong.  Matrix eigenvalue solver based on reconfigurable photonic neural network.  11(17):4089–4099.  ISSN 2192-8614.  doi: 10.1515/nanoph-2022-0109.  URL https://www.degruyter.com/document/doi/10.1515/nanoph-2022-0109/html?lang=en.  Publisher: De Gruyter.
  * [8] Farshid Ashtiani, Alexander J. Geers, and Firooz Aflatouni.  An on-chip photonic deep neural network for image classification.  606(7914):501–506.  ISSN 1476-4687.  doi: 10.1038/s41586-022-04714-0.  URL https://www.nature.com/articles/s41586-022-04714-0.  Publisher: Nature Publishing Group.
  * Tait et al. [a] Alexander N. Tait, Thomas Ferreira de Lima, Ellen Zhou, Allie X. Wu, Mitchell A. Nahmias, Bhavin J. Shastri, and Paul R. Prucnal.  Neuromorphic photonic networks using silicon photonic weight banks.  7(1):7430, a.  ISSN 2045-2322.  doi: 10.1038/s41598-017-07754-z.  URL https://www.nature.com/articles/s41598-017-07754-z.  Number: 1 Publisher: Nature Publishing Group.
  * [10] Lorenzo De Marinis and Nicola Andriolli.  Photonic integrated neural network accelerators.  In _Photonics in Switching and Computing 2021 (2021), paper W3B.1_ , page W3B.1. Optica Publishing Group.  doi: 10.1364/PSC.2021.W3B.1.  URL https://opg.optica.org/abstract.cfm?uri=PSC-2021-W3B.1.
  * [11] Zhongjin Lin, Bhavin J. Shastri, Shangxuan Yu, Jingxiang Song, Yuntao Zhu, Arman Safarnejadian, Wangning Cai, Yanmei Lin, Wei Ke, Mustafa Hammood, Tianye Wang, Mengyue Xu, Zibo Zheng, Mohammed Al-Qadasi, Omid Esmaeeli, Mohamed Rahim, Grzegorz Pakulski, Jens Schmid, Pedro Barrios, Weihong Jiang, Hugh Morison, Matthew Mitchell, Xiaogang Qiang, Xun Guan, Nicolas A. F. Jaeger, Leslie A. n Rusch, Sudip Shekhar, Wei Shi, Siyuan Yu, Xinlun Cai, and Lukas Chrostowski.  65 GOPS/neuron photonic tensor core with thin-film lithium niobate photonics.  URL http://arxiv.org/abs/2311.16896.
  * [12] Paul R. Prucnal and Bhavin J. Shastri.  _Neuromorphic Photonics_.  CRC Press.  ISBN 978-1-4987-2524-8.  Google-Books-ID: VbvODgAAQBAJ.
  * [13] Alexander Norman Tait.  Silicon photonic neural networks.  URL https://dataspace.princeton.edu/handle/88435/dsp01vh53wz43k.  Accepted: 2018-04-26T18:46:40Z Publisher: Princeton, NJ : Princeton University.
  * [14] Peter L. McMahon.  The physics of optical computing.  5(12):717–734.  ISSN 2522-5820.  doi: 10.1038/s42254-023-00645-5.  URL https://www.nature.com/articles/s42254-023-00645-5.  Publisher: Nature Publishing Group.
  * Killoran et al. [a] Nathan Killoran, Thomas R. Bromley, Juan Miguel Arrazola, Maria Schuld, Nicolás Quesada, and Seth Lloyd.  Continuous-variable quantum neural networks.  1(3):033063, a.  doi: 10.1103/PhysRevResearch.1.033063.  URL https://link.aps.org/doi/10.1103/PhysRevResearch.1.033063.  Publisher: American Physical Society.
  * [16] Sophie Choe.  Continuous variable quantum MNIST classifiers.  URL http://arxiv.org/abs/2204.01194.
  * [17] Shikha Bangar, Leanto Sunny, Kübra Yeter-Aydeniz, and George Siopsis.  Experimentally realizable continuous-variable quantum neural networks.  108(4):042414.  doi: 10.1103/PhysRevA.108.042414.  URL https://link.aps.org/doi/10.1103/PhysRevA.108.042414.  Publisher: American Physical Society.
  * [18] Lars S. Madsen, Fabian Laudenbach, Mohsen Falamarzi Askarani, Fabien Rortais, Trevor Vincent, Jacob F. F. Bulmer, Filippo M. Miatto, Leonhard Neuhaus, Lukas G. Helt, Matthew J. Collins, Adriana E. Lita, Thomas Gerrits, Sae Woo Nam, Varun D. Vaidya, Matteo Menotti, Ish Dhand, Zachary Vernon, Nicolás Quesada, and Jonathan Lavoie.  Quantum computational advantage with a programmable photonic processor.  606(7912):75–81.  ISSN 1476-4687.  doi: 10.1038/s41586-022-04725-x.  URL https://www.nature.com/articles/s41586-022-04725-x.  Number: 7912 Publisher: Nature Publishing Group.
  * [19] Saumil Bandyopadhyay, Alexander Sludds, Stefan Krastanov, Ryan Hamerly, Nicholas Harris, Darius Bunandar, Matthew Streshinsky, Michael Hochberg, and Dirk Englund.  Single chip photonic deep neural network with accelerated training.  URL http://arxiv.org/abs/2208.01623.
  * [20] Maria Schuld and Nathan Killoran.  Quantum machine learning in feature hilbert spaces.  122(4):040504.  doi: 10.1103/PhysRevLett.122.040504.  URL https://link.aps.org/doi/10.1103/PhysRevLett.122.040504.  Publisher: American Physical Society.
  * [21] Vojtěch Havlíček, Antonio D. Córcoles, Kristan Temme, Aram W. Harrow, Abhinav Kandala, Jerry M. Chow, and Jay M. Gambetta.  Supervised learning with quantum-enhanced feature spaces.  567(7747):209–212.  ISSN 1476-4687.  doi: 10.1038/s41586-019-0980-2.  URL https://www.nature.com/articles/s41586-019-0980-2.  Publisher: Nature Publishing Group.
  * [22] Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, Jake Vanderplas, Alexandre Passos, David Cournapeau, Matthieu Brucher, Matthieu Perrot, and Édouard Duchesnay.  Scikit-learn: Machine learning in python.  12(85):2825–2830.  ISSN 1533-7928.  URL http://jmlr.org/papers/v12/pedregosa11a.html.
  * [23] Gavin E. Crooks.  Gradients of parameterized quantum gates using the parameter-shift rule and gate decomposition.  URL http://arxiv.org/abs/1905.13311.
  * [24] V. D. Vaidya, B. Morrison, L. G. Helt, R. Shahrokshahi, D. H. Mahler, M. J. Collins, K. Tan, J. Lavoie, A. Repingon, M. Menotti, N. Quesada, R. C. Pooser, A. E. Lita, T. Gerrits, S. W. Nam, and Z. Vernon.  Broadband quadrature-squeezed vacuum and nonclassical photon number correlations from a nanophotonic device.  6(39):eaba9186.  doi: 10.1126/sciadv.aba9186.  URL https://www.science.org/doi/10.1126/sciadv.aba9186.  Publisher: American Association for the Advancement of Science.
  * Zhang et al. [b] Y. Zhang, M. Menotti, K. Tan, V. D. Vaidya, D. H. Mahler, L. G. Helt, L. Zatti, M. Liscidini, B. Morrison, and Z. Vernon.  Squeezed light from a nanophotonic molecule.  12(1):2233, b.  ISSN 2041-1723.  doi: 10.1038/s41467-021-22540-2.  URL https://www.nature.com/articles/s41467-021-22540-2.  Number: 1 Publisher: Nature Publishing Group.
  * [26] A. I. Lvovsky and M. G. Raymer.  Continuous-variable optical quantum-state tomography.  81(1):299–332.  ISSN 0034-6861, 1539-0756.  doi: 10.1103/RevModPhys.81.299.  URL https://link.aps.org/doi/10.1103/RevModPhys.81.299.
  * [27] Pieter Kok and Brendon W. Lovett.  _Introduction to Optical Quantum Information Processing_.  Cambridge University Press, 1 edition.  ISBN 978-0-521-51914-4 978-1-139-19365-8.  doi: 10.1017/CBO9781139193658.  URL https://www.cambridge.org/core/product/identifier/9781139193658/type/book.
  * [28] Amira Abbas, David Sutter, Christa Zoufal, Aurelien Lucchi, Alessio Figalli, and Stefan Woerner.  The power of quantum neural networks.  1(6):403–409.  ISSN 2662-8457.  doi: 10.1038/s43588-021-00084-1.  URL https://www.nature.com/articles/s43588-021-00084-1.  Publisher: Nature Publishing Group.
  * [29] Martín Abadi, Paul Barham, Jianmin Chen, Zhifeng Chen, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Geoffrey Irving, Michael Isard, Manjunath Kudlur, Josh Levenberg, Rajat Monga, Sherry Moore, Derek G. Murray, Benoit Steiner, Paul Tucker, Vijay Vasudevan, Pete Warden, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng.  TensorFlow: A system for large-scale machine learning.  URL http://arxiv.org/abs/1605.08695.
  * [30] Ville Bergholm, Josh Izaac, Maria Schuld, Christian Gogolin, Shahnawaz Ahmed, Vishnu Ajith, M. Sohaib Alam, Guillermo Alonso-Linaje, B. AkashNarayanan, Ali Asadi, Juan Miguel Arrazola, Utkarsh Azad, Sam Banning, Carsten Blank, Thomas R. Bromley, Benjamin A. Cordier, Jack Ceroni, Alain Delgado, Olivia Di Matteo, Amintor Dusko, Tanya Garg, Diego Guala, Anthony Hayes, Ryan Hill, Aroosa Ijaz, Theodor Isacsson, David Ittah, Soran Jahangiri, Prateek Jain, Edward Jiang, Ankit Khandelwal, Korbinian Kottmann, Robert A. Lang, Christina Lee, Thomas Loke, Angus Lowe, Keri McKiernan, Johannes Jakob Meyer, J. A. Montañez-Barrera, Romain Moyard, Zeyue Niu, Lee James O’Riordan, Steven Oud, Ashish Panigrahi, Chae-Yeun Park, Daniel Polatajko, Nicolás Quesada, Chase Roberts, Nahum Sá, Isidor Schoch, Borun Shi, Shuli Shu, Sukin Sim, Arshpreet Singh, Ingrid Strandberg, Jay Soni, Antal Száva, Slimane Thabet, Rodrigo A. Vargas-Hernández, Trevor Vincent, Nicola Vitucci, Maurice Weber, David Wierichs, Roeland Wiersema, Moritz Willmann, Vincent Wong, Shaoming Zhang, and Nathan Killoran.  PennyLane: Automatic differentiation of hybrid quantum-classical computations.  URL http://arxiv.org/abs/1811.04968.
  * Killoran et al. [b] Nathan Killoran, Josh Izaac, Nicolás Quesada, Ville Bergholm, Matthew Amy, and Christian Weedbrook.  Strawberry fields: A software platform for photonic quantum computing.  3:129, b.  ISSN 2521-327X.  doi: 10.22331/q-2019-03-11-129.  URL http://arxiv.org/abs/1804.03159.
  * [32] Thomas Ferreira de Lima, Alexander N. Tait, Hooman Saeidi, Mitchell A. Nahmias, Hsuan-Tung Peng, Siamak Abbaslou, Bhavin J. Shastri, and Paul R. Prucnal.  Noise analysis of photonic modulator neurons.  26(1):1–9.  ISSN 1558-4542.  doi: 10.1109/JSTQE.2019.2931252.  URL https://ieeexplore.ieee.org/document/8782580.
  * Tait et al. [b] Alexander N. Tait, Mitchell A. Nahmias, Bhavin J. Shastri, and Paul R. Prucnal.  Broadcast and weight: An integrated network for scalable photonic spike processing.  32(21):4029–4041, b.  ISSN 1558-2213.  doi: 10.1109/JLT.2014.2345652.  URL https://ieeexplore.ieee.org/document/6872524.  Conference Name: Journal of Lightwave Technology.
  * [34] Viraj Bangari, Bicky Marquez, Heidi Miller, A.N. Tait, Mitchell Nahmias, Thomas Ferreira de Lima, Hsuan-Tung Peng, Paul Prucnal, and Bhavin Shastri.  Digital electronics and analog photonics for convolutional neural networks (DEAP-CNNs).  PP:1–1.  doi: 10.1109/JSTQE.2019.2945540.
  * Zhang et al. [c] Weipeng Zhang, Chaoran Huang, Hsuan-Tung Peng, Simon Bilodeau, Aashu Jha, Eric Blow, Thomas Ferreira de Lima, Bhavin J. Shastri, and Paul Prucnal.  Silicon microring synapses enable photonic deep learning beyond 9-bit precision.  9(5):579–584, c.  ISSN 2334-2536.  doi: 10.1364/OPTICA.446100.  URL https://opg.optica.org/optica/abstract.cfm?uri=optica-9-5-579.
  * [36] Olivier Pfister.  Continuous-variable quantum computing in the quantum optical frequency comb.  53(1):012001.  ISSN 0953-4075.  doi: 10.1088/1361-6455/ab526f.  URL https://dx.doi.org/10.1088/1361-6455/ab526f.  Publisher: IOP Publishing.
  * [37] Vicente Martin, Juan Pedro Brito, Carmen Escribano, Marco Menchetti, Catherine White, Andrew Lord, Felix Wissel, Matthias Gunkel, Paulette Gavignet, Naveena Genay, Olivier Le Moult, Carlos Abellán, Antonio Manzalini, Antonio Pastor-Perales, Victor López, and Diego López.  Quantum technologies in the telecommunications industry.  8(1):19.  ISSN 2662-4400, 2196-0763.  doi: 10.1140/epjqt/s40507-021-00108-9.  URL https://epjquantumtechnology.springeropen.com/articles/10.1140/epjqt/s40507-021-00108-9.
  * [38] Koji Azuma, Sophia E. Economou, David Elkouss, Paul Hilaire, Liang Jiang, Hoi-Kwong Lo, and Ilan Tzitrin.  Quantum repeaters: From quantum networks to the quantum internet.  95(4):045006.  doi: 10.1103/RevModPhys.95.045006.  URL https://link.aps.org/doi/10.1103/RevModPhys.95.045006.  Publisher: American Physical Society.
  * [39] Alexander Sludds, Saumil Bandyopadhyay, Zaijun Chen, Zhizhen Zhong, Jared Cochrane, Liane Bernstein, Darius Bunandar, P. Ben Dixon, Scott A. Hamilton, Matthew Streshinsky, Ari Novack, Tom Baehr-Jones, Michael Hochberg, Manya Ghobadi, Ryan Hamerly, and Dirk Englund.  Delocalized photonic deep learning on the internet’s edge.  378(6617):270–276.  doi: 10.1126/science.abq8271.  URL https://www.science.org/doi/10.1126/science.abq8271.  Publisher: American Association for the Advancement of Science.
  * [40] Yunping Bai, Xingyuan Xu, Mengxi Tan, Yang Sun, Yang Li, Jiayang Wu, Roberto Morandotti, Arnan Mitchell, Kun Xu, and David J. Moss.  Photonic multiplexing techniques for neuromorphic computing.  12(5):795–817.  ISSN 2192-8614.  doi: 10.1515/nanoph-2022-0485.  URL https://www.degruyter.com/document/doi/10.1515/nanoph-2022-0485/html?lang=en.  Publisher: De Gruyter.