# Foveal-pit inspired filtering of DVS spike response
††thanks:

1st Shriya T.P. Gupta  Microsoft Research and Development Center
Microsoft Corporation Pvt. Ltd.
Bangalore, India
<EMAIL_ADDRESS>2nd Pablo Linares-Serrano  3rd Basabdatta Sen
Bhattacharya  Instituto de Microelectrónica de Sevilla
IMSE-CNM (CSIC and Universidad de Sevilla)
Sevilla, Spain
<EMAIL_ADDRESS>Department of Computer Science
BITS Pilani Goa Campus
Goa, India
<EMAIL_ADDRESS>4th Teresa Serrano-Gotarredona  Instituto
de Microelectrónica de Sevilla
IMSE-CNM (CSIC and Universidad de Sevilla)
Sevilla, Spain
<EMAIL_ADDRESS>
###### Abstract

In this paper, we present results of processing Dynamic Vision Sensor (DVS)
recordings of visual patterns with a retinal model based on foveal-pit
inspired Difference of Gaussian (DoG) filters. A DVS sensor was stimulated
with varying number of vertical white and black bars of different spatial
frequencies moving horizontally at a constant velocity. The output spikes
generated by the DVS sensor were applied as input to a set of DoG filters
inspired by the receptive field structure of the primate visual pathway. In
particular, these filters mimic the receptive fields of the midget and parasol
ganglion cells (spiking neurons of the retina) that sub-serve the photo-
receptors of the foveal-pit. The features extracted with the foveal-pit model
are used for further classification using a spiking convolutional neural
network trained with a backpropagation variant adapted for spiking neural
networks.

###### Index Terms:

dynamic vision sensor, neural filtering, spiking neural network,
classification, difference of gaussian, convolution, foveal-pit

##  I Introduction

Recent advances in deep learning [1, 2] have led to state-of-the-art
performance for varied classification tasks in natural language processing,
computer vision and speech recognition. Traditional Artificial Neural Networks
(ANN) use idealized computing units which have a differentiable, non-linear
activation function allowing stacking of such neurons in multiple trainable
layers. The existence of derivatives makes it possible to carry out large
scale training of these architectures with gradient based optimization methods
[3] using high computing resources like Graphic Processing Units (GPU).
However, this prevents the use of such deep learning models for essential
real-life applications like mobile devices and autonomous systems that have
limited compute power.

Spiking Neural Networks (SNN) have been proposed as an energy-efficient
alternative to ANNs as they simulate the event-based information processing of
the brain [4]. These bio-inspired SNNs follow an asynchronous method of event
processing using spiking neurons. The internal state of a spiking neuron is
updated when it receives an action potential and consequently an output spike
is fired when the membrane voltage crosses a pre-defined threshold. Further,
improvements in neuromorphic engineering allow the implementation of SNNs on
neuromorphic hardware platforms [5] that lead to a much higher efficiency in
terms of power and speed compared to conventional GPU based computing systems.

Although SNNs are considered as the third generation of neural networks
holding the potential for sparse and low-power computation, their
classification performance is considerably lower than those of ANNs. This can
be attributed to the fact that gradient optimization techniques like the
backpropagation algorithm can’t be implemented in SNNs due to the discrete
nature of spiking neurons. A common technique for training SNN models is the
Hebbian learning inspired Spike Timing Dependent Plasticity (STDP) that is
used in several state-of-the-art approaches [6, 7]. Other works like [8, 9]
have adapted the gradient descent algorithm for SNNs using a differentiable
approximation of spiking neurons. Our approach also employs a similar modified
backpropagation algorithm proposed by Hunsberger et al. [10] that is
implemented in the Nengo-DL library [11].

As shown by Camunas-Mesa et al. [12], the efficiency gain of SNNs from event-
based processing can be further improved through the use of inputs from event-
based sensors like a neuromorphic Dynamic Vision Sensor (DVS) [13]. Event
driven sensors represent the information dynamically by asynchronously
transmitting the address event of each pixel and hence avoid processing
redundant data. However, the classification accuracy drops drastically when
using real sensory data from a physical spiking silicon retina, since the
spike events are no longer Poissonian [14].

Figure 1: The two-stage architecture of the proposed DVS based spiking
convolutional neural network. In the first stage, the DVS is stimulated with a
pattern of vertical white and black bars and the generated spike responses are
processed using foveal-pit inspired DoG filters. The second stage is composed
of the convolutional and pooling layers which are used to classify the
features extracted from the first stage.

In our previous work [15], we had demonstrated the effect of foveal-pit
inspired filtering for synthetically generated datasets like MNIST [3] and
Caltech [16]. In this work, we present the results of applying similar neural
filtering to data generated by the DVS. In our proposed model, we process DVS
outputs using bio-inspired filters that simulate receptive fields of the
midget and parasol ganglion cells of the primate retina. The DVS is stimulated
with vertical black and white bars having a constant displacement of 2 pixels
from frame to frame. The foveal-pit informed Difference of Gaussian (DoG)
filters are applied to the DVS recordings in order to capture the most
perceptually important information from the input data. The use of DoG
functions to model retinal filters was originally proposed by Rullen et al.
[17] and the receptive fields of the foveal-pit are implemented as in
Bhattacharya et al. [18].

The processed features are then used to perform the classification using a
Spiking Convolutional Neural Network (SCNN). The SCNN architecture is inspired
by two previous works viz. Diehl et al. [7] and Kheradpisheh et al. [6], while
the model is implemented as in Gupta et al. [15]. Each input is presented to
the network for a total duration of 60 timesteps and the predictions are
assigned based on the voltages measured from the output neurons. The empirical
results demonstrate that the application of neural filtering to DVS recordings
leads to an improvement of 35% in classification accuracy compared to the
unfiltered DVS spike responses. Out of the filtered scenarios, the highest
performance of 100% is achieved using the off-center parasol ganglion cells.

The rest of the paper is organized as follows: Section II describes the
architecture of the proposed model including the response generation and
filtering, Section III provides the results of the experiments and Section IV
contains the conclusion and future directions.

##  II Methodology

The overall architecture of our model consists of two main stages: the first
stage is made up of the DVS response generation and neural filtering of output
spikes; the second stage consists of performing classification using the SCNN.
The proposed model is shown in Fig. 1 and each of the individual stages are
covered in detail in the following sub-sections.

###  II-A Dynamic Vision Sensor Responses

We have used a $128\times 128$ sized neuromorphic DVS developed by Serrano-
Gotarredona et al. [13] to capture non-synthetic visual data. Each pixel of
the DVS processes the input continuously and emits a spike based on the
variation in the illumination impinging upon it [19]. A sample illustration is
provided in Fig. 2 using a sinusoidal input stimulus having a frequency of
10Hz. The first row depicts the pixel’s illumination over time whereas the
remaining two rows capture the emission of spikes over the same duration
corresponding to changes in illumination. An increase in illumination leads to
a positive spike whereas a decrease in illumination leads to a negative spike
as seen in the last row of Fig. 2.

Figure 2: Illustration of the DVS spike response generated using a sinusoidal
input stimulus. The first row represents the pixel’s illumination over time.
The second row depicts the positive spikes whereas the last row represents the
negative spikes corresponding to a decrease in illumination [19].

For our experiments, the DVS was placed in front of a monitor displaying a
pattern of equally wide black and white vertical bars as shown in Fig. 3. The
bars were moved horizontally across the screen such that a displacement of 2
pixels is applied from frame to frame. The number of bars were varied from 2,
4, 8, 16, 32, 64 to 128 and these $K=7$ categories correspond to the final
labels for our multiclass classification problem. The events generated by the
DVS were captured in the Address Event Representation (AER) format using the
jAER software [20].

Figure 3: The DVS setup to record spike responses when presented with a simple
pattern moving across a computer screen as visual stimulus.

###  II-B Retina-inspired filtering

The DVS recordings generated from the first stage of our model are passed to a
set of neural filters simulating the primate visual system. As proposed by
Kheradpisheh et al. [6], we have used DoG functions to implement these
biologically inspired filters sub-serving the retinal foveal pit. The foveal
pit is a circular region of 200 $\mu m$ diameter that lies at the center of
the foveola. This region also has the highest visual acuity in the primate
retina and is most accessible to incoming light. The fovea is sub-served by a
ganglion cell layer composed of midget and parasol cells. The retinal ganglion
cells are the only spiking neurons of the primate visual system and their
axons transmit the received information from the retina to other parts of the
brain.

1: kernel = dog_func(mat_dim, cent_dev, circ_shift)

2: [ma, na] = size(input)

3: [mb, nb] = size(kernel)

4: filt_out = zeros(ma, na)

5: r1 = ceil(mb/2)

6: s1 = ceil(nb/2)

7: for i = 1 to i = ma do

8: for j =1 to j = na do

9: i1 = max(0, i-r1);

10: for r = max(1, r1-i+1) to r = mb do

11: i1 = i1 + 1

12: j1 = max(0, j-s1)

13: for s = max(1, s1-j+1) to s = nb do

