# SASA: saliency-aware self-adaptive snapshot compressive imaging

###### Abstract

The ability of snapshot compressive imaging (SCI) systems to efficiently
capture high-dimensional (HD) data depends on the advent of novel optical
designs to sample the HD data as two-dimensional (2D) compressed measurements.
Nonetheless, the traditional SCI scheme is fundamentally limited, due to the
complete disregard for high-level information in the sampling process. To
tackle this issue, in this paper, we pave the first mile toward the advanced
design of adaptive coding masks for SCI. Specifically, we propose an efficient
and effective algorithm to generate coding masks with the assistance of
saliency detection, in a low-cost and low-power fashion. Experiments
demonstrate the effectiveness and efficiency of our approach. Code is
available at: https://github.com/IndigoPurple/SASA.

figure While traditional snapshot compressive imaging (SCI) directly generates
random binary matrices as coding masks to obtain compressed measurements, we
propose a framework elegantly with saliency detection to perform self-adaptive
sampling, and thus boosting the quality of compressed measurement and the
performance of reconstruction result.

††∗Corresponding author.††This work is supported in part by the Research
Grants Council (GRF 17201822), by the Research Postgraduate Student Innovation
Award (The University of Hong Kong), and by ACCESS –- AI Chip Center for
Emerging Smart Systems, Hong Kong SAR.

Index Terms—  snapshot compressive imaging, coded image sensing, compressed
sensing, imaging systems, computational imaging

Notation | Description
---|---
${\bf X}\in\mathbb{R}^{H\times W\times C}$ | The video frames we aim to compress and reconstruct.
$\forall c=1,\dots,C$, ${\bf X}_{c}={\bf X}(:,:,c)\in\mathbb{R}^{H\times W}$
is the c-th video frame.
Here, $H$, $W$ and $C$ are the image height, image width and compression rate,
respectively.
${{\bf E}}\in\mathbb{R}^{H\times W}$ | The measurement noise.
${{\bf Y}}\in\mathbb{R}^{H\times W}$ | The compressed measurement.
${\bf A}\in\\{0,1\\}^{H\times W\times C}$ | The sensing matrix we use to compress the video frames.
$\forall h=1,\dots,H$, $\forall w=1,\dots,W$, $\forall c=1,\dots,C$, ${\bf
A}_{c}={\bf A}(:,:,c)\in\\{0,1\\}^{H\times W}$ is the c-th coding mask,
$a_{hwc}={\bf A}(h,w,c)\in\\{0,1\\}$ is the element on the h-th row and w-th
column of the c-th coding mask.
$Bernoulli(p)$ | The Bernoulli distribution, which is the discrete probability distribution of a random variable which
takes the value 1 with probability $p$ and the value 0 with probability
$q=1-p$.
${{\bf S}}\in\mathbb{R}^{H\times W\times D}$ | The saliency maps, where ${{\bf S}}_{d}={{\bf S}}(:,:,d)\in\\{0,1\\}^{H\times W}$, is the d-th saliency map.
Here, value 1 indicates saliency, and value 0 indicates no saliency.
$D$ is the maximum number of detections to examine.
${{\bf P}}\in\mathbb{R}^{H\times W}$ | The sampling probability, where $p_{hw}={{\bf P}}(h,w)\in[0,1]$ is a probability parameter.
Table 1: Key notations and descriptions.

##  1 introduction

Snapshot compressive imaging (SCI) systems address the challenges associated
with high-dimensional (HD) visual signal acquisition and analysis [1]. These
systems capture HD data as compressed two-dimensional (2D) measurements, and
then utilize compressed sensing algorithms to reconstruct the HD data [2, 3,
4, 5]. By employing optical designs, SCI systems enable the acquisition of HD
visual signals in a single snapshot, as opposed to conventional methods. This
capability proves advantages in scenarios where capturing dynamic or fast-
moving scenes is crucial. Therefore, SCI systems offer promising solutions for
real-world applications in fields such as surveillance, biomedical imaging,
remote sensing

The key principle behind SCI lies in the compression of the acquired HD data.
Rather than directly sampling the entire HD signal, SCI systems capture a
compressed version of the data, resulting in compressed measurements. However,
the existing research in this area has predominantly focused on the
reconstruction algorithms, while the study of compression sampling methods
remains lacking, let alone optimizing the mask design. Additionally, the
current sampling methods treat all pixels equally by employing random binary
masks, completely disregarding high-level information such as objects and
saliency. As a result, the importance of different regions of the image is not
properly considered.

An even more serious issue is that the “sampling - reconstruction - tasks”
scheme of traditional SCI is fundamentally limited. On one hand, the complete
ignorance of high-level information in the sampling process leads to wasteful
resource allocation for content-irrelevant computation. On another hand, the
independence of subsequent computer vision tasks hampers the sampling and
design efficiency of SCI [6].

To solve those problems, as Figure SASA: saliency-aware self-adaptive snapshot
compressive imaging shows, we propose an efficient and effective algorithm to
generate adaptive coding masks with the assistance of saliency detection, in a
low-cost and low-power fashion. Our key insight lies in that the sampling
process to acquire compressed measurements should be performed adaptively and
intelligently, which assigns a higher sampling probability to salient image
regions and a lower probability to non-salient regions.

Our main contributions are as follows:

  * •

To the best of our knowledge, this work is the first to design adaptive coding
masks specifically for SCI.

  * •

Our proposed method introduces a low-cost and efficient approach by performing
saliency detection directly on the measurement. This strategy requires only a
few atomic operations, making it highly suitable for real-time applications.
Notably, the processing speed is on average 250fps on a single laptop CPU.

  * •

Quatitative and qualitative experiments on classical datasets demonstrate the
effectiveness of our method.

Fig. 1: The pipeline of traditional SCI and our proposed saliency-aware self-
adaptive SCI, respectively.

##  2 related work

SCI systems consist of two components: an encoding process that compresses the
HD data into a 2D measurement by hardware [7, 8, 1], and a decoding process
that reconstructs the desired signal by algorithms [3, 9, 5, 10, 11, 12, 3,
13, 14, 15, 16, 17, 18, 19, 20, 21, 22]. Nevertheless, the existing research
in this area focused on the reconstruction algorithms, leaving a gap in the
study of compression sampling methods.

