# Synchronizing Chaos using Reservoir Computing

Amirhossein Nazerian  Mechanical Engineering Department, University of New
Mexico, Albuquerque, NM, 87131  Chad Nathe  Mechanical Engineering Department,
University of New Mexico, Albuquerque, NM, 87131  Joseph D. Hart  US Naval
Research Laboratory, Washington, DC 20375  Francesco Sorrentino  Mechanical
Engineering Department, University of New Mexico, Albuquerque, NM, 87131

###### Abstract

We attempt to achieve isochronal synchronization between a drive system
unidirectionally coupled to a response system, under the assumption that
limited knowledge on the states of the drive is available at the response.
Machine learning techniques have been previously implemented to estimate the
states of a dynamical system from limited measurements. We consider situations
in which knowledge of the non-measurable states of the drive system is needed
in order for the response system to synchronize with the drive. We use a
reservoir computer to estimate the non-measurable states of the drive system
from its measured states and then employ these measured states to synchronize
the response system with the drive.

A large literature has investigated synchronization of chaos, see e.g., [1, 2,
3, 4, 5, 6, 7] and the review paper [8]. However, little attention has been
devoted to the usage of machine learning techniques to enable synchronization.
A typical problem is that of synchronizing two identical chaotic systems, a
drive system that is unidirectionally coupled to a response system. This is
typically achieved by having the drive communicate part of its state to the
response. Here we use a reservoir observer within a control loop to
reconstruct other states of the drive at the response system and use these
reconstructed states to aid synchronization. We also show how this proposed
scheme can be used to control the time evolution of the response system on an
unstable periodic orbit embedded in the chaotic attractor. The robustness of
our proposed approach is studied against measurement noise affecting the
information transmitted from the drive to the response system.

##  I Introduction

The need of regulating the behavior of nonlinear systems is a common
requirement in many physical, social, and biological applications [9, 10, 11,
12]. Observers are classically used in controls applications to estimate
states of a dynamical system that cannot be directly measured; feedback
control can then be performed using the states that are directly available and
those that are reconstructed using the observer. An important problem in the
control of nonlinear systems, then, is the construction of a sufficiently
accurate observer.

In the absence of a physical model of the system to be observed, one must turn
to a data-driven approach. One such approach that is particularly effective
when limited data is available is that of reservoir computing. A reservoir
computer is a type of recurrent neural network that is designed to be easy to
train [13, 14]. While reservoir computers are most often used in an autonomous
mode [15], they have also been found to be an efficient and effective observer
of dynamical systems such as Lorenz oscillators [16, 17, 18], semiconductor
lasers [19], ecological models [20], and spatiotemporal systems [16, 21].

There is ample observation of synchronization in the natural world [22], yet
it is unclear whether this observed synchronization may be the result of an
underlying learning process. Autonomous reservoir computing models have been
found to be capable of synchronizing with the nonlinear systems upon which
they are trained [23, 24] as well as with other identical reservoir computing
models [25, 26, 27]. A modification of reservoir computing termed deep
reservoir computing has been used for direct learning of a control rule [28].
Additionally, autonomous reservoir computing has been used for the estimation
of unstable periodic orbits of an unknown dynamical system and for the pinning
control of a network on to the estimated unstable periodic orbit [29].

Despite these successes, and even though the combination of an observer and a
controller is a well established paradigm in control theory, little work has
been devoted to study how reservoir computing can be used as an observer for
control applications. In this work, we train a reservoir computer as an
observer, which is used to mediate the pinning control of a response system by
a drive system. The reservoir observer, driven by the measured variable of the
drive system, is used to estimate unmeasured variables of the drive system.
The measured variable and the estimated variables are then used to drive the
response system to a controlled trajectory. We find that the introduction of
the reservoir computer can lead to a dramatic reduction in the coupling
strength required.

More formally, we consider a problem in which a drive system is
unidirectionally coupled to a response system. The drive system produces a
reference trajectory for the response system with the goal of controlling or
synchronizing the time evolution of the response system on that of the drive
system. Another reason why this particular setting may be of interest is that
by choosing the initial condition of the drive system to lie on an unstable
periodic orbit (UPO) it may be possible to control the time evolution of the
response to converge on that particular UPO.

We take the response and the drive to be described by the same equations and
we assume that when these systems are uncoupled, their dynamics evolve on the
same chaotic attractor. We consider the following set of generic dynamical
equations that describe the time evolution of the drive, $\pmb{x}_{D}(t)$, and
the response, $\pmb{x}_{R}(t)$,

$\displaystyle\dot{\pmb{x}}_{D}(t)$ $\displaystyle=\pmb{F}(\pmb{x}_{D}(t))$
(1) $\displaystyle\dot{\pmb{x}}_{R}(t)$
$\displaystyle=\pmb{F}(\pmb{x}_{R}(t))+\kappa
H(\pmb{x}_{D}(t)-\pmb{x}_{R}(t))$

where $\pmb{F}:\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}$ and $m$ is the number
of states in the dynamics system. The matrix $H$ is size $m\times m$ and
describes the coupling scheme. The scalar $k>0$ is the coupling strength.

In the rest of this paper, without loss of generality we take $m=3$ and write
${\pmb{x}}_{D}(t)=[x_{D}(t),y_{D}(t),z_{D}(t)]$,
${\pmb{x}}_{R}(t)=[x_{R}(t),y_{R}(t),z_{R}(t)]$. We assume the entries of the
matrix $H=\\{H_{ij}\\}$ to be either zeros or ones, and consider different
coupling schemes, as a result of different choices of the matrix $H$. For
example, by setting all the entries of the matrix $H$ equal to zero except for
entry $H_{12}=1$, we indicate that the state $y_{D}(t)$ from the drive system
is available to the response system and that $y_{D}(t)$ appears in the
equation of the first state of the response system. We will also refer to this
coupling scheme as $y_{D}(t)\rightarrow x_{R}(t)$. We will instead refer to a
coupling scheme as $\hat{y}_{D}(t)\rightarrow x_{R}(t)$ when $y_{D}(t)$ is not
a measured variable from the response system, but an estimate of
$\hat{y}_{D}(t)$.

The rest of this paper is organized as follows. In Sec. II, reservoir
computing and the general procedure of the drive-response system are
discussed. Examples of the Chen and Rössler systems are also provided, and the
effect of measurement noise on the performance of the reservoir computer is
discussed. Conclusions are provided in Sec. III. In Appendix A, we discuss how
the integration of the continuous time system and the evolution of the
discrete time RC are performed.

##  II Reservoir Computing

We begin by introducing the reservoir equation,

$\pmb{r}(t+\Delta
t)=(1-\alpha)\pmb{r}(t)+\alpha\tanh{(A\pmb{r}(t)+s(t)\pmb{w}_{\text{in}})}$
(2)

where $A$ is the coupling matrix, $s(t)$ is the drive signal, $\Delta t$ is
the discrete time step, and $\pmb{w}_{\text{in}}$ is a vector of random
elements drawn from a standard Gaussian distribution, i.e.,
$\mathcal{N}(1,0)$. We set $\Delta t=0.001$ throughout this paper unless
stated otherwise. The matrix $A$ is the adjacency matrix of a directed Erdos
Renyi network with $N$ nodes and connectivity probability $p=0.1$. We make $A$
to be have a spectral radius of $\rho$ via the operation, $A\leftarrow
A\rho/\lambda$, where $\lambda$ is the largest real eigenvalue and $0<\rho\leq
1$.

From the time evolution of the reservoir equation, Eq. (2), one can construct
the readout matrix,

$\Omega=\begin{bmatrix}r_{1}(1)&r_{2}(1)&\cdots&r_{N}(1)&1\\\
r_{1}(2)&r_{2}(2)&\cdots&r_{N}(2)&1\\\ \vdots&\vdots&\ddots&\vdots&\vdots\\\
r_{1}(T_{1})&r_{2}(T_{1})&\cdots&r_{N}(T_{1})&1\\\ \end{bmatrix}$ (3)

where $r_{i}(t)$ is the readout of node $i$ at time $t$ and $t=T_{1}$
indicates the end of the training phase. The last column of the matrix
$\Omega$ is set to 1 to account for any constant offset in the fit. We then
relate the readouts to the training signal, $y(t)$, with additive noise, via
the unknown coefficients contained in the vector,
$\pmb{w}_{\textnormal{out}}$,

$\Omega\pmb{w}_{\textnormal{out}}=\mathbf{g}$ (4)

where $\mathbf{g}(t)$ is the training signal. We then compute the unknown
coefficients vector $\pmb{w}_{\textnormal{out}}$ via the equation,

$\pmb{w}_{\textnormal{out}}=\mathbf{\Omega}^{\dagger}\mathbf{g}.$ (5)

Here, $\mathbf{\Omega}^{\dagger}$ is given as,

$\mathbf{\Omega}^{\dagger}=\mathbf{\left(\Omega^{T}\Omega+\beta
I\right)^{-1}\Omega^{T}}.$ (6)

In the above equation, $\beta$ is the ridge-regression parameter used to avoid
overfitting [16] and $\mathbf{I}$ is the identity matrix. For all the
simulations in this work, $\beta=10^{-9}$ is used. For $\beta=0$,
$\mathbf{\Omega}^{\dagger}$ is the pseudo inverse matrix of $\mathbf{\Omega}$.
Next, we define the training fit signal as,

$\mathbf{h}=\Omega\pmb{w}_{\textnormal{out}}.$ (7)

Lastly, the training error is computed as,

$\Delta_{\text{tr}}=\frac{\textnormal{std}(\mathbf{h}-{\mathbf{g}})}{\textnormal{std}(\mathbf{g})}$
(8)

where the notation $\textnormal{std}(\cdot)$ denotes the standard deviation.

The testing phase is carried out in the same manner, except using the
previously computed coefficients contained in, $\pmb{w}_{\textnormal{out}}$,
to compute the fit signal, $\mathbf{h}_{\text{ts}}$ according to the
$\Omega_{\text{ts}}$ matrix generated by the new drive signal, i.e., the
testing signal $\mathbf{g}_{\text{ts}}(t)$. We call the length of the testing
phase $T_{s}$. The testing error is computed as

