# Zero-delay Consistent Signal Reconstruction from Streamed Multivariate Time
Series††thanks: This work was supported by the SFI Offshore Mechatronics grant
237896/O30 and the IKTPLUSS INDURB grant 270730/O70 from the Norwegian Science
Foundation, and the VALIDATE project grant 101057263 from the EU HORIZON-RIA.

Emilio Ruiz-Moreno1, Luis Miguel López-Ramos12, ,
and Baltasar Beferull-Lozano12,

1WISENET Center, Department of ICT, University of Agder, Grimstad, Norway
2Simula Metropolitan Center for Digital Engineering, Oslo, Norway

###### Abstract

Digitalizing real-world analog signals typically involves sampling in time and
discretizing in amplitude. Subsequent signal reconstructions inevitably incur
an error that depends on the amplitude resolution and the temporal density of
the acquired samples. From an implementation viewpoint, consistent signal
reconstruction methods have proven a profitable error-rate decay as the
sampling rate increases. Despite that, these results are obtained under
offline settings. Therefore, a research gap exists regarding methods for
consistent signal reconstruction from data streams. This paper presents a
method that consistently reconstructs streamed multivariate time series of
quantization intervals under a zero-delay response requirement. On the other
hand, previous work has shown that the temporal dependencies within univariate
time series can be exploited to reduce the roughness of zero-delay signal
reconstructions. This work shows that the spatiotemporal dependencies within
multivariate time series can also be exploited to achieve improved results.
Specifically, the spatiotemporal dependencies of the multivariate time series
are learned, with the assistance of a recurrent neural network, to reduce the
roughness of the signal reconstruction on average while ensuring consistency.
Our experiments show that our proposed method achieves a favorable error-rate
decay with the sampling rate compared to a similar but non-consistent
reconstruction.

###### Index Terms:

Analog-to-digital conversion, consistent signal reconstruction, stream
learning, multivariate time series, spline interpolation, sequential decision
making, recurrent neural networks.

##  I Introduction

(a) Snapshot at time $x_{t}$ of a multivariate signal being converted from
analog to digital. The acquired data samples form a streamed multivariate time
series of quantization intervals.

(b) Snapshot at time $x_{t}$ of the multivariate time series from Fig. 1(a)
being reconstructed. The consistent signal reconstruction method is
represented visually as a microprocessor. The green shaded area indicates the
currently reconstructed portion of the signal, while the red represents the
future.

Figure 1: Illustration of the considered signal acquisition-reconstruction
process spanning from time $0$ to $X$. Stationary spatial relationships are
represented as edges in a graph with black dashed lines.

Most real-world physical quantities can only be observed by sampling their
continuous-time waveforms (analog signals). If the analog signal is
bandlimited, i.e., its spectral density has bounded support, we can ensure a
reversible discretization in time by sampling at, or above, its Nyquist
sampling rate [1, 2]. This is a universal sufficient condition but not
necessary for all signals: some non-bandlimited signals can be perfectly
reconstructed after uniformly sampled at, or above, its finite rate of
innovation. Examples of such signals are streams of Dirac delta functions and
_splines_ [3]. Nonetheless, we can only observe a modified version of the
sampled analog signal in practice because of the physical limitations of
acquisition systems. For instance, acquisition systems such as analog-to-
digital converters (ADCs) introduce both a sampling and a quantization step.
The consequent discretization in amplitude causes an irreversible loss of
information, making perfect signal reconstruction no longer possible.

Even though quantization is a deterministic operation, the signal
reconstruction error-rate decay dependencies can be analyzed effectively from
a stochastic point of view. For example, the reconstruction error caused by a
uniform quantizer can be well approximated and modeled as an additive uniform
white noise source independent of the input signal when certain assumptions
are met [4, 5]. Under these assumptions, the quantization noise spectral
density is constant. Thus, if an analog signal $\psi$ is bandlimited, we can
reduce the amount of in-band quantization noise by sampling above its Nyquist
sampling rate (oversampling). Accordingly, we can reduce the reconstruction
error by digitally filtering out the quantization noise above the original
signal bandwidth. In this scenario, it is well known that low-pass filtering
reduces the mean squared error (MSE) of a reconstructed signal $f$ by a factor
proportional to the squared quantization step size $\Delta^{2}$ and inversely
proportional to the oversampling ratio $R$. Formally,

$\mathbb{E}\left[\left\|f-\psi\right\|^{2}_{2}\right]\propto\frac{\Delta^{2}}{R},$
(1)

where the expectation is taken over the amplitude of the reconstruction error,
modeled as a uniform random variable [4]. The result shown in (1) is useful
for identifying the main factors that influence the signal reconstruction
error-rate decay under an oversampling regime. However, it does not reflect
the effects of any quantization scheme or reconstruction method other than
uniform quantization and linear decoding, i.e., low-pass filtering of the
quantized signal at a cut-off frequency.

In contrast, deterministic analyses, based on deterministic bounds over the
quantization levels, show that it is possible to reconstruct quantized signals
with a squared norm error that is asymptotically proportional to the square of
the quantization step size and inversely proportional to the square of the
oversampling ratio, that is

$\left\|f-\psi\right\|^{2}_{2}=\mathcal{O}\left(\frac{\Delta^{2}}{R^{2}}\right).$
(2)

The trend in (2) can be achieved by performing signal reconstruction via
_consistent signal estimates_ , where the term “consistent” here means
coherent with all available knowledge about the original analog signal and the
acquisition system. This has been shown for bandlimited signals [6, 7], and
non-bandlimited signals with finite rate of innovation [8]. The behavior
presented in (2) suggests that when a uniformly quantized and oversampled
signal is reconstructed using consistent estimates, the reconstruction error
can be asymptotically reduced at the same rate by either increasing the
oversampling ratio or decreasing the quantization step size. This balanced
signal reconstruction error-rate decay has important practical implications.
For instance, in practice, it is more convenient to increase the oversampling
ratio because the implementation cost of the analog circuitry necessary for
halving the quantization step size is much higher than that for doubling the
oversampling ratio [9].

Notwithstanding the above, consistent signal reconstructions usually require
sophisticated methods far more computationally complex than linear decoding,
thus hindering their practical implementation. Only a few previous studies
have explored how to alleviate the computational complexity of yielding
consistent reconstructions while maintaining the asymptotic behavior
established in (2). For example, in [10], linear decoding computational
complexity is achieved by using different quantization step sizes, while the
difficulty lies in the design of the quantization scheme. As another example,
the authors of [11] design a subtractive dithering method aimed at recursively
solving an overdetermined system of linear equations set up from quantized
data described by frames111Frames provide redundant and usually non-unique
representations of vectors [12, 13].. However, the resulting signal estimates
are only guaranteed to be consistent with the most recently analyzed frame
representation in the recursion. In summary, there seems to be a research gap
regarding low complexity consistent signal reconstruction methods under
uniform quantization schemes. On the other hand, and to the best of our
knowledge, there are no consistent signal reconstruction methods designed for
sequentially reconstructing streamed data samples. Consequently, there are no
studies about their error-rate decay dependencies either. Guaranteeing a
consistent signal reconstruction from streamed data requires a different
approach than the previously discussed (offline) consistent signal
reconstruction methods. Unlike them, and as an additional research gap,
consistent signal reconstruction approaches from streaming data require the
ability to be executed online in addition to low complexity.

This paper aims to cover the aforementioned gaps. To this end, we focus on a
class of multivariate, non-bandlimited, temporally and spatially
interdependent smooth signals. This is justified because most signals of
interest are smooth, usually elapse within a finite timespan, and rarely occur
in complete isolation. Then, we consider an ADC that uniformly samples and
quantizes these multivariate signals, resulting in streamed multivariate time
series of quantization intervals, as illustrated in Fig. 1(a). From here, our
goal is to design a consistent signal reconstruction approach from streaming
data that satisfies two key requirements: signal smoothness and a zero-delay
response. These requirements arise naturally in a wide range of signal
reconstruction problems, ranging from high-speed digital to analog conversion
[14] to online trajectory planning [15, 16], among others. However, despite
the practical importance of these requirements, very few signal reconstruction
methods have been proposed to deal with them jointly. To the best of our
knowledge, our previous work [17] is the only approach that expressly
addresses such a smooth and zero-delay signal reconstruction challenge,
although only for the univariate case. Its main idea is to exploit, in a data-
driven manner, the temporal dynamics of the streaming data to reduce a given
trade-off between the squared residuals and the roughness, i.e., a derivative-
based measure of the model complexity [18], of the reconstructed signal. This
is achieved by training a recurrent neural network within a stationary
parametric policy, i.e., a time-invariant (zero-delay) decision strategy based
on tuneable parameters, on top of a recurrent neural network (RNN) [19], over
a set of representative streaming data.