14: j1 = j1 + 1

15: filt_out(i, j) += kernel(r, s) * input(i1, j1)

16: end for

17: end for

18: end for

19: end for

20: return filt_out

Algorithm 1 Algorithm for filtering the DVS spike response using foveal-pit
inspired DoG functions.

The midget and parasol ganglion cells have two types of centre surround
receptive fields — on-centre-off-surround and off-centre-on-surround. We have
modelled these receptive fields using DoG functions as specified in
Bhattacharya et al. [18]. The off-center midget cells have a matrix size of
$5\times 5$ with standard deviation of 0.8 whereas the on-center midget cells
are of size $11\times 11$ with standard deviation of 1.04. Similarly, the off-
center parasol cells have a size of $61\times 61$ with a standard deviation of
8 while the on-center parasol cells are of size $243\times 243$ with a
standard deviation of 10.4. These DoG functions are then applied to the DVS
spike responses using Algorithm 1.

###  II-C Convolutional Network Architecture

The asynchronous DVS recordings generated from the previous stage are split
into individual frames for training our frame-based classifier composed of
convolutional layers. This modified dataset is created following the procedure
of Stromatias et al. [21] to produce an analog vector representation. The SCNN
architecture used in our work consists of three convolutional and pooling
layers which are made up of Leaky Integrate and Fire (LIF) neurons.

TABLE I: Dimensions of the SCNN layers. Layer | No. of filters | Input size | Kernel size
---|---|---|---
Conv1 | 8 | (128, 128) | 3
Pool1 | - | (128, 128) | 2
Conv2 | 16 | (64, 64) | 3
Pool2 | - | (64, 64) | 2
Conv3 | 32 | (32, 32) | 3
Pool3 | - | (32, 32) | 2
Flatten | - | (16, 16) | -
Dense | - | (1, 8192) | -
Outputs | - | (1, 7) | -

Traditional deep learning architectures use sigmoid neurons which are
differentiable non-linearities, whereas the spiking neurons used in SCNNs are
non-differentiable. Hence, we use a differentiable approximation of the
spiking neurons during training and the actual spiking neurons during
inference as proposed by Hunsberger et al. [10]. Since we use a rate-based
approximation during training, the model is run only for a single timestep
whereas during testing with the spiking neurons, the model is run for 60
timesteps to collect the cumulative spike output over time.

The convolution is carried out on the $128\times 128$ input arrays using
filters of size $3\times 3$. The first, second and third convolutional layers
of the SCNN are made up of $2^{3}$, $2^{4}$ and $2^{5}$ filters respectively,
followed by a pooling operation after each convolution. The synaptic
connections between the neurons of these layers are modelled as the trainable
weights of the network which are optimized by minimizing the loss function of
the overall SCNN. The exact dimensions of the individual layers are provided
in Table I.

###  II-D Training and Inference

For our multiclass classification problem with $K=7$ categories, we convert
the outputs of the last pooling layer into a 1-D vector using a flatten
operation. This is followed by a dense layer with all-to-all connectivity
having K neurons which generates a $K\times 1$ output vector. A softmax
classifier is used to transform these output values into a set of K
probabilities:

$\mathbf{Y}(X,W)=\frac{e^{w_{i}x_{i}}}{\sum_{j=1}^{K}e^{w_{j}x_{j}}}\forall
i=1,\cdots,K$ (1)

where $x_{i}\in X$ and $w_{i}\in W$ are the inputs and weights of the dense
layer respectively, and $\mathbf{Y}$ is the prediction probabilities that sum
to 1. The Negative Log Likelihood (NLL) loss for the overall network is
computed using the one-hot encoded output labels $\mathbf{L}$ and the
softmaxed probabilities $\mathbf{Y}$ with NLL defined as:

$\mathbf{O}(X,W)=-\frac{1}{M}\sum_{i}^{M}\sum_{j}^{K}F_{i}(j)*\log(\mathbf{Y}(X,W))$
(2)

where M is the mini-batch size, $F_{i}(j)=1$ when $j=L_{i}$ and zero
otherwise. The SCNN is trained end-to-end using a spiking approximation of the
backpropogation algorithm adapted for SNNs. This is done by minimizing the NLL
loss using the procedure described in Gupta et al. [15] with a duration of 3
epochs and a mini-batch size of 20.

For the inference stage, we pass the input images from the testing corpora and
measure the voltages (mV) of the output layer neurons. These values are
generated using the probe function of the Nengo-DL library and represent the
progressively increasing membrane potentials. Thus, the neuron having the
highest voltage over a 60 ms simulation time period is assigned as the
predicted class for that epoch.

##  III Experimental Methods and Results

The filtered spikes responses from the 128$\times$128 sized DVS sensor was
split into individual frames for each recording to be passed as input to the
subsequent convolutional network. This resulted in a total of 3552 images for
the unfiltered scenario and a collection of 3503 images for the filtered
recordings. In each case, the images were then partitioned in the ratio of 9:1
to create the corresponding training and testing corpora.

Generation of the DVS spikes responses along with the filtering was
implemented entirely in Matlab, while the SCNN and its various layers were
coded in Python using the Nengo-DL library [11]. The generated .mat files of
the dataset were loaded into the Python network using the Scipy library [22]
and the experiments were carried out on the GPU accessed via Google
Colaboratory [23].

###  III-A Quantitative Effects of Filtering

To assess the effects of incorporating the neural filtering on DVS recordings,
we ran two experiments with the SCNN for a total duration of 60 timesteps. The
empirical results are summarized in Table II. For the first scenario of using
unfiltered DVS frames, the model achieves an accuracy of 65% which is
significantly lower than the values in the remaining rows that correspond to
the filtered DVS inputs. This demonstrates that introducing the foveal-pit
inspired neural filtering into our retinal model leads to a considerable
improvement of 35% even for simplistic visual patterns such as those in our
dataset.

TABLE II: Accuracies (%) for the frame-based DVS input Scenario | Cell - Type | CircShift | Accuracy
---|---|---|---
Unfiltered | - | - | 65.0 %
| off-center midget |  | 77.5 %
Filtered | on-center midget | 0 | 85.0 %
| off-center parasol |  | 92.5 %
| on-center parasol |  | 87.5 %
| off-center midget |  | 77.5 %
Filtered | on-center midget | 1 | 85.0 %
| off-center parasol |  | 100.0 %
| on-center parasol |  | 85.0 %

Amongst the filtered outputs, the parasol ganglion cells have a comparatively
higher increase in accuracy compared to the midget cells. Since the parasol
cells have larger dimensions and capture the overall background information,
they lead to a significant improvement in classification of distinct patterns
without intricate details. Thus, the parasol cells lead to a larger
performance gain achieving a highest accuracy in the shifted case, as our
dataset is composed of only vertical black and white bars. On the other hand,
midget cells have smaller dimensions which allows them to capture only the
finer details of an image and hence they contribute lesser to the overall
increase in classification accuracy.

From Table II, we can also observe that the variations in accuracy for
filtering with different ganglion cell types is almost comparable for both the
cases of with and without any circular-shift; circular-shift refers to the
case where the DoG filters are circular shifted to ‘wrap’ on the raster at all
four edges. This reduces artifacts due to edges. The alternative scenario is
zero padding at all four edges. The circular-shift value in the Table II is
set to 1 indicating all cases where filtering was performed using circular
shift at edges, and is set to 0 otherwise.

###  III-B Qualitative Effects of Filtering

For analysing the qualitative effects of neural filtering, we generated raster
plots using the analog vector representation of DVS responses as shown in Fig.
4. The neuron numbers range from a value of 0 to 16384 as they represent the
pixels of the 128$\times$128 sized electronic retina used in our experiments.
The blue markers depict a positive event corresponding to an increase in
illumination as the moving edges go from black to white. On the contrary, the
red markers represent a negative event and indicate a decrease in illumination
as the edges go from white to black.

Figure 4(a) illustrates the unfiltered scenario which has the least
distinction between positive and negative events of the input stimulus. This
lack of differentiation between black and white bars of the visual pattern
also leads to a drop in the classification accuracy which was previously
observed in Table II. Further, in the filtered raster plots of Fig. 4, we note
that all the ganglion cell filters capture edges more effectively compared to
the unfiltered case as there is a clear distinction in the positive and
negative events of the filtered raster plots. This improved distinction also
leads to a higher classification performance as seen in Section III-A.

(a) Unfiltered scenario

(b) Filtered: Unshifted midget

(c) Filtered: Unshifted parasol

(d) Filtered: Circular-shifted midget

(e) Filtered: Circular-shifted parasol

Figure 4: Raster plot of the (a) unfiltered DVS spike response and filtered
DVS response using the off-center (b) unshifted midget cell (c) unshifted
parasol cell (d) circular-shifted midget cell (e) circular-shifted parasol
cell.

