# Atomically Sharp 318nm Gd:AlGaN Ultraviolet Light Emitting Diodes on Si with
Low Threshold Voltage

Thomas F. Kent  Department of Materials Science and Engineering, The Ohio
State University, Columbus, Ohio 43210, USA  Santino D. Carnevale  Department
of Materials Science and Engineering, The Ohio State University, Columbus,
Ohio 43210, USA  Roberto C. Myers  Department of Materials Science and
Engineering, The Ohio State University, Columbus, Ohio 43210, USA  Deparment
of Electrical and Computer Engineering, The Ohio State University, Columbus,
Ohio 43210, USA

###### Abstract

Self-assembled AlxGa1-xN polarization-induced nanowire light emitting diodes
(PINLEDs) with Gd-doped AlN active regions are prepared by plasma-assisted
molecular beam epitaxy on Si substrates. Atomically sharp electroluminescence
(EL) from Gd intra-f-shell electronic transitions at 313nm and 318nm is
observed under forward biases above 5V. The intensity of the Gd 4f EL scales
linearly with current density and increases at lower temperature. The low
field excitation of Gd 4f EL in PINLEDs is contrasted with high field
excitation in metal/Gd:AlN/polarization-induced n-AlGaN devices; PINLED
devices offer over a three fold enhancement in 4f EL intensity at a given
device bias.

Keywords: AlGaN, Nanowires, Rare Earth, Electroluminescent Devices

In this letter we report on ultraviolet light emitting diodes based on self-
assembled AlGaN nanowire heterojunctions doped with Gd, which emit ultraviolet
(UV) radiation at 318nm and operate at lower device bias compared with
existing Gd electroluminescent technology. Following the realization of the
ruby laser over half a century agomaiman , many materials have been developed
as hosts for elements with optically-active energetic transitions to produce
lasers and electroluminescent devices with emission in the visible and
infrared part of the electromagnetic spectrum. In the last twenty years,
electrically driven devices utilizing rare earth phosphors in wide or medium
gap semiconductors have begun to make an appearancerack . These thin film
electroluminescent devices (TFED) have found applications as visible and IR
emitters, and with the increasing use of wide gap materials, such as GaN and
AlN are currently being developed to take of advantage of phosphors in the
blue and UV.

Figure 1: Optically active Gd3+ 4f levels in the active region of a III-N
nanowire polarization induced light emitting diode. a. Calculated device band
diagram showing location of optically active Gd ions. Insert shows device
schematic as measured. b. Device IV exhibiting turn on in forward bias at 5V
Figure 2: Optically active Gd ions in the depletion region of a polarization-
induced light emitting diode (PINLED). a. Room temperature device EL spectra
for multiple currents. b. Detail of spectral region of interest for the Gd 4f
shell EL peaks, showing first excited and second excited state transitions.
Insert shows behavior of peak intensity with increasing current density. c.
Detail of spectra region of interest for the Gd 4f EL at multiple
temperatures. Insert shows variation of peak intensity with temperature.

Design of optoelectronic devices which utilize atomic transitions in the rare
earths for EL can offer a number of advantages over EL produced from band-to-
band recombination. First, the transition energy is dictated by the energy
level scheme of the 4f orbital, which is relatively unperturbed by the
crystalline environment due to the fact that the 5d and 5s,5p orbitals extend
to further radial distances, and fill before the f shell of lower principal
quantum numberatkins . This shielding by the earlier filled 5d and 5s orbitals
causes the energies of the 4f transitions to be relatively insensitive to
crystalline imperfections, unlike transitions based on band-to-band
transitions in semiconductors. Band-to-band optical transitions are well
knowndavies02  to be sensitive to deep levels, exciton-phonon interaction, and
crystalline disorder which can lead to broadening of emission or parasitic
emission at an unintended energy. Additionally, due to the decoupling of the
4f orbital with the lattice, emission from rare earth centers is spectrally
pure with common FWHM of less than 30meV.

In recent years, the family of III-nitrides, particularly the pseudobinary
AlxGa1-xN, has become attractive as a potential host for rare earth elemental
phosphorswakahara . This material system offers a high breakdown field
(12MV/cm for AlN) and sufficiently wide bandgaps (3.4 to 6.1eV) to accommodate
phosphors with emission from the UV to IR. Equally significant is the way in
which rare earths incorporate in the wurtzite structure of GaN and AlN alloys.
Rare earth (RE) ions commonly exhibit the RE3+ ionization state, which makes
them isovalent with the Al3+ and Ga3+ cations of AlxGa1-xN. RE atoms exhibit
high solubility (in excess of 1at% has been reportedwang ) and regularly
incorporate at the cation site. Due to the wide range of energy level schemes
of the lanthanides, EL from RE centers in AlN and GaN have been reported for
Er3+ (IR, green)Lee , Tm3+(blue)Lee , Eu3+(Red)wakahara , and Gd3+(UV)kita .
Additonally, room temperature optically pumped lasing in Eu-doped GaNpark  has
been successfully demonstrated.

Although most RE phosphors have been developed for optical transitions in the
visible and infrared parts of the electromagnetic spectrum, the energy scheme
of Gd3+ in wurtzite AlN offers an energy difference between ground and first
excited state of 3.90eV (318nm). The spectrally-narrow and energetically-
stable nature of the Gd3+ fluorescence emission make it a potential candidate
for spectroscopic and lithographic applications in the UV. This led to
exploration of dilutely Gd doped AlxGa1-xN in the form of fluorescencezavada
and cathodoluminescence experimentsvetter ; gruber ; kita .Although the 4f
levels in the RE3+ are typically thought not to interact with the surrounding
lattice, cathodoluminescence data for Gd:AlN thin films show phonon replica
satellite peaks of the Gd3+ 6P7/2$\rightarrow$8S7/2 (318nm) transitionsvetter
. These data suggest that the 4f electrons in Gd3+ in AlN are not completely
decoupled from the host lattice.

Although there have been a number of reportskita ; zavada ; vetter ; gruber
on the spectroscopy of Gd:AlGaN compounds, less work has focused on
development of active optoelectronic devices that utilize Gd3+ 4f transitions.
This is likely due to the difficulty of achieving electrical contact to uid-
AlN. Reportskita ; kitayama  have been made of a “field emission device”
consisting of a reactive ion sputtered AlxGd1-xN film with metal contacts,
forming a MIS structure whereby a high voltage 270V to $>$1kV driven across
the device produces fluorescence of the Gd3+ ions, likely by the process of
impact excitation.

Figure 3: Optically active Gd ions in a purely unipolar metal-insulator-
semiconductor heterostructure a. Device heterostructure and calculated band
diagram. Insert shows device schematic. b. Device IV curve exhibiting
rectification and a turn on voltage above 10V. c. EL spectra showing weak EL
of the Gd 4f 6P7/2$\rightarrow$8S7/2 transition. Insert shows detail of the 4f
EL peak with an overlay of the more intense spectrum from the PINLED device
under comparable conditions.

The dominant mechanism of RE3+ 4f excitation in solids depends on the electric
field regime where the device functions. MIS devices function largely in the
high field regime, causing the dominant mechanism to be direct transfer of
kinetic energy from hot electrons to the RE3+ center by collision, known as
impact excitation. At lower fields, the excitation mechanism becomes more
complex, involving typically a multi-step, defect assisted Auger process, or
exciton localization by the RE3+ center.godlewski ; bodiou  In this work, we
study two seperate heterostructures which are designed to generate EL under
low electric field conditions (pn-diode) and high electric field conditions
(Gd:AlGaN MIS structure). In the PINLED, the initially large built in electric
field in the depletion region is reduced with forward bias of the device. In
contrast, the MIS structure has initially flat bands and thus no electric
field in the active region and when biased, the electric field is increased in
the Gd doped region.

Self-assembled III-nitride nanowire heterostructures grown by plasma assisted
molecular beam epitaxy have recently gained popularity for applications
requiring high crystalline quality and largely mismatched, complex
heterostructures that would otherwise be difficult to form in thin films due
to strain considerationscarnevale1 ; carnevale2 ; carnevale3 . In addition,
they have been shown to function at high current densitiescarnevale3 .
Additionally, UV-LEDs based on III-nitride nanowire heterostructures have been
demonstrated to accommodate active regions spanning from high %Al AlGaN to
GaNcarnevale3 .

Polarization-induced nanowire diodes (PINLEDs)carnevale2  containing Gd doping
in their active region are prepared by plasma-assisted molecular beam epitaxy
on n-Si(111) in the III-limited growth regime, the details of which are
discussed elsewherecarnevale2 . These devices consist of a GaN nucleation
layer followed by a linear grade in composition from GaN to AlN over 100nm.
This is followed by an active region consisting of 2.4ML of GdN deposited
between two 5nm uid-AlN spacers. The structure is then linearly graded in
composition from AlN back to GaN. This structure forms a pn-diode as shown in
Fig. 1a, the band diagram of which is calculated with a self consistent
Schrödinger-Poisson solverbandeng , as described in refcarnevale2 . The
spontaneous polarization present in wurtzite AlxGa1-xN combined with a
gradient in composition give rise to polarization-induced hole and electron
doping, respectivelysimon . In addition to the built in polarization doping,
wires are supplementary doped with Mg and Si in p and n regions, assuming
c-axis,Ga-polar orientation of the nanowires.

From these heterostructures, electrical devices are fashioned by depositing
semitransparent 10nm/20nm Ni/Au contacts with an electron beam evaporator for
the top contact, which connect the vertical ensemble in parallel. The back
contact is fashioned by mechanically removing the nanowires adjacent to the
top contact with a diamond scribe and thermally diffusing In metal directly to
the n-Si with a soldering ion. Current-Voltage (IV) behavior of the devices
are measured with a probe station and an Agilent B1500 semiconductor parameter
analyzer. Device IV’s show rectifying behavior with device turn on at 5V, as
shown in Fig. 1b.

After the IV behavior of the devices is characterized, they are transferred to
a variable temperature UV-VIS-NIR spectroscopy system consisting of a closed
cycle ARS DMX20-OM cryostat and a Princeton instruments SP2500i 0.5m
spectrometer equipped with a Princeton instruments PIXIS100 UV-VIS-NIR CCD.
Devices are connected to a Yokogawa DC constant current source and Keithley
2700 data acquisition system. Prior to the collection of any spectra, a
background spectrum is collected with no current injection in the device.
Constant currents from 10mA to 120mA, corresponding to device current
densities from 1.32A/cm2 to 14.5A/cm2 (the current density through any given
nanowire is unknown since not all individual nanowires give ELlimbach ) are
sourced at room temperature and the resulting EL is collected through a 50mm,
f/2 uv-fused silica singlet lens, collimated and subsequently focused onto the
entrance slit of the 0.5m spectrometer.

Room temperature EL spectra, shown for multiple current densities in Fig. 2a.,
exhibit multiple emission peaks from the UV through the visible parts of the
spectrum. The sharp peak at 318nm (Fig. 2b.), corresponds to the
6P7/2$\rightarrow$8S7/2 first excited state to ground state transition of the
Gd3+ iongruber . More careful inspection of spectra around the Gd atomic line
region reveals an additional peak at the correct energy for the
6P5/2$\rightarrow$8S7/2 second excited state to ground state transitiongruber
. In addition, the intensity of these peaks scale linearly with current
density, within the range of currents investigated. The Gd emission exhibits
FWHM of 23.1meV. Additional broad EL peaks in the 400nm-700nm range are
attributed to below band gap defects in AlGaN, due to observation of identical
emission spectra in non-Gd containing devices, though the graded nature of
structure prevents precise identification of the deep levels responsible.

In order to further investigate the mechanism by which the UV emission occurs,
variable temperature EL measurements are conducted from room temperature to
30K, the results of which are shown in Fig. 2c. Spectral intensity of the main
6P7/2$\rightarrow$8S7/2 Gd 4f peak is observed to be invariant at temperatures
above 75K, below which the intensity increases dramatically (Fig. 2c.,
insert). Similar behavior has been previously observed in Eu-doped GaNnyein
as well as Er-doped InGaPneuhalfen  which is attributed to thermal quenching
of the multi-step excitation mechanism of the RE3+ ion. At low temperature
another peak becomes distinct at 324nm. This peak has been previously
identified as a phonon replica of the primary 318nm peak in
cathodoluminescence experimentsvetter . From analysis of the peak positions, a
phonon mode energy of 72.2meV is measured, which is smaller than the LO phonon
energy of the surrounding AlN matrix (110meVbergman ) as well as GaN
(92meVcingolani ). This phonon energy agrees within the resolution of the
spectrometer used with results from cathodoluminescence experiments, which
report 72.9meVvetter .

Many RE electroluminescent devices rely on the impact excitation mechanism to
drive intra-f-shell EL. In the interest of investigating EL under these
conditions in nanowire based devices, a structure consisting of an n-type
nanowire graded from GaN to AlN, with a 200nm uid-AlN layer doped with Gd ions
(1E18cm-3) and capped with a small amount of n++ AlN for a top contact is
prepared. Growth conditions for the heterostructure are identical to the
n-region and depletion region of the PINLED, thus similar defect content are
expected in both structures. An identical device contact scheme to that of the
heterojunction diode is used and is shown in Fig. 3a. This device is again
rectifying, as shown in Fig. 3b., but is less conductive than the PINLED
device, producing 5.9A/cm2 compared to 11.2A/cm2 when forward biased to 15V.
At 15V it is calculated that the active region of the MIS device develops an
electric field in excess of 0.75MV/cm, where as the PINLED should have
approximately flat band conditions in the active region, due to reduction of
band bending in forward bias. The lower conductivity of the MIS devices can be
attributed to the uid-AlN center region as well as a large Schottky barrier
between the n++AlN and the Ti/Au top contact which produce additional series
resistance in the device over the PINLED. Electroluminescence spectroscopy
(Fig. 3c.) reveals a weak, but detectable peak at 318nm among a large
background of defect luminescence, indicating that some Gd ions are being
excited by hot electrons passing through the structure as well as impact
excitation of band to defect luminescence. Comparison of emission from PINLED
devices and Gd:AlGaN MIS structures (Fig. 3c., insert) shows a 372%
enhancement of intensity of the 6P7/2$\rightarrow$8S7/2 transition in the
PINLED devices over the Gd:AlGaN MIS devices at 15V bias. Additionally, no
peak corresponding to the 6P5/2$\rightarrow$8S7/2 higher order transition is
present in the emission spectra from the Gd:AlGaN MIS structure. It is noted
that due to the variety of possible mechanisms for RE 4f excitation in III-V
materials, the difference in performance between the MIS device and the PINLED
device could be affected by phenomena which are extrinsic to the E-field
regime, such as preferential interaction of one carrier type with the RE ion.