Concretely, this paper presents the first approach to zero-delay consistent
signal reconstruction from a streamed multivariate time series of quantization
intervals, see Fig. 1(b). Specifically, we devise a stationary parametric RNN-
guided policy that exploits the spatiotemporal dependencies of the streaming
data to reduce the roughness of the signal reconstruction while ensuring
consistency. Moreover, we experimentally investigate its error-rate decay
behavior.

The main contributions of this paper can be enumerated as follows:

  * •

We formalize the concept of consistent signal reconstruction from streaming
data under zero-delay response requirements. As a result, our formulation
generalizes the concept of consistency beyond offline settings.

  * •

To the best of our knowledge, we devise the first trainable method for zero-
delay signal reconstruction from multivariate time series of quantization
intervals, which can also enforce consistency.

  * •

We show experimentally that the reconstruction error incurred by our proposed
method decays at the same rate by decreasing the quantization step size or
increasing the oversampling ratio. Moreover, the error-rate decay slope with
respect to the oversampling ratio doubles (in logarithmic scale) the one
obtained with a similar zero-delay smooth signal reconstruction method that
does not enforce consistency.

The rest of the paper is organized as follows. Sec. II introduces the main
notation used throughout the paper. Sec. III formalizes the notion of
consistency and presents the concept of consistent signal reconstruction from
streaming data under zero-delay requirements. Then, Sec. IV describes our
proposed method from a sequential decision-making perspective, and Sec. V
analyzes its error-rate decay behavior experimentally. Finally, Sec. VI
concludes the paper.

##  II Mathematical notation

This section introduces the mathematical notation most recurrently used
throughout the paper.

Vectors and matrices are represented in bold lowercase and bold capital
letters, respectively. Given a vector $\bm{v}=[v_{1},\dots,v_{C}]^{\top}$ its
cth component is indicated as $[\bm{v}]_{c}\triangleq v_{c}$. Similarly, given
a matrix $\bm{M}\in\mathbb{R}^{R\times C}$, the element in the rth row and cth
column is indicated as $[\bm{M}]_{r,c}$. The notation $[\bm{v}]_{i:j}$ refers
to the sliced vector $[v_{i},\dots,v_{j}]\in\mathbb{R}^{j-i+1}$. We use
Euler’s notation for the derivative operator; thus, $D^{k}_{x}$ denotes the
kth derivative over the variable $x$. Lastly, we refer to functions $f$ as
elements of a function space or explicitly allude to its domain and codomain
with the compact notation $f(x)$, depending on the context.

##  III Consistency

This section describes and models the function space comprising the signals
considered in this work and the acquisition system. Accordingly, we identify
and define the set of consistent signal estimates and introduce the related
concept of consistent signal reconstruction from streamed acquired samples.

###  III-A Function space of multivariate smooth signals

Most real-world physical processes are bounded and smooth due to energy
conservation [20]. Moreover, most of them occur within a certain finite time
lapse. Consequently, we assume any real-world process-associated measurable
quantity to be described by smooth and non-bandlimited signals.

Let us define $\mathcal{W}_{\rho}$ as the function space of real functions
defined over the domain $(0,X]\subseteq\mathbb{R}$ with $\rho-1$ absolutely
continuous derivatives and with the $\rho$th derivative square integrable.
Then, any of the aforementioned real-world signals can be accurately modeled
by a function $\psi\in\mathcal{W}_{\rho}$. In practice, $\psi$ can be
approximated, or reconstructed, from a set of $T_{*}$ uniformly sampled
measurements $\\{\psi(tX/T_{*})\\}^{T_{*}}_{t=1}$ as follows

$\displaystyle
f_{*}\triangleq\text{arg}\underset{f\in\mathcal{W}_{\rho}}{\text{ min }}$
$\displaystyle\int^{X}_{0}\left(D^{\rho}_{x}f(x)\right)^{2}\,dx$ (3a) subject
to: $\displaystyle
f\left(\frac{tX}{T_{*}}\right)=\psi\left(\frac{tX}{T_{*}}\right),\,\,\forall
t\in\mathbb{N}^{[1,T_{*}]},$ (3b)

where the metric used as objective in (3a) is usually referred to as
_roughness_ and $f_{*}$ is unique as long as $T_{*}\geq\rho$ [21, 22].
Actually, any function $\psi\in\mathcal{W}_{\rho}$ with roughness bounded by
an arbitrary constant is known to satisfy

$\underset{x\in[X/T_{*},X)}{\text{ sup
}}\left|D^{k}_{x}\psi(x)-D^{k}_{x}f_{*}(x)\right|\leq\text{const}\,\cdot\,\left(\frac{X}{T_{*}}\right)^{\rho-k},$
(4)

for all $k\in\\{0,\dots,\rho-1\\}$ [23]. That is, the approximation error, up
to the $(\rho-1)$th derivative, is upper bounded by a term inversely
proportional to the number $T_{*}$ of measurements.

Since the approximation error in (4) can be made arbitrarily small by
increasing the number of samples $T_{*}$, we restrict ourselves without loss
of generality to those signals in $\mathcal{W}_{\rho}$ whose approximation, as
in (3), is exact, i.e., the left-hand term in (4) is zero for all allowed kth
derivatives. The function space where such signals of interest belong is
denoted as $\mathcal{V}_{\rho}\subseteq\mathcal{W}_{\rho}$, which in turn is
constituted by splines, i.e., piecewise polynomial functions, of order
$2\rho-1$ with $T_{*}$ knots (joint points) at
$\mathcal{X}_{*}=\\{tX/T_{*}\\}^{T_{*}}_{t=1}$ [23].

On the other hand, more than one signal often concurrently stems from the same
physical process. Given that, one can expect spatial interdependencies among
them. For this reason, this work deals with multivariate signals of the form
$\bm{\Psi}=[\psi^{(1)},\dots,\psi^{(N)}]^{\top}$ with
$\psi^{(n)}\in\mathcal{V}_{\rho_{n}}$ for all $n\in\mathbb{N}^{[1,N]}$,
namely,
$\bm{\Psi}\in\mathcal{V}_{\rho_{1}}\times\dots\times\mathcal{V}_{\rho_{N}}$.
However, and for the sake of simplicity in the notation, we rather focus on
those multivariate signals in
$\mathcal{V}_{\rho}\times\dots\times\mathcal{V}_{\rho}=\mathcal{V}^{N}_{\rho}$,
i.e., same order $\rho$ across all $N$ dimensions, since the considered
formulation can be easily adapted to the presented general case.

###  III-B Acquisition system and irreversibility

Thanks to the reconstruction procedure presented in (3), a multivariate signal
$\bm{\Psi}\in\mathcal{V}^{N}_{\rho}$ can be reversibly sampled despite being
non-bandlimited. That is, we can recover the multivariate signal $\bm{\Psi}$
from finitely many measurements as long as they contain its set of knots,
namely, the signal values $\bm{\Psi}(x)$ sampled at $x\in\mathcal{X}_{*}$. For
example, we can measure the knots by sampling uniformly at a rate
$\nu_{*}=T_{*}/X$. In fact, any uniform sampling rate $\nu=R\nu_{*}$, with $R$
being a positive integer oversampling ratio, measures the $T_{*}$ knots
because $\mathcal{X}=\\{t\nu^{-1}\\}^{T}_{t=1}\supseteq\mathcal{X}_{*}$.

Ideally, these measurements would be acquired through an infinite amplitude-
resolution sampling mechanism, namely, an ideal sampler. However, this is
impossible in practice due to the physical limitations of measuring devices
[24]. Instead, the usual procedure is to assume an ideal sampler and then
model the resulting finite amplitude-resolution measurements as the outcome of
a quantization mapping. In this work, we assume an ideal sampler working at a
uniform sampling rate $\nu$ followed by a uniform midtread quantizer with a
quantization step size $\Delta^{(n)}$ for each of the $n$ signals.

The resulting acquired data consists of a streamed multivariate time series of
$N$ time series of length $T=RT_{*}$. We refer to every tth term of the
streamed multivariate time series as an _observation_ , $\bm{o}_{t}$. Every
tth observation is constructed from $N$ quantization intervals described by
their time stamp $x_{t}=tX/T\in\mathcal{X}$, common for all intervals of the
tth observation, their center $y^{(n)}_{t}\in\mathbb{R}$, and half
quantization step size $\epsilon^{(n)}=\Delta^{(n)}/2\in\mathbb{R}_{+}$.
Formally, $\bm{o}_{t}=[x_{t},\bm{y}_{t}^{\top},\bm{\epsilon}^{\top}]^{\top}$
where $\bm{y}_{t}=[y^{(1)}_{t},\dots,y^{(N)}_{t}]^{\top}$ and
$\bm{\epsilon}=[\epsilon^{(1)},\dots,\epsilon^{(N)}]$.