Amongst the outputs of the filtered DVS responses, Fig. 4(b) and Fig. 4(c)
represent the cases having a circular-shift value of 0 while Fig. 4(d) and
Fig. 4(e) correspond to a circular-shift value of 1. The variation in filtered
responses are almost similar for these two cases as seen in Table II. However,
the result obtained using the unshifted parasol cell in Fig. 4(c) has lesser
clarity between the positive and negative events and hence contains less
distinguishable edges than in Fig. 4(e) using the circular-shifted filter.
Thus, the best classification accuracy of 100% in Section III-A was achieved
using the circular-shifted parasol filter as it captures the changes in
illumination more effectively.

Additionally, we observe that the plots generated using midget ganglion cells
leads to a noisier qualitative output than those using the parasol ganglion
cells. This is because the midget cells have smaller dimensions and are able
to pick up only the finer details of an image. But since the simplistic visual
patterns used in our experiments lack any intricate details, the midget cells
contribute lesser to the overall improvement in the classification
performance. Thus, the highest accuracies in Table II are obtained using the
parasol cells as they capture the larger and more significant information
contained in the input data.

##  IV Conclusion and Future Work

In this paper, we have presented a novel method for processing the DVS spike
responses of a visual pattern with foveal-pit inspired DoG filters that
simulate the primate retinal system. The pattern was composed of varying
number of vertical white and black bars of different spatial frequencies
moving at a fixed velocity. The outputs from the sensor are applied as input
to the bio-inspired neural filters that model the receptive field structure of
midget and parasol ganglion cells of the foveal-pit. These processed features
are passed as input to our spiking convolutional neural network architecture
which classifies the frame-based version of the filtered responses into seven
corresponding categories. The SCNN is composed of convolutional and pooling
layers and is trained with a modified backpropogation algorithm using a
differentiable approximation of spiking neurons [10].

The proposed model demonstrates the effect of applying neural filtering to
real DVS data generated from a neuromorphic vision sensor. This builds upon
our previous work [15] that depicted the results of foveal-pit inspired
filtering for synthetically generated datasets like MNIST [3] and Caltech
[16]. Our model achieves a promising performance of 92.5% using the unshifted
off-center parasol ganglion cell and an accuracy of 100% in the circular-
shifted scenario, which is an improvement of 35% over the classification using
unfiltered DVS responses. The empirical results indicate the importance of the
foveal-pit inspired neural filtering in redundancy reduction of the DVS inputs
and in discarding irrelevant background information.

For our proposed network, the asynchronous DVS recordings generated from the
first stage of the model were converted to an analog vector representation for
training the frame-based classifier composed of convolution layers. As future
work, we plan to adapt our spiking convolutional network architecture to
directly process event-based data and evaluate the effects of the bio-inspired
neural filtering on continuous outputs of a neuromorphic DVS. Also, the
dataset used in this work is limited in terms of variation in the inputs as
well as the size of the training and testing corpora. Hence, we would like to
further verify the effect of the DoG filters on DVS spike responses of larger
and more complex datasets.

##  V Acknowledgements

This work was funded by the EU grant PCI2019-111826-2 “APROVIS3D”, by Spanish
grant from the Ministry of Science and Innovation PID2019-105556GB-C31
“NANOMIND” (with support from the European Regional Development Fund) and by
the CSIC 2018-50E008 AVE project. BSB is supported by the Science and
Engineering Research Board (SERB) of India Fund CRG/2019/003534.

## References

  * [1] Y. Bengio, Y. LeCun _et al._ , “Scaling learning algorithms towards ai,” _Large-scale kernel machines_ , vol. 34, no. 5, pp. 1–41, 2007.
  * [2] Y. LeCun, Y. Bengio, and G. Hinton, “Deep learning,” _nature_ , vol. 521, no. 7553, p. 436, 2015.
  * [3] Y. LeCun, L. Bottou, Y. Bengio, P. Haffner _et al._ , “Gradient-based learning applied to document recognition,” _Proceedings of the IEEE_ , vol. 86, no. 11, pp. 2278–2324, 1998.
  * [4] W. Gerstner and W. M. Kistler, _Spiking neuron models: Single neurons, populations, plasticity_. Cambridge university press, 2002.
  * [5] A. R. Young, M. E. Dean, J. S. Plank, and G. S. Rose, “A review of spiking neuromorphic hardware communication systems,” _IEEE Access_ , vol. 7, pp. 135 606–135 620, 2019.
  * [6] S. R. Kheradpisheh, M. Ganjtabesh, S. J. Thorpe, and T. Masquelier, “Stdp-based spiking deep convolutional neural networks for object recognition,” _Neural Networks_ , vol. 99, pp. 56–67, 2018.
  * [7] P. U. Diehl, D. Neil, J. Binas, M. Cook, S.-C. Liu, and M. Pfeiffer, “Fast-classifying, high-accuracy spiking deep networks through weight and threshold balancing,” in _2015 International joint conference on neural networks (IJCNN)_. ieee, 2015, pp. 1–8.
  * [8] E. O. Neftci, C. Augustine, S. Paul, and G. Detorakis, “Event-driven random back-propagation: Enabling neuromorphic deep learning machines,” _Frontiers in neuroscience_ , vol. 11, p. 324, 2017.
  * [9] H. Mostafa, “Supervised learning based on temporal coding in spiking neural networks,” _IEEE transactions on neural networks and learning systems_ , vol. 29, no. 7, pp. 3227–3235, 2017.
  * [10] E. Hunsberger and C. Eliasmith, “Training spiking deep networks for neuromorphic hardware,” _arXiv preprint arXiv:1611.05141_ , 2016.
  * [11] T. Bekolay, J. Bergstra, E. Hunsberger, T. DeWolf, T. C. Stewart, D. Rasmussen, X. Choo, A. Voelker, and C. Eliasmith, “Nengo: a python tool for building large-scale functional brain models,” _Frontiers in neuroinformatics_ , vol. 7, p. 48, 2014.
  * [12] L. Camunas-Mesa, C. Zamarreno-Ramos, A. Linares-Barranco, A. J. Acosta-Jimenez, T. Serrano-Gotarredona, and B. Linares-Barranco, “An event-driven multi-kernel convolution processor module for event-driven vision sensors,” _IEEE Journal of Solid-State Circuits_ , vol. 47, no. 2, pp. 504–517, 2011\.
  * [13] T. Serrano-Gotarredona and B. Linares-Barranco, “A 128$\times$128 1.5% contrast sensitivity 0.9% fpn 3 $\mu$s latency 4 mw asynchronous frame-free dynamic vision sensor using transimpedance preamplifiers,” _IEEE Journal of Solid-State Circuits_ , vol. 48, no. 3, pp. 827–838, 2013.
  * [14] G. Orchard, A. Jayawant, G. K. Cohen, and N. Thakor, “Converting static image datasets to spiking neuromorphic datasets using saccades,” _Frontiers in neuroscience_ , vol. 9, p. 437, 2015.
  * [15] S. T. Gupta and B. S. Bhattacharya, “Implementing a foveal-pit inspired filter in a spiking convolutional neural network: a preliminary study,” in _2020 International Joint Conference on Neural Networks (IJCNN)_. IEEE, 2020, pp. 1–8.
  * [16] L. Fei-Fei, R. Fergus, and P. Perona, “Learning generative visual models from few training examples: An incremental bayesian approach tested on 101 object categories,” in _2004 conference on computer vision and pattern recognition workshop_. IEEE, 2004, pp. 178–178.
  * [17] R. V. Rullen and S. J. Thorpe, “Rate coding versus temporal order coding: what the retinal ganglion cells tell the visual cortex,” _Neural computation_ , vol. 13, no. 6, pp. 1255–1283, 2001.
  * [18] B. S. Bhattacharya and S. B. Furber, “Biologically inspired means for rank-order encoding images: A quantitative analysis,” _IEEE transactions on neural networks_ , vol. 21, no. 7, pp. 1087–1099, 2010.
  * [19] B. S. Bhattacharya, T. Serrano-Gotarredona, L. Balassa, A. Bhattacharya, A. B. Stokes, A. Rowley, I. Sugiarto, and S. Furber, “A spiking neural network model of the lateral geniculate nucleus on the spinnaker machine,” _Frontiers in neuroscience_ , vol. 11, p. 454, 2017.
  * [20] “Jaer project,” 2007. [Online]. Available: http://jaerproject.net/
  * [21] E. Stromatias, M. Soto, T. Serrano-Gotarredona, and B. Linares-Barranco, “An event-driven classifier for spiking neural networks fed with synthetic or dynamic vision sensor data,” _Frontiers in neuroscience_ , vol. 11, p. 350, 2017.
  * [22] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python,” _Nature Methods_ , vol. 17, pp. 261–272, 2020.
  * [23] “Colaboratory.” [Online]. Available: https://colab.research.google.com/


abbr

# Non-crossing convex quantile regression

Sheng Dai  Timo Kuosmanen  Xun Zhou

(April 2022)

###### Abstract