In conclusion, polarization induced light emitting diodes (PINLEDS) doped with
Gd in an AlN active region are prepared by plasma assisted molecular beam
epitaxy on n-Si substrates. These devices function at one to two orders of
magnitude lower biases than previously reported Gd:AlN electroluminescent
devices, making them attractive for low power, UV EL applications,
particularly portable devices. When forward biased, devices emit sharp peaks
at 318nm and 313nm, which correspond to the Gd intra-f-shell
6P7/2$\rightarrow$8S7/2 and 6P5/2$\rightarrow$8S7/2 transitions, respectively
and scale linearly with current density. Emission intensity is shown to be
temperature independent above 75K, below which it increases strongly. By
studying two different devices, designed to produced Gd 4f EL under both low
and high electric field conditions, we observe a significant improvement in
emission intensity for PINLED devices which function under low-field operation
conditions over hot electron MIS devices. Although this device has been
applied to Gd, it would be possible in principle to dope with any of the 4f
phosphor rare earths to achieve spectrally stable electrically driven emission
at a variety of wavelengths.

This work is supported by the Center of Emergent Materials at OSU under NSF
DMR-0820414 and National Science Foundation CAREER award (DMR-1055164). S. D.
Carnevale acknowledges support from the National Science Foundation Graduate
Research Fellowship Program (2011101708).

## References

  * (1) Maimn, T. H. “Stimulated Optical Radiation in Ruby” Nature, 187 4736, pp. 493-494 (1960)
  * (2) P. D. Rack, P. H. Holloway “The structure, device physics and material properties of thin film electroluminescent displays” Mat. Sci. and Eng., R21(19988)171-219(1998)
  * (3) P. W. Atkins Molecular Quantum Mechanics 2nd Ed. Oxford, New York. 1992 pg. 223
  * (4) J. H. Davies The Physics of Low-dimensional Semiconductors: An Introduction Cambridge University Press (1997)
  * (5) A. Wakahara, H. Sekiguchi, H. Okada, Y. Takagi “Current status for light-emitting diode with Eu-doped GaN active layer grown by MBE” J. Lumin. (2012)
  * (6) R. Wang and A. J. Steckl “Effect of growth conditions on Eu${}^{3}+$ luminescence in GaN” J. Crys. Grow. 312 (2010) 680-684
  * (7) D. S. Lee and A. J. Steckl “Enhanced blue and green emission in rare-earth-doped GaN electroluminescent devices by optical photopumping” Appl. Phys. Lett. 81, 2331(2002);
  * (8) T. Kita, et. al. “Narrow-band deep-ultraviolet light emitting device using Al1-xGdxN” Appl. Phys. Lett. 93, 211901 (2008)
  * (9) J. H. Park. and A. J. Steckl “Laser action in Eu-doped GaN thin-film cavity at room temperature” Appl. Phys. Lett. 85, 4588 (2004)
  * (10) J. M . Zavada et. al., “Ultraviolet photoluminescence from Gd-implanted AlN epilayers” Appl. Phys. Lett. 89, 152107 (2006)
  * (11) U. Vetter, J. Zenneck and Hofsass “Intense ultraviolet cathodoluminescence at 318nm from Gd3+-doped AlN” Appl. Phys. Lett. 83, 2145 (2003)
  * (12) J. B. Gruber, U. Vetter, H. Hofsass, B. Zandi and M. F. Reid “Spectra and energy levels of Gd3+ (4f7) in AlN” Phys. Rev. B 69, 195202 (2004)
  * (13) M.A. Scarpulla, C. S. Gallinat, S. Mack, J. S. Speck, A. C. Gossard, J. Crys. Grow. 311, 1239 (2009)
  * (14) S. Kitayama, et. al. “Influence of local atomic configuration in AlGdN phosphor thin films on deep ultra-violet luminescence intensity” J. Appl. Phys. 110, 093108 (2011)
  * (15) M. Godlewski, M. Leskel’́a, CRC Critical Reviews in Solid State and Materials Sciences 19 (1994) 199 “Excitation and recombination processes during elecroluminescence of rare earth-activated materials”
  * (16) L. Bodiou, A. Braud “Direc evidence of trap-mediated excitation in GaN:Er3+ with a two-color experiment” Appl. Phys. Lett 93, 151107 (2008)
  * (17) S. D. Carnevale , J. Yang , Patrick J. Phillips , M. J. Mills, and R. C. Myers “Three-Dimensional GaN/AlN Nanowire Heterostructures by Separating Nucleation and Growth Processes” Nano. Lett.,2011, 11(2), pp 866-867
  * (18) S. D. Carnevale, T. F. Kent, P. J. Phillips, M. J. Mills, S. Rajan, and R. C. Myers “Polarization-Induced pn Diodes in Wide-Band-Gap Nanowires with Ultraviolet Electroluminescence” Nano. Lett.,2012, 12(2), pp 915-920
  * (19) S. D. Carnevale, C. Marginean, P. J. Phillips, T. F. Kent, A. T. M. G. Sarwar, M. J. Mills, and R. C. Myers “Coaxial nanowire resonant tunneling diodes from non-polar AlN/GaN on silicon” Appl. Phys. Lett. 100, 142115 (2012)
  * (20) M. Grundman (2005) BandEng Available online: http://my.ece.ucsb.edu/mgrundmann/bandeng.htm
  * (21) J. Simon, V. Protasenko, C. Lian, H. Xing, D. Jena “Polarization-induced hole doping in wide-band-gap uniaxial semiconductor heterostructures” Science 2010, 327, 60-64
  * (22) F. Limbach, C. Hauswalf, J. Lähnemann, M. Wölz, O. Brandt, A. Trampert, M Hanke,et. al., Nanotechnology 23 (2012) 465301
  * (23) E. Nyein, U. Hömmerich, J. Heikenfeld, D. Lee, A. J. Steckl, et. al. Appl. Phys. Lett. 82, 1655 (2003)
  * (24) A. Neuhalfen, B. Wessels Appl. Phys. Lett. 60, 2657 (1992)
  * (25) L. Bergman, M. Dutta, C. Balkas, R. Davis, J. Christman, et. al. J. Appl. Phys. 85, 3535 (1999)
  * (26) 3A. Cingolani, M. Ferrara, M. Lugara, and G. Scamarcio, Solid State Commun. 58, 823 (1986)


# Feedback Design for Devising Optimal Epidemic Control Policies