Unfortunately, the acquisition system under consideration is irreversible.
That is, any of the multivariate signals in $\mathcal{V}^{N}_{\rho}$ may not
be perfectly recovered from a finite set of measurements anymore if acquired
as described above. This drawback is two-sourced. First, the sampler works
sequentially. Thus, intermediate acquisition stages are irreversible because
not all knots, i.e., the multivariate signal values at $\mathcal{X}_{*}$, have
been measured yet. Second, the quantization step is a many-to-one mapping,
i.e., non-injective and hence, not invertible. Consequently, incurring a
reconstruction error in this context is unavoidable.

###  III-C Consistent signal estimate

(a) Intersection of hyperslabs as a Venn diagram. At time stamp $t$, the
functions within the area enclosed by the dotted black line,
$\mathcal{H}^{(n)}_{1:t}$, describe current consistent signal estimates.

(b) Comparison between the ongoing reconstructed signals using a non-
consistent (NC) and a consistent (C) method. Fig. 2(a) illustrates how both
methods produce a sequence of consistent signal estimates, namely,
$\\{\tilde{f}^{(n)}_{t}\\}^{T}_{t=1}$ and $\\{f^{(n)}_{t}\\}^{T}_{t=1}$.
However, the sequence $\\{\tilde{f}^{(n)}_{t}\\}^{T}_{t=1}$ does not lead to a
consistent signal reconstruction, while the sequence
$\\{f^{(n)}_{t}\\}^{T}_{t=1}$ does.

Figure 2: Visual representation of a consistent signal estimate and a
consistent signal reconstruction.

As stated in the previous Sec. III-B, it is generally impossible to recover a
multivariate signal $\bm{\Psi}\in\mathcal{V}^{N}_{\rho}$ from its stream of
observations. Nevertheless, we can select a relatively reduced set of
functions where the multivariate signal $\bm{\Psi}$ lies. Specifically, for
every nth time series and for any already received ith quantization interval
we can construct its associated hyperslab
$\mathcal{H}^{(n)}_{i}\subseteq\mathcal{V}_{\rho}$ as

$\mathcal{H}^{(n)}_{i}=\left\\{h\in\mathcal{V}_{\rho}:\left|h(x_{i})-y^{(n)}_{i}\right|\leq\epsilon^{(n)}\right\\}.$
(5)

Then, at time step $t$, we can state that each signal $\psi^{(n)}$ belongs to
the intersection of its so-far constructed hyperslabs. That is,
$\psi^{(n)}\in\mathcal{H}^{(n)}_{1:t}$ with
$\mathcal{H}^{(n)}_{1:t}\triangleq\bigcap^{t}_{i=1}\mathcal{H}^{(n)}_{i}$, as
illustrated in Fig. 2(a). In this context, any signal estimate described by a
function drawn from $\mathcal{H}^{(n)}_{1:t}$ is called a consistent signal
estimate at time step $t$, or simply consistent signal estimate, when clear by
context. This is because any reacquisition, i.e., processing the newly
estimated signal through the same acquisition system until the same time step
$t$, leads to the same set of already received $t$ quantization intervals, and
conventionally, we refer to any signal estimate as consistent if and only if:
i) it passes through all the (so-far) received quantization intervals, and ii)
it is described by a function in $\mathcal{V}_{\rho}$ [9].

Recall from Sec. III-A that $\mathcal{V}_{\rho}$ is constituted by splines of
order $2\rho-1$, with $\rho-1$ continuous derivatives and $T_{*}$ knots
located at time stamps in $\mathcal{X}_{*}$. For this reason, we construct
every nth signal estimate as a spline composed of $T$ pieces, or _function
sections_ , as

$f^{(n)}(x)=\begin{cases}g^{(n)}_{1}(x),&\text{if }0<x\leq x_{1}\\\
g^{(n)}_{2}(x),&\text{if }x_{1}<x\leq x_{2}\\\ \vdots\\\
g^{(n)}_{T}(x),&\text{if }x_{T-1}<x\leq x_{T}\end{cases}$ (6)

where every $t$th function section $g^{(n)}_{t}:(x_{t-1},x_{t}]\to\mathbb{R}$
is a linear combination of polynomials of the form

$g^{(n)}_{t}(x)={\bm{a}^{(n)}_{t}}^{\top}\bm{p}_{t}(x),$ (7)

with combination coefficients $\bm{a}^{(n)}_{t}\in\mathbb{R}^{2\rho}$ and
basis vector function $\bm{p}_{t}:(x_{t-1},x_{t}]\to\mathbb{R}^{2\rho}$
defined as

$\bm{p}_{t}(x)=\left[1,(x-x_{t-1}),\dots,(x-x_{t-1})^{2\rho-1}\right]^{\top}.$
(8)

Then, the remaining requirements for $f^{(n)}$ to belong to
$\mathcal{V}_{\rho}$ are that the location of its knots must cover all time
stamps in $\mathcal{X}_{*}$ and it has to be continuous up to the $(\rho-1)$th
derivative. By design, the knots of $f^{(n)}$ lie at $\mathcal{X}$ hence,
covering all time stamps in $\mathcal{X}_{*}$ since
$\mathcal{X}_{*}\subseteq\mathcal{X}$. Regarding its continuity, for $f^{(n)}$
to be $(\rho-1)$-smooth, the function sections must satisfy

$\underset{x\to x_{t-1}^{-}}{\text{lim
}}D^{k}_{x}\,g^{(n)}_{t-1}(x)=\underset{x\to x_{t-1}^{+}}{\text{lim
}}D^{k}_{x}\,g^{(n)}_{t}(x),$ (9)

for every $k\in[0,\rho-1]$ and $t\in[2,T]$. This can be ensured by imposing
the following equality constraint

$\left[\bm{a}^{(n)}_{t}\right]_{1:\rho}=\bm{e}^{(n)}_{t-1},$ (10)

for each $t\in\mathbb{N}^{[1,T]}$, where
$\bm{e}^{(n)}_{t}\in\mathbb{R}^{\rho}$ is a vector with elements given by

$[\bm{e}^{(n)}_{t}]_{i}=\frac{1}{(i-1)!}\sum^{2\rho}_{j=1}\left[\bm{a}^{(n)}_{t}\right]_{j}\,u_{t-1}^{j-i}\prod^{i-1}_{k=1}(j-k),$
(11)

where $u_{t}\triangleq x_{t}-x_{t-1}$, except $u_{1}=x_{1}$, and where
$\bm{e}^{(n)}_{0}$ determines the initial boundary conditions of the signal
estimate that the user can set ([17, Proposition 1]).

In this way, any signal estimate $f^{(n)}$ constructed as in (6), and
satisfying (9), belongs to $\mathcal{V}_{\rho}$. Additionally, and as stated
before, if at time step $t$ it passes through all the so-far received
quantization intervals, i.e., $f^{(n)}$ satisfies
$|f^{(n)}(x_{i})-y^{(n)}_{i}|\leq\epsilon^{(n)}$ for all
$i\in\\{1,\dots,t\\}$, or equivalently $f^{(n)}\in\mathcal{H}^{(n)}_{1:t}$,
then it is a consistent signal estimate at time step $t$.

###  III-D Zero-delay consistent signal reconstruction

Recall that the acquisition system modeled in Sec. III-B generates a stream of
observations. In this paper, we pursue its converse reconstruction approach;
that is, our goal is to sequentially reconstruct a signal
$\bm{\Psi}\in\mathcal{V}^{N}_{\rho}$ as its observations are received. To this
end, we rely on zero-delay multivariate signal reconstruction methods. In this
context, they yield a sequence of signal estimates in
$\mathcal{V}^{N}_{\rho}$, where the tth multivariate signal estimate
$\bm{f}_{t}=[f^{(1)}_{t},\dots,f^{(N)}_{t}]^{\top}$ is used to reconstruct the
multivariate time series between the last two received observations, i.e.,
$\bm{o}_{t-1}$ and $\bm{o}_{t}$. Moreover, every reconstruction step is
performed as soon as the tth observation is received and without knowing the
to-be-received observations222Implementing a zero-delay signal reconstruction
method also requires a reduced constant complexity per iteration to guarantee
the required execution speed and scalability [25]..