$\Delta_{\text{ts}}=\frac{\textnormal{std}(\mathbf{h}_{\text{ts}}-{\mathbf{g}_{\text{ts}}})}{\textnormal{std}(\mathbf{g}_{\text{ts}})}.$
(9)

In Fig. 1, we show the process of training the reservoir computer in (a) and
then using it for control in (b). I/R is the input-to-reservoir function and
is represented by the term $s(t)\pmb{w}_{\text{in}}$ in Eq. (2). R/O is the
reservoir-to-output function represented by Eq. (7).

Figure 1: (a) Training phase of reservoir. The reservoir input signal is
$x_{D}(t)$ and it is trained on $y_{D}(t)$. The error, $E(t)$, is calculated
with the fit signal, $\hat{y}(t)$, and the training signal. (b) Control
configuration. The response system takes as an input, $E(t)$, where $E(t)$ is
$\pmb{x}_{R}(t)-\pmb{x}_{D}(t)$. The matrix, $H$, describes what state
variables are used in coupling.

###  II.1 Chen system

We take both the drive and the response systems to be described by the Chen
equation $\dot{\pmb{x}}(t)=\pmb{F}(\pmb{x}(t))$ where $\pmb{x}=[x,\ y,\
z]^{\top}$ and

$\pmb{F}(\pmb{x}(t))=\begin{bmatrix}a(y(t)-x(t))\\\ (c-a-z(t))+cy(t)\\\
x(t)y(t)-\beta z(t)\end{bmatrix},$ (10)

where here we have set $a=35,c=28$ and $\beta=8/3$. The goal is to synchronize
the trajectories of the response system and the drive system, using the scheme
shown in Fig. 1. It is assumed that the only component of the drive system
accessible to the response system is the time evolution of the state
$x_{D}(t)$.

The master stability function [30] predicts that by coupling $x_{D}\rightarrow
y_{R}$, the two systems synchronize when the coupling strength $\kappa>10.62$
[31]. However, if the information about $y_{D}$ is available, the two systems
synchronize by the coupling $y_{D}\rightarrow y_{R}$ when $\kappa>3.54$ [31].
Hence, we would like to estimate $\hat{y}_{D}$ from $x_{D}$ in order to
synchronize the two systems with a lower coupling strength $3.54<\kappa<10.62$
through coupling $\hat{y}_{D}\rightarrow y_{R}$. Note that coupling
$x_{D}\rightarrow x_{R}$ does not result in synchronization for any value of
$k$, as shown in [31].

(A) (B)
---
(C)
Figure 2:  Chen system. (A) Coupled $x_{D}\rightarrow y_{R}$,
$\kappa_{\min}\simeq 10.78$ (B) Coupled $\hat{y}_{D}\rightarrow
y_{R},\kappa_{\min}\simeq 3.96$ (C) Coupled both $x_{D}\rightarrow y_{R}$ and
$\hat{y}_{D}\rightarrow y_{R},\kappa_{\min}\simeq 3.00$. The errorbars
represent standard deviations.

We introduce the instantaneous and average synchronization errors,

$\displaystyle\begin{split}E(t):=&\left[\left(\frac{x_{D}(t)-x_{R}(t)}{\textnormal{std}(x_{D}(t))_{t}}\right)^{2}+\left(\frac{y_{D}(t)-y_{R}(t)}{\textnormal{std}(y_{D}(t))_{t}}\right)^{2}\right.\\\
&+\left.\left(\frac{z_{D}(t)-z_{R}(t)}{\textnormal{std}(z_{D}(t))_{t}}\right)^{2}\right]^{\frac{1}{2}}\\\
\bar{E}=&<E(t)>_{t},\end{split}$ (11)

respectively, where $<\cdot>_{t}$ indicates an average over the time interval
$t=[500,1000]$ seconds, and $\textnormal{std}(\cdot)_{t}$ returns the standard
deviation (a positive scalar value) of the time-series in the argument over
the time interval $t$.

Figure 2 shows the synchronization error $\bar{E}$ for different coupling
schemes and different coupling strengths, when the drive and the response
system are initialized from randomly chosen initial conditions from the Chen
system attractor. The results we obtain are in agreement with the predictions
of the master stability function. From each plot in Fig. 2 we see that there
is a critical coupling strength $\kappa_{\min}$ above which a transition from
asynchrony to synchrony is observed. Fig. 2 (A) corresponds to coupling
$x_{D}\rightarrow y_{R}$, for which $\kappa_{\min}\simeq 10.78$; Fig. 2 (B)
corresponds to coupling $y_{D}\rightarrow y_{R}$, for which
$\kappa_{\min}\simeq 3.96$. Finally, Fig. 2 (C) corresponds to coupling in
both $x_{D}\rightarrow y_{R}$ and $y_{D}\rightarrow y_{R},$ for which
$\kappa_{\min}\simeq 3.00$.

By comparing Figs. 2 (B) and 2 (C) with Fig. 2 (A) we see that a much lower
coupling strength is needed for synchronization when coupling involves the
state $y_{D}(t)$, hence in what follows we attempt to reconstruct
$\hat{y}_{D}(t)\approx y_{D}(t)$ from knowledge of $x_{D}(t)$.

Our approach illustrated in Fig.1 is to indirectly couple the response system
to the drive system through the estimate $\hat{y}_{D}(t)$ produced by the
reservoir observer. Hence, the coupled system equations are,

$\displaystyle\dot{\pmb{x}}_{D}(t)$ $\displaystyle=\pmb{F}(\pmb{x}_{D}(t))$
(12) $\displaystyle\dot{\pmb{x}}_{R}(t)$
$\displaystyle=\pmb{F}(\pmb{x}_{R}(t))+\kappa\begin{bmatrix}0\\\
\hat{y}_{D}-y_{R}+x_{D}-x_{R}\\\ 0\end{bmatrix}$

If the two systems are coupled ‘ideally’, then $\hat{y}_{D}$ is replaced by
the true value $y_{D}$.

As an example, Fig. 3 shows the time evolutions of the $x$ components of the
drive and response system as they are coupled through the RC and ideally. We
tentatively set the RC parameters equal to $\alpha=0.61$, $\rho=0.9$, 500
nodes, and the input weights are drawn randomly from a uniform distribution
between $-0.5$ and $0.5$. The two systems are initialized from different
initial conditions on the chaotic attractor of the Chen system. We see that
the two coupling schemes produce comparable levels of synchronization between
the drive and response systems, which appear qualitatively similar to the
human eye. However, a calculation shows that the two synchronization errors
differ by several orders of magnitude. This can be seen in Figure 4, which
compares the synchronization errors $E_{\textnormal{RC}}(t)$ and
$E_{\textnormal{Ideal}}(t)$ through a 1000$s$ simulations. Here, the
instantanous synchronization errors are $E_{\textnormal{RC}}(t)$ when coupled
through the RC, and $E_{\textnormal{Ideal}}(t)$ when coupled ideally. We see
that while the synchronization error attains a small value throughout the
simulation when coupling uses the estimated state, the error attains a much
smaller value (equal to the numerical precision of the computer) when coupling
uses the true state. This is due to the fact that the estimate produced by the
RC observer is not exactly the same as the true signal, but just an
approximation. It is also important to point out that the ideal case shown in
Fig. 4 does not consider the presence of noise in the signal transmitted from
the drive system to the response system. This irrealistic assumption is
removed in Sec. II.4.

Figure 3: The plots show the time evolution of the $x$ component of the drive
and the response systems when coupled through the RC in $x_{D}\rightarrow
y_{R}$ and $\hat{y}_{D}\rightarrow y_{R}$ (top) and when coupled ideally in
$x_{D}\rightarrow y_{R}$ and $y_{D}\rightarrow y_{R}$ (bottom.) Here, the
coupling strength $\kappa=3.1$ for both cases. Figure 4: Synchronization error
of Chen’s drive and response systems from time series in Fig. 3. The plot
shows $E_{\textnormal{RC}}(t)$, the error when the coupling $x_{D}\rightarrow
y_{R}$ and $\hat{y}_{D}\rightarrow y_{R}$ is through estimation by RC, and
$E_{\textnormal{Ideal}}(t)$, the error in the ideal case when all the states
of the drive system are known, $x_{D}\rightarrow y_{R}$ and $y_{D}\rightarrow
y_{R}$.

The results in Figs. 2 and 3 show the feasibility of our approach. An
important step is to properly choose the hyperparameters of the RC observer.
Figure 5 (A) shows the training, testing, and synchronization errors,
$\Delta_{\textnormal{tr}}$, $\Delta_{\textnormal{ts}}$, and $\bar{E}$,
respectively, as the parameter $\alpha$ in (2) is varied. We see that as
$\alpha$ grows to 0.4, the training and synchronization errors are
substantially reduced. However, for $0.4\leq\alpha\leq 1$, the testing errors
attain higher standard deviations, especially for larger values of $\alpha$.
Figure 5 (B) shows the effect of varying the spectral radius $\rho$ over the
training, testing, and synchronization errors, for $\alpha=0.5$. We see that
the performance we achieve is quite robust to variations in $\rho$, as long as
$\rho\geq 0.1$. We set $\rho=0.9$ in our simulations.

(A)
---
(B)
Figure 5: Chen system. Coupling is in $x_{D}\rightarrow x_{R}$ and
$\hat{y}_{D}\rightarrow y_{R}$ with $\kappa=10$ (A) We plot the training
$\Delta_{\textnormal{tr}}$, testing $\Delta_{\textnormal{ts}}$, and
synchronization $E$ errors as a function of the leakage parameter, $\alpha$
for a fixed value of the spectral radius $\rho=0.9$. (B) We plot the training,
testing, and synchronization error as a function of the spectral radius,
$\rho$, of the $A$ matrix, for a fixed value of the leakage parameter
$\alpha=0.5$.

###  II.2 Rössler

The dynamical equation of a Rössler system is
$\dot{\pmb{x}}(t)=\pmb{F}(\pmb{x}(t))$ where $\pmb{x}=[x,\ y,\ z]^{\top}$ and

$\pmb{F}(\pmb{x}(t))=\begin{bmatrix}-y(t)-z(t))\\\ x(t)+ay(t)\\\
b+(x(t)-c)z(t)\end{bmatrix}$ (13)