Muhammad Umar B. Niazi  Philip E. Paré  Karl H. Johansson  Laboratory for
Information and Decision Systems, Massachusetts Institute of Technology, 77
Massachusetts Avenue, Cambridge, MA 02139, USA (Email<EMAIL_ADDRESS>Division of Decision and Control Systems, Digital Futures, School of
Electrical Engineering and Computer Science, KTH Royal Institute of
Technology, Stockholm, Sweden (Email<EMAIL_ADDRESS> Elmore Family School of
Electrical and Computer Engineering, Purdue University, IN, USA (Email:
<EMAIL_ADDRESS>
###### Abstract

This paper proposes a feedback design that effectively copes with
uncertainties for reliable epidemic monitoring and control. There are several
optimization-based methods to estimate the parameters of an epidemic model by
utilizing past reported data. However, due to the possibility of noise in the
data, the estimated parameters may not be accurate, thereby exacerbating the
model uncertainty. To address this issue, we provide an observer design that
enables robust state estimation of epidemic processes, even in the presence of
uncertain models and noisy measurements. Using the estimated model and state,
we then devise optimal control policies by minimizing a predicted cost
functional. To demonstrate the effectiveness of our approach, we implement it
on a modified SIR epidemic model. The results show that our proposed method is
efficient in mitigating the uncertainties that may arise in epidemic
monitoring and control.

###### keywords:

Estimation and control of epidemics, output feedback control, nonlinear
observer design, optimal control.

††thanks: This work is supported the European Union’s Horizon Research and
Innovation Programme under Marie Skłodowska-Curie grant agreement No.
101062523 and the US National Science Foundation, grant NSF-ECCS #2032258.

##  1 Introduction

Since the onset of the SARS-CoV-2 outbreak, there has been a surge of interest
in epidemic processes from many fields including the controls community. These
works typically consider analysis, parameter identification, state estimation,
forecasting, and/or control of a particular compartmental model that may or
may not be networked, e.g., (Hota et al., 2021). Further, there is also a rich
body of literature from the controls field prior to the COVID-19 Pandemic (Wan
et al., 2007; Nowzari et al., 2016; Mei et al., 2017; Paré et al., 2018). In
this work, we present a unified framework for parameter estimation, state
estimation, and optimal control on a generic class of nonlinear models that
includes most of the deterministic epidemiological spreading models in the
literature.

The parameter and state estimation problems are questions of identifiability
and observability, respectively. Conventionally, differential geometric
techniques were employed for obtaining sufficient conditions that verify these
notions for nonlinear systems (Grewal and Glover, 1976; Hermann and Krener,
1977). On the other hand, to obtain necessary and sufficient conditions, Diop
and Fliess (1991) introduced differential algebraic methods for
identifiability and observability. The concepts were further developed and
applied to biological models by Audoly et al. (2001) and Saccomani et al.
(2003).

Once identifiability has been verified, the parameters must be estimated, for
which several methods exist in the literature. Most common among them are the
gradient-based and Newton-type methods like Levenberg–Marquardt and trust
region reflective algorithms; see (Bard, 1974) and (Ljung, 1999). Barz et al.
(2015) present some useful regularization techniques to ensure that the
parameter estimation problem is well-posed. For some basic epidemic models,
like SIR, explicit expressions for parameters are derived by Hadeler (2011)
and Magal and Webb (2018). These ideas have also been explored for networked
epidemics (Paré et al., 2020). However, if the state variables in those
expressions cannot be directly measured, these techniques cannot be employed.

On the other hand, constrained optimal control is also a rich area of research
with two techniques typically employed in practice. The first one is to use
Pontryagin’s minimum principle for computing the optimal control trajectory
(Kirk, 2004, Chapter 5). However, in general, this principle is only a
necessary condition of optimality. It is sufficient only in the case when the
Hamiltonian functional is convex in the state variable (Sethi, 2019, Chapter
2). A practically superior method to solve these types of problems is to
convert them to a constrained nonlinear optimization problem (Betts, 2010) and
then use a numerical solver. Optimal control has been employed for epidemic
mitigation in several works, e.g., (Köhler et al., 2021; Acemoglu et al.,
2021).

Unlike parameter estimation and optimal control, the literature on state
estimation of epidemic processes is lacking. Designing robust observers to
accurately estimate the current state of the epidemic model is the missing
component in the feedback optimal control. However, due to being nonlinear,
observer design for epidemic processes is quite challenging. Extended Kalman
filtering techniques (Rajaei et al., 2021; Gomez-Exposito et al., 2021; Azimi
et al., 2022) are based on linearization and can only provide local
guarantees. Therefore, one must know the initial state quite accurately to
obtain a good state estimate when using these techniques. On the other hand,
observer design techniques for general nonlinear systems with global
guarantees turn out to be very conservative for epidemic models (Niazi and
Johansson, 2022).

Our main contributions include the extension of the observer proposed by Niazi
and Johansson (2022) and provide robust guarantees under model and data
uncertainties. Moreover, devising optimal epidemic policies using state
estimates from the observer is another novelty considered in this paper. Our
results demonstrate that incorporating a robust observer in the feedback loop
yields more reliable epidemic control policies.

The rest of the paper is organized as follows. Section 2 presents the generic
class of nonlinear epidemic models and formulates the problem of interest.
Section 3 outlines our proposed framework and Section 4 provides a necessary
background. Section 5 presents our proposed algorithm for robust state
estimation and Section 6 demonstrates our method on a modified SIR epidemic
model.

Notations. The Euclidean norm of $x\in\mathbb{R}^{n}$ is denoted as
$\|x\|\doteq\sqrt{x^{\top}x}$. For a function $w\in
L^{\infty}(\mathbb{R};\mathbb{R}^{n})$, the essential supremum norm
$\|w\|_{\infty}\doteq\operatorname*{ess\,sup}_{t\in\mathbb{R}}\|w(t)\|$. By
$w_{[t_{0},t_{1}]}$, we denote the restriction of $w$ to $[t_{0},t_{1}]$ for
some $t_{1}>t_{0}$. The maximum singular value of $M\in\mathbb{R}^{n\times m}$
is denoted as $\sigma_{\max}(M)$. An identity matrix of size $n\times n$ is
$I_{n}$. For $M\in\mathbb{R}^{n\times n}$, $\mathtt{sym}(M)\doteq M+M^{\top}$,
and $M\geq 0$ ($M>0$) means that $M$ is positive semi-definite (resp.,
definite).

##  2 Problem Definition

We consider a class of deterministic epidemic models that can be described by

$\displaystyle\dot{x}(t)$ $\displaystyle=Ax(t)+Gf(Hx(t),u(t))$ (1a)
$\displaystyle y(t)$ $\displaystyle=Cx(t)$ (1b)

where $x(t)\in\mathcal{X}\subset\mathbb{R}^{n_{x}}$ is the state,
$u(t)\in\mathcal{U}\subset\mathbb{R}^{n_{u}}$ is the control input, and
$y(t)\in\mathbb{R}^{n_{y}}$ is the measured output. Each element $x_{i}(t)$ of
the state vector corresponds to a different epidemic variable or compartment,
and each element $y_{i}(t)$ of the output vector corresponds to a certain
measurable epidemic variable. On the other hand, each element $u_{i}(t)$ of
the input vector corresponds to a certain pharmaceutical or non-pharmaceutical
interventions like improving medical facilities, testing and isolation,
quarantining, vaccination, lockdown, social distancing, and travel
restrictions, which can be enforced by a public authority.

We have $A\in\mathbb{R}^{n_{x}\times n_{x}},G\in\mathbb{R}^{n_{x}\times
n_{g}},C\in\mathbb{R}^{n_{y}\times n_{x}}$ with

$\displaystyle A$ $\displaystyle\doteq A(\theta),\quad G\doteq G(\theta),\quad
C\doteq C(\theta)$

where $\theta$ is the vector of epidemic parameters. The matrix
$H\in\\{0,1\\}^{n_{H}\times n_{x}}$ is known and specifies the state variables
involved in the nonlinear function
$f:\mathcal{X}\times\mathcal{U}\rightarrow\mathbb{R}^{n_{f}}$, where $f$ is
smooth and thus Lipschitz continuous on a compact domain
$\mathcal{X}\times\mathcal{U}$. That is, for every $x,\hat{x}\in\mathcal{X}$
and $u\in\mathcal{U}$, there exists $\ell\geq 0$ such that

$\|f(Hx,u)-f(H\hat{x},u)\|\leq\ell\|Hx-H\hat{x}\|$ (2)

where

$\ell=\sup_{(x,u)\in\mathcal{X}\times\mathcal{U}}\sigma_{\max}\left(\frac{\partial
f}{\partial x}(Hx,u)\right).$ (3)

Note that $f$ depends only on the state $x(t)$ and the input $u(t)$, and not
on the parameters $\theta$.

###### Remark 1

The class of nonlinear systems (1) captures a variety of epidemic models in
the literature. For instance, all the basic SIS, SIR, SEIR models (Hethcote,
1989, 1994, 2000) and their variants (Arino and Van den Driessche, 2003;
Giordano et al., 2020, 2021; Niazi et al., 2021) can be written in the form of
(1a). The networked epidemic models (Nowzari et al., 2016; Mei et al., 2017;
Paré et al., 2020) can also written compactly as (1a). $\triangle$

The reported data on a time interval $[t_{0},t_{1}]$ is given by

$\displaystyle\bar{u}(t)$ $\displaystyle=u(t)+\delta_{u}(t)$ (4a)
$\displaystyle\bar{y}(t)$ $\displaystyle=y(t)+\delta_{y}(t)$ (4b)

where $t_{0}\geq 0$ is the time of epidemic onset, $t_{1}>t_{0}$ is the
current time, and $\delta_{u}(t)$ and $\delta_{y}(t)$ represent the
uncertainties in the input-output data. The uncertainties $\delta_{u}(t)$ and
$\delta_{y}(t)$ are unknown and account for clerical errors and delays in
recording and reporting the data.

Problem statement. Given the input-output data $(\bar{u},\bar{y})$ for the
past time interval $[t_{0},t_{1}]$, we first aim to estimate the parameters
$\theta$ and the current state $x(t_{1})$ of (1). Then, based on the estimated
model, we devise an optimal control policy $u(t)$ for a future time interval
$[t_{1},t_{2}]$, $t_{2}>t_{1}$, by minimizing a given cost functional

$J(x_{[t_{1},t_{2}]},u_{[t_{1},t_{2}]})\doteq\int_{t_{1}}^{t_{2}}q(x,u,t)\text{d}t$
(5)

subject to a set of specified constraints

$\begin{array}[]{rclrcl}r_{i}(x_{[t_{1},t_{2}]},u_{[t_{1},t_{2}]})&=&0,&\quad
i&=&1,2,\dots,k\\\ s_{j}(x_{[t_{1},t_{2}]},u_{[t_{1},t_{2}]})&\leq&0,&\quad
j&=&1,2,\dots,l\end{array}$ (6)

where $x_{[t_{1},t_{2}]}\doteq x(t,u;\hat{x}_{t_{1}},\hat{\theta})$ is the
predicted state trajectory obtained by integrating (1a) for
$t\in[t_{1},t_{2}]$ using the values of estimated parameters
$\theta=\hat{\theta}$ and choosing the initial condition as the estimated
state $x(t_{1})=\hat{x}_{t_{1}}$.

##  3 Outline of the Proposed Method

Given the past input-output data (4) and the cost functional (5) with
constraints (6), the proposed feedback design for (1) has three main
constituents.

  1. 1.

Parameter estimation. Given the past input-output data
$(\bar{u}(t),\bar{y}(t))$, for $t\in[t_{0},t_{1}]$, we estimate the model
parameters $\theta$ by solving

$\hat{\theta}=\arg\min_{\theta\in\Theta}\int_{t_{0}}^{t_{1}}\|\bar{y}(t)-y(t,\bar{u};\theta)\|\text{d}t$
(7)

where $y(t,\bar{u};\theta)$ is the output of (1) at time $t$ given the
parameters $\theta$.

  2. 2.

State estimation. Given the past input-output data $(\bar{u}(t),\bar{y}(t))$,
for $t\in[t_{0},t_{1}]$, we design a state observer that estimates the current
state $x(t_{1})$ by solving

$\hat{x}(t_{1})=\arg\min_{x_{t_{1}}\in\mathcal{X}}\int_{t_{0}}^{t_{1}}\|\bar{y}(t)-y(t,\bar{u};x_{t_{1}},\hat{\theta})\|\text{d}t$
(8)

where $y(t,\bar{u};x_{t_{1}},\hat{\theta})$ is the output of (1) at time $t$
when the state trajectory goes through $x_{t_{1}}$ at time $t_{1}$ and the
parameters $\theta=\hat{\theta}$.

  3. 3.

Optimal control: Given the estimated model (1) with parameters $\hat{\theta}$
and the state $\hat{x}(t_{1})$, we obtain optimal control policies by solving

$\displaystyle u^{*}(t)=$ $\displaystyle\arg\min_{u\in\mathcal{U}}J(x,u)$ (9a)
subject to (6) and
$\displaystyle\left\\{\begin{array}[]{ccl}\dot{x}&=&A(\hat{\theta})x+G(\hat{\theta})f(Hx,u)\\\\[5.0pt]
t&\in&[t_{1},t_{2}],\;x(t_{1})=\hat{x}(t_{1})\end{array}\right.$ (9d)

where the cost functional $J(x,u)$ is defined in (5).

Optimal Control Epidemic Process State Estimation Parameter Estimation
$+$$+$$\delta_{y}$$\delta_{u}$$\bar{y}$$\bar{u}$$u$$y$$\hat{\theta}$$\hat{x}$$\hat{\theta}$
Figure 1: Block scheme of optimal feedback control.

##  4 Background Material

In this section, we introduce notions of identifiability and observability,
and describe techniques to solve the parameter estimation and optimal control
problems.

###  4.1 Verifying identifiability and observability

Identifiability is necessary for estimating model parameters from the input-
output data $(\bar{u},\bar{y})$. If the model has non-identifiable parameters,
then it is not possible to estimate them uniquely (Saccomani et al., 2003).

###### Definition 2

System (1) is locally identifiable if, for almost any parameter vector
$\theta\in\Theta$, there exists a neighborhood
$\mathcal{N}(\theta)\subseteq\Theta$ such that, for every
$\hat{\theta}\in\mathcal{N}(\theta)$,

$y(t;\theta)=y(t;\hat{\theta}),\forall t\geq
0\Leftrightarrow\theta=\hat{\theta}.$

On the other hand, the notion of observability guarantees whether or not
input-output data contains sufficient information about the system’s state
(Bernard et al., 2022).

###### Definition 3

System (1) is locally observable if, for any $\tau\geq 0$ and almost any
$x(\tau)\in\mathcal{X}$, there exists a neighborhood
$\mathcal{N}(x(\tau))\in\mathcal{X}$ such that, for every
$\hat{x}(\tau)\in\mathcal{N}(x(\tau))$,

$y(t,u;x(\tau))=y(t,u;\hat{x}(\tau)),\forall t\in[\tau,\infty)\Leftrightarrow
x(\tau)=\hat{x}(\tau).$

Identifiability and observability of (1) are required for the well-posedness
of the parameter and state estimation problems, respectively. To verify these
notions for (1), we employ the GenSSI algorithm proposed by Ligon et al.
(2018), which tests the injectivity of the observation map obtained by taking
Lie derivatives of the output.

###  4.2 Solving the parameter estimation problem

Once identifiability of the system has been verified, the parameters must be
estimated. To carry out this process, optimization algorithms like Levenberg-
Marquardt (Moré, 1978) and trust region (Byrd et al., 1987) can be used.

###  4.3 Solving the optimal control problem

In order to solve (9), we convert it to a constrained nonlinear optimization
problem (Betts, 2010) and use interior point (Byrd et al., 1999) or trust
region reflective (Coleman and Li, 1994) methods.

##  5 Algorithm for Robust State Estimation

The main challenge in the proposed scheme (Figure 1) is the state estimation
problem. Existing observer design techniques for the state estimation of
epidemic processes are quite conservative and often turn out to be infeasible
(Niazi and Johansson, 2022). Moreover, observers for epidemic models are often
designed for specific compartmental models and cannot be adapted to other
models. Here, we extend an observer for general epidemic models proposed by
Niazi and Johansson (2022) to include model uncertainties and measurement
noise.

After estimating the parameters from the data (4), (1) can be written as an
uncertain nonlinear system

$\displaystyle\dot{x}(t)$
$\displaystyle=\hat{A}x(t)+\hat{G}f(Hx(t),\bar{u}(t))+w(t)$ (10a)
$\displaystyle\bar{y}(t)$ $\displaystyle=\hat{C}x(t)+v(t)$ (10b)

where $w(t)\in\mathbb{R}^{n_{x}}$ is the model uncertainty,
$v(t)\in\mathbb{R}^{n_{y}}$ is the measurement noise, and

$\hat{A}\doteq A(\hat{\theta}),\quad\hat{G}\doteq
G(\hat{\theta}),\quad\hat{C}\doteq C(\hat{\theta}).$

Notice that the model uncertainty and measurement noise result from the
uncertainties in the input-output data and the parameter estimation error
$\theta-\hat{\theta}$.

Consider the observer proposed by Niazi and Johansson (2022):

$\displaystyle\dot{z}(t)$
$\displaystyle=Mz(t)+(ML+J)\bar{y}(t)+N\hat{G}f(q(t),\bar{u}(t))$ (11a)
$\displaystyle\hat{x}(t)$ $\displaystyle=z(t)+L\bar{y}(t)$ (11b)
$\displaystyle\hat{y}(t)$ $\displaystyle=\hat{C}\hat{x}(t)$ (11c)

where $q(t)\doteq H\hat{x}(t)+K(\bar{y}(t)-\hat{y}(t))$,
$J,L\in\mathbb{R}^{n_{x}\times n_{y}}$ and $K\in\mathbb{R}^{n_{H}\times
n_{y}}$ are matrices to be designed, and

$\displaystyle M=\hat{A}-L\hat{C}\hat{A}-J\hat{C},\quad N=I_{n_{x}}-L\hat{C}.$

Here, $z(t)\in\mathbb{R}^{n_{x}}$ is the observer’s state, and
$\hat{x}(t)\in\mathbb{R}^{n_{x}}$ and $\hat{y}(t)\in\mathbb{R}^{n_{y}}$ are
the state and output estimate of (10).

Consider the semidefinite programming (SDP) problem:

$\displaystyle\text{minimize}~{}\mu~{}\text{subject to}$ (12a)
$\displaystyle\left[\begin{array}[]{cc}\mathtt{sym}(P\hat{A}-R\hat{C}\hat{A}-S\hat{C})+Q&(P-R\hat{C})\hat{G}\\\
G^{\top}(P-R\hat{C})^{\top}&-I_{n_{f}}\end{array}\right]<0$ (12d)
$\displaystyle\left[\begin{array}[]{cc}-Q&(H-K\hat{C})^{\top}\\\
H-K\hat{C}&-\frac{1}{\ell^{2}}I_{n_{H}}\end{array}\right]\leq 0$ (12g)
$\displaystyle\left[\begin{array}[]{cc}-\mu I_{n_{x}}&R\\\
R^{\top}&-I_{n_{y}}\end{array}\right]\leq 0$ (12j) $\displaystyle
P=P^{\top}>0~{}\text{and}~{}Q=Q^{\top}>0$ (12k)

where $\ell$ is the Lipschitz constant obtained from (3), and

$J=P^{-1}S,\quad L=P^{-1}R.$

Now, we show that feasibility of (12) is sufficient for the existence of a
robust observer (11).

###### Theorem 4

If the SDP problem (12) is feasible, then there exist a $\mathcal{KL}$
function $\beta$ and $\mathcal{K}_{\infty}$ functions $\alpha_{1},\alpha_{2}$
such that the estimation error satisfies

$\|x(t_{1})-\hat{x}(t_{1})\|\leq\beta(\|x(t_{0})-\hat{x}(t_{0})\|,t_{1})\\\
+\alpha_{1}(\|w_{[t_{0},t_{1}]}\|_{\infty})+\alpha_{2}(\|v_{[t_{0},t_{1}]}\|_{\infty}).$

[Sketch] Define the estimation error

$e\doteq x-\hat{x}=x-(z+L\bar{y})=\eta-Lv$ (13)

where $x,\bar{y}$ are given in (10), $z,\hat{x}$ are given in (11), and
$\eta\doteq Nx-z.$ Then, by taking the time derivative of $\eta$, we obtain an
error system

$\displaystyle\dot{\eta}$ $\displaystyle=M\eta+NG\tilde{f}+Nw-(ML+J)v,\quad
e=\eta-Lv$

where $\tilde{f}\doteq
f(Hx,\bar{u})-f(H\hat{x}+K(\bar{y}-\hat{C}\hat{x}),\bar{u})$ and the
dependence on $t$ is omitted for brevity. Here, $\eta(t)$ is the state, $e(t)$
is the output, and $w(t),v(t)$ are unknown inputs. Finally, subject to the
feasibility of (12), the proof is completed by showing input-to-output
stability (Sontag and Wang, 2000) of the above system. ∎

The above theorem states that the state estimation error (13) is stable with
respect to the data and parameter uncertainties. Particularly, the noise
attenuation is directly related to the parameter $\mu$ in (12). Moreover, in
the absence of uncertainties, the estimate $\hat{x}(t_{1})$ asymptotically
converges to the true $x(t_{1})$ given that $t_{1}$ is sufficiently large. As
a result of the $\mathcal{KL}$ function $\beta$, the transient error resulting
from the poor choice of $\hat{x}(t_{0})$ also converges to zero
asymptotically.

##  6 Application to SIDHER Epidemic Model

In this paper, we demonstrate the proposed method on an SIDHER epidemic model.
After providing the model design, we implement the proposed method step-by-
step and provide the simulation results.

###  6.1 Model design

We consider an SIDHER epidemic model (Susceptible, Infected, Detected,
Hospitalized, Extinct, and Recovered), which is illustrated in Figure 2 and
given by

$\displaystyle\dot{\mathtt{S}}(t)$
$\displaystyle=\lambda\mathtt{R}(t)-\beta\mathtt{S}(t)\mathtt{I}(t)(1-u_{1}(t))-\nu\mathtt{S}(t)u_{4}(t)$
(14a) $\displaystyle\dot{\mathtt{I}}(t)$
$\displaystyle=-\gamma\mathtt{I}(t)+\beta\mathtt{S}(t)\mathtt{I}(t)(1-u_{1}(t))-\tau\mathtt{I}(t)u_{3}(t)$
(14b) $\displaystyle\dot{\mathtt{D}}(t)$
$\displaystyle=-(\rho+\phi)\mathtt{D}(t)+\tau\mathtt{I}(t)u_{3}(t)$ (14c)
$\displaystyle\dot{\mathtt{H}}(t)$
$\displaystyle=-\xi\mathtt{H}(t)(1-u_{2}(t))+\phi\mathtt{D}(t)-\sigma\mathtt{H}(t)u_{2}(t)$
(14d) $\displaystyle\dot{\mathtt{E}}(t)$
$\displaystyle=\xi\mathtt{H}(t)(1-u_{2}(t))$ (14e)
$\displaystyle\dot{\mathtt{R}}(t)$
$\displaystyle=-\lambda\mathtt{R}(t)+\gamma\mathtt{I}(t)+\rho\mathtt{D}(t)+\nu\mathtt{S}(t)u_{4}(t)+\sigma\mathtt{H}(t)u_{2}(t)$
(14f)

where all the state variables
$\mathtt{S}(t),\mathtt{I}(t),\mathtt{D}(t),\mathtt{H}(t),\mathtt{E}(t),\mathtt{R}(t)\in[0,1]$
and control input
$u(t)=[\begin{array}[]{cccc}u_{1}(t)&u_{2}(t)&u_{3}(t)&u_{4}(t)\end{array}]^{\top}$.
Note that, for every $t\geq 0$,

$\mathtt{S}(t)+\mathtt{I}(t)+\mathtt{D}(t)+\mathtt{H}(t)+\mathtt{E}(t)+\mathtt{R}(t)=1.$
(15)

$\mathtt{S}$ $\mathtt{I}$ $\mathtt{D}$ $\mathtt{H}$ $\mathtt{E}$ $\mathtt{R}$ $\beta$$\nu$$\tau$$\gamma$$\phi$$\rho$$\sigma$$\xi$$\lambda$ Figure 2: Block diagram of SIDHER epidemic model. Table 1: Description of control inputs, which are unitless. Control | Description
---|---
$u_{1}$ | Stringency of NPIs
$u_{2}$ | Proportion of medical resources dedicated
$u_{3}$ | Testing capacity per population
$u_{4}$ | Vaccination capacity per population
Table 2: Description of parameters. Parameter | Description | Unit
---|---|---
$\beta$ | Infection rate | 1/day
$\gamma$ | Recovery rate of undetected cases | 1/day
$\rho$ | Recovery rate of detected cases | 1/day
$\sigma$ | Recovery rate of hospitalized cases | 1/day
$\xi$ | Mortality rate of hospitalized cases | 1/day
$\lambda$ | Rate at which people lose immunity | 1/day
$\phi$ | Hospitalization rate | 1/day
$\tau$ | Testing rate of undetected cases | 1/day
$\nu$ | Vaccination rate of susceptible cases | 1/day

Control inputs. The control inputs, summarized in Table 1, are described
below:

  * •

$u_{1}(t)\in[0,1]$, which is the stringency of the non-pharmaceutical
interventions (NPIs) like lockdown, social distancing, travel restriction,
etc.

  * •

$u_{2}(t)\in[0,1]$, which is a measure of resources invested in improving
medical resources and infrastructure, implementing better treatment methods of
hospitalized cases, and medical facilities dedicated to the epidemic.

  * •

$u_{3}(t)\in[0,1]$, which is the testing capacity in terms of the fraction of
population that can be tested.

  * •

$u_{4}(t)\in[0,1]$, which is the vaccination capacity in terms of the fraction
of population that can be vaccinated.

Parameters. Summarized in Table 2, the parameters are described below:

  * •

The infection rate $\beta\geq 0$ is the product of contact rate (average
number of contacts each person makes per day) and the infection probability
(probability that a susceptible person gets infected after coming in contact
with an infected person).

  * •

The recovery rates $\gamma,\rho,\sigma\in[0,1]$ are the inverses of average
infection periods of infected, detected, and hospitalized cases, respectively.
The mortality rate $\xi\in[0,1]$ is the inverse of average number of days
after which a typical non-surviving hospitalized case dies. The rates
$\lambda,\phi\in[0,1]$ are the inverses of the average number of days after
which recovered or vaccinated cases lose their immunity and detected cases are
hospitalized, respectively.

  * •

The testing and vaccination rates $\tau,\nu\in[0,1]$ are the fractions of
infected individuals that get tested and susceptible individuals that get
vaccinated every day on average, respectively.

Measured outputs. The model outputs are the following:

  * •

$y_{1}=\nu\mathtt{S}$: Since we measure the proportion of population
vaccinated per day $\nu\mathtt{S}u_{4}$ and the vaccination capacity $u_{4}$
is known, we obtain the output $y_{1}$ by dividing both.

  * •

$y_{2}=\tau\mathtt{I}$: Since we measure the proportion of population tested
per day $\tau\mathtt{I}u_{3}$ and the testing capacity $u_{3}$ is known, we
obtain the output $y_{2}$ by dividing both.

  * •

$y_{3}=\mathtt{D}$: Active number of detected infected cases.

  * •

$y_{4}=\rho\mathtt{D}$: Daily number of cases recovering after being detected.
This may not be directly measured if people are not tested again after their
infectious period. In such a case, this output variable can be inferred from
the daily reported data on active detected cases $y_{3}$, active hospitalized
cases $y_{6}$, and daily number of people recovering from hospitals $y_{7}$
(Niazi et al., 2021).

  * •

$y_{5}=\phi\mathtt{D}$: Daily number of cases hospitalized after being
detected.

  * •

$y_{6}=\mathtt{H}$: Active number of hospitalized cases.

  * •

$y_{7}=\sigma\mathtt{H}$: Daily number of cases recovering after
hospitalization.

  * •

$y_{8}=\xi\mathtt{H}$: Daily number of deaths. We assume that non-surviving
cases are first detected and then hospitalized as their symptoms worsen.

  * •

$y_{9}=\mathtt{E}$: Total number of deaths.

  * •

$y_{10}=\mathtt{S}+\mathtt{I}+\mathtt{R}$: Since (15) holds and
$\mathtt{D},\mathtt{H},\mathtt{E}$ are measured, we can obtain the output
$y_{10}$ by subtracting $\mathtt{D}+\mathtt{H}+\mathtt{E}$ from $1$.

Model in vector form. We can write the model (14) in the form (1) with
$f(Hx,u)=[\begin{array}[]{ccccc}\mathtt{S}\mathtt{I}&\mathtt{S}\mathtt{I}u_{1}&\mathtt{H}u_{2}&\mathtt{I}u_{3}&\mathtt{S}u_{4}\end{array}]^{\top}$
and

$\displaystyle A$
$\displaystyle={\scriptsize\left[\begin{array}[]{rrcrrr}0&0&0&0&0&\lambda\\\
0&-\gamma&0&0&0&0\\\ 0&0&-(\rho+\phi)&0&0&0\\\ 0&0&\phi&-\xi&0&0\\\
0&0&0&\xi&0&0\\\ 0&\gamma&\rho&0&0&-\lambda\end{array}\right]},\quad
H={\scriptsize\left[\begin{array}[]{cccccc}1&0&0&0&0&0\\\ 0&1&0&0&0&0\\\
0&0&0&1&0&0\end{array}\right]},$ $\displaystyle G$
$\displaystyle={\scriptsize\left[\begin{array}[]{rrrrr}-\beta&\beta&0&0&-\nu\\\
\beta&-\beta&0&-\tau&0\\\ 0&0&0&\tau&0\\\ 0&0&-\sigma+\xi&0&0\\\
0&0&-\xi&0&0\\\ 0&0&\sigma&0&\nu\end{array}\right]},\quad
C={\scriptsize\left[\begin{array}[]{cccccc}\nu&0&0&0&0&0\\\ 0&\tau&0&0&0&0\\\
0&0&1&0&0&0\\\ 0&0&\rho&0&0&0\\\ 0&0&\phi&0&0&0\\\ 0&0&0&1&0&0\\\
0&0&0&\sigma&0&0\\\ 0&0&0&\xi&0&0\\\ 0&0&0&0&1&0\\\
1&1&0&0&0&1\end{array}\right]}.$

Note that $(A,C)$ is an observable pair.

Identifiability and observability. Identifiability and observability of (14)
is verified by the GenSSI software (Ligon et al., 2018) in MATLAB. Fig. 3
illustrates the resulting observability and identifiability tableau, which is
the zero-pattern structure of the Jacobian of the observability and
identifiability map. The black boxes in the figure represent non-zero terms of
the Jacobian, whereas white area represents zero terms. From the figure, it
can be seen that the Jacobian has full generic rank, thus implying structural
identifiability and observability. Therefore, it can be inferred that for
almost all values of the initial states and parameters, (14) is at least
locally identifiable and observable.

Figure 3: Observability and identifiability tableau.

###  6.2 Simulation results

Data generation. We generate noisy synthetic data by simulating the model (14)
for $t\in[t_{0},t_{1}]$, where $t_{0}=0$ and $t_{1}=30$ days, with “true”
parameters provided in Table 3 to illustrate the proposed method. The input-
output data $(\bar{u},\bar{y})$ is corrupted with white Gaussian noise that is
sampled from $\mathcal{N}(0,10^{-6})$. The nominal input is chosen as

$u(t)=\left[\begin{array}[]{c}0.01\lceil\sin(t/2)\rfloor+0.015\\\
0.01\lceil\cos(t/2)\rfloor+0.015\\\ 0.01\lceil\sin(t/3)\rfloor+0.015\\\
0.01\lceil\cos(t/3)\rfloor+0.015\end{array}\right].$

The true initial state is chosen to be

$x_{t_{0}}=[\begin{array}[]{cccccc}0.999&0.0005&0.0005&0&0&0\end{array}]^{\top}$

which is only used for data generation and is not known by the parameter and
state estimation algorithms.

Parameter estimation. The parameters $\rho,\phi,\sigma,\xi$ can be estimated
using the least square solution by

$\displaystyle\hat{\rho}$
$\displaystyle=\int_{t_{0}}^{t_{1}}\frac{\bar{y}_{3}(t)^{\top}\bar{y}_{4}(t)}{\bar{y}_{3}(t)^{\top}\bar{y}_{3}(t)}\text{d}t,$
$\displaystyle\hat{\phi}$
$\displaystyle=\int_{t_{0}}^{t_{1}}\frac{\bar{y}_{3}(t)^{\top}\bar{y}_{5}(t)}{\bar{y}_{3}(t)^{\top}\bar{y}_{3}(t)}\text{d}t$
$\displaystyle\hat{\sigma}$
$\displaystyle=\int_{t_{0}}^{t_{1}}\frac{\bar{y}_{6}(t)^{\top}\bar{y}_{7}(t)}{\bar{y}_{6}(t)^{\top}\bar{y}_{6}(t)}\text{d}t,$
$\displaystyle\hat{\xi}$
$\displaystyle=\int_{t_{0}}^{t_{1}}\frac{\bar{y}_{6}(t)^{\top}\bar{y}_{8}(t)}{\bar{y}_{6}(t)^{\top}\bar{y}_{6}(t)}\text{d}t.$

The data $\bar{y}$ is discrete and can be interpolated for computing the above
integrals. By fixing these estimated parameters, the remaining parameters are
estimated using the nlgreyest function from the System Identification Toolbox
in MATLAB employing the Trust-Region-Reflective Algorithm. See Table 3 for the
estimated parameters. Notice that, for parameter estimation, we do not know
the true initial state $x_{t_{0}}$. Instead, we guess the initial state
appropriately, where $\mathtt{S}_{t_{0}}$ is chosen uniformly at random from
$[0.95,1]$; $\mathtt{I}_{t_{0}}$ from $[0,0.05]$;
$\mathtt{D}_{t_{0}},\mathtt{H}_{t_{0}},\mathtt{E}_{t_{0}}$ are obtained from
$y_{3}(t_{0}),y_{6}(t_{0}),y_{9}(t_{0})$; and $\mathtt{R}_{t_{0}}$ is obtained
from (15).

Table 3: True and estimated values of the model parameters. Parameter | True value | Estimated value
---|---|---
$\beta$ | 0.3500 | 0.3530
$\gamma$ | 0.1000 | 0.0981
$\rho$ | 0.0500 | 0.0501
$\sigma$ | 0.0400 | 0.0399
$\xi$ | 0.0200 | 0.0202
$\lambda$ | 0.0167 | 0.0383
$\phi$ | 0.1429 | 0.1428
$\tau$ | 0.3000 | 0.2757
$\nu$ | 0.0100 | 0.0100

###### Remark 5

Since the system (14) is locally identifiable, and the parameter estimation
problem (7) is non-convex, the solution obtained in Table 3 might be local.
That is, the solution depends on the appropriate initialization of these
parameters for the optimization algorithm. To address this issue, one can
adopt global optimization techniques to ensure that the true parameters can be
efficiently estimated. Moreover, the observer can be extended to
simultaneously estimate the state and parameters of the model. However, this
is an interesting prospect that will be addressed in the future. $\triangle$

Figure 4: State estimation by the proposed observer. Figure 5: Online output
tracking by the observer.

State estimation. After obtaining the estimated parameters, we employ the
observer (11) for estimating the state of (14). The observer is designed by
solving (12), which yields the design matrices $J,K,L$111Due to their large
size, it is not possible to provide these matrices in the paper. Kindly
contact the authors for the code and data.. The result of the state estimation
method is demonstrated in Fig. 4. It can be seen that the observer converges
very closely to the true state within three days of the epidemic outbreak.
Secondly, Fig. 5 shows the output prediction/tracking by the observer at every
time instant. Notice that, even though it was initialized arbitrarily, the
observer converges very fast and minimizes the difference between the measured
output and the predicted output.

Optimal control. In the optimal control problem, we choose the cost functional
(5) as

$J(x,u)=\int_{t_{1}}^{t_{2}}\left(x(t)^{\top}\Gamma x(t)+u(t)^{\top}\Lambda
u(t)\right)\text{d}t$

where $\Gamma\geq 0$ and $\Lambda>0$ are chosen to be

$\displaystyle\Gamma$ $\displaystyle=\mathtt{diag}(0.01,1,0,2,10,0)$
$\displaystyle\Lambda$ $\displaystyle=\mathtt{diag}(0.01,0.01,0.01,0.01).$

In short, we would like to minimize the susceptible, infected, hospitalized,
and extinct cases in the time interval $[t_{1},t_{2}]$. The control inputs are
bounded as $0\leq u_{1}(t)\leq 1$, $0\leq u_{2}(t)\leq 0.9$, $0.1\leq
u_{3}(t)\leq 0.7$, and $0\leq u_{4}(t)\leq 0.7$. Notice that, in general, it
may be physically impossible to make $u_{2}(t)=1$ because that would mean that
inflow to $\mathtt{E}$ is zero, i.e., no one dies due to the epidemic.
Moreover, we have $u_{3}(t)>0$ to allow a baseline diagnosis of the infected
cases. If $u_{3}(t)=0$, then the inflow to $\mathtt{D}$ is zero, i.e., no one
is diagnosed with the disease.

The equality constraint in (6) is (14) with estimated parameters (Table 3) and
initialized at estimated state $\hat{x}(t_{1})$. The inequality constraints in
(6) are

$\mathtt{I}(t)-\bar{\mathtt{I}}\leq 0,\;\mathtt{H}(t)-\bar{\mathtt{H}}\leq
0,\;\mathtt{E}(t)-\bar{\mathtt{E}}\leq 0$

where we choose $\bar{I}=0.5$, $\bar{\mathtt{H}}=0.05$, and
$\bar{\mathtt{E}}=0.005$.

Figure 6: State trajectory with optimal control input. Figure 7: Sequence of
optimal control input.

We use the fmincon solver in MATLAB for obtaining the optimal solution of (9).
We consider a piecewise constant control input $u(t)$, where $u(t)$ remains
constant for a period of 14 days. The optimal control sequence is obtained for
five such periods, i.e., the solution recommends how the policy should be
varied in the future after every 14 days. The obtained optimal control input
is shown in Fig. 7 and the resulting predicted state trajectory in Fig. 6.
Notice from Fig. 6 that the constraints on $\mathtt{I}$, $\mathtt{H}$, and
$\mathtt{E}$ are satisfied. From the forecast of $\mathtt{E}$, using the
polynomial fit of the estimated trajectory, it can be seen that the number of
deaths are significantly reduced under the optimal control algorithm. This
prediction is reliable only when, in addition to the parameters, the current
state $x(t_{1})$ is accurately estimated. The results may vary significantly
if the unmeasured states are not chosen accurately. From Fig. 7, we can see
that the optimal controller recommends to enact a full lockdown for the first
four weeks to suppress the epidemic growth, i.e., to make
$\dot{\mathtt{I}}<0$. The lockdown is then lifted gradually by increasing the
testing and vaccination capacities. On the other hand, to avoid the number of
deaths and to satisfy the hard constraint $\mathtt{E}(t)\leq\bar{\mathtt{E}}$,
the optimal controller recommends to employ full medical resources throughout
the finite future horizon.

##  7 Concluding Remarks

We presented a unified framework for feedback optimal control of epidemic
processes via a general class of nonlinear compartmental models. For the
parameter estimation and constrained optimal control, we employed existing
methods and techniques from system identification and optimal control theory.
For the state estimation, however, the proposed robust observer design
criteria is a novel contribution of this paper. By considering a realistic
epidemic model, we demonstrated how optimal control policies can be devised by
employing the proposed feedback design.

The method proposed in this paper could be helpful in devising optimal control
policies during an epidemic outbreak. Robustness and fast convergence of the
observer guarantees accurate state estimation of an epidemic model even when
few samples of the data are collected. Extending the current method to a
moving horizon estimation and control is a straightforward task and will be
addressed in the future. Another interesting direction is the observer design
for simultaneous parameter and state estimation. Recently proposed learning-
based Luenberger observer design method by Niazi et al. (2022) could be
helpful in this regard. Such an observer is crucial because the existing
parameter estimation techniques require an accurate guess of the initial
state, which may not be possible in real epidemic scenarios.

## References

  * Acemoglu et al. (2021) Acemoglu, D., Fallah, A., Giometto, A., Huttenlocher, D., Ozdaglar, A., Parise, F., and Pattathil, S. (2021).  Optimal adaptive testing for epidemic control: Combining molecular and serology tests.  _arXiv preprint arXiv:2101.00773_.
  * Arino and Van den Driessche (2003) Arino, J. and Van den Driessche, P. (2003).  A multi-city epidemic model.  _Mathematical Population Studies_ , 10(3), 175–193.
  * Audoly et al. (2001) Audoly, S., Bellu, G., D’Angio, L., Saccomani, M.P., and Cobelli, C. (2001).  Global identifiability of nonlinear models of biological systems.  _IEEE Transactions on Biomedical Engineering_ , 48(1), 55–65.
  * Azimi et al. (2022) Azimi, V., Sharifi, M., Fakoorian, S., Nguyen, T., and Van Huynh, V. (2022).  State estimation-based robust optimal control of influenza epidemics in an interactive human society.  _Information Sciences_ , 592, 340–360.
  * Bard (1974) Bard, Y. (1974).  _Nonlinear Parameter Estimation_.  Academic Press, New York, USA.
  * Barz et al. (2015) Barz, T., Körkel, S., Wozny, G., et al. (2015).  Nonlinear ill-posed problem analysis in model-based parameter estimation and experimental design.  _Computers & Chemical Engineering_, 77, 24–42.
  * Bernard et al. (2022) Bernard, P., Andrieu, V., and Astolfi, D. (2022).  Observer design for continuous-time dynamical systems.  _Annual Reviews in Control_.
  * Betts (2010) Betts, J.T. (2010).  _Practical Methods for Optimal Control and Estimation using Nonlinear Programming_.  SIAM.
  * Byrd et al. (1999) Byrd, R.H., Hribar, M.E., and Nocedal, J. (1999).  An interior point algorithm for large-scale nonlinear programming.  _SIAM Journal on Optimization_ , 9(4), 877–900.
  * Byrd et al. (1987) Byrd, R.H., Schnabel, R.B., and Shultz, G.A. (1987).  A trust region algorithm for nonlinearly constrained optimization.  _SIAM Journal on Numerical Analysis_ , 24(5), 1152–1170.
  * Coleman and Li (1994) Coleman, T.F. and Li, Y. (1994).  On the convergence of interior-reflective Newton methods for nonlinear minimization subject to bounds.  _Mathematical Programming_ , 67(1), 189–224.
  * Diop and Fliess (1991) Diop, S. and Fliess, M. (1991).  Nonlinear observability, identifiability, and persistent trajectories.  In _Proc. 30th IEEE Conference on Decision and Control_ , 714–719.
  * Giordano et al. (2020) Giordano, G., Blanchini, F., Bruno, R., Colaneri, P., Di Filippo, A., Di Matteo, A., and Colaneri, M. (2020).  Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy.  _Nature Medicine_ , 26(6), 855–860.
  * Giordano et al. (2021) Giordano, G., Colaneri, M., Di Filippo, A., Blanchini, F., Bolzern, P., De Nicolao, G., Sacchi, P., Colaneri, P., and Bruno, R. (2021).  Modeling vaccination rollouts, SARS-CoV-2 variants and the requirement for non-pharmaceutical interventions in Italy.  _Nature Medicine_ , 27(6), 993–998.
  * Gomez-Exposito et al. (2021) Gomez-Exposito, A., Rosendo-Macias, J.A., and Gonzalez-Cagigal, M.A. (2021).  Monitoring and tracking the evolution of a viral epidemic through nonlinear kalman filtering: Application to the COVID-19 case.  _IEEE Journal of Biomedical and Health Informatics_.
  * Grewal and Glover (1976) Grewal, M. and Glover, K. (1976).  Identifiability of linear and nonlinear dynamical systems.  _IEEE Transactions on Automatic Control_ , 21(6), 833–837.
  * Hadeler (2011) Hadeler, K. (2011).  Parameter identification in epidemic models.  _Mathematical Biosciences_ , 229(2), 185–189.
  * Hermann and Krener (1977) Hermann, R. and Krener, A. (1977).  Nonlinear controllability and observability.  _IEEE Transactions on Automatic Control_ , 22(5), 728–740.
  * Hethcote (1989) Hethcote, H.W. (1989).  Three basic epidemiological models.  In _Applied Mathematical Ecology_ , 119–144. Springer.
  * Hethcote (1994) Hethcote, H.W. (1994).  A thousand and one epidemic models.  In _Frontiers in Mathematical Biology_ , 504–515. Springer.
  * Hethcote (2000) Hethcote, H.W. (2000).  The mathematics of infectious diseases.  _SIAM review_ , 42(4), 599–653.
  * Hota et al. (2021) Hota, A.R., Godbole, J., and Paré, P.E. (2021).  A closed-loop framework for inference, prediction, and control of SIR epidemics on networks.  _IEEE Transactions on Network Science and Engineering_ , 8(3), 2262–2278.
  * Kirk (2004) Kirk, D.E. (2004).  _Optimal Control Theory: An Introduction_.  Courier Corporation.
  * Köhler et al. (2021) Köhler, J., Schwenkel, L., Koch, A., Berberich, J., Pauli, P., and Allgöwer, F. (2021).  Robust and optimal predictive control of the covid-19 outbreak.  _Annual Reviews in Control_ , 51, 525–539.
  * Ligon et al. (2018) Ligon, T.S., Fröhlich, F., Chiş, O.T., Banga, J.R., Balsa-Canto, E., and Hasenauer, J. (2018).  GenSSI 2.0: Multi-experiment structural identifiability analysis of SBML models.  _Bioinformatics_ , 34(8), 1421–1423.
  * Ljung (1999) Ljung, L. (1999).  _System Identification: Theory for the User_.  Prentice-Hall, Upper Saddle River, NJ, 2 edition.
  * Magal and Webb (2018) Magal, P. and Webb, G. (2018).  The parameter identification problem for SIR epidemic models: Identifying unreported cases.  _Journal of Mathematical Biology_ , 77(6), 1629–1648.
  * Mei et al. (2017) Mei, W., Mohagheghi, S., Zampieri, S., and Bullo, F. (2017).  On the dynamics of deterministic epidemic propagation over networks.  _Annual Reviews in Control_ , 44, 116–128.
  * Moré (1978) Moré, J.J. (1978).  The Levenberg-Marquardt algorithm: Implementation and theory.  In _Numerical Analysis_ , 105–116. Springer.
  * Niazi et al. (2022) Niazi, M.U.B., Cao, J., Sun, X., Das, A., and Johansson, K.H. (2022).  Learning-based design of Luenberger observers for autonomous nonlinear systems.  _arXiv preprint arXiv:2210.01476_.
  * Niazi and Johansson (2022) Niazi, M.U.B. and Johansson, K.H. (2022).  Observer design for the state estimation of epidemic processes.  _arXiv preprint arXiv:2207.11977_.
  * Niazi et al. (2021) Niazi, M.U.B., Kibangou, A., Canudas-de Wit, C., Nikitin, D., Tumash, L., and Bliman, P.A. (2021).  Modeling and control of epidemics through testing policies.  _Annual Reviews in Control_ , 52, 554–572.
  * Nowzari et al. (2016) Nowzari, C., Preciado, V.M., and Pappas, G.J. (2016).  Analysis and control of epidemics: A survey of spreading processes on complex networks.  _IEEE Control Systems Magazine_ , 36(1), 26–46.
  * Paré et al. (2020) Paré, P.E., Beck, C.L., and Başar, T. (2020).  Modeling, estimation, and analysis of epidemics over networks: An overview.  _Annual Reviews in Control: Special Issue on Systems and Control Research Efforts Against COVID-19 and Future Pandemics_ , 50, 345–360.
  * Paré et al. (2018) Paré, P.E., Beck, C.L., and Nedić, A. (2018).  Epidemic processes over time-varying networks.  _IEEE Trans. on Control of Network Systems_ , 5(3), 1322–1334.
  * Rajaei et al. (2021) Rajaei, A., Raeiszadeh, M., Azimi, V., and Sharifi, M. (2021).  State estimation-based control of COVID-19 epidemic before and after vaccine development.  _Journal of Process Control_ , 102, 1–14.
  * Saccomani et al. (2003) Saccomani, M.P., Audoly, S., and D’Angiò, L. (2003).  Parameter identifiability of nonlinear systems: The role of initial conditions.  _Automatica_ , 39(4), 619–632.
  * Sethi (2019) Sethi, S.P. (2019).  _Optimal Control Theory: Applications to Management Science and Economics_.  Springer Nature Switzerland AG, 3 edition.
  * Sontag and Wang (2000) Sontag, E. and Wang, Y. (2000).  Lyapunov characterizations of input to output stability.  _SIAM Journal on Control and Optimization_ , 39(1), 226–249.
  * Wan et al. (2007) Wan, Y., Roy, S., and Saberi, A. (2007).  Network design problems for controlling virus spread.  In _Proc. 46th IEEE Conference on Decision and Control_ , 3925–3932.


* Blocki et al. (2012) Jeremiah Blocki, Avrim Blum, Anupam Datta, and Or Sheffet.  The johnson-lindenstrauss transform itself preserves differential privacy.  In _Proceedings of the 53rd Annual IEEE Symposium on Foundations of Computer Science (FOCS)_ , pages 410–419, New Brunswick, NJ, 2012\.
  * Blum et al. (2005) Avrim Blum, Cynthia Dwork, Frank McSherry, and Kobbi Nissim.  Practical privacy: the SuLQ framework.  In _Proceedings of the Twenty-fourth ACM SIGACT-SIGMOD-SIGART Symposium on Principles of Database Systems (PODS)_ , pages 128–138, Baltimore, MD, 2005.
  * Boufounos and Baraniuk (2008) Petros Boufounos and Richard G. Baraniuk.  1-bit compressive sensing.  In _Proceedings of the 42nd Annual Conference on Information Sciences and Systems (CISS)_ , pages 16–21, Princeton, NJ, 2008.
  * Candès et al. (2006) Emmanuel J. Candès, Justin K. Romberg, and Terence Tao.  Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information.  _IEEE Trans. Inf. Theory_ , 52(2):489–509, 2006\.
  * Chang and Lin (2011) Chih-Chung Chang and Chih-Jen Lin.  Libsvm: a library for support vector machines.  _ACM Transactions on Intelligent Systems and Technology (TIST)_ , 2(3):27, 2011.
  * Charikar et al. (2004) Moses Charikar, Kevin Chen, and Martin Farach-Colton.  Finding frequent items in data streams.  _Theor. Comput. Sci._ , 312(1):3–15, 2004.
  * Charikar (2002) Moses S Charikar.  Similarity estimation techniques from rounding algorithms.  In _Proceedings of the Thiry-Fourth Annual ACM Symposium on Theory of Computing (STOC)_ , pages 380–388, Montreal, Canada, 2002.
  * Chaudhuri and Monteleoni (2008) Kamalika Chaudhuri and Claire Monteleoni.  Privacy-preserving logistic regression.  In _Advances in Neural Information Processing Systems (NIPS)_ , pages 289–296, Vancouver, Canada, 2008.
  * Chaudhuri et al. (2011) Kamalika Chaudhuri, Claire Monteleoni, and Anand D. Sarwate.  Differentially private empirical risk minimization.  _J. Mach. Learn. Res._ , 12:1069–1109, 2011.
  * Chen et al. (2017) Danqi Chen, Adam Fisch, Jason Weston, and Antoine Bordes.  Reading wikipedia to answer open-domain questions.  In _Proceedings of the 55th Annual Meeting of the Association for Computational Linguistics (ACL)_ , pages 1870–1879, Vancouver, Canada, 2017\.
  * Chen et al. (2015) Wenlin Chen, James Wilson, Stephen Tyree, Kilian Weinberger, and Yixin Chen.  Compressing Neural Networks with the Hashing Trick.  In _Proceedings of the 32nd International Conference on Machine Learning (ICML)_ , pages 2285–2294, Lille, France, 2015.
  * Cormode et al. (2018a) Graham Cormode, Somesh Jha, Tejas Kulkarni, Ninghui Li, Divesh Srivastava, and Tianhao Wang.  Privacy at scale: Local differential privacy in practice.  In _Proceedings of the 2018 International Conference on Management of Data (SIGMOD)_ , pages 1655–1658, Houston, TX, 2018a.
  * Cormode et al. (2018b) Graham Cormode, Tejas Kulkarni, and Divesh Srivastava.  Marginal release under local differential privacy.  In _Proceedings of the 2018 International Conference on Management of Data (SIGMOD)_ , pages 131–146, Houston, TX, 2018b.
  * Cortes and Vapnik (1995) Corinna Cortes and Vladimir Vapnik.  Support-vector networks.  _Mach. Learn._ , 20(3):273–297, 1995.
  * Dahl et al. (2013) George E. Dahl, Jack W. Stokes, Li Deng, and Dong Yu.  Large-scale malware classification using random projections and neural networks.  In _Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , pages 3422–3426, Vancouver, Canada, 2013.
  * Dasgupta (2000) Sanjoy Dasgupta.  Experiments with random projection.  In _Proceedings of the 16th Conference in Uncertainty in Artificial Intelligence (UAI)_ , pages 143–151, Stanford, CA, 2000.
  * Dasgupta and Freund (2008) Sanjoy Dasgupta and Yoav Freund.  Random projection trees and low dimensional manifolds.  In _Proceedings of the 40th Annual ACM Symposium on Theory of Computing (STOC)_ , pages 537–546, Victoria, Canada, 2008.
  * Datar et al. (2004) Mayur Datar, Nicole Immorlica, Piotr Indyk, and Vahab S Mirrokni.  Locality-sensitive hashing scheme based on p-stable distributions.  In _Proceedings of the Twentieth Annual Symposium on Computational Geometry (SCG)_ , pages 253–262, Brooklyn, NY, 2004.
  * Desfontaines and Pejó (2020) Damien Desfontaines and Balázs Pejó.  Sok: Differential privacies.  _Proc. Priv. Enhancing Technol._ , 2020(2):288–313, 2020.
  * Dickens et al. (2022) Charlie Dickens, Justin Thaler, and Daniel Ting.  Order-invariant cardinality estimators are differentially private.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , New Orleans, LA, 2022.
  * Ding et al. (2017) Bolin Ding, Janardhan Kulkarni, and Sergey Yekhanin.  Collecting telemetry data privately.  In _Advances in Neural Information Processing Systems (NIPS)_ , pages 3571–3580, Long Beach, CA, 2017.
  * Dong et al. (2022) Jinshuo Dong, Aaron Roth, and Weijie J Su.  Gaussian differential privacy.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 84(1):3–37, 2022.
  * Dong et al. (2008) Wei Dong, Moses Charikar, and Kai Li.  Asymmetric distance estimation with sketches for similarity search in high-dimensional spaces.  In _Proceedings of the 31st Annual International ACM SIGIR Conference on Research and Development in Information Retrieval (SIGIR)_ , pages 123–130, 2008.
  * Donoho (2006) David L. Donoho.  Compressed sensing.  _IEEE Trans. Inf. Theory_ , 52(4):1289–1306, 2006.
  * Dwork and Lei (2009) Cynthia Dwork and Jing Lei.  Differential privacy and robust statistics.  In _Proceedings of the 41st Annual ACM Symposium on Theory of Computing (STOC)_ , pages 371–380, Bethesda, MD, 2009.
  * Dwork and Roth (2014) Cynthia Dwork and Aaron Roth.  The algorithmic foundations of differential privacy.  _Found. Trends Theor. Comput. Sci._ , 9(3-4):211–407, 2014.
  * Dwork and Rothblum (2016) Cynthia Dwork and Guy N Rothblum.  Concentrated differential privacy.  _arXiv preprint arXiv:1603.01887_ , 2016.
  * Dwork et al. (2006) Cynthia Dwork, Frank McSherry, Kobbi Nissim, and Adam D. Smith.  Calibrating noise to sensitivity in private data analysis.  In _Proceedings of the Third Theory of Cryptography Conference (TCC)_ , pages 265–284, New York, NY, 2006.
  * Erlingsson et al. (2014) Úlfar Erlingsson, Vasyl Pihur, and Aleksandra Korolova.  RAPPOR: randomized aggregatable privacy-preserving ordinal response.  In _Proceedings of the 2014 ACM SIGSAC Conference on Computer and Communications Security (CCS)_ , pages 1054–1067, Scottsdale, AZ, 2014.
  * Fan and Li (2022) Chenglin Fan and Ping Li.  Distances release with differential privacy in tree and grid graph.  In _IEEE International Symposium on Information Theory (ISIT)_ , pages 2190–2195, 2022.
  * Fan et al. (2022) Chenglin Fan, Ping Li, and Xiaoyun Li.  Private graph all-pairwise-shortest-path distance release with improved error rate.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , New Orleans, LA, 2022.
  * Fang et al. (2023) Huang Fang, Xiaoyun Li, Chenglin Fan, and Ping Li.  Improved convergence of differential private sgd with gradient clipping.  In _Proceedings of the Eleventh International Conference on Learning Representations (ICLR)_ , Kigali, Rwanda, 2023.
  * Feldman et al. (2009) Dan Feldman, Amos Fiat, Haim Kaplan, and Kobbi Nissim.  Private coresets.  In _Proceedings of the 41st Annual ACM Symposium on Theory of Computing (STOC)_ , pages 361–370, Bethesda, MD, 2009.
  * Fern and Brodley (2003) Xiaoli Zhang Fern and Carla E. Brodley.  Random projection for high dimensional data clustering: A cluster ensemble approach.  In _Proceedings of the Twentieth International Conference (ICML)_ , pages 186–193, Washington, DC, 2003.
  * Freund et al. (2007) Yoav Freund, Sanjoy Dasgupta, Mayank Kabra, and Nakul Verma.  Learning the structure of manifolds using random projections.  In _Advances in Neural Information Processing Systems (NIPS)_ , pages 473–480, Vancouver, Canada, 2007.
  * Friedman et al. (1975) Jerome H. Friedman, F. Baskett, and L. Shustek.  An algorithm for finding nearest neighbors.  _IEEE Transactions on Computers_ , 24:1000–1006, 1975.
  * Ge et al. (2018) Jason Ge, Zhaoran Wang, Mengdi Wang, and Han Liu.  Minimax-optimal privacy-preserving sparse PCA in distributed systems.  In _Proceedings of the International Conference on Artificial Intelligence and Statistics (AISTATS)_ , pages 1589–1598, Playa Blanca, Lanzarote, Canary Islands, Spain, 2018.
  * Goemans and Williamson (1995) Michel X. Goemans and David P. Williamson.  Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming.  _J. ACM_ , 42(6):1115–1145, 1995.
  * Gupta et al. (2010) Anupam Gupta, Katrina Ligett, Frank McSherry, Aaron Roth, and Kunal Talwar.  Differentially private combinatorial optimization.  In _Proceedings of the Twenty-First Annual ACM-SIAM Symposium on Discrete Algorithms (SODA)_ , pages 1106–1125, Austin, TX, 2010.
  * Haddadpour et al. (2020) Farzin Haddadpour, Belhal Karimi, Ping Li, and Xiaoyun Li.  Fedsketch: Communication-efficient and private federated learning via sketching.  _arXiv preprint arXiv:2008.04975_ , 2020.
  * Haeberlen and Khanna (2014) Justin Hsu Marco Gaboardi Andreas Haeberlen and Sanjeev Khanna.  Differential privacy: An economic method for choosing epsilon.  _arXiv preprint arXiv:1402.3329_ , 2014.
  * Indyk and Motwani (1998) Piotr Indyk and Rajeev Motwani.  Approximate nearest neighbors: Towards removing the curse of dimensionality.  In _Proceedings of the Thirtieth Annual ACM Symposium on the Theory of Computing (STOC)_ , pages 604–613, Dallas, TX, 1998.
  * Jain et al. (2018) Prateek Jain, Om Dipakbhai Thakkar, and Abhradeep Thakurta.  Differentially private matrix completion revisited.  In _Proceedings of the 35th International Conference on Machine Learning (ICML)_ , pages 2220–2229, Stockholmsmässan, Stockholm, Sweden, 2018.
  * Jayaraman and Evans (2019) Bargav Jayaraman and David Evans.  Evaluating differentially private machine learning in practice.  In _Proceedings of the 28th USENIX Security Symposium (USENIX Security)_ , pages 1895–1912, Santa Clara, CA, 2019.
  * Johnson and Lindenstrauss (1984) William B. Johnson and Joram Lindenstrauss.  Extensions of Lipschitz mapping into Hilbert space.  _Contemporary Mathematics_ , 26:189–206, 1984.
  * Kairouz et al. (2014) Peter Kairouz, Sewoong Oh, and Pramod Viswanath.  Extremal mechanisms for local differential privacy.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , pages 2879–2887, Montreal, Canada, 2014.
  * Kasiviswanathan et al. (2013) Shiva Prasad Kasiviswanathan, Kobbi Nissim, Sofya Raskhodnikova, and Adam D. Smith.  Analyzing graphs with node differential privacy.  In _Proceedings of the 10th Theory of Cryptography Conference (TCC)_ , pages 457–476, Tokyo, Japan, 2013.
  * Kenny et al. (2021) Christopher T Kenny, Shiro Kuriwaki, Cory McCartan, Evan TR Rosenman, Tyler Simko, and Kosuke Imai.  The use of differential privacy for census data and its impact on redistricting: The case of the 2020 us census.  _Science advances_ , 7(41):eabk3283, 2021.
  * Kenthapadi et al. (2013) Krishnaram Kenthapadi, Aleksandra Korolova, Ilya Mironov, and Nina Mishra.  Privacy via the johnson-lindenstrauss transform.  _J. Priv. Confidentiality_ , 5(1), 2013.
  * Knudson et al. (2016) Karin Knudson, Rayan Saab, and Rachel Ward.  One-bit compressive sensing with norm estimation.  _IEEE Trans. Inf. Theory_ , 62(5):2748–2758, 2016.
  * Krizhevsky and Hinton (2009) Alex Krizhevsky and Geoffrey Hinton.  Learning multiple layers of features from tiny images.  _Technical Report, University of Toronto_ , 2009.
  * Laurent and Massart (2000) Béatrice Laurent and Pascal Massart.  Adaptive estimation of a quadratic functional by model selection.  _The Annals of Statistics_ , pages 1302–1338, 2000.
  * LeCun et al. (1998) Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner.  Gradient-based learning applied to document recognition.  _Proc. IEEE_ , 86(11):2278–2324, 1998.
  * Leng et al. (2014) Cong Leng, Jian Cheng, and Hanqing Lu.  Random subspace for binary codes learning in large scale image retrieval.  In _Proceedings of the 37th International ACM SIGIR Conference on Research and Development in Information Retrieval (SIGIR)_ , pages 1031–1034, Gold Coast, Australia, 2014.
  * Li (2019) Ping Li.  Sign-full random projections.  In _Proceedings of the Thirty-Third AAAI Conference on Artificial Intelligence (AAAI)_ , pages 4205–4212, Honolulu, HI, 2019.
  * Li and Li (2023a) Ping Li and Xiaoyun Li.  OPORP: One permutation + one random projection.  _arXiv preprint arXiv:2302.03505_ , 2023a.
  * Li and Zhao (2022) Ping Li and Weijie Zhao.  GCWSNet: Generalized consistent weighted sampling for scalable and accurate training of neural networks.  In _Proceedings of the 31st ACM International Conference on Information and Knowledge Management (CIKM)_ , Atlanta, GA, 2022.
  * Li et al. (2006) Ping Li, Trevor J Hastie, and Kenneth W Church.  Very sparse random projections.  In _Proceedings of the 12th ACM SIGKDD international conference on Knowledge discovery and data mining (KDD)_ , pages 287–296, Philadelphia, PA, 2006.
  * Li et al. (2014) Ping Li, Michael Mitzenmacher, and Anshumali Shrivastava.  Coding for random projections.  In _Proceedings of the 31th International Conference on Machine Learning (ICML)_ , pages 676–684, Beijing, China, 2014.
  * Li and Li (2019a) Xiaoyun Li and Ping Li.  Generalization error analysis of quantized compressive learning.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , pages 15124–15134, Vancouver, Canada, 2019a.
  * Li and Li (2019b) Xiaoyun Li and Ping Li.  Random projections with asymmetric quantization.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , pages 10857–10866, Vancouver, Canada, 2019b.
  * Li and Li (2021) Xiaoyun Li and Ping Li.  Quantization algorithms for random Fourier features.  In _Proceedings of the 38th International Conference on Machine Learning (ICML)_ , pages 6369–6380, Virtual Event, 2021.
  * Li and Li (2023b) Xiaoyun Li and Ping Li.  Differentially private one permutation hashing and bin-wise consistent weighted sampling.  _arXiv preprint_ , 2023b.
  * Mironov (2017) Ilya Mironov.  Rényi differential privacy.  In _Proceedings of the 30th IEEE Computer Security Foundations Symposium (CSF)_ , pages 263–275, Santa Barbara, CA, 2017.
  * Nissim et al. (2007) Kobbi Nissim, Sofya Raskhodnikova, and Adam D. Smith.  Smooth sensitivity and sampling in private data analysis.  In _Proceedings of the 39th Annual ACM Symposium on Theory of Computing (STOC)_ , pages 75–84, San Diego, CA, 2007.
  * Owen (1980) Donald Bruce Owen.  A table of normal integrals: A table.  _Communications in Statistics-Simulation and Computation_ , 9(4):389–419, 1980.
  * Rabanser et al. (2019) Stephan Rabanser, Stephan Günnemann, and Zachary C. Lipton.  Failing loudly: An empirical study of methods for detecting dataset shift.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , pages 1394–1406, Vancouver, Canada, 2019.
  * Rothchild et al. (2020) Daniel Rothchild, Ashwinee Panda, Enayat Ullah, Nikita Ivkin, Ion Stoica, Vladimir Braverman, Joseph Gonzalez, and Raman Arora.  FetchSGD: Communication-efficient federated learning with sketching.  In _Proceedings of the 37th International Conference on Machine Learning (ICML)_ , pages 8253–8265, Virtual Event, 2020.
  * Shrivastava and Li (2014) Anshumali Shrivastava and Ping Li.  In defense of minhash over simhash.  In _Proceedings of the Seventeenth International Conference on Artificial Intelligence and Statistics (AISTATS)_ , pages 886–894, Reykjavik, Iceland, 2014.
  * Singhal et al. (2021) Karan Singhal, Hakim Sidahmed, Zachary Garrett, Shanshan Wu, John Rush, and Sushant Prakash.  Federated reconstruction: Partially local federated learning.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , virtual, 2021.
  * Slawski and Li (2018) Martin Slawski and Ping Li.  On the trade-off between bit depth and number of samples for a basic approach to structured signal recovery from b-bit quantized linear measurements.  _IEEE Trans. Inf. Theory_ , 64(6):4159–4178, 2018.
  * Smith et al. (2020) Adam D. Smith, Shuang Song, and Abhradeep Thakurta.  The flajolet-martin sketch itself preserves differential privacy: Private counting with minimal space.  In _Advances in Neural Information Processing Systems_ , virtual, 2020\.
  * Soria-Comas et al. (2017) Jordi Soria-Comas, Josep Domingo-Ferrer, David Sánchez, and David Megías.  Individual differential privacy: A utility-preserving formulation of differential privacy guarantees.  _IEEE Trans. Inf. Forensics Secur._ , 12(6):1418–1429, 2017.
  * Stausholm (2021) Nina Mesing Stausholm.  Improved differentially private euclidean distance approximation.  In _Proceedings of the 40th ACM SIGMOD-SIGACT-SIGAI Symposium on Principles of Database Systems (PODS)_ , pages 42–56, Virtual Event, China, 2021.
  * Tomita et al. (2020) Tyler M. Tomita, James Browne, Cencheng Shen, Jaewon Chung, Jesse Patsolic, Benjamin Falk, Carey E. Priebe, Jason Yim, Randal C. Burns, Mauro Maggioni, and Joshua T. Vogelstein.  Sparse projection oblique randomer forests.  _J. Mach. Learn. Res._ , 21:104:1–104:39, 2020.
  * Vempala (2005) Santosh S Vempala.  _The random projection method_ , volume 65.  American Mathematical Soc., 2005.
  * Warner (1965) Stanley L Warner.  Randomized response: A survey technique for eliminating evasive answer bias.  _Journal of the American Statistical Association_ , 60(309):63–69, 1965.
  * Wei et al. (2020) Kang Wei, Jun Li, Ming Ding, Chuan Ma, Howard H. Yang, Farhad Farokhi, Shi Jin, Tony Q. S. Quek, and H. Vincent Poor.  Federated learning with differential privacy: Algorithms and performance analysis.  _IEEE Trans. Inf. Forensics Secur._ , 15:3454–3469, 2020\.
  * Weinberger et al. (2009) Kilian Q. Weinberger, Anirban Dasgupta, John Langford, Alexander J. Smola, and Josh Attenberg.  Feature hashing for large scale multitask learning.  In _Proceedings of the 26th Annual International Conference on Machine Learning (ICML)_ , pages 1113–1120, Montreal, Canada, 2009.
  * Wu et al. (2019) Jun Wu, Jingrui He, and Jiejun Xu.  DEMO-Net: Degree-specific graph neural networks for node and graph classification.  In _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining (KDD)_, pages 406–415, Anchorage, AK, 2019.
  * Xu et al. (2013) Jia Xu, Zhenjie Zhang, Xiaokui Xiao, Yin Yang, Ge Yu, and Marianne Winslett.  Differentially private histogram publication.  _VLDB J._ , 22(6):797–822, 2013.
  * Xu et al. (2021) Zhaozhuo Xu, Beidi Chen, Chaojian Li, Weiyang Liu, Le Song, Yingyan Lin, and Anshumali Shrivastava.  Locality sensitive teaching.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , pages 18049–18062, virtual, 2021.
  * Yang et al. (2015) Bin Yang, Issei Sato, and Hiroshi Nakagawa.  Bayesian differential privacy on correlated data.  In _Proceedings of the 2015 ACM SIGMOD International Conference on Management of Data (SIGMOD)_ , pages 747–762, Melbourne, Australia, 2015.
  * Zhang and Li (2020) Hang Zhang and Ping Li.  Optimal estimator for unlabeled linear regression.  In _Proceedings of the 37th International Conference on Machine Learning (ICML)_ , pages 11153–11162, Virtual Event, 2020.
  * Zhang et al. (2012) Jun Zhang, Zhenjie Zhang, Xiaokui Xiao, Yin Yang, and Marianne Winslett.  Functional mechanism: Regression analysis under differential privacy.  _Proc. VLDB Endow._ , 5(11):1364–1375, 2012\.
  * Zhang et al. (2014) Jun Zhang, Graham Cormode, Cecilia M. Procopiuc, Divesh Srivastava, and Xiaokui Xiao.  Privbayes: private data release via bayesian networks.  In _International Conference on Management of Data (SIGMOD)_ , pages 1423–1434, Snowbird, UT, 2014.
  * Zhang et al. (2022) Shan Zhang, Lei Wang, Naila Murray, and Piotr Koniusz.  Kernelized few-shot object detection with efficient integral aggregation.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 19185–19194, New Orleans, LA, 2022.
  * Zhang et al. (2021) Shijie Zhang, Hongzhi Yin, Tong Chen, Zi Huang, Lizhen Cui, and Xiangliang Zhang.  Graph embedding for recommendation against attribute inference attacks.  In _Proceedings of the Web Conference (WWW)_ , pages 3002–3014, Virtual Event / Ljubljana, Slovenia, 2021.
  * Zhang et al. (2020) Zhaoqi Zhang, Panpan Qi, and Wei Wang.  Dynamic malware analysis with feature engineering and feature learning.  In _Proceedings of the Thirty-Fourth AAAI Conference on Artificial Intelligence (AAAI)_ , pages 1210–1217, New York, NY, 2020.
  * Zhao et al. (2022) Fuheng Zhao, Dan Qiao, Rachel Redberg, Divyakant Agrawal, Amr El Abbadi, and Yu-Xiang Wang.  Differentially private linear sketches: Efficient implementations and applications.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , 2022.
  * Zymnis et al. (2010) Argyrios Zymnis, Stephen P. Boyd, and Emmanuel J. Candès.  Compressed sensing with quantized measurements.  _IEEE Signal Process. Lett._ , 17(2):149–152, 2010.

##  Appendix A Deferred Proofs

###  A.1 Proof of Lemma 3.6

Lemma 3.6 (Half-normal tail bound). Let $X_{1},...,X_{n}$ be iid $N(0,1)$
random variables and let $Y_{i}=|X_{i}|$, $\forall i$. Denote
$Z=\sum_{i=1}^{n}Y_{i}$. Then for any $t>0$, it holds that

$\displaystyle Pr(Z\geq\sqrt{2n^{2}\log 2+2nt})\leq\exp(-t).$

###### Proof.

By the Markov inequality, we write (with some $\zeta>0$)

$\displaystyle Pr(Z\geq t)=Pr(\exp(\zeta\sum_{i=1}^{n}Y_{i})\geq\exp(\zeta
t))\leq\frac{\mathbb{E}[\exp(\zeta\sum_{i=1}^{n}Y_{i}]}{\exp(\zeta t)}.$

By the independence of $Y_{i}$ and the moment generating function of half-
normal distribution, we have

$\displaystyle\mathbb{E}[\exp(\zeta\sum_{i=1}^{n}Y_{i}]=\prod_{i=1}^{n}\mathbb{E}[\exp(\zeta
Y_{i})]=\exp(\frac{n\zeta^{2}}{2})(2\Phi(\zeta))^{n}$
$\displaystyle=\exp(\frac{n\zeta^{2}}{2}+n\log 2\Phi(\zeta))$
$\displaystyle\leq\exp(\frac{n\zeta^{2}}{2}+n\log 2).$

Hence, we have

$\displaystyle Pr(Z\geq t)$
$\displaystyle\leq\exp(\frac{n\zeta^{2}}{2}-t\zeta+n\log 2),$

which is minimized at $\zeta=\frac{t}{n}$, leading to

$\displaystyle Pr(Z\geq t)$ $\displaystyle\leq\exp(-\frac{t^{2}}{2n}+n\log
2).$

Taking $t=\sqrt{2n^{2}\log 2+2nt}$ completes the proof. ∎

###  A.2 Proof of Lemma 4.2

###### Proof.

The bivariate normal density function is

$\displaystyle f(x,y)$
$\displaystyle=\frac{1}{2\pi\sigma_{x}\sigma_{y}\sqrt{1-\rho^{2}}}\exp\left(-\frac{\frac{x^{2}}{\sigma_{x}^{2}}-\frac{2\rho
xy}{\sigma_{x}\sigma_{y}}+\frac{y^{2}}{\sigma_{y}^{2}}}{2(1-\rho^{2})}\right)$
$\displaystyle=\frac{1}{2\pi\sigma_{x}\sigma_{y}\sqrt{1-\rho^{2}}}\exp\left(-\frac{x^{2}}{2\sigma_{x}^{2}}\right)\exp\left(-\frac{(\frac{y}{\sigma_{y}}-\rho\frac{x}{\sigma_{x}})^{2}}{2(1-\rho^{2})}\right).$

Therefore, we have $\mathbb{E}\left[|X|\big{|}|X|>|Y|\right]=\frac{A}{P}$,
with

$\displaystyle A$
$\displaystyle=\int_{-\infty}^{\infty}\frac{|x|}{\sqrt{2\pi}\sigma_{x}\sqrt{1-\rho^{2}}}\exp\left(-\frac{x^{2}}{2\sigma_{x}^{2}}\right)dx\int_{-|x|}^{|x|}\frac{1}{\sqrt{2\pi}\sigma_{y}}\exp\left(-\frac{(\frac{y}{\sigma_{y}}-\rho\frac{x}{\sigma_{x}})^{2}}{2(1-\rho^{2})}\right)dy,$
$\displaystyle P$
$\displaystyle=\int_{-\infty}^{\infty}\frac{1}{\sqrt{2\pi}\sigma_{x}\sqrt{1-\rho^{2}}}\exp\left(-\frac{x^{2}}{2\sigma_{x}^{2}}\right)dx\int_{-|x|}^{|x|}\frac{1}{\sqrt{2\pi}\sigma_{y}}\exp\left(-\frac{(\frac{y}{\sigma_{y}}-\rho\frac{x}{\sigma_{x}})^{2}}{2(1-\rho^{2})}\right)dy.$

Note that $P=Pr(|X|>|Y|)$ in the first statement of the theorem. Our
calculation will use the following two identities involving the Gaussian
functions (Owen, 1980):

$\displaystyle\int_{0}^{\infty}\phi(ax)\Phi(bx)dx=\frac{1}{2\pi|a|}\left(\frac{\pi}{2}+\tan^{-1}\left(\frac{b}{|a|}\right)\right),$
(33)
$\displaystyle\int_{0}^{\infty}x\phi(ax)\Phi(bx)dx=\frac{1}{2\sqrt{2\pi}}\left(1+\frac{b}{\sqrt{1+b^{2}}}\right),$
(34)

where $\phi(x)$ and $\Phi(x)$ are the pdf and cdf of the standard Gaussian
distribution.

With a proper change of random variables, we can compute $A$ as

$\displaystyle A=$
$\displaystyle\int_{-\infty}^{\infty}\frac{|x|}{\sqrt{2\pi}\sigma_{x}\sqrt{1-\rho^{2}}}\exp\left(-\frac{x^{2}}{2\sigma_{x}^{2}}\right)dx\int_{\frac{\frac{-|x|}{\sigma_{y}}-\rho\frac{x}{\sigma_{x}}}{\sqrt{1-\rho^{2}}}}^{\frac{\frac{|x|}{\sigma_{y}}-\rho\frac{x}{\sigma_{x}}}{\sqrt{1-\rho^{2}}}}\sqrt{1-\rho^{2}}\frac{1}{\sqrt{2\pi}}e^{-s^{2}}ds$
$\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\frac{|x|}{\sqrt{2\pi}\sigma_{x}}\exp\left(-\frac{x^{2}}{2\sigma_{x}^{2}}\right)\left[\Phi\left(\frac{\frac{|x|}{\sigma_{y}}-\rho\frac{x}{\sigma_{x}}}{\sqrt{1-\rho^{2}}}\right)-\Phi\left(\frac{\frac{-|x|}{\sigma_{y}}-\rho\frac{x}{\sigma_{x}}}{\sqrt{1-\rho^{2}}}\right)\right]dx$
$\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\frac{\sigma_{x}|t|}{\sqrt{2\pi}}\exp\left(-\frac{t^{2}}{2}\right)\left[\Phi\left(\frac{\frac{\sigma_{x}}{\sigma_{y}}|t|-\rho
t}{\sqrt{1-\rho^{2}}}\right)-\Phi\left(-\frac{\frac{\sigma_{x}}{\sigma_{y}}|t|-\rho
t}{\sqrt{1-\rho^{2}}}\right)\right]dt$ $\displaystyle\mathrel{\mathop{:}}=$
$\displaystyle A_{1}-A_{2}.$

For the first term we have

$\displaystyle A_{1}$
$\displaystyle=\sigma_{x}\left[\int_{0}^{\infty}\frac{t}{\sqrt{2\pi}}\exp\left(-\frac{t^{2}}{2}\right)\Phi\left(\frac{\frac{\sigma_{x}}{\sigma_{y}}-\rho}{\sqrt{1-\rho^{2}}}t\right)dt+\int_{-\infty}^{0}\frac{-t}{\sqrt{2\pi}}\exp\left(-\frac{t^{2}}{2}\right)\Phi\left(-\frac{\frac{\sigma_{x}}{\sigma_{y}}+\rho}{\sqrt{1-\rho^{2}}}t\right)dt\right]$
$\displaystyle=\sigma_{x}\left[\int_{0}^{\infty}\frac{t}{\sqrt{2\pi}}\exp\left(-\frac{t^{2}}{2}\right)\Phi\left(\frac{\frac{\sigma_{x}}{\sigma_{y}}-\rho}{\sqrt{1-\rho^{2}}}t\right)dt+\int_{0}^{\infty}\frac{s}{\sqrt{2\pi}}\exp\left(-\frac{s^{2}}{2}\right)\Phi\left(\frac{\frac{\sigma_{x}}{\sigma_{y}}+\rho}{\sqrt{1-\rho^{2}}}s\right)ds\right]$
$\displaystyle=\sigma_{x}\left[\frac{1}{2\sqrt{2\pi}}\left(1+\frac{\frac{\frac{\sigma_{x}}{\sigma_{y}}-\rho}{\sqrt{1-\rho^{2}}}}{\sqrt{1+\frac{(\frac{\sigma_{x}}{\sigma_{y}}-\rho)^{2}}{1-\rho^{2}}}}\right)+\frac{1}{2\sqrt{2\pi}}\left(1+\frac{\frac{\frac{\sigma_{x}}{\sigma_{y}}+\rho}{\sqrt{1-\rho^{2}}}}{\sqrt{1+\frac{(\frac{\sigma_{x}}{\sigma_{y}}+\rho)^{2}}{1-\rho^{2}}}}\right)\right]$
$\displaystyle=\sigma_{x}\left[\frac{1}{\sqrt{2\pi}}+\frac{1}{2\sqrt{2\pi}}\left(\frac{r-\rho}{\sqrt{1+r^{2}-2r\rho}}+\frac{r+\rho}{\sqrt{1+r^{2}+2r\rho}}\right)\right],$

where we denote $r=\frac{\sigma_{x}}{\sigma_{y}}$ and use (34). Similarly, we
have that

$\displaystyle A_{2}$
$\displaystyle=\sigma_{x}\left[\int_{0}^{\infty}\frac{t}{\sqrt{2\pi}}\exp\left(-\frac{t^{2}}{2}\right)\Phi\left(-\frac{r+\rho}{\sqrt{1-\rho^{2}}}t\right)dt+\int_{-\infty}^{0}\frac{-t}{\sqrt{2\pi}}\exp\left(-\frac{t^{2}}{2}\right)\Phi\left(\frac{r-\rho}{\sqrt{1-\rho^{2}}}t\right)dt\right]$
$\displaystyle=\sigma_{x}\left[\int_{0}^{\infty}\frac{t}{\sqrt{2\pi}}\exp\left(-\frac{t^{2}}{2}\right)\Phi\left(-\frac{r+\rho}{\sqrt{1-\rho^{2}}}t\right)dt+\int_{0}^{\infty}\frac{s}{\sqrt{2\pi}}\exp\left(-\frac{s^{2}}{2}\right)\Phi\left(-\frac{r-\rho}{\sqrt{1-\rho^{2}}}s\right)ds\right]$
$\displaystyle=\sigma_{x}\left[\frac{1}{\sqrt{2\pi}}-\frac{1}{2\sqrt{2\pi}}\left(\frac{r-\rho}{\sqrt{1+r^{2}-2r\rho}}+\frac{r+\rho}{\sqrt{1+r^{2}+2r\rho}}\right)\right].$

Therefore, we obtain

$\displaystyle
A(\rho,r)=A_{1}-A_{2}=\frac{\sigma_{x}}{\sqrt{2\pi}}\left(\frac{r-\rho}{\sqrt{1+r^{2}-2r\rho}}+\frac{r+\rho}{\sqrt{1+r^{2}+2r\rho}}\right).$
(35)

To compute $P$, by doing a similar change of variables, we have

$\displaystyle P=$
$\displaystyle\int_{-\infty}^{\infty}\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{t^{2}}{2}\right)\left[\Phi\left(\frac{r|t|-\rho
t}{\sqrt{1-\rho^{2}}}\right)-\Phi\left(-\frac{r|t|-\rho
t}{\sqrt{1-\rho^{2}}}\right)\right]dt$ $\displaystyle\mathrel{\mathop{:}}=$
$\displaystyle P_{1}-P_{2}.$

Using (33), we obtain

$\displaystyle P_{1}$
$\displaystyle=\int_{0}^{\infty}\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{t^{2}}{2}\right)\Phi\left(\frac{r-\rho}{\sqrt{1-\rho^{2}}}t\right)dt+\int_{-\infty}^{0}\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{t^{2}}{2}\right)\Phi\left(-\frac{r+\rho}{\sqrt{1-\rho^{2}}}t\right)dt$
$\displaystyle=\int_{0}^{\infty}\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{t^{2}}{2}\right)\Phi\left(\frac{r-\rho}{\sqrt{1-\rho^{2}}}t\right)dt+\int_{0}^{\infty}\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{s^{2}}{2}\right)\Phi\left(\frac{r+\rho}{\sqrt{1-\rho^{2}}}s\right)ds$
$\displaystyle=\frac{1}{2\pi}\left(\frac{\pi}{2}+\tan^{-1}\left(\frac{r-\rho}{\sqrt{1-\rho^{2}}}\right)\right)+\frac{1}{2\pi}\left(\frac{\pi}{2}+\tan^{-1}\left(\frac{r+\rho}{\sqrt{1-\rho^{2}}}\right)\right)$
$\displaystyle=\frac{1}{2}+\frac{1}{2\pi}\left[\tan^{-1}\left(\frac{r-\rho}{\sqrt{1-\rho^{2}}}\right)+\tan^{-1}\left(\frac{r+\rho}{\sqrt{1-\rho^{2}}}\right)\right],$
$\displaystyle P_{2}$
$\displaystyle=\frac{1}{2}-\frac{1}{2\pi}\left[\tan^{-1}\left(\frac{r-\rho}{\sqrt{1-\rho^{2}}}\right)+\tan^{-1}\left(\frac{r+\rho}{\sqrt{1-\rho^{2}}}\right)\right],$

which leads to

$\displaystyle P(\rho,r)=P_{1}-P_{2}=$
$\displaystyle\frac{1}{\pi}\left[\tan^{-1}\left(\frac{r-\rho}{\sqrt{1-\rho^{2}}}\right)+\tan^{-1}\left(\frac{r+\rho}{\sqrt{1-\rho^{2}}}\right)\right]$
(36)

Therefore, we know that

$\displaystyle\mathbb{E}\left[|X|\big{|}|X|>|Y|\right]=\frac{A(\rho,r)}{P(\rho,r)}=\sigma_{x}\sqrt{\frac{\pi}{2}}\cdot\frac{\frac{r-\rho}{\sqrt{1+r^{2}-2r\rho}}+\frac{r+\rho}{\sqrt{1+r^{2}+2r\rho}}}{\tan^{-1}\left(\frac{r-\rho}{\sqrt{1-\rho^{2}}}\right)+\tan^{-1}\left(\frac{r+\rho}{\sqrt{1-\rho^{2}}}\right)},$

with $r=\sigma_{x}/\sigma_{y}$. We now investigate the derivative of $P$. By
some algebra, we can show that

$\displaystyle\frac{\partial
P(\rho,r)}{\partial\rho}=\frac{2r\rho(r^{2}-1)}{(1+r^{2}-2r\rho)(1+r^{2}+2r\rho)\sqrt{1-\rho^{2}}}.$

When $0<r\leq 1$, $\frac{\partial P(\rho,r)}{\partial\rho}\geq 0$ when
$\rho\leq 0$ and $\frac{\partial P(\rho,r)}{\partial\rho}\leq 0$ when
$\rho>0$. Therefore, $\max_{\rho}P(\rho,r)=P(0,r)=\frac{2}{\pi}\tan^{-1}(r)$.

Tail bound. By our previous calculations, the conditional distribution of $X$
given $|X|>|Y|$ is

$\displaystyle
f(x\big{|}|X|>|Y|)=\frac{\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{x^{2}}{2}\right)\left[\Phi\left(\frac{r|x|-\rho
x}{\sqrt{1-\rho^{2}}}\right)-\Phi\left(-\frac{r|x|-\rho
x}{\sqrt{1-\rho^{2}}}\right)\right]}{P},\quad x\in\mathbb{R},$