From the discussion in Sec. III-C, one may expect to obtain a consistent
signal reconstruction using a zero-delay signal reconstruction method that
generates a sequence of consistent signal estimates. That is, for each $n$th
time series every $t$th signal estimate $f^{(n)}_{t}$ belongs to
$\mathcal{H}^{(n)}_{1:t}$. However, it is important to notice that a method
producing consecutive consistent signal estimates does not necessarily achieve
a consistent reconstruction. This is because guaranteeing a consistent
reconstruction requires the sequence of signal estimates to fulfill additional
requisites. Specifically, every two consecutive consistent signal estimates,
namely $f^{(n)}_{t-1}$ and $f^{(n)}_{t}$, must ensure smoothness at the
takeover time stamp $x_{t-1}$, as shown in Fig. 2(b). Formally, a consistent
signal reconstruction is attainable from a sequence of consistent signal
estimates $\\{f^{(n)}_{t}\\}^{T}_{t=1}$ if they satisfy that
$c^{(n)}_{t}\in\mathcal{H}^{(n)}_{1:t}$, where

$c^{(n)}_{t}(x)=\begin{cases}f^{(n)}_{t-1}(x)&\text{if }0<x\leq x_{t-1},\\\
f^{(n)}_{t}(x)&\text{if }x_{t-1}<x\leq X,\end{cases}$ (12)

for all $t\in\mathbb{N}^{[2,T]}$.

Zero-delay reconstruction methods can yield a consistent signal reconstruction
by sequentially updating the spline-based signal estimate coefficients.
Concretely, if at every time step $t$, we update the remaining $(T-t)$ spline
coefficients, i.e., $\bm{a}^{(n)}_{i}$ for $i=t,\dots,T$ subject to satisfying
the corresponding continuity constraints in (10), and subject to passing
through the last received quantization interval, i.e., the tth spline
coefficients obey
$|{\bm{a}^{(n)}_{t}}^{\top}\bm{p}_{t}(x_{t})-y^{(n)}_{t}|\leq\epsilon^{(n)}$,
then the resulting signal estimate update $f^{(n)}_{t}$ is part of a sequence
satisfying (12). In this way, the overall signal reconstruction, i.e., the
reconstructed signal at the last time stamp $T$, is completely described by
the last consistent signal estimate $f^{(n)}_{T}$. In practice, updating just
the tth spline coefficients at every time step $t$ eventually leads to the
same consistent signal reconstruction.

###  III-E Reconstruction error-rate decay behavior

To the best of our knowledge, and as stated in the Introduction, the problem
of reconstructing a quantized time series using consistent signal estimates
has only been addressed when all samples are available in advance. Notice that
in that scenario, the signal estimate and the signal reconstruction coincide
because only one signal estimate is needed. However, in this work, we
explicitly account for the sequential nature of the acquisition system, so it
is conceptually more challenging since we cannot achieve a consistent
reconstruction with just a consistent signal estimate but with a sequence of
them, subject to (12), instead. This is a noticeable difference with respect
to previously discussed (offline) consistent signal reconstruction methods,
where consistency entails a profitable error-rate decay from an implementation
viewpoint [6, 7, 10, 8]. These facts open the question of whether a method for
zero-delay consistent signal reconstructions has the corresponding effects in
the reconstruction error-rate decay behavior.

We anticipate a positive answer to the posed question, but before presenting
such an outcome, we introduce our zero-delay consistent signal reconstruction
methods.

##  IV Zero-delay spline interpolation

This section introduces a novel zero-delay smooth signal reconstruction
approach based on policy training that can accommodate multivariate quantized
data and enforce consistency as defined in Sec. III.

###  IV-A A trainable multivariate approach to zero-delay spline interpolation

The problem of reconstructing a multivariate time series of quantization
intervals using a sequence of spline-based signal estimates under smoothness
and zero-delay requirements can be formalized from a dynamic programming
perspective [26].

At time step $t$, we encode the condition of the so-far reconstructed signal,
and the last received quantization interval in a vector-valued variable
referred to as the tth state $\bm{s}_{t}\in\mathcal{S}$, being $\mathcal{S}$
the state space. Every tth state is constructed as
$\bm{s}_{t}=[x_{t-1},\bm{o}^{\top}_{t},\bm{e}^{\top}_{t-1}]^{\top}$ with
$\bm{e}_{t-1}=[{\bm{e}^{(1)}_{t-1}}^{\top},\dots,{\bm{e}^{(N)}_{t-1}}^{\top}]^{\top}$
where every nth component $\bm{e}^{(n)}_{t-1}$ represents the right hand term
at $\eqref{eq:cc_vf}$. In fact, every tth state is fully determined once all
nth spline coefficients $\bm{a}^{(n)}_{t}$ are fixed, and the tth observation
$\bm{o}_{t}$ is received. This can be explicitly described by a state update
mechanism through a deterministic mapping, as
$\bm{s}_{t+1}=F(\bm{s}_{t},\bm{a}_{t},\bm{o}_{t+1})$, allowing us to identify
all visitable states handily. On the other hand, selecting all nth spline
coefficients at time step t can be understood as an action. This is because
every nth function section $g^{(n)}_{t}$, as defined in (7), is fully
determined as soon as $\bm{a}^{(n)}_{t}$ is chosen and hence so is the
resulting piece of the subsequent signal reconstruction. Explicitly, we
construct every tth action by concatenating all nth spline coefficients, i.e.,
$\bm{a}_{t}=[{\bm{a}^{(1)}_{t}}^{\top},\dots,{\bm{a}^{(N)}_{t}}^{\top}]^{\top}$,
and denote the action space as $\mathcal{A}\subseteq\mathbb{R}^{2N\rho}$. We
distinctively indicate the set of actions accessible from a given state
$\bm{s}_{t}$ as the admissible action set
$\mathcal{A}(\bm{s}_{t})\subseteq\mathcal{A}$. Particularly, the tth action is
selected through a stationary parametric policy
$\bm{\mu}_{\bm{\theta}}:\mathcal{S}\to\mathcal{A}$. Stationary policies do not
change over time, and they are suitable for making decisions in problems
governed by stationary spatiotemporal dynamics and with varying time steps. On
the other hand, policy parametrization helps to reduce the pool of candidate
policies and allows, with some flexibility, to accommodate the problem
structure, e.g., spatiotemporal dynamics and constraints, into the policy.

In this context, any reconstructed signal can be assessed and therefore
compared by computing a cumulative cost
$K:\mathcal{S}\times\mathcal{A}\to\mathbb{R}$ through the traversed state-
action pairs by following a given policy $\bm{\mu}_{\bm{\theta}}$. From here,
we can tune the policy parameters via policy training by solving the following
optimization problem

$\displaystyle\text{arg}\underset{\bm{\theta}\in\mathbb{R}^{P}}{\text{ min }}$
$\displaystyle\sum_{m=1}^{M}\sum^{T}_{t=1}K\left(\bm{s}_{m,t},\bm{\mu}_{\bm{\theta}}(\bm{s}_{m,t})\right)$
(13a) s. to:
$\displaystyle\bm{s}_{m,t}=F\left(\bm{s}_{m,t-1},\bm{\mu}_{\bm{\theta}}(\bm{s}_{m,t-1}),\bm{o}_{m,t}\right),\forall
m,t,$ (13b)
$\displaystyle\bm{\mu}_{\bm{\theta}}(\bm{s}_{m,t})\in\mathcal{A}(\bm{s}_{m,t}),\forall
m,t,$ (13c)

where the integer $M$ denotes the number of example multivariate time series,
indexed by m, and where all the stream of observations $\bm{o}_{m,t}$ and
initial states $\bm{s}_{m,0}$ are given. Thus, if a sufficiently large
collection of representative multivariate time series is available, then a
trained policy will outperform, on average, any other policy with different
parameter values when tested over a previously-unseen stream of observations.

###  IV-B Policy architectures and variations

TABLE I: Per time-series cost $\kappa(\cdot,\cdot)$ and admissible action set
$\mathcal{A}(\cdot)$ for the consistent
(${\color[rgb]{.9,0.9,0.9}\blacksquare}$) and smoothing ($\square$) policy
variations.

Per time-series cost evaluated at the tth state-action pair
---
${\bm{a}^{(n)}_{t}}^{\top}\bm{M}_{t}\bm{a}^{(n)}_{t}$
$\left(\bm{p}_{t}^{\top}(x_{t})\bm{a}^{(n)}_{t}-y^{(n)}_{t}\right)^{2}+\eta\,{\bm{a}^{(n)}_{t}}^{\top}\bm{M}_{t}\bm{a}^{(n)}_{t}$,
with $\eta>0$

Admissible action set for the tth state
---
$\begin{aligned}
\left\\{\bm{a}\in\mathcal{A}\,:\right.\,&\left|\bm{p}^{\top}_{t}(x_{t})\bm{a}^{(n)}-y^{(n)}_{t}\right|\leq\epsilon^{(n)},\\\
&\left.[\bm{a}^{(n)}]_{1:\rho}=\bm{e}^{(n)}_{t-1},\,\forall
n\in[1,N]\right\\}\end{aligned}$
$\left\\{\bm{a}\in\mathcal{A}\,:\,[\bm{a}^{(n)}]_{1:\rho}=\bm{e}^{(n)}_{t-1},\,\forall
n\in[1,N]\right\\}$