where we set $a=b=0.2$ and $c=9$. In what follows, we set the initial
conditions of the drive and response system to be randomly chosen points on
the chaotic attractor. We examine the effect of the gain, $\kappa$, on the
synchronization error in the scenarios where we couple $x_{D}\rightarrow
x_{R}$, $y_{D}\rightarrow y_{R}$, and both simultaneously. We then move on to
show the performance of the RC when $x_{D}\rightarrow x_{R}$ and
$\hat{y}_{D}\rightarrow y_{R}$ are the couplings.

In Fig. 6 we perform a similar study as what previously shown in Fig. 5 to
select the hyperparameters of the RC for the case of the Rössler system.

(A)
---
(B)
Figure 6: Rössler system. Coupling is in $x_{D}\rightarrow x_{R}$ and
$\hat{y}_{D}\rightarrow y_{R}$ with $\kappa=0.15$. (A) We plot the training,
testing, and synchronization error as a function of the leakage parameter,
$\alpha$. We set $\rho=0.6$ (B) We plot the training, testing, and
synchronization error as a function of the spectral radius, $\rho$, of the $A$
matrix. We set $\alpha=0.05$ (A) (B)
---
(C)
Figure 7: Rössler system. (A) Coupling $x_{D}\rightarrow x_{R}$ (B) Coupling
$\hat{y}_{D}\rightarrow y_{R}$ (C) Coupling is both $x_{D}\rightarrow x_{R}$
and $\hat{y}_{D}\rightarrow y_{R}$. The errorbars represent the standard
deviations.

Figure 7 shows a clear advantage of incorporating the RC into the control
loop. In this figure, we set $\alpha=0.002$ and $\rho=0.9$. In plot (A) we
show the synchronization error versus the coupling strength for the accessible
coupling, $x_{D}\rightarrow x_{R}$. In plot (B) we show that $\hat{y}_{D}$,
generated by the RC, is sufficient enough to drive the synchronization error
close to zero. Plot (C) shows that when both couplings are used, a lower value
for the coupling strength is needed to achieve synchronization.

###  II.3 Unstable periodic orbits

Reservoir computing has been previously used to detect unstable periodic
orbits (UPOs) of chaotic systems [29]. In this section, we aim to synchronize
a response system on the trajectory of a drive system, evolving on a UPO. The
time-evolutions of the drive and the response systems are still described by
Eq. (1), but we consider the special case that the drive system is on a UPO
embedded in the chaotic attractor.

Without loss of generality, we focus on the case of the Rössler system,
obeying Eq. (13), with the same parameters used before, $a=b=0.2$ and $c=9$.
First, we find the trajectory of an unstable periodic orbit of the Rössler
system by using the MatCont toolbox [32] for MATLAB. We used the period-1
trajectory as the drive signal, as this is expected to be the most difficult
to synchronize [33]. The reservoir is trained on the UPO signal using 100
nodes, $\alpha=9.7(10^{-4}),\rho=0.9$, and input weights $\pmb{w}_{\text{in}}$
randomly drawn from a uniform distribution between $-0.5$ and $0.5$; the
reservoir then produces an estimate of the unavailable $y$ state from
knowledge of the available $x$ state of the UPO signal.

(A)
---
(B)
Figure 8: Time trajectories of the $x$ components of the drive and the
response systems when the drive Rössler system evolves on an unstable periodic
orbit. In (A), the two systems are coupled with $x_{D}\rightarrow x_{R}$ with
the coupling strength $\kappa=0.3$. In (B), the two systems are coupled with
$x_{D}\rightarrow x_{R}$ and $\hat{y}_{D}\rightarrow y_{R}$ with the same
$\kappa=0.3$. The estimation $\hat{y}_{D}$ is done through a reservoir
computer, as explained in the text.

Figure 8 shows that by using the RC we are successful in achieving
synchronization on the period-1 UPO with a coupling strength ($\kappa=0.3$)
for which the drive-response system could not synchronize only through
coupling in the $x$ component.

###  II.4 Noise

In order to compare the results of the RC fairly with ideal knowledge, we
consider the situation in which there may be measurement noise present in data
acquisition. We proceed under the assumption that all the information obtained
from the drive is noise corrupted. This holds true both during the training of
the RC as well as implementing it in the control configuration. Following Ref
[34], we set the magnitude of noise present in the training and control
configurations to be equal. Our noise corrupted signals are then,

$\tilde{\pmb{x}}_{D}(t)=\pmb{x}_{D}(t)+\epsilon\sqrt{\frac{\Delta
t}{T}}\pmb{\zeta}(t),$ (14)

$\tilde{\pmb{x}}_{D}(t)=[\tilde{x}_{D}(t),\tilde{y}_{D}(t),\tilde{z}_{D}(t)]$,
where $\epsilon$ is the magnitude of noise, $T$ is the approximate period of
oscillation, and $\pmb{\zeta}(t)$ is a vector with the same dimension of
$\pmb{x}_{D}(t)$ composed of elements randomly drawn from a standard normal
distribution.

We consider the same case of the Rössler system described in Sec. II.2 and
adopt the same trained RC used in Fig. 7 with coupling strength $\kappa=0.2$.
In Fig. 9 we plot $\bar{E}$ as we vary the magnitude of the noise $\epsilon$.
In (A) we compare the case of coupling $\tilde{y}_{D}\rightarrow y_{R}$ (noise
corrupted) with the case $\hat{y}_{D}\rightarrow y_{R}$ (noise corrupted+RC).
In (B) we compare the case of coupling $\tilde{y}_{D}\rightarrow y_{R}$ and
$\tilde{x}_{D}\rightarrow x_{R}$ (noise corrupted) with the case
$\hat{y}_{D}\rightarrow y_{R}$ and $\tilde{x}_{D}\rightarrow x_{R}$ (noise
corrupted+RC).

We see that for all values of $\epsilon$ in the plotted range, the error is
higher in the case in which $\hat{y}_{D}$ is reconstructed using an RC
observer than in the case in which the signal is actually available at the
receiver but corrupted with noise. The presence of noise results in a
deterioration of the RC performance but overall our proposed strategy based on
RC observation appears to be quite robust to the presence of noise, even for
large values of $\epsilon$. This is true for both the coupling schemes shown
in Fig. 9(A) and (B).

(A)
---
(B)
Figure 9: Rössler system. We plot the synchronization error as a function of
the magnitude of noise, $\epsilon$, as seen in Eq. (14). (A) compares the case
of coupling $\tilde{y}_{D}\rightarrow y_{R}$ (noise corrupted) with the case
$\hat{y}_{D}\rightarrow y_{R}$ (noise corrupted+RC). (B) compares the case of
coupling $\tilde{y}_{D}\rightarrow y_{R}$ and $\tilde{x}_{D}\rightarrow x_{R}$
(noise corrupted) with the case $\hat{y}_{D}\rightarrow y_{R}$ and
$\tilde{x}_{D}\rightarrow x_{R}$ (noise corrupted+RC).

##  III Conclusions

In this paper, we used reservoir computing as an observer within a control
loop to estimate the unmeasurable states of a system from its measurable
states. We consider two identical chaotic systems: a drive system
undirectionally coupled to a response system. We have successfully applied
this approach to reconstruct unavailable states of the drive system at the
response system and to isochronally synchronize the response to the drive.

In both the cases of the Chen system and of the Rössler system, we have shown
that usage of the reservoir observer allows us to achieve synchronization for
lower values of the coupling strength than it would be possible otherwise.
This is particularly relevant to practical situations in which physical
constraints may affect our ability to freely choose the strength of the
coupling. We showed successful implementation of our proposed scheme to
control the state of the response system on an unstable periodic orbit
embedded within the chaotic attractor. The simulations with noisy measured
data also demonstrated that the reservoir computer can very well estimate the
state trajectory even in the presence of measurement noise.

##  Appendix A Numerical integration

Here, we discuss the numerical method used to integrate the continuous-time
dynamics of the drive-response system when they are coupled through the
discrete-time reservoir computer.

After the training and testing phases are done, a random point from the
testing trajectory is selected to be set as the initial condition for the
drive system. We denote this by $\pmb{x}_{0}=\pmb{x}(t_{s})$ where $t_{s}$ is
the randomly chosen time-point from the testing signal. The corresponding
$\pmb{r}(t_{s})$ from the testing phase is set as the initial condition for
Eq. (2). A different random point on the training signal is chosen as the
initial condition for the response system.

Algorithm 1 The modified 4th order Runge-Kutta

1:$\tilde{\pmb{F}}(t,\pmb{y},\hat{y}_{D}),\pmb{y}(t_{0}),\pmb{r}(t_{0}),t,h,K,\alpha,A,\pmb{w}_{\text{in}},\pmb{w}_{\text{out}}$

2:$\pmb{y},\pmb{r}$

3:$k\leftarrow K$

4:for $j=0,\ldots,n-1$ do

5: $\pmb{y}_{j}\leftarrow\pmb{y}(t_{j})$

6: $\hat{y}_{D}\leftarrow\pmb{r}(t_{j})^{\top}\pmb{w}_{\text{out}}$

7: $k_{1}\leftarrow\tilde{\pmb{F}}(t_{j},\pmb{y}_{j},\hat{y}_{D})$

8:
$k_{2}\leftarrow\tilde{\pmb{F}}(t_{j}+\frac{h}{2},\pmb{y}_{j}+\frac{h}{2}k_{1},\hat{y}_{D})$

9:
$k_{3}\leftarrow\tilde{\pmb{F}}(t_{j}+\frac{h}{2},\pmb{y}_{j}+\frac{h}{2}k_{2},\hat{y}_{D})$

10: $k_{4}\leftarrow\tilde{\pmb{F}}(t_{j}+h,\pmb{y}_{j}+hk_{1},\hat{y}_{D})$

11:
$\pmb{y}(t_{j+1})\leftarrow\pmb{y}_{j}+\frac{h}{6}(k_{1}+2k_{2}+2k_{3}+k_{4})$

12: if $k=K$ then

