# Efficient Sound Field Reconstruction with Conditional Invertible Neural
Networks

Xenofon Karakonstantis, Efren Fernandez-Grande, Peter Gerstoft

###### Abstract

In this study, we introduce a method for estimating sound fields in
reverberant environments using a conditional invertible neural network (CINN).
Sound field reconstruction can be hindered by experimental errors, limited
spatial data, model mismatches, and long inference times, leading to
potentially flawed and prolonged characterizations. Further, the complexity of
managing inherent uncertainties often escalates computational demands or is
neglected in models. Our approach seeks to balance accuracy and computational
efficiency, while incorporating uncertainty estimates to tailor
reconstructions to specific needs. By training a CINN with Monte Carlo
simulations of random wave fields, our method reduces the dependency on
extensive datasets and enables inference from sparse experimental data. The
CINN proves versatile at reconstructing Room Impulse Responses (RIRs), by
acting either as a likelihood model for maximum a posteriori estimation or as
an approximate posterior distribution through amortized Bayesian inference.
Compared to traditional Bayesian methods, the CINN achieves similar accuracy
with greater efficiency and without requiring its adaptation to distinct sound
field conditions.

###### Index Terms:

##  I Introduction

With spatial audio technologies promptly developing, there is an escalating
need for sophisticated sound field reconstruction algorithms. This demand is
driven by the need to replicate, analyze, and control auditory experiences.
This surge of interest spans many fields, from architectural acoustics to
audio signal processing, where researchers are driven to explore a wide range
of computational methods and modeling techniques in pursuit of accurately
represented sound fields.

Sound field reconstruction involves recreating and characterizing sound field
distributions within an enclosed space conditioned on limited pressure
measurements. It finds extensive utility across multiple domains. Sound field
analysis, facilitates the exploration of acoustic environments, providing
insights into propagation characteristics, composite wave interactions, and
sound decay within confined spaces [1, 2]. In virtual reality, this
methodology assists in crafting immersive auditory experiences, enhancing the
realism of virtual environments [3, 4]. Moreover, applications in sound field
control [5, 6], personalized sound zones [7, 8] or reproduction methods [9,
10] benefit from prior knowledge of sound fields, enabling the creation of
optimized auditory spaces to accommodate diverse preferences or needs within
shared environments.

Traditional sound field reconstruction techniques often exploit spatial
regularization methods to enforce smoothness or sparsity constraints on the
reconstructed field. [11, 12, 13, 14, 15] These methods aim to mitigate noise
amplification using prior knowledge about the sound field’s spatial
composition and spectro-temporal properties. Additionally, spatial covariance-
based techniques, such as kernel methods using Bessel functions, capture
spatial interactions and propagate acoustic information across sensor arrays
[16, 17, 18] Further advancements include exploring localized structures
within the data via dictionary learning approaches [19] or spatial
decomposition models that incorporate the physics-based priors over the signal
properties [20, 21]. However, these methods often show limited adaptability,
requiring adjustments whenever the characteristics of a sound field change.
They can also be time-consuming, making drawing inferences less efficient.
Additionally, the analytical tools describing wave propagation usually focus
on the frequency domain [22]. In practical terms, this means that to achieve
tasks such as sound field reproduction, it’s necessary to use a separate model
for each discrete frequency.

The emergence of deep learning, particularly the integration of neural
networks, presents a promising avenue for significant advancement in this
field. Deep learning can capture intricate patterns and statistical
dependencies within sound field data [23, 24, 25], thereby enabling more
accurate and efficient reconstruction processes. Recent years have seen a
notable increase in deep learning within the sound field reconstruction
literature [26, 27, 28, 29]. A pivotal aspect of deep learning for sound field
reconstruction involves integrating principles of physics into the learning
framework. The addition of physical constraints act as a form of inductive
bias, guiding the learning process toward solutions that align with the
underlying physical phenomena governing sound propagation. This incorporation
can manifest explicitly by directly solving the equation describing wave
propagation [30, 31], or implicitly by constraining the deep learning model
solutions within the range of feasible solutions [32, 28].

Within the landscape of deep learning for parameter estimation, invertible
neural networks (INNs) stand out based on the principles and applications of
normalizing flows [33]. When conditioned on an auxiliary variable in the form
of complementary data or measurements—thereby becoming CINNs—these networks
have proven effective in inferring parameters across high-dimensional spaces
[34, 35, 36, 37]. They can approximate complex probability distributions by
employing several bijective transformations, an architecture that presents
multiple advantages: It enables Bayesian inference, allows for efficient
sampling, and facilitates the estimation of uncertainties. This combination of
flexibility and computational efficiency enables real-time inference, making
them very suitable for swiftly processing large-scale spatio-temporal data.

The principle that allows for this efficiency is the concept of “amortization”
over the computational expenses incurred during the initial phase of training
the CINN. Amortized Bayesian inference efficiently delineates the association
between observations and their posterior distributions using a conditional
density estimator, notably a trained CINN [38]. This approach amortizes the
initial training cost, facilitating expedited inference for new observations
using the pre-trained model. Consequently, the method capitalizes on the
initial computational cost, enabling efficient posterior estimation for
subsequent observations without retraining.

This study proposes a CINN architecture with a series of normalizing flows to
address sound field reconstruction. We propose assimilating the forward
process of plane wave propagation from random wave fields satisfying the
Helmholtz equation so that the network is poised for broad applicability
across acoustic environments. Unlike other methods, this methodological choice
circumvents the complexities associated with explicit regularization and
selection of specific models, as the network’s training spans an extensive
range of sound fields and frequencies. Moreover, the versatility of the
proposed approach enables both accurate reconstruction of sound fields and a
framework for rapid prediction of a sound field distribution, without
compromising precision.

The bijective nature of CINNs also enables the model to not only expedite the
sound field parameters but also to assess uncertainty, thereby providing a
measure of confidence in its predictions. This is achieved by drawing samples
from a known analytic base distribution and processing these through the CINN,
thus deriving the parameter distribution post-inference. Consequently, our
approach affords a robust measure of the predictive confidence, integral to
sound field analysis or any desired post-processing application.

To validate the approach, we conduct experiments utilizing data acquired from
a fully furnished auditorium and compare the results with a principled method
for Bayesian inference, specifically a hierarchical plane wave model employing
Markov Chain Monte Carlo (MCMC). Furthermore, we interpret the uncertainty
provided in the invertible networks’ predictions and evaluate its capabilities
at RIR estimation.

The subsequent sections are structured as follows. Section II provides an
overview of the methods employed, which includes the acoustic model,
normalizing flows, details about the architecture and training of the CINN,
and details on a hierarchical plane wave model baseline. Section III outlines
the experimental dataset utilized for model validation, while Section IV
presents details specific to the training and inference procedures, along with
the results and analysis. This includes discussions on a small-scale example,
baseline comparisons, uncertainty quantification, and RIR reconstruction.
Lastly, Section V offers concluding remarks.

##  II Methods

###  II-A Plane wave acoustic model

In an acoustic domain devoid of sources, the linear sound field is governed by
the homogeneous Helmholtz equation, exhibiting time-harmonic dependence. This
implies that at any position $\mathbf{r}_{m}=[x_{m},y_{m},z_{m}]^{T}$ the
measured sound pressure within a spatially extended region
$\Omega\in\mathbb{R}^{3}$ can be expressed as[22]

$\mathbf{p}=\mathbf{H}\mathbf{x}+\mathbf{n},$ (1)

where $\mathbf{p}=\left[p_{1},\ldots,p_{M}\right]^{T}$ denotes the measured
sound pressure vector. The dictionary $\mathbf{H}\in\mathbb{C}^{M\times N}$
comprises $N$ plane waves, with each element
$\mathbf{H}_{mn}=\mathrm{e}^{j\mathbf{k}_{n}\cdot\mathbf{r}_{m}}$ reflecting
the influence of the $n$-th plane wave on the $m$-th measurement position. The
vector $\mathbf{x}=\left[x_{1},\ldots,x_{N}\right]^{T}$ contains the unknown
plane wave coefficients, and $\mathbf{n}=\left[n_{1},\ldots,n_{M}\right]^{T}$
represents the measurement noise at each observation point. The plane wave
coefficients, associated with the direction $\mathbf{k}_{n}$, are crucial for
the acoustic field’s evaluation. The additive noise term $\mathbf{n}$
encompasses various sources of unwanted disturbances or measurement errors
that may affect the accuracy of the measured sound pressure. Only propagating
waves are included in this study, implying far-field assumptions.

The task of determining these plane wave coefficients is critical for
interpolating or extrapolating the sound pressure at positions
$\mathbf{r}_{R}=\left[x_{R},y_{R},z_{R}\right]^{T}$, where measurements are
not available. This process is complicated by the inherent challenges of ill-
posedness, rank deficiency, and underdetermination, primarily due to the
discrepancy between the number of measured pressures $M$ and the number of
plane waves $N$ used to reconstruct the sound field (i.e. $N>M$).

###  II-B Normalizing flows

Normalizing flows are parameterized bijective transformations, designed to
establish an inverse mapping between two distributions. Specifically, a
function $\boldsymbol{g}$ and its exact inverse $\boldsymbol{g}^{-1}$,
facilitate the conversion of $N$-dimensional latent variables
$\mathbf{z}\sim\pi_{z}(\mathbf{z})$ into $N$-dimensional representations of
target data $\mathbf{w}\sim\pi^{*}_{w}(\mathbf{w})$, and vice-versa such that

$\displaystyle\boldsymbol{g}:$
$\displaystyle\,\mathbf{z}\rightarrow\mathbf{w},$ (2)
$\displaystyle\boldsymbol{g}^{-1}:$
$\displaystyle\,\mathbf{w}\rightarrow\mathbf{z}.$ (3)