Despite recent efforts to bridge the gap between compressive imaging and
computer vision, such as directly performing vision tasks on the compressed
measurements rather than the reconstructed data [23, 6], the fundamental issue
remains unresolved: the complete disregard of high-level information during
the sampling process. As a result, there is a pressing need to conduct
research on optimizing the mask design. Traditional SCI systems [1] treat all
pixels equally, employing random binary masks that fail to consider crucial
high-level information, such as objects and saliency. Consequently, the
importance of different regions within an image is not appropriately taken
into account.

Therefore, in this paper, we aim to pave the first mile towards an adaptive
SCI system that leverages high-level information during the sampling process.

##  3 method

In this section, we first formulate the sampling process in video SCI. Next,
we propose an efficient and effective algorithm to generate adaptive coding
masks, in a low-cost and low-power fashion. Table 1 outlines key notations
used in this paper.

Figure 1 provides an intuitive understanding of our pipeline, jointly with the
comparison to traditional SCI. While traditional SCI directly generates random
binary matrices as coding masks to obtain compressed measurements, our
framework elegantly utilizes saliency detection to perform self-adaptive
sampling.

To begin with, the acquisition of compressed measurements in video SCI could
be mathematically modeled as:

$\displaystyle{{\bf Y}}=\sum^{C}_{c=1}{\bf A}_{c}\odot{\bf X}_{c}+{{\bf E}},$
(1)

where ${{\bf Y}}$ is the compressed measurement, ${\bf A}$ is the sensing
matrix, ${\bf X}$ is the video frames we aim to compress and reconstruct,
${{\bf E}}$ is the measurement noise, and $\odot$ denotes the Hadamard
(element-wise) product.

In traditional SCI [1], the sensing matrix ${\bf A}$ is always generated by
randomly sampling elements from a Bernoulli distribution with probability
$p=0.5$. In our SCI system, we initialize the first sensing matrix by the
traditional method:

$\displaystyle\begin{aligned} {\bf A}^{(0)}=\big{[}a_{hwc}^{(0)}\big{]},\qquad
a_{hwc}^{(0)}\sim Bernoulli(0.5),\\\ \forall h=1,\dots,H,\ \forall
w=1,\dots,W,\ \forall c=1,\dots,C,\end{aligned}$ (2)

which is used to preliminary compress the dynamic scene and capture the first
measurement:

$\displaystyle{{\bf Y}}^{(0)}=\sum^{C}_{c=1}{\bf A}^{(0)}_{c}\odot{\bf
X}^{(0)}_{c}+{{\bf E}}^{(0)}.$ (3)

Once the first measurement is obtained, we perform advanced vision processing
technology, i.e., saliency detection directly on the measurement in stead of
on the reconstructed video, to boost inference speed and reduce bandwidth
occupation, memory footprint and energy consumption. Specifically, we adopt
the light-weight algorithm from [24], denoted as $f$, which takes the
measurement ${{\bf Y}}^{(t)}$ corresponding to the $t$-th video sequence as
input and outputs saliency maps ${{\bf S}}^{(t)}$:

$\displaystyle{{\bf S}}^{(t)}=f({{\bf Y}}^{(t)};D),\qquad t=0,1,2,\dots,$ (4)

where $D$ is a parameter that determines the maximum number of detections to
examine. Based on the estimated saliency maps, we calculate the sampling
probability for the next video sequence by retrieving saliency in the
compressed domain with low bandwidth:

$\displaystyle{{\bf P}}^{(t+1)}=\frac{1}{D}\sum^{D}_{d=1}{{\bf
S}}_{d}^{(t)},\qquad t=0,1,2,\dots,$ (5)

where ${{\bf P}}$ is the sampling probability used to guide sensing matrix
generation for the next video sequence. It assigns higher sampling
probabilities to salient image regions and lower to non-salient, ensuring an
adaptive and effective sampling process. To ensure sampling coverage in
regions without salient events, we then assign a probability of $1/D$ to areas
with zero probability. According to ${{\bf P}}=[p_{hw}]$, the sensing matrix
is updated by:

$\displaystyle\begin{aligned} {\bf
A}^{(t+1)}=\big{[}a^{(t+1)}_{hwc}\big{]},\qquad a^{(t+1)}_{hwc}\sim
Bernoulli\big{(}p_{hw}^{(t+1)}\big{)},\\\ \forall h=1,\dots,H,\ \forall
w=1,\dots,W,\ \forall c=1,\dots,C,\\\ t=0,1,2,\dots.\end{aligned}$ (6)

Here, our key insight lies in that the sampling process of acquiring
compressed measurements should be performed adaptively, which assigns a higher
sampling probability to salient image regions and lower to non-salient
regions.

Finally, the measurement of the next video sequence is captured using the
dynamically updated sensing matrix:

$\displaystyle\begin{aligned} {{\bf Y}}^{(t+1)}=\sum^{C}_{c=1}{\bf
A}^{(t+1)}_{c}\odot{\bf X}^{(t+1)}_{c}+{{\bf E}}^{(t+1)},\\\
t=0,1,2,\dots.\end{aligned}$ (7)

It is worth noting that our method performs saliency detection directly on the
measurement, and requires only a few atomic operations to generate sensing
matrices. In our experiments, the processing speed is on average 250fps on a
single laptop CPU, which means that it enables a low-speed camera to capture
high-speed scenes with $C\times 250$ fps while intelligently and adaptively
varying the coding masks.

Fig. 2: Comparisons of reconstruction results on the dataset Kobe.
Reconstruction is by ADMM-TV [2].  Fig. 3: Comparisons of reconstruction
results on the dataset Vehicle. Reconstruction is by DEQSCI [4].

