# Non-linear Functional Modeling using Neural Networks

Aniruddha Rajendra Rao
Department of Statistics
Pennsylvania State University, USA
Matthew Reimherr*
Department of Statistics
Pennsylvania State University, USA

###### Abstract

We introduce a new class of non-linear models for functional data based on
neural networks. Deep learning has been very successful in non-linear
modeling, but there has been little work done in the functional data setting.
We propose two variations of our framework: a functional neural network with
continuous hidden layers, called the Functional Direct Neural Network (FDNN),
and a second version that utilizes basis expansions and continuous hidden
layers, called the Functional Basis Neural Network (FBNN). Both are designed
explicitly to exploit the structure inherent in functional data. To fit these
models we derive a functional gradient based optimization algorithm. The
effectiveness of the proposed methods in handling complex functional models is
demonstrated by comprehensive simulation studies and real data examples.

Keywords: Functional Data Analysis, Deep Learning, Functional Regression,
Functional Neural Network.

##  1 Introduction

In today’s world, repeated measurements taken over time or space are collected
across a variety of different areas of industry, government, and academia. For
instance, health monitoring of patients, daily weather information, and stock
prices, to name only a few. In such cases, the data of interest changes
continuously over time. One approach to modeling such data is functional data
analysis (Ramsay and Silverman,, 1997; Kokoszka and Reimherr,, 2018).

Functional Data Analysis (FDA) is a branch of statistics that enables
researchers to analyze and model relationships between variables measured over
a particular domain, such as time or space (Ramsay and Silverman,, 1997;
Kokoszka and Reimherr,, 2018; Ferraty and Romain,, 2011; Ferraty and Vieu,,
2006; Horváth and Kokoszka,, 2012). FDA has been a rich area of research in
the last few decades. The advantage FDA has over multivariate data analysis is
that FDA enables smoothing and interpolation procedures that can yield
functional representations of a finite set of observations. Also, FDA has
applicability to data that are recorded irregularly over time or some other
domain.

Functional regression can mainly be classified into three major categories
according to the role played by the functional data in each model: scalar
response and functional predictors (scalar-on-function regression); functional
response and scalar predictors (function-on-scalar regression); and functional
response and functional predictors (function-on-function regression). This
paper concentrates on the first case, that is the scalar-on-function
regression. Even though Functional Regression, especially scalar-on-function
regression, is perhaps one of the most thoroughly researched topics within the
broader literature on FDA, there is still a great need for superior non-linear
models. It is not only a very active area for research, but also is very
challenging with functional data due to the curse of dimensionality (Geenens,,
2011).

Some common ways of dealing with non-linearity in functional data are by using
the generalized additive model, also known as Continuously Additive Model
(CAM) (McLean et al.,, 2014; Fan et al.,, 2015; Ma and Zhu,, 2016; Maller et
al.,, 2013; Wang and Ruppert,, 2013; Reimherr et al.,, 2017), Quadratic Model
(Yao and Muller,, 2010), Single Index Model (Stoker,, 1986; Eilers et al.,,
2009; Ferraty,, 2011), and Multiple Index Model (James and Silverman,, 2004;
Chen et al.,, 2012; Ferraty et al.,, 2013). The single index model has the
outcome related to a linear functional of the predictors through a non-linear
transformation, whereas the multiple index model is analogous but with
multiple linear projections. In these models, you can either choose to capture
the complex non-linear relations or the higher order interactions. In
particular, the multiple index model allows for complex interactions but non-
linearity is weak because we are only using a few projections whereas in CAM,
we can capture the non-linearity but interactions are completely ignored, and
the quadratic model assumes a linear structure with interactions without any
provision for non-linearity. In general, these models assume a very particular
structure, which limits the kind of non-linear patterns that can be learned,
especially in terms of interactions within the time points for the predictors.
Most of these approaches have been discussed in the review papers by Reiss et
al., (2016) and Morris, (2014).

Neural Network (NN) and Deep Learning (DL) have attracted a great deal of
attention and research due to their application in areas like time series
(Dorffner,, 1996; Qin et al.,, 2017; Khashei and Bijari,, 2010), image
processing/computer vision (Gregor et al.,, 2015; Shi et al.,, 2016;
Kanellopoulos and Wilkinson,, 1997; Wang and Yang,, 2019; Zhang et al.,,
2017), speech processing/natural language processing (NLP) (Seide et al.,,
2011; Graves et al.,, 2013; Chan et al.,, 2016; Morchid,, 2018), genetics
(Singh and Singh,, 2017; Asadi et al.,, 2014), self-driving cars (Maqueda et
al.,, 2018; Chen and Huang,, 2017), and video learning (Yue-Hei Ng et al.,,
2015; Karpathy et al.,, 2014). Neural network models by combining multiple
linear combinations of the predictors after passing them through
activation/link functions (called hidden layers) to learn complex non-linear
relationships. Even with so much success, neural network and deep learning are
still relatively undeveloped for functional data. A traditional neural network
includes no special provision for the temporal dependency seen in FDA.
Recurrent Neural Network (RNN), especially Long Short-Term Memory (LSTM), have
an edge over conventional feed-forward neural networks when it comes to time
series or sequential data. This is due to their ability to selectively
remember patterns for a long duration of time compared to regular RNN (Joshi
et al.,, 2020). The main drawback of deep learning techniques is that they
can’t handle curves directly, especially if the curves are measured at
different time points or with a different number of observations, that is if
the data is irregular.

Rossi et al., (2002); Rossi and Conan-Guez, (2006); Conan-Guez and Rossi,
(2002) developed and demonstrated the efficacy of Functional Neural Network
(FNN). FNN was then further investigated by Wang et al., (2019, 2019), where
they demonstrated the superiority of FNN to state-of-the-art alternatives
(CNN, RNN, LSTM). The extension of the neural network to functional data is
intriguing, as it enables efficient non-linear learning when the predictors
are curves or surfaces. Rossi et al., (2002) proposed a way to feed functional
data as an input into a neural network and proved the validity by theoretical
arguments and simulation experiments. But their initial work (Rossi et al.,,
2002) required the data to be dense and regular. This requirement was overcome
with a slightly different FNN version (Conan-Guez and Rossi,, 2002) where they
used curve projections. Wang et al., (2019) extended this work to deal with
sparsity. The basic approach is that it only requires the entire time interval
to be well covered by pooled observations across all the subjects. They pool
measurements among all the subjects to estimate a set of basis functions and
then use these basis functions to build a sparse functional neuron that
extracts features for each subject. This method is developed for sparse FDA
but it is equivalent to the original FNN with dense data. The main idea of FNN
is that only the first layer is a functional hidden layer with functional
neurons that give a scalar output, which is then connected to regular hidden
layers and finally the output layer. In short, the functional layer is
followed by a traditional neural network. This is a straight forward approach
but it does not embrace the functional data completely. Other extensions of
FNN for Spatio-Temporal modeling and Function-on-Function regression were
explored in Rao et al., (2020); Wang et al., (2020). We use the idea of FNN as
an inspiration to develop our novel approach where we extend the initial
functional hidden layer with functional neurons to multiple continuous hidden
layers with continuous neurons with no regular hidden layers.

In this work, we study the problem of non-linear modeling from functional data
using neural networks. There are mainly two challenges: 1) How to design a
continuous hidden layer consisting of continuous neurons for a neural network
that can learn functional representation from the functional inputs; and 2)
how to train such a model. We propose a novel approach that takes the
functional input and gives a scalar output, which can be solved using either
Functional Direct Neural Network (FDNN) or Functional Basis Neural Network
(FBNN). The main contributions of this work are as follows:

  * •

We introduce a new structure for non-linear modeling of functional data with a
scalar response.

  * •

We give two strategies, FDNN and FBNN, to solve the problem.

  * •

We derive the functional gradient based learning algorithm for both strategies
to optimize the parameters.

  * •

We conduct simulation studies to demonstrate the effectiveness of FDNN and
FBNN in varying settings. We also apply the models to several real-world
examples.

The rest of this paper is organized as follows: In Section 2, we first present
some background on FDA and FNN. Then we describe our approach and give two
strategies to solve it, Functional Direct Neural Network (FDNN) and Functional
Basis Neural Network (FBNN). We follow this with simulations and several
benchmark applications in Section 3. We consider multiple simulation settings
depending on how the non-linear relation is defined between the response and
functional predictor. We demonstrate the efficiency of our methods over
current state-of-the-art approaches. In the final section, we present our
concluding remarks and future research directions which pertain to deeper
statistical theory and extension to a variety of functional data settings.

##  2 Methods

We begin by providing the necessary notation used in the paper, followed by a
summary of a few of the major non-linear models for scalar-on-function
regression before describing our methodology. Suppose we have $N\in\mathbb{N}$
subjects observed over a compact time interval $\mathcal{T}$, and for the
$i^{th}$ subject, the inputs are $R$ random functions that can be denoted as
{$X_{i,1},\ldots,X_{i,R}\\}\in\mathcal{L}^{2}(\mathcal{T})$. The output is a
scalar variable $Y_{i}\in\mathbb{R}$ for $i=1,\ldots,N$. We assume
$X_{i,r}:=\\{X_{i,r}(t):t\in[0,1]\\}$ for $r=1,\ldots,R$ are fully observed
functions, which is also known as dense functional data analysis (Kokoszka and
Reimherr,, 2018). The main objective is to learn the non-linear mapping,
$F:\mathcal{L}^{2}(\mathcal{T})\times\dots\times\mathcal{L}^{2}(\mathcal{T})\to\mathbb{R}$
from the functional inputs $X_{i,r}$ to the scalar output $Y_{i}$

$\text{E}[Y_{i}|X_{i,1},\ldots,X_{i,R}]=F\left(X_{i,1},\ldots,X_{i,R}\right).$
(1)

###  2.1 Functional Regression

There is vast literature in FDA when it comes to scalar-on-function regression
and we briefly describe the most popular ones, which will be later used in
simulations for comparison of different methods. In these models, the
integrals without limits are taken over the entire domain.

The scalar-on-function regression model under different relation between the
functional predictors and the scalar response is as follows:

  * •

Linear Model:

$\text{E}[Y_{i}|X_{i,1},\ldots,X_{i,R}]=\alpha+\sum_{r=1}^{R}\int\beta_{r}(t)X_{i,r}(t)dt$
(2)

where $\beta_{r}(t)$ is the parameter function (Kokoszka and Reimherr,, 2018;
Ramsay and Silverman,, 1997). This is also known as Funcitonal Linear Model
(FLM).

  * •

Continuously Additive Model (CAM):

$\text{E}[Y_{i}|X_{i,1},\ldots,X_{i,R}]=\sum_{r=1}^{R}\int f(X_{i,r}(t),t)dt$
(3)

where bivariate parameter function $(x,t)\rightarrow f(x,t)$ is smooth (McLean
et al.,, 2014; Maller et al.,, 2013; Wang and Ruppert,, 2013; Fan et al.,,
2015; Ma and Zhu,, 2016).

  * •

Single Index Model:

$\text{E}[Y_{i}|X_{i,1},\ldots,X_{i,R}]=g\left(\sum_{r=1}^{R}<\beta_{r},X_{i,r}>\right)$
(4)

where the function $g()$ is any smooth function defined on the real line
(Eilers et al.,, 2009; Chen et al.,, 2012).

  * •

Multiple Index Model:

$\text{E}[Y_{i}|X_{i,1},\ldots,X_{i,R}]=g\left(\sum_{r=1}^{R}<\beta_{1r},X_{i,r}>,...,\sum_{r=1}^{R}<\beta_{Pr},X_{i,r}>\right)$
(5)

Multiple Index models are extensions of the single index models (Chen et al.,,
2012; Ferraty et al.,, 2013; Muller and Yao,, 2008; Zhu et al.,, 2013).

  * •

Quadratic Model:

$\displaystyle\text{E}[Y_{i}|X_{i,1},\ldots,X_{i,R}]$
$\displaystyle=\sum_{r=1}^{R}\int\beta_{r}(t)X_{i,r}(t)dt+\sum_{r=1}^{R}\int\int\beta(s,t)X_{i,r}(t)X_{i,r}(s)dtds,$
(6) $\displaystyle\text{E}[Y_{i}|X_{i,1},\ldots,X_{i,R}]$
$\displaystyle=\sum_{r=1}^{R}\int f(X_{i,r}(t),t)dt+\sum_{r=1}^{R}\int\int
f(X_{i,r}(t),X_{i,r}(s),s,t)dtds$ (7)