with $P=Pr(|X|>|Y|)$ in (36) the normalizing constant to make the integral
equal to 1.

The conditional tail probability can be computed as follows. For some $t>0$,
by symmetry,

$\displaystyle Pr(|X|>t,|X|>|Y|)$ $\displaystyle=$ $\displaystyle
2\int_{t}^{\infty}\frac{1}{\sqrt{2\pi}\sigma_{x}}\exp\left(-\frac{x^{2}}{2\sigma_{x}^{2}}\right)\left[\Phi\left(\frac{r|x|-\rho
x}{\sqrt{1-\rho^{2}}}\right)-\Phi\left(-\frac{r|x|-\rho
x}{\sqrt{1-\rho^{2}}}\right)\right]dx$ $\displaystyle=$ $\displaystyle
2\int_{\frac{t}{\sigma_{x}}}^{\infty}\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{x^{2}}{2}\right)\left[\Phi\left(\frac{r|x|-\rho
x}{\sqrt{1-\rho^{2}}}\right)-\Phi\left(-\frac{r|x|-\rho
x}{\sqrt{1-\rho^{2}}}\right)\right]dx$ $\displaystyle\mathrel{\mathop{:}}=$
$\displaystyle 2(\tilde{P}_{1}-\tilde{P}_{2}).$