Dataset | Aerial | Vehicle | Kobe | Traffic | Average
---|---|---|---|---|---
Recon. | ADMM-TV [2]
Trad. SCI | 24.54, 0.877 | 23.88, 0.822 | 25.38, 0.821 | 20.15, 0.740 | 23.49, 0.815
SASA | 24.61, 0.882 | 24.44, 0.862 | 26.35, 0.877 | 20.19, 0.746 | 23.90, 0.842
Recon. | GAP-TV [3]
Trad. SCI | 24.69, 0.861 | 24.51, 0.867 | 25.97, 0.863 | 20.44, 0.763 | 23.90, 0.839
SASA | 24.71, 0.864 | 24.57, 0.872 | 26.80, 0.895 | 20.48, 0.771 | 24.13, 0.851

Table 2: The results in terms of PSNR (dB) and SSIM by different compressed
measurements on classical four datasets.

$D$ | Aerial | Vehicle | Kobe | Traffic | Average
---|---|---|---|---|---
10 | 24.77, 0.871 | 23.69, 0.867 | 25.16, 0.845 | 19.92, 0.740 | 23.39, 0.831
20 | 24.55, 0.859 | 24.13, 0.856 | 25.89, 0.857 | 20.39, 0.759 | 23.74, 0.833
30 | 24.71, 0.864 | 24.57, 0.871 | 26.80, 0.895 | 20.48, 0.770 | 24.13, 0.850
40 | 24.56, 0.865 | 24.53, 0.873 | 26.62, 0.896 | 20.26, 0.754 | 23.99, 0.847
50 | 24.42, 0.862 | 24.41, 0.871 | 26.48, 0.896 | 20.16, 0.745 | 23.87, 0.843

Table 3: Ablation study on different max detection numbers in saliency
detection. Reconstruction is by GAP-TV [3].

##  4 experiment

Experimental Setting.  We initialize the first sensing matrix and compressed
measurement following the traditional SCI [1], and then use the proposed
method to generate adaptive sensing matrices for capturing sequential
compressed measurements. For evaluations, four of the six classical SCI
datasets [9] including Aerial, Vehicle, Kobe and Traffic are adopted, while
two (Drop and Runner) are excluded because the videos are too short and they
contain only a single measurement. Then two classical SCI reconstruction
algorithms ADMM-TV [2] and GAP-TV[3], and a state-of-the-art method DEQSCI [4,
25] are conducted, on different measurements captured by traditional SCI and
our proposed framework (SASA), to reconstruct the videos.

Quantitative and Qualitative Evaluations.  Quantitative comparison results of
different compressed measurements are provided in Table 2, where our method
achieves approximately 0.2$\sim$0.5 dB improvement in PSNR and 0.01$\sim$0.03
in SSIM [26] on average. The improvement indicates our method could
reconstruct images with relatively fine structure, which is confirmed by
qualitative evaluations in Figures 2 and 3. As shown in these figures,
reconstruction results from traditional SCI have more artifacts and
distortions around margins. Our method can maintain a clear and accurate image
structure, thus leading to higher performance. Ours adjusts the coding masks
during the sampling process, seamlessly integrates with existing
reconstruction algorithms, and enhances performance across all baselines.
Furthermore, the reconstruction quality can be further improved using various
video processing techniques [27, 28, 29, 30, 31].

Ablation Study.  In Table 3, we experimented with different values for $D$,
the maximum number of detections to examine in saliency detection, and found
$D=30$ is the best setting.

##  5 conclusion

This work addresses the limitations of existing SCI systems by an efficient
and effective framework to adjust coding masks, paving the first mile towards
adaptive SCI systems.