Quantile crossing is a common phenomenon in shape constrained nonparametric
quantile regression. A recent study by Wang2014c has proposed to address this
problem by imposing non-crossing constraints to convex quantile regression.
However, the non-crossing constraints may violate an intrinsic quantile
property. This paper proposes a penalized convex quantile regression approach
that can circumvent quantile crossing while better maintaining the quantile
property. A Monte Carlo study demonstrates the superiority of the proposed
penalized approach in addressing the quantile crossing problem.
Keywords: Quantile function, Quantile crossing, Convex quantile regression,
Simultaneous estimation, Regularization
JEL Codes: C1, C6, C13

##  1 Introduction

Quantile estimation has been widely applied in various fields of economics and
econometrics (see, e.g., Wang2014c, *Wang2014c; Jradi2019, *Jradi2019;
Tsionas2020b, *Tsionas2020b; Kuosmanen2020b, *Kuosmanen2020b; Zhao2021,
*Zhao2021). However, when multiple quantiles are separately estimated to
obtain a family of conditional quantile functions, two or more quantile curves
may cross on the condition that the distribution functions and their
associated inverse functions are not monotone increasing (He1997, *He1997).
Such quantile crossing is a longstanding problem in quantile regression.

To our knowledge, there are three commonly seen approaches to avoid quantile
crossing: post-processing, stepwise estimation, and simultaneous estimation.
In the post-processing procedure, a non-crossing assumption is usually
enforced via a sorting or monotonic rearrangement of the original estimated
non-monotone functions (e.g., Dette2008, *Dette2008; Chernozhukov2010,
*Chernozhukov2010). This indirect approach is effective in estimating the
conditional quantile, but lacks the ability to quantify the effects of the
predictors (Bondell2010, *Bondell2010). In the stepwise procedure, it prevents
an estimated quantile function from crossing the previously estimated one by
adding an extra set of non-crossing constraints iteratively to the regression
model (e.g., Wu2009a, *Wu2009a); but this approach cannot offer simultaneous
estimates. In the simultaneous estimation, non-crossing constraints are
imposed to ensure that the estimated conditional quantile functions are
monotone nondecreasing, with all quantiles being estimated simultaneously
(e.g., Takeuchi2006, *Takeuchi2006; Bondell2010, *Bondell2010). More recently,
Wang2014c extend this simultaneous estimation technique to convex quantile
regression (sCQR). However, the non-crossing constraints may disturb the
quantile property (Takeuchi2006, *Takeuchi2006).

This paper develops a new non-crossing approach to nonparametric quantile
function estimation. Compared with sCQR (Wang2014c, *Wang2014c), our approach
based on the penalized convex quantile regression (pCQR) independently
estimates multiple non-crossing quantiles but can better satisfy the intrinsic
quantile property. Furthermore, the proposed pCQR approach can fit the true
quantile functions more reliably and robustly.

##  2 Penalized convex quantile regression

Consider a general nonparametric regression model with observations
$\\{{\boldsymbol{x}}_{i},y_{i}\\}_{i=1}^{n}$ satisfying

$\displaystyle y_{i}=f({\boldsymbol{x}}_{i})+\varepsilon_{i},\qquad\mbox{for
}i=1,\ldots,n,$ (1)

where $y_{i}\in{\mathbb{R}}$ and ${\boldsymbol{x}}_{i}\in{\mathbb{R}}^{d}$ are
output and inputs variables, and $\varepsilon_{i}$ is a error term with zero
mean. Accordingly, for a given quantile $\tau\in(0,1)$, the nonparametric
quantile function $Q_{y}(\tau\,|\,{\boldsymbol{x}})$ is defined as

$\displaystyle
Q_{y_{i}}(\tau\,|\,{\boldsymbol{x}}_{i})=f({\boldsymbol{x}}_{i})+F_{\varepsilon_{i}}^{-1}(\tau),$
(2)

where $F_{\varepsilon_{i}}$ is the distribution function of the error term
$\varepsilon_{i}$.

To estimate quantiles empirically, we resort to convex quantile regression
(CQR) that does not require any assumptions about the functional form of the
regression function $f$ or its smoothness, but imposes the shape constraints
such as monotonicity and concavity. Specifically, CQR estimates the quantile
function (2) by solving the following linear programming problem (Wang2014c,
*Wang2014c)

$\displaystyle\underset{\alpha,{\boldsymbol{\beta}},\varepsilon^{\text{+}},\varepsilon^{-}}{\mathop{\min}}$
$\displaystyle\,\tau\sum\limits_{i=1}^{n}{\varepsilon_{i}^{+}}+(1-\tau)\sum\limits_{i=1}^{n}{\varepsilon_{i}^{-}}$
(3) s.t. $\displaystyle
y_{i}=\mathbf{\alpha}_{i}+{\boldsymbol{\beta}}_{i}^{{}^{\prime}}{\boldsymbol{x}}_{i}+\varepsilon_{i}^{+}-\varepsilon_{i}^{-}$
$\displaystyle\forall i$
$\displaystyle\alpha_{i}+{\boldsymbol{\beta}}_{i}^{{}^{\prime}}{\boldsymbol{x}}_{i}\leq\alpha_{h}+{\boldsymbol{\beta}}_{h}^{{}^{\prime}}{\boldsymbol{x}}_{i}$
$\displaystyle\forall i,h$
$\displaystyle{\boldsymbol{\beta}}_{i}\geq{\boldsymbol{0}}$
$\displaystyle\forall i$ $\displaystyle\varepsilon_{i}^{+}\geq 0,\
\varepsilon_{i}^{-}\geq 0$ $\displaystyle\forall i$

where the first set of constraints can be interpreted as a multivariate
regression equation, the second set of constraints imposes concavity on the
quantile function, the third set of constraints guarantees monotonicity, and
the last refers to sign constraints of the error terms. Note that there exists
an intrinsic quantile property in terms of the optimal solutions to problem
(3), $\hat{\varepsilon}_{i}^{+}$and $\hat{\varepsilon}_{i}^{-}$.

######  Theorem 1.

For any $\tau\in(0,1)$, the number of strict positive residuals
($\hat{\varepsilon}_{i}^{+}>0$) by $n_{\tau}^{+}$ and the number of strict
negative residuals ($\hat{\varepsilon}_{i}^{-}>0$) by $n_{\tau}^{-}$ always
satisfy the inequalities:

$\frac{n_{\tau}^{+}}{n}\leq
1-\tau\quad\mbox{and}\quad\frac{n_{\tau}^{-}}{n}\leq\tau.$

###### Proof.

See proofs in Wang2014c and Kuosmanen2020b. ∎

Compared with the conventional full frontier estimation, the quantile function
estimation is more robust to random noise, heteroscedasticity, and the choice
of direction vectors. However, when separately estimating each conditional
quantile function $Q_{y}(\tau\,|\,{\boldsymbol{x}})$, CQR is likely to violate
the assumption that the distribution functions and their associated inverse
functions should be monotone nondecreasing; see Fig. A1 for an example of the
quantile crossing problem detected in our empirical application of CQR.

We notice that the quantile crossing problem could be addressed by
simultaneous estimation, which imposes an extra set of linear non-crossing
constraints in the CQR approach (see, e.g., Takeuchi2006, *Takeuchi2006;
Wang2014c, *Wang2014c). Following Wang2014c, the simultaneous convex quantile
regression (sCQR) estimator of $j$ conditional quantile functions at
$0<\tau_{1}<\tau_{2}<\cdots<\tau_{j}<1$ is formulated as

$\displaystyle\underset{\alpha,{\boldsymbol{\beta}},\varepsilon^{\text{+}},\varepsilon^{-},C}{\mathop{\min}}$
$\displaystyle\,\sum\limits_{j=1}^{{\mathcal{J}}}\Big{(}\tau_{j}\sum\limits_{i=1}^{n}{\varepsilon_{i,j}^{+}}+(1-\tau_{j})\sum\limits_{i=1}^{n}{\varepsilon_{i,j}^{-}}\Big{)}$
(4) s.t. $\displaystyle
y_{i,j}=\mathbf{\alpha}_{i,j}+{\boldsymbol{\beta}}_{i,j}^{{}^{\prime}}{\boldsymbol{x}}_{i}+\varepsilon_{i,j}^{+}-\varepsilon_{i,j}^{-}$
$\displaystyle\forall i,j$
$\displaystyle\alpha_{i,j}+{\boldsymbol{\beta}}_{i,j}^{{}^{\prime}}{\boldsymbol{x}}_{i}\leq\alpha_{h,j}+{\boldsymbol{\beta}}_{h,j}^{{}^{\prime}}{\boldsymbol{x}}_{i}$
$\displaystyle\forall i,h,j$
$\displaystyle\alpha_{i,j}+{\boldsymbol{\beta}}_{i,j}^{{}^{\prime}}{\boldsymbol{x}}_{i}+C_{i,j}\leq\alpha_{i,j+1}+{\boldsymbol{\beta}}_{i,j+1}^{{}^{\prime}}{\boldsymbol{x}}_{i}$
$\displaystyle\forall i,1\leq j\leq{\mathcal{J}}-1$
$\displaystyle{\boldsymbol{\beta}}_{i,j}\geq{\boldsymbol{0}}$
$\displaystyle\forall i,j$ $\displaystyle\varepsilon_{i,j}^{+}\geq 0,\
\varepsilon_{i,j}^{-}\geq 0$ $\displaystyle\forall i,j$