Given parameters $\boldsymbol{\theta}$ of the transformations $\boldsymbol{g}$
and $\boldsymbol{g}^{-1}$, the ultimate aim of a normalizing flow is to
closely approximate a data distribution $\pi^{*}_{w}(\mathbf{w})$ through a
probability density $q_{\theta}(\mathbf{w})$ derived from the change-of-
variable formula

$\displaystyle q_{\theta}(\mathbf{w})$
$\displaystyle=\pi_{z}\left(\boldsymbol{g}^{-1}(\mathbf{w};\boldsymbol{\theta})\right)\left|\left(\frac{d\boldsymbol{g}^{-1}(\mathbf{w};\boldsymbol{\theta})}{d\mathbf{w}}\right)\right|,$
(4)

where $\left|\cdot\right|$ denotes the determinant of the Jacobian matrix of
$\boldsymbol{g}^{-1}$, crucial for maintaining the computational efficiency of
the parameterized transformation. An invertible neural network $\mathbf{g}$
combines several normalizing flows that are parameterized by neural networks.
This is done to achieve the mapping of a simple tractable distribution
$\mathbf{z}_{0}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$ to the data space
$\mathbf{w}$ using $N_{f}$ such transformations. This compositional nature,
structured as

$\mathbf{w}=\mathbf{g}(\mathbf{z}_{0})=\boldsymbol{g}_{N_{f}}\circ\cdots\circ\boldsymbol{g}_{0}(\mathbf{z}_{0}),$
(5)

ensures that the final variable $\mathbf{w}$ is derived through an intricate
yet computationally tractable process, effectively bridging simple analytical
forms and the multifaceted nature of the data.

###  II-C CINNs for modeling plane wave propagation

The Bayesian framework offers a structured approach for estimating the
coefficients $\mathbf{x}$ condition on the measured pressure $\mathbf{p}$.
These coefficients $\mathbf{x}$ are estimated via the posterior probability
distribution $\pi(\mathbf{x}|\mathbf{p})$ derived from Bayes’ theorem [39]

$\pi(\mathbf{x}|\mathbf{p})=\frac{\pi(\mathbf{p}|\mathbf{x})\pi(\mathbf{x})}{\pi(\mathbf{p})},$
(6)

where $\pi(\mathbf{p}|\mathbf{x})$ represents the likelihood,
$\pi(\mathbf{x})$ is the prior and $\pi(\mathbf{p})$ refers to the evidence.

This research aims to develop a CINN that can closely approximate the
posterior distribution $\pi(\mathbf{x}\mid\mathbf{p})$ over plane wave
coefficients $\mathbf{x}$. This involves an invertible function
$\mathbf{g}:\mathbb{C}^{N}\rightarrow\mathbb{C}^{N}$, parameterized by
$\boldsymbol{\theta}$ and conditioned on the pressure vector $\mathbf{p}$.
Following [34] and [37] we define the approximation of the true posterior by

$q_{\theta}(\mathbf{x}\mid\mathbf{p})\approx\pi(\mathbf{x}\mid\mathbf{p})$ (7)

where a model aims to match the distribution of $\mathbf{x}$ and $\mathbf{p}$
across all possible values. The model reparameterizes $q_{\theta}$ using
network $\mathbf{g}$, to transform between coefficients $\mathbf{x}$ and a
complex Gaussian latent variable $\mathbf{z}_{0}$ as follows

$\mathbf{x}\sim
q_{\theta}(\mathbf{x}\mid\mathbf{p})\Longleftrightarrow\mathbf{x}=\mathbf{g}(\mathbf{z}_{0};\mathbf{p},\boldsymbol{\theta})\text{
with
}\mathbf{z}_{0}\sim\mathcal{N}_{c}(\mathbf{z}_{0}\mid\mathbf{0},\mathbf{I}).$
(8)

To ensure that the outputs of
$\mathbf{g}(\mathbf{z}_{0};\mathbf{p},\boldsymbol{\theta})$ follow the target
posterior $\pi(\mathbf{x}\mid\mathbf{p})$, we minimize the Kullback-Leibler
(KL) divergence between the true posterior $\pi(\mathbf{x}\mid\mathbf{p})$ and
$q_{\theta}$ for any possible observation of pressure $\mathbf{p}$, targeting
neural network parameters $\hat{\boldsymbol{\theta}}$ that fulfill the
objective

$\displaystyle\widehat{\boldsymbol{\theta}}$
$\displaystyle=\underset{\boldsymbol{\theta}}{\operatorname{argmin}}\,\mathbb{E}_{\pi(\mathbf{p})}\left[\mathbb{K}\mathbb{L}\left(\pi(\mathbf{x}\mid\mathbf{p})\|q_{\boldsymbol{\theta}}(\mathbf{x}\mid\mathbf{p})\right)\right]$
(9)
$\displaystyle=\underset{\boldsymbol{\theta}}{\operatorname{argmin}}\,\mathbb{E}_{\pi(\mathbf{p})}\left[\mathbb{E}_{\pi(\mathbf{x}\mid\mathbf{p})}\left[\log\pi(\mathbf{x}\mid\mathbf{p})-\log
q_{\boldsymbol{\theta}}(\mathbf{x}\mid\mathbf{p})\right]\right]$ (10)

Since the $\log$ posterior density $\log\pi(\mathbf{x}\mid\mathbf{p})$ does
not depend on $\boldsymbol{\theta}$, (10) can be further simplified to

$\displaystyle\widehat{\boldsymbol{\theta}}$
$\displaystyle=\underset{\boldsymbol{\theta}}{\operatorname{argmax}}\,\mathbb{E}_{\pi(\mathbf{p})}\left[\mathbb{E}_{\pi(\mathbf{x}\mid\mathbf{p})}\left[\log
q_{\boldsymbol{\theta}}(\mathbf{x}\mid\mathbf{p})\right]\right]$ (11)
$\displaystyle=\underset{\boldsymbol{\theta}}{\operatorname{argmax}}\,\iint\pi(\mathbf{p},\mathbf{x})\log
q_{\boldsymbol{\theta}}(\mathbf{x}\mid\mathbf{p})d\mathbf{p}d\mathbf{x}.$ (12)

This optimization seeks parameters that maximize the log-likelihood of
observing $\mathbf{x}$ given $\mathbf{p}$ for all combinations of $\mathbf{x}$
and $\mathbf{p}$. Since
$\mathbf{g}^{-1}(\mathbf{x};\mathbf{p},\boldsymbol{\theta})=\mathbf{z}$ by
design, using the change of variables formula leads to

$q_{\theta}(\mathbf{x}\mid\mathbf{p})=\pi\left(\mathbf{z}=\mathbf{g}^{-1}(\mathbf{x};\mathbf{p},\boldsymbol{\theta})\right)\left|\left(\frac{d\mathbf{g}^{-1}(\mathbf{x};\mathbf{p},\boldsymbol{\theta})}{d\mathbf{x}}\right)\right|.$
(13)

Thus, we can rewrite our objective as

$\displaystyle\widehat{\boldsymbol{\theta}}$
$\displaystyle=\underset{\boldsymbol{\theta}}{\operatorname{argmax}}\iint\pi(\mathbf{p},\mathbf{x})\log
q_{\boldsymbol{\theta}}(\mathbf{x}\mid\mathbf{p})d\mathbf{p}d\mathbf{x}$ (14)
$\displaystyle=\underset{\boldsymbol{\theta}}{\operatorname{argmax}}\iint\pi(\mathbf{p},\mathbf{x})\left(\log\pi\left(\mathbf{g}^{-1}(\mathbf{x};\mathbf{p},\boldsymbol{\theta})\right)\right.$
$\displaystyle\qquad\qquad\qquad\left.+\log\left|\boldsymbol{J}_{\mathbf{g}^{-1}}\right|\right)d\mathbf{p}d\mathbf{x},$
(15)

where we have abbreviated
$d\mathbf{g}^{-1}(\mathbf{x};\mathbf{p},\boldsymbol{\theta})/d\mathbf{x}$ (the
Jacobian of $\mathbf{g}^{-1}$ evaluated at $\mathbf{x}$ and $\mathbf{p}$ ) as
$\boldsymbol{J}_{\mathbf{g}^{-1}}$. Due to the architecture of the CINN, the
$\log\left|\boldsymbol{J}_{\mathbf{g}^{-1}}\right|$ is easy to compute (see
Sect. II-D).

Given a joint set of $D$ sound field samples
$\left\\{\mathbf{x}_{d},\mathbf{p}_{d}\right\\}_{d=1}^{D}$ we can approximate
the expectations by minimizing the Monte Carlo estimate of the negative of
(15), as delineated by

$\displaystyle\widehat{\boldsymbol{\theta}}$
$\displaystyle=\underset{\boldsymbol{\theta}}{\operatorname{argmin}}\frac{1}{D}\sum_{d=1}^{D}-\log
q_{\theta}\left(\mathbf{x}_{d}\mid\mathbf{p}_{d}\right)$ (16)
$\displaystyle=\underset{\boldsymbol{\theta}}{\operatorname{argmin}}\frac{1}{D}\sum_{d=1}^{D}\left(-\log\pi\left(\mathbf{g}^{-1}(\mathbf{x}_{d};\mathbf{p}_{d},\boldsymbol{\theta})\right)-\log\left|\boldsymbol{J}_{\mathbf{g}^{-1}}^{(d)}\right|\right)$
(17)
$\displaystyle=\underset{\boldsymbol{\theta}}{\operatorname{argmin}}\frac{1}{D}\sum_{d=1}^{D}\left(\frac{\left\|\mathbf{g}^{-1}\left(\mathbf{x}_{d};\mathbf{p}_{d},\boldsymbol{\theta}\right)\right\|_{2}^{2}}{2}-\log\left|\mathbf{J}_{\mathbf{g}^{-1}}^{(d)}\right|\right).$
(18)