In this work, we are interested in stationary parametric policies of the form

$\bm{\mu}_{\bm{\theta}}(\bm{s}_{t})=\text{arg}\underset{\bm{a}\in\mathcal{A}(\bm{s}_{t})}{\text{
min
}}\left\\{K(\bm{s}_{t},\bm{a})+J_{\bm{\theta}}(\bm{s}_{t},\bm{a};\bm{h}_{t})\right\\},$
(14)

where $K$ denotes the same cost used in the objective (13a), and the map
$J_{\bm{\theta}}:\mathcal{S}\times\mathcal{A}\times\mathbb{R}^{L}\to\mathbb{R}$
is a parametric approximation of the cost-to-go [27]. Here,
$\bm{h}_{t}\in\mathbb{R}^{L}$ represents a latent state, common for all $N$
time series at the tth time step, and it plays the role of policy memory [28,
29, 30]. Notice that the choice of $J_{\bm{\theta}}$ determines the
architecture of the policy, i.e., the relation between the $P$ parameters
within $\bm{\theta}$. In this work, we adopt three possible policy
architectures.

First, we consider a _myopic_ policy consisting of a parameterless policy
architecture where $J_{\bm{\theta}}(\bm{s}_{t},\bm{a};\bm{h}_{t})=0$ for all
$t$. The myopic policy is arguably the simplest architecture-wise but ignores
the cost-effect of current actions in future decisions. Second, we propose an
_RNN-based_ policy which corresponds to the following parametric cost-to-go
approximation

$J_{\bm{\theta}}(\bm{s}_{t},\bm{a};\bm{h}_{t})=\sum^{N}_{n=1}\lambda_{n}\left\|\bm{a}^{(n)}-\begin{bmatrix}\bm{0}_{\rho}\\\
\bm{r}^{(n)}_{t}\end{bmatrix}\right\|^{2}_{2},$ (15)

where the vectors $\bm{r}^{(n)}_{t}\in\mathbb{R}^{\rho}$, for all
$n\in\mathbb{N}^{[1,N]}$, and $\bm{h}_{t}$ represent the output value and
latent state of an RNN
$R_{\bm{\theta}^{\prime}}:\mathcal{S}\times\mathbb{R}^{L}\to\mathbb{R}^{N\rho}\times\mathbb{R}^{L}$,
respectively. They are obtained through the following relation

$R_{\bm{\theta}^{\prime}}(\bm{s}_{t};\bm{h}_{t})=\begin{bmatrix}\bm{r}_{t}\\\
\bm{h}_{t+1}\end{bmatrix},$ (16)

where
$\bm{r}_{t}=[{\bm{r}^{(1)}_{t}}^{\top},\dots,{\bm{r}^{(N)}_{t}}^{\top}]^{\top}$,
$\bm{\theta}^{\prime}\in\mathbb{R}^{P-N}$ and
$\bm{\theta}=[\lambda_{1},\dots,\lambda_{N},{\bm{\theta}^{\prime}}^{\top}]^{\top}\in\mathbb{R}^{N}_{+}\times\mathbb{R}^{P-N}$.
The vector $\bm{0}_{\rho}$ is used in (15) to disregard the components of
every $\bm{a}^{(n)}$ that are already decided by imposing the continuity
constraints in (10). The RNN-based policy promotes actions close, in the sense
of a Euclidean distance, to the output of a certain RNN, in this case,
$R_{\bm{\theta}^{\prime}}$. Thus, an RNN that successfully captures the
spatiotemporal dynamics of the process generating the multivariate time series
can lead to admissible actions with a reduced cost-effect in future decisions.
Finally, we utilize a _batch_ policy involving an exact cost-to-go
approximation. That is, at time step $t$,
$J_{\bm{\theta}}(\bm{s}_{t},\bm{a};\bm{h}_{t})=J^{*}_{t+1}(F(\bm{s}_{t},\bm{a},\bm{o}_{t+1}))$
where the mapping $F$ refers to the state update mechanism described in Sec.
IV-A, and $J^{*}_{t+1}(\bm{s}_{t+1})$ is the optimal cost for a tail
subproblem that starts at state $\bm{s}_{t+1}$, and time step $t+1$, and ends
at time step $T$. It is computed as

$\displaystyle
J^{*}_{t}(\bm{s}_{t})=\underset{\\{\bm{a}_{\tau}\in\mathcal{A}(\bm{s}_{\tau})\\}}{\text{
min }}$
$\displaystyle\sum^{T}_{\tau=t}K\left(\bm{s}_{\tau},\bm{a}_{\tau}\right)$
(17a) subject to:
$\displaystyle\bm{s}_{\tau+1}=F\left(\bm{s}_{\tau},\bm{a}_{\tau},\bm{o}_{\tau+1}\right),\,\forall\tau,$
(17b)

for all $t\in\mathbb{N}^{[1,T]}$, being $K$ the same cost used in (13a). The
name of the policy comes from the fact that you need the whole batch of
streamed observations in advance to construct the exact cost-to-go, i.e.,
solving (17). Thus, in practice, it cannot be used under zero-delay
requirements. The batch policy is parameterless, and by definition, it yields
the signal reconstruction with the lowest possible cost.

From the presented policy architectures, notice that we can build different
policy variations depending on the choice of the cost $K$ and admissible set
$\mathcal{A}(\cdot)$. In this work, we devise a policy yielding consistent
signal reconstructions and present another policy variation designed to
reconstruct noisy data smoothly but not consistently. Both presented policy
variations, called from now on consistent and smoothing, respectively, are
based on the fact that the roughness, as expressed in (3a), can be
equivalently computed from a trajectory of state-action pairs, as we have
shown in our previous work ([17, Proposition 2]). Specifically, the roughness
of a spline-based signal estimate
$\bm{f}_{T}=[f^{(1)}_{T},\dots,f^{(N)}_{T}]^{\top}$ over the domain $(0,X]$
can be computed as

$\sum^{N}_{n=1}\int^{X}_{0}\left(D^{\rho}_{x}f^{(n)}_{T}(x)\right)^{2}\,dx=\sum^{T}_{t=1}\sum^{N}_{n=1}{\bm{a}^{(n)}_{t}}^{\top}\bm{M}_{t}\bm{a}^{(n)}_{t},$
(18)

where the elements of $\bm{M}_{t}\in\bm{S}^{2\rho}_{+}$ are obtained by

$\left[\bm{M}_{t}\right]_{i,j}=\begin{cases}0&\text{if }i\leq\rho\text{ or
}j\leq\rho,\\\
u_{t}^{i+j-2\rho-1}\frac{\prod^{\rho}_{k=1}(i-k)(j-k)}{i+j-2\rho-1}&\text{otherwise},\end{cases}$
(19)

with $u_{t}\triangleq x_{t}-x_{t-1}$, except $u_{1}=x_{1}$. Similarly, we can
construct a roughness-aware penalization from a cost of the form

$K(\bm{s}_{t},\bm{a}_{t})=\sum^{N}_{n=1}\kappa\left(\bm{s}^{(n)}_{t},\bm{a}^{(n)}_{t}\right),$
(20)

where $\bm{s}^{(n)}_{t}$ is an auxiliar per time-series state constructed from
the tth state,
$\bm{s}_{t}=[x_{t-1},\bm{o}^{\top}_{t},\bm{e}^{\top}_{t-1}]^{\top}$, as
$\bm{s}^{(n)}_{t}=[x_{t-1},x_{t},y^{(n)}_{t},\epsilon^{(n)},{\bm{e}^{(n)}_{t-1}}^{\top}]^{\top}$
and where $\kappa$ denotes a per time-series cost. Table I summarizes our
proposed policy variations.

Briefly, the consistent policy variation ensures smoothness while enforcing
that the reconstructed signal passes through the quantization intervals, thus
yielding a consistent signal reconstruction. On the other hand, the smoothing
policy variation extends the smoothing interpolation approach, introduced in
our previous work [17], to multivariate time series. Specifically, it aims to
reduce the squared residuals with the center of the intervals while promoting
smoothness, but it does not enforce that the ongoingly reconstructed signal
passes through the quantization intervals, as illustrated in Fig. 3, hence
neither consistency.

Figure 3: Snapshot of an ongoing signal reconstruction using the smoothing
policy variation.

###  IV-C Policy evaluation