where $C_{i,j}\geq 0$ are small nonnegative constants for quantiles, which are
introduced in sCQR to ensure that $Q_{y}(\tau_{j}\,|\,{\boldsymbol{x}})\leq
Q_{y}(\tau_{j+1}\,|\,{\boldsymbol{x}}),\forall i\text{ }\text{and}\text{
}j\in{\mathcal{J}}$. For the purpose of non-crossing, $C$ can simply be given
by zero; that is, there may exist touching rather than crossing between two
neighboring quantiles (see Fig. 1(a) for an illustration). In practice,
however, after enforcing the non-crossing constraints, sCQR may violate the
quantile property (Theorem 1) due to the fact that the approach simultaneously
optimizes for both the quantile property and the non-crossing property
(Takeuchi2006, *Takeuchi2006).

This paper proposes an alternative to sCQR to address the quantile crossing
problem. By using the $L_{2}$-norm regularization on subgradients
${\boldsymbol{\beta}}_{i}$, we formulate penalized convex quantile regression
(pCQR) as

$\displaystyle\underset{\alpha,{\boldsymbol{\beta}},\varepsilon^{+},\varepsilon^{-}}{\mathop{\min}}$
$\displaystyle\,\tau\sum\limits_{i=1}^{n}{\varepsilon_{i}^{+}}+(1-\tau)\sum\limits_{i=1}^{n}{\varepsilon_{i}^{-}}+\gamma\sum\limits_{i=1}^{n}||{\boldsymbol{\beta}}_{i}||^{2}_{2}$
(5) s.t. $\displaystyle
y_{i}=\alpha_{i}+{\boldsymbol{\beta}}_{i}^{{}^{\prime}}{\boldsymbol{x}}_{i}+\varepsilon_{i}^{+}-\varepsilon_{i}^{-}$
$\displaystyle\forall i$
$\displaystyle\alpha_{i}+{\boldsymbol{\beta}}_{i}^{{}^{\prime}}{\boldsymbol{x}}_{i}\leq\alpha_{h}+{\boldsymbol{\beta}}_{h}^{{}^{\prime}}{\boldsymbol{x}}_{i}$
$\displaystyle\forall i,h$
$\displaystyle{\boldsymbol{\beta}}_{i}\geq{\boldsymbol{0}}$
$\displaystyle\forall i$ $\displaystyle\varepsilon_{i}^{+}\geq 0,\
\varepsilon_{i}^{-}\geq 0$ $\displaystyle\forall i$

where $\gamma\geq 0$ is the tuning parameter and $||\cdot||_{2}$ denotes the
standard Euclidean norm. As $\gamma$ approaches to zero, pCQR (5) collapses to
the original CQR problem (3). The rationale behind the avoidance of quantile
crossing in pCQR lies in that as $\gamma\xrightarrow{}\infty$, the
regularization will dominate the minimization and then all estimated
subgradients ${\boldsymbol{\beta}}_{i}$ “flatten out” to 0. In this case, the
estimated quantile functions will be horizontal lines (for $d=1$) or planes
($d>1$).

As $\gamma$ increases, the quantile property may also be disturbed in pCQR.
But for sufficiently small $\gamma$, the quantile property can be guaranteed
in theory. Therefore, we design the following Algorithm 1 to obtain the
supremum $\gamma^{*}$ so as to avoid quantile crossing and ensure the quantile
property as well as possible in estimating multiple quantile functions.

Data:
$\\{{\boldsymbol{x}}_{i},y_{i}\\}_{i=1}^{n}\in{\mathbb{R}}^{d}\times{\mathbb{R}}$,
$\tau_{1}$ and $\tau_{2}$ ($\tau_{1}<\tau_{2}$)

1 $\mathrm{out}=0$ and $\gamma=0$;

2 while  _$\mathrm{out}=0$_ do

3 Solve problem (5) with quantiles $\tau_{1}$ and $\tau_{2}$, separately, to
calculate $\hat{Q}_{y}(\tau_{1}\,|\,{\boldsymbol{x}}_{i})$ and
$\hat{Q}_{y}(\tau_{2}\,|\,{\boldsymbol{x}}_{i})$;

4 if
_$\sum\limits_{i=1}^{n}\mathbf{1}_{\\{\hat{Q}_{y}(\tau_{1}\,|\,{\boldsymbol{x}}_{i})-\hat{Q}_{y}(\tau_{2}\,|\,{\boldsymbol{x}}_{i})\leq
0\\}}\neq n$_ then

5 Re-solve problem (5) using the updated $\gamma$;

6

7 else

8 $\mathrm{out}=1$;

9

10 $\gamma=\gamma+0.01$;

11

Result: $\gamma^{*}$

Algorithm 1 Searching the optimal tuning parameter $\gamma^{*}$.

We proceed to illustrate how non-crossing quantile functions look like with a
real dataset used in Kuosmanen2020b. It contains plant-level data on 130 U.S.
electric power plants operating in 2014; see Kuosmanen2020b for a more
detailed description of the data. For the sake of demonstration, we simply
consider a univariate case of one input and one output. The input is the total
cost involved in electricity production and the output is the net electricity
generation of each power plant. Both variables are in natural logarithm.

An application of CQR to the empirical data finds that the 15${}^{\text{th}}$
quantile curve crosses the 25${}^{\text{th}}$ quantile curve twice (see Fig.
A1). We then demonstrate how the sCQR and pCQR approaches can address this
problem. It is evident from Fig. 1 that both approaches manage to circumvent
the quantile crossing problem; that is, we observe that
$\widehat{Q}_{y}(0.25\,|\,{\boldsymbol{x}}_{i})$ is greater than or equal to
$\widehat{Q}_{y}(0.15\,|\,{\boldsymbol{x}}_{i})$ in both approaches. However,
the shapes of the estimated quantile functions in Figs. 1(a) and 1(b) (see
particularly the upper right corner) are slightly different. As mentioned
earlier, this difference arises because sCQR tries to simultaneously optimize
for the non-crossing property, the quantile property, and the production
axioms, whereas the pCQR approach independently estimates the quantile
production functions. Further, the difference affects which approach can
better retain the quantile property.

##  3 Monte Carlo study

We perform a Monte Carlo study to examine whether pCQR or sCQR can better
satisfy the quantile property while addressing quantile crossing. Consider the
following data generating process (Dai2021a, *Dai2021a)

$y_{i}=\prod\limits_{d=1}^{D}{\boldsymbol{x}}^{\frac{0.8}{d}}_{d,i}+v_{i}-u_{i},$

where the input matrix ${\boldsymbol{x}}_{i}\in{\mathbb{R}}^{n\times d}$ is
generated independently from $U[1,10]$, and noise $v_{i}$ and inefficiency
$u_{i}$ are drawn independently from $N(0,\sigma_{v}^{2})$ and
$N^{+}(0,\sigma_{u}^{2})$, respectively. To investigate the robustness of the
quantile approaches, following Aigner1977 we use different combinations of
noise ($\sigma$) and signal to noise ratio ($\lambda$), that is,
($\sigma^{2}$, $\lambda$) = (1.88, 1.66), (1.63, 1.24), and (1.35, 0.83).

We consider 54 scenarios with $n\in\\{99,199,499\\}$, $d\in\\{2,3,4\\}$, and
$\tau\in\\{0.85,0.90,0.95\\}$. Each scenario is replicated 500 times using the
pyStoNED package (Dai2021b, *Dai2021b) on Python with the standard solver
Mosek (9.3). We then compute the ramp loss
($\text{RL}=|\frac{1}{n}\sum_{i=1}^{n}\textbf{1}_{y_{i}>{{\widehat{Q}}_{y}}(\tau\,|\,{\boldsymbol{x}}_{i})}-\tau|$)
(Takeuchi2006, *Takeuchi2006) to examine the quantile property and the mean
squared error (MSE) to evaluate the finite-sample performance. Replications
for which no quantile crossing happen (i.e., $\gamma^{*}=0$) are excluded from
the calculations of RL and MSE. Note that the smaller the ramp loss, the
better the quantile performance.

Tables 1 and A1 present the estimated ramp loss and MSE statistics across
different scenarios. The results clearly show that compared with sCQR, pCQR
has lower ramp loss in virtually all the scenarios and lower MSE in all the
scenarios. This is because the optimal tuning parameter $\gamma^{*}$ used in
the pCQR simulations mainly locates in the interval $[0.01,0.03]$ (see Fig.
2), suggesting that the quantile property can be guaranteed to a certain
extent and thus pCQR can better fit the true quantile functions. Several other
findings are summarized as follows:

  * •

The higher the dimension or the noise in data space, the lower the ramp loss.
As $d$ or $\sigma$ increases, the data space becomes more sparse, thereby
indicating that the probability of crossing between two neighboring quantiles
is relatively small.

  * •

The differences in the ramp loss among quantiles in pCQR are smaller than
those in sCQR due to the different estimation strategies, i.e., independent
and simultaneous estimation, respectively.

  * •