By minimizing (18), the CINN acts as a a conditional density estimator and can
efficiently parameterize any sound field through amortized Bayesian inference
(see Sec. II-F). The graphical representation of the CINN is outlined in Fig.
1, which displays both forward and inverse directions of the CINN.

Figure 1: Probabilistic graph of a CINN composed of multiple normalizing flows
$\mathbf{g}_{0}\dots\mathbf{g}_{N_{f}-1}$ and their inverses
$\mathbf{g}_{1}^{-1}\dots\mathbf{g}_{N_{f}}^{-1}$ where the variables
$\mathbf{z}_{0},\mathbf{z}_{1},\dots,\mathbf{z}_{N_{f}-1}$ are latent
variables, $\mathbf{x}$ are the plane wave coefficients, and $\mathbf{p}$ is a
conditional variable, used as an auxiliary input of sound pressure to
normalizing flows

To facilitate interactions among input variables in normalizing flows, each
flow is randomly shuffled through permutations, an invertible linear
transformation with an absolute determinant equal to 1.[40] Specifically,
these transformations are implemented through an LU decomposition for each
flow

$\mathbf{z}_{n+1}=\mathbf{PL}\mathbf{U}\boldsymbol{g}_{n}(\mathbf{z}_{n};\mathbf{p},\boldsymbol{\theta})$
(19)

where $\mathbf{L}$ is lower triangular with ones on the diagonal, $\mathbf{U}$
is upper triangular with non-zero diagonal entries, and $\mathbf{P}$ is a
permutation matrix. The determinant of the transformation is readily obtained
by multiplying the diagonal elements of the upper matrix $\mathbf{U}$ of the
decomposition. Moreover, the inverse transformation can be executed by solving
two triangular systems, one for $\mathbf{U}$ and one for $\mathbf{L}$. More
details are found in [40].

We adopt a strategy where the model is fitted to a Monte-Carlo simulation of
random wave fields [41, 28]. Simulating the forward model of wave propagation
allows the network to learn the summary statistics of the underlying process
[38]. In this particular instance and similar to (1), the sound fields are a
linear combination of $N$ plane waves with randomly sampled coefficients, such
that the pressure at fixed positions is given by

$\mathbf{p}_{d}=\mathbf{H}\mathbf{x}_{d}+\mathbf{n}_{d},$ (20)

where $\mathbf{x}_{i}$ and $\mathbf{n}_{i}$ are random variables given by

$\pi_{x}(\mathbf{x}_{d}|\sigma_{x})=\mathcal{N}_{c}(\mathbf{x}_{i}|0,\sigma_{x}^{2}\mathbf{I}),$
(21)

and

$\pi_{n}(\mathbf{n}_{d}|\sigma_{p})=\mathcal{N}_{c}\left(\mathbf{n}_{d}|0.\sigma_{p}^{2}\mathbf{I}\right),$
(22)

The hyperprior $\sigma_{p}$ is determined by

$\sigma_{p}^{2}=\frac{\mathbb{E}\left[\check{\mathbf{p}}_{d}^{2}\right]}{\operatorname{SNR}},$
(23)

where $\check{\mathbf{p}}_{i}=\mathbf{H}\mathbf{x}_{d}$ and the variable SNR
is uniformly distributed with lower bound $\text{SNR}_{l}$ and upper bound
$\text{SNR}_{h}$, both in decibels, given by

$\pi(\operatorname{SNR})=\mathcal{U}\left(10^{\left(\text{SNR}_{l}/10\right)},10^{\left(\text{SNR}_{h}/10\right)}\right).$
(24)

###  II-D Conditional coupling transform

In practice, CINNs can use a variety of bijective normalizing flows—such as
linear, element-wise, planar, radial, and autoregressive transformations—to
model complex distributions [42]. Among these, the coupling transformation
stands out for its straightforward yet versatile approach to handling high-
dimensional data [43]. This method begins by dividing the input vector
$\mathbf{z}\in\mathbb{C}^{N}$ into two segments
$\mathbf{z}^{A}\in\mathbb{C}^{\mathrm{n}}$ and
$\mathbf{z}^{B}\in\mathbb{C}^{N-\mathrm{n}}$. The segment $\mathbf{z}^{B}$ is
then merged with a corresponding pressure vector
$\mathbf{p}\in\mathbb{C}^{M}$, forming a new vector
$\mathbf{u}^{B}\in\mathbb{C}^{N+M-\mathrm{n}}$. Through a neural network
$\mathrm{NN}$, parameterized by $\boldsymbol{\psi}$, a set of parameters
$\boldsymbol{\phi}$ is derived for the invertible function
$\boldsymbol{h}:\mathbb{C}^{n}\rightarrow\mathbb{C}^{n}$. The transformation,
defined as
$\mathbf{y}=\boldsymbol{g}\left(\mathbf{z};\mathbf{p},\boldsymbol{\theta}\right)$,
is given by

$\displaystyle\boldsymbol{\phi}_{1},\dots,\boldsymbol{\phi}_{n}$
$\displaystyle=\mathrm{NN}\left(\mathbf{u}^{B};\boldsymbol{\psi}\right),$ (25)
$\displaystyle\mathbf{y}$
$\displaystyle=\left[\left[h_{1}\left(z_{1}^{A};\boldsymbol{\phi}_{1}\right),\ldots,h_{\mathrm{n}}\left(z_{\mathrm{n}}^{A};\boldsymbol{\phi}_{\mathrm{n}}\right)\right]^{T},\mathbf{z}^{B}\right]^{T},$
(26)

where each
$h_{i}\left(\cdot;\boldsymbol{\phi}_{i}\right):\mathbb{C}\rightarrow\mathbb{C}$
represents a scalar bijection.

The resulting Jacobian matrix of this transform is block triangular,

$\displaystyle\mathbf{J}_{\boldsymbol{g}}$
$\displaystyle=\left[\begin{array}[]{ll}\partial\boldsymbol{y}^{A}/\partial\boldsymbol{z}^{A}&\partial\boldsymbol{y}^{A}/\partial\boldsymbol{z}^{B}\\\
\partial\boldsymbol{y}^{B}/\partial\boldsymbol{z}^{A}&\partial\boldsymbol{y}^{B}/\partial\boldsymbol{z}^{B}\end{array}\right]$
(29)
$\displaystyle=\left[\begin{array}[]{cc}\mathbf{J}_{\boldsymbol{h}}&\partial\boldsymbol{y}^{A}/\partial\boldsymbol{z}^{B}\\\
\mathbf{0}&\mathbf{I}\end{array}\right],$ (32)

where $\mathbf{J}_{\boldsymbol{h}}$ denotes the Jacobian of the transformed
$\mathbf{z}^{A}$ components, and $\mathbf{I}$ is the identity matrix
corresponding to the unchanged $\mathbf{z}^{B}$ components, while
$\boldsymbol{y}^{A}$ and $\boldsymbol{y}^{B}$ are the respective partioned
outputs of the transform. Achieving exact invertibility in a single pass
involves reversing the described process and utilizing $\boldsymbol{g}^{-1}$,
$\boldsymbol{h}^{-1}$ and $\mathbf{y}$ in place of $\boldsymbol{g}$,
$\boldsymbol{h}$, and $\mathbf{z}$ [44].

###  II-E Rational quadratic splines

In principle, the normalizing flows $\boldsymbol{g}$ have the flexibility to
assume any parametric and invertible form. Among these, rational quadratic
splines (RQS) are particularly noteworthy. They offer a balance of gradient
stability throughout the learning process and versatility in modeling [45].
Rational quadratic splines are functions composed of segments (bins)
representing the ratio of two quadratic polynomials. The segments are
delimited by $K+1$ coordinates known as knots
$\left\\{\left(z^{(k)},y^{(k)}\right)\right\\}_{k=0}^{K}$. For instance, in
the $i$-th dimension of the input $\mathbf{z}^{A}$ which is binned in the
$k$-th segment of the spline, the element-wise functions of (26) take the form

$h_{i}\left(z_{i};\phi_{i}\right)=\frac{\alpha^{(k)}(\xi)}{\beta^{(k)}(\xi)},$
(33)

where $\alpha^{(k)}$ and $\beta^{(k)}$ are two polynomials, whose parameters
are determined by a neural network and
$\xi(z_{i})=\left(z_{i}-z_{i}^{k}\right)/\left(z_{i}^{k+1}-z_{i}^{k}\right)$
represents the normalized distance along the $z$ axis within the
$k^{\text{th}}$ bin. Therefore, $\xi(z_{i})\in[0,1]$. The spline is
constructed within the specified bounding box $[-B,B]\times[-B,B]$, with the
identity function applied elsewhere, while $B$ is typically determined before
training and set to a value that covers the range of the input data
$\mathbf{x}$. More details on the definition of RQS can be found in Appendix
A.

Following [45], the implementation of the monotonic rational quadratic
coupling transform can be outlined by the process below.

1: A neural network takes the combined vector of pressure and coefficients
$\mathbf{u}^{B}$ as input (see (25) Sec. II-D) and generates an unconstrained
parameter vector $\boldsymbol{\phi}_{i}$ for each $i=1,\ldots,n$. Each
$\boldsymbol{\phi}_{i}$ has a length of $3K-1$.

2: Divide $\boldsymbol{\phi}_{i}$ into three parts:
$\boldsymbol{\phi}_{i}^{w}$, $\boldsymbol{\phi}_{i}^{h}$, and
$\boldsymbol{\phi}_{i}^{d}$. $\boldsymbol{\phi}_{i}^{w}$ and
$\boldsymbol{\phi}_{i}^{h}$ have a length of $K$, while
$\boldsymbol{\phi}_{i}^{d}$ has a length of $K-1$.