Equation (6) represents a Quadratic Model given by Yao and Muller, (2010)
whereas Equation (7) represents a Complex Quadratic Model defined for
simulation purposes.

The above models are the most common approaches used in functional regression
for linear and non-linear cases. They are commonly fit using basis expansions,
Functional Principal Component Analysis (FPCA), or penalization methods. The
non-linear models are analogous to their multivariate counterparts, which help
to overcome the curse of dimensionality. In practice, it is difficult to know
the actual form of the relation between the response and the functional
predictors. Therefore, having a method that can deal with any kind of
relationship without knowing the truth will be very useful.

###  2.2 Functional Neural Network

A neural network consists of three components: an input layer, hidden layers,
and an output layer. The input layer and output layer consist of neurons/nodes
that are linked by multiple hidden layers consisting of neurons. In the
forward propagation, that is to construct a predicted value, each neuron sums
the weighted signal it receives from the input connections and produces an
output signal which is a non-linear function of these inputs. Back propagation
(Rumelhart et al.,, 1986) (that is gradient descent) is used to determine the
correct weights that minimize the overall cost in the neural network.

Figure 1 shows the architecture of the FNN model developed by Rossi et al.,
(2002) using the idea of a generalized neuron (Stinchcombe,, 1999). In FNN,
the first hidden layer is a functional hidden layer which consists of
functional neurons that give scalar outputs. These scalar values are then fed
into a regular neural network framework. A functional neuron is defined as
follows:

$H(X,W_{r})=\sigma\Big{(}b+\sum_{r=1}^{R}\int W_{r}(t)X_{r}(t)dt\Big{)},$ (8)

where $b\in\mathbb{R}$ is the unknown intercept, $W_{r}\left(t\right)$ is the
parameter function ($r=1,\ldots,R$) and $\sigma(\cdot)$ is a non-linear
activation function from $\mathbb{R}$ to $\mathbb{R}$.

Figure 1: The architecture of a FNN with two functional neurons in the first
layer and two numerical neurons in the second layer.

In forward propagation, as seen in Figure 1, we first pass through the
functional layer having functional neurons. For the value of $H(X,W_{r}(t))$,
we use Equation (8), where numerical approximation is used to get the integral
value. The forward propagation calculation in the subsequent hidden layers is
straightforward because of the regular neural network layers. In the
backpropagation step, we can go from the output layer to the second hidden
layer using the partial derivatives that can be easily computed. It is
critical to make sure that the partial derivatives at the first hidden layer
for the parameters $W_{r}(t)$ exist. Under certain assumptions (Rossi et al.,,
2002; Wang et al.,, 2019), we can compute the partial derivatives required in
the first hidden layer as:

$\frac{\partial H(X,W_{r})}{\partial
W_{r}}(t)=\sigma^{\prime}\Big{(}b+\sum_{r^{\prime}=1}^{R}\int
W_{r^{\prime}}(t^{\prime})X_{r^{\prime}}(t^{\prime})dt^{\prime}\Big{)}\times
X_{r}(t)$ (9)

where $\sigma^{\prime}(\cdot)$ represents the first derivative of the
activation function $\sigma(\cdot)$.

###  2.3 Functional Direct Neural Network and Functional Basis Neural Network

In our proposed approach, we develop a continuous mapping from layer to layer.
To do this, we make use of what we call continuous neurons and refer to a
hidden layer with multiple continuous neurons as a continuous hidden layer. In
our neural network framework, we have three types of layers: an input layer,
which is the first layer that takes in the functional predictors, the
continuous hidden layers containing continuous neurons, and a functional
output layer that provides the final fitted value using functional neurons.
Unlike neural network, which learns using scalar weights, we learn using
weight functions that are continuous over time (or some other continuum). The
$l^{th}$ continuous hidden layer and its $k^{th}$ continuous neuron is defined
as:

$\begin{split}H^{(l)}_{(k)}(s)&=\sigma\Big{(}b^{(l)}_{(k)}(s)+\sum_{j=1}^{J}\int
w^{(l)}_{(j,k)}(s,t)H^{(l-1)}_{(j)}(t)dt\Big{)}\\\
\widehat{Y}=H^{(L+1)}&=\sigma\Big{(}b^{(L+1)}+\sum_{j=1}^{J}\int
w^{(L+1)}_{(j)}(t)H^{(L)}_{(j)}(t)dt\Big{)}\end{split}$ (10)

where $l=1,2,3,...,L$, $H^{(0)}_{(j)}(t):=X_{j}(t)$, $b^{(L+1)}\in\mathbb{R}$
is the unknown intercept (also called the bias in machine learning),
$b^{(l)}_{(k)}(s)\in\mathcal{L}^{2}(T)$ is the unknown intercept function,
$w^{(l)}_{(j,k)}(s,t)$ is the bivariate parameter function for the $k^{th}$
continuous neuron in the $l^{th}$ hidden layer coming from the $j^{th}$
continuous neuron of the $(l-1)^{th}$ hidden layer, $w^{(L+1)}_{(j)}(t)$ is a
parameter function for output layer coming from the $j^{th}$ continuous neuron
of the $(L)^{th}$ hidden layer and $\sigma(\cdot)$ is a non-linear activation
function.

Comparing Figure 1 to Figure 2, and Equation (8) to Equation (10), we can
clearly see the difference between our approach and the functional neural
network defined by Rossi et al., (2002). Our approach uses continuous hidden
layers consisting of continuous neurons that output a function and learn
bivariate functional weights, unlike the FNN where only the first layer uses
functional neurons that output a scalar value. We use the functional neurons
only in the last layer to get the final output value. Also, unlike FNN where
the functional layer is connected to multiple regular hidden layers of
classical NN, our approach has multiple continuous hidden layers connected to
each other and a functional layer at the end, which gives the scalar output
value. A major motivation for this approach is to preserve the functional
nature of the data for as long as possible to provide a richer structure for
the prediction that also exploits the domain information and continuity of the
predictor functions.

We implement our proposed approach using two strategies: (1) a more direct
solution when the data is observed on a common grid or (2) using basis
expansions to allow for more irregular sampling of the functional predictors.
The former approach we call as Functional Direct Neural Network (FDNN) while
the latter approach we refer to as Functional Basis Neural Network (FBNN).
FDNN works well with densely observed data on a common grid. We can apply
standard tools from the neural network in FDNN directly, making it
straightforward to implement. We are learning the weight functions similar to
the way we learn the weights in the neural network and also have the
flexibility to change the grid length (s) between the continuous hidden
layers. FBNN can handle irregular data with the help of basis expansions by
pre-imputation and also allows for greater parsimony if the data is sampled on
a very dense grid, as we need to only learn the weights of the basis
expansion.

Figure 2: General architecture for our proposed approach with multiple
continuous hidden layers.

Functional Direct Neural Network (FDNN): For simplicity, we assume that Figure
2 has $R$ inputs (functional predictors), $L$ continuous hidden layers, each
with $J$ incoming and $K$ outgoing connections. FDNN has the same flexibility
as a neural network, where we can specify the number of continuous hidden
layers and the number of continuous neurons in each of them. Commonly used
activation functions in the continuous neurons include ReLU, tanh,
sigmoid/logistic and linear. The activation function used in the functional
output layer depends on the response. In our simulations, for the continuous
response, we use the linear activation function in the functional output
layer, whereas for the binary response, we use the sigmoid activation
function. In the forward propagation phase, we reach from the functional
predictors to the scalar response via a chain of multiple CAMs, and in the
backpropagation phase, we want to learn the weight functions and the intercept
functions.

We derive a functional gradient based optimization algorithm to optimize the
network parameters. The gradient of the weight function measures the change in
a functional to a change in a function on which the functional depends. The
needed assumptions and mathematical tools can be found in Rossi et al.,
(2002); Wang et al., (2019). We extend the gradient descent optimization
technique to the functional setting using calculus of variations (Olver,,
2019), in particular by taking Fréchet derivatives. In the backpropagation
step, we pass through the network in reverse, as we do in a neural network,
but we calculate the partial derivative for the weight functions, the
intercept functions and the intercept. Let the loss function be
$\mathcal{L}(\theta)$, where $\theta$ is the collection of all functional and
scalar parameters. The loss function is defined based on the response,
quadratic loss for a continuous response and a cross entropy loss for a binary
response. The partial derivatives required for the parameters of the
continuous neurons in the continuous hidden layers of our proposed
architecture (as shown in Figure 2) are as follows:

$\begin{split}\frac{\partial H^{(l)}_{(k)}}{\partial
b^{(l)}_{(k)}}(s)&=\sigma^{\prime}\Big{(}b^{(l)}_{(k)}(s)+\sum_{j=1}^{J}\int
w^{(l)}_{(j,k)}(s,t)H^{(l-1)}_{(j)}(t)dt\Big{)}\\\ \frac{\partial
H^{(l)}_{(k)}}{\partial
w^{(l)}_{(j,k)}}(s,t)&=\sigma^{\prime}\Big{(}b^{(l)}_{(k)}(s)+\sum_{j^{\prime}=1}^{J}\int
w^{(l)}_{(j^{\prime},k)}(s,t^{\prime})H^{(l-1)}_{(j^{\prime})}(t^{\prime})dt^{\prime}\Big{)}\times
H^{(l-1)}_{(j)}(t)\end{split}$ (11)

where $l=1,2,3,...,L$, $H^{(0)}_{(j)}(t):=X_{j}(t)$, and
$\sigma^{\prime}(\cdot)$ represents the first derivative of $\sigma(\cdot)$.

Next, the partial derivatives required for the parameters of the functional
neurons in the functional output layer are,

$\begin{split}\frac{\partial H^{(L+1)}}{\partial
b^{(L+1)}}&=\sigma^{\prime}\Big{(}b^{(L+1)}+\sum_{j=1}^{J}\int
w^{(L+1)}_{(j)}(t)H^{(L)}_{(j)}(t)dt\Big{)}\\\ \frac{\partial
H^{(L+1)}}{\partial
w^{(L+1)}_{(j)}}(t)&=\sigma^{\prime}\Big{(}b^{(L+1)}+\sum_{j^{\prime}=1}^{J}\int
w^{(L+1)}_{(j^{\prime})}(t^{\prime})H^{(L)}_{(j^{\prime})}(t^{\prime})dt^{\prime}\Big{)}\times
H^{(L)}_{(j)}(t)\end{split}$ (12)

Appendix A contains in-depth information regarding the partial derivatives. We
can consider the values of the grid/domain ($s$), the number of continuous
hidden layers, and the number of continuous neurons in each continuous hidden
layer as hyperparameters. The general way of optimizing for the parameters is
still the same in FDNN, using backpropagation with the help of the gradients
given above.

Functional Basis Neural Network (FBNN): Similar to FDNN, we have continuous
hidden layers with continuous neurons in FBNN, but the continuous neurons are
expanded using basis expansion. Unlike FDNN, we learn the weights for the
basis functions rather than the weight function itself. As we are using basis
functions, we can deal with irregular data more easily (Kokoszka and
Reimherr,, 2018). We get the forward feed in FBNN by using basis expansion
where we replace the weight functions ($w^{(l)}_{(j,k)}(s,t)$) in Equation
(10) by a linear combination of basis functions as seen below. We will again
see the details of our approach using the architecture shown in Figure 2. Let
subscripts $j,k$ represent the continuous neurons of the neural network and
superscript $l$ denote the $l^{th}$ layer for the weight parameter
$(W^{(l)}_{(j,k)})$. We define $v$ to be some standard collection of basis
functions, like splines, wavelets, or sine and cosine functions. Also, let us
assume that we use the same basis function and the same number of basis in
each of the continuous neurons in the continuous hidden layer across the
network for simplicity. This gives us the following Equation:

$\begin{split}H^{(l)}_{(k)}(s)&=\sigma\Big{(}b^{(l)}_{(k)}(s)+\sum_{j=1}^{J}\int
w^{(l)}_{(j,k)}(s,t)H^{(l-1)}_{(j)}(t)dt\Big{)}\\\
&=\sigma\Big{(}\sum_{b=1}^{B}b^{(l)}_{(k)}v^{*}_{b}(s)+\sum_{j=1}^{J}\sum_{c=1}^{C}\sum_{d=1}^{D}w^{(l)}_{(j,k),c,d}v_{c}(s)\int
v_{d}(t)H^{(l-1)}_{(j)}(t)dt\Big{)}\\\
&=\sigma\Big{(}\sum_{b=1}^{B}b^{(l)}_{(k)}v^{*}_{b}(s)+\sum_{j=1}^{J}\sum_{c=1}^{C}\sum_{d=1}^{D}w^{(l)}_{(j,k)c,d}v_{c}(s)A_{(j),d}^{(l)}\Big{)}\\\
&=\sigma\Big{(}B^{(l)}_{(k)}v^{*}(s)+\sum_{j=1}^{J}W^{(l)}_{(j,k)}v(s)A_{(j)}^{(l)}\Big{)}\\\
\end{split}$ (13)

$\begin{split}\widehat{Y}=H^{(L+1)}&=\sigma\Big{(}b^{(L+1)}+\sum_{j=1}^{J}\int
w^{(L+1)}_{(j)}(t)H^{(L)}_{(j)}(t)dt\Big{)}\\\
&=\sigma\Big{(}b^{(L+1)}+\sum_{j=1}^{J}\sum_{c=1}^{C}w^{(L+1)}_{(j),c}\int
v_{c}(t)H^{(L)}_{(j)}(t)dt\Big{)}\\\
&=\sigma\Big{(}b^{(L+1)}+\sum_{j=1}^{J}\sum_{c=1}^{C}w^{(L+1)}_{(j),c}A_{(j),c}^{(L+1)}\Big{)}=\sigma\Big{(}b^{(L+1)}+\sum_{j=1}^{J}W^{(L+1)}_{(j)}A_{(j)}^{(L+1)}\Big{)}\end{split}$
(14)

where $l=1,2,3,...,L$, $H^{(0)}_{(j)}(t):=X_{j}(t)$, and
$b^{(L+1)}\in\mathbb{R}$ is the unknown intercept. $B^{(l)}_{(k)}$,
$W^{(l)}_{(j,k)}$ and $W^{(L+1)}_{(j)}$ are the unknown basis weights
corresponding to their respective basis functions $v^{*}(s),$ $v(s),$ and
$v(t)$. $B$, $C$ and $D$ specify the number of basis used in each of the
continuous neuron of the continuous hidden layers. $\sigma(\cdot)$ is a non-
linear activation function.

For FBNN, we can again extend to as many continuous hidden layers as we want
with any number of continuous neurons. The activation function selection is
the same as FDNN. We can consider values of grid/domain ($s$), the number of
continuous hidden layers, the number of continuous neurons in each continuous
hidden layer, the selection of the basis functions, and the number of basis as
hyperparameters. We learn using the basis weights of the basis functions, to
reach from the functional predictors to the scalar response via a chain of
multiple CAMs in the forward propagation phase. The partial derivatives
required for the parameters of the continuous neurons in the continuous hidden
layers under FBNN in our proposed architecture (as seen in Figure 2) are as
follows:

$\begin{split}\frac{\partial H^{(l)}_{(j)}(s)}{\partial
B^{(l)}_{(k)}}&=\sigma^{\prime}\Big{(}B^{(l)}_{(k)}v^{*}(s)+\sum_{j=1}^{J}W^{(l)}_{(j,k)}v(s)A_{(j)}^{(l)}\Big{)}\times
v^{*}(s)\\\ \frac{\partial H^{(l)}_{(j)}(s)}{\partial
W^{(l)}_{(j,k)}}&=\sigma^{\prime}\Big{(}B^{(l)}_{(k)}v^{*}(s)+\sum_{j^{\prime}=1}^{J}W^{(l)}_{(j^{\prime},k)}v(s)A_{(j^{\prime})}^{(l)}\Big{)}\times
v(s)A_{(j)}^{(l)}\end{split}$ (15)

where l=1,2,3,…,L, $H^{(0)}_{(j)}(t):=X_{j}(t)$, and $\sigma^{\prime}(\cdot)$
represents the first derivative of $\sigma(\cdot)$.

Next, the partial derivatives required for the parameters of the functional
neurons in the functional output layer are,

$\begin{split}\frac{\partial H^{(L+1)}}{\partial
b^{(L+1)}}&=\sigma^{\prime}\Big{(}b^{(L+1)}+\sum_{j=1}^{J}W^{(L+1)}_{(j)}A_{(j)}^{(L+1)}\Big{)}\\\
\frac{\partial H^{(L+1)}}{\partial
W^{(L+1)}_{(j)}}&=\sigma^{\prime}\Big{(}b^{(L+1)}+\sum_{j^{\prime}=1}^{J}W^{(L+1)}_{(j^{\prime})}A_{(j^{\prime})}^{(L+1)}\Big{)}\times
A_{(j)}^{(L+1)}\end{split}$ (16)

The general way of optimizing for the parameters is the same in FBNN, using
backpropagation, but the gradients are found using the above partial
derivatives.

Regularization: Over-fitting is a common problem with machine learning and
nonparametric statistical methods. They are often over-parameterized and thus,
care must be taken to ensure that the results still produce reliable out of
sample predictions. Regularization in neural network to deal with over-fitting
is an active area of research (Yang et al.,, 2020; Ludwig et al.,, 2014;
Ayinde et al.,, 2019). Common methods used to address this issue include
penalization, drop out, data augmentation, and early stopping method. We
concentrate on early stopping, which utilizes a validation set to stop the
model fitting process early before it starts over-fitting. Naturally, both
FDNN and FBNN tend to over-fit if the models are too deep or run for a large
number of iterations. Early stopping is easy to implement and, as we will see
in Section 3, it is effective at improving the performance of FDNN and FBNN.

##  3 Empirical Results

In this section, we present multiple simulations to evaluate the effectiveness
of FDNN and FBNN, as well as benchmarking with multiple real world data sets.
We demonstrate the following objectives through our results: 1) FDNN and FBNN
are effectively trained by the derived optimization method 2) FDNN and FBNN
capture the relation between the functional predictors and the scalar response
3) FBNN and FDNN outperform compared to other methods.

###  3.1 Simulation Study

For simplicity, we consider a single predictor function (that is R=1) which is
dense and regular. We simulate $n=1500$ iid random curves
$\left\\{X_{1}(t),\cdots,X_{n}(t)\right\\}$ from a Gaussian process with mean
0 and covariance

$C_{X}(t,s)=\frac{\sigma^{2}}{\Gamma(\nu)2^{\nu-1}}\left(\frac{\sqrt{2\nu}|t-s|}{\rho}\right)^{\nu}K_{\nu}\left(\frac{\sqrt{2\nu}|t-s|}{\rho}\right),$

which is the Matérn covariance function and $K_{\nu}$ is the modified Bessel
function of the second kind. We set $\rho=0.5,$ $\nu=5/2$ and $\sigma^{2}=1.$
These curves are evaluated at $m=200$ equally-spaced time points from $[0,1]$.

The data generating models for the continuous and binary response are given in
Tables 1 and 3 respectively. We add $\epsilon_{i}$ (random noise) to our
continuous response model to adjust the signal to noise ratio. For the binary
response, we assume that we have an outcome $Y_{i}$ related to function
$X_{i}$ via a link function $g:$ $E(Y_{i}|X_{i})=g^{-1}\left(\eta_{i}\right)$
where $\eta_{i}$ is defined as the response in any of the models mentioned in
Section (2.2). Since we have a binary response, the link function used is the
logistic function, which is also the canonical link function. Each of these
settings is simulated 100 times and we divide the 1500 samples into 1000 for
training and 500 for testing. For the Deep Learning methods, we use early
stopping which requires a validation set. We use 500 samples from the training
set as the validation set in such cases. We are primarily interested in the
prediction performance and therefore, we report the Root Mean Square Error
(RMSE) of the prediction for the continuous response, and the classification
error as well as the log-likelihood for the binary response.

For FDNN and FBNN, we consider multiple combinations of continuous hidden
layers (l=1,2) and continuous neurons (j=1,2,4). The value for the grid (s) in
the continuous neurons is in the range (30,300). For FBNN, we use B-splines as
the basis functions and the number of basis used is in the range (5,10).

Simulation Scenarios
---
1\. Linear:
$Y_{i}=\alpha+\sum_{r=1}^{R}\int\beta_{r}(t)X_{i,r}(t)dt+\epsilon_{i}\qquad\quad\text{
where }\beta(t)=5\times\sin(2\pi t)\text{ and }\alpha=0$
2\. CAM: $Y_{i}=\sum_{r=1}^{R}\int
f(X_{i,r}(t),t)dt+\epsilon_{i}\qquad\quad\quad\hskip 11.09654pt\text{ where
}f(X_{i,r}(t),t)=(X_{i,r}(t)^{2})$
3\. Single Index:$\text{
}Y_{i}=g\left(\sum_{r=1}^{R}<\beta_{r},X_{i,r}>\right)+\epsilon_{i}\quad\hskip
6.25958pt\text{ where }g(a)=a^{2},\text{ }\beta(t)=5\times\sin(2\pi t)$
4\. Multiple Index:$\text{
}Y_{i}=g\left(\sum_{r=1}^{R}<\beta_{1,r},X_{i,r}>,...,\sum_{r=1}^{R}<\beta_{P,r},X_{i,r}>\right)+\epsilon_{i}\quad$
$\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\hskip
2.84526pt\text{ where }g(a,b)=a^{2}+b^{2},\text{
}\beta_{1}(t)=5\times\sin(2\pi t)$,
$\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\quad\hskip
3.55658pt\beta_{2}(t)=5\times\sin(3\pi t)$
5\. Quadratic:$\text{
}Y_{i}=\sum_{r=1}^{R}\int\beta_{r}(t)X_{i,r}(t)dt+\sum_{r=1}^{R}\int\int\beta(s,t)X_{i,r}(t)X_{i,r}(s)dtds+\epsilon_{i}\qquad\qquad$
$\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\hskip
2.84526pt\text{ where }\beta(t,s)=\beta_{1}(t)\otimes\beta_{2}(s),\text{
}\beta_{1}(t)=5\times\sin(3\pi t),$
$\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\quad\hskip
4.2679pt\beta_{2}(s)=5\times\sin(\pi s),\text{ }\otimes\text{is tensor
products}$
6\. Complex Quadratic:$\text{ }Y_{i}=\sum_{r=1}^{R}\int
f(X_{i,r}(t),t)dt+\sum_{r=1}^{R}\int\int
f(X_{i,r}(t),X_{i,r}(s),s,t)dtds+\epsilon_{i}\qquad\qquad$
$\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\hskip
2.84526pt\text{ where }f(X_{i,r}(t),t)=(X_{i,r}(t)^{2}),$
$\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\quad\hskip
4.2679ptf(X_{i,r}(t),X_{i,r}(s),s,t)=(X_{i,r}(t)\times X_{i,r}(s))^{2}$

Table 1: Different simulation scenarios with continuous response and a single
predictor function.

In Table 2, we summarize the performance of each method under different
simulation settings with a continuous response. We observe that when the true
model is Linear, four models do equally well: FLM, CAM, and both of our
methods. When the truth is CAM, then our approaches along with CAM do equally
well. More substantial gains are realized as the true model becomes more
complex and interactions are considered. In such situations, our methods
dominate others. Interestingly, under the Complex Quadratic settings, all the
methods struggle substantially compared to our approach, but there is still
potential for improvement (say, through deeper networks) as the optimal RMSE
is one (since that is the variance of the error). The gain from early stopping
can be seen in Appendix B Table 6 and Appendix B Table 7. Table 4 shows a
similar trend with the binary response. When the truth is Logistic or CAM,
then FLM, CAM, and both of our methods perform equally well. Our methods
outperform others as the true model becomes more complex and interactions are
considered. We see similar results when comparing the log-likelihood values as
seen in Appendix B Table 8. Note that all the table values for the deep
learning approaches are obtained using early stopping and the values
highlighted in the tables are within two standard error from the lowest RMSE
or classification error for each mapping.