13:
$\pmb{r}(t_{j+1})\leftarrow(1-\alpha)\pmb{r}(t_{j})+\alpha\tanh(A\pmb{r}(t_{j})+\allowbreak{\qquad\qquad
x_{D}}(t_{j})\pmb{w}_{\text{in}})$

14: $k\leftarrow 1$

15: else

16: $\pmb{r}(t_{j+1})\leftarrow\pmb{r}(t_{j})$

17: $k\leftarrow k+1$

18: end if

19:end for

20:return $\pmb{y},\pmb{r}$

The integration of Eq. (1) is done using the 4th-order fixed step Runge-Kutta
method. We assume the $x_{D}(t)$ component of the drive signal is known, and
we try to estimate $\hat{y}_{D}(t)$. We use the compact notation for Eq. (1)
as $\pmb{y}(t)=\tilde{\pmb{F}}(t,\pmb{y}(t),\hat{y}_{D}(t))$. The integration
interval is $0\leq t\leq T_{f}$ with time points $t_{j}$, $j=0,\ldots,n$
sampled at a fixed sampling time $h$. For simplicity, we assume $\Delta t$,
the time step of the RC in Eq. (2), to be an integer multiple of $h$, i.e.,
$K=\Delta t/h$ where $K\geq 1$ is an integer. For large enough $K$, the
integration procedure using the Runge-Kutta method has an acceptable numerical
error. For our simulations in this paper, we set $K=1$, since $h=\Delta
t=0.001$ showed to provide good accuracy for the integration process. See
Algorithm 1 for the pseudo-code of the modified 4th order Runge-Kutta for the
case that $x_{D}(t)$ of the drive system is known and the goal is to estimate
$\hat{y}_{D}(t)$ of the drive system.

## Acknowledgement

This work was partly funded by NIH Grant No. 1R21EB028489-01A1 and by the
Naval Research Lab’s Basic Research Program.

## Data Availability

The data that support the findings of this study are available within the
article.

## References

  * [1] I. Aranson and N. Rul’kov, “Nontrivial structure of synchronization zones in multidimensional systems,” Physics Letters A, vol. 139, no. 8, pp. 375–378, 1989.
  * [2] A. S. Pikovskii, “Synchronization and stochastization of array of self-excited oscillators by external noise,” Radiophysics and Quantum Electronics, vol. 27, no. 5, pp. 390–395, 1984.
  * [3] H. Fujisaka and T. Yamada, “Stability theory of synchronized motion in coupled-oscillator systems,” Prog. Theor. Phys., vol. 69, no. 1, pp. 32–47, 1983.
  * [4] V. Afraimovich, N. Verichev, and M. Rabinovich, “Stochastic synchronization of oscillations in dissipative systems,” Radiofizika, vol. 29, no. 9, pp. 1050–1060, 1986.
  * [5] L. M. Pecora and T. L. Carroll, “Synchronization in chaotic systems,” Physical review letters, vol. 64, no. 8, p. 821, 1990.
  * [6] L. M. Pecora and T. L. Carroll, “Master stability functions for synchronized coupled systems,” Physical review letters, vol. 80, no. 10, p. 2109, 1998\.
  * [7] A. Nazerian, S. Panahi, I. Leifer, D. Phillips, H. A. Makse, and F. Sorrentino, “Matryoshka and disjoint cluster synchronization of networks,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 32, no. 4, p. 041101, 2022.
  * [8] L. M. Pecora and T. L. Carroll, “Synchronization of chaotic systems,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 25, no. 9, p. 097611, 2015.
  * [9] P. A. Iglesias and B. P. Ingalls, Control theory and systems biology.  MIT press, 2010.
  * [10] D. Del Vecchio, A. J. Dy, and Y. Qian, “Control theory meets synthetic biology,” Journal of The Royal Society Interface, vol. 13, no. 120, p. 20160380, 2016.
  * [11] Y.-Y. Liu and A.-L. Barabási, “Control principles of complex systems,” Reviews of Modern Physics, vol. 88, no. 3, p. 035006, 2016.
  * [12] I. Klickstein, A. Shirin, and F. Sorrentino, “Energy scaling of targeted optimal control of complex networks,” Nature Communications, vol. 8, 2017\.
  * [13] H. Jaeger, “The “echo state” approach to analysing and training recurrent neural networks-with an erratum note,” Bonn, Germany: German National Research Center for Information Technology GMD Technical Report, vol. 148, no. 34, p. 13, 2001.
  * [14] W. Maass, T. Natschläger, and H. Markram, “Real-time computing without stable states: A new framework for neural computation based on perturbations,” Neural computation, vol. 14, no. 11, pp. 2531–2560, 2002\.
  * [15] J. Pathak, Z. Lu, B. R. Hunt, M. Girvan, and E. Ott, “Using machine learning to replicate chaotic attractors and calculate lyapunov exponents from data,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 27, no. 12, p. 121102, 2017.
  * [16] Z. Lu, J. Pathak, B. Hunt, M. Girvan, R. Brockett, and E. Ott, “Reservoir observers: Model-free inference of unmeasured variables in chaotic systems,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 27, no. 4, p. 041102, 2017.
  * [17] T. L. Carroll and J. D. Hart, “Time shifts to reduce the size of reservoir computers,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 32, no. 8, p. 083122, 2022.
  * [18] J. D. Hart, F. Sorrentino, and T. L. Carroll, “Time-shift selection for reservoir computing using a rank-revealing qr algorithm,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 33, no. 4, 2023.
  * [19] A. Cunillera, M. C. Soriano, and I. Fischer, “Cross-predicting the dynamics of an optically injected single-mode semiconductor laser using reservoir computing,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 29, no. 11, p. 113113, 2019.
  * [20] L.-W. Kong, Y. Weng, B. Glaz, M. Haile, and Y.-C. Lai, “Reservoir computing as digital twins for nonlinear dynamical systems,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 33, no. 3, 2023.
  * [21] R. S. Zimmermann and U. Parlitz, “Observing spatio-temporal dynamics of excitable media using reservoir computing,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 28, no. 4, p. 043118, 2018.
  * [22] S. Strogatz, Sync: The Emerging Science of Spontaneous Order.  Hyperion. New York, 2003.
  * [23] P. Antonik, M. Gulina, J. Pauwels, and S. Massar, “Using a reservoir computer to learn chaotic attractors, with applications to chaos synchronization and cryptography,” Physical Review E, vol. 98, no. 1, p. 012215, 2018.
  * [24] T. Weng, H. Yang, C. Gu, J. Zhang, and M. Small, “Synchronization of chaotic systems and their machine-learning models,” Physical Review E, vol. 99, no. 4, p. 042203, 2019.
  * [25] W. Hu, Y. Zhang, R. Ma, Q. Dai, and J. Yang, “Synchronization between two linearly coupled reservoir computers,” Chaos, Solitons & Fractals, vol. 157, p. 111882, 2022.
  * [26] T. Weng, X. Chen, Z. Ren, H. Yang, J. Zhang, and M. Small, “Synchronization of machine learning oscillators in complex networks,” Information Sciences, vol. 630, pp. 74–81, 2023.
  * [27] J. D. Hart, “Estimating the master stability function from the time series of one oscillator via machine learning,” arXiv preprint arXiv:2304.13125, 2023\.
  * [28] D. Canaday, A. Pomerance, and D. J. Gauthier, “Model-free control of dynamical systems with deep reservoir computing,” Journal of Physics: Complexity, vol. 2, no. 3, p. 035025, 2021.
  * [29] Q. Zhu, H. Ma, and W. Lin, “Detecting unstable periodic orbits based only on time series: When adaptive delayed feedback control meets reservoir computing,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 29, no. 9, p. 093125, 2019.
  * [30] L. Pecora and T. Carroll, “Master stability functions for synchronized coupled systems,” Phys. Rev. Lett., vol. 80, pp. 2109–2112, 1998.
  * [31] L. Huang, Q. Chen, Y.-C. Lai, and L. M. Pecora, “Generic behavior of master-stability functions in coupled nonlinear dynamical systems,” Phys. Rev. E, vol. 80, p. 036204, 2009.
  * [32] A. Dhooge, W. Govaerts, and Y. A. Kuznetsov, “Matcont: A matlab package for numerical bifurcation analysis of odes,” ACM Trans. Math. Softw., vol. 29, p. 141–164, jun 2003.
  * [33] B. R. Hunt and E. Ott, “Optimal periodic orbits of chaotic systems,” Physical review letters, vol. 76, no. 13, p. 2254, 1996.
  * [34] C. Nathe, C. Pappu, N. A. Mecholsky, J. Hart, T. Carroll, and F. Sorrentino, “Reservoir computing with noise,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 33, 04 2023.  041101\.