3: Apply softmax functions to $\boldsymbol{\phi}_{i}^{w}$ and
$\boldsymbol{\phi}_{i}^{h}$, followed by scaling by $2B$. These outputs
represent the widths and heights of the $K$ bins, ensuring positivity and
coverage of the interval $[-B,B]$.

4: The cumulative sums of the widths and heights yield the $K+1$ knots.

5: $\boldsymbol{\phi}_{i}^{d}$ undergoes a softplus transformation,
representing the values of the derivatives at the internal knots.

Evaluating the rational-quadratic spline transform at a given location $x$
involves efficiently determining the corresponding bin, typically achieved
through binary search due to sorted bins. The Jacobian determinant is computed
as a closed-form product of quotient derivatives, while the inversion process
entails solving a quadratic equation, the coefficients of which depend on the
target value.

Unlike simpler transformations, which offer limited flexibility, a
differentiable RQS with sufficient bin count can approximate any monotonic
function within the specified interval $[-B,B]^{2}$. Moreover, due to the
closed-form, tractable Jacobian determinant it can be analytically inverted.

###  II-F Reconstructing posterior sound fields via amortized Bayesian
inference

In order for sound field models to be useful in practice, their parameter
estimation should be feasible within reasonable time limits. Amortized
Bayesian inference aims to efficiently map the relationship between
observations and posterior distributions through a conditional density
estimator such as the trained CINN [38]. The initial cost of training the CINN
is ‘amortized’, as inference on new observations of pressure $\mathbf{p}$
becomes more efficient. This efficiency is due to the ability to achieve
posterior estimates $\mathbf{x}\sim
q_{\boldsymbol{\theta}}(\mathbf{x}\mid\mathbf{p})$ via a single forward pass
through the CINN.

Therefore, sampling from an approximate posterior involves first drawing
samples from the latent distribution
$\mathbf{z}_{0}\sim\mathcal{N}(\mathbf{z}_{0}\mid\mathbf{0},\mathbf{I})$ and
then passing them through the function
$\mathbf{g}\left(\mathbf{z};\mathbf{p},\widehat{\boldsymbol{\theta}}\right)$,
as depicted in Fig. 1. This sampling mechanism is highly efficient, typically
requiring only a fraction of a second to complete.

Since the prior over coefficients $\mathbf{x}$ follows a normal distribution,
e.g., (21), one can estimate the exact MAP coefficients $\widehat{\mathbf{x}}$
by first obtaining the latent variable $\widehat{\mathbf{z}}$ which maximizes
the posterior probability $\pi(\mathbf{x}|\mathbf{p})$. The likelihood is
given by

$\displaystyle\pi(\mathbf{p}|\mathbf{x})$
$\displaystyle=\pi_{z}\left(\mathbf{p}|\mathbf{z};\widehat{\boldsymbol{\theta}}\right)$
$\displaystyle=\frac{1}{\sqrt{2\pi\sigma_{p}^{2}}}\exp\left[-\frac{\left(\mathbf{p}-\mathbf{H}\cdot\mathbf{g}\left(\mathbf{z};\mathbf{p},\widehat{\boldsymbol{\theta}}\right)\right)^{2}}{2\sigma_{p}^{2}}\right],$
(34)

where $\sigma_{p}^{2}$ represents the variance in the measurement noise.
Therefore, The MAP coefficients are determined as the solution to the
optimization problem

$\displaystyle\widehat{\mathbf{x}}$
$\displaystyle=\mathbf{g}\left(\widehat{\mathbf{z}};\mathbf{p},\widehat{\boldsymbol{\theta}}\right)=\underset{\mathbf{x}}{\operatorname{argmax}}\left(\operatorname{log}\pi\left(\mathbf{x}|\mathbf{p}\right)\right)$
$\displaystyle=\underset{\mathbf{z}}{\operatorname{argmin}}\Bigg{[}\frac{\left\lVert\mathbf{p}-\mathbf{H}\cdot\mathbf{g}\left(\mathbf{z};\mathbf{p},\widehat{\boldsymbol{\theta}}\right)\right\rVert_{2}^{2}}{2\sigma_{p}^{2}}$
$\displaystyle+\frac{1}{2\sigma_{z}^{2}}\left\lVert\mathbf{z}\right\rVert_{2}^{2}+\frac{1}{2\sigma_{x}^{2}}\left\lVert\mathbf{g}\left(\mathbf{z};\mathbf{p},\widehat{\boldsymbol{\theta}}\right)\right\rVert_{2}^{2}\Bigg{]},$
(35)

with $\sigma_{z}=1$.

Algorithm 1 outlines the training and inference process of the CINN. The
procedure is divided into two main phases: training and inference. During the
training phase, the algorithm iteratively samples wave coefficients and noise
from their respective prior distributions, simulates pressure fields, and
updates the neural network parameters by minimizing a loss function until
convergence is achieved. The second phase describes the inference stage which,
depending on the specific goal, can either compute the MAP estimates as
described in (II-F) or generate samples from the approximate posterior
distribution (7), providing insights into the underlying uncertainty.

Algorithm 1 CINN for amortized Bayesian inference

1: Training Phase (Batch Size $D$):

2: repeat

3: for $d=1$ to $D$ do

4: Sample coefficients from prior:
$\mathbf{x}_{d}\sim\pi_{x}(\mathbf{x}_{d}|\sigma_{x})$.

5: Sample signal-to-noise ratio (SNR) from its prior:
$\text{SNR}\sim\pi(\text{SNR})$.

6: Define noise variance based on SNR:
$\sigma_{p}^{2}=\frac{\mathbb{E}\left[\|\mathbf{p}_{d}\|^{2}\right]}{\text{SNR}}$.

7: Sample noise from its distribution:
$\mathbf{n}_{d}\sim\pi_{n}(\mathbf{n}_{d}|\sigma_{p})$.

8: Simulate the pressure field:
$\mathbf{p}_{d}=\mathbf{H}\mathbf{x}_{d}+\mathbf{n}_{d}$.

9: Pass the data $\\{\mathbf{x}_{d},\mathbf{p}_{d}\\}$ through the inverse
CINN:
$\mathbf{z}_{d}=\mathbf{g}^{-1}(\mathbf{x}_{d};\mathbf{p}_{d},\boldsymbol{\theta})$.

10: Compute the loss according to (18) and obtain gradients w.r.t.
$\boldsymbol{\theta}$.

11: Update $\boldsymbol{\theta}$ using backpropagation.

12: end for

13: until convergence to $\widehat{\boldsymbol{\theta}}$

14: Inference Phase (Experimental pressure $\mathbf{p}$):

15: if estimate
$\underset{\mathbf{x}}{\text{argmax}}\;\pi(\mathbf{x}|\mathbf{p})$ then

16: Sample a latent variable:
$\mathbf{z}\sim\mathcal{N}(\mathbf{z}\mid\mathbf{0},\mathbf{I})$.

17: repeat

18: Pass $\mathbf{p}$ and $\mathbf{z}$ through the forward CINN:
$\mathbf{x}=\mathbf{g}(\mathbf{z};\mathbf{p},\widehat{\boldsymbol{\theta}})$.

19: Compute loss according to (II-F) and obtain gradients w.r.t. $\mathbf{z}$.

20: Update $\mathbf{z}$ using backpropagation

21: until convergence to $\widehat{\mathbf{x}}$

22: Return
$\widehat{\mathbf{x}}=\underset{\mathbf{x}}{\text{argmax}}\;\pi(\mathbf{x}|\mathbf{p})$.

23: else if sample from $q_{\boldsymbol{\theta}}(\mathbf{x}|\mathbf{p})$ then

24: for $l=1$ to $L$ samples do

25: Sample a latent variable:
$\mathbf{z}_{l}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$.

26: Pass $\mathbf{p}$ and $\mathbf{z}_{l}$ through the forward CINN:
$\mathbf{x}_{l}=\mathbf{g}(\mathbf{z}_{l};\mathbf{p},\widehat{\boldsymbol{\theta}})$.

27: end for

28: Return $\\{\mathbf{x}_{l}\\}_{l=1}^{L}$ as samples from the approximate
posterior $q_{\boldsymbol{\theta}}(\mathbf{x}|\mathbf{p})$.

29: end if

###  II-G Hierarchical Bayes for inference and inversion

We design a hierarchical plane wave model with sparsity-inducing priors to
compare with the proposed CINN architecture. The hierarchical nature of the
model allows for the incorporation of multiple layers of abstraction, each
layer refining and focusing on the information passed from the previous one,
thus enabling a more nuanced understanding of the sound field. The assignment
of sparsity-inducing priors is not merely a computational convenience but a
reflection of the physical reality that many sound fields are inherently
sparse in their representation, particularly at low frequencies [11, 12].

Assuming the noise ${\bf n}\sim\mathcal{N}_{c}\left({\bf
n}|\mathbf{0},\sigma_{p}^{2}\mathbf{I}\right)$ with $\sigma_{p}^{2}$
expressing variance of model, then the likelihood is a product of complex
Gaussians

$\pi(\mathbf{p}|\mathbf{x};\sigma_{p})=\prod_{m=1}^{M}\mathcal{N}_{c}\left(p_{m}|\mathbf{H}_{m}\mathbf{x},\sigma_{p}^{2}\mathbf{I}\right).$
(36)

Prior probabilities over the coefficients can be assigned hierarchically,
emphasizing the selection of weakly informative priors. In particular, a
complex normal distribution is selected to represent each heteroscedastic
complex coefficient such that

$\pi(\mathbf{x}|\boldsymbol{\sigma}_{x})=\prod_{n=1}^{N}\mathcal{N}_{c}(x_{n}|0,\sigma_{x_{n}}^{2}).$
(37)