Figure 4: Conceptualization of the hyperplane projection $P_{\hat{H}^{(n)}_{t}}$ in the inner product space $\left(\mathbb{R}^{\rho},\langle\cdot,\cdot\rangle_{\bm{D}^{(n)}_{t}}\right)$. TABLE II: Terms of the quadratic form for the smoothing policy variation. This notation is shared with the rest of the paper with the incorporation of $\bm{P}_{t}\triangleq\bm{p}_{t}(x_{t})\bm{p}_{t}(x_{t})^{\top}$ and $\bm{v}^{(n)}_{t}\triangleq[\bm{0}^{\top}_{\rho},{\bm{r}^{(n)}_{t}}^{\top}]^{\top}$. | $\bm{A}^{(n)}_{t}$ | $\bm{b}^{(n)}_{t}$
---|---|---
Myopic | $\bm{P}_{t}+\eta\bm{M}_{t}$ | $-2y^{(n)}_{t}\bm{p}_{t}(x_{t})$
RNN | $\bm{P}_{t}+\eta\bm{M}_{t}+\lambda_{n}\bm{I}_{2\rho}$ | $-2(y^{(n)}_{t}\bm{p}_{t}(x_{t})+\lambda_{n}\bm{v}^{(n)}_{t})$
TABLE III: Terms of the quadratic form for the consistent policy variation. Same notation as Table II. | $\bm{A}^{(n)}_{t}$ | $\bm{b}^{(n)}_{t}$
---|---|---
Myopic | $\bm{M}_{t}$ | $\bm{0}_{2\rho}$
RNN | $\bm{M}_{t}+\lambda_{n}\bm{I}_{2\rho}$ | $-2\lambda_{n}\bm{v}^{(n)}_{t}$
Figure 5: Representation of a tth forward pass of our RNN-based policy. CAT
refers to a vector concatenation layer. Evaluating the policy consists of
solving a convex optimization problem with respect to $\bm{a}$, see (14) and
Sec. IV-C. Therefore, its last step can be seen as a convex optimization layer
[31].

By design, evaluating the considered policies involves solving optimization
problem (14). The myopic and RNN-based policy architectures, in any of their
variations given in Table I, can be equivalently evaluated by solving the
following quadratic convex problem

$\bm{\mu}_{\bm{\theta}}(\bm{s}_{t})=\text{arg}\underset{\bm{a}\in\mathcal{A}(\bm{s}_{t})}{\text{min}}\left\\{\sum^{N}_{n=1}{\bm{a}^{(n)}}^{\top}\bm{A}^{(n)}_{t}\bm{a}^{(n)}+{\bm{b}^{(n)}_{t}}^{\top}\bm{a}^{(n)}\right\\}$
(21)

where the terms $\bm{A}^{(n)}_{t}\in\bm{S}_{+}^{2\rho}$ and
$\bm{b}^{(n)}_{t}\in\mathbb{R}^{2\rho}$ depend on the chosen policy
architecture and variation, as displayed in Tables II and III. The policy
evaluation in (21) is presented as an intermediate step for clarity. Then,
after some algebraic steps, it can be equivalently reformulated as

$\bm{\mu}_{\bm{\theta}}(\bm{s}_{t})=\text{arg}\underset{\bm{a}\in\mathcal{A}(\bm{s}_{t})}{\text{min}}\left\\{\sum^{N}_{n=1}\|\bm{a}^{(n)}-\bm{\alpha}^{(n)}_{t}\|^{2}_{\bm{D}^{(n)}_{t}}\right\\}$
(22)

where every $\bm{D}^{(n)}_{t}\in\bm{S}^{\rho}_{++}$ and
$\bm{\alpha}^{(n)}_{t}\in\mathbb{R}^{\rho}$ are computed from the elements of
the tth state $\bm{s}_{t}$, through $\bm{A}^{(n)}_{t}$ and $\bm{b}^{(n)}_{t}$,
as

$\bm{D}^{(n)}_{t}=\bm{B}_{2}^{\top}\bm{A}^{(n)}_{t}\bm{B}_{2},$ (23)

and

$\bm{\alpha}^{(n)}_{t}=-{\bm{D}^{(n)}_{t}}^{-1}\left(\bm{B}_{2}\bm{A}^{(n)}_{t}\bm{B}_{1}\bm{e}^{(n)}_{t-1}+\frac{1}{2}\bm{B}_{2}^{\top}\bm{b}^{(n)}_{t}\right),$
(24)

for all $n\in\mathbb{N}^{[1,N]}$ and $t\in\mathbb{N}^{[1,T]}$. The auxiliary
matrices
$\bm{B}_{1}\triangleq[\bm{I}_{\rho},\bm{0}_{\rho\times\rho}]^{\top}\in\mathbb{R}^{2\rho\times\rho}$
and
$\bm{B}_{2}\triangleq[\bm{0}_{\rho\times\rho},\bm{I}_{\rho}]^{\top}\in\mathbb{R}^{2\rho\times\rho}$
are defined for notational brevity. The evaluation form in (22) is
analytically preferable since it corresponds to a minimum norm problem or a
hyperslab projection problem [32, 17], for the smoothing or consistent policy
variations, respectively. This is because the aforementioned problems admit
the following closed-form solution

$\bm{\mu}_{\theta}(\bm{s}_{t})=\begin{bmatrix}{\bm{a}^{(1)}_{t}}^{*}\\\
\vdots\\\ {\bm{a}^{(N)}_{t}}^{*}\end{bmatrix},\text{ with
}\,{\bm{a}^{(n)}_{t}}^{*}=\begin{bmatrix}\bm{e}^{(n)}_{t-1}\\\
\zeta^{(n)}_{t}\left(\bm{\alpha}^{(n)}_{t}\right)\end{bmatrix},\,\forall n,$
(25)

where $\zeta^{(n)}_{t}:\mathbb{R}^{\rho}\to\mathbb{R}^{\rho}$ is the identity
operator for the smoothing policy variation and the projection onto the
hyperslab

$\tilde{\mathcal{H}}^{(n)}_{t}\triangleq\left\\{\bm{\beta}\in\mathbb{R}^{\rho}:\left|\bm{\beta}^{\top}\bm{q}_{t}-w^{(n)}_{t}\right|\leq\epsilon^{(n)}\right\\},$
(26)

for the consistent one, where
$\bm{q}_{t}\triangleq[\bm{p}_{t}(x_{t})]_{\rho+1:2\rho}$ and
$w^{(n)}_{t}\triangleq
y^{(n)}_{t}-{\bm{e}^{(n)}_{t-1}}^{\top}[\bm{p}_{t}(x_{t})]_{1:\rho}$ for all
$t\in\mathbb{N}^{[1,T]}$ and $n\in\mathbb{N}^{[1,N]}$. Every hyperslab
$\tilde{\mathcal{H}}^{(n)}_{t}$ can be visualized as the set of all points
which belong between and onto the hyperplanes

$\displaystyle\hat{H}^{(n)}_{t}\triangleq\left\\{\bm{\beta}\in\mathbb{R}^{\rho}:\bm{\beta}^{\top}\bm{q}_{t}=w^{(n)}+\epsilon^{(n)}\right\\},$
(27a)
$\displaystyle\check{H}^{(n)}_{t}\triangleq\left\\{\bm{\beta}\in\mathbb{R}^{\rho}:\bm{\beta}^{\top}\bm{q}_{t}=w^{(n)}-\epsilon^{(n)}\right\\}.$
(27b)

Based on this observation, the projection onto the hyperslab
$\tilde{\mathcal{H}}^{(n)}_{t}$ can be computed as

$\zeta^{(n)}_{t}(\bm{\beta})=\begin{cases}P_{\hat{H}^{(n)}_{t}}(\bm{\beta}),&\text{if
}\bm{\beta}^{\top}\bm{q}_{t}>w^{(n)}_{t}+\epsilon^{(n)},\\\
\bm{\beta},&\text{if
}\left|\bm{\beta}^{\top}\bm{q}_{t}-w^{(n)}_{t}\right|\leq\epsilon^{(n)},\\\
P_{\check{H}^{(n)}_{t}}(\bm{\beta}),&\text{if
}\bm{\beta}^{\top}\bm{q}_{t}<w^{(n)}_{t}-\epsilon^{(n)},\end{cases}$ (28)

where $P_{\hat{H}^{(n)}_{t}}:\mathbb{R}^{\rho}\to\hat{H}^{(n)}_{t}$ and
$P_{\check{H}^{(n)}_{t}}:\mathbb{R}^{\rho}\to\check{H}^{(n)}_{t}$ denote the
projections onto the hyperplanes $\hat{H}^{(n)}_{t}$ and
$\check{H}^{(n)}_{t}$, respectively. These projections can be computed in
closed form as