* Foss-Feig  _et al._ [2010b] M. Foss-Feig, M. Hermele, V. Gurarie, and A. M. Rey, Heavy fermions in an optical lattice, Physical Review A 82, 053624 (2010b).
  * Riegger  _et al._ [2018] L. Riegger, N. Darkwah Oppong, M. Höfer, D. Fernandes, I. Bloch, and S. Fölling, Localized Magnetic Moments with Tunable Spin Exchange in a Gas of Ultracold Fermions, Physical Review Letters 120, 143601 (2018).
  * Darkwah Oppong  _et al._ [2019] N. Darkwah Oppong, L. Riegger, O. Bettermann, M. Höfer, J. Levinsen, M. Parish, I. Bloch, and S. Fölling, Observation of Coherent Multiorbital Polarons in a Two-Dimensional Fermi Gas, Physical Review Letters 122, 193604 (2019).
  * Darkwah Oppong  _et al._ [2022] N. Darkwah Oppong, G. Pasqualetti, O. Bettermann, P. Zechmann, M. Knap, I. Bloch, and S. Fölling, Probing Transport and Slow Relaxation in the Mass-Imbalanced Fermi-Hubbard Model, Physical Review X 12, 031026 (2022).
  * Zhang  _et al._ [2014] X. Zhang, M. Bishof, S. L. Bromley, C. V. Kraus, M. S. Safronova, P. Zoller, A. M. Rey, and J. Ye, Spectroscopic observation of SU( N )-symmetric interactions in Sr orbital magnetism, Science 345, 1467 (2014).
  * Duan  _et al._ [2003] L.-M. Duan, E. Demler, and M. D. Lukin, Controlling Spin Exchange Interactions of Ultracold Atoms in Optical Lattices, Physical Review Letters 91, 090402 (2003).
  * García-Ripoll  _et al._ [2004] J. J. García-Ripoll, M. A. Martin-Delgado, and J. I. Cirac, Implementation of Spin Hamiltonians in Optical Lattices, Physical Review Letters 93, 250405 (2004).
  * Köhler  _et al._ [2006] T. Köhler, K. Góral, and P. S. Julienne, Production of cold molecules via magnetically tunable Feshbach resonances, Reviews of Modern Physics 78, 1311 (2006).
  * Bloch  _et al._ [2008] I. Bloch, J. Dalibard, and W. Zwerger, Many-body physics with ultracold gases, Reviews of Modern Physics 80, 885 (2008).
  * Sonderhouse  _et al._ [2020] L. Sonderhouse, C. Sanner, R. B. Hutson, A. Goban, T. Bilitewski, L. Yan, W. R. Milner, A. M. Rey, and J. Ye, Thermodynamics of a deeply degenerate SU(N)-symmetric Fermi gas, Nature Physics 16, 1216 (2020).
  * Milner  _et al._ [2023] W. R. Milner, L. Yan, R. B. Hutson, C. Sanner, and J. Ye, High-fidelity imaging of a band insulator in a three-dimensional optical lattice clock (2023), arXiv:2301.03343 [cond-mat, physics:physics].
  * Cheuk  _et al._ [2015] L. W. Cheuk, M. A. Nichols, M. Okan, T. Gersdorf, V. V. Ramasesh, W. S. Bakr, T. Lompe, and M. W. Zwierlein, Quantum-Gas Microscope for Fermionic Atoms, Physical Review Letters 114, 193001 (2015).
  * Parsons  _et al._ [2015] M. F. Parsons, F. Huber, A. Mazurenko, C. S. Chiu, W. Setiawan, K. Wooley-Brown, S. Blatt, and M. Greiner, Site-Resolved Imaging of Fermionic Li 6 in an Optical Lattice, Physical Review Letters 114, 213002 (2015).
  * Sherson  _et al._ [2010] J. F. Sherson, C. Weitenberg, M. Endres, M. Cheneau, I. Bloch, and S. Kuhr, Single-atom-resolved fluorescence imaging of an atomic Mott insulator, Nature 467, 68 (2010).
  * Haller  _et al._ [2015] E. Haller, J. Hudson, A. Kelly, D. A. Cotta, B. Peaudecerf, G. D. Bruce, and S. Kuhr, Single-atom imaging of fermions in a quantum-gas microscope, Nature Physics 11, 738 (2015).
  * Bohrdt  _et al._ [2018] A. Bohrdt, D. Greif, E. Demler, M. Knap, and F. Grusdt, Angle-resolved photoemission spectroscopy with quantum gas microscopes, Phys. Rev. B 97, 125117 (2018).
  * Greiner  _et al._ [2001] M. Greiner, I. Bloch, O. Mandel, T. W. Hänsch, and T. Esslinger, Exploring phase coherence in a 2d lattice of bose-einstein condensates, Phys. Rev. Lett. 87, 160405 (2001).
  * Albash and Lidar [2018] T. Albash and D. A. Lidar, Adiabatic quantum computation, Reviews of Modern Physics 90, 015002 (2018).
  * Trebst  _et al._ [2006] S. Trebst, U. Schollwöck, M. Troyer, and P. Zoller, $d$-wave resonating valence bond states of fermionic atoms in optical lattices, Phys. Rev. Lett. 96, 250402 (2006).
  * Sørensen  _et al._ [2010] A. S. Sørensen, E. Altman, M. Gullans, J. V. Porto, M. D. Lukin, and E. Demler, Adiabatic preparation of many-body states in optical lattices, Physical Review A 81, 061603 (2010).
  * Anderson [1967] P. W. Anderson, Infrared Catastrophe in Fermi Gases with Local Scattering Potentials, Physical Review Letters 18, 1049 (1967).
  * Gebert  _et al._ [2014] M. Gebert, H. Küttler, and P. Müller, Anderson’s Orthogonality Catastrophe, Communications in Mathematical Physics 329, 979 (2014).
  * Abramowitz and Stegun [1964] M. Abramowitz and I. A. Stegun,  _Handbook of Mathematical Functions_ (National Bureau of Standards, 1964).
  * Busch and Penson [1987] U. Busch and K. A. Penson, Tight-binding electrons on open chains: Density distribution and correlations, Physical Review B 36, 9271 (1987).
  * Noschese  _et al._ [2013] S. Noschese, L. Pasquini, and L. Reichel, Tridiagonal Toeplitz matrices: properties and novel applications, Numerical Linear Algebra with Applications 20, 302 (2013).


${\mathbf{C}}_{{\mathbf{x}}}=\frac{1}{N}\boldsymbol{\Sigma}^{2}$ as computed
by PCA. These principal components were also used as the relevant information
to be reconstructed, i.e., ${\mathbf{T}}={\mathbf{I}}$. The sensor noise was
assumed to be isotropic with covariance
${\mathbf{C}}_{{\mathbf{n}}_{\mathscr{S}}}=\sigma^{2}{\mathbf{I}}_{d_{\mathscr{S}}}$
with $\sigma=0.02$. We tried many other values of $\sigma$, yielding different
sensor locations, none of which could be used for nonlinear reconstruction.
The ones we show at $\sigma=0.02$ are representative.

###  A.4. Maximum Likelihood D-Optimal Selection

We used the maximum likelihood D-optimal selection technique based on convex
relaxation found in [26] in order to choose sensors to try to reconstruct only
the $3$rd and $4$th principal components of the shock-mixing layer snapshots.
That is, if
${\mathbf{U}}=\begin{bmatrix}{\mathbf{u}}_{1}&{\mathbf{u}}_{2}&\cdots\end{bmatrix}$
is the matrix of principal components, we model the state as a linear
combination of ${\mathbf{u}}_{3}$ and ${\mathbf{u}}_{4}$ together with
isotropic Gaussian noise. We try to find the sensors so that the correct
coefficients on ${\mathbf{u}}_{3}$ and ${\mathbf{u}}_{4}$ can be recovered
with high confidence from the measurements. The rationale for doing so is the
fact that these two components are sufficient to nonlinearly reconstruct the
state of the system if they can be measured. As in Section A.3 above, we use a
direct Python translation of a MATLAB code written by S. Joshi and S. Boyd,
which may be found at
https://web.stanford.edu/~boyd/papers/matlab/sensor_selection/.

###  A.5. Pivoted QR Factorization

For the pivoted QR factorization method [18, 8] applied to the shock-mixing
layer flow, we represent the state approximately as a linear combination of
the leading three principal components. Scipy’s implementation of pivoted QR
factorization found at
https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.qr.html was
used to select among the $2210$ allowable sensors those that allow robust
reconstruction of these first three principal components. We also tried
representing the state using more principal components and taking the first
three sensor locations chosen via pivoted QR factorization. As with the case
when only three principal components are used, these sensors do not enable
nonlinear reconstruction of the state.

###  A.6. Secant-Based Detectable Differences

The secant-based detectable difference method was implemented using the
accelerated greedy algorithm of M. Minoux [35] to optimize the objective
computed over all secants between points in the training data set consisting
of $N=750$ snapshots of the shock-mixing layer velocity field. We select among
the $2210$ sensor measurements of $u$ and $v$ velocity components on a grid of
$1105$ spatial locations taken directly from the shock-mixing layer snapshot
data. The target variables were chosen to be the leading two Isomap
coordinates
${\mathbf{g}}({\mathbf{x}})=\left(\phi_{1}({\mathbf{x}}),\phi_{2}({\mathbf{x}})\right)$,
which reveal the phase angle $\theta$. The greedy algorithm first reveals the
two sensor locations marked by green stars and then the black star in Figure 1
over the range of $0.02\leq\gamma\leq 0.06$, which can be used to reveal the
exact phase of the system. Choosing smaller values of $\gamma$ produce
different sensors that can also be used to reveal the phase, but with reduced
robustness to measurement perturbations. Gaussian process regression [46] was
used to reconstruct the leading $100$ principal components of the flowfields
from the sensor measurements. We used scikit learn’s implementation which can
be found at https://scikit-
learn.org/stable/modules/generated/sklearn.gaussian_process.GaussianProcessRegressor.html
together with a Matérn and white noise kernel whose parameters were optimized
during the fit.

For the torus example, the relevant information we wish to reconstruct are the
leading $100$ Isomap eigen-coordinates
${\mathbf{g}}({\mathbf{x}})=\left(\phi_{1}({\mathbf{x}}),\ldots,\phi_{100}({\mathbf{x}})\right)$
computed from $2000$ points sampled from the torus according to Eq. 11. The
objective function was evaluated using secants between $\\#(\mathcal{B})=100$
randomly sampled base points and the original set of $N=2000$ points. The
correct three coordinates $\phi_{1},\phi_{2},\phi_{7}$ are selected from among
the first $100$ consistently across a wide range of measurement separation
values $0.05\leq\gamma\leq 3.0$. We note that these values vary slightly with
the selected base points and these particular values hold only for one
instance.

###  A.7. Secant-Based Amplification Tolerance

Like the secant-based detectable difference method described above, the
secant-based amplification tolerance method was implemented using the same
data, secant vectors, and target variables with the accelerated greedy
algorithm. A bisection search was used to find the smallest Lipschitz constant
$L=1868$ to within a tolerance of $1$ for which the algorithm selects three
sensors on the shock-mixing layer flow. Three (different) sensors that
correctly reveal the state of the flow are selected by this algorithm over a
range $1868\leq L\leq 47624$, above which only two sensors that cannot reveal
the state are selected. We also find that with $L=129$, the minimum possible
number of sensors exceeds $\\#(\mathscr{S}_{K})/(1+\ln{\kappa})=3.18>3$.
Therefore, the minimum possible reconstruction Lipschitz constant using three
sensors that one might find by an exhaustive combinatorial search must be
greater than $129$. We admit that this is likely a rather pessimistic bound,
but we cannot check it as there are $\binom{2210}{3}\approx 1.8\times 10^{9}$
possible choices for three sensors in this problem.