The hyper-priors $\sigma_{x}^{2}$ are selected to be weakly informative, as
they put more probability mass on lower values, and are selected such that
there is an individual variance associated independently with every
coefficient

$\pi(\boldsymbol{\sigma}_{x}^{2};\alpha,\beta)=\prod_{n=1}^{N}\mathcal{IG}(\sigma_{x_{n}}|\alpha,\beta),$
(38)

where $\mathcal{IG}(\cdot)$ is an inverse Gamma distribution, with $\alpha$
and $\beta$ referring to the shape and scale parameters respectively. This
selection of variances also aligns with the sparse Bayesian learning
principles, such that the selection of priors is independent of the scale of
the measured data (e.g., scale-invariance) [46]. The noise in the model is
expressed via the hyper-prior $\sigma_{p}$, which in turn is estimated by the
standard deviation of the noise floor, obtained during the experimental
process. Therefore, the posterior distribution is given by

$\pi(\mathbf{x}|\mathbf{p},\boldsymbol{\sigma}_{x},\sigma_{p})=\frac{\pi(\mathbf{p}|\mathbf{x};\sigma_{p})\pi(\mathbf{x}|\boldsymbol{\sigma}_{x})\pi(\boldsymbol{\sigma}_{x}^{2};\alpha,\beta)}{\pi(\mathbf{p})}.$
(39)

The complete hierarchical model is depicted in Fig. 2. The posterior
probability (39) is often approximated by using MCMC sampling methods, [39] a
class of algorithms for numerically approximating high-dimensional integrals.

Figure 2: Probabilistic graph of hierarchical plane wave model

The maximum a posteriori (MAP) estimate of the coefficients
$\widehat{\mathbf{x}}$ is determined by

$\widehat{\mathbf{x}}\propto\operatorname{argmax}{\pi(\mathbf{p}|\mathbf{x};\sigma_{p})\pi(\mathbf{x}|\boldsymbol{\sigma}_{x})\pi(\boldsymbol{\sigma}_{x}^{2};\alpha,\beta)},$
(40)

where the evidence $\pi(\mathbf{p})$ remains constant across the unknown
coefficients.

Challenges arise when one requires the complete distribution as described by
(39), due to the computational demands of evaluating the likelihood function,
heightened by the complexity of the underlying physical system which requires
numerous evaluations of the forward model. Moreover, the high-dimensional
nature of the parameters $\mathbf{x}$ complicates the derivation of posterior
samples, underscoring the need for efficient computational strategies in
Bayesian inference. These are common issues with sampling from high-
dimensional posteriors using MCMC algorithms.

##  III Experimental dataset

Acoustic measurements were conducted in Auditorium A of the Niels Bohr
Institute in Copenhagen, Denmark, a renowned Danish scientific facility.
Largely preserved in its original state, Auditorium A has been a venue for
numerous historical lectures and pioneering talks on quantum mechanics. Figure
3 illustrates the auditorium layout, along with images depicting the
measurement process.

The auditorium is equipped with large windows on the wall on the southeast
side (at $x=0$ m), complemented by light drapery, while concrete walls enclose
the remaining sections. Wooden benches, arranged for seating during lectures,
occupy a significant portion of the auditorium’s length (spanning
$y\in[0,6.3]$ m). The floor of the auditorium is slanted, progressively
reducing the room height at each successive bench section toward the rear.
Table I provides details on the auditorium’s reverberation time and volume.
Room impulse responses (RIRs) were measured using a Universal Robots UR5
robotic arm equipped with a 1/2” omnidirectional Brüel & Kjær pressure
microphone. The measurements were conducted across a uniform cuboid grid
comprising 4913 positions (17 per dimension), with each dimension measuring
$|x|=0.5$ m, $|y|=0.5$ m, and $|z|=0.44$ m. RIRs were obtained using the
exponential sine sweep and inverse filter method, employing 5-second sweeps. A
Dynaudio BM6 2-way passive loudspeaker served as the sound source, positioned
at the lecturer’s location (front) of the auditorium, while the measurement
aperture was situated at the rear, above the last seating row.

For fitting the invertible neural network, only a subset of this data was
used, by uniformly subsampling the complete grid to 4 positions per dimension,
and also limiting the aperture size to finally obtain an array of dimensions
$\left|x\right|=0.28125$ m, $\left|y\right|=0.28125$ m,
$\left|z\right|=0.2475$ m composed of 64 microphones. The rest of the measured
RIRs were used for experimental validation.

TABLE I: Size and reverberation times in octave bands from the studied room
(Niels Bohr Institute - Auditorium A)

Volume (m3) | Reverberation Time - $T_{60}$ (s)
---|---
$134$ | 125 Hz | 250 Hz | 500 Hz | 1 kHz | 2 kHz | 4 kHz
0.572 | 0.728 | 0.797 | 0.874 | 1.03 | 1.06

Figure 3: Experimental dataset layout in the Niels Bohr Institue - Auditorium
A (top) and photographs of the auditorium and the robotic arm used for
measuring in the microphone array configuration (bottom)

##  IV Results

###  IV-A CINN training

The training process of the CINN involves conducting Bayesian hyperparameter
selection using the Optuna package [47] to ensure optimal performance. In our
experiments, we set $\sigma_{x}=0.3$ to determine the prior over the
coefficients of the sound fields, while specifying signal-to-noise ratio (SNR)
levels ranging from $\text{SNR}_{l}=25$ to $\text{SNR}_{h}=45$ for the random
wave field data. An initial learning rate of $\eta=10^{-5}$ was employed,
coupled with an exponential decay learning rate scheduler using the Adam
optimizer [48]. The training dataset (and the later experimental validation)
consisted of $N=642$ plane waves, with a batch size of 256 sound fields. The
coefficient vectors were split in half when processed by each normalizing flow
layer (i.e. $\mathrm{n}=321$). Finally, to accommodate complex-valued data, we
adopt a strategy that involves independently treating the real and imaginary
parts of a complex variable. The real and imaginary components of any
arbitrary complex variable $\mathbf{v}\in\mathbb{C}^{N_{v}}$ are concatenated
along the same axis, effectively doubling the dimensionality of the variable
such that $\mathbf{u}\in\mathbb{R}^{2N_{v}}$

The CINN was trained jointly on all frequencies of interest, which allows for
generalization across different room configurations and frequencies. Each RQS
comprised $K=8$ bins, and the CINN was composed of $N_{f}=12$ flows in total,
with random permutations incorporated between flows (except for the last
flow). The feed-forward neural networks parameterizing the spline flows were
each composed of two layers with 1024 neurons per layer and employed
Exponential Linear Unit (ELU) activation functions with a hyperparameter
$\alpha_{\text{ELU}}$ value of 1. The implementation is carried out using the
publicly available Framework for Easily Invertible Architectures Python
package [36].

Training the CINN involved $100000$ iterations on a single NVIDIA A100 80GB
GPU, utilizing only 18% of the system memory. The training duration spanned a
total of 16 hours, ensuring thorough convergence. Training validation was
conducted by evaluating the averaged normalized mean square error
($\text{NMSE}_{\text{oct}}$) across $N_{b}=6$ octave bands between 125 Hz and
4 kHz, of the 64-channel experimental dataset (see Sec. III). This evaluation
entailed sampling from the CINN at periodic intervals and subsequently
averaging the sampled coefficients, yielding the posterior mean for each
frequency, and then reconstructing the pressure at the microphone positions.
The NMSE and $\text{NMSE}_{\text{oct}}$ are given by

NMSE
$\displaystyle=\frac{||\mathbf{p}-\mathbf{p}_{\diamond}||^{2}}{||\mathbf{p}||^{2}},$
(41) $\displaystyle\text{NMSE}_{\text{oct}}$
$\displaystyle=\frac{1}{N_{b}}\sum_{n=1}^{N_{b}}\text{NMSE}_{n},$ (42)

where $\mathbf{p}$ refers to the experimental (measured) pressure and
$\mathbf{p}_{\diamond}$ to the estimated pressure, respectively. Figure 4
presents both the NMSEoct (in logarithmic scale) and the negative log-
likelihood (NLL), or the logarithm of (16), observed during the training
process. The $\text{NMSE}_{\text{oct}}$ demonstrates the CINN’s capability to
fit out-of-distribution data, in this case, that of the experimental
observations. Meanwhile, the NLL indicates that the CINN quickly learns to
model the forward propagation of plane waves when conditioned on the pressure,
even in the early stages of training using the simulated data.

Figure 4: Validation using $10\log_{10}\left(\text{NMSE}_{\text{oct}}\right)$
during CINN training (top) and negative log-likelihood (logarithm of (16)) as
training loss on simulated random wave fields (bottom)

###  IV-B Single plane wave posterior estimation with CINN

To assess the performance of the proposed invertible neural network, we
conduct an investigation focused on estimating the posterior distribution of
single-plane wave sound fields. A CINN model is trained using a dataset
consisting of sound fields generated by random plane waves emanating from 600
uniformly sampled directions within a 2D circular domain. These sound fields
are captured by an $8\times 8$ microphone array spanning a $5\times 5$ m2
aperture, with amplitudes and phases drawn from distributions outlined in Sec.
II-C. Training proceeds with the same hyperparameters detailed in Sec. IV-A,
except for the learning rate and number of iterations, adjusted to
$\eta=2\cdot 10^{-4}$ and 2500, respectively. For MAP inference, as discussed
in (II-F), we employ 1000 iterations with a learning rate of $\eta_{z}=0.05$
while utilizing the Adam optimizer. The rest of the parameters are determined
by the training data.

Following training, we subject the model to an out-of-distribution scenario,
exposing it to a simulated single-plane wave sound field characterized by
unseen amplitude, phase, and direction, yet sampled using the same microphone
array configuration. Figure 5 illustrates the true simulated sound field
alongside the reconstructed sound field (MAP estimate) and the element-wise
standard deviation of the invertible network.