## References

  * [1] X. Yuan, D. J. Brady, and A. K. Katsaggelos, “Snapshot compressive imaging: Theory, algorithms, and applications,” _IEEE Signal Processing Magazine_ , vol. 38, no. 2, pp. 65–88, 2021.
  * [2] S. H. Chan, X. Wang, and O. A. Elgendy, “Plug-and-play admm for image restoration: Fixed-point convergence and applications,” _IEEE Transactions on Computational Imaging_ , 2016.
  * [3] X. Yuan, “Generalized alternating projection based total variation minimization for compressive sensing,” in _2016 IEEE International Conference on Image Processing (ICIP)_. IEEE, 2016, pp. 2539–2543.
  * [4] Y. Zhao, S. Zheng, and X. Yuan, “Deep equilibrium models for snapshot compressive imaging,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 37, no. 3, 2023, pp. 3642–3650.
  * [5] Q. Yang and Y. Zhao, “Revisit dictionary learning for video compressive sensing under the plug-and-play framework,” in _SPIE Proceedings_ , vol. 12166. SPIE, 2022, pp. 2018–2025.
  * [6] Z. Zhang, B. Zhang, X. Yuan, S. Zheng, X. Su, J. Suo, D. J. Brady, and Q. Dai, “From compressive sampling to compressive tasking: retrieving semantics in compressed domain with low bandwidth,” _PhotoniX_ , vol. 3, no. 1, pp. 1–22, 2022.
  * [7] A. Wagadarikar, R. John, R. Willett, and D. Brady, “Single disperser design for coded aperture snapshot spectral imaging,” _Applied optics_ , vol. 47, no. 10, pp. B44–B51, 2008.
  * [8] P. Llull, X. Liao, X. Yuan, J. Yang, D. Kittle, L. Carin, G. Sapiro, and D. J. Brady, “Coded aperture compressive temporal imaging,” _Optics express_ , vol. 21, no. 9, pp. 10 526–10 545, 2013.
  * [9] X. Yuan, Y. Liu, J. Suo, and Q. Dai, “Plug-and-play algorithms for large-scale snapshot compressive imaging,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2020, pp. 1447–1457.
  * [10] Y. Liu, X. Yuan, J. Suo, D. J. Brady, and Q. Dai, “Rank minimization for snapshot compressive imaging,” _IEEE transactions on pattern analysis and machine intelligence_ , 2018.
  * [11] J. Yang, X. Liao, X. Yuan, P. Llull, D. J. Brady, G. Sapiro, and L. Carin, “Compressive sensing by learning a Gaussian mixture model from measurements,” _IEEE Transaction on Image Processing_ , vol. 24, no. 1, pp. 106–119, January 2015.
  * [12] J. Yang, X. Yuan, X. Liao, P. Llull, G. Sapiro, D. J. Brady, and L. Carin, “Video compressive sensing using Gaussian mixture models,” _IEEE Transaction on Image Processing_ , 2014.
  * [13] M. Qiao, Z. Meng, J. Ma, and X. Yuan, “Deep learning for video compressive sensing,” _APL Photonics_ , vol. 5, no. 3, p. 030801, 2020.
  * [14] S. Zheng, C. Wang, X. Yuan, and H. L. Xin, “Super-compression of large electron microscopy time series by deep compressive sensing learning,” _Patterns_ , vol. 2, no. 7, p. 100292, 2021.
  * [15] L. Wang, M. Cao, Y. Zhong, and X. Yuan, “Spatial-temporal transformer for video snapshot compressive imaging,” _arXiv preprint arXiv:2209.01578_ , 2022.
  * [16] Z. Cheng, B. Chen, R. Lu, Z. Wang, H. Zhang, Z. Meng, and X. Yuan, “Recurrent neural networks for snapshot compressive imaging,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 2022.
  * [17] Z. Meng and X. Yuan, “Perception inspired deep neural networks for spectral snapshot compressive imaging,” in _2021 IEEE International Conference on Image Processing (ICIP)_. IEEE, 2021, pp. 2813–2817.
  * [18] Z. Meng, S. Jalali, and X. Yuan, “Gap-net for snapshot compressive imaging,” _arXiv preprint arXiv:2012.08364_ , 2020.
  * [19] Z. Wu, J. Zhang, and C. Mou, “Dense deep unfolding network with 3d-cnn prior for snapshot compressive imaging,” _arXiv preprint arXiv:2109.06548_ , 2021.
  * [20] Y. Zhao, “Mathematical cookbook for snapshot compressive imaging,” _arXiv preprint arXiv:2202.07437_ , 2022.
  * [21] Y. Zhao, Q. Zeng, and E. Y. Lam, “Adaptive compressed sensing for real-time video compression, transmission, and reconstruction,” in _IEEE International Conference on Data Science and Advanced Analytics (DSAA)_. IEEE, 2023.
  * [22] Z. Y. Chong, Y. Zhao, Z. Wang, and E. Y. Lam, “Solving inverse problems in compressive imaging with score-based generative models,” in _IEEE International Conference on Data Science and Advanced Analytics (DSAA)_. IEEE, 2023.
  * [23] S. Lu, X. Yuan, and W. Shi, “Edge compression: An integrated framework for compressive imaging processing on cavs,” in _2020 IEEE/ACM Symposium on Edge Computing (SEC)_. IEEE, 2020, pp. 125–138.
  * [24] M.-M. Cheng, Z. Zhang, W.-Y. Lin, and P. Torr, “Bing: Binarized normed gradients for objectness estimation at 300fps,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2014, pp. 3286–3293.
  * [25] Y. Zhao, S. Zheng, and X. Yuan, “Deep equilibrium models for video snapshot compressive imaging,” _arXiv preprint arXiv:2201.06931_ , 2022.
  * [26] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: from error visibility to structural similarity,” _IEEE transactions on image processing_ , vol. 13, no. 4, pp. 600–612, 2004.
  * [27] Y. Zhao, M. Ji, R. Huang, B. Wang, and S. Wang, “EFENet: Reference-based video super-resolution with enhanced flow estimation,” in _Artificial Intelligence: First CAAI International Conference, CICAI 2021, Hangzhou, China, June 5–6, 2021, Proceedings, Part I 1_. Springer, 2021, pp. 371–383.
  * [28] Y. Zhao, H. Zheng, Z. Wang, J. Luo, and E. Y. Lam, “MANet: improving video denoising with a multi-alignment network,” in _2022 IEEE International Conference on Image Processing (ICIP)_. IEEE, 2022, pp. 2036–2040.
  * [29] Y. Zhao, G. Li, and E. Y. Lam, “Cross-camera human motion transfer by time series analysis,” _IEEE International Conference on Acoustics, Speech and Signal Processing_ , 2024.
  * [30] Y. Zhao, H. Zheng, M. Ji, and R. Huang, “Cross-camera deep colorization,” in _CAAI International Conference on Artificial Intelligence_. Springer, 2022, pp. 3–17.
  * [31] Y. Zhao, H. Zheng, J. Luo, and E. Y. Lam, “Improving video colorization by test-time tuning,” in _IEEE International Conference on Image Processing (ICIP)_ , 2023.


# Cubature Kalman Filter Based Training of Hybrid Differential Equation
Recurrent Neural Network Physiological Dynamic Models