When applied to select from among the leading $100$ Isomap eigen-coordinates
on the torus example with the same setup as the secant-based detectable
differences method, the amplification tolerance method selects the appropriate
collection $\phi_{1},\phi_{2},\phi_{7}$ over the range $7.1\leq L\leq 25$. We
note that these value vary slightly with the selected base points and these
particular values hold only for one instance.

##  Appendix B Submodularity of Objectives

We will need the definition of a modular function given below.

######  Definition B.1 (Modular Function).

Denote the set of all subsets of $\mathscr{M}$ by $2^{\mathscr{M}}$. A real-
valued function of the subsets $f:2^{\mathscr{M}}\to\mathbb{R}$ is called
“modular” when it can be written as a sum

(101) $f(\mathscr{S})=\sum_{j\in\mathscr{S}}a_{j}$

of constants $a_{j}$, $j\in\mathscr{M}$.

The key ingredient needed to prove submodularity for the objectives described
in Section 4 is the following lemma.

######  Lemma B.2 (Concave Composed with Modular is Submodular).

Let $h:\mathbb{R}\to\mathbb{R}$ be a concave function and let
$a:2^{\mathscr{M}}\to\mathbb{R}$ defined by

(102) $a(\mathscr{S})=\sum_{j\in\mathscr{S}}a_{j}$

be a modular function (Def. B.1) of subsets $\mathscr{S}\subseteq\mathscr{M}$
with $a_{j}\geq 0$ for all $j\in\mathscr{M}$. Then the function
$f:2^{\mathscr{M}}\to\mathbb{R}$ defined by

(103) $f(\mathscr{S})=h(a(\mathscr{S}))$

is submodular.

###### Proof.

Suppose that
$\mathscr{S}\subseteq\mathscr{S}^{\prime}\subseteq{\mathscr{M}}\setminus\\{j\\}$.
By concavity of $h$ we have

(104)
$h_{\alpha}=h((1-\alpha)a(\mathscr{S})+\alpha(a(\mathscr{S}^{\prime})+a_{j}))\geq(1-\alpha)h_{0}+\alpha
h_{1}$

for every $\alpha\in[0,1]$, where we note that $h_{0}=f(\mathscr{S})$ and
$h_{1}=f(\mathscr{S}^{\prime}\cup\\{j\\})$.

Since $\\{a_{l}\\}$ are non-negative we have $a(\mathscr{S})\leq
a(\mathscr{S})+a_{j}\leq a(\mathscr{S}^{\prime})+a_{j}$ and
$a(\mathscr{S})\leq a(\mathscr{S}^{\prime})\leq
a(\mathscr{S}^{\prime})+a_{j}$. We can therefore find

(105)
$\alpha_{1}=\frac{a_{j}}{a(\mathscr{S}^{\prime})+a_{j}-a(\mathscr{S})},\quad\alpha_{2}=\frac{a(\mathscr{S}^{\prime})-a(\mathscr{S})}{a(\mathscr{S}^{\prime})+a_{j}-a(\mathscr{S})}$

so that $h_{\alpha_{1}}=f(\mathscr{S}\cup\\{j\\})$ and
$h_{\alpha_{2}}=f(\mathscr{S}^{\prime})$. Note that $\alpha_{1}+\alpha_{2}=1$.

We now use Eq. 104 at $\alpha_{1}$ and $\alpha_{2}$ to bound the increments of
$f$:

(106)
$f(\mathscr{S}\cup\\{j\\})-f(\mathscr{S})=h_{\alpha_{1}}-h_{0}\geq\alpha_{1}(h_{1}-h_{0}),$
(107)
$f(\mathscr{S}^{\prime}\cup\\{j\\})-f(\mathscr{S}^{\prime})=h_{1}-h_{\alpha_{2}}\leq(1-\alpha_{2})(h_{1}-h_{0})$

Combining the bounds Eq. 106 and Eq. 107 on the increments using
$1-\alpha_{2}=\alpha_{1}$ we conclude that $f$ is submodular

(108) $f(\mathscr{S}\cup\\{j\\})-f(\mathscr{S})\geq
f(\mathscr{S}^{\prime}\cup\\{j\\})-f(\mathscr{S}^{\prime}).$

∎

Using Lemma B.2 it suffices to observe that each of the objectives described
in Section 4 can be written as the composition of a concave function and a
modular function. We carry this out below in addition to proving normalization
and monotonicity for these objectives.

######  Lemma B.3 (Detectable Difference Objective is Submodular).

Suppose that the target variables ${\mathbf{g}}$ and measurements
${\mathbf{m}}_{j}$, $j\in\mathscr{M}$ are measurable functions. If $\mu$ and
$\nu$ are measures on $\mathcal{X}$, then the function defined by

(109)
$f(\mathscr{S})=\int_{\begin{subarray}{c}({\mathbf{x}},{\mathbf{x}}^{\prime})\in\mathcal{X}\times\mathcal{X}:\\\
\|{\mathbf{g}}({\mathbf{x}})-{\mathbf{g}}({\mathbf{x}}^{\prime})\|_{2}\geq\varepsilon\end{subarray}}w_{\gamma,{\mathbf{x}},{\mathbf{x}}^{\prime}}(\mathscr{S})\|{\mathbf{g}}({\mathbf{x}})-{\mathbf{g}}({\mathbf{x}}^{\prime})\|_{2}^{2}\
d\mu({\mathbf{x}})\nu(d{\mathbf{x}}^{\prime}),$

for any $\varepsilon\geq 0$ with

(110)
$w_{\gamma,{\mathbf{x}},{\mathbf{x}}^{\prime}}(\mathscr{S})=\min\left\\{\frac{1}{\gamma^{2}}\|{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}})-{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}^{\prime})\|_{2}^{2},\
1\right\\},$

is normalized so that $f(\emptyset)=0$, monotone non-decreasing so that
$\mathscr{S}\subseteq\mathscr{S}^{\prime}\ \Rightarrow\ f(\mathscr{S})\leq
f(\mathscr{S}^{\prime})$, and submodular (Def. 4.2).

###### Proof.

Normalization is obvious. It suffices to prove that the function
$w_{{\mathbf{x}},{\mathbf{x}}^{\prime}}(\mathscr{S})$ is monotone and
submodular for any fixed ${\mathbf{x}},{\mathbf{x}}^{\prime}\in\mathcal{X}$.
For if we suppose that

(111)
$\mathscr{S}\subseteq\mathscr{S}^{\prime}\subseteq{\mathscr{M}}\setminus\\{j\\}\quad\Rightarrow\quad
w_{\gamma,{\mathbf{x}},{\mathbf{x}}^{\prime}}(\mathscr{S}\cup\\{j\\})-w_{\gamma,{\mathbf{x}},{\mathbf{x}}^{\prime}}(\mathscr{S})\\\
\geq
w_{\gamma,{\mathbf{x}},{\mathbf{x}}^{\prime}}(\mathscr{S}^{\prime}\cup\\{j\\})-w_{\gamma,{\mathbf{x}},{\mathbf{x}}^{\prime}}(\mathscr{S}^{\prime}),$

then multiplying both sides of the inequality by
$\|{\mathbf{g}}({\mathbf{x}})-{\mathbf{g}}({\mathbf{x}}^{\prime})\|_{2}^{2}$
and integrating proves that $f$ is submodular. The same argument also proves
monotonicity.

Let ${\mathbf{x}},{\mathbf{x}}^{\prime}\in\mathcal{X}$ be fixed. The squared
separation between the measurements is given by a modular (Def. B.1) sum

(112) $\mathscr{S}\ \mapsto\
\|{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}})-{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}^{\prime})\|_{2}^{2}=\sum_{j\in\mathscr{S}}\|{\mathbf{m}}_{j}({\mathbf{x}})-{\mathbf{m}}_{j}({\mathbf{x}}^{\prime})\|_{2}^{2}$

of non-negative constants
$\|{\mathbf{m}}_{j}({\mathbf{x}})-{\mathbf{m}}_{j}({\mathbf{x}}^{\prime})\|_{2}^{2}$
over each $j\in\mathscr{S}$. Since $x\mapsto\min\\{x/\gamma^{2},\ 1\\}$ is a
non-decreasing function, it follows that
$\mathscr{S}\subseteq\mathscr{S}^{\prime}\ \Rightarrow\
w_{{\mathbf{x}},{\mathbf{x}}^{\prime}}(\mathscr{S})\leq
w_{{\mathbf{x}},{\mathbf{x}}^{\prime}}(\mathscr{S}^{\prime})$, proving
monotonicity.

Submodularity of $w_{{\mathbf{x}},{\mathbf{x}}^{\prime}}(\mathscr{S})$ follows
from Lemma B.2 since $w_{{\mathbf{x}},{\mathbf{x}}^{\prime}}(\mathscr{S})$ is
the composition of a concave function $x\mapsto\min\\{x/\gamma^{2},\ 1\\}$
with the modular function in Eq. 112. ∎

######  Lemma B.4 (Lipschitz Objective is Submodular).

Suppose that the target variables ${\mathbf{g}}$ and measurements
${\mathbf{m}}_{j}$, $j\in\mathscr{M}$ are measurable functions. If $\mu$ and
$\nu$ are measures on $\mathcal{X}$, then the function defined by

(113)
$f(\mathscr{S})=\int_{\begin{subarray}{c}({\mathbf{x}},{\mathbf{x}}^{\prime})\in\mathcal{X}\times\mathcal{X}:\\\
{\mathbf{g}}({\mathbf{x}})\neq{\mathbf{g}}({\mathbf{x}}^{\prime})\end{subarray}}g_{{\mathbf{x}},{\mathbf{x}}^{\prime}}(\mathscr{S})\
d\mu({\mathbf{x}})\nu(d{\mathbf{x}}^{\prime}),$

with

(114)
$g_{{\mathbf{x}},{\mathbf{x}}^{\prime}}(\mathscr{S})=\min\left\\{\frac{\|{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}})-{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}^{\prime})\|_{2}^{2}}{\|{\mathbf{g}}({\mathbf{x}})-{\mathbf{g}}({\mathbf{x}}^{\prime})\|_{2}^{2}},\
\frac{1}{L^{2}}\right\\},$