The model yields promising outcomes, particularly in areas adjacent to the
aperture, where the predicted sound field closely matches the ground truth
observations. However, discrepancies between the predicted and true sound
fields become more evident with increasing distance from the measurement
aperture. As expected, this divergence coincides with heightened uncertainty
estimates, as indicated by elevated element-wise standard deviations in
regions where the model’s predictions deviate from observed data. The CINN
demonstrates adeptness in quantifying uncertainty, offering valuable insights
into areas of prediction discrepancy. While apparent in this straightforward
scenario, such findings suggest that CINNs can contribute significantly to
enhancing our comprehension of both model and data uncertainties in sound
field analysis and signal processing applications, such as sensor placement.

Figure 5: Single plane wave reference sound field, MAP inference, and point-
wise standard deviation using a CINN; the microphone positions are
superimposed over the reference sound field

###  IV-C Quantifying uncertainty in CINN sound field reconstruction

Quantifying uncertainty in sound field reconstruction using wave-based models
is crucial for ensuring the reliability of the results obtained. Given the
absence of ground truth sound field data in most practical scenarios, relying
solely on a single point estimate may obscure the validity of the solution,
underscoring the necessity of incorporating the inherent uncertainty inherent
in Bayesian models. This subsection emphasizes the importance of uncertainty
quantification in the context of CINN sound field reconstruction when using
experimental data, supported by qualitative insights.

Figure 6 illustrates extrapolated sound fields (volume of 0.11 m3) beyond the
array bounds (volume of 0.02 m3) at frequencies of 400 Hz, 800 Hz, 1200 Hz,
and 1600 Hz. It juxtaposes the ground truth with the CINN MAP prediction and
the point-wise standard deviation, as defined by (7), employing a 64-channel
microphone array detailed in Sec. III. The CINN exhibits the ability to
accurately reconstruct the sound fields, even when extrapolating to dimensions
nearly twice the size of the measurement array aperture while capturing
underlying spatial characteristics across frequencies. Notably, the CINN
demonstrates higher confidence in its extrapolated predictions at lower
frequencies, with standard deviation values increasing as frequency increases.
This observation aligns with the fact that lower frequencies with longer
wavelengths exhibit greater spatial coherence, allowing for extrapolation over
large distances relative to their wavelength. Conversely, higher frequencies
have reduced spatial coherence, leading to increased variability in
extrapolation due to the smaller scale of spatial features relative to the
wavelength. Additionally, the discernible periodic-like pattern in the
standard deviation also suggests this. As the wavelength approaches or becomes
smaller than the inter-microphone spacing, the array’s ability to
unambiguously capture and reconstruct the wavefront diminishes, evident in the
array’s sensitivity to certain directions.

Figure 6: (Color online) Measured and CINN (MAP) extrapolated sound fields
with point-wise standard deviation for frequencies $f=400,\;800,\;1200$ and
$1600$ Hz; the 64-channel validation microphone array is super-imposed on the
first (ground truth) sound field at $f=300$ Hz

###  IV-D Comparison of CINN with hierarchical Bayes MCMC for experimental
sound field estimation

The hierarchical Bayes plane wave model parameters are obtained using the
gradient-based No-U-Turn sampler (NUTS), which belongs to the class of
Hamiltonian Monte-Carlo (HMC) algorithms, chosen for its efficiency,
particularly in handling high-dimensional data [39]. The Python package Pyro
is utilized for MCMC sampling [49]. The plane wave dictionary $\mathbf{H}$ is
also composed of 642 plane waves and the hyperprior $\sigma_{p}$ is also
obtained from the noise floor in the measurements. The rest of the hyperpriors
are set to $\alpha=3$ and $\beta=1$ accordingly. Model selection was based on
the Pareto-smoothed importance sampling leave-one-out cross-validation
criteria [50] at $f=500$ Hz, with the same model used for all other
experimental cases. Two chains are generated for each model, each consisting
of 1000 warm-up samples and a further 2000 posterior samples, which are then
mixed to remove potential non-representative samples. The hierarchical Bayes
model is hereby denoted as MCMC for ease of comparison. Each hierarchical
model posterior (per frequency) requires approximately 7 minutes per chain.

On the other hand, sampling from the CINN posterior is very efficient,
typically taking only a few seconds to generate thousands of samples. This
process involves sampling from the latent distribution followed by a
straightforward forward pass through the network. Additionally, recovering the
MAP plane wave coefficients, as detailed in (II-F), involves iterating 1500
times with a learning rate of $\eta_{z}=0.05$, a task that takes less than a
minute on a GPU. Hyperparameter selection, including the aforementioned
parameters, is conducted using the Optuna package for Bayesian hyperparameter
selection, yielding fairly consistent results across minor variations in
parameter settings.

To evaluate the performance of either method we use the NMSE in decibels as
denoted in (41) and the Spatial Coherence (SC) between the measured pressure
$\mathbf{p}$ and the predicted pressure $\mathbf{p}_{\diamond}$ given by

$\textsc{SC}=\frac{\left|\mathbf{p}^{H}\mathbf{p}_{\diamond}\right|^{2}}{(\mathbf{p}^{H}\mathbf{p})(\mathbf{p}_{\diamond}^{H}\mathbf{p}_{\diamond})},$
(43)

where $H$ denotes the Hermitian transpose operation, the SC metric shows the
degree of similarity between the pressure fields with values between 0 and 1.

Figure 7 presents the NMSE and SC across frequencies for both the MCMC and
CINN models and the 64 microphone array as shown in Fig. 3. Specifically for
the CINN, results are showcased for both the MAP prediction and the posterior
mean prediction (i.e.,
$\mathbb{E}\left[q_{\boldsymbol{\theta}}(\mathbf{x}\mid\mathbf{p})\right]$).
Additionally, the corresponding metric-wise standard deviations for each
method are denoted as $\sigma_{\text{CINN}}$ and $\sigma_{\text{MCMC}}$. The
CINN MAP prediction outperforms the respective CINN posterior mean prediction,
which is expected given that the CINN models a multivariate distribution
across multiple frequencies and is likely not unimodal. Nonetheless, the CINN
yields favorable results considering its inference time is negligible, and its
adaptability to various types of sound fields.

Both the CINN MAP estimate and the MCMC MAP estimate exhibit similar
performance, with noticeable divergence occurring primarily beyond 1 kHz. This
divergence is anticipated, as the hierarchical model can often favor smaller
coefficient variances (i.e., $\left\lVert\mathbf{x}\right\rVert\approx 0$) due
to the inverse gamma prior assigned over them. Consequently, a degree of
sparsity is observed at lower frequencies, but this tendency diminishes at
higher frequencies. This behavior is also reflected in the uncertainty in the
SC metric for the MCMC method, which is more pronounced at lower frequencies
compared to higher frequencies. Employing a less informative prior may enhance
performance at higher frequencies, albeit potentially at the expense of
performance at lower frequencies.

The CINN circumvents this model selection issue by relying on data-driven
priors. Notably, both CINN MAP and posterior mean point estimates show to be
more robust at higher frequencies, reflected in the improvements of both NMSE
and SC metrics. Furthermore, the uncertainty associated with the CINN
estimates exhibits a decreasing trend in standard deviation with increasing
frequency. This trend can be attributed to the CINN’s training methodology,
which involved random wave models without any sparsity constraints.

Figure 7: Performance of CINN and MCMC MAP estimates, as well as the CINN
posterior mean and the respective 95 % confidence intervals of the NMSE (top)
and the spatial coherence (bottom), using the experimental 64 channel
microphone array

We investigate the robustness of each model concerning the decimation of the
microphone array, where the arrays undergo uniform decimation in each
dimension, reducing the total number of microphones from 216 to 27. It is
important to highlight that the aperture size is adjusted proportionally to
ensure a Nyquist rate well above the frequencies under examination.

Figure 8 depicts the performance of the posterior means for both the CINN and
MCMC methods, with their respective 95% confidence interval (approximately
$2\sigma$). In the case of MCMC, this corresponds to the MAP estimate. The
comparison reveals the comparable performance of both methods at low
frequencies in terms of NMSE, with the MCMC slightly outperforming the
proposed CINN; however, the opposite occurs at high frequencies.

The CINN demonstrates a slight decrease in NMSE with an increasing number of
microphones. Conversely, the MCMC method does not exhibit this behavior, as a
larger dataset assigns more weight to the likelihood and less to the prior,
particularly when utilizing sparsity-imposing hyperpriors. Consequently,
additional data may necessitate recalibration of the hierarchical model.

At low frequencies, the CINN saturates for a large number of microphones,
evident at $f=250$ Hz for 96 and 216 microphones. Notably, the standard
deviation, indicated by the error bars, is significantly larger than the rest
of the predictions. This disparity likely stems from variance in the
measurement noise present in the experimental data, which may differ from the
training data used for the CINN, and is more pronounced at low frequencies.

Figure 8: Perfomance of CINN and MCMC methods over frequency, each subfigure
displays the posterior mean NMSE and 95% confidence interval at different
microphone array configurations

###  IV-E Room impulse response reconstruction

Figure 9 demonstrates an example of an interpolated RIR (i.e., within the
array bounds), with the top subfigure depicting the CINN MAP prediction and
the bottom showing the posterior mean prediction within its respective 99.7%
confidence interval, both superimposed over the measured RIR. It is noteworthy
that while the MAP requires an optimization step as outlined in (II-F), its
balance of accuracy and efficiency renders it an appealing method for RIR
reconstruction. On the other hand, even if not as accurate, obtaining the
posterior mean RIR is nearly instantaneous. As mentioned, it merely
necessitates sampling from the latent distribution and performing as many
forward passes of the network as there are discrete frequencies. This is
significant as RIR reconstruction entails a heavy computational load and is a
time-consuming task, an aspect often overlooked in the literature.