For $\tilde{P}_{1}$, using polar coordinates we have

$\displaystyle\tilde{P}_{1}=$
$\displaystyle\frac{1}{2\pi}\int_{\frac{t}{\sigma_{x}}}^{\infty}e^{-\frac{x^{2}}{2}}dx\int_{-\infty}^{\frac{r-\rho}{\sqrt{1-\rho^{2}}}x}e^{-\frac{y^{2}}{2}}dy$
$\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int_{-\frac{\pi}{2}}^{\tan^{-1}(\frac{r-\rho}{\sqrt{1-\rho^{2}}})}d\theta\int_{\frac{t}{\sigma_{x}\cos(\theta)}}^{\infty}e^{-\frac{r^{2}}{2}}rdr$
$\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int_{-\frac{\pi}{2}}^{\tan^{-1}(\frac{r-\rho}{\sqrt{1-\rho^{2}}})}\exp\left(-\frac{t^{2}}{2\sigma_{x}^{2}\cos^{2}(\theta)}\right)d\theta.$

Similarly,

$\displaystyle\tilde{P}_{2}$
$\displaystyle=\frac{1}{2\pi}\int_{-\frac{\pi}{2}}^{\tan^{-1}(-\frac{r+\rho}{\sqrt{1-\rho^{2}}})}\exp\left(-\frac{t^{2}}{2\sigma_{x}^{2}\cos^{2}(\theta)}\right)d\theta.$