Ahmet Demirkaya1, Tales Imbiriba1, Kyle Lockwood1, Sumientra Rampersad1,
Elie Alhajjar2, Giovanna Guidoboni3, Zachary Danziger4, Deniz Erdoğmuş1
1Department of Electrical and Computer Engineering, Northeastern University,
Boston, MA 02115, USA; e-mail: {demirkaya.a, t.imbiriba, lockwood.k,
s.rampersad<EMAIL_ADDRESS>2The United States Military
Academy, West Point, NY; e-mail<EMAIL_ADDRESS>3Department of
Electrical Engineering and Computer Science, Department of Mathematics,
University of Missouri, Columbia, MO; e-mail<EMAIL_ADDRESS>4Department of Biomedical Engineering, Florida International University,
Miami, FL; e-mail<EMAIL_ADDRESS>work was supported in part by NIH
(OT2OD030524) and NSF (DMS-1853222/1853303, IIS-1715858, M3X-20040457).

###### Abstract

Modeling biological dynamical systems is challenging due to the
interdependence of different system components, some of which are not fully
understood. To fill existing gaps in our ability to mechanistically model
physiological systems, we propose to combine neural networks with physics-
based models. Specifically, we demonstrate how we can approximate missing
ordinary differential equations (ODEs) coupled with known ODEs using Bayesian
filtering techniques to train the model parameters and simultaneously estimate
dynamic state variables. As a study case we leverage a well-understood model
for blood circulation in the human retina and replace one of its core ODEs
with a neural network approximation, representing the case where we have
incomplete knowledge of the physiological state dynamics. Results demonstrate
that state dynamics corresponding to the missing ODEs can be approximated well
using a neural network trained using a recursive Bayesian filtering approach
in a fashion coupled with the known state dynamic differential equations. This
demonstrates that dynamics and impact of missing state variables can be
captured through joint state estimation and model parameter estimation within
a recursive Bayesian state estimation (RBSE) framework. Results also indicate
that this RBSE approach to training the NN parameters yields better outcomes
(measurement/state estimation accuracy) than training the neural network with
backpropagation through time in the same setting.

###### Index Terms:

Recurrent neural networks, Bayesian filtering, differential equations, retina
circulation.

##  I INTRODUCTION

The dynamic behavior of biological systems is often modeled via sets of
ordinary differential equations (ODEs) that describe the evolution of
biological quantities, such as blood pressure or molecular concentrations,
over time. Such models are based on known physics of the described systems and
data acquired from experiments with humans, animals, cell cultures, or
combinations of such data in multi-scale models. Complex biological systems
consist of many interconnecting processes. Generally, the ODEs for some of
these processes are not known and it may be difficult to acquire data required
to fit model parameters and formulate a mathematical expression underlying the
biophysics. Physics-informed machine learning has emerged as a method to fill
the gaps in models based on differential equations in general. For example,
Raissi et al. developed a deep learning framework to solve problems involving
nonlinear partial differential equations [1], and Weinan et al. used deep
learning to solve high-dimensional parabolic equations [2]. While these
methods are promising, training a physics-informed deep neural network can be
cumbersome due to long cycles and unknown initial states, which reduce the
effectiveness of backpropagation through time. Here, we present the use of
recursive Bayesian estimation and cubature Kalman filters to efficiently train
a physics-informed recurrent neural network.

Recursive Bayesian estimation is a probabilistic method to estimate a
probability density function (PDF) recursively by extracting information about
parameters, or states, of a dynamical system in real time from measurements of
the system output and a system model. The estimate of a state in a given time
instant is based on the old estimate of that state and the new measurement
data. Kalman filters are recursive Bayesian filters for multivariate normal
distributions, which are commonly used in signal processing, navigation and
automatic control. The use of Kalman filters to train recurrent neural
networks (RNNs) has been proposed and implemented successfully in the past [3,
4]. Here we present the first implementation of Kalman-filter-trained RNNs for
the estimation of missing ODEs in models of complex biological systems.

We use a validated model for the blood circulation in the human retina to test
our method [5]. We compare our proposed method to the conventional method,
i.e., backpropagation through time (BPTT), for various levels of additive
noise. Results show that our method was successful in estimating the system
behavior with an NRMSE of 0.038 and was superior to backpropagation for all
tested parameters.

##  II Methods

###  II-A Cubature Kalman Filter

In this paper we train a hybrid ODE and neural network model using recursive
Bayesian state estimation. Specifically, when the predictive and observation
models are linear and Gaussian, the state posterior recursion integrals can be
solved analytically leading to the well-known Kalman time and measurement
update equations [6]. When nonlinear models are in play, required integrals
often become intractable and numerical strategies must be sought. Alternatives
include linearization of nonlinear functions (extended Kalman filters, EKF) or
sampling methods such as particle filters [7], unscented Kalman filters (UKF)
[8], or cubature Kalman filters (CKF) [9], which assume different levels of
system simplicity. For instance, EKF, UKF and CKF assume Gaussianity of the
measurement and transitional models while handling the integration exploiting
this Gaussianity in different ways. While EKF linearizes the models using a
first-order Taylor expansion, UKF and CKF use unscented and cubature rules to
compute integrals such as

$\displaystyle
I(\ell)=\int_{\mathcal{D}}\ell(\boldsymbol{x})p(\boldsymbol{x})d\boldsymbol{x}$
(1)

where $\ell$ is a nonlinear function of $\boldsymbol{x}\in\mathbb{R}^{d_{x}}$
and $p(\boldsymbol{x})=\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$ is a
Gaussian PDF with mean $\boldsymbol{\mu}$ and covariance
$\boldsymbol{\Sigma}$, as a weighted sum of function evaluations of a finite
number of deterministic points. For the third-degree cubature rule, the
integral in (1) can be approximated as

$\displaystyle
I(\ell)\approx\frac{1}{2d_{x}}\sum_{j=1}^{2d_{x}}\ell(\boldsymbol{S}^{\top}{\boldsymbol{\xi}}_{j}+\boldsymbol{\mu})$
(2)

where ${\boldsymbol{\xi}}_{j}=[{\boldsymbol{1}}]_{j}\sqrt{2d_{x}/2}$ are
deterministic points [9], and $\boldsymbol{S}$ is the lower triangular
Cholesky decomposition such that
$\boldsymbol{\Sigma}=\boldsymbol{S}\boldsymbol{S}^{\top}$. It is important to
highlight that the cubature rule demands only two points per dimension of
$\boldsymbol{x}$, i.e., $2d_{x}$, to evaluate the sum in (2), making it more
suitable when working in high-dimensional state-spaces. Assuming Gaussianity
of state posteriors, CKF can solve the integrals required in the Bayesian
recursion as well as the moments (mean and covariance) of the new state
posterior [9].

In contrast, particle filters do not assume any particular distribution,
instead approximating the distribution as a linear combination of Dirac
deltas. Thus, moments of propagated particles can be easily computed. One
drawback of particle filters is the high number of particles needed to
accurately represent distributions. This issue is profoundly aggravated if the
state-space dimension is large, making this filtering strategy unfeasible in
such scenarios [10].

###  II-B Training Recurrent Neural Networks with Augmented Latent Space

Bayesian filtering (BF) has been considered to train neural networks applied
to dynamical systems [11], using EKF [12] and more recently applied to target
tracking in indoor navigation [13], where the state transition and measurement
models can be arbitrary feed-forward neural networks such as:

$\displaystyle\dot{\boldsymbol{s}}$
$\displaystyle=g(\boldsymbol{s},\boldsymbol{u};\tilde{{\boldsymbol{\omega}}})+\tilde{\boldsymbol{n}}$
$\displaystyle\boldsymbol{y}$
$\displaystyle=q(\boldsymbol{s};{\boldsymbol{\phi}})+\boldsymbol{r}$ (3)

where
$\tilde{\boldsymbol{n}}\sim\mathcal{N}(0,\sigma_{\tilde{n}}^{2}{\boldsymbol{I}})$,
$\boldsymbol{s}$ is the vector of latent states, $\boldsymbol{u}$ are the
external inputs, $\boldsymbol{r}$ is the additive measurement noise with
$\boldsymbol{r}\sim\mathcal{N}(0,\sigma_{\boldsymbol{r}}^{2}{\boldsymbol{I}})$,
and $\tilde{{\boldsymbol{\omega}}}$ and ${\boldsymbol{\phi}}$ are the model
parameters of $g$ and $q$, respectively.

The discretized model can be obtained by approximating the derivative via
finite differences leading to

$\displaystyle\boldsymbol{s}_{k}$
$\displaystyle=\boldsymbol{s}_{k-1}+g(\boldsymbol{s}_{k-1},\boldsymbol{u}_{k-1};{\boldsymbol{\omega}})+\boldsymbol{n}_{k-1}$
$\displaystyle\boldsymbol{y}_{k}$
$\displaystyle=q(\boldsymbol{s}_{k};{\boldsymbol{\phi}})+\boldsymbol{r}_{k}.$
(4)

where we assumed that the time-step $\Delta_{t}$ is incorporated in the NN
weights, $\boldsymbol{n}_{k}$ is process noise with
$\boldsymbol{n}_{k}\sim\mathcal{N}(0,\sigma_{n}^{2}{\boldsymbol{I}})$, where
$\sigma_{n}={\Delta_{t}}\sigma_{\tilde{n}}$, and $\boldsymbol{r}_{k}$ is the
additive measurement noise with
$\boldsymbol{r}_{k}\sim\mathcal{N}(0,\sigma_{\boldsymbol{r}}^{2}{\boldsymbol{I}})$.

In order to learn the model parameters using a BF approach, we augment the
latent space as:

$\displaystyle{\boldsymbol{\omega}}_{k}$
$\displaystyle={\boldsymbol{\omega}}_{k-1}$
$\displaystyle\boldsymbol{\phi}_{k}$ $\displaystyle=\boldsymbol{\phi}_{k-1}$
$\displaystyle\boldsymbol{s}_{k}$
$\displaystyle=\boldsymbol{s}_{k-1}+g(\boldsymbol{s}_{k-1},\boldsymbol{u}_{k-1};{\boldsymbol{\omega}}_{k-1})+\boldsymbol{n}_{k-1}$
(5) $\displaystyle\boldsymbol{y}_{k}$
$\displaystyle=q(\boldsymbol{s}_{k};{\boldsymbol{\phi}}_{k})+\boldsymbol{r}_{k}$

which could be rewritten as

$\displaystyle\boldsymbol{x}_{k}$
$\displaystyle=f(\boldsymbol{x}_{k-1},\boldsymbol{u}_{k-1})+{\boldsymbol{\varepsilon}}_{k-1}$
$\displaystyle\boldsymbol{y}_{k}$
$\displaystyle=h(\boldsymbol{x}_{k})+\boldsymbol{r}_{k}$ (6)

where
$\boldsymbol{x}_{k}=[{\boldsymbol{\omega}}_{k}^{\top},\boldsymbol{\phi}_{k}^{\top},\boldsymbol{s}_{k}^{\top}]^{\top}$,
$f(\boldsymbol{x}_{k-1},\boldsymbol{u}_{k-1})=[{\boldsymbol{\omega}}_{k-1},\boldsymbol{\phi}_{k-1},g(\boldsymbol{s}_{k-1},\boldsymbol{u}_{k-1};{\boldsymbol{\omega}}_{k-1})]^{\top}$,
and $h(\boldsymbol{x}_{k})=q(\boldsymbol{s}_{k};\boldsymbol{\phi}_{k})$. The
above model can then be tackled using filtering techniques [12, 13] as
discussed in Section II-A.

###  II-C Application to the Retinal Circulation

The proposed approach was tested on a validated model for the retinal
circulation proposed by Guidoboni et al. [5]. The model consists of four
nonlinear ODEs describing the blood flow through the retinal vasculature and
the central retinal vessels. Leveraging the electric analogy to fluid flow,
electric currents and electric potentials represent blood flow rates and blood
pressures ($P$), while resistors and capacitors represent hydraulic resistance
and vessel compliance [14]. The model includes five vascular compartments,
whose internal pressures are denoted from 1 to 5: central retinal artery
(CRA): $P_{1}$; arterioles: $P_{2}$; capillaries: $P_{3}$; venules: $P_{4}$;
central retinal vein (CRV): $P_{5}$. Pressures $P_{in}$ and $P_{out}$ are
inlet and outlet pressures that drive the network. By applying Kirchhoff’s
laws of currents and voltages and leveraging the constitutive equations of
each resistor and capacitor, the following ODE system is obtained:

$\dot{P_{1}}=\mathcal{F}_{1}(P_{1},P_{2};P_{in})\\\ $ (7)

$\dot{P_{2}}=\mathcal{F}_{2}(P_{1},P_{2},P_{4})$ (8)

$\dot{P_{4}}=\mathcal{F}_{4}(P_{2},P_{4},P_{5})$ (9)

$\dot{P_{5}}=\mathcal{F}_{5}(P_{4},P_{5};P_{out})$ (10)

The explicit expressions of the functional $\mathcal{F}_{i}$ representing the
time rate of change of the pressures $P_{i}$ can be found in [5]. To test our
proposed method, we treat $P_{4}$ as though its dynamics were unknown and
replace $\mathcal{F}_{4}$ in equation (15) with a neural network with a built-
in integrator as follows

$\dot{P_{4}}=q_{4}(P_{2},P_{4},P_{5};{\boldsymbol{\omega}}),$ (11)

where $q_{4}$ is a feed-forward neural network with parameters
${\boldsymbol{\omega}}$, while equations (7), (8) and (10) are used for the
dynamics of pressures $P_{1}$, $P_{2}$ and $P_{5}$. Note that equations
(7)–(10) have recursive relations that become apparent when time is
discretized:
$\boldsymbol{s}_{k+1}=\boldsymbol{s}_{k}+\Delta_{t}\mathbf{f}(\boldsymbol{s}_{k},\boldsymbol{u}_{k})$,
where $\boldsymbol{s}$ consists of four pressure values $P_{i}$, and
$\mathbf{f}$ represents the functions on the right side of equations (7), (8),
(11) and (10). Therefore, the neural network given in equation (11) is trained
using an RNN that encapsulates this discrete-time approximation for
integration over time. In the RNN, state variables are used to represent the
pressures in each retinal vascular compartment. Both backpropagation through
time and cubature Kalman filters can be used as training methodologies. In the
next section, their performances are compared. The feed-forward neural network
$q_{4}$ is encapsulated in the RNN that couples known differential equations
with the NN model for the unknown dynamics. The $q_{4}$ used in the following
numerical analyses of both BPTT and CKF training consists of a single hidden
layer with 20 hidden units and a total of 103 parameters. For the BPTT
training, we used Adam optimizer with a learning rate of $5\times 10^{-3}$.

For the models in (II-B), we set the vector of latent states as discretized
pressure values at time $k$, i.e.,
$\boldsymbol{s}_{k}=[P_{1,k},\ldots,P_{5,k}]^{\top}$,
$g(\boldsymbol{s}_{k};{\boldsymbol{\omega}}_{k})=[\mathcal{F}_{1}(P_{1,k},P_{2,k};P_{in,k}),\;\mathcal{F}_{2}(P_{1,k},P_{2,k},P_{4,k}),\;\,q_{4}(P_{2,k},P_{4,k},\newline
P_{5,k};{\boldsymbol{\omega}}_{k}),\allowbreak\mathcal{F}_{5}(P_{4,k},P_{5,k};P_{out,k})]$,
and $h(\boldsymbol{s}_{k})=\boldsymbol{H}\boldsymbol{s}_{k}$, where
$\boldsymbol{H}$ is a matrix with zeros and ones selecting measured pressures
and blocking unobservable states. For credible practices of modeling and
simulation, we followed "Ten Simple Rules" suggested by Erdemir et. al [15].

##  III RESULTS

We used the retinal circulation model from [5] to demonstrate the
effectiveness of CKF training of a hybrid ODE and RNN model. In order to
compare the performance of CKF and BPTT training, we performed Monte Carlo
simulations with 100 runs each, with random initialization of the model
parameters at each run. As error metrics we consider the mean absolute
percentage error (MAPE) and normalized root-mean-square error (NRMSE), which
are computed as
$\text{MAPE}=100d_{y}^{-1}\allowbreak\sum_{i=1}^{d_{y}}\left(\sum_{k=1}^{n_{t}}|\hat{y}_{i,k}-y_{i,k}|/\sum_{k=1}^{n_{t}}|y_{i,k}|\right)$
and
$\text{NRMSE}=\allowbreak\\{\sum_{i=1}^{d_{y}}[\sum_{k=1}^{n_{t}}(\hat{y}_{i,k}-y_{i,k})^{2}/(n_{t}(y^{\max}_{i}-y^{\min}_{i}))]/d_{y}\\}^{1/2}$,
respectively, with
$\boldsymbol{y}_{k}=[P_{1,k},\ldots,P_{5,k}]^{\top}\in\mathbb{R}^{d_{y}}$, and
$n_{t}$ is the total number of time samples used.

Pressures $P_{in}$ and $P_{out}$ were used as inputs, $P_{1}$, $P_{2}$ and
$P_{5}$ were used as observable states to train the model and $P_{4}$ was
assumed to be not observable and modeled as a latent space. To take potential
noise in the system into account, we introduced white noise to all inputs and
observable states in our experiments. We used five different noise levels with
signal-to-noise ratios (SNR) given in Table I, where resulting quantitative
error metrics are also reported. Both MAPE and NRMSE increased with decreasing
SNR for CKF as well as BPTT. Across different noise levels, the errors for
BPTT were 1.9–5.4 times higher than for CKF.

Figs. 1–3 show the ground truth values, i.e., values obtained using the ODEs
provided by Guidoboni et al. [5], and estimations of the observable and
unobserved states after training the model with noisy pressure values. Fig. 1
presents a single simulation for the CKF-trained model, demonstrating the
noise on the inputs and estimated values. Despite the existence of noise, the
estimated values converge quickly to the ground truth values. Figs. 2 and 3
demonstrate the mean and 95% confidence intervals for the Monte Carlo
simulations using CKF (Fig. 2) and BPTT (Fig. 3) training for simulations with
an SNR of 22.56. For both methods, the estimated values converge to the ground
truth values. However, the CKF-trained model is able to estimate the
unobserved state in fewer steps with higher confidence than the same model
trained with BPTT. Fig. 4 presents learning curves for both methods, showing
that CKF training reaches an error below the noise level faster and converge
in fewer steps than BPTT.

Figure 1: An example simulation for the CKF-trained model with white noise
(SNR of 22.56) added to the inputs and observable states.

Figure 2: Result of CKF training on 100 Monte Carlo simulations with white
noise (SNR of 22.56). Ground truth values (solid lines) and estimated values
(dashed lines) are shown, with 95% confidence intervals for all estimated
states.

Figure 3: Result of BPTT training on 100 Monte Carlo simulations with white
noise (SNR of 22.56). Ground truth values (solid lines) and estimated values
(dashed lines) are shown, with 95% confidence intervals for all estimated
states.

Figure 4: Learning curves of CKF and BPTT training for a simulation with white noise (SNR of 22.56). Results on the training set (solid lines) and test set (dashed lines) are shown in comparison with the noise level (black solid line). TABLE I: MAPE and NRMSE averaged across 100 Monte Carlo experiments SNR |  CKF MAPE |  BPTT MAPE |  CKF NRMSE |  BPTT NRMSE
---|---|---|---|---
49.53 |  1.47 |  5.36 |  0.038 |  0.204
39.52 |  1.82 |  5.85 |  0.044 |  0.206
32.58 |  2.32 |  5.74 |  0.054 |  0.207
29.51 |  2.59 |  5.92 |  0.061 |  0.208
22.56 |  3.54 |  6.61 |  0.093 |  0.208

##  IV Discussion

The experimental results indicate that RNN models can fill modeling gaps in
ODE-based physiological models given that the set of inputs and outputs of the
model can be observed. In the retinal circulation model we obtained error
metrics in terms of mean absolute value reaching thresholds below the noise
floor. When comparing CKF and BPTT results, we observed that CKF provides
significant improvement for both training convergence and test performance as
can be attested by the results presented in Fig. 4 and Table I. This is due to
the fact that BF strategies keep updating the states continuously whenever a
new measurement is available. Although we opted to stop RNN parameter updates
during the testing phase, in practice we could keep updating the parameters
indefinitely, allowing the model to possibly adapt to physiological states not
contained in the training set.

A common problem in solving ODEs is how to determine the initial conditions.
This issue persists in the hybrid ODE-RNN scheme discussed in this paper,
especially for the model trained with BPTT. In theory, this could be solved by
including the initial states as extra model parameters that need to be
estimated in both training and testing phases. However, when testing this
approach we were never able to obtain accurate estimations for the initial
physiological states. On the other hand, by adopting the BF approach we
changed the problem’s perspective from _estimating the previous state_ with
BPTT to _estimating the current state given the past set of measurements_ with
CKF. Thus, in the BF formulation we may disregard the initial states since the
filtering updates quickly move the states near the correct value. This
difference is key in modeling complex systems since hundreds or even thousands
of initial state conditions might have to be found.

The implementation aspects of the BF strategy is also important, since
nonlinear systems often lead to numerical approaches that scale poorly with
the dimension of the state space. Here, we considered third-degree cubature
rules to solve the filtering integrals under the assumption of Gaussianity of
the states. Cubature rules are particularly interesting since the number of
required particles are relatively small (twice the state-space dimension) and
deterministic, and the solution is exact for monomials of degree three or less
[9]. Propagating such particles through the state transition and observation
models, however, leads to multiple independent model evaluations, which
drastically increases the computational cost of the optimization procedure
during the training phase. Nevertheless, independent evaluation of such
particles allows for easy parallelization of the algorithm.

##  V CONCLUSIONS

In this paper we presented a strategy to complement ODE-based models with
RNNs, where we compared two different learning methodologies: backpropagation
and Bayesian filtering with CKF. Using a validated human retinal circulation
model, we showed that both methodologies were capable of approximating the
missing ODE, and that CKF greatly outperformed BPTT due to the time updates of
latent states even during the testing phase. These results show that BF-based
strategies hold great potential for approximating the dynamics of unknown
states in ODE systems.

## References

  * [1] M. Raissi, P. Perdikaris, and G. Karniadakis, “Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations,” _Journal of Computational Physics_ , vol. 378, pp. 686–707, 2019. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0021999118307125
  * [2] W. E, J. Han, and A. Jentzen, “Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations,” _Communications in Mathematics and Statistics_ , vol. 5, no. 4, p. 349–380, Nov 2017. [Online]. Available: http://dx.doi.org/10.1007/s40304-017-0117-6
  * [3] S. Singhal and L. Wu, “Training multilayer perceptrons with the extended Kalman algorithm,” in _Advances in Neural Information Processing Systems_ , 1988, pp. 133–140.
  * [4] G. Puskorius and L. Feldkamp, “Decoupled extended Kalman filter training of feedforward layered networks,” in _Proceedings of the International Joint Conference on Neural Networks_ , 1991, pp. 771–777.
  * [5] G. Guidoboni, A. Harris, S. Cassani, J. Arciero, B. Siesky, A. Amireskandari, L. Tobe, P. Egan, I. Januleviciene, and J. Park, “Intraocular pressure, blood pressure, and retinal blood flow autoregulation: a mathematical model to clarify their relationship and clinical relevance,” _Investigative Ophthalmology & Visual Science_, vol. 55, no. 7, pp. 4105–4118, 2014.
  * [6] S. Särkkä, _Bayesian filtering and smoothing_. Cambridge University Press, 2013, no. 3.
  * [7] M. S. Arulampalam, S. Maskell, N. Gordon, and T. Clapp, “A tutorial on particle filters for online nonlinear/non-gaussian bayesian tracking,” _IEEE Transactions on signal processing_ , vol. 50, no. 2, pp. 174–188, 2002\.
  * [8] E. A. Wan, R. Van Der Merwe, and S. Haykin, “The unscented Kalman filter,” _Kalman filtering and neural networks_ , vol. 5, no. 2007, pp. 221–280, 2001\.
  * [9] I. Arasaratnam and S. Haykin, “Cubature Kalman filters,” _IEEE Transactions on automatic control_ , vol. 54, no. 6, pp. 1254–1269, 2009.
  * [10] T. Imbiriba and P. Closas, “Enhancing particle filtering using gaussian processes,” in _2020 IEEE 23rd International Conference on Information Fusion (FUSION)_. IEEE, 2020, pp. 1–7.
  * [11] S. Haykin, _Kalman filtering and neural networks_. John Wiley & Sons, 2004, vol. 47.
  * [12] D. Erdogmus, J. C. Sanchez, and J. C. Principe, “Modified Kalman filter based method for training state-recurrent multilayer perceptrons,” in _Proceedings of the 12th IEEE Workshop on Neural Networks for Signal Processing_. IEEE, 2002, pp. 219–228.
  * [13] P. Wu, T. Imbiriba, G. LaMountain, J. Vilà-Valls, and P. Closas, “Wifi fingerprinting and tracking using neural networks,” in _Proceedings of the 32nd International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS+ 2019)_ , 2019, pp. 2314–2324.
  * [14] R. Sacco, G. Guidoboni, and A. Mauri, _A comprehensive physically-based approach to modeling in bioengineering and life sciences_. Academic Press. Elsevier, 2019, iSBN: 9780128125182.
  * [15] A. Erdemir, L. Mulugeta, J. P. Ku, A. Drach, M. Horner, T. M. Morrison, G. C. Y. Peng, R. Vadigepalli, W. W. Lytton, and J. G. Myers, “Credible practice of modeling and simulation in healthcare: ten rules from a multidisciplinary perspective,” _Journal of Translational Medicine_ , vol. 18, no. 1, p. 369, Dec. 2020.