Figure 9: Normalized CINN MAP estimated RIR (top) and posterior mean RIR and
99.7 % confidence interval (bottom) at an interpolated position both
superimposed over the measured reference RIR

Figure 10 displays the reconstructed MAP and posterior mean RIRs at a distance
of 0.12 m from the aperture. The MAP estimate demonstrates noteworthy
predictive capability, adeptly representing the numerous dense reflections
inherent in this dataset, with most peaks aligning well in both magnitude and
phase. Conversely, the posterior mean exhibits less precision in terms of
magnitude, as the peaks extend beyond the reference RIR when normalized.
However, it can correctly predict the alignment of individual reflections with
respect to the reference RIR.

The majority of uncertainty in either reconstructed RIR is concentrated at the
beginning of the responses, before the direct part, but also encompasses the
early reflections. This phenomenon can be attributed to pre-ringing, a
characteristic of frequency domain methods, which aligns with the hypothesis
of plane wave underrepresentation, potentially inducing this phenomenon.

Figure 10: Normalized CINN MAP estimated RIR (top) and posterior mean RIR and
99.7 % confidence interval (bottom) at an extrapolated position both
superimposed over the measured reference RIR

##  V Conclusions

A conditional invertible neural network, capable of learning the forward
operation of plane wave propagation, has been proposed for sound field
reconstruction. This approach offers several advantages over existing
methodologies. Firstly, it does not necessitate exhaustive datasets for
training, as the training data was simulated and obtained using random wave
fields. This reduces the burden of data collection and allows for more
efficient training processes.

Furthermore, the model provides interpretable uncertainty estimates,
delineating mainly uncertainties apparent in the acoustic model. This allows
for a better understanding of the reconstructed sound fields’ reliability and
the potential need for model refinement when no ground truth data is
available.

Moreover, the proposed neural network can serve as an invertible prior for MAP
estimation of plane wave coefficients, or be directly sampled to obtain
posterior mean coefficients, conditioned on the measured pressure and
parameters obtained during training. This flexibility in estimation methods
broadens the versatility of the model and enables one to choose the approach
that best suits their specific requirements.

In comparative evaluations and under challenging acoustic conditions, the
proposed neural network outperforms the hierarchical Bayes model at high
frequencies, showing improved performance in the experimental validation.
Additionally, it exhibits greater consistency than the hierarchical Bayes
model in terms of robustness to array decimation and variability of data.
Lastly, the proposed approach is well suited for RIR reconstruction, offering
either the accuracy of MAP-estimated RIRs or the rapid estimation of posterior
mean RIRs.

## References

  * [1] M. Nolan, S. A. Verburg, J. Brunskog, and E. Fernandez-Grande, “Experimental characterization of the sound field in a reverberation rooma),” _The Journal of the Acoustical Society of America_ , vol. 145, no. 4, pp. 2237–2246, Apr. 2019.
  * [2] M. Berzborn and M. Vorländer, “Directional sound field decay analysis in performance spaces,” _Building Acoustics_ , vol. 28, no. 3, pp. 249–263, Sep. 2021.
  * [3] T. Deppisch, S. V. A. Garí, P. Calamia, and J. Ahrens, “Direct and Residual Subspace Decomposition of Spatial Room Impulse Responses,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 31, pp. 927–942, 2023.
  * [4] L. Birnie, T. Abhayapala, V. Tourbabin, and P. Samarasinghe, “Mixed Source Sound Field Translation for Virtual Binaural Application With Perceptual Validation,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 29, pp. 1188–1203, 2021.
  * [5] F. M. Heuchel, D. Caviedes-Nozal, J. Brunskog, F. T. Agerkvist, and E. Fernandez-Grande, “Large-scale outdoor sound field control,” _The Journal of the Acoustical Society of America_ , vol. 148, no. 4, pp. 2392–2402, Oct. 2020.
  * [6] T. Abe, S. Koyama, N. Ueno, and H. Saruwatari, “Amplitude Matching for Multizone Sound Field Control,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 31, pp. 656–669, 2023.
  * [7] L. Vindrola, M. Melon, J.-C. Chamard, and B. Gazengel, “Use of the filtered-x least-mean-squares algorithm to adapt personal sound zones in a car cabin,” _The Journal of the Acoustical Society of America_ , vol. 150, no. 3, pp. 1779–1793, Sep. 2021.
  * [8] J. Cadavid, M. B. Møller, C. S. Pedersen, S. Bech, T. van Waterschoot, and J. Østergaard, “Performance of low-frequency sound zones with very fast room impulse response measurements,” _The Journal of the Acoustical Society of America_ , vol. 155, no. 1, pp. 757–768, 01 2024.
  * [9] E. Erdem, Z. Cvetković, and H. Hacıhabiboğlu, “3D Perceptual Soundfield Reconstruction via Virtual Microphone Synthesis,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 31, pp. 1305–1317, 2023.
  * [10] H. Zuo, T. D. Abhayapala, and P. N. Samarasinghe, “Particle Velocity Assisted Three Dimensional Sound Field Reproduction Using a Modal-Domain Approach,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 28, pp. 2119–2133, 2020.
  * [11] R. Mignot, G. Chardon, and L. Daudet, “Low frequency interpolation of room impulse responses using compressed sensing,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 22, no. 1, pp. 205–216, 2013.
  * [12] S. A. Verburg and E. Fernandez-Grande, “Reconstruction of the sound field in a room using compressive sensing,” _The Journal of the Acoustical Society of America_ , vol. 143, no. 6, pp. 3770–3779, 2018.
  * [13] S. Koyama and L. Daudet, “Sparse Representation of a Spatial Sound Field in a Reverberant Environment,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 13, no. 1, pp. 172–184, Mar. 2019.
  * [14] N. Antonello, E. De Sena, M. Moonen, P. A. Naylor, and T. Van Waterschoot, “Room impulse response interpolation using a sparse spatio-temporal representation of the sound field,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 25, no. 10, pp. 1929–1941, 2017.
  * [15] J. M. Schmid, E. Fernandez-Grande, M. Hahmann, C. Gurbuz, M. Eser, and S. Marburg, “Spatial reconstruction of the sound field in a room in the modal frequency range using Bayesian inference,” _The Journal of the Acoustical Society of America_ , vol. 150, no. 6, pp. 4385–4394, Dec. 2021.
  * [16] N. Ueno, S. Koyama, and H. Saruwatari, “Kernel ridge regression with constraint of helmholtz equation for sound field interpolation,” in _2018 16th International Workshop on Acoustic Signal Enhancement (IWAENC)_. IEEE, 2018, pp. 1–440.
  * [17] D. Caviedes-Nozal, N. A. Riis, F. M. Heuchel, J. Brunskog, P. Gerstoft, and E. Fernandez-Grande, “Gaussian processes for sound field reconstruction,” _The Journal of the Acoustical Society of America_ , vol. 149, no. 2, pp. 1107–1119, 2021.
  * [18] J. G. C. Ribeiro, S. Koyama, and H. Saruwatari, “Kernel Interpolation of Acoustic Transfer Functions with Adaptive Kernel for Directed and Residual Reverberations,” in _ICASSP 2023 - 2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , Jun. 2023, pp. 1–5, iSSN: 2379-190X.
  * [19] M. Hahmann, S. A. Verburg, and E. Fernandez-Grande, “Spatial reconstruction of sound fields using local and data-driven functions,” _The Journal of the Acoustical Society of America_ , vol. 150, no. 6, pp. 4417–4428, 2021.
  * [20] M. Pezzoli, M. Cobos, F. Antonacci, and A. Sarti, “Sparsity-Based Sound Field Separation in the Spherical Harmonics Domain,” in _ICASSP 2022 - 2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , May 2022, pp. 1051–1055, iSSN: 2379-190X.
  * [21] E. Fernandez-Grande, D. Caviedes-Nozal, M. Hahmann, X. Karakonstantis, and S. A. Verburg, “Reconstruction of room impulse responses over extended domains for navigable sound field reproduction,” in _2021 Immersive and 3D Audio: from Architecture to Automotive (I3DA)_. IEEE, 2021, pp. 1–8.
  * [22] E. G. Williams, _Fourier acoustics: sound radiation and nearfield acoustical holography_. Elsevier, 1999.
  * [23] M. Cobos, J. Ahrens, K. Kowalczyk, and A. Politis, “An overview of machine learning and other data-based methods for spatial audio capture, processing, and reproduction,” _EURASIP Journal on Audio, Speech, and Music Processing_ , vol. 2022, no. 1, p. 10, May 2022.
  * [24] G. Götz, S. J. Schlecht, and V. Pulkki, “Common-Slope Modeling of Late Reverberation,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 31, pp. 3945–3957, 2023.
  * [25] E. Zea, E. Brandão, M. Nolan, J. Cuenca, J. Andén, and U. P. Svensson, “Sound absorption estimation of finite porous samples with deep residual learning,” _The Journal of the Acoustical Society of America_ , vol. 154, no. 4, pp. 2321–2332, 10 2023.
  * [26] F. Lluis, P. Martinez-Nuevo, M. Bo Møller, and S. Ewan Shepstone, “Sound field reconstruction in rooms: Inpainting meets super-resolution,” _The Journal of the Acoustical Society of America_ , vol. 148, no. 2, pp. 649–659, 2020.
  * [27] M. Pezzoli, D. Perini, A. Bernardini, F. Borra, F. Antonacci, and A. Sarti, “Deep prior approach for room impulse response reconstruction,” _Sensors_ , vol. 22, no. 7, p. 2710, Jan. 2022.
  * [28] X. Karakonstantis and E. Fernandez-Grande, “Generative adversarial networks with physical sound field priors,” _The Journal of the Acoustical Society of America_ , vol. 154, no. 2, pp. 1226–1238, 2023.
  * [29] F. Miotello, L. Comanducci, M. Pezzoli, A. Bernardini, F. Antonacci, and A. Sarti, “Reconstruction of Sound Field through Diffusion Models,” _arXiv preprint arXiv:2312.08821_ , Feb. 2024.
  * [30] K. Shigemi, S. Koyama, T. Nakamura, and H. Saruwatari, “Physics-informed convolutional neural network with bicubic spline interpolation for sound field estimation,” in _2022 International Workshop on Acoustic Signal Enhancement (IWAENC)_. IEEE, 2022, pp. 1–5.
  * [31] X. Karakonstantis, D. Caviedes-Nozal, A. Richard, and E. Fernandez-Grande, “Room impulse response reconstruction with physics-informed deep learning,” _The Journal of the Acoustical Society of America_ , vol. 155, no. 2, pp. 1048–1059, Feb. 2024.
  * [32] M. Olivieri, M. Pezzoli, F. Antonacci, and A. Sarti, “A physics-informed neural network approach for nearfield acoustic holography,” _Sensors_ , vol. 21, no. 23, p. 7834, 2021.
  * [33] G. Papamakarios, E. Nalisnick, D. J. Rezende, S. Mohamed, and B. Lakshminarayanan, “Normalizing flows for probabilistic modeling and inference,” _The Journal of Machine Learning Research_ , vol. 22, no. 1, pp. 57:2617–57:2680, Jan. 2021.
  * [34] L. Ardizzone, C. Lüth, J. Kruse, C. Rother, and U. Köthe, “Guided image generation with conditional invertible neural networks,” _arXiv preprint arXiv:1907.02392_ , 2019.
  * [35] G. A. Padmanabha and N. Zabaras, “Solving inverse problems using conditional invertible neural networks,” _Journal of Computational Physics_ , vol. 433, p. 110194, 2021.
  * [36] L. Ardizzone, T. Bungert, F. Draxler, U. Köthe, J. Kruse, R. Schmier, and P. Sorrenson, “Framework for Easily Invertible Architectures (FrEIA),” 2018-2022. [Online]. Available: https://github.com/vislearn/FrEIA
  * [37] S. T. Radev, U. K. Mertens, A. Voss, L. Ardizzone, and U. Köthe, “BayesFlow: Learning Complex Stochastic Models With Invertible Neural Networks,” _IEEE Transactions on Neural Networks and Learning Systems_ , vol. 33, no. 4, pp. 1452–1466, Apr. 2022.
  * [38] K. Cranmer, J. Brehmer, and G. Louppe, “The frontier of simulation-based inference,” _Proceedings of the National Academy of Sciences_ , vol. 117, no. 48, pp. 30 055–30 062, 2020.
  * [39] J. B. Rubin, H. S. Carlin, D. B. Stern, and A. Gelman, _Bayesian Data Analysis_. New York: Chapman and Hall/CRC, Jun. 1995.
  * [40] D. P. Kingma and P. Dhariwal, “Glow: Generative Flow with Invertible 1x1 Convolutions,” in _Advances in Neural Information Processing Systems_ , vol. 31. Curran Associates, Inc., 2018.
  * [41] E. Fernandez-Grande, X. Karakonstantis, D. Caviedes-Nozal, and P. Gerstoft, “Generative models for sound field reconstruction,” _The Journal of the Acoustical Society of America_ , vol. 153, no. 2, pp. 1179–1190, 2023.
  * [42] I. Kobyzev, S. D. Prince, and M. A. Brubaker, “Normalizing flows: An introduction and review of current methods,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 43, no. 11, pp. 3964–3979, nov 2021.
  * [43] L. Dinh, J. Sohl-Dickstein, and S. Bengio, “Density estimation using Real NVP,” _arXiv preprint arXiv:1605.08803_ , Feb. 2017, arXiv:1605.08803 [cs, stat] type: article.
  * [44] L. Dinh, D. Krueger, and Y. Bengio, “NICE: Non-linear Independent Components Estimation,” _arXiv preprint arXiv:1410.8516_ , Apr. 2015, arXiv:1410.8516 [cs] type: article.
  * [45] C. Durkan, A. Bekasov, I. Murray, and G. Papamakarios, “Neural spline flows,” in _Advances in Neural Information Processing Systems_ , H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, Eds., vol. 32. Curran Associates, Inc., 2019.
  * [46] M. E. Tipping, “Sparse bayesian learning and the relevance vector machine,” _Journal of machine learning research_ , vol. 1, no. Jun, pp. 211–244, 2001.
  * [47] T. Akiba, S. Sano, T. Yanase, T. Ohta, and M. Koyama, “Optuna: A Next-generation Hyperparameter Optimization Framework,” _arXiv preprint arXiv:1907.10902_ , Jul. 2019, arXiv:1907.10902 [cs, stat] type: article.
  * [48] D. P. Kingma and J. Ba, “Adam: A Method for Stochastic Optimization,” _arXiv preprint arXiv:1412.6980_ , Jan. 2017, arXiv:1412.6980 [cs] type: article.
  * [49] E. Bingham, J. P. Chen, M. Jankowiak, F. Obermeyer, N. Pradhan, T. Karaletsos, R. Singh, P. Szerlip, P. Horsfall, and N. D. Goodman, “Pyro: Deep universal probabilistic programming,” _The Journal of Machine Learning Research_ , vol. 20, no. 1, pp. 973–978, 2019.
  * [50] A. Vehtari, D. Simpson, A. Gelman, Y. Yao, and J. Gabry, “Pareto Smoothed Importance Sampling,” _arXiv preprint arXiv:1507.02646_ , Aug. 2022.