Mapping | FLM | CAM | NN | FNN | FDNN | FBNN
---|---|---|---|---|---|---
Linear | 1.011 | 1.015 | 1.185 | 1.150 | 1.019 | 1.011
CAM | 1.705 | 1.053 | 1.451 | 1.298 | 1.076 | 1.058
Single Index | 1.711 | 1.532 | 1.515 | 1.225 | 1.067 | 1.061
Multiple Index | 1.808 | 1.350 | 1.609 | 1.290 | 1.100 | 1.060
Quadratic | 1.646 | 1.211 | 1.480 | 1.273 | 1.092 | 1.092
Complex Quadratic | 1.807 | 1.414 | 1.905 | 1.632 | 1.232 | 1.202
Table 2: Comparison of the RMSE for the methods under different setting
between the functional predictors and the continuous response. Simulation
Scenarios
---
1\. Logistic:$\text{
logit}(\eta_{i})=\alpha+\sum_{r=1}^{R}\int\beta_{r}(t)X_{i,r}(t)dt\qquad\quad\quad\text{
where }\beta(t)=5\times\sin(2\pi t)\text{ and }\alpha=0$
2\. CAM:$\text{ logit}(\eta_{i})=\sum_{r=1}^{R}\int
f(X_{i,r}(t),t)dt\qquad\quad\quad\quad\hskip 15.649pt\text{ where
}f(X_{i,r}(t),t)=sin(X_{i,r}(t))$
3\. Single Index:$\text{
logit}(\eta_{i})=g\left(\sum_{r=1}^{R}<\beta_{r},X_{i,r}>\right)\quad\quad\hskip
9.95845pt\text{ where }g(a)=sin(a),\text{ }\beta(t)=5\times\sin(2\pi t)$
4\. Multiple Index:$\text{
logit}(\eta_{i})=g\left(\sum_{r=1}^{R}<\beta_{1,r},X_{i,r}>,...,\sum_{r=1}^{R}<\beta_{P,r},X_{i,r}>\right)\quad$
$\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\hskip
3.55658pt\text{ where }g(a,b)=sin(sin(a)+b),\text{
}\beta_{1}(t)=5\times\sin(2\pi t),$
$\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\hskip
4.2679pt\beta_{2}(t)=5\times\sin(3\pi t)$
5\. Quadratic:$\text{
logit}(\eta_{i})=\sum_{r=1}^{R}\int\beta_{r}(t)X_{i,r}(t)dt+\sum_{r=1}^{R}\int\int\beta(s,t)X_{i,r}(t)X_{i,r}(s)dtds\qquad\qquad$
$\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\hskip
3.55658pt\text{ where }\beta(t,s)=\beta_{1}(t)\otimes\beta_{2}(s),\text{
}\beta_{1}(t)=5\times\sin(3\pi t),$
$\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\hskip
4.2679pt\beta_{2}(s)=5\times\sin(\pi s),\text{ }\otimes\text{is tensor
products}$
6\. Complex Quadratic:$\text{ logit}(\eta_{i})=\sum_{r=1}^{R}\int
f(X_{i,r}(t),t)dt+\sum_{r=1}^{R}\int\int
f(X_{i,r}(t),X_{i,r}(s),s,t)dtds\qquad\qquad$
$\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\hskip
3.55658pt\text{ where }f(X_{i,r}(t),t)=sin(X_{i,r}(t)),$
$\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\hskip
4.2679ptf(X_{i,r}(t),X_{i,r}(s),s,t)=sin(X_{i,r}(t)\times X_{i,r}(s))$

Table 3: Different simulation scenario with binary response and a single predictor function. Mapping | FLM | CAM | NN | FNN | FDNN | FBNN
---|---|---|---|---|---|---
Logistic | 0.276 | 0.276 | 0.291 | 0.288 | 0.278 | 0.277
CAM | 0.390 | 0.392 | 0.402 | 0.398 | 0.393 | 0.392
Single Index | 0.365 | 0.366 | 0.373 | 0.372 | 0.357 | 0.355
Multiple Index | 0.374 | 0.376 | 0.392 | 0.410 | 0.370 | 0.368
Quadratic | 0.363 | 0.328 | 0.344 | 0.338 | 0.322 | 0.324
Complex Quadratic | 0.394 | 0.392 | 0.403 | 0.398 | 0.392 | 0.393
Table 4: Comparison of the classification error for the methods under
different setting between the functional predictors and the binary response.

###  3.2 Real Data

We analyze three different data sets as seen in Figure 3. The fat content
spectrometric data (Ferraty and Vieu,, 2006) consists of food samples that
contain finely chopped pure meat with different fat content, the growth curves
from Berkeley Growth Study (Tuddenham and Snyder,, 1954) contains information
about the heights (in cm) of the children at different age and the speech
recognition data from TIMIT (available at http://statweb.stanford.edu/
tibs/Ele-mStatLearn/) where we have speech signals for different phonemes. The
spectrometric data is recorded on a Tecator Infractec Food and Feed Analyzer
working in the wavelength range 850-1050nm by the near-infrared (NIR)
transmission principle. For each meat sample, the data consists of a 100
channel spectrum of absorbances (-log10 of the transmittance measured by a
spectrometer) and a fat content measure derived from an analytic chemical
processing. We predict the fat content of a meat sample on the basis of its
NIR absorbance spectrum. The Growth data contains the heights of 39 boys and
54 girls from age 1 to 18 and the ages at which they were collected. The goal
is to build a classifier for the gender using the growth curves. In the speech
recognition data, the predictor is a voice signal as a function of time and
the response is the name of the vowel sounds pronounced, which can be regarded
as a discrete random variable. We only take a subset of the data, similar to
Li and Song, (2017) and only consider two phonemes (response) transcribed as
follows: “aa” as the vowel in “dark” and “ao” as the first vowel in “water”.
The distinction between the two groups is not straightforward as seen in
Figure 3, which shows the phoneme curves for each of the two groups computed
by a log-periodogram of 150 (Epifanio,, 2008). We use the speech signals as
functional predictors to classify the different phonemes.

Figure 3: Sample curves for each of the three data set. Data | FLM | CAM | NN | FNN | FDNN | FBNN
---|---|---|---|---|---|---
Spectrometric | 2.517 | 5.103 | 1.249 | 2.777 | 0.822 | 1.197
Growth | 0.150 | 0.150 | 0.450 | 0.450 | 0.100 | 0.100
Phonemes | 0.187 | 0.231 | 0.237 | 0.232 | 0.156 | 0.160
Table 5: Comparison of different methods on the real data using RMSE (Fat
Spectrum) and Classification Error (Growth and Phonemes).

In the spectrometric data, we have 216 samples of which 176 are for training
(40 samples for validation) and 40 samples are kept aside for testing. The
response is the fat content value. It should be noted that the spectra are
finely sampled, leading to very smooth curves as seen in Figure 3. Next, we
look into the growth curves from Berkeley Growth Study. The heights (in cm) of
the children in the data set are measured from 1 to 18 years at 31 (m=31)
timepoints. We divide the data into 73 and 20 samples for the training (20
samples for validation) and testing respectively. We take the response as
Gender. In the last data set, we divide the data into 640 and 160 samples for
training (160 samples for validation) and testing respectively and the
response is the phoneme. We can see from Table 5 that FDNN and FBNN are
performing much better compared to all the traditional methods in the
Spectrometric data for predicting the fat content of the meat sample. For the
two classification problems, we see a similar trend where our methods have
lower classification error compared to other methods whether the sample size
is small or large. The true nature of the relation between the functional
predictors and the response is never known in real world settings but our
methods can learn any kinds of patterns and excel in such situations as seen
in Table 5.

##  4 Discussion and Conclusions

In this paper, a novel approach for complex functional regression with a
scalar response has been developed. It is a natural extension of a neural
network to functional data. We use continuous hidden layers with continuous
neurons, allowing us to capture temporal dependency and complicated
relationships. We have devised two ways to implement this approach, FDNN and
FBNN. Our approach is flexible as we don’t assume a particular structure that
limits the kinds of non-linear patterns learned, especially, interactions
within time points for the predictors.

Understanding both the methods further will help us determine when it is more
beneficial to use FDNN or FBNN. Both, in general, seem to be working on par
with the best existing methods under simple cases but are outperforming in
more complicated non-linear settings. The proposed approach is also appealing
in practice as it provides a way for researchers to systematically conduct
deep learning tasks when they have functional predictors. The main advantage
of FDNN over FBNN is that it is straightforward to apply on dense and nicely
gridded data, works on the same ideas as a neural network and requires less
tuning, so it is a good approach to start with. FBNN is more flexible as it
can deal with irregular samples but is dependent on the imputation quality as
well as the selection of the basis functions. FBNN learns basis weights
allowing for greater parsimony with very dense grids. Further theoretical
investigation of these approaches is needed. Hence, an avenue for future work
is developing statistical theory to give these methods sound backing. Also,
the performance of FBNN in a sparse setting can be of interest but needs
nontrivial extensions.

Lastly, Functional regression has been studied extensively (especially the
scalar-on-function models) and non-linear modeling for functional data has
received more attention recently. But there are still many questions that
remain open, in particular function-on-function models. The applicability of
such models is of huge significance as they can be used in many areas
including time series, speech analysis, and computer vision. There is a need
for a more sophisticated and flexible function-on-function model in such
domains, as functional data becomes more prevalent in today’s world. Hence,
extending our work to the function-on-function case will enable us to model
functional predictors with functional response and help capture their complex
relation.

## Acknowledgement

This research was supported in part by the following grant to Pennsylvania
State University: NSF SES-1853209