is normalized so that $f(\emptyset)=0$, monotone non-decreasing so that
$\mathscr{S}\subseteq\mathscr{S}^{\prime}\ \Rightarrow\ f(\mathscr{S})\leq
f(\mathscr{S}^{\prime})$, and submodular (Def. 4.2).

###### Proof.

Normalization is obvious. It suffices to prove that the function
$g_{{\mathbf{x}},{\mathbf{x}}^{\prime}}(\mathscr{S})$ is monotone and
submodular for any fixed ${\mathbf{x}},{\mathbf{x}}^{\prime}\in\mathcal{X}$.
For if we suppose that

(115)
$\mathscr{S}\subseteq\mathscr{S}^{\prime}\subseteq{\mathscr{M}}\setminus\\{j\\}\quad\Rightarrow\quad
g_{{\mathbf{x}},{\mathbf{x}}^{\prime}}(\mathscr{S}\cup\\{j\\})-g_{{\mathbf{x}},{\mathbf{x}}^{\prime}}(\mathscr{S})\geq
g_{{\mathbf{x}},{\mathbf{x}}^{\prime}}(\mathscr{S}^{\prime}\cup\\{j\\})-g_{{\mathbf{x}},{\mathbf{x}}^{\prime}}(\mathscr{S}^{\prime}),$

then integrating both sides of the inequality proves that $f$ is submodular.
The same argument also proves monotonicity.

Let ${\mathbf{x}},{\mathbf{x}}^{\prime}\in\mathcal{X}$ be fixed. The squared
separation between the measurements is given by a modular (Def. B.1) sum

(116) $\mathscr{S}\ \mapsto\
\|{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}})-{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}^{\prime})\|_{2}^{2}=\sum_{j\in\mathscr{S}}\|{\mathbf{m}}_{j}({\mathbf{x}})-{\mathbf{m}}_{j}({\mathbf{x}}^{\prime})\|_{2}^{2}$

of non-negative constants
$\|{\mathbf{m}}_{j}({\mathbf{x}})-{\mathbf{m}}_{j}({\mathbf{x}}^{\prime})\|_{2}^{2}$
over each $j\in\mathscr{S}$. Since

(117) $x\ \mapsto\
\min\left\\{\frac{x}{\|{\mathbf{g}}({\mathbf{x}})-{\mathbf{g}}({\mathbf{x}}^{\prime})\|_{2}^{2}},\
\frac{1}{L^{2}}\right\\}$

is a non-decreasing function, it follows that
$\mathscr{S}\subseteq\mathscr{S}^{\prime}\ \Rightarrow\
g_{{\mathbf{x}},{\mathbf{x}}^{\prime}}(\mathscr{S})\leq
g_{{\mathbf{x}},{\mathbf{x}}^{\prime}}(\mathscr{S}^{\prime})$, proving
monotonicity.

Submodularity of $g_{{\mathbf{x}},{\mathbf{x}}^{\prime}}(\mathscr{S})$ follows
from Lemma B.2 since $g_{{\mathbf{x}},{\mathbf{x}}^{\prime}}(\mathscr{S})$ is
the composition of the concave function in Eq. 117 with the modular function
in Eq. 116. ∎

##  Appendix C Proofs

###### Proposition 4.5: Separation Guarantee on Underlying Set.

The result follows immediately from the triangle inequality. Let
${\mathbf{x}},{\mathbf{x}}^{\prime}\in\mathcal{X}$ and
${\mathbf{x}}_{i},{\mathbf{x}}_{j}\in\mathcal{X}_{N}$ so that
$\|{\mathbf{x}}-{\mathbf{x}}_{i}\|_{2}<\varepsilon_{0}$ and
$\|{\mathbf{x}}^{\prime}-{\mathbf{x}}_{j}\|_{2}<\varepsilon_{0}$. Then
$\varepsilon+2\varepsilon_{0}\|{\mathbf{g}}\|_{\text{lip}}\leq\|{\mathbf{g}}({\mathbf{x}})-{\mathbf{g}}({\mathbf{x}}^{\prime})\|_{2}$
implies that

(118)
$\begin{split}\varepsilon+2\varepsilon_{0}\|{\mathbf{g}}\|_{\text{lip}}&\leq\|{\mathbf{g}}({\mathbf{x}})-{\mathbf{g}}({\mathbf{x}}^{\prime})\|_{2}\\\
&\leq\|{\mathbf{g}}({\mathbf{x}})-{\mathbf{g}}({\mathbf{x}}_{i})\|_{2}+\|{\mathbf{g}}({\mathbf{x}}^{\prime})-{\mathbf{g}}({\mathbf{x}}_{j})\|_{2}+\|{\mathbf{g}}({\mathbf{x}}_{i})-{\mathbf{g}}({\mathbf{x}}_{j})\|_{2}\\\
&<\|{\mathbf{g}}({\mathbf{x}}_{i})-{\mathbf{g}}({\mathbf{x}}_{j})\|_{2}+2\varepsilon_{0}\|{\mathbf{g}}\|_{\text{lip}},\end{split}$

hence,
$\|{\mathbf{g}}({\mathbf{x}}_{i})-{\mathbf{g}}({\mathbf{x}}_{j})\|_{2}\geq\varepsilon$.
By assumption, this implies that
$\|{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}_{i})-{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}_{j})\|_{2}\geq\gamma$
and

(119)
$\begin{split}\gamma&\leq\|{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}_{i})-{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}_{j})\|_{2}\\\
&\leq\|{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}_{i})-{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}})\|_{2}+\|{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}^{\prime})-{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}_{j})\|_{2}+\|{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}})-{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}^{\prime})\|_{2}\\\
&<2\varepsilon_{0}\|{\mathbf{m}}_{\mathscr{S}}\|_{\text{lip}}+\|{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}})-{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}^{\prime})\|_{2},\end{split}$

hence,
$\|{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}})-{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}^{\prime})\|_{2}>\gamma-2\varepsilon_{0}\|{\mathbf{m}}_{\mathscr{S}}\|_{\text{lip}}$
as claimed. ∎

###### Proposition 4.6: Amplification Guarantee on Underlying Set.

The result follows immediately from the triangle inequality. Let
${\mathbf{x}},{\mathbf{x}}^{\prime}\in\mathcal{X}$ and
${\mathbf{x}}_{i},{\mathbf{x}}_{j}\in\mathcal{X}_{N}$ so that
$\|{\mathbf{x}}-{\mathbf{x}}_{i}\|_{2}<\varepsilon_{0}$ and
$\|{\mathbf{x}}^{\prime}-{\mathbf{x}}_{j}\|_{2}<\varepsilon_{0}$, then

(120)
$\begin{split}\|{\mathbf{g}}({\mathbf{x}})-{\mathbf{g}}({\mathbf{x}}^{\prime})\|_{2}&\leq\|{\mathbf{g}}({\mathbf{x}})-{\mathbf{g}}({\mathbf{x}}_{i})\|_{2}+\|{\mathbf{g}}({\mathbf{x}}^{\prime})-{\mathbf{g}}({\mathbf{x}}_{j})\|_{2}+\|{\mathbf{g}}({\mathbf{x}}_{i})-{\mathbf{g}}({\mathbf{x}}_{j})\|_{2}\\\
&<2\varepsilon_{0}\|{\mathbf{g}}\|_{\text{lip}}+L\|{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}_{i})-{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}_{j})\|_{2}\\\
&\leq
2\varepsilon_{0}\|{\mathbf{g}}\|_{\text{lip}}+L\|{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}})-{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}_{i})\|_{2}+L\|{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}^{\prime})-{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}_{j})\|_{2}\\\
&\hskip
56.9055pt+L\|{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}})-{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}^{\prime})\|_{2}\\\
&<2\varepsilon_{0}\|{\mathbf{g}}\|_{\text{lip}}+2L\varepsilon_{0}\|{\mathbf{m}}_{\mathscr{S}}\|_{\text{lip}}+L\|{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}})-{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}^{\prime})\|_{2}.\end{split}$

Gathering terms on $\varepsilon_{0}$ completes the proof. ∎

###### Proposition 6.1: Noisy Separation Guarantee.

Choose ${\mathbf{x}}_{i},{\mathbf{x}}_{j}\in\mathcal{X}_{N}$ and suppose that

(121)
$\|{\mathbf{g}}({\mathbf{x}}_{i})-{\mathbf{g}}({\mathbf{x}}_{j})\|_{2}\geq\varepsilon+2\delta_{v}.$

Then we have

(122)
$\begin{split}\|\left({\mathbf{g}}({\mathbf{x}}_{i})+{\mathbf{v}}_{i}\right)-\left({\mathbf{g}}({\mathbf{x}}_{j})+{\mathbf{v}}_{j}\right)\|_{2}&\geq\|{\mathbf{g}}({\mathbf{x}}_{i})-{\mathbf{g}}({\mathbf{x}}_{j})\|_{2}-\|{\mathbf{v}}_{i}\|-\|{\mathbf{v}}_{j}\|\\\
&\geq\|{\mathbf{g}}({\mathbf{x}}_{i})-{\mathbf{g}}({\mathbf{x}}_{j})\|_{2}-2\delta_{v}\\\
&\geq\varepsilon\end{split}$

By our assumption, this implies

(123)
$\|\left({\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}_{i})+{\mathbf{u}}_{i,\mathscr{S}}\right)-\left({\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}_{j})+{\mathbf{u}}_{j,\mathscr{S}}\right)\|_{2}\geq\gamma,$

and so we have

(124)
$\begin{split}\|{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}_{i})-{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}_{j})+{\mathbf{u}}_{j,\mathscr{S}}\|_{2}&\geq\|\left({\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}_{i})+{\mathbf{u}}_{i,\mathscr{S}}\right)-\left({\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}_{j})+{\mathbf{u}}_{j,\mathscr{S}}\right)\|_{2}-\|{\mathbf{u}}_{i,\mathscr{S}}\|-\|{\mathbf{u}}_{j,\mathscr{S}}\|\\\
&\geq\gamma-2\delta_{u}.\end{split}$