For both approaches, the MSE increases as more inputs are included and
decreases as the sample size gets larger.

  * •

The performance of both approaches in terms of MSE becomes worse as the signal
to noise ratio increases.

Overall, the pCQR approach can better satisfy the quantile property and fit
the true quantile functions while at the same time addressing the quantile
crossing problem. Regularizing the quantile function, instead of imposing an
extra set of non-crossing constraints, proves a better remedy to quantile
crossing according to our simulations.

(a) sCQR ($C=0$)

(b) pCQR ($\gamma^{*}=0.01$)

Figure 1: Empirical illustration of estimated non-crossing quantile functions using the U.S. power plant data. Table 1: Ramp loss and MSE of two neighboring quantiles with $n=499$. $d$ | ($\sigma^{2},\lambda$) | $\tau_{1}$ | $\tau_{2}$ |  RL${}_{\tau_{1}}$ |  MSE${}_{\tau_{1}}$ |  RL${}_{\tau_{2}}$ |  MSE${}_{\tau_{2}}$
---|---|---|---|---|---|---|---
pCQR | sCQR | pCQR | sCQR | pCQR | sCQR | pCQR | sCQR
2 | (1.35, 0.83) | 0.85 | 0.90 | 0.928 | 0.993 | 0.065 | 0.093 | 0.968 | 1.035 | 0.075 | 0.115
|  | 0.90 | 0.95 | 0.998 | 1.078 | 0.074 | 0.117 | 1.035 | 1.108 | 0.095 | 0.176
| (1.63, 1.24) | 0.85 | 0.90 | 0.899 | 0.953 | 0.075 | 0.108 | 0.933 | 1.028 | 0.085 | 0.135
|  | 0.90 | 0.95 | 0.972 | 1.061 | 0.081 | 0.134 | 0.991 | 1.051 | 0.108 | 0.208
| (1.88, 1.66) | 0.85 | 0.90 | 0.898 | 0.967 | 0.079 | 0.117 | 0.984 | 1.069 | 0.095 | 0.155
|  | 0.90 | 0.95 | 0.890 | 0.963 | 0.093 | 0.153 | 0.923 | 1.012 | 0.117 | 0.228
3 | (1.35, 0.83) | 0.85 | 0.90 | 0.811 | 0.942 | 0.103 | 0.205 | 0.826 | 0.947 | 0.115 | 0.265
|  | 0.90 | 0.95 | 0.775 | 0.902 | 0.122 | 0.253 | 0.745 | 0.860 | 0.166 | 0.402
| (1.63, 1.24) | 0.85 | 0.90 | 0.832 | 0.950 | 0.110 | 0.226 | 0.832 | 0.940 | 0.127 | 0.296
|  | 0.90 | 0.95 | 0.827 | 0.959 | 0.138 | 0.281 | 0.751 | 0.878 | 0.175 | 0.435
| (1.88, 1.66) | 0.85 | 0.90 | 0.792 | 0.913 | 0.127 | 0.250 | 0.827 | 0.960 | 0.144 | 0.325
|  | 0.90 | 0.95 | 0.789 | 0.914 | 0.172 | 0.332 | 0.749 | 0.899 | 0.214 | 0.512
4 | (1.35, 0.83) | 0.85 | 0.90 | 0.658 | 0.837 | 0.163 | 0.368 | 0.671 | 0.833 | 0.187 | 0.484
|  | 0.90 | 0.95 | 0.627 | 0.794 | 0.207 | 0.487 | 0.481 | 0.601 | 0.269 | 0.750
| (1.63, 1.24) | 0.85 | 0.90 | 0.587 | 0.801 | 0.168 | 0.397 | 0.627 | 0.795 | 0.193 | 0.526
|  | 0.90 | 0.95 | 0.696 | 0.861 | 0.217 | 0.535 | 0.539 | 0.700 | 0.279 | 0.828
| (1.88, 1.66) | 0.85 | 0.90 | 0.584 | 0.793 | 0.189 | 0.446 | 0.622 | 0.793 | 0.218 | 0.589
|  | 0.90 | 0.95 | 0.646 | 0.795 | 0.285 | 0.581 | 0.504 | 0.640 | 0.381 | 0.920
Figure 2: Empirical distribution of $\gamma^{*}$ with $n=499$,
$(\sigma^{2},\lambda)=(1.35,0.83)$, $\tau_{1}=0.85$, and $\tau_{2}=0.90$.

##  4 Conclusions

In this paper, a penalized convex quantile regression approach has been
developed to address the quantile crossing problem. The proposed algorithm can
search the optimal tuning parameter such that the occurrences of quantile
crossing are avoided and the quantile property is ensured as well as possible.
A Monte Carlo study confirms the superiority of the proposed approach compared
to previous work in addressing the quantile crossing problem. We believe the
proposed approach can be readily extended from convex quantile regression to
other nonparametric and parametric quantile regression techniques. We leave
such extensions as fascinating avenues for future research.

## References

  * [1] [Aigner et al.]Aigner, Lovell Schmidt1977Aigner1977 Aigner, D., Lovell, C. A. K. Schmidt, P. 1977, Formulation and estimation of stochastic frontier production function models, Journal of Econometrics 6, 21–37.
  * [2] [Bondell et al.]Bondell, Reich Wang2010Bondell2010 Bondell, H. D., Reich, B. J. Wang, H. 2010, Noncrossing quantile regression curve estimation, Biometrika 97, 825–838.
  * [3] [Chernozhukov et al.]Chernozhukov, FernÁndez-Val Galichon2010Chernozhukov2010 Chernozhukov, V., FernÁndez-Val, I. Galichon, A. 2010, Quantile and probability curves without crossing, Econometrica 78, 1093–1125.
  * [4] Dai2021Dai2021a Dai, S. 2021, Variable selection in convex quantile regression: L1-norm or L0-norm regularization?, arxiv preprint arxiv:2107.03119.
  * [5] [Dai et al.]Dai, Fang, Lee Kuosmanen2021Dai2021b Dai, S., Fang, Y. H., Lee, C. Y. Kuosmanen, T. 2021, pyStoNED: A Python package for convex regression and frontier estimation, arxiv preprint arxiv:2109.12962.
  * [6] Dette Volgushev2008Dette2008 Dette, H. Volgushev, S. 2008, Non-crossing non-parametric estimates of quantile curves, Journal of the Royal Statistical Society. Series B (Statistical Methodology) 70, 609–627.
  * [7] He1997He1997 He, X. 1997, Quantile curves without crossing, American Statistician 51, 186–192.
  * [8] [Jradi et al.]Jradi, Parmeter Ruggiero2019Jradi2019 Jradi, S., Parmeter, C. F. Ruggiero, J. 2019, Quantile estimation of the stochastic frontier model, Economics Letters 182, 15–18.
  * [9] Kuosmanen Zhou2021Kuosmanen2020b Kuosmanen, T. Zhou, X. 2021, Shadow prices and marginal abatement costs: Convex quantile regression approach, European Journal of Operational Research 289, 666–675.
  * [10] Liu Wu2009Wu2009a Liu, Y. Wu, Y. 2009, Stepwise multiple quantile regression estimation using non-crossing constraints, Statistics and Its Interface 2, 299–310.
  * [11] [Takeuchi et al.]Takeuchi, Le, Sears Smola2006Takeuchi2006 Takeuchi, I., Le, Q. V., Sears, T. D. Smola, A. J. 2006, Nonparametric quantile estimation, Journal of Machine Learning Research 7, 1231–1264.
  * [12] [Tsionas et al.]Tsionas, Assaf Andrikopoulos2020Tsionas2020b Tsionas, M. G., Assaf, A. G. Andrikopoulos, A. 2020, Quantile stochastic frontier models with endogeneity, Economics Letters 188\.
  * [13] [Wang et al.]Wang, Wang, Dang Ge2014Wang2014c Wang, Y., Wang, S., Dang, C. Ge, W. 2014, Nonparametric quantile frontier estimation under shape restriction, European Journal of Operational Research 232, 671–678.
  * [14] Zhao2021Zhao2021 Zhao, S. 2021, Quantile estimation of stochastic frontier models with the normal-half normal specification: A cumulative distribution function approach, Economics Letters 206, 109998.
  * [15]

## Appendix

###  A Supplementary tables and figures

Figure A1: Illustration of the quantile crossing problem in CQR estimation.