## References

  * Asadi et al., (2014) Asadi, E., Silva, M. G., Antunes, C. H., Dias, L., and Glicksman, L. (2014).  Multi-objective optimization for building retrofit: A model using genetic algorithm and artificial neural network and an application.  Energy and Buildings, 81:444 – 456.
  * Ayinde et al., (2019) Ayinde, B. O., Inanc, T., and Zurada, J. M. (2019).  Regularizing deep neural networks by enhancing diversity in feature extraction.  IEEE Transactions on Neural Networks and Learning Systems, 30(9):2650–2661.
  * Chan et al., (2016) Chan, W., Jaitly, N., Le, Q., and Vinyals, O. (2016).  Listen, attend and spell: A neural network for large vocabulary conversational speech recognition.  In 2016 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 4960–4964.
  * Chen et al., (2012) Chen, D., Hall, P., and Müller, H. (2012).  Single and multiple index functional regression models with nonparametric link.  Annals of Statistics - ANN STATIST, 39.
  * Chen and Huang, (2017) Chen, Z. and Huang, X. (2017).  End-to-end learning for lane keeping of self-driving cars.  In 2017 IEEE Intelligent Vehicles Symposium (IV), pages 1856–1860.
  * Conan-Guez and Rossi, (2002) Conan-Guez, B. and Rossi, F. (2002).  Multi-layer perceptrons for functional data analysis: A projection based approach.  In ICANN.
  * Dorffner, (1996) Dorffner, G. (1996).  Neural networks for time series processing.  Neural Network World, 6:447–468.
  * Eilers et al., (2009) Eilers, P. H., Li, B., and Marx, B. D. (2009).  Multivariate calibration with single-index signal regression.  Chemometrics and Intelligent Laboratory Systems, 96(2):196 – 202\.  Chimiometrie 2007, Lyon, France, 29-30 November 2007.
  * Epifanio, (2008) Epifanio, I. (2008).  Shape descriptors for classification of functional data.  Technometrics, 50(3):284–294.
  * Fan et al., (2015) Fan, Y., James, G. M., and Radchenko, P. (2015).  Functional additive regression.  Ann. Statist., 43(5):2296–2325.
  * Ferraty, (2011) Ferraty, F. (2011).  Recent Advances in Functional Data Analysis and Related Topics.
  * Ferraty et al., (2013) Ferraty, F., Goia, A., Salinelli, E., and Vieu, P. (2013).  Functional projection pursuit regression.  TEST, 22.
  * Ferraty and Romain, (2011) Ferraty, F. and Romain, Y. (2011).  The Oxford Handbook of Functional Data Analysis.
  * Ferraty and Vieu, (2006) Ferraty, F. and Vieu, P. (2006).  Nonparametric Functional Data Analysis: Theory and Practice (Springer Series in Statistics).  Springer-Verlag, Berlin, Heidelberg.
  * Geenens, (2011) Geenens, G. (2011).  Curse of dimensionality and related issues in nonparametric functional regression.  Statistics Surveys, 5.
  * Graves et al., (2013) Graves, A., Mohamed, A., and Hinton, G. (2013).  Speech recognition with deep recurrent neural networks.  In 2013 IEEE International Conference on Acoustics, Speech and Signal Processing, pages 6645–6649.
  * Gregor et al., (2015) Gregor, K., Danihelka, I., Graves, A., Rezende, D. J., and Wierstra, D. (2015).  Draw: A recurrent neural network for image generation.
  * Horváth and Kokoszka, (2012) Horváth, L. and Kokoszka, P. (2012).  Inference for Functional Data with Applications.  Springer Series in Statistics. Springer New York.
  * James and Silverman, (2004) James, G. and Silverman, B. (2004).  Functional adaptive model estimation.  Journal of the American Statistical Association, 100.
  * Joshi et al., (2020) Joshi, A., Khosravy, M., and Gupta, N. (2020).  Machine Learning for Predictive Analysis: Proceedings of ICTIS 2020.  Lecture Notes in Networks and Systems. Springer Singapore.
  * Kanellopoulos and Wilkinson, (1997) Kanellopoulos, I. and Wilkinson, G. G. (1997).  Strategies and best practice for neural network image classification.  International Journal of Remote Sensing, 18(4):711–725.
  * Karpathy et al., (2014) Karpathy, A., Toderici, G., Shetty, S., Leung, T., Sukthankar, R., and Fei-Fei, L. (2014).  Large-scale video classification with convolutional neural networks.  In Proceedings of the IEEE conference on Computer Vision and Pattern Recognition, pages 1725–1732.
  * Khashei and Bijari, (2010) Khashei, M. and Bijari, M. (2010).  An artificial neural network (p,d,q) model for timeseries forecasting.  Expert Systems with Applications, 37(1):479 – 489.
  * Kokoszka and Reimherr, (2018) Kokoszka, P. and Reimherr, M. (2018).  Introduction to Functional Data Analysis.  New York: Chapman and Hall/CRC.
  * Li and Song, (2017) Li, B. and Song, J. (2017).  Nonlinear sufficient dimension reduction for functional data.  Ann. Statist., 45(3):1059–1095.
  * Ludwig et al., (2014) Ludwig, O., Nunes, U., and Araujo, R. (2014).  Eigenvalue decay: A new method for neural network regularization.  Neurocomputing, 124:33 – 42.
  * Ma and Zhu, (2016) Ma, H. and Zhu, Z. (2016).  Continuously dynamic additive models for functional data.  Journal of Multivariate Analysis, 150:1 – 13.
  * Maller et al., (2013) Maller, H., Wu, Y., and Yao, F. (2013).  Continuously additive models for nonlinear functional regression.  Biometrika, 100(3):607–622.
  * Maqueda et al., (2018) Maqueda, A. I., Loquercio, A., Gallego, G., García, N., and Scaramuzza, D. (2018).  Event-based vision meets deep learning on steering prediction for self-driving cars.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR).
  * McLean et al., (2014) McLean, M. W., Hooker, G., Staicu, A., Scheipl, F., and Ruppert, D. (2014).  Functional generalized additive models.  Journal of Computational and Graphical Statistics, 23(1):249–269.
  * Morchid, (2018) Morchid, M. (2018).  Parsimonious memory unit for recurrent neural networks with application to natural language processing.  Neurocomputing, 314:48 – 64.
  * Morris, (2014) Morris, J. (2014).  Functional regression.  Annual Review of Statistics and Its Application, 2.
  * Muller and Yao, (2008) Muller, H. and Yao, F. (2008).  Functional additive models.  Journal of the American Statistical Association, 103:1534–1544.
  * Olver, (2019) Olver, P. J. (2019).  Introduction to the calculus of variations.
  * Qin et al., (2017) Qin, Y., Song, D., Chen, H., Cheng, W., Jiang, G., and Cottrell, G. (2017).  A dual-stage attention-based recurrent neural network for time series prediction.
  * Ramsay and Silverman, (1997) Ramsay, J. and Silverman, B. (1997).  Functional Data Analysis.  Springer series in statistics. Springer.
  * Rao et al., (2020) Rao, A. R., Wang, Q., Wang, H., Khorasgani, H., and Gupta, C. (2020).  Spatio-temporal functional neural networks.  In 2020 IEEE 7th International Conference on Data Science and Advanced Analytics (DSAA), pages 81–89.
  * Reimherr et al., (2017) Reimherr, M., Sriperumbudur, B., and Taoufik, B. (2017).  Optimal prediction for additive function-on-function regression.
  * Reiss et al., (2016) Reiss, P., Goldsmith, J., Shang, H. L., and Ogden, R. (2016).  Methods for scalar‐on‐function regression.  International Statistical Review, 85.
  * Rossi and Conan-Guez, (2006) Rossi, F. and Conan-Guez, B. (2006).  Theoretical properties of projection based multilayer perceptrons with functional inputs.  Neural Processing Letters, 23.
  * Rossi et al., (2002) Rossi, F., Conan-Guez, B., and Fleuret, F. (2002).  Functional data analysis with multi layer perceptrons.  In Proceedings of the 2002 International Joint Conference on Neural Networks. IJCNN’02 (Cat. No.02CH37290), volume 3, pages 2843–2848 vol.3.
  * Rumelhart et al., (1986) Rumelhart, D. E., Hinton, G. E., and Williams, R. J. (1986).  Learning Internal Representations by Error Propagation, page 318–362.  MIT Press, Cambridge, MA, USA.
  * Seide et al., (2011) Seide, F., Li, G., and Yu, D. (2011).  Conversational speech transcription using context-dependent deep neural networks.  pages 437–440.
  * Shi et al., (2016) Shi, W., Caballero, J., Huszar, F., Totz, J., Aitken, A. P., Bishop, R., Rueckert, D., and Wang, Z. (2016).  Real-time single image and video super-resolution using an efficient sub-pixel convolutional neural network.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR).
  * Singh and Singh, (2017) Singh, A. and Singh, K. K. (2017).  Satellite image classification using genetic algorithm trained radial basis function neural network, application to the detection of flooded areas.  Journal of Visual Communication and Image Representation, 42:173 – 182.
  * Stinchcombe, (1999) Stinchcombe, M. (1999).  Neural network approximation of continuous functionals and continuous functions on compactifications.  Neural Networks, 12(3):467 – 477.
  * Stoker, (1986) Stoker, T. M. (1986).  Consistent estimation of scaled coefficients.  Econometrica, 54(6):1461–1481.
  * Tuddenham and Snyder, (1954) Tuddenham, R. and Snyder, M. (1954).  Physical growth of california boys and girls from birth to eighteen years.  Publications in child development. University of California, Berkeley, 1:183–364.
  * Wang et al., (2020) Wang, Q., Wang, H., Gupta, C., Rao, A. R., and Khorasgani, H. (2020).  A non-linear function-on-function model for regression with time series data.
  * Wang et al., (2019) Wang, Q., Zheng, S., Farahat, A., Serita, S., and Gupta, C. (2019).  Remaining useful life estimation using functional data analysis.
  * Wang et al., (2019) Wang, Q., Zheng, S., Farahat, A., Serita, S., Saeki, T., and Gupta, C. (2019).  Multilayer perceptron for sparse functional data.  In 2019 International Joint Conference on Neural Networks (IJCNN), pages 1–10.
  * Wang and Yang, (2019) Wang, W. and Yang, Y. (2019).  Development of convolutional neural network and its application in image classification: A survey.  Optical Engineering, 58:1.
  * Wang and Ruppert, (2013) Wang, X. and Ruppert, D. (2013).  Optimal prediction in an additive functional model.  Statistica Sinica, 25.
  * Yang et al., (2020) Yang, Y., Yuan, Y., Chatzimichailidis, A., Sloun, R., Lei, L., and Chatzinotas, S. (2020).  Proxsgd: Training structured neural networks under regularization and constraints.  In International Conference on Learning Representations.
  * Yao and Muller, (2010) Yao, F. and Muller, H. (2010).  Functional quadratic regression.  Biometrika, 97(1):49–64.
  * Yue-Hei Ng et al., (2015) Yue-Hei Ng, J., Hausknecht, M., Vijayanarasimhan, S., Vinyals, O., Monga, R., and Toderici, G. (2015).  Beyond short snippets: Deep networks for video classification.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 4694–4702.
  * Zhang et al., (2017) Zhang, Z., Wang, H., Xu, F., and Jin, Y. (2017).  Complex-valued convolutional neural network and its application in polarimetric sar image classification.  IEEE Transactions on Geoscience and Remote Sensing, 55(12):7177–7188.
  * Zhu et al., (2013) Zhu, H., Yao, F., and Zhang, H. (2013).  Structured functional additive regression in reproducing kernel hilbert spaces.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 76.

##  5 Appendix

## Appendix A. Mathematical Details for FDNN and FBNN

The partial derivatives used in FDNN are as follows,

$\begin{split}\frac{\partial H^{(l)}_{(k)}}{\partial
w^{(l)}_{(j,k)}}(s,t)&=\sigma^{\prime}\Big{(}b^{(l)}_{(k)}(s)+\sum_{j^{\prime}=1}^{J}\int
w^{(l)}_{(j^{\prime},k)}(s,t^{\prime})H^{(l-1)}_{(j^{\prime})}(t^{\prime})dt^{\prime}\Big{)}\\\
&\times\int\frac{\partial}{\partial
w^{(l)}_{(j,k)}(s,t)}w^{(l)}_{(j,k)}(s,t)H^{(l-1)}_{(j)}(t)dt\\\
&=\sigma^{\prime}\Big{(}b^{(l)}_{(k)}(s)+\sum_{j^{\prime}=1}^{J}\int
w^{(l)}_{(j^{\prime},k)}(s,t^{\prime})H^{(l-1)}_{(j^{\prime})}(t^{\prime})dt^{\prime}\Big{)}\times
H^{(l-1)}_{(j)}(t)\\\ \frac{\partial H^{(L+1)}}{\partial
w^{(L+1)}_{(j)}}(t)&=\sigma^{\prime}\Big{(}b^{(L+1)}+\sum_{j^{\prime}=1}^{J}\int
w^{(L+1)}_{(j^{\prime})}(t^{\prime})H^{(L)}_{(j^{\prime})}(t^{\prime})dt^{\prime}\Big{)}\\\
&\times\int\frac{\partial}{\partial
w^{(L+1)}_{(j)}(t)}w^{(L+1)}_{(j)}(t)H^{(L)}_{(j)}(t)dt\\\
&=\sigma^{\prime}\Big{(}b^{(L+1)}+\sum_{j^{\prime}=1}^{J}\int
w^{(L+1)}_{(j^{\prime})}(t^{\prime})H^{(L)}_{(j^{\prime})}(t^{\prime})dt^{\prime}\Big{)}\times
H^{(L)}_{(j)}(t)\end{split}$ (17)

where $\sigma^{\prime}(\cdot)$ represents the first derivative of the
activation function $\sigma(\cdot)$.

## Appendix B. Early Stopping and Log-likelihood

Mapping | FDNN | FBNN | FDNN (2,2) | FDNN (4,4) | FBNN (2,2) | FBNN (4,4)
---|---|---|---|---|---|---
Linear | $\mathbf{1.069}$ | ${1.015}$ | 1.452 | 1.440 | 1.244 | 1.271
CAM | $\mathbf{1.118}$ | ${1.084}$ | 1.460 | 1.459 | 1.217 | 1.202
Single Index | $\mathbf{1.151}$ | ${1.359}$ | 1.457 | 1.461 | 1.243 | 1.200
Multiple Index | ${1.154}$ | $\mathbf{1.091}$ | 1.495 | 1.473 | 1.207 | 1.219
Quadratic Linear | $\mathbf{1.136}$ | ${1.205}$ | 1.495 | 1.465 | 1.284 | 1.263
Complex Quadratic | $\mathbf{1.287}$ | ${1.345}$ | 1.586 | 1.591 | 1.386 | 1.309
Table 6: Comparison of Deep FDNN and FBNN under different settings between the functional predictors and the continuous response. Mapping | FDNN | FBNN | FDNN (2,2) | FDNN (4,4) | FBNN (2,2) | FBNN (4,4)
---|---|---|---|---|---|---
Linear | $\mathbf{1.019}$ | $\mathbf{1.011}$ | 1.036 | 1.043 | 1.034 | 1.028
CAM | 1.076 | $\mathbf{1.062}$ | 1.107 | 1.072 | 1.105 | $\mathbf{1.058}$
Single Index | $\mathbf{1.067}$ | 1.346 | 1.069 | 1.122 | 1.078 | $\mathbf{1.061}$
Multiple Index | 1.100 | 1.089 | 1.120 | 1.092 | 1.089 | $\mathbf{1.060}$
Quadratic Linear | $\mathbf{1.092}$ | 1.184 | 1.121 | 1.112 | 1.119 | $\mathbf{1.092}$
Complex Quadratic | 1.232 | 1.337 | 1.283 | 1.301 | 1.277 | $\mathbf{1.202}$
Table 7: Comparison of Deep FDNN and FBNN with early stopping under different
settings between the functional predictors and the continuous response.