Therefore, we have established that

(125)
$\forall{\mathbf{x}}_{i},{\mathbf{x}}_{j}\in\mathcal{X}_{N}\qquad\|{\mathbf{g}}({\mathbf{x}}_{i})-{\mathbf{g}}({\mathbf{x}}_{j})\|_{2}\geq\varepsilon+2\delta_{v}\\\
\Rightarrow\quad\|{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}_{i})-{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}_{j})+{\mathbf{u}}_{j,\mathscr{S}}\|_{2}\geq\gamma-2\delta_{u}.$

The conclusion follows immediately by Proposition 4.5. ∎

###### Theorem 5.5: Down-Sampled Amplification Guarantee.

For simplicity, we will drop $L$ from the subscript on our objective since the
threshold $L$ for the Lipschitz constant remains fixed throughout the proof.
Let us begin by fixing a set $\mathscr{S}\subseteq\mathscr{M}$ and define the
random variables

(126)
$Z_{\mathscr{S}}({\mathbf{b}}_{i})=\max_{{\mathbf{x}}\in\mathcal{X}_{N}}\mathbbm{1}\big{\\{}\|{\mathbf{g}}({\mathbf{b}}_{i})-{\mathbf{g}}({\mathbf{x}})\|_{2}>L\|{\mathbf{m}}_{\mathscr{S}}({\mathbf{b}}_{i})-{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}})\|_{2}\big{\\}},$

for $i=1,\ldots,m$. If $Z_{\mathscr{S}}({\mathbf{b}}_{i})=0$ then every secant
between ${\mathbf{b}}_{i}$ and points of $\mathcal{X}_{N}$ satisfies the
desired bound on the amplification. Otherwise, there is some point
${\mathbf{x}}\in\mathcal{X}$ for which

(127)
$\|{\mathbf{g}}({\mathbf{b}}_{i})-{\mathbf{g}}({\mathbf{x}})\|_{2}>L\|{\mathbf{m}}_{\mathscr{S}}({\mathbf{b}}_{i})-{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}})\|_{2}$

and so $Z_{\mathscr{S}}({\mathbf{b}}_{i})=1$. We observe that
$Z_{\mathscr{S}}({\mathbf{b}}_{1}),\ldots,Z_{\mathscr{S}}({\mathbf{b}}_{m})$
are independent, identically distributed Bernoulli random variables whose
expectation

(128)
$\mathbb{E}[Z_{\mathscr{S}}({\mathbf{b}}_{i})]=\mu\big{(}\big{\\{}{\mathbf{x}}\in\mathcal{X}\
:\ \exists{\mathbf{x}}_{j}\in\mathcal{X}_{N}\quad\mbox{s.t.}\\\
\|{\mathbf{g}}({\mathbf{x}})-{\mathbf{g}}({\mathbf{x}}_{j})\|_{2}>L\|{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}})-{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}_{j})\|_{2}\big{\\}}\big{)}$

is the $\mu$-measure of points in $\mathcal{X}$ that are not adequately
separated from points in the $\varepsilon_{0}$-net $\mathcal{X}_{N}$ by the
measurements ${\mathbf{m}}_{\mathscr{S}}$. Suppose that for a fixed
${\mathbf{x}}\in\mathcal{X}$ we have

(129) $\|{\mathbf{g}}({\mathbf{x}})-{\mathbf{g}}({\mathbf{x}}_{j})\|_{2}\leq
L\|{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}})-{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}_{j})\|_{2}$

for every ${\mathbf{x}}_{j}\in\mathcal{X}_{N}$. By definition of
$\mathcal{X}_{N}$, for any ${\mathbf{x}}^{\prime}\in\mathcal{X}$, there is an
${\mathbf{x}}_{j}\in\mathcal{X}_{N}$ with
$\|{\mathbf{x}}^{\prime}-{\mathbf{x}}_{j}\|_{2}<\varepsilon_{0}$ and so we
have

(130)
$\begin{split}\|{\mathbf{g}}({\mathbf{x}})-{\mathbf{g}}({\mathbf{x}}^{\prime})\|_{2}&\leq\|{\mathbf{g}}({\mathbf{x}})-{\mathbf{g}}({\mathbf{x}}_{j})\|_{2}+\|{\mathbf{g}}({\mathbf{x}}_{j})-{\mathbf{g}}({\mathbf{x}}^{\prime})\|_{2}\\\
&<L\|{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}})-{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}_{j})\|_{2}+\varepsilon_{0}\|{\mathbf{g}}\|_{\text{lip}}\\\
&\leq
L\|{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}})-{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}^{\prime})\|_{2}+L\|{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}^{\prime})-{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}_{j})\|_{2}+\varepsilon_{0}\|{\mathbf{g}}\|_{\text{lip}}\\\
&<L\|{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}})-{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}^{\prime})\|_{2}+\left(\|{\mathbf{g}}\|_{\text{lip}}+L\|{\mathbf{m}}_{\mathscr{S}}\|_{\text{lip}}\right)\varepsilon_{0}.\end{split}$

It follows that $\mathbb{E}[Z_{\mathscr{S}}({\mathbf{b}}_{i})]$ is an upper
bound on the $\mu$-measure of points in $\mathcal{X}$ for which the relaxed
amplification threshold is exceeded, that is,

(131)
$\mathbb{E}[Z_{\mathscr{S}}({\mathbf{b}}_{i})]\geq\mu\big{(}\big{\\{}{\mathbf{x}}\in\mathcal{X}\
:\
\exists{\mathbf{x}}^{\prime}\in\mathcal{X}\quad\mbox{s.t.}\quad\|{\mathbf{g}}({\mathbf{x}})-{\mathbf{g}}({\mathbf{x}}^{\prime})\|_{2}\\\
\geq
L\|{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}})-{\mathbf{m}}_{\mathscr{S}}({\mathbf{x}}^{\prime})\|_{2}+\left(\|{\mathbf{g}}\|_{\text{lip}}+L\|{\mathbf{m}}_{\mathscr{S}}\|_{\text{lip}}\right)\varepsilon_{0}\big{\\}}\big{)}.$

By assumption, we have a set $\mathscr{S}\subseteq\mathscr{M}$ so that
$Z_{\mathscr{S}}({\mathbf{b}}_{i})=0$ for each $i=1,\ldots,m$. And so it
remains to bound the difference between the empirical and true expectation of
$Z_{\mathscr{S}}({\mathbf{b}}_{i})$ uniformly over every subset
$\mathscr{S}\subseteq\mathscr{M}$. For fixed $\mathscr{S}$, the one-sided
Hoeffding inequality gives

(132)
$\mathbb{P}\Big{\\{}\frac{1}{m}\sum_{i=1}^{m}\left(\mathbb{E}[Z_{\mathscr{S}}({\mathbf{b}}_{i})]-Z_{\mathscr{S}}({\mathbf{b}}_{i})\right)\geq\delta\Big{\\}}\leq
e^{-2m\delta^{2}}.$

Unfixing $\mathscr{S}$ via the union bound over all
$\mathscr{S}\subseteq\mathscr{M}$ and applying our assumption about the number
of base points $m$ yields

(133)
$\mathbb{P}\bigcup_{\mathscr{S}\subseteq\mathscr{M}}\Big{\\{}\frac{1}{m}\sum_{i=1}^{m}\left(\mathbb{E}[Z_{\mathscr{S}}({\mathbf{b}}_{i})]-Z_{\mathscr{S}}({\mathbf{b}}_{i})\right)\geq\delta\Big{\\}}\leq
e^{\\#(\mathscr{M})\ln{2}-2m\delta^{2}}\leq p.$

Since our assumed choice of $\mathscr{S}$ has
$f_{m}(\mathscr{S})=f_{m}(\mathscr{M})$ it follows that all
$Z_{\mathscr{S}}({\mathbf{b}}_{i})=0$, $i=1,\ldots,m$, hence we have

(134) $\mathbb{E}[Z_{\mathscr{S}}({\mathbf{b}}_{i})]<\delta$

with probability at least $1-p$. Combining this with Eq. 131 completes the
proof. ∎

##  Appendix D Description of the Accelerated Greedy Algorithm

Since each objective function $f$ presented in Section 4 is submodular, it is
possible to use an “accelerated greedy” (AG) algorithm to obtain the same
solution as the naive greedy algorithm with a provably minimal number of
objective function evaluations compared to a broad class of algorithms [35].
Let the increase in the objective function obtained by adding the sensor $j$
to the set $\mathscr{S}$ be called
$\Delta_{j}(\mathscr{S})=f(\mathscr{S}\cup\\{j\\})-f(\mathscr{S})$. Instead of
evaluating $\Delta_{j}(\mathscr{S}_{k-1})$ for every measurement in
$\mathscr{M}\setminus\mathscr{S}_{k-1}$, AG keeps track of an upper bound
$\hat{\Delta}_{j}\geq\Delta_{j}(\mathscr{S}_{k-1})$ on the increments for each
sensor. Since submodularity of $f$ means that the increments
$\Delta_{j}(\mathscr{S})$ can only decrease as the size of $\mathscr{S}$
increases, it is sufficient to have the maximum upper bound
$\hat{\Delta}_{j^{*}}\geq\hat{\Delta}_{j}$, $\forall
j\in\mathscr{M}\setminus\mathscr{S}_{k-1}$ be tight
$\hat{\Delta}_{j^{*}}=\Delta_{j^{*}}(\mathscr{S}_{k-1})$ in order to conclude
that $\Delta_{j^{*}}(\mathscr{S}_{k-1})$ is the largest increment. The rest of
the upper bounds on the increments can remain loose since they are smaller
than the tight maximum upper bound. The AG algorithm finds largest upper bound
$\hat{\Delta}_{j^{*}}$ and updates it so that it is tight. If
$\hat{\Delta}_{j^{*}}$ is still the greatest upper bound, then $j^{*}=j_{k}$
achieves the largest increment and is added to $\mathscr{S}_{k-1}$. Otherwise
if $\hat{\Delta}_{j^{*}}$ is no longer the largest upper bound, the new
largest upper bound is selected and process repeated until a tight maximum
upper bound is obtained.