$\displaystyle
P_{\hat{H}^{(n)}_{t}}(\bm{\beta})=\bm{\beta}-\left(\bm{\beta}^{\top}\bm{q}_{t}-w^{(n)}_{t}-\epsilon^{(n)}\right)\frac{\left(\bm{D}^{(n)}_{t}\right)^{-1}\bm{q}_{t}}{\bm{q}_{t}^{\top}\left(\bm{D}^{(n)}_{t}\right)^{-1}\bm{q}_{t}},$
(29a) $\displaystyle
P_{\check{H}^{(n)}_{t}}(\bm{\beta})=\bm{\beta}-\left(\bm{\beta}^{\top}\bm{q}_{t}-w^{(n)}_{t}+\epsilon^{(n)}\right)\frac{\left(\bm{D}^{(n)}_{t}\right)^{-1}\bm{q}_{t}}{\bm{q}_{t}^{\top}\left(\bm{D}^{(n)}_{t}\right)^{-1}\bm{q}_{t}}.$
(29b)

See Fig. 4 for some intuition behind the closed-form projection in (29a). The
same reasoning applies to the other hyperplane.

###  IV-D Policy training

Unlike the myopic and batch policy architectures, which are parameterless, the
RNN-based architecture requires to be tuned, in this case, through (13), for
adequate performance.

Evaluating the RNN-based policy can be seen as performing a forward pass from
a deep learning point of view [33]. From this perspective, the policy training
problem (13) can be tackled via backpropagation through time [34]. Here, the
closed-form evaluation described in Sec. IV-C allows us to compute and
propagate the gradient of the tth action $\bm{a}_{t}$ with respect to the
parameters contained in $\bm{\theta}$ without the need of unrolling numerical
optimizers [35] or other specialized tools [31]. An overall scheme of our
chosen RNN-based policy and its components can be found in Fig. 5.

##  V Experimental reconstruction error-rate decay

In this section, we experimentally analyze the reconstruction error-rate decay
incurred by the consistent variation of the RNN-based policy (see Sec. IV-B),
with respect to the oversampling ratio and quantization step size of the
acquisition system.

To this end, we first set a specific function space of multivariate smooth
signals, as described in Sec. III-A. For simplicity, we fix the dimensionality
of the multivariate signal to $N=2$, and specify the roughness by choosing
$\rho=2$. Then, we describe a form to generate multivariate signals from
$\mathcal{V}^{2}_{2}$. These multivariate signals, along with an acquisition
system with settled oversampling ratio and quantization stepsize values, shape
the time series data of the experiment. Separately, we use the myopic policy
architecture as a benchmark, or base method, which does not exploit the
spatiotemporal dependencies among time series, and the batch policy
architecture as a baseline, or best possible performance indicator. One
central goal of this section is to compare the performance of the RNN-based
policy architecture that considers the possible spatial relations between the
time series (multivariate) against a policy based on the same RNN architecture
but ignoring the spatial relations (univariate). This is done for both policy
variations: consistent and smoothing.

###  V-A Data generation

Figure 6: Scheme of the acyclic VAR(1) process generating the $2$ series of
knots for any mth sequence. The thickness of the arrows represents the
magnitude of the autoregressive parameters.

Unlike real-world acquired signals, which usually include uncertainty due to
unidentified external noise sources or lack of knowledge about their
underlying generating process, synthetic signals can be generated in a
controlled environment. On the one hand, synthetic data is usually a
simplified version of real-world data, hence not always accurately
representing the real world. On the other hand, the mechanisms used to
generate synthetic data are fully known. This is especially convenient in
studies involving irreversible transformations such as quantization since
there is no other unidentified source of error than the quantization and
reconstruction methods at work. In this paper, we use synthetically generated
multivariate time series to delimit, as much as possible, the dependencies
between the reconstruction error and the reconstruction method used. Moreover,
the synthetic data we use is complex enough to accommodate the smooth behavior
as well as the spatiotemporal dependencies assumed in Sec. III-A.

First, we generate $M=288$ sequences composed of $N=2$ series of $T_{*}=100$
knots each, where every knot has been computed recurrently from a stable
vector autoregressive VAR(1) process
$\bm{z}_{m,t}=\bm{\Phi}\bm{z}_{m,t-1}+\textbf{w}_{m,t}$ with autoregressive
parameters $\bm{\Phi}\in\mathbb{R}^{2}$, being $[\bm{\Phi}]_{1,2}=0$ to avoid
self feedback, Gaussian innovation
$\textbf{w}_{m,t}\sim\mathcal{N}\left(\bm{0},\text{diag}([0.1,0.1]^{\top})\right)$
and where $\bm{z}_{m,0}=\bm{0}_{2}$ for all $m\in\mathbb{N}^{[1,M]}$. In this
way, the resulting sequences of knots inherit the spatiotemporal relations of
the VAR(1) generator process, as illustrated in Fig. 6. Second, each series of
knots is arranged uniformly, with a unitary period $1/\nu^{*}=1$, and is
interpolated with cubic natural splines, which are the smoothest interpolators
in terms of the $\rho=2$ roughness [36]. Lastly, the resulting continuous
signals are uniformly sampled at a frequency $\nu=R\nu^{*}=R$ with a natural
oversampling ratio $R$ spanning from $1$ to $5$ and quantized by a uniform
midtread quantizer, with a half-step size
$\epsilon\in\\{0.1,0.08,0.06,0.04,0.02\\}$. As a result, we generate a dataset
consisting of $288$ multivariate time series with elements of the form
$\tilde{y}^{(n)}_{m,t}=\text{round}(z^{(n)}_{m,t}/2\epsilon)2\epsilon$ for
every $n\in\\{1,2\\}$, $m\in\mathbb{N}^{[1,288]}$ and
$t\in\mathbb{N}^{[1,T]}$, where $T=100R$.

###  V-B Experimental setup

The dataset is divided into $196$ samples for training, $64$ for validation,
and $32$ for testing. Then, to avoid data leaking, the whole dataset is
standardized along each time series with the mean $\gamma^{(n)}$ and standard
deviation $\sigma^{(n)}$ of the training partition. Explicitly,
$y^{(n)}_{m,t}=(\tilde{y}^{(n)}_{m,t}-\gamma^{(n)})/\sigma^{(n)}$ and
$\epsilon^{(n)}=\epsilon/\sigma^{(n)}$ for all $n\in\\{1,2\\}$,
$m\in\mathbb{N}^{[1,288]}$ and $t\in\mathbb{N}^{[1,100R]}$, resulting in
quantization intervals as described in Sec. III-B. Regarding the RNN-based
policy architecture, it consists of $2$ stacked gated recurrent units (GRU)
layers [37], with a latent state of size $L=48$ and an input of size $32$. The
input and output layers are set as linear layers matching the required
dimensionality(see Fig. 5 for a visual representation), which accordingly
varies depending on whether we reconstruct multivariate time series as a
multivariate signal or as separate individual univariate signals. Lastly, the
RNN-based policies have been trained using the adaptive moments (Adam)
optimizer [38], with $\beta_{1}=0.9$, $\beta_{2}=0.999$, without weight decay,
a learning rate of $0.002$ and with a gradient norm clip value of $0.1$ over
mini-batches of $32$ samples.

###  V-C Results and discussion

Figure 7: Training-validation curves for the consistent variation of the RNN-
based policy. The loss metric is the average total cost per function section
computed from the per time-series cost presented in Table I.

(a) Consistent policy variation.

(b) Smoothing policy variation.

Figure 8: Average MSE decay curves over the test set in log-log scale. The
fitting curves are log-linear functions adjusted from the multivariate average
MSE values. The shaded area represents the confidence in the slope of the fit.
The values in Fig. 8(b) have been computed for an $\eta=0.001$.

The training-validation curves obtained for the consistent variation of the
RNN-based policy, from data acquired with $R=1$ and $\epsilon=0.1$, are
presented in Fig. 7. As expected, the consistent RNN-based policy outperforms
the myopic policy for both the multivariate and univariate reconstruction
cases. Moreover, as expected, the multivariate reconstruction outperforms the
univariate one, illustrating the model capabilities to successfully exploit
spatial dependencies (and not only the temporal ones) from multivariate time
series.