Note: (2,2) means two continuous hidden layer with two continuous neurons each
and (4,4) means two continuous hidden layer with four continuous neurons each.

Mapping | FLM | CAM | NN | FNN | FDNN | FBNN
---|---|---|---|---|---|---
Logistic | 0.544 | 0.544 | 0.591 | 0.562 | 0.548 | 0.546
CAM | 0.663 | 0.662 | 0.698 | 0.693 | 0.665 | 0.663
Single Index | 0.665 | 0.666 | 0.680 | 0.673 | 0.663 | 0.662
Multiple Index | 0.661 | 0.660 | 0.685 | 0.668 | 0.658 | 0.658
Quadratic | 0.639 | 0.595 | 0.652 | 0.615 | 0.584 | 0.596
Complex Quadratic | 0.664 | 0.663 | 0.693 | 0.670 | 0.661 | 0.661
Table 8: Comparison of log-likelihood for the methods under different settings
between the functional predictors and the binary response.


# Two-photon interference from a quantum emitter in hexagonal boron nitride

Clarisse Fournier1, Sébastien Roux1,2, Kenji Watanabe3, Takashi Taniguchi4,
Stéphanie Buil1, Julien Barjon1, Jean-Pierre Hermier1, Aymeric Delteil1 1
Université Paris-Saclay, UVSQ, CNRS, GEMaC, 78000, Versailles, France.
2 Université Paris-Saclay, ONERA, CNRS, Laboratoire d’étude des
microstructures, 92322, Châtillon, France.
3 Research Center for Functional Materials, National Institute for Materials
Science, 1-1 Namiki, Tsukuba 305-0044, Japan
4 International Center for Materials Nanoarchitectonics, National Institute
for Materials Science, 1-1 Namiki, Tsukuba 305-0044, Japan
<EMAIL_ADDRESS>
###### Abstract

Recently discovered quantum emitters in two-dimensional (2D) materials have
opened new perspectives of integrated photonic devices for quantum
information. Most of these applications require the emitted photons to be
indistinguishable, which has remained elusive in 2D materials. Here, we
investigate two-photon interference of a quantum emitter generated in
hexagonal boron nitride (hBN) using an electron beam. We measure the
correlations of zero-phonon-line photons in a Hong-Ou-Mandel (HOM)
interferometer under non-resonant excitation. We find that the emitted photons
exhibit a partial indistinguishability of $0.44\pm 0.11$ in a 3 ns time
window, which corresponds to a corrected value of $0.56\pm 0.11$ after
accounting for imperfect emitter purity. The dependence of the HOM visibility
on the width of the post-selection time window allows us to estimate the
dephasing time of the emitter to be $\sim 1.5$ ns, about half the limit set by
spontaneous emission. A visibility above 90 % is under reach using Purcell
effect with up-to-date 2D material photonics.

Two-photon interference is essential for many photonic implementations of
quantum information protocols, from linear optical quantum computing klm01  to
distant entanglement generation Cabrillo99 ; Bernien13 ; Delteil15  and
quantum communication Gisin07 . The indistinguishability of two single-photon
pulses – which quantifies their ability to interfere – results in the so-
called Hong-Ou-Mandel (HOM) effect hom87 . The latter refers to the fact that
perfectly indistinguishable photons simultaneously reaching the two input
ports of a beamsplitter always exit the beamsplitter from the same output port
Bouchard21 . Experimental observation of the HOM effect between consecutive
photons from a quantum emitter constitutes an important milestone in the use
of a physical system for the generation of scalable photonic qubits. Among the
physical systems able to generate indistinguishable photons, solid-state
single-photon emitters (SPEs) have been widely investigated due to their
potential for integration in photonic devices Aharonovich16 . Thus, photon
indistinguishability has been experimentally demonstrated with III-V
semiconductor quantum dots Santori02 ; Somaschi16 ; Thoma17 ; Gao13  and color
centers in three-dimensional wide bandgap crystals Sipahigil12 ; Sipahigil14 ;
Morioka20 .

In turn, recently discovered quantum emitters in 2D materials, comprising
trapped excitons in transition metal dichalcogenides chakraborty15 ; he15 ;
koperski15 ; srivastava15 ; tonndorf15  and color centers in hBN tran16 ;
bourrelier16 ; martinez16 , have raised a growing interest owing to the
perspectives of extreme miniaturization and integration into complex
heterostructures Geim13  – yet without demonstration of two-photon
interference to date. Among these systems, a recently discovered family of hBN
SPEs stands out – a class of blue emitting color centers (abbreviated
B-centers in the following) that can be generated at controlled locations
using an electron beam. Their zero-phonon-line (ZPL) center wavelength is
consistently found within 3 meV around 436 nm shevitski19 ; Fournier21 ;
Gale22 . Several studies have already demonstrated their spectral stability,
narrow linewidth, brightness and single-photon emission up to room temperature
Fournier21 ; Gale22 ; Horder22 .

In this letter, we characterize two-photon interference of light emitted by an
individual B-center. Our sample consists of a single hBN crystal grown using
high pressure, high temperature conditions Taniguchi07 , that we exfoliated on
a SiO2(285 nm)/Si substrate. We generate a SPE ensemble in a commercial
scanning electron microscope (SEM) using a slightly defocused electron beam
(diameter $\sim$300 nm) under 15 kV acceleration voltage and 10 nA current,
following Fournier21  (figure 1a). We subsequently characterize the sample in
a confocal microscope operating in a helium closed-cycle cryostat, keeping the
sample at 4 K. The sample is optically excited by a pulsed diode laser of 405
nm wavelength, 850 $\mu$W power and 80 MHz repetition rate (figure 1b) that is
focused on the sample using a microscope objective of numerical aperture 0.8.
The SPE luminescence is collected through the same objective, and coupled into
a single-mode fiber. In the following, we focus on an individual SPE with a
ZPL centered at 436.24 nm. Figure 1c shows a low resolution spectrum of the
SPE, exhibiting the usual spectral shape of the B-centers, which comprises a
narrow ZPL (40 % of the emission) and an acoustic phonon sideband (60 %). We
ensure that spectral diffusion of the ZPL is limited, as shown figure 1d,
where the wavelength fluctuations are contained below 20 pm.

Figure 1: (a) Irradiation by a 15 kV electron beam generates B-centers in a
multilayer hBN crystal. (b) Energy levels of the SPE: A 405 nm laser excites
the emitter. Non-radiative relaxation occurs, followed by emission of a photon
at 436 nm. (c) Low resolution spectrum of the SPE, where the ZPL and the
acoustic phonon sideband can be observed. Inset: High resolution spectrum,
limited by the spectrometer resolution of about 100 $\mu$eV. The phonon
pedestal of the main panel is no more visible due to the higher resolution.
(d) High resolution spectra as a function of time measured with 2 s
integration time during 5 min, showing the stability of the SPE.
Figure 2: (a) Experimental setup for the HOM experiment. LCR: Liquid crystal
retarder; BS: Beamsplitter; APD: Avalanche photodiode; FPC: Fiber polarization
controller. TCSPC: time-correlated single photon counter. (b) Histogram of the
photon detection events at the output port. The instrument response function
is indicated by the gray shadow. The orange shading of width $\Delta t=3$ ns
indicates the post-selection time window. A fit to the data (not shown)
provides the spontaneous emission time $\Gamma_{\mathrm{sp}}^{-1}=1.9$ ns. (c)
Second-order correlation function of the photons emitted during $\Delta t$,
yielding $g^{(2)}(0)=0.14\pm 0.03$. The horizontal axis is in units of
repetition periods. (d) Two-photon coincidences in Hong-Ou-Mandel
configuration. The orange (blue) dots denote the normalized coincidence rate
for parallel (orthogonal) polarization configuration of the second
beamsplitter input ports. The center peak provides the values
$g^{(2)}_{\mathrm{HOM,\parallel}}(0)=0.32\pm 0.05$ (parallel case) and
$g^{(2)}_{\mathrm{HOM,\perp}}(0)=0.58\pm 0.07$ (orthogonal case). The light
gray bars marks the theoretical values for distinguishable photons of
$g^{(2)}(0)=0.14$. The dark gray bar indicates the theoretical value for fully
indistinguishable photons of $g^{(2)}(0)=0.14$.

Figure 2a depicts the experimental setup used for two-photon interference
characterization. The photoluminescence is collected in a single-mode fiber
that channels the photons to a delayed Mach-Zehnder interferometer. The ZPL is
filtered using a transmission grating of 1379 grooves per millimeter. Together
with subsequent coupling to single-mode fibers, it implements a narrow
bandpass filter of 100 GHz bandwidth. This spectral width is larger than the
time-averaged linewidth but much narrower than the width of the acoustic
phonon sideband of about 7.8 meV (1.9 THz). A polarizer ensures that the input
photons have a well-defined linear polarization at the input port of the
delayed Mach-Zehnder interferometer. One of the arms is delayed by the same
amount as our repetition period (12.5 ns) using a 2.6 m fiber, such that two
consecutively emitted photons can simultaneously impinge on the beamsplitter.
A liquid crystal retarder is inserted in the other arm to rotate the photon
polarization by 90 degrees when suited, allowing photon polarization at the
two input ports of the second beamsplitter to be either identical (parallel)
or orthogonal. The total count rate of the ZPL photons at the output ports is
about 1200 counts per second. Figure 2b shows a histogram of the photon
detection times at the output ports. We first consider photon detections
occurring during the $\Delta t=$ 3 ns time window highlighted by the orange
shadowing on figure 2b, that is located after the laser pulse of width 550 ps
(gray shadowing on figure 2b). Figure 2c shows the second-order photon
correlations measured in a Hanbury Brown and Twiss configuration of the
interferometer (i.e. in a single arm). The relative height of the center
period allows to infer the photon purity of $g^{(2)}_{\mathrm{HBT}}(0)=0.14\pm
0.03$, limited by reminiscent background signal and dark counts.

We then measure the photon coincidences in the HOM configuration of the
interferometer. Figure 2d gives the normalized coincidences measured during 36
hours while alternating between parallel and orthogonal polarizations using
the liquid crystal retarder, considering photons detected during the same time
window of width $\Delta t$. The significant reduction of the center period
value $g^{(2)}_{\mathrm{HOM,\parallel}}(0)=0.32\pm 0.05$ in the parallel
polarization case as compared with the orthogonal case
$g^{(2)}_{\mathrm{HOM,\perp}}(0)=0.58\pm 0.07$ is a signature of photon
coalescence. The raw (uncorrected) degree of indistinguishability of the
emitted photons is then given by the interference visibility defined as
$V_{\mathrm{HOM}}=1-g^{(2)}_{\mathrm{HOM,\parallel}}(0)/g^{(2)}_{\mathrm{HOM,\perp}}(0)$.
In our case, we find $V_{\mathrm{HOM}}=0.44\pm 0.11$. In the case of single
photons with ideal purity (i.e. $g^{(2)}_{\mathrm{HBT}}(0)=0$), this quantity
ranges between $V_{\mathrm{HOM}}=0$ (perfectly distinguishable photons) and
$V_{\mathrm{HOM}}=1$ (perfectly indistinguishable photons). When multiple
detection events are not negligible, the theoretical bounds of
$g^{(2)}_{\mathrm{HOM,\parallel}}(0)$ and $g^{(2)}_{\mathrm{HOM,\perp}}(0)$
are offset upwards Santori02 , such that the corrected visibility reads
$V_{\mathrm{HOM}}^{\mathrm{corr}}=\left(1+2g^{(2)}_{\mathrm{HBT}}(0)\right)V_{\mathrm{HOM}}$.
In the case of the studied SPE, the theoretical upper (resp. lower) bounds
accounting for our finite value of $g^{(2)}_{\mathrm{HBT}}(0)$ are indicated
by light (resp. dark) gray bars on the center period of the histogram shown
figure 2d. Accordingly, we find a corrected HOM visibility of
$V_{\mathrm{HOM}}^{\mathrm{corr}}=0.56\pm 0.11$. This number is comparable
with the indistinguishability of non-resonantly excited quantum dots Santori02
; Thoma17  or color centers in wide gap 3D semiconductors Morioka20 . This
observation of HOM interference from single photons emitted by a 2D material
quantum emitter constitutes the main result of our study, and suggests
possible practical applications of hBN for optical quantum information, upon
further improvement of the visibility.