Therefore, we obtain

$\displaystyle Pr(|X|>t,|X|>|Y|)$
$\displaystyle=\frac{1}{\pi}\int_{\tan^{-1}(-\frac{r+\rho}{\sqrt{1-\rho^{2}}})}^{\tan^{-1}(\frac{r-\rho}{\sqrt{1-\rho^{2}}})}\exp\left(-\frac{t^{2}}{2\sigma_{x}^{2}\cos^{2}(\theta)}\right)d\theta$
$\displaystyle=\frac{1}{\pi}\int_{-\tan^{-1}(\frac{r+\rho}{\sqrt{1-\rho^{2}}})}^{\tan^{-1}(\frac{r-\rho}{\sqrt{1-\rho^{2}}})}\exp\left(-\frac{t^{2}}{2\sigma_{x}^{2}\cos^{2}(\theta)}\right)d\theta$
$\displaystyle\leq
e^{-\frac{t^{2}}{2\sigma_{x}^{2}}}\frac{1}{\pi}\int_{-\tan^{-1}(\frac{r+\rho}{\sqrt{1-\rho^{2}}})}^{\tan^{-1}(\frac{r-\rho}{\sqrt{1-\rho^{2}}})}d\theta,$

since $\cos^{2}(\theta)\in[0,1]$. Notice that $P$ in (36) can be written as
$P=\frac{1}{\pi}\int_{-\tan^{-1}(\frac{r+\rho}{\sqrt{1-\rho^{2}}})}^{\tan^{-1}(\frac{r-\rho}{\sqrt{1-\rho^{2}}})}d\theta$.
Hence, we know that the conditional tail probability is

$\displaystyle
Pr(|X|>t\big{|}|X|>|Y|)=\frac{Pr(|X|>t,|X|>|Y|)}{Pr(|X|>|Y|)}\leq\exp\left(-\frac{t^{2}}{2\sigma_{x}^{2}}\right),\quad\forall
r>0,\rho\in(-1,1).$

At the boundaries $\rho=1$, $\rho=-1$, one can verify $Pr(|X|>|Y|)=0$. This
concludes the proof. ∎