Regarding the reconstruction error-rate decay behavior, we have computed the
average MSE over all time series in the test set. This has been done for every
considered policy architecture (once they are trained, if applicable) and
policy variation. As shown in Fig. 8, we have presented the results by
grouping the consistent and smoothing policy variations versus the
quantization half-step size $\epsilon$ (with fixed $R=1$) or versus the
oversampling ratio $R$ (with fixed $\epsilon=0.1$). The slope of the fitting
curve reflects the average error-rate decay in logarithmic scale incurred by
the multivariate reconstructions of the RNN-based policies. In the case of the
consistent policy variation, the module of the quotient between the slope
obtained by varying the quantization step size ($0.99\pm 0.09$) over the slope
obtained by increasing the oversampling ratio ($-0.812\pm 0.008$) is $1.22\pm
0.11$. This quotient ratio is only (at most within) two standard deviations
from a unit value, suggesting that the balanced error-rate decay property can
be achieved for consistent reconstructions even under zero-delay requirements.
As to the smoothing policy variation (with slope $-0.435\pm 0.002$), we are
interested in its relative performance with respect to the consistent one.
Here, the quotient ratio between the slopes of the fits is $1.87\pm 0.02$,
meaning that the error-rate decay for the consistent policy variation is
nearly doubled with respect to the smoothing one. Again as expected, a
property of consistent reconstructions manifests even under zero-delay
requirements.

##  VI Conclusion

This paper designs a method for zero-delay signal reconstruction from streamed
multivariate time series of quantization intervals, both for the smoothing and
consistent discussed variants. The error-rate decay behavior of the proposed
methods has been empirically analyzed showing that the error incurred by the
consistent variant decreases nearly twice (in logarithmic scale) as fast as
the error incurred by the (non-consistent) smoothing variant, as the
oversampling ratio is increased. Moreover, the former error rate is
practically the same as the rate the error decreases when the quantization
step size is reduced, which is advantageous from an implementation
perspective. Finally, we also observe that successfully exploiting the
spatiotemporal dependencies of the streaming quantized data, in this case via
policy training, helps to reduce the average signal reconstruction error.

## References

  * [1] A. J. Jerri, “The shannon sampling theorem—its various extensions and applications: A tutorial review,” Proceedings of the IEEE, vol. 65, no. 11, pp. 1565–1596, 1977.
  * [2] H. D. Luke, “The origins of the sampling theorem,” IEEE Communications Magazine, vol. 37, no. 4, pp. 106–108, 1999.
  * [3] M. Vetterli, P. Marziliano, and T. Blu, “Sampling signals with finite rate of innovation,” IEEE transactions on Signal Processing, vol. 50, no. 6, pp. 1417–1428, 2002.
  * [4] W. R. Bennett, “Spectra of quantized signals,” The Bell System Technical Journal, vol. 27, no. 3, pp. 446–472, 1948.
  * [5] R. M. Gray, “Quantization noise spectra,” IEEE Transactions on information theory, vol. 36, no. 6, pp. 1220–1244, 1990.
  * [6] N. T. Thao and M. Vetterli, “Deterministic analysis of oversampled a/d conversion and decoding improvement based on consistent estimates,” IEEE Transactions on signal processing, vol. 42, no. 3, pp. 519–531, 1994.
  * [7] N. T. Thao and M. Vetterli, “Reduction of the mse in r-times oversampled a/d conversion o (1/r) to o (1/r/sup 2/),” IEEE Transactions on Signal Processing, vol. 42, no. 1, pp. 200–203, 1994.
  * [8] I. Jovanovic and B. Beferull-Lozano, “Oversampled a/d conversion and error-rate dependence of nonbandlimited signals with finite rate of innovation,” IEEE Transactions on Signal Processing, vol. 54, no. 6, pp. 2140–2154, 2006.
  * [9] T.-T. Nguyen, Deterministic analysis of oversampled A/D conversion and Sigma Delta modulation, and decoding improvements using consistent estimates.  PhD thesis, Columbia University, 1993.
  * [10] B. Beferull-Lozano and A. Ortega, “Efficient quantization for overcomplete expansions in/spl ropf//sup n,” IEEE Transactions on Information Theory, vol. 49, no. 1, pp. 129–150, 2003.
  * [11] S. Rangan and V. K. Goyal, “Recursive consistent estimation with bounded noise,” IEEE Transactions on Information Theory, vol. 47, no. 1, pp. 457–464, 2001.
  * [12] V. K. Goyal, J. Kovačević, and J. A. Kelner, “Quantized frame expansions with erasures,” Applied and Computational Harmonic Analysis, vol. 10, no. 3, pp. 203–233, 2001.
  * [13] J. Kovacevic, V. K. Goyal, and M. Vetterli, “Fourier and wavelet signal processing,” Fourier Wavelets. org, pp. 1–294, 2013.
  * [14] C. Schmidt, C. Kottke, V. Jungnickel, and R. Freund, “High-speed digital-to-analog converter concepts,” in Next-Generation Optical Communication: Components, Sub-Systems, and Systems VI, vol. 10130, pp. 133–141, SPIE, 2017.
  * [15] S. A. Bazaz and B. Tondu, “Minimum time on-line joint trajectory generator based on low order spline method for industrial manipulators,” Robotics and Autonomous Systems, vol. 29, no. 4, pp. 257–268, 1999.
  * [16] T. Marcucci, P. Nobel, R. Tedrake, and S. Boyd, “Fast path planning through large collections of safe boxes,” arXiv preprint arXiv:2305.01072, 2023\.
  * [17] E. Ruiz-Moreno, L. M. López-Ramos, and B. Beferull-Lozano, “A trainable approach to zero-delay smoothing spline interpolation,” 2023.
  * [18] P. Craven and G. Wahba, “Smoothing noisy data with spline functions: estimating the correct degree of smoothing by the method of generalized cross-validation,” Numerische mathematik, vol. 31, no. 4, pp. 377–403, 1978.
  * [19] H. Salehinejad, S. Sankar, J. Barfett, E. Colak, and S. Valaee, “Recent advances in recurrent neural networks,” arXiv preprint arXiv:1801.01078, 2017.
  * [20] O. Kosheleva and V. Kreinovich, “Why physical processes are smooth or almost smooth: A possible physical explanation based on intuitive ideas behind energy conservation,” Mathematical Structures and Modeling, 2021.
  * [21] I. J. Schoenberg, “On interpolation by spline functions and its minimal properties,” in On Approximation Theory/Über Approximationstheorie, pp. 109–129, Springer, 1964.
  * [22] P. J. Green and B. W. Silverman, Nonparametric regression and generalized linear models: a roughness penalty approach.  Crc Press, 1993.
  * [23] G. Wahba, Spline models for observational data.  SIAM, 1990.
  * [24] B. Widrow and I. Kollár, Quantization noise: roundoff error in digital computation, signal processing, control, and communications.  Cambridge University Press, 2008.
  * [25] H. S. Wilf, Algorithms and complexity.  AK Peters/CRC Press, 2002.
  * [26] R. Bellman, “Dynamic programming,” Science, vol. 153, no. 3731, pp. 34–37, 1966.
  * [27] D. Bertsekas, Dynamic programming and optimal control: Volume I, vol. 1.  Athena scientific, 2012.
  * [28] L. Peshkin, N. Meuleau, and L. Kaelbling, “Learning policies with external memory,” arXiv preprint cs/0103003, 2001.
  * [29] A. M. Schäfer, Reinforcement learning with recurrent neural networks.  PhD thesis, Osnabrück, Univ., Diss., 2008, 2008.
  * [30] M. Zhang, Z. McCarthy, C. Finn, S. Levine, and P. Abbeel, “Learning deep neural network policies with continuous memory states,” in 2016 IEEE international conference on robotics and automation (ICRA), pp. 520–527, IEEE, 2016.
  * [31] A. Agrawal, B. Amos, S. Barratt, S. Boyd, S. Diamond, and J. Z. Kolter, “Differentiable convex optimization layers,” Advances in neural information processing systems, vol. 32, 2019.
  * [32] K. Slavakis, P. Bouboulis, and S. Theodoridis, “Online learning in reproducing kernel hilbert spaces,” in Academic Press Library in Signal Processing, vol. 1, pp. 883–987, Elsevier, 2014.
  * [33] I. Goodfellow, Y. Bengio, and A. Courville, Deep learning.  MIT press, 2016.
  * [34] P. J. Werbos, “Backpropagation through time: what it does and how to do it,” Proceedings of the IEEE, vol. 78, no. 10, pp. 1550–1560, 1990.
  * [35] V. Monga, Y. Li, and Y. C. Eldar, “Algorithm unrolling: Interpretable, efficient deep learning for signal and image processing,” IEEE Signal Processing Magazine, vol. 38, no. 2, pp. 18–44, 2021.
  * [36] S. N. Wood, Generalized additive models: an introduction with R.  chapman and hall/CRC, 2006.
  * [37] K. Cho, B. Van Merriënboer, D. Bahdanau, and Y. Bengio, “On the properties of neural machine translation: Encoder-decoder approaches,” arXiv preprint arXiv:1409.1259, 2014.
  * [38] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” arXiv preprint arXiv:1412.6980, 2014.