##  Appendix A Mapping variables with rational quadratic splines

The splines defined by [45] employ $K$ distinct rational-quadratic functions
(bins), delimited by $K+1$ coordinates known as knots
$\left\\{\left(z^{(k)},y^{(k)}\right)\right\\}_{k=0}^{K}$. These knots
increase monotonically from $\left(z^{(0)},y^{(0)}\right)=(-B,-B)$ to
$\left(z^{(K)},y^{(K)}\right)=(B,B)$. Internally the spline derivatives take
on $K-1$ arbitrary positive values $\left\\{\phi^{(k)}\right\\}_{k=1}^{K-1}$,
while the boundary derivatives are given by $\phi^{(0)}=\phi^{(K)}=1$ which
act as linear ‘tails’.

The method constructs a monotonic, continuously-differentiable, rational-
quadratic spline passing through knots with prescribed derivatives. The
rational-quadratic function $h^{(k)}$ within the $k^{\text{th}}$ bin is
expressed as

$\displaystyle h(z^{(k)})$
$\displaystyle=\frac{\alpha^{(k)}(\xi)}{\beta^{(k)}(\xi)}$
$\displaystyle=y^{(k)}+\frac{\left(y^{(k+1)}-y^{(k)}\right)\left[s^{(k)}\xi^{2}+\phi^{(k)}\xi(1-\xi)\right]}{s^{(k)}+\left[\phi^{(k+1)}+\phi^{(k)}-2s^{(k)}\right]\xi(1-\xi)},$
(44)

with $s_{k}=$ $\left(y^{k+1}-y^{k}\right)/\left(z^{k+1}-z^{k}\right)$
corresponding to the slopes of the lines connecting the coordinates and
$\xi(z)=\left(z-z^{k}\right)/\left(z^{k+1}-z^{k}\right)$ represents the
normalized distance along the $z$ axis within the $k^{\text{th}}$ bin.

The logarithm of the determinant of its Jacobian can be computed by summing
the logarithm of derivatives of (A) with respect to each transformed $z$ value
of the input vector. This derivative is given by

$\displaystyle\frac{\mathrm{d}h}{\mathrm{d}z}$
$\displaystyle=\frac{\mathrm{d}}{\mathrm{d}z}\left[\frac{\alpha^{(k)}(\xi)}{\beta^{(k)}(\xi)}\right]$
$\displaystyle=\frac{\left(s^{(k)}\right)^{2}\left[\phi^{(k+1)}\xi^{2}+2s^{(k)}\xi(1-\xi)+\phi^{(k)}(1-\xi)^{2}\right]}{\left[s^{(k)}+\left[\phi^{(k+1)}+\phi^{(k)}-2s^{(k)}\right]\xi(1-\xi)\right]^{2}}.$
(45)

Analytical computation of the inverse rational-quadratic function involves
inverting (A), which in turn leads to solving the following quadratic equation

$q(z)=\alpha(\xi(z))-y\beta(\xi(z))=a\xi(z)^{2}+b\xi(z)+c=0.$ (46)

Due to monotonicity, only one root is obtained. The solution is given by
$\xi(z)=2c/\left(-b-\sqrt{b^{2}-4ac}\right)$, where

$\displaystyle
a=\left(y^{(k+1)}-y^{(k)}\right)\left[s^{(k)}-\phi^{(k)}\right]$
$\displaystyle\quad+\left(y-y^{(k)}\right)\left[\phi^{(k+1)}+\phi^{(k)}-2s^{(k)}\right]$
(47) $\displaystyle b=\left(y^{(k+1)}-y^{(k)}\right)\phi^{(k)}$
$\displaystyle\quad-\left(y-y^{(k)}\right)\left[\phi^{(k+1)}+\phi^{(k)}-2s^{(k)}\right]$
(48) $\displaystyle c=-s^{(k)}\left(y-y^{(k)}\right)$ (49)

where the coefficients $a$, $b$, and $c$ depend on the target output $y$.
These expressions allow us to determine $\xi(z)$, subsequently enabling the
computation of the inverse transformation of $h$.