Table A1: Ramp loss and MSE of two neighboring quantiles with $n=99$ and $199$. $n$ | $d$ | ($\sigma^{2},\lambda$) | $\tau_{1}$ | $\tau_{2}$ |  RL${}_{\tau_{1}}$ |  MSE${}_{\tau_{1}}$ |  RL${}_{\tau_{2}}$ |  MSE${}_{\tau_{2}}$
---|---|---|---|---|---|---|---|---
pCQR | sCQR | pCQR | sCQR | pCQR | sCQR | pCQR | sCQR
99 | 2 | (1.35, 0.83) | 0.85 | 0.90 | 0.824 | 0.882 | 0.239 | 0.350 | 0.869 | 0.884 | 0.276 | 0.444
|  |  | 0.90 | 0.95 | 0.849 | 0.899 | 0.300 | 0.459 | 0.835 | 0.794 | 0.396 | 0.683
|  | (1.63, 1.24) | 0.85 | 0.90 | 0.806 | 0.882 | 0.248 | 0.378 | 0.857 | 0.904 | 0.287 | 0.484
|  |  | 0.90 | 0.95 | 0.888 | 0.935 | 0.324 | 0.481 | 0.867 | 0.845 | 0.427 | 0.727
|  | (1.88, 1.66) | 0.85 | 0.90 | 0.831 | 0.898 | 0.272 | 0.435 | 0.891 | 0.916 | 0.327 | 0.556
|  |  | 0.90 | 0.95 | 0.857 | 0.918 | 0.336 | 0.534 | 0.792 | 0.809 | 0.439 | 0.803
| 3 | (1.35, 0.83) | 0.85 | 0.90 | 0.719 | 0.753 | 0.349 | 0.619 | 0.667 | 0.702 | 0.400 | 0.777
|  |  | 0.90 | 0.95 | 0.626 | 0.687 | 0.442 | 0.807 | 0.446 | 0.381 | 0.618 | 1.262
|  | (1.63, 1.24) | 0.85 | 0.90 | 0.671 | 0.721 | 0.400 | 0.688 | 0.632 | 0.657 | 0.469 | 0.878
|  |  | 0.90 | 0.95 | 0.629 | 0.695 | 0.516 | 0.875 | 0.393 | 0.364 | 0.713 | 1.392
|  | (1.88, 1.66) | 0.85 | 0.90 | 0.677 | 0.741 | 0.439 | 0.782 | 0.622 | 0.676 | 0.500 | 0.989
|  |  | 0.90 | 0.95 | 0.611 | 0.662 | 0.541 | 1.046 | 0.437 | 0.404 | 0.742 | 1.651
| 4 | (1.35, 0.83) | 0.85 | 0.90 | 0.412 | 0.491 | 0.591 | 1.041 | 0.334 | 0.342 | 0.715 | 1.384
|  |  | 0.90 | 0.95 | 0.382 | 0.421 | 0.798 | 1.381 | 0.239 | 0.194 | 1.118 | 2.148
|  | (1.63, 1.24) | 0.85 | 0.90 | 0.455 | 0.550 | 0.688 | 1.139 | 0.403 | 0.427 | 0.833 | 1.497
|  |  | 0.90 | 0.95 | 0.381 | 0.395 | 0.874 | 1.491 | 0.229 | 0.182 | 1.264 | 2.384
|  | (1.88, 1.66) | 0.85 | 0.90 | 0.502 | 0.606 | 0.664 | 1.180 | 0.450 | 0.468 | 0.813 | 1.567
|  |  | 0.90 | 0.95 | 0.412 | 0.449 | 0.931 | 1.624 | 0.239 | 0.231 | 1.339 | 2.567
199 | 2 | (1.35, 0.83) | 0.85 | 0.90 | 0.837 | 0.916 | 0.139 | 0.201 | 0.914 | 0.964 | 0.160 | 0.251
|  |  | 0.90 | 0.95 | 0.914 | 0.997 | 0.159 | 0.250 | 0.934 | 0.981 | 0.206 | 0.372
|  | (1.63, 1.24) | 0.85 | 0.90 | 0.869 | 0.957 | 0.140 | 0.210 | 0.935 | 0.978 | 0.166 | 0.276
|  |  | 0.90 | 0.95 | 0.944 | 1.016 | 0.178 | 0.291 | 0.959 | 1.003 | 0.233 | 0.427
|  | (1.88, 1.66) | 0.85 | 0.90 | 0.868 | 0.931 | 0.150 | 0.236 | 0.894 | 0.949 | 0.186 | 0.311
|  |  | 0.90 | 0.95 | 0.945 | 1.012 | 0.186 | 0.322 | 0.948 | 0.968 | 0.223 | 0.459
| 3 | (1.35, 0.83) | 0.85 | 0.90 | 0.656 | 0.777 | 0.203 | 0.374 | 0.703 | 0.779 | 0.231 | 0.485
|  |  | 0.90 | 0.95 | 0.733 | 0.830 | 0.256 | 0.472 | 0.561 | 0.621 | 0.338 | 0.721
|  | (1.63, 1.24) | 0.85 | 0.90 | 0.696 | 0.815 | 0.221 | 0.436 | 0.708 | 0.790 | 0.248 | 0.568
|  |  | 0.90 | 0.95 | 0.748 | 0.881 | 0.266 | 0.550 | 0.631 | 0.707 | 0.364 | 0.848
|  | (1.88, 1.66) | 0.85 | 0.90 | 0.713 | 0.837 | 0.249 | 0.461 | 0.740 | 0.839 | 0.277 | 0.595
|  |  | 0.90 | 0.95 | 0.657 | 0.716 | 0.335 | 0.653 | 0.516 | 0.553 | 0.438 | 1.013
| 4 | (1.35, 0.83) | 0.85 | 0.90 | 0.584 | 0.703 | 0.282 | 0.663 | 0.568 | 0.597 | 0.328 | 0.876
|  |  | 0.90 | 0.95 | 0.452 | 0.570 | 0.409 | 0.918 | 0.350 | 0.364 | 0.537 | 1.430
|  | (1.63, 1.24) | 0.85 | 0.90 | 0.503 | 0.670 | 0.331 | 0.736 | 0.507 | 0.605 | 0.380 | 0.969
|  |  | 0.90 | 0.95 | 0.453 | 0.569 | 0.476 | 0.974 | 0.321 | 0.354 | 0.617 | 1.533
|  | (1.88, 1.66) | 0.85 | 0.90 | 0.428 | 0.622 | 0.354 | 0.806 | 0.427 | 0.554 | 0.423 | 1.074
|  |  | 0.90 | 0.95 | 0.472 | 0.630 | 0.515 | 1.084 | 0.332 | 0.375 | 0.657 | 1.687


graph $G_{\text{2D}}$ has a bulk node (black circle) for each $Z$ stabilizer
generator, and a bulk edge (black) for each data qubit. A bulk edge coming
from a bulk node corresponds to the edge’s data qubit being in the support of
the node’s stabilizer. The graph also contains boundary nodes (white boxes)
and boundary edges (blue), which do not correspond to stabilizers or data
qubits. Each bulk and boundary edge is assigned weight one and zero
respectively. The minimum weight decoder is then implemented as follows. After
the error $E$ is applied, the nodes corresponding to unsatisfied stabilizers
are highlighted. If an odd number of stabilizers was unsatisfied, one of the
boundary nodes is also highlighted. Highlighted nodes are then efficiently
paired together by the minimum weight connections in the graph, by Edmonds’
algorithm [Edmonds(1965), Kolmogorov(2009)]. The correction $C$ is applied to
the edges in the connection. Note that any single $\mathcal{O}(p)$ fault in
this noise model corresponds to a weight one edge on the graph.

(a)

(b)

Figure 20: Circuits for measuring (a) $Z$-type, and (b) $X$-type generators.
Identity gates (black rectangles) are inserted in the $Z$-type stabilizer
measurement circuits to ensure that all measurements are synchronized. Note
that unlike in [Fowler et al.(2012)Fowler, Mariantoni, Martinis, and Cleland],
to be consistent with the other schemes in this paper, we assume that we can
prepare and measure in both the $X$ and $Z$ basis.

For circuit noise, we introduce a measurement qubit for each stabilizer
generator, as represented by gray circles in Fig. 19(a), and circuits must be
specified to implement the measurements, such as those in Fig. 20. The
performance of the code is sensitive to the choice of circuit [Tomita and
Svore(2014)], for example a poor choice could allow a single fault to cause a
logical failure for $d=3$ for any choice of decoder.

To implement the decoder, first construct a new three dimensional graph
$G_{\text{3D}}$ by stacking $d$ copies of the planar graph $G_{\text{2D}}$
that was shown in Fig. 19(b), and adding new bulk (boudnary) edges to connect
bulk (boudnary) nodes in neighboring layers. We also add additional diagonal
edges such that any single $\mathcal{O}(p)$ fault in the measurement circuits
corresponds to a weight-one edge in $G_{\text{3D}}$ (see Fig. 21). For
simplicity, we do not involve further possible optimizations such as setting
edge weights based on precise probabilities and including $X$-$Z$ correlations
[Wang et al.(2011)Wang, Fowler, and Hollenberg].

(a)

(b)

(c)

(d)

(e)

(f)