The limited value of the corrected HOM visibility can be possibly attributed
to fast dephasing of the optical dipole. The total dephasing rate of an
optical transition can be written $\gamma=\Gamma_{\mathrm{sp}}/2+\gamma^{*}$,
where $\Gamma_{\mathrm{sp}}=1/T_{1}$ is the spontaneous emission rate, and
$\gamma^{*}=1/T_{2}^{*}$ denotes the rate of dephasing caused by reservoirs
other than the vacuum electromagnetic field. In the pulsed regime, only
$\gamma^{*}$ causes a reduction of the photon indistinguishability kiraz04 ;
bylander03 . An expected consequence of dephasing is that extending the
integration window $\Delta t$ would degrade the HOM visibility by allowing a
larger delay time between detected pairs Bouchard21 ; Martinez22 ; Proux15 ;
Kim16 ; Nawrath19 . It is then in principle possible to estimate $\gamma^{*}$
by observing the influence of the postselection time window on the
interference visibility. Figure 3a (blue dots) plots the measured HOM
visibility $V_{\mathrm{HOM}}^{\mathrm{corr}}$ as a function of the post-
selection window size $\Delta t$. We can observe an overall decrease of the
HOM visibility as $\Delta t$ increases, which is consistent with the effect of
dephasing. The decay of $V_{\mathrm{HOM}}^{\mathrm{corr}}(\Delta t)$ can be
fitted by an exponential function, of which we have fixed the intercept to be
1. This result is compatible with the assumption that finite purity and
dephasing are the main sources of imperfect HOM visibility. The fit yields the
decay time $\tau_{V}=2.0$ ns. To relate this timescale to an estimation of
$\gamma^{*}$, we simulate the two-photon interference visibility of a two-
level atom with dephasing based on numerical integration of a master equation
loudon ; Fischer16 ; qutip1 ; qutip2  accounting for spontaneous emission and
dephasing in the Lindblad form. Our simulation accounts for a time post-
selection of the coincidences in the same way as experimentally realized. The
details of the simulation are exposed in the Supplemental Material
supplemental . We compute the two-photon interference visibility while varying
the post-selection time window, from which we extract the visibility decay
timescale $\tau_{V}$. By repeating this process for different values of the
dephasing rate $\gamma^{*}$, we obtain a one-to-one relation between
$\tau_{V}$ and $\gamma^{*}$, which we plot on figure 3b. The dephasing rate
associated with our experimental data corresponds to $T_{2}^{*}=2.4\pm 0.7$
ns. This value is in agreement with prior estimations of B-center dephasing
time based on the decay of Rabi oscillations in resonant excitation Horder22 ,
where values of $T_{2}^{*}$ between 0.7 and more than 2 ns were estimated,
depending on the emitter and on the laser power. In turn, this places the
total dephasing time $T_{2}=(\Gamma_{\mathrm{sp}}/2+\gamma^{*})^{-1}=1.5$ ns,
about halfway from the perfect indistinguishability limit $T_{2}=2T_{1}=3.8$
ns.

Figure 3: (a) Blue dots: HOM visibility $V_{\mathrm{HOM}}^{\mathrm{corr}}$ as
a function of the post-selection time window width $\Delta t$ . Orange line:
exponential fit to the data, yielding $\tau_{V}=1.7$ ns. (b) Blue line:
Simulated value of $\tau_{V}$ as a function of the dephasing rate
$\gamma^{*}$, in units of radiative lifetime $\Gamma_{\mathrm{sp}}^{-1}$.
Orange dot: Value corresponding to the experimental data (a).

The value of indistinguishability we measured could be further improved by
using resonant excitation, which selectively addresses the transition of
interest. This leads to a decrease of the saturation power by several orders
of magnitude and therefore of the environment noise, yielding up to Fourier-
limited linewidths Fournier23 . Additionally, embedding the emitter in a
cavity would increase the ratio between the rates of spontaneous emission and
dephasing. We calculate that a Purcell factor of 7 (15) would increase the
indistinguishability to 80 % (90 %), while improving at the same time the
collection efficiency. Similar or higher spontaneous emission enhancements
have already been achieved in 2D materials Froch21 ; Iff21 . Integrating
B-centers into electrical gates for fine tuning of the optical resonance would
then allow interference between distinct emitters Zhigulin23 . The
indistinguishability of photons emitted by hBN B-centers, together with their
already-established advantageous photophysical properties and the possibility
to locate them at a prechosen position, opens the way to a broad range of
applications in integrated quantum photonics and optical quantum information
based on 2D materials.

##  I Acknowledgments

The authors acknowledge Aurélie Pierret and Michael Rosticher for the flake
exfoliation. This work is supported by the French Agence Nationale de la
Recherche (ANR) under reference ANR-21-CE47-0004-01 (E$-$SCAPE project). This
work also received funding from the European Union’s Horizon 2020 research and
innovation program under Grant No. 881603 (Graphene Flagship Core 3). K.W. and
T.T. acknowledge support from JSPS KAKENHI (Grant Numbers 19H05790, 20H00354
and 21H05233).

## References

  * (1) E. Knill, R. Laflamme, and G. J. Milburn, A scheme for efficient quantum computation with linear optics, Nature 409, 46 (2001).
  * (2) C. Cabrillo, J. I. Cirac, P. García-Fernández, and P. Zoller, Creation of entangled states of distant atoms by interference, Phys. Rev. A 59, 1025 (1999).
  * (3) H. Bernien, B. Hensen, W. Pfaff, G. Koolstra, M. S. Blok, L. Robledo, T. H. Taminiau, M. Markham, D. J. Twitchen, L. Childress, and R. Hanson, Heralded entanglement between solid-state qubits separated by three metres, Nature 497, 86 (2013).
  * (4) A. Delteil, Z. Sun, W.-B. Gao, E. Togan, S. Fält, and A. Imamoğlu, Generation of heralded entanglement between distant hole spins, Nature Phys. 12, 218 (2016).
  * (5) N. Gisin and R. Thew, Quantum communication, Nat. Photon. 1, 165 (2007).
  * (6) C. K. Hong, Z. Y. Ou, and L. Mandel, Measurement of subpicosecond time intervals between two photons by interference, Phys. Rev. Lett. 59, 2044 (1987).
  * (7) F. Bouchard, A. Sit, Y. Zhang, R. Fickler, F. M. Miatto, Y. Yao, F. Sciarrino, and E. Karimi, Two-photon interference: the Hong–Ou–Mandel effect, Rep. Prog. Phys. 84 012402 (2021).
  * (8) I. Aharonovich, D. Englund, and M. Toth, Solid-state single-photon emitters, Nat. Photon. 10 631 (2016).
  * (9) C. Santori, D. Fattal, J. Vučković, G. S. Solomon, and Y. Yamamoto, Indistinguishable photons from a single-photon device, Nature 419, 594 (2002).
  * (10) N. Somaschi, V. Giesz, L. De Santis, J. C. Loredo, M. P. Almeida, G. Hornecker, S. L. Portalupi, T. Grange, C. Antón, J. Demory et al., Near-optimal single-photon sources in the solid state, Nat. Photonics 10, 340 (2016).
  * (11) A. Thoma, P. Schnauber, J. Böhm, M. Gschrey, J.-H. Schulze, A. Strittmatter, S. Rodt, T. Heindel, and S. Reitzenstein, Two-photon interference from remote deterministic quantum dot microlenses, Appl. Phys. Lett. 110, 011104 (2017).
  * (12) W.-B. Gao, P. Fallahi, E. Togan, A. Delteil, Y.S. Chin, J. Miguel-Sanchez, and A. Imamoğlu, Quantum teleportation from a propagating photon to a solid-state spin qubit, Nat. Commun. 4, 2744 (2013).
  * (13) A. Sipahigil, M. L. Goldman, E. Togan, Y. Chu, M. Markham, D. J. Twitchen, A. S. Zibrov, A. Kubanek, and M. D. Lukin, Quantum Interference of Single Photons from Remote Nitrogen-Vacancy Centers in Diamond, Phys. Rev. Lett. 108, 143601 (2012).
  * (14) A. Sipahigil, K. D. Jahnke, L. J. Rogers, T. Teraji, J. Isoya, A. S. Zibrov, F. Jelezko, and M. D. Lukin, Indistinguishable Photons from Separated Silicon-Vacancy Centers in Diamond, Phys. Rev. Lett. 113, 113602 (2014).
  * (15) N. Morioka, N. Morioka, C. Babin, R. Nagy, I. Gediz, E. Hesselmeier, D. Liu, M. Joliffe, M. Niethammer, D. Dasari, V. Vorobyov, et al., Spin-controlled generation of indistinguishable and distinguishable photons from silicon vacancy centres in silicon carbide, Nature Commun. 11, 2516 (2020).
  * (16) C. Chakraborty, L. Kinnischtzke, K. M. Goodfellow, R. Beams, and A. N. Vamivakas, Voltage-controlled quantum light from an atomically thin semiconductor, Nat. Nanotechnol. 10, 507 (2015).
  * (17) Y.-M. He, G. Clark, J. R. Schaibley, Y. He, M.-C. Chen, Y.-J. Wei, X. Ding, Q. Zhang, W. Yao, X. Xu, Single quantum emitters in monolayer semiconductors, et al, Nat. Nanotechnol. 10, 497 (2015).
  * (18) M. Koperski, K. Nogajewski, A. Arora, V. Cherkez, P. Mallet, J.-Y. Veuillen, J. Marcus, P. Kossacki, and M. Potemski, Single photon emitters in exfoliated WSe2 structures , Nat. Nanotechnol. 10, 503 (2015).
  * (19) A. Srivastava, M. Sidler, A. V. Allain, D. S. Lembke, A. Kis, and A. Imamoğlu, Optically active quantum dots in monolayer WSe2 Nat. Nanotechnol. 10, 491 (2015).
  * (20) P. Tonndorf, R. Schmidt, R. Schneider, J. Kern, M. Buscema, G. A. Steele, A. Castellanos-Gomez, H. S. J. van der Zant, S; Michaelis de Vasconcellos, and R. Bratschitsch, Single-photon emission from localized excitons in an atomically thin semiconductor, Optica 2, 347 (2015).
  * (21) T. T. Tran, K. Bray, M. J. Ford, M. Toth, and I. Aharonovich, Quantum emission from hexagonal boron nitride monolayers, Nat. Nanotechnol. 11, 37 (2016).
  * (22) R. Bourrelier, S. Meuret, A. Tararan, O. Stéphan, M. Kociak, L. H. G. Tizei, and A. Zobelli, Bright UV Single Photon Emission at Point Defects in h-BN, Nano Lett. 16, 4317 (2016).
  * (23) L. J. Martínez, T. Pelini, V. Waselowski, J. R. Maze, B. Gil, G. Cassabois, and V. Jacques, Efficient single photon emission from a high-purity hexagonal boron nitride crystal, Phys. Rev. B 94, 121405 (2016).
  * (24) A. K. Geim, and I. V. Grigorieva, Van der Waals heterostructures, Nature 499, 419 (2013).
  * (25) C. Fournier, A. Plaud, S. Roux, A. Pierret, M. Rosticher, K. Watanabe, T. Taniguchi, S. Buil, X. Quélin, J. Barjon, J.-P. Hermier, A. Delteil, Position-controlled quantum emitters with reproducible emission wavelength in hexagonal boron nitride, Nat. Commun. 12, 3779 (2021).
  * (26) A. Gale, C. Li, Y. Chen, K. Watanabe, T. Taniguchi, I. Aharonovich, and M. Toth, Site-Specific Fabrication of Blue Quantum Emitters in Hexagonal Boron Nitride, ACS Photonics 9, 2170 (2022).
  * (27) B. Shevitski, M. Gilbert, C. T. Chen, C. Kastl, E. S. Barnard, E. Wong, D. F. Ogletree, K. Watanabe, T. Taniguchi, A. Zettl, and S. Aloni, Phys. Rev. B 100, 155419 (2019).
  * (28) J. Horder, S. White, A. Gale, C. Li, K. Watanabe, T. Taniguchi, M. Kianinia, I. Aharonovich, and M. Toth, Coherence Properties of Electron-Beam-Activated Emitters in Hexagonal Boron Nitride Under Resonant Excitation, Phys. Rev. Appl. 18, 064021 (2022).
  * (29) T. Taniguchi, and K. Watanabe, Synthesis of High-Purity Boron Nitride Single Crystals under High Pressure by using Ba–BN Solvent, J. Cryst. Growth 303, 525 (2007).
  * (30) J. Bylander, I. Robert-Philip, and I. Abram, Interference and correlation of two independent photons, Eur. Phys. J. D 22, 295 (2003).
  * (31) A. Kiraz, M. Atatüre, and A. Imamoğlu, Quantum-dot single-photon sources: Prospects for applications in linear optics quantum-information processing, Phys. Rev. A 69, 032305 (2004).
  * (32) J. A. Martínez, R. A. Parker, K. C. Chen, . M. Purser, L. Li, C. P. Michaels, A. M. Stramma, R. Debroux, I. B. Harris, M. H. Appel, et al., Photonic Indistinguishability of the Tin-Vacancy Center in Nanostructured Diamond, Phys. Rev. Lett. 129, 173603 (2022).
  * (33) R. Proux, M. Maragkou, E. Baudin, C. Voisin, P. Roussignol, and C. Diederichs, Measuring the Photon Coalescence Time Window in the Continuous-Wave Regime for Resonantly Driven Semiconductor Quantum Dots, Phys. Rev. Lett. 114, 067401 (2015).
  * (34) J.-H. Kim, T. Cai, C. J. K. Richardson, R. P. Leavitt, and E. Waks, Two-photon interference from a bright single-photon source at telecom wavelengths, Optica 3, 577 (2016).
  * (35) C. Nawrath, F. Olbrich, M. Paul, S. L. Portalupi, M. Jetter, and P. Michler, Coherence and indistinguishability of highly pure single photons from non-resonantly and resonantly excited telecom C-band quantum dots, Appl. Phys. Lett. 115, 023103 (2019).
  * (36) R. Loudon, The Quantum Theory of Light, 2nd ed. (Oxford University Press, New York, 1983).
  * (37) K. A. Fischer, K. Müller, K. G. Lagoudakis, and J. Vučković, Dynamical modeling of pulsed two-photon interference, New J. Phys. 18, 113053 (2016).
  * (38) J. R. Johansson, P. D. Nation, and F. Nori, QuTiP: An opensource Python framework for the dynamics of open quantum systems, Comput. Phys. Commun. 183, 1760 (2012).
  * (39) J. R. Johansson, P. D. Nation, and F. Nori, QuTiP 2: A Python framework for the dynamics of open quantum systems, Comput. Phys. Commun. 184, 1234 (2013).
  * (40) See Supplemental Material at [URL will be inserted by publisher] for the detail of the simulations of time-gated two-photon interference.
  * (41) C. Fournier, K. Watanabe, T. Taniguchi, S. Buil, J. Barjon, J.-P. Hermier, and A. Delteil, Investigating the fast spectral diffusion of a quantum emitter in hBN using resonant excitation and photon correlations, arXiv:2303.05315 (2023).
  * (42) J. E. Fröch, C. Li, Y. Chen, M. Toth, M. Kianinia, S. Kim, and I. Aharonovich, Purcell Enhancement of a Cavity-Coupled Emitter in Hexagonal Boron Nitride, Small 18, 2104805 (2021).
  * (43) O. Iff, Q. Buchinger, M. Moczała-Dusanowska, M. Kamp, S. Betzold, M. Davanco, K. Srinivasan, S. Tongay, C. Antń-Solanas, S. Höfling, and C. Schneider, Purcell-Enhanced Single Photon Source Based on a Deterministically Placed WSe2 Monolayer Quantum Dot in a Circular Bragg Grating Cavity, Nano Lett. 21, 4715 (2021).
  * (44) I. Zhigulin, J. Horder, V. Ivady, S. J. U. White, A. Gale, C. Li, C. J. Lobo, M. Toth, I. Aharonovich, and M. Kianinia, Stark effect of quantum blue emitters in hBN, arXiv:2208.00600 (2022).