Figure 21: Examples of a single fault leading to diagonal edges in
$G_{\text{3D}}$. Dark arrows represent the CNOT sequence. (a) An $X$ error
occurs during the third time step in the CNOT gate acting on the central data
qubit. (b) During the fifth time step of this round, the $X$ error is detected
by the $Z$ type measurement qubit to the top right. (c) The $X$ error is not
detected by the bottom left $Z$ type stabilizer until the following round. (d)
An $XX$ error occurs on the third CNOT of an $X$ measurement circuit, which is
detected by the $Z$ measurement to the right. (e) Detection by the left $Z$
stabilizer does not occur until the next round. (f) The corresponding edges in
$G_{\text{3D}}$, green for (a-c), and blue for (d-e). Here we show two rounds
of the graph ignoring boundary edges.

All simulations of the surface code are performed using the circuit noise
model in Section 1.2, with the graph $G_{\text{3D}}$ described above as
follows (to correct $X$ errors):

  1. 1.

Data acquisition: Stabilizer outcomes are stored over $d$ rounds of noisy
error correction, followed by one round of perfect error correction. The net
error $E$ applied over all $d$ rounds is recorded.

  2. 2.

Highlight nodes: Nodes in the graph $G_{\text{3D}}$ are highlighted if the
corresponding $Z$-type stabilizer outcome changes in two consecutive rounds.
555For an odd number of highlighted vertices, highlight the boundary vertex.

  3. 3.

Minimum weight matching: Find a minimal edge set forming paths that terminate
on highlighted nodes. Highlight the edge set.

  4. 4.

Vertical collapse: The highlighted edges in $G_{\text{3D}}$ are mapped edges
in the planar graph $G_{\text{2D}}$, and are then added modulo 2.

  5. 5.

Correction: The $X$-type correction $C_{X}$ is applied to highlighted edges in
$G_{\text{2D}}$.

The $Z$ correction $C_{Z}$ is found analogously. Finally, if the residual
Pauli $R=EC_{X}C_{Z}$ is a logical operator, we say the protocol succeeded,
otherwise we say it failed.

FTEC scheme | Noise model | Number of qubits | Time steps ($T_{\mathrm{time}}$) | Pseudo-threshold
---|---|---|---|---
Flag-EC $[\\![7,1,3]\\!]$ | $\tilde{p}=p$ | 9 | $36\leq T_{\mathrm{time}}\leq 108$ | $p_{\mathrm{pseudo}}=(3.39\pm 0.10)\times 10^{-5}$
Flag-EC $[\\![7,1,3]\\!]$ |  | 11 | $34\leq T_{\mathrm{time}}\leq 104$ | $p_{\mathrm{pseudo}}=(2.97\pm 0.01)\times 10^{-5}$
Table 6: Pseudo-thresholds and circuit depth for flag-EC protocols using two
and four ancilla qubits applied to the $[\\![7,1,3]\\!]$ code. The results are
presented for the noise models where $\tilde{p}=p$ and $\tilde{p}=p/100$.

##  Appendix G Compact implementation of flag error correction

$[\\![5,1,3]\\!]$ | $[\\![7,1,3]\\!]$ |  $[\\![19,1,5]\\!]$ code |  $[\\![17,1,5]\\!]$ code
---|---|---|---
$X_{1}Z_{2}Z_{3}X_{4}$ | $Z_{4}Z_{5}Z_{6}Z_{7}$ |  $Z_{1}Z_{2}Z_{3}Z_{4}$, $X_{1}X_{2}X_{3}X_{4}$ |  $Z_{1}Z_{2}Z_{3}Z_{4}$, $X_{1}X_{2}X_{3}X_{4}$
$X_{2}Z_{3}Z_{4}X_{5}$ | $Z_{2}Z_{3}Z_{6}Z_{7}$ |  $Z_{1}Z_{3}Z_{5}Z_{7}$, $X_{1}X_{3}X_{5}X_{7}$ |  $Z_{1}Z_{3}Z_{5}Z_{6}$, $X_{1}X_{3}X_{5}X_{6}$
$X_{1}X_{3}Z_{4}Z_{5}$ | $Z_{1}Z_{3}Z_{5}Z_{7}$ |  $Z_{12}Z_{13}Z_{14}Z_{15}$, $X_{12}X_{13}X_{14}X_{15}$ |  $Z_{5}Z_{6}Z_{9}Z_{10}$, $X_{5}X_{6}Z_{9}Z_{10}$
$Z_{1}X_{2}X_{4}Z_{5}$ | $X_{4}X_{5}X_{6}X_{7}$ |  $Z_{1}Z_{2}Z_{5}Z_{6}Z_{8}Z_{9}$, $X_{1}X_{2}X_{5}X_{6}X_{8}X_{9}$ |  $Z_{7}Z_{8}Z_{11}Z_{12}$, $X_{7}X_{8}X_{11}X_{12}$
| $X_{2}X_{3}X_{6}X_{7}$ |  $Z_{6}Z_{9}Z_{16}Z_{19}$, $X_{6}X_{9}X_{16}X_{19}$ |  $Z_{9}Z_{10}Z_{13}Z_{14}$, $X_{9}X_{10}X_{13}X_{14}$
| $X_{1}X_{3}X_{5}X_{7}$ |  $Z_{16}Z_{17}Z_{18}Z_{19}$, $X_{16}X_{17}X_{18}X_{19}$ |  $Z_{11}Z_{12}Z_{15}Z_{16}$, $X_{11}X_{12}X_{15}X_{16}$
|  |  $Z_{10}Z_{11}Z_{12}Z_{15}$, $X_{10}X_{11}X_{12}X_{15}$ |  $Z_{8}Z_{12}Z_{16}Z_{17}$, $X_{8}X_{12}X_{16}X_{17}$
|  | $Z_{8}Z_{9}Z_{10}Z_{11}Z_{16}Z_{17}$ | $Z_{3}Z_{4}Z_{6}Z_{7}Z_{10}Z_{11}Z_{14}Z_{15}$
|  | $Z_{5}Z_{7}Z_{8}Z_{11}Z_{12}Z_{13}$ | $X_{3}X_{4}X_{6}X_{7}X_{10}X_{11}X_{14}X_{15}$
|  | $X_{5}X_{7}X_{8}X_{11}X_{12}X_{13}$ |
|  | $X_{8}X_{9}X_{10}X_{11}X_{16}X_{17}$ |
$\overline{X}=X^{\otimes 5}$,$\overline{Z}=Z^{\otimes 5}$ |  $\overline{X}=X^{\otimes 7}$,$\overline{Z}=Z^{\otimes 7}$ |  $\overline{X}=X^{\otimes 19}$,$\overline{Z}=Z^{\otimes 19}$ |  $\overline{X}=X^{\otimes 17}$,$\overline{Z}=Z^{\otimes 17}$
Table 7: Stabilizer generators for the 5-qubit code [Laflamme et
al.(1996)Laflamme, Miquel, Paz, and Zurek], $d=3$ (Steane code)
[Steane(1996)], $d=5$ ($[\\![19,1,5]\\!]$ code) and $d=5$ ($[\\![17,1,5]\\!]$
code) family of color codes [Bombin and Martin-Delgado(2006)]. The last row
illustrates representatives of the codes logical operators. Figure 22: Circuit
for measuring the $Z$ stabilizer generators of the $[\\![7,1,3]\\!]$ code
using one flag qubit and three measurement qubits. The circuit is constructed
such that any single fault at a bad location leading to an error of weight
greater than one will cause the circuit to flag. Moreover, any error that
occurs when the circuit flags due to a single fault has a unique syndrome.
Figure 23: Logical failure rates of the flag 1-FTEC protocols using two and
four ancilla qubits applied to the $[\\![7,1,3]\\!]$ Steane code.

In [Chao and Reichardt(2017a)], it was shown that by using extra ancilla
qubits in the flag-EC protocol, it is possible to measure multiple stabilizer
generators during one measurement cycle which could reduce the circuit depth.
Note that for the Steane code, measuring the $Z$ stabilizers using Fig. 2(b)
requires only one extra time step. In this section we compare logical failure
rates of the $[\\![7,1,3]\\!]$ code using the flag-EC method of Section 2.1
which requires only two ancilla qubits and a flag-EC method which uses four
ancilla qubits but that can measure all $Z$ stabilizer generators in one cycle
(see Fig. 22). All $X$ stabilizers are measured in a separate cycle.

Logical failure rates for $\tilde{p}=p$ are shown in Fig. 23. Pseudo-
thresholds and the number of time steps required to implement the protocols
are given in Table 6. Note that measuring stabilizers using two ancilla’s
requires at most two extra time steps. Furthermore, the extra ancilla’s for
measuring multiple stabilizers result in more idle qubit locations compared to
using only two ancilla qubits. With the added locations for errors to be
introduced, the flag error correction protocol using only two ancilla’s
achieves a higher pseudo-threshold compared to the protocol using more
ancilla’s. Thus assuming that reinitializing qubits can be done without
introducing many errors into the system, FTEC using fewer qubits could achieve
lower logical failure rates compared to certain schemes using more qubits.

##  Appendix H Stabilizer generators of various codes.

In Table 7 we provide stabilizer generators for the $[\\![5,1,3]\\!]$ code,
$[\\![7,1,3]\\!]$ Steane code, $[\\![19,1,5]\\!]$ and $[\\![17,1,5]\\!]$ color
codes.