Supplementary Material

Two-photon interference from a quantum emitter in hexagonal boron nitride
Clarisse Fournier1, Sébastien Roux1,2, Kenji Watanabe3, Takashi Taniguchi4,
Julien Barjon1, Stéphanie Buil1, Jean-Pierre Hermier1, Aymeric Delteil1

1 Université Paris-Saclay, UVSQ, CNRS, GEMaC, 78000, Versailles, France.
2 Université Paris-Saclay, ONERA, CNRS, Laboratoire d’étude des
microstructures, 92322, Châtillon, France.
3 Research Center for Functional Materials, National Institute for Materials
Science, 1-1 Namiki, Tsukuba 305-0044, Japan
4 International Center for Materials Nanoarchitectonics, National Institute
for Materials Science, 1-1 Namiki, Tsukuba 305-0044, Japan
——————–<EMAIL_ADDRESS>

Influence of the post-selection on the HOM visibility

It is well known that dephasing degrades the visibility of two-photon
interference when the time difference between coincidences is larger than the
dephasing time [33-35]. Time postselection of photon detection events, allowed
by suitable detector resolution, can thus increase the Hong-Ou-Mandel (HOM)
visibility and has been used to boost the fidelity of quantum protocols based
on HOM interference – at the price of a lower count rate [3,4]. In this
section, we calculate the evolution of the HOM visibility upon modification of
the post-selection time window based on numerical integration of the master
equation.

Our system is modeled by a two-level system with computational states
$|g\rangle$ and $|e\rangle$, and is described in the rotating frame, in the
absence of coherent laser drive. As a consequence, the time evolution is
purely dissipative, governed by the Lindblad master equation
$\dot{\rho}=\sum_{c_{i}}c_{i}\rho
c_{i}^{\dagger}-\frac{1}{2}\left\\{c_{i}^{\dagger}c_{i},\rho\right\\}$, with
$c_{i}$ the relevant collapse operators. The dissipation terms are spontaneous
emission at rate $\Gamma_{\mathrm{sp}}$ and pure dephasing at rate
$\gamma^{*}$. The corresponding collapse operators read, respectively,
$\sqrt{\Gamma_{\mathrm{sp}}}\sigma^{-}$ and $\sqrt{\gamma^{*}/2}\sigma_{z}$
(in the convention where $\gamma^{*}=1/T_{2}^{*}$ [32]), where
$\sigma_{z}=|e\rangle\langle e|-|g\rangle\langle g|$, and
$\sigma^{-}=|g\rangle\langle e|$. The effect of the non-resonant laser pulse
is accounted for by taking $|e\rangle$ as the initial state. In the far field,
the field operators are locked to the source operators [36], such that
$a(t)\propto\sigma^{-}\left(t-r/c\right)$, with $a(t)$ the photon annihilation
operator. The field observables are therefore expressed in terms of source
operators in the following.

The HOM visibility can be expressed as a function of the first- and second-
order coherence of the light field [37]:

$g^{(2)}_{\mathrm{HOM}}(0)=\frac{1}{2}g^{(2)}_{\mathrm{HBT}}(0)+\frac{1}{2}\left(1-\frac{1}{N^{2}}\int\int\mathrm{d}t\mathrm{d}t^{\prime}\left|G^{(1)}(t,t^{\prime})\right|^{2}\right)$

with
$G^{(1)}(t,t^{\prime})=\langle\sigma^{+}(t)\sigma^{-}(t^{\prime})\rangle$, $N$
a normalization factor
$N=\int\mathrm{d}t\langle\sigma^{+}(t)\sigma^{-}(t)\rangle=\int\mathrm{d}t\langle\sigma_{ee}(t)\rangle$
and $g^{(2)}_{\mathrm{HBT}}(0)=0$ for an ideal single-photon source.

$G^{(1)}(t,t^{\prime})$ is numerically calculated based on the quantum
regression theorem using the Qutip toolbox [38,39]. To account for the
temporal post-selection, the integrals run over $t,t^{\prime}\in[0,\Delta t]$,
where $\Delta t$ is the post-selection window size.
$g^{(2)}_{\mathrm{HOM}}(0)$ is calculated as a function of $\Delta t$ for
various values of $\gamma^{*}$.

Figure S1 shows the result of the calculation for six values of $\gamma^{*}$.
It can be seen that the visibility
$V_{\mathrm{HOM}}=1-2g^{(2)}_{\mathrm{HOM}}(0)$ decays with increasing $\Delta
t$. The decay time depends on the dephasing, and is plotted as a function of
$\gamma^{*}$ in main text figure 3b. The visibility for $\Delta t\approx 0$ is
always high since the photons are always indistinguishable at short delays
[33], and decays faster for larger dephasing rates. The long-time asymptote,
corresponding to integration of all detection events, coincides with the bare
(non-postselected) indistinguishability
$T_{2}/2T_{1}=\Gamma_{\mathrm{sp}}/(\Gamma_{\mathrm{sp}}+2\gamma^{*})$ [30].

Figure S1: $V_{\mathrm{HOM}}$ as a function of $\Delta t$ for different values
of $\gamma^{*}$.

Influence of the Purcell factor on the HOM visibility

As mentioned in the previous section, the integrated HOM visibility (without
any post-selection) for an ideal emitter is $V_{\mathrm{HOM}}=T_{2}/2T_{1}$
[30] – the ideal case $T_{2}=2T_{1}$ providing $V_{\mathrm{HOM}}=1$. In a
cavity, the Purcell effect modifies the lifetime according to
$T_{1}^{\prime}=T_{1}/F_{p}$, where $F_{p}$ is the Purcell factor. Given our
estimation of $T_{2}^{*}=2.4$ ns, an indistinguishability of 0.80 (resp 0.90)
corresponds to $F_{p}=7$ (resp. $F_{p}=15$).

## References

  * (1) [3] H. Bernien, B. Hensen, W. Pfaff, G. Koolstra, M. S. Blok, L. Robledo, T. H. Taminiau, M. Markham, D. J. Twitchen, L. Childress, and R. Hanson, Heralded entanglement between solid-state qubits separated by three metres, Nature 497, 86 (2013).
  * (2) [4] A. Delteil, Z. Sun, W.-B. Gao, E. Togan, S. Fält, and A. Imamoğlu, Generation of heralded entanglement between distant hole spins, Nature Phys. 12, 218 (2016).
  * (3) [30] J. Bylander, I. Robert-Philip, and I. Abram, Interference and correlation of two independent photons, Eur. Phys. J. D 22, 295 (2003).
  * (4) [32] J. A. Martínez, R. A. Parker, K. C. Chen, . M. Purser, L. Li, C. P. Michaels, A. M. Stramma, R. Debroux, I. B. Harris, M. H. Appel, et al., Photonic Indistinguishability of the Tin-Vacancy Center in Nanostructured Diamond, Phys. Rev. Lett. 129, 173603 (2022).
  * (5) [33] R. Proux, M. Maragkou, E. Baudin, C. Voisin, P. Roussignol, and C. Diederichs, Measuring the Photon Coalescence Time Window in the Continuous-Wave Regime for Resonantly Driven Semiconductor Quantum Dots, Phys. Rev. Lett. 114, 067401 (2015).
  * (6) [34] J.-H. Kim, T. Cai, C. J. K. Richardson, R. P. Leavitt, and E. Waks, Two-photon interference from a bright single-photon source at telecom wavelengths, Optica 3, 577 (2016).
  * (7) [35] C. Nawrath, F. Olbrich, M. Paul, S. L. Portalupi, M. Jetter, and P. Michler, Coherence and indistinguishability of highly pure single photons from non-resonantly and resonantly excited telecom C-band quantum dots, Appl. Phys. Lett. 115, 023103 (2019).
  * (8) [36] R. Loudon, The Quantum Theory of Light, 2nd ed. (Oxford University Press, New York, 1983).
  * (9) [37] K. A. Fischer, K. Müller, K. G. Lagoudakis, and J. Vučković, Dynamical modeling of pulsed two-photon interference, New J. Phys. 18, 113053 (2016).
  * (10) [38] J. R. Johansson, P. D. Nation, and F. Nori, QuTiP: An opensource Python framework for the dynamics of open quantum systems, Comput. Phys. Commun. 183, 1760 (2012).
  * (11) [39] J. R. Johansson, P. D. Nation, and F. Nori, QuTiP 2: A Python framework for the dynamics of open quantum systems, Comput. Phys. Commun. 184, 1234 (2013).