On the Approximation of Bi-Lipschitz Maps by Invertible Neural Networks
Invertible neural networks (INNs) represent an important class of deep neural
network architectures that have been widely used in several applications. The
universal approximation properties of INNs have also been established recently.
However, the approximation rate of INNs is largely missing. In this work, we
provide an analysis of the capacity of a class of coupling-based INNs to
approximate bi-Lipschitz continuous mappings on a compact domain, and the
result shows that it can well approximate both forward and inverse maps
simultaneously. Furthermore, we develop an approach for approximating
bi-Lipschitz maps on infinite-dimensional spaces that simultaneously
approximate the forward and inverse maps, by combining model reduction with
principal component analysis and INNs for approximating the reduced map, and we
analyze the overall approximation error of the approach. Preliminary numerical
results show the feasibility of the approach for approximating the solution
operator for parameterized second-order elliptic problems.
Bangti Jin, Zehui Zhou, Jun Zou
# On the Approximation of Bi-Lipschitz Maps by Invertible Neural Networks+

###### Abstract

Invertible neural networks (INNs) represent an important class of deep neural network architectures that have been widely used in several applications. The universal approximation properties of INNs have also been established recently. However, the approximation rate of INNs is largely missing. In this work, we provide an analysis of the capacity of a class of coupling-based INNs to approximate bi-Lipschitz continuous mappings on a compact domain, and the result shows that it can well approximate both forward and inverse maps simultaneously. Furthermore, we develop an approach for approximating bi-Lipschitz maps on infinite-dimensional spaces that simultaneously approximate the forward and inverse maps, by combining model reduction with principal component analysis and INNs for approximating the reduced map, and we analyze the overall approximation error of the approach. Preliminary numerical results show the feasibility of the approach for approximating the solution operator for parameterized second-order elliptic problems.

**Keywords**: invertible neural network; bi-Lipschitz map; error estimate; operator approximation.

## 1 Introduction

Invertible neural networks (INNs) are a class of neural network (NN) architectures with invertibility by design, via special invertible layers called the flow layers. INNs often enjoy tractable numerical algorithms to compute the inverse map and Jacobian determinant, e.g., with explicit inversion formulas. These distinct features have made them very attractive for a variety of machine learning tasks, e.g., generative modeling [16, 31, 29], probabilistic modeling [38, 17, 23, 6], solving inverse problems [2, 1, 3], modeling nonlinear dynamics [9] and point cloud generation [44].

There are several different classes of INNs, including invertible residual networks (iResNet) [7, 43], neural ordinary differential equations (NODEs) [11, 13, 18] and coupling-based neural networks [16, 17, 25, 31, 2]. For iResNet, Behrmann et al [7] leveraged the viewpoint of ResNets as an Euler discretization of ODEs and proved the standard ResNet architecture can be made invertible by adding a simple normalization step to control the Lipschitz constant of the NN during training. The inverse is not available in closed form but can be obtained through a fixed-point iteration. Chen et al [13] proposed using black-box ODE solvers as a model component, and developed a class of new models, i.e., NODEs, for time-series modeling, supervised learning, and density estimation etc. NODEs indirectly models an invertible function by transforming an input vector through an ordinary differential equation (ODE). Dupont and Doucet [18] introduced a class of more expressive and empirically stable models, augmented neural ODEs (ANODEs), which have a lower computational cost. There is another family of popular INNs using coupling-based frameworks, including non-linear independent component estimation (NICE) [16], real-valued non-volume preserving (RealNVPs) [17], invertible reversible residual network
(i-RevNet) [25], and generative flow (Glow) [31]. In coupling-based NNs, one employs a highly restricted NN architecture where only some of the input variables undergo some transformations, and the rest of the input variables become the output as-is so that the splitting facilitates the inversion process.

The great empirical successes of INNs in diverse applications have sparked intensive interest in theoretical study of INNs, e.g., approximation theory and training. The construction of INNs imposes certain restrictions on the choice of admissible layers, and hence their approximation properties are more limited and far less well understood than that for standard fully connected NNs (see [15] for an overview of the latter). Indeed, only very recently, Teschima et al [41] proved a first universal approximation property in \(\mathcal{L}^{p}\) spaces for coupling based INNs. Lyu et al [39] extended the results to \(C^{k}\) universality for \(C^{k}\) diffeomorphism and also established the parametric version. More recently, Ishikawa et al [24] developed a more general theoretical framework, covering also other types of INNs, e.g., NODEs. However, approximation rates of INNs appear unknown so far.

In this work, we continue the line of research on the approximation theory for INNs. We explicitly construct a coupling-based INN for simultaneously approximating the forward and inverse maps of a bi-Lipschitz map on a compact domain, by combining a shallow NN with a deep coupling-based NN, and provide an analysis of the approximation capacity of the constructed INNs. We show that quantitatively it can approximate the forward and inverse processes simultaneously to a certain tolerance; see Theorem 2.1 for the precise statement. This result complements the existing results [41, 39, 24] by providing an explicit construction and also quantitative estimates. To the best of our knowledge, this result is the first of the kind in the literature. Methodologically, the analysis is inspired by the techniques for establishing approximation capabilities of INNs / NODEs [4, 35, 41].

To illustrate the use of the approximation result, we apply INNs to approximate bi-Lipschitz maps in infinite-dimensional spaces, which belongs to operator learning and represents a fundamental task in the study of many scientific problems [33]. We focus on the situation where the ambient dimensionality is high but the intrinsic dimensionality is low, so that model reduction can be applied effectively. There are many methods for feature projection to reduce the dimensionality, e.g., principal component analysis (PCA) [22, 8], nonnegative matrix factorization [34] and Johnson-Lindenstrauss embedding [27]. We propose to couple INN with PCA, i.e., first applying PCA encoder and decoder to reduced the bi-Lipschitz map in infinite-dimensional spaces to finite-dimensional subspaces, and then approximating the reduced map with INNs. We provide an analysis of the overall approximation error. Further, we provide preliminary numerical illustrations with model problems in parametric PDEs.

The idea of coupling of DNNs with model reduction has been explored in several works (see, e.g., [42, 28, 40, 10, 36, 12]). Wojcik and Kurdziel [42] proposed using Johnson-Lindenstrauss embedding to reduce the input dimensionality, and Karnik et al [28] developed complexity bounds for approximating a Holder (or uniformly) continuous function under a low-complexity structure assumption. Bhattacharya et al [10] proposed an approach using PCA for dimensionality reduction, and then employing a fully connected feedforward NN to approximate the reduced nonlinear mapping. Pineda and Petersen [40] investigated restricting the domain of the forward operator to finite-dimensional spaces, and proved the existence of a robust-to-noise NN approximation of the operator. Our work continues the active line of research on developing operator approximations, by extending the encoder-decoder type approach in [10] to simultaneously approximating the forward and inverse maps using INNs. Liu et al [36, 12] thoroughly analyzed general encoder-decoder based approaches in terms of approximation and generalization errors, and applied the results to canonical PDEs. In this work, we use the approximation result of INNs developed herein instead of more well-understood feedforward fully connected NNs, in order to derive the estimation error on the INN based operator approximations.

The rest of the paper is organized as follows. In Section 2, we develop an INN architecture for approximating the bi-Lipschitz map, and analyze the approximation error of the constructed INN. In Section 3, we develop a novel operator approximation based on model reduction (PCA) and INN, and discuss the error estimation for the overall scheme. Finally, in Section 4 we present some preliminary numerical experiments to illustrate the approach.
## 2 Invertible neural networks

In this section, we construct a coupling-based invertible neural network (INN) to approximate the bi-Lipschitz forward operator \(F:\mathbb{R}^{d}\to\mathbb{R}^{d}\), \(d\geq 2\).

**Assumption 2.1**.: _The operator \(F\) is bi-Lipschitz continuous with Lipschitz constants \(L_{F}\) and \(L_{F^{-1}}\) for \(F\) and its inverse \(F^{-1}\), respectively._

Coupling-based NNs split the input \(x\in\mathbb{R}^{d}\) into two parts and then employ an affine coupling layer whose inverse can be found directly. Specifically, fix \(k\in\mathbb{N}\) with \(k<d\), we split \(x\in\mathbb{R}^{d}\) into \(x=(x_{\leq k},x_{>k})\), with the subvectors \(x_{\leq k}:=(x_{1},\ldots,x_{k})\in\mathbb{R}^{k}\) and \(x_{>k}:=(x_{k+1},\ldots,x_{d})\in\mathbb{R}^{d-k}\), and set

\[x_{\leq k}^{\prime}= x_{\leq k}\odot\exp(s_{1}(x_{>k}))+t_{1}(x_{>k}),\] \[x_{>k}^{\prime}= x_{>k}\odot\exp(s_{2}(x_{\leq k}^{\prime}))+t_{2}(x_{\leq k}^{ \prime}),\]

with operators \(s_{1},t_{1}:\mathbb{R}^{d-k}\longrightarrow\mathbb{R}^{k}\) and \(s_{2},t_{2}:\mathbb{R}^{k}\longrightarrow\mathbb{R}^{d-k}\), which are all taken to be NNs. The inverse of this process can be obtained directly as

\[x_{>k}= (x_{>k}^{\prime}-t_{2}(x_{\leq k}^{\prime}))\odot\exp(-s_{2}(x_{ \leq k}^{\prime})),\] \[x_{\leq k}= (x_{\leq k}^{\prime}-t_{1}(x_{>k}))\odot\exp(-s_{1}(x_{>k})).\]

Intuitively, it reduces the spatial resolution of the coefficients and potentially increasing the layer size may improve the expressivity of the resulting NN. Numerically, it has been shown that the restriction to the affine coupling layers will not noticeably reduce the expressivity of the NN [2]. The universal approximation properties of (affine) coupling based INNs were studied in [41, 39, 24]. We provide an explicit construction of INNs to approximate bi-Lipschitz maps, which yields an explicit bound on the INN approximation.

### Construction

Now we construct an INN, denoted by \(F_{\mathrm{nn}}\), to approximate the bi-Lipschitz map \(F:\mathbb{R}^{d}\to\mathbb{R}^{d}\), over any compact set \(K\subset\mathbb{R}^{d}\). Throughout we assume that \(K=[0,1]^{d}\) and \(\{x^{i}\}_{i=1}^{N}\) are taken from a uniform grid in \(K\). For a general compact set \(K\subset\mathbb{R}^{d}\), we can always find an affine mapping \(A_{f}:\;x\to W_{a}x+b_{a}\) with \(W_{a}\in\mathbb{R}^{d\times d}\) and \(b_{a}\in\mathbb{R}^{d}\) such that \(A_{f}(K)\subset[0,1]^{d}\). Let \(\{x^{i}\}_{i=1}^{N}=\{x^{\alpha}\}_{\alpha\in[n-1]^{d}}\), where \([n]\) denotes the set \([n]=\{0,\cdots,n\}\), \(N=n^{d}\) and \(x^{\alpha}=\frac{\alpha}{n}\in\mathbb{R}^{d}\). We also denote the \(j\)th entry of a vector \(x\in\mathbb{R}^{d}\) by \(x_{j}\). For a given set \(K\), we denote by \(|K|\) its Lebesgue measure. We denote the standard Euclidean norm of a vector (and also the spectral norm of a matrix) by \(\|\cdot\|_{2}\). For any \(f:K\to\mathbb{R}^{d}\), we define

\[\|f\|_{\mathcal{L}^{p}(K)}= \left\{\begin{aligned} \Big{(}\int_{K}\|f(x)\|_{2}^{p} \mathrm{d}x\Big{)}^{\frac{1}{p}},& 1\leq p<\infty,\\ \sup_{x\in K}\|f(x)\|_{2},& p=\infty.\end{aligned}\right.\]

Following the idea in [35, Section 4.3], we first approximate the map \(F\) by a composition mapping without using NNs. This is achieved in the following lemma. The function \(h^{r}\) (depending on the choice of \(n\) and \(r\)) serves a piecewise linear approximation to the following piecewise constant function:

\[\lim_{r\to 1^{-}}h^{r}(x)=\left\{\begin{array}{ll}\frac{j}{n},&x\in[ \frac{j}{n},\frac{j+1}{n}),\;j\in[n-1],\\ x,&x\in\mathbb{R}\setminus[0,1).\end{aligned}\right.\]

For any \(r\in(0,1)\), \(h^{r}\) is strictly increasing, and can be realized by a rectified linear unit (ReLU) NN exactly; see Section 2.1.1. Then by choosing \(r\in(0,1)\) sufficiently close to \(1\), we obtain the final estimate.
**Lemma 2.1**.: _For any \(r\in(0,1)\), let \(H^{r}:\mathbb{R}^{d}\to\mathbb{R}^{d}\) be defined by \(H^{r}(x)=(h^{r}(x_{1}),\cdots,h^{r}(x_{d}))^{t}\), with the component function \(h^{r}:\mathbb{R}\to\mathbb{R}\) given by_

\[h^{r}(x)=\left\{\begin{array}{ll}\frac{j}{n},&x=\frac{j}{n},\,j \in[n-1],\\ \frac{j+(1-r)}{n},&x=\frac{j+r}{n},\,j\in[n-1],\\ x,&x\in\mathbb{R}\setminus[0,1],\end{array}\right.\]

_being continuous over \(\mathbb{R}\) and linear in the intervals \([\frac{j}{n},\frac{j+r}{n}]\) and \([\frac{j+r}{n},\frac{j+1}{n}]\) for \(j\in[n-1]\). Then the following estimates hold_

\[\|F\circ H^{r}-F\|_{\mathcal{L}^{2}(K)}^{2} \leq L_{F}^{2}(2r-1)^{2}dn^{-2},\] \[\|H^{r}\circ F^{-1}-F^{-1}\|_{\mathcal{L}^{2}(F(K))}^{2} \leq L_{F}^{d}(2r-1)^{2}dn^{-2}.\]

Proof.: First, since \(|K|=1\), we obtain the following estimate

\[\|F\circ H^{r}-F\|_{\mathcal{L}^{2}(K)}^{2}= \|F\circ H^{r}-F\circ\mathrm{id}\|_{\mathcal{L}^{2}(K)}^{2}=\|F \circ(H^{r}-\mathrm{id})\|_{\mathcal{L}^{2}(K)}^{2}\] \[\leq |K|\|F\circ(H^{r}-\mathrm{id})\|_{\mathcal{L}^{\infty}(K)}^{2}= \|F\circ(H^{r}-\mathrm{id})\|_{\mathcal{L}^{\infty}(K)}^{2}.\]

Then, by the definitions of the functions \(H^{r}\) and \(h^{r}\), we can derive

\[\|F\circ(H^{r}-\mathrm{id})\|_{\mathcal{L}^{\infty}(K)}\leq\omega_{F}(\|H^{r} -\mathrm{id}\|_{\mathcal{L}^{\infty}(K)})\leq\omega_{F}(d^{\frac{1}{2}}\max_{ x\in[0,n^{-1}]}(h^{r}(x)-x))=\omega_{F}(|2r-1|d^{\frac{1}{2}}n^{-1}),\]

where \(\omega_{F}\) denotes the modulus of continuity of the mapping \(F\), defined by \(\omega_{F}(s)=\sup_{\|x-\tilde{x}\|_{2}\leq s}\|F(x)-F(\tilde{x})\|_{2}\leq L _{F}s\). Thus, we obtain

\[\|F\circ H^{r}-F\|_{\mathcal{L}^{2}(K)}^{2}\leq L_{F}^{2}(2r-1)^{2}dn^{-2}.\]

Similarly, by the inequality \(|F(K)|\leq L_{F}^{d}|K|=L_{F}^{d}\), we can bound \(F^{-1}-H^{r}\circ F^{-1}\). 

**Remark 2.1**.: _The estimates in Lemma 2.1 are derived in the \(\mathcal{L}^{2}\) norm, and the same results can be extended to the \(\mathcal{L}^{p}\) space, with \(1\leq p\leq\infty\). By the construction of \(h^{r}\), \(\lim_{r\to 1^{-}}h^{r}\) is a piecewise constant function. Since \(F\in\mathcal{L}^{2}(K)\), we can use piecewise constant mappings to approximate \(F\) in \(K\). The definition of \(H^{r}\) implies that \(\lim_{r\to 1^{-}}F\circ H^{r}\) may serve as such a piecewise constant mapping._

Next we construct two INNs \(H^{r}_{\mathrm{nn}}\) and \(\tilde{F}_{\mathrm{nn}}\) to approximate \(H^{r}\) and \(F\), respectively. The detailed procedures are given in the following two parts. Throughout, for any \(s\in\mathbb{N}\), the notation \(\mathbf{1}_{s}\) and \(\mathbf{0}_{s}\) denote a length-\(s\) column / row vector of ones and zeros, respectively. The notation \(O_{d\times}\) denotes a zero matrix of size \(d\times s\).

#### 2.1.1 INN for \(H^{r}\)

First we construct an INN \(H^{r}_{\mathrm{nn}}\) to realize the map \(H^{r}\). For any \(r\in(0,1)\), the function \(h^{r}:\mathbb{R}\to\mathbb{R}\) is strictly increasing and hence bijective and continuous piecewise linear with the non-differentiable points given by \((2n+1)\) tuples \(\{(a_{i},a_{i})\}_{i=0}^{n}\cup\{(p_{i},q_{i})\}_{i=1}^{n}\), with \(a_{i}=\frac{i}{n}\), \(p_{i}=\frac{i-1+r}{n}\) and \(q_{i}=\frac{i-r}{n}\). It can be exactly represented by a two-layer NN with the ReLU activation function \(\sigma(t)\)[4] (see also [21]). Specifically, let \(\sigma(t)=\max(0,t)\). Then it follows from the elementary identities

\[\frac{a_{i}-q_{i}}{a_{i}-p_{i}}-\frac{q_{i}-a_{i-1}}{p_{i}-a_{i-1 }}=\frac{\frac{i}{n}-\frac{i-r}{n}}{\frac{i}{n}-\frac{i-1+r}{n}}-\frac{\frac{ i-r}{n}-\frac{i-1}{n}}{\frac{i-1+r}{n}-\frac{i-1}{n}}=\frac{r}{1-r}-\frac{1-r}{r}= \frac{2r-1}{r(1-r)},\] \[\frac{q_{i+1}-a_{i}}{p_{i+1}-a_{i}}-\frac{a_{i}-q_{i}}{a_{i}-p_{i }}=\frac{\frac{i+1-r}{n}-\frac{i}{n}}{\frac{i+r}{n}-\frac{i}{n}}-\frac{\frac{i}{ n}-\frac{i-r}{n}}{\frac{i-1+r}{n}}=\frac{1-r}{r}-\frac{r}{1-r}=-\frac{2r-1}{r(1-r)},\]

for \(i=1,\ldots,n-1\), that we can express \(h^{r}\) as

\[h^{r}(x)= -\sigma(-x)+\frac{q_{1}-a_{0}}{p_{1}-a_{0}}\sigma(x-a_{0})+\Big{(} \frac{a_{1}-q_{1}}{a_{1}-p_{1}}-\frac{q_{1}-a_{0}}{p_{1}-a_{0}}\Big{)}\sigma(x- p_{1})+\Big{(}\frac{q_{2}-a_{1}}{p_{2}-a_{1}}-\frac{a_{1}-q_{1}}{a_{1}-p_{1}} \Big{)}\sigma(x-a_{1})\]
\[\|\tilde{h}^{r}-h^{r}\|_{\mathcal{F}^{\infty}([0,1])}\leq\omega_{h^{r}}( \Delta)=\omega_{h^{r}}(n^{-1})\leq 2n^{-1},\]

_where \(\Delta\) is the maximum of the minimal distance of each data point with all other data points and \(\omega_{h^{r}}\) is the modulus of continuity, i.e., \(\omega_{h^{r}}(\Delta)=\sup_{|x-\tilde{x}|\leq|\Delta|}(h^{r}(x)-h^{r}(\tilde{ x}))\). Then for \(H^{r}_{\mathrm{nn}}(x):=(\tilde{h}^{r}(x_{1}),\cdots,\tilde{h}^{r}(x_{d}))\), we can derive_

\[\|H^{r}_{\mathrm{nn}}-H^{r}\|_{\mathcal{F}^{\infty}(K)}\leq\sqrt{d}\|\tilde{h} ^{r}-h^{r}\|_{\mathcal{F}^{\infty}([0,1])}\leq 2\sqrt{d}n^{-1}. \tag{2.1}\]

#### 2.1.2 INN for \(F\)

Now we explicitly construct an INN \(\tilde{F}_{\mathrm{nn}}\) to approximate the bi-Lipschitz map \(F\). The overall construction is quite lengthy and technical and is divided into several steps. Note that the construction is non-unique: there is an alternative construction by lifting the vector into \(\mathbb{R}^{2d+2}\); see the appendix for details. Throughout, we assume that we have access to a finite collection of evaluations \(\{(x^{\alpha},y^{\alpha})\}_{\alpha\in[n-1]^{d}}\) on a uniform grid. The construction proceeds in the following four steps:

1. make a small perturbation \(\eta\) to \(\{x^{\alpha}\}_{\alpha\in[n-1]^{d}}\) such that \(\{\eta(x^{\alpha})_{d}\}_{\alpha\in[n-1]^{d}}\) are distinct, i.e., the last component of all the mapped points \(\{\eta(x^{\alpha})_{d}\}_{\alpha\in[n-1]^{d}}\) are distinct;
2. construct a mapping \(\varphi^{N}\) that satisfies, for any \(\alpha\in[n-1]^{d}\), that \((\varphi^{N}\circ\eta(x^{\alpha}))_{j}=y^{\alpha}_{j}\) for any \(j=1,\cdots,d-1\) and \((\varphi^{N}\circ\eta(x^{\alpha}))_{d}=\eta(x^{\alpha})_{d}\), i.e., keeping the last entry of all the points unchanged;
3. choose an index \(j_{0}\in\{1,\cdots,d-1\}\), and then make a small perturbation \(\tilde{\eta}\), which ensures that \(\{[\tilde{\eta}\circ\varphi^{N}\circ\eta(x^{\alpha})]_{j_{0}}\}_{\alpha\in[n- 1]^{d}}\) are all distinct;
4. derive some \(\tilde{\varphi}^{N}\) so that \([\tilde{\varphi}^{N}\circ\tilde{\eta}\circ\varphi^{N}\circ\eta(x^{\alpha})]_{ d}=y^{\alpha}_{d}\) for any \(\alpha\in[n-1]^{d}\) and \([\tilde{\varphi}^{N}\circ\tilde{\eta}\circ\varphi^{N}\circ\eta(x^{\alpha})]_{ j}=[\tilde{\eta}\circ\varphi^{N}\circ\eta(x^{\alpha})]_{j}\) for any \(j=1,\cdots,d-1\).
The key is to ensure that each of these steps can be realized by a coupling-based INN, and Lipschitz constants of these INNs and their inverses are under precise control. The explicit constructions of these invertible networks follow that of coupling based INNs, and the details are given in the following lemmas. The construction represents the main technical contribution of the work.

Now we introduce the main technical tool, i.e., control family, in the construction. This concept was used in the dynamical system approach to universal approximation property of DNNs [35]. We employ a control family \(\mathcal{F}\) of three-layer fully connected NN architecture with ReLU \(\sigma\) defined by

\[\mathcal{F}:=\big{\{}V\sigma\big{(}W^{(2)}\sigma(W^{(1)}\cdot+b^{(1)})+b^{(2)} \big{)}\ |\ V,W^{(1)},W^{(2)}\in\mathbb{R}^{d\times d}.b^{(1)},b^{(2)}\in\mathbb{R}^{d }\big{\}}. \tag{2.2}\]

Accordingly, we define the attainable set \(\mathcal{A}_{\mathcal{F}}\) corresponding to \(\mathcal{F}\) by

\[\mathcal{A}_{\mathcal{F}}:=\bigcup_{T>0}\Big{\{}\psi_{\tau_{k}}^{f_{k}}\circ \psi_{\tau_{k-1}}^{f_{k-1}}\circ\cdots\circ\psi_{\tau_{1}}^{f_{1}}\ |\ k\geq 1,f_{j}\in \mathcal{F},\tau_{j}>0,\forall j=1,\cdots,k\ \text{and}\ \sum_{j=1}^{k}\tau_{j}=T\Big{\}}, \tag{2.3}\]

with \(\psi_{\tau}^{f}=y(\tau)\), where \(y(t)\) is the unique solution of the following initial value problem of the systems of first-order autonomous ordinary differential equations (ODEs):

\[\begin{cases}\frac{\mathrm{d}y}{\mathrm{d}t}=f(y),\quad t>0,\\ y(0)=x.\end{cases} \tag{2.4}\]

Since \(f(y)\) is Lipschitz continuous in \(y\), by the classical Picard-Lindelof theorem for ODEs (see, e.g., [32, Theorem 2.2.1]), problem (2.4) has a unique continuously differentiable solution \(y(t)\). Further, we have

\[\psi_{\tau}^{f}(x)=y(\tau;x)=x+\int_{0}^{\tau}f(y(t;x))\mathrm{d}t. \tag{2.5}\]

The map \(\psi_{\tau}^{-f}=z(\tau)\) given by the unique solution \(z(t)\) of the following initial value problem

\[\begin{cases}\frac{\mathrm{d}z}{\mathrm{d}t}=-f(z),\quad t>0,\\ z(0)=x,\end{cases} \tag{2.6}\]

is the inverse of the flow map \(\psi_{\tau}^{f}\), since \(\psi_{\tau}^{-f}\circ\psi_{\tau}^{f}(x)=\psi_{\tau}^{f}\circ\psi_{\tau}^{-f}( x)=x\). Thus, the functions in the attainable set \(\mathcal{A}_{\mathcal{F}}\) are invertible, which serve as the natural candidates for constructing various INNs. Moreover, if \(f(x)\in\mathcal{F}\), we have \(Df(Ax+b)\in\mathcal{F}\), for any \(D,A\in\mathbb{R}^{d\times d}\) and \(b\in\mathbb{R}^{d}\). Indeed, the constructions below use these elementary facts heavily.

It follows from the identity (2.5) and direct computation with the chain rule that the Jacobian \(J_{\psi_{\tau}^{f}}(x)\) of the flow map \(\psi_{\tau}^{f}\) (with respect to the input \(x\)) is given by

\[J_{\psi_{\tau}^{f}}(x)=I+\int_{0}^{\tau}J_{f}(y(t))J_{y(t)}(x)\mathrm{d}t=I+ \int_{0}^{\tau}J_{f}(\psi_{t}^{f}(x))J_{\psi_{t}^{f}}(x)\mathrm{d}t, \tag{2.7}\]

where \(J_{f}\) denotes the Jacobian of \(f\) with respect to its argument. Moreover, the solution \(J_{\psi_{\tau}^{f}}(x)\) to (2.7) is given by

\[J_{\psi_{\tau}^{f}}(x)=e^{\int_{0}^{\tau}J_{f}(\psi_{t}^{f}(x))\mathrm{d}t}. \tag{2.8}\]

Let \(\Delta_{j}(x^{\alpha})=\min_{\alpha_{1}\neq\alpha_{2}}|x_{j}^{\alpha_{1}}-x_ {j}^{\alpha_{2}}|\) be the minimal distance in the \(j\)th coordinate of the collection of points \(\{x^{\alpha}\}_{\alpha\in[n-1]^{d}}\). The construction of the INN \(\tilde{F}_{\mathrm{nn}}\) below requires \(\Delta_{j}(x^{\alpha})>0\) for some \(j\) at the beginning. The mapping \(\eta\in\mathcal{A}_{\mathcal{F}}\) constructed in the following lemma ensures \(\Delta_{d}\big{(}\eta(x^{\alpha})\big{)}=N^{-1}\). Thus, it fulfills step (i) of the overall procedure.

**Lemma 2.2**.: _For the uniform grids \(\{x^{\alpha}=\frac{\alpha}{n}\}_{\alpha\in[n-1]^{d}}\) of the hypercube \(K=[0,1]^{d}\), there exists an invertible mapping \(\eta\in\mathcal{A}_{\mathcal{F}}\) such that_

\[(\eta(x))_{j}=x_{j},\quad\forall j=1,\cdots,d-1,\ x\in\mathbb{R} ^{d},\quad\Delta_{d}\big{(}\eta(x^{\alpha})\big{)}=N^{-1},\] \[\max_{x\in\mathbb{R}^{d}}\|J_{\eta}(x)\|_{2}\leq\frac{n}{n-1}\quad \text{and}\quad\max_{x\in\mathbb{R}^{d}}\|J_{\eta^{-1}}(x)\|_{2}\leq\frac{n}{n- 1}.\]
Proof.: Let \(f_{j}(x)=\operatorname{diag}(\mathbf{0}_{d-1},1)\sigma(\operatorname{diag}( \mathbf{0}_{d-1},x_{j}))\in\mathcal{F}\), \(j=1,\cdots,d-1\), and define the mapping \(\eta:=\psi_{n^{-(d-1)}}^{f_{d-1}}\circ\cdots\circ\psi_{n^{-2}}^{f_{2}}\circ\psi _{n^{-1}}^{f_{1}}\in\mathcal{A}_{\mathcal{F}}\). Then there holds \(\{\eta(x^{\alpha})_{d}\}_{\alpha\in[n-1]^{d}}=\{iN^{-1}\}_{i=0}^{N-1}\). In fact, by the definitions of \(\psi_{\tau}^{f}\) and \(f_{j}\), we have

\[\eta(x^{\alpha})_{d}=x_{d}^{\alpha}+\sum_{j=1}^{d-1}n^{-j}x_{j}^{\alpha}. \tag{2.9}\]

Thus, \(\eta\) maps the set \(\{x^{\alpha}=\frac{\alpha}{n}\}_{\alpha\in[n-1]^{d}}\) to a set of \(N\) points with distinct values for the last component, and \(\Delta_{d}\big{(}\eta_{d}(x^{\alpha})\big{)}=N^{-1}.\) Meanwhile, for any \(j=1,\cdots,d-1\), by the identity (2.8), we have

\[J_{\psi_{n^{-j}}^{f_{j}}}(x)=\exp\Big{(}\int_{0}^{n^{-j}}J_{f_{j}}(\psi_{t}^{ f_{j}}(x))\mathrm{d}t\Big{)}:=\exp\bigl{(}n^{-j}A_{j}(x_{j})\big{)}=I+n^{-j}A_{j}(x_ {j}),\]

since the matrix \(A_{j}(x_{j})\in\mathbb{R}^{d\times d}\) has only one nonzero element \([A_{j}(x_{j})]_{d,j}=\frac{\mathrm{d}\sigma(x_{j})}{\mathrm{d}x_{j}}:=a_{j}\), which is also independent of \(x_{d}\), and thus \(A_{j}(x_{j})^{k}=0\) for any \(k\geq 2\). Furthermore, we have \(|a_{j}|\leq 1\). Also note that the components \(x_{j}\), \(j=1,\cdots,d-1\), stay unchanged under any of the mappings \(\psi_{n^{-j}}^{f_{j}}\). Thus, there holds

\[J_{\eta}(x)= J_{\psi_{n^{-(d-1)}}^{f_{d-1}}\circ\cdots\circ\psi_{n^{-1}}^{f_{1}} }(x)=J_{\psi_{n^{-d+1}}^{f_{d-1}}}(x)\cdots J_{\psi_{n^{-1}}^{f_{1}}}(x)\] \[= (I+n^{-(d-1)}A_{d-1}(x_{d-1}))\cdots(I+n^{-1}A_{1}(x_{1}))=I+\sum _{j=1}^{d-1}n^{-j}A_{j}(x_{j}),\]

in view of the following identity for any \(i,j=1,\cdots,d-1\): \(A_{i}(x_{i})A_{j}(x_{j})=0\). Thus, we can bound the spectral norm \(\|J_{\eta}(x)\|_{2}\) of the Jacobian \(J_{\eta}(x)\) by

\[\|J_{\eta}(x)\|_{2}\leq 1+\Big{\|}\sum_{j=1}^{d-1}n^{-j}A_{j}(x_{j})\Big{\|}_{2}\leq 1 +\left(\operatorname{tr}\Bigl{(}\bigl{(}\sum_{j=1}^{d-1}n^{-j}A_{j}(x_{j}) \bigr{)}^{t}\bigl{(}\sum_{j=1}^{d-1}n^{-j}A_{j}(x_{j})\bigr{)}\Bigr{)}\right)^ {\frac{1}{2}}\] \[\leq 1+\left(\sum_{j=1}^{d-1}n^{-2j}\right)^{\frac{1}{2}}\leq 1+ \left(\frac{n^{-2}}{1-n^{-2}}\right)^{\frac{1}{2}}\leq 1+\frac{1}{n-1}=\frac{n}{n-1}.\]

Since \(\det(J_{\eta}(x))=1\), the mapping \(\eta\) is invertible, and

\[\max_{x\in\mathbb{R}^{d}}\|J_{\eta^{-1}}(x)\|_{2}= \max_{x\in\mathbb{R}^{d}}\|(J_{\eta}(x))^{-1}\|_{2}=\Big{\|}I- \sum_{j=1}^{d-1}n^{-j}A_{j}(x_{j})\Big{\|}_{2}\leq\frac{n}{n-1}.\]

This shows the last estimate and completes the proof of the lemma. 

**Remark 2.3**.: _The forward process of the mapping \(\eta:x\to y\) in Lemma 2.2 is equivalent to a coupling-based INN \((\)with one weight layer and one added layer, i.e., identity mapping, and \(n\) neurons\()\) and the inverse process \(\eta^{-1}\) likewise [17, 2, 41]_

\[y_{j}=\begin{cases}x_{j},&j=1,\ldots,d-1,\\ x_{d}+\sum_{j=1}^{d-1}n^{-j}\sigma(x_{j}),&j=d,\end{cases}\quad\text{and}\quad x _{j}=\begin{cases}y_{j},&j=1,\cdots,d-1,\\ y_{d}-\sum_{j=1}^{d-1}n^{-j}\sigma(x_{j}),&j=d.\end{cases}\]

For the construction below, it is convenient to introduce the hypercube

\[S_{x^{\alpha}}(s):=(x_{1}^{\alpha}-s,x_{1}^{\alpha}+s)\times\cdots\times(x_{d}^ {\alpha}-s,x_{d}^{\alpha}+s),\]

which is a hypercube centered at \(x^{\alpha}\), each side of length \(2s\). Also we decompose a vector \(x\in\mathbb{R}^{d}\) into \(x=\begin{pmatrix}x^{\prime}\\ x_{d}\end{pmatrix}\), where the subvector \(x^{\prime}\in\mathbb{R}^{d-1}\) contains the first \(d-1\) components of the vector \(x\). Then we construct an NN, for any \(\alpha\in[n-1]^{d}\), to transport \((x^{\alpha})^{\prime}\) to \((y^{\alpha})^{\prime}(=\bigl{(}F(x^{\alpha})\bigr{)}^{\prime})\), and keep the last component \(x_{d}\) for any \(x\in\mathbb{R}^{d}\) and \((x)^{\prime}\) for any \(x\in(S_{x^{\alpha}}(\frac{1}{2N}))^{c}=\mathbb{R}^{d}\setminus S_{x^{\alpha}}( \frac{1}{2N})\) unchanged. Note that the set \(\{S_{x^{\alpha}}(\frac{1}{2N})\}_{\alpha\in[n-1]^{d}}\) consists of a collection of disjoint open hypercubes in \(\mathbb{R}^{d}\).
**Lemma 2.3**.: _Let \(\{(x^{\alpha})_{d}\}_{\alpha\in[n-1]^{d}}\) be distinct with the smallest distance \(\frac{1}{N}\). Then for any \(\alpha_{1}\in[n-1]^{d}\), there exists an invertible mapping \(\varphi\in\mathcal{A}_{\mathcal{F}}\), such that_

\[\varphi(x^{\alpha_{1}})^{\prime}=(y^{\alpha_{1}})^{\prime},\quad(\varphi(x))_{ d}= x_{d},\;\forall x\in\mathbb{R}^{d}\quad\text{and}\quad\varphi(x)=x,\; \forall x\in(S_{x^{\alpha_{1}}}(\tfrac{1}{2N}))^{c}.\]

_Furthermore, the following estimates hold_

\[\max_{x\in\mathbb{R}^{d}}\|J_{\varphi}(x)\|_{2}\leq 1+6N\|(y^{\alpha_{1}}-x^{ \alpha_{1}})^{\prime}\|_{2}\quad\text{and}\quad\max_{x\in\mathbb{R}^{d}}\|J_{ \varphi^{-1}}(x)\|_{2}\leq 1+6N\|(y^{\alpha_{1}}-x^{\alpha_{1}})^{\prime}\|_{2}.\]

Proof.: First, we define an intermediate mapping

\[\ell(x)=\sigma\big{(}-\sigma(\tfrac{1}{2}W^{0}x)+b^{0}\big{)}-\sigma\big{(}- \sigma(W^{0}x)+b^{0}\big{)}:=\ell^{(1)}(x)+\ell^{(2)}(x),\]

with the weight matrix \(W^{0}=\big{(}O_{d\times(d-1)}\quad\mathbf{1}_{d}\big{)}\in\mathbb{R}^{d\times d}\) and bias vector \(b^{0}=\mathbf{1}_{d}\in\mathbb{R}^{d}\). The mapping \(\ell(x)\) can be equivalently written as \(\ell(x)=\ell_{0}(x_{d})\mathbf{1}_{d}\), where the function \(\ell_{0}:\mathbb{R}\longrightarrow\mathbb{R}\) is defined by

\[\ell_{0}(x)=\sigma\big{(}-\sigma(\tfrac{1}{2}x)+1\big{)}-\sigma\big{(}-\sigma (x)+1\big{)}. \tag{2.10}\]

Note that \(\ell_{0}(x)\) is a hat function supported on \(x\in[0,2]\):

\[\ell_{0}(x)=\left\{\begin{array}{ll}0,&x\in(-\infty,0)\cup(2,+\infty),\\ \frac{x}{2},&x\in[0,1],\\ 1-\frac{x}{2},&x\in[1,2].\end{array}\right.\]

Next let \(A=\text{diag}(\mathbf{0}_{d-1},2\Delta_{d}^{-1}(x^{\alpha}))\in\mathbb{R}^{d \times d}\) and \(b^{\alpha_{1}}=(\mathbf{0}_{d-1},1-2\Delta_{d}^{-1}(x^{\alpha})x_{d}^{\alpha_ {1}})^{t}\in\mathbb{R}^{d}\). Note that the minimal distance \(\Delta_{d}(x^{\alpha})=\min_{\alpha_{2}\neq\alpha_{3}}|x_{d}^{\alpha_{2}}-x_{d }^{\alpha_{3}}|=\frac{1}{N}\) is independent of \(\alpha\), so is the matrix \(A\). Then for any \(\alpha_{1}\in[n-1]^{d}\), the following properties hold

\[(Ax^{\alpha_{1}}+b^{\alpha_{1}})_{d}=1\quad\text{and}\quad(Ax+b^{\alpha_{1}}) _{d}=2\Delta_{d}^{-1}(x^{\alpha})(x_{d}-x_{d}^{\alpha_{1}})+1\notin(0,2),\quad \forall x\in(S_{x^{\alpha_{1}}}(\tfrac{1}{2N}))^{c},\]

since \(|x_{d}-x_{d}^{\alpha_{1}}|\geq\frac{1}{2N}\) for any \(x\in(S_{x^{\alpha_{1}}}(\tfrac{1}{2N}))^{c}\) implies \(|2\Delta_{d}^{-1}(x^{\alpha})(x_{d}-x_{d}^{\alpha_{1}})|\geq 1\). Further, it implies that for any \(\alpha_{1}\in[n-1]^{d}\),

\[\ell(Ax^{\alpha_{1}}+b^{\alpha_{1}})=\tfrac{1}{2}\mathbf{1}_{d}\quad\text{and }\quad\ell(Ax+b^{\alpha_{1}})=0,\quad\forall x\in(S_{x^{\alpha_{1}}}(\tfrac{1 }{2N}))^{c}.\]

Next, we define, for some fixed \(\alpha_{1}\in[n-1]^{d}\), that

\[D_{\alpha_{1}}= \text{diag}(2(y_{1}^{\alpha_{1}}-x_{1}^{\alpha_{1}}),\cdots,2(y_{ d-1}^{\alpha_{1}}-x_{d-1}^{\alpha_{1}}),0)\in\mathbb{R}^{d\times d}.\]

Then we claim that \(\varphi:=\psi_{1}^{f}\) with \(f(x)=D_{\alpha_{1}}\ell(Ax+b^{\alpha_{1}})\) is the desired mapping if \(\varphi\in\mathcal{A}_{\mathcal{F}}\) and it is invertible. In fact, we can decompose \(\ell\) into \(\ell=\ell^{(1)}+\ell^{(2)}\), with \(\ell^{(1)},\ell^{(2)}\in\mathcal{F}\). Note that these two mappings actually depend only on the unchangeable value \(x_{d}\) along the time. Let \(f^{(i)}(x)=D_{\alpha_{1}}\ell^{(i)}(Ax+b^{\alpha_{1}})\) for \(i=1,2\), then \(\varphi=\psi_{1}^{f}=\psi_{1}^{f^{(1)}}\circ\psi_{1}^{f^{(2)}}\in\mathcal{A}_{ \mathcal{F}}\). By the definitions of \(f\) and \(\ell\), we have

\[f(x)= D_{\alpha_{1}}\ell(Ax+b^{\alpha_{1}})=2\Big{[}\sigma\Big{(}-\sigma \big{[}\tfrac{1}{2}\big{(}2\Delta_{d}^{-1}(x^{\alpha})(x_{d}-x_{d}^{\alpha_{1}}) +1\big{)}\big{]}+1\Big{)}\] \[-\sigma\Big{(}-\sigma\big{[}2\Delta_{d}^{-1}(x_{\alpha})(x_{d}-x_ {d}^{\alpha_{1}})+1\big{]}+1\Big{)}\Big{]}\begin{pmatrix}(y^{\alpha_{1}}-x^{ \alpha_{1}})^{\prime}\\ 0\end{pmatrix}:=h_{\alpha_{1}}(x_{d})\begin{pmatrix}(y^{\alpha_{1}}-x^{\alpha_ {1}})^{\prime}\\ 0\end{pmatrix}. \tag{2.11}\]

Then there holds

\[J_{\varphi}(x)= I+\frac{\mathrm{d}h_{\alpha_{1}}(x_{d})}{\mathrm{d}x_{d}}\begin{pmatrix} O_{(d-1)\times(d-1)}&(y^{\alpha_{1}}-x^{\alpha_{1}})^{\prime}\\ O_{1\times(d-1)}&0\end{pmatrix}:=I+H_{\alpha_{1}}(x).\]

It follows from direct computation that

\[\Big{|}\frac{\mathrm{d}h_{\alpha_{1}}(x_{d})}{\mathrm{d}x_{d}}\Big{|}\leq 2( \Delta_{d}^{-1}(x^{\alpha})+2\Delta_{d}^{-1}(x^{\alpha}))=6\Delta_{d}^{-1}(x^{ \alpha})=6N\quad\text{and}\quad\det(J_{\varphi}(x))=1.\]
This directly implies that \(\varphi\in\mathcal{A}_{\mathcal{F}}\) is invertible and further,

\[\|J_{\varphi}(x)\|_{2}\leq 1+\|H_{\alpha_{1}}(x)\|_{2}\leq 1+\sqrt{\operatorname{ tr}\bigl{(}(H_{\alpha_{1}}(x))^{t}H_{\alpha_{1}}(x)\bigr{)}}\leq 1+6N\|(y^{ \alpha_{1}}-x^{\alpha_{1}})^{\prime}\|_{2}.\]

Similarly, we deduce

\[\max_{x\in\mathbb{R}^{d}}\|J_{\varphi^{-1}}(x)\|_{2}=\max_{x\in\mathbb{R}^{d}} \|(J_{\varphi}(x))^{-1}\|_{2}=\max_{x\in\mathbb{R}^{d}}\|I-H_{\alpha_{1}}(x)\| _{2}\leq 1+6N\|(y^{\alpha_{1}}-x^{\alpha_{1}})^{\prime}\|_{2}.\]

This completes the proof of the lemma. 

**Remark 2.4**.: _The forward process of the mapping \(\varphi_{\alpha}:x\to y\) in Lemma 2.3 is equivalent to the following coupling-based NN \((\)with six weight layers and two added layers, i.e., identity mapping, and \(d\) neurons on each layer\()\) with the mapping \(h_{\alpha}:\mathbb{R}^{d}\to\mathbb{R}^{d}\) defined as (2.11) and likewise the inverse process \(\varphi^{-1}\):_

\[\begin{cases}y^{\prime}=x^{\prime}+h_{\alpha}(x_{d})(y^{\alpha}-x^{\alpha})^{ \prime},\\ y_{d}=x_{d},\end{cases}\quad\text{and}\quad\begin{cases}x_{d}=y_{d},\\ x^{\prime}=y^{\prime}-h_{\alpha}(x_{d})(y^{\alpha}-x^{\alpha})^{\prime}.\end{cases}\]

_Equivalently, \(y=x+D_{\alpha}\ell(Ax+b^{\alpha})\)._

So far we have already successfully transported the set \(\{x^{\alpha}\}_{\alpha\in[n-1]^{d}}\) to \(\{z^{\alpha}\}_{\alpha\in[n-1]^{d}}\) with \(z^{\alpha}=\bigl{(}(y^{\alpha})^{\prime t},\eta(x_{d}^{\alpha})\bigr{)}^{t}\). Following the ideas of Lemmas 2.2 and 2.3, we can further transport the collection \(\{z^{\alpha}\}_{\alpha\in[n-1]^{d}}\) of points to \(\{y^{\alpha}\}_{\alpha\in[n-1]^{d}}\). To this end, we first perturb the set \(\{z^{\alpha}\}_{\alpha\in[n-1]^{d}}\) by \(\tilde{\eta}\) so that \(\Delta_{j}(\tilde{\eta}(z^{\alpha}))\geq\epsilon\) for some \(j\neq d\) and \(\epsilon>0\). This distinctness at the \(j\)th coordinate is needed for performing the last step of the construction, and it can be achieved as follows. If there exists some \(j\neq d\) such that \(\Delta_{j}(z^{\alpha})>0\), we take \(\epsilon=\Delta_{j}(z^{\alpha})\) and \(\tilde{\eta}=\operatorname{id}\). Otherwise, we construct the desired mapping \(\tilde{\eta}\in\mathcal{A}_{\mathcal{F}}\) in Lemma 2.5 below. The next lemma analyzes the fundamental case of two points \(z^{1},z^{2}\in\mathbb{R}^{d}\).

**Lemma 2.4**.: _Let \(z^{1},z^{2}\in\mathbb{R}^{d}\), \(z^{1}\neq z^{2}\) and \(z_{1}^{1}=z_{1}^{2}\). Then for any \(a\in\mathbb{R}\), \(\epsilon:=a-z_{1}^{1}\) and \(0<\Delta\leq\max_{j=1,\cdots,d}|z_{j}^{1}-z_{j}^{2}|\), there exists an invertible mapping \(\tilde{\eta}_{0}\in\mathcal{A}_{\mathcal{F}}\) such that_

\[(\tilde{\eta}_{0}(z^{1}))_{1}=a,\quad(\tilde{\eta}_{0}(z^{1}))_{j }=z_{j}^{1},\quad j=2,\cdots,d,\] \[\tilde{\eta}_{0}(z)=z,\ \forall z\in(S_{z^{1}}(\Delta))^{c}\quad\text{and} \quad\|\tilde{\eta}_{0}(z)-z\|_{\mathcal{L}^{\infty}(\mathbb{R}^{d})}\leq| \epsilon|.\]

Proof.: The construction of \(\tilde{\eta}_{0}\) is inspired by the argument of [35, Lemma 4.14]. Since \(z^{1}\neq z^{2}\), there exists an index \(j\neq 1\) such that \(z_{j}^{1}\neq z_{j}^{2}\), and we may assume \(z_{2}^{1}\neq z_{2}^{2}\). Then for any \(0<\Delta\leq|z_{2}^{1}-z_{2}^{2}|\), consider the following initial value problem

\[\begin{cases}\begin{aligned} \frac{\operatorname{d}y}{ \operatorname{d}t}&=g(W_{22}y+e^{2}):=f(y),\quad t>0,\\ y(0)&=z,\end{aligned}\end{cases}\]

with \(W_{22}\in\mathbb{R}^{d\times d}\), \(e^{2}\in\mathbb{R}^{d}\) and \(g\) given by

\[W_{22}=\operatorname{diag}(0,\Delta^{-1},\mathbf{0}_{d-2})\in \mathbb{R}^{d\times d},\quad e^{2}=(0,1-\Delta^{-1}z_{2}^{1},\mathbf{0}_{d-2}) ^{t}\in\mathbb{R}^{d},\] \[g(z):=\sigma\big{(}-\sigma\big{(}\tfrac{1}{2}W_{12}z)+e^{1} \big{)}-\sigma\big{(}-\sigma(W_{12}z)+e^{1}\big{)},\]

where the matrix \(W_{12}\in\mathbb{R}^{d\times d}\) has only one nonzero entry \([W_{12}]_{12}=1\) and the vector \(e^{1}\in\mathbb{R}^{d}\) has one nonzero entry \([e^{1}]_{1}=1\). By the definitions of \(\Delta\) and \(g\), we have \(f(z^{1})=\frac{1}{2}e^{1}\) and \(f(z)=0\), for any \(z\in\mathbb{R}^{d}\) with \(|z_{2}-z_{2}^{1}|\geq\Delta\), since for \(z\in\mathbb{R}^{d}\) with \(|z_{2}-z_{2}^{1}|\geq\Delta\), \((W_{22}z+e^{2})_{2}=\Delta^{-1}(z_{2}-z_{2}^{1})+1\not\in(0,2)\), and \(g\) is supported on the interval \([0,2]\). Thus, for any \(z\in\mathbb{R}^{d}\) satisfying \(|z_{2}-z_{2}^{1}|\geq\Delta\), \(z_{j}^{1}\) for any \(j=2,\cdots,d\) stay unchanged. Moreover, since \(f(z^{1})=\frac{e^{1}}{2}\), for any \(\epsilon:=a-z_{1}^{1}\geq 0\), when \(t=2\epsilon\), we have

\[[\psi_{t}^{f}(z^{1})]_{1}-z_{1}^{1}=[\psi_{t}^{f}(z^{1})]_{1}-z_{1}^{2}=\epsilon,\]
\[\|\psi_{t}^{f}(z)-z\|_{\mathcal{L}^{\infty}(\mathbb{R}^{d})}\leq|[\psi_{t}^{f}(z^{ 1})]_{1}-z_{1}^{1}|=\epsilon,\]

since by assumption \(z_{1}^{1}=z_{1}^{2}\), and \(|f(z)|\leq\frac{1}{2}\). Similarly, for any \(\epsilon=a-z_{1}^{1}<0\), when \(t=-2\epsilon\), we have

\[[\psi_{t}^{-f}(z^{1})]_{1}-z_{1}^{1}=[\psi_{t}^{-f}(z^{1})]_{1}-z_ {1}^{2}=\epsilon,\] \[\|\psi_{t}^{-f}(z)-z\|_{\mathcal{L}^{\infty}(\mathbb{R}^{d})}\leq| [\psi_{t}^{-f}(z^{1})]_{1}-z_{1}^{1}|=-\epsilon.\]

Let

\[\tilde{\eta}_{0}(z):=\left\{\begin{array}{ll}\psi_{2\epsilon}^{f}(z),&\text { if }\epsilon\geq 0,\\ \psi_{-2\epsilon}^{-f}(z),&\text{ otherwise.}\end{array}\right.\]

Then repeating the argument of Lemma 2.3 shows \(\tilde{\eta}_{0}\in\mathcal{A}_{\mathcal{F}}\) by decomposing \(g=g^{(1)}+g^{(2)}\), where \(g^{(1)},g^{(2)}\in\mathcal{F}\) depend only on the unchangeable value. Note that

\[J_{\tilde{\eta}_{0}}(z)=\exp\bigl{(}2\epsilon\frac{\mathrm{d}h(z_{2})}{ \mathrm{d}z_{2}}W_{12}\bigr{)}:=\exp\bigl{(}2\epsilon A(z_{2})\bigr{)}=I+2 \epsilon A(z_{2}),\]

with \(A(z_{2})=\frac{\mathrm{d}h(z_{2})}{\mathrm{d}z_{2}}W_{12}\), where

\[h(z_{2})=\left\{\begin{array}{ll}\frac{1}{2}\bigl{(}\Delta^{-1}(z_{2}-z_{2} ^{1})+1\bigr{)},&\text{ if }\Delta^{-1}(z_{2}-z_{2}^{1})+1\in(0,1],\\ 1-\frac{1}{2}\bigl{(}\Delta^{-1}(z_{2}-z_{2}^{1})+1\bigr{)},&\text{ if }\Delta^{-1}(z_{2}-z_{2}^{1})+1\in(1,2],\\ 0,&\text{ otherwise.}\end{array}\right.\]

Thus, \(\det(J_{\tilde{\eta}_{0}}(z))=1\), and the mapping \(\tilde{\eta}\) is invertible and the inverse \((J_{\tilde{\eta}_{0}}(z))^{-1}\) of \(J_{\tilde{\eta}_{0}}(z)\) is given by \((J_{\tilde{\eta}_{0}}(z))^{-1}=I-2\epsilon A(z_{2})\). Then with the choice \(\tilde{\eta}_{0}\), the asserted properties in the lemma follow. 

**Remark 2.5**.: _The forward process of the mapping \(\eta:x\to y\) in Lemma 2.4 is equivalent to a coupling-based NN \((\)with four weight layers, two added layers, i.e., identity mapping and \(d\) neurons on each layer\()\), and likewise the inverse process \(\eta^{-1}\):_

\[y_{j}=\begin{cases}x_{j}+2\epsilon h(x_{2}),&j=1,\\ x_{j},&j=2,\cdots,d,\end{cases}\quad\text{and}\quad x_{j}=\begin{cases}y_{j},&j= 2,\cdots,d,\\ y_{j}+2\epsilon h(x_{2}),&j=1.\end{cases}\]

The next lemma generalizes the transformation in Lemma 2.4 to multiple points.

**Lemma 2.5**.: _Let the set of points \(\{z^{\alpha}=\bigl{(}(y^{\alpha})^{\prime t},x_{d}^{\alpha}\bigr{)}^{t}\}_{ \alpha\in[n-1]^{d}}\) be such that \(\Delta_{d}(z^{\alpha})=N^{-1}\). Then for any \(\epsilon\in(0,1)\), there exists an INN \(\tilde{\eta}\in\mathcal{A}_{\mathcal{F}}\) and an index \(j_{0}<d\) such that_

\[\Delta_{j_{0}}\bigl{(}\tilde{\eta}(z^{\alpha})\bigr{)}\geq N^{-1}\epsilon,\quad \max_{z\in\mathbb{R}^{d}}|(\tilde{\eta}(z)-z)_{j_{0}}|<\epsilon\quad\text{and} \quad(\tilde{\eta}(z^{\alpha}))_{j}=z_{j}^{\alpha},\;\forall j\neq j_{0}. \tag{2.12}\]

_Moreover, the following estimates hold_

\[\max_{x\in\mathbb{R}^{d}}\|J_{\tilde{\eta}}(x)\|_{2}\leq 1+2N\epsilon\quad \text{and}\quad\max_{x\in\mathbb{R}^{d}}\|J_{\tilde{\eta}^{-1}}(x)\|_{2}\leq 1+2N\epsilon.\]

Proof.: For any \(\epsilon\in(0,1)\) and any fixed index \(j_{0}<d\), there exists a sequence of disjoint closed intervals \(\{I_{j}\}_{j=1}^{J}\), each of length \(\epsilon\), with \(\mathbb{N}\ni J\leq N\), such that the inclusion \(\{z_{j_{0}}^{\alpha}\}_{\alpha\in[n-1]^{d}}\subset\bigcup_{j=1}^{J}I_{j}\) holds. Consider an interval \(I_{j}:=[b_{l},b_{r}]\) and assume that there are \(K\) components of \(\{z_{j_{0}}^{\alpha}\}_{\alpha\in[n-1]^{d}}\), denoted by \(\{z_{j_{0}}^{\alpha_{i}}\}_{i=1}^{K}\), lie in \(I_{j}\). Clearly we have \(K\leq N\). Then there exists a set \(\{a^{\alpha_{i}}\}_{i=1}^{K}\subset I_{j}\) such that

\[\min\big{(}\min_{1\leq i<j\leq K}|a^{\alpha_{i}}-a^{\alpha_{j}}|,\min_{1\leq i \leq K}|a^{\alpha_{i}}-b_{r}|\big{)}=K^{-1}\epsilon\geq N^{-1}\epsilon.\]

By Lemma 2.4, for any index \(1\leq i\leq K\), we can construct an invertible mapping \(\tilde{\eta}_{\alpha_{i}}\in\mathcal{A}_{\mathcal{F}}\) with \(\Delta=(2N)^{-1}\) such that

\[(\tilde{\eta}_{\alpha_{i}}(z^{\alpha_{i}}))_{j_{0}}=a^{\alpha_{i}},\;\;\|\tilde{ \eta}_{\alpha_{i}}(z)-z\|_{\mathcal{L}^{\infty}(\mathbb{R}^{d})}=|a^{\alpha_{i} }-z_{j_{0}}^{\alpha_{i}}|\leq\epsilon,\]
and \(\tilde{\eta}_{\alpha_{i}}(z)=z\) for any \(z\in(S_{z^{\alpha_{i}}}(\frac{1}{2N}))^{c}\). Similarly, there exist a set \(\{a^{\alpha_{i}}\}_{i=1}^{N}\subset\bigcup_{j=1}^{J}I_{j}\) and a sequence of invertible mapping \(\{\tilde{\eta}_{\alpha_{i}}\}_{i=1}^{N}\) such that

\[\min_{1\leq i<j\leq N}|a^{\alpha_{i}}-a^{\alpha_{j}}|\geq N^{-1}\epsilon,\ \ (\tilde{\eta}_{\alpha_{i}}(z^{\alpha_{i}}))_{j_{0}}=a^{\alpha_{i}},\ \ \|\tilde{\eta}_{\alpha_{i}}(z)-z\|_{\mathcal{L}^{\infty}( \mathbb{R}^{d})}=|a^{\alpha_{i}}-z_{j_{0}}^{\alpha_{i}}|\leq\epsilon,\]

and \(\tilde{\eta}_{\alpha_{i}}(z)=z\) for any \(z\in(S_{z^{\alpha_{i}}}(\frac{1}{2N}))^{c}\). Furthermore, we have \(J_{\tilde{\eta}_{\alpha_{i}}}(z)=I+2\epsilon A_{\alpha_{i}}(z_{d})\), where the matrix \(A_{\alpha_{i}}(z_{d})\in\mathbb{R}^{d\times d}\) has only one nonzero entry at \([A_{\alpha_{i}}(z_{d})]_{j_{0},d}\) with \(|[A_{\alpha_{i}}(z_{d})]_{j_{0},d}|\leq\frac{1}{2}\Delta^{-1}\). Next we define \(\tilde{\eta}:=\tilde{\eta}_{\alpha_{N}}\circ\cdots\circ\tilde{\eta}_{\alpha_ {1}}\in\mathcal{A}_{\mathcal{F}}\). It can be shown that \(\tilde{\eta}\) satisfies the properties in (2.12), since

\[\Delta_{j_{0}}\big{(}\tilde{\eta}(z^{\alpha})\big{)}=\min_{1\leq i<j\leq N}|a ^{\alpha_{i}}-a^{\alpha_{j}}|\geq N^{-1}\epsilon.\]

Further, since the supports of \(J_{\tilde{\eta}_{\alpha_{i}}}\) are disjoint from each other, there hold

\[\|J_{\tilde{\eta}}(x)\|_{2} \leq\max_{i=1,\cdots,N}\|J_{\tilde{\eta}_{\alpha_{i}}}(x)\|_{2} \leq\max_{i=1,\cdots,N}\|I+2\epsilon A_{\alpha_{i}}(z_{d})\|_{2}\leq 1+2 \epsilon\max_{i=1,\cdots,N}\|A_{\alpha_{i}}(z_{d})\|_{2}\] \[\leq 1+2\epsilon\max_{i=1,\cdots,N}\sqrt{\operatorname{tr}(A_{ \alpha_{i}}(z_{d}))^{t}(A_{\alpha_{i}}(z_{d}))}\leq 1+\epsilon\Delta^{-1}=1+2N\epsilon,\] \[\|(J_{\tilde{\eta}}(x))^{-1}\|_{2} \leq\max_{i=1,\cdots,N}\|I-2\epsilon A_{\alpha_{i}}(z_{d})\|_{2} \leq 1+2N\epsilon.\]

This completes the proof of the lemma. 

Under the condition \(\Delta_{j_{0}}\big{(}\tilde{\eta}(z^{\alpha})\big{)}\geq N^{-1}\epsilon\), we can now construct an INN \(\tilde{\varphi}\) for any \(\alpha_{1}\in[n-1]^{d}\), to transport \(z_{d}^{\alpha_{1}}\) to \(y_{d}^{\alpha_{1}}=(\tilde{\varphi}(z^{\alpha_{1}}))_{d}\), and keep \((z)^{\prime}\) for any \(z\in\mathbb{R}^{d}\) and \(z_{d}\) for any \(z\in(S_{z^{\alpha_{1}}}(\frac{\epsilon}{2N}))^{c}\) unchanged.

**Lemma 2.6**.: _Let \(\{z_{j_{0}}^{\alpha}\}_{\alpha\in[n-1]^{d}}\) be distinct with the smallest distance equal or larger than \(\Delta_{j_{0}}=\frac{\epsilon}{N}\) for some \(j_{0}\neq d\) and \(\epsilon>0\). Then for any \(\alpha_{1}\in[n-1]^{d}\), there exists an INN \(\tilde{\varphi}\in\mathcal{A}_{\mathcal{F}}\) such that_

\[(\tilde{\varphi}(z^{\alpha_{1}}))_{d}=y_{d}^{\alpha_{1}},\quad(\tilde{\varphi} (z))_{j}=z_{j},\;\forall z\in\mathbb{R}^{d},\;j\neq d\quad\text{and}\quad \tilde{\varphi}(z)=z,\;\forall z\in(S_{z^{\alpha_{1}}}(\frac{\epsilon}{2N}))^{c}.\]

_Furthermore, the following estimates hold_

\[\max_{z\in\mathbb{R}^{d}}\|J_{\tilde{\varphi}}(z)\|_{2}\leq 1+6N\epsilon^{-1}|y_{d}^ {\alpha_{1}}-z_{d}^{\alpha_{1}}|\quad\text{and}\quad\max_{x\in\mathbb{R}^{d}} \|J_{\tilde{\varphi}^{-1}}(x)\|_{2}\leq 1+6N\epsilon^{-1}|y_{d}^{\alpha_{1}}-z_{d}^{ \alpha_{1}}|.\]

Proof.: We define an intermediate map \(\ell(z)=\sigma\big{(}-\sigma\big{(}\frac{1}{2}Pz\big{)}+e^{1}\big{)}-\sigma \big{(}-\sigma\big{(}Pz\big{)}+e^{1}\big{)}\), where \(P\in\mathbb{R}^{d\times d}\) has only one nonzero entry \([P]_{d,j_{0}}=1\) and \(e^{1}=\mathbf{1}_{d}\in\mathbb{R}^{d}\). Let \(A\in\mathbb{R}^{d\times d}\) with only one nonzero entry \([A]_{j_{0},j_{0}}=2\Delta_{j_{0}}^{-1}\) and \(b^{\alpha_{1}}\in\mathbb{R}^{d}\) with only one nonzero entry \(b_{j_{0}}^{\alpha_{1}}=1-2\Delta_{j_{0}}^{-1}z_{j_{0}}^{\alpha_{1}}\). Then there hold

\[(Az^{\alpha_{1}}+b^{\alpha_{1}})_{j_{0}}=1,\quad\quad\forall\alpha _{1}\in[n-1]^{d},\] \[(Az+b^{\alpha_{1}})_{j_{0}}\notin(0,2),\quad\forall z\in(S_{z^{ \alpha_{1}}}(\frac{\epsilon}{2N}))^{c}.\]

These two properties together imply

\[[\ell(Az^{\alpha_{1}}+b^{\alpha_{1}})] =(\mathbf{0}_{d-1},\tfrac{1}{2})^{t},\quad\forall\alpha_{1}\in[n-1 ]^{d}, \tag{2.13}\] \[\ell(Az+b^{\alpha_{1}}) =0,\quad\quad\quad\quad\forall z\in(S_{z^{\alpha_{1}}}(\frac{ \epsilon}{2N}))^{c}. \tag{2.14}\]

For a fixed index \(\alpha_{1}\in[n-1]^{d}\), define \(D_{\alpha_{1}}=\operatorname{diag}(\mathbf{0}_{d-1},2(y_{d}^{\alpha_{1}}-z_{d}^{ \alpha_{1}}))\in\mathbb{R}^{d\times d}.\) Then \(\tilde{\varphi}:=\psi_{1}^{f}\in\mathcal{A}_{\mathcal{F}}\) with \(f(z)=D_{\alpha_{1}}\ell(Az+b_{\alpha_{1}})\) is the desired mapping. Indeed, by the definition of \(f\), we have

\[f(x)= 2\Big{[}\sigma\Big{(}-\sigma\big{[}\tfrac{1}{2}\big{(}2\Delta_{j_ {0}}^{-1}(z_{j_{0}}-z_{j_{0}}^{\alpha_{1}})+1\big{)}\big{]}+1\Big{)}\] \[\qquad-\sigma\Big{(}-\sigma\big{[}2\Delta_{j_{0}}^{-1}(z_{j_{0}}-z_{j _{0}}^{\alpha_{1}})+1\big{]}+1\Big{)}\Big{]}\begin{pmatrix}\mathbf{0}_{d-1} \\ y_{d}^{\alpha_{1}}-z_{d}^{\alpha_{1}}\end{pmatrix}:=h_{\alpha_{1}}(z_{j_{0}}) \begin{pmatrix}\mathbf{0}_{d-1}\\ y_{d}^{\alpha_{1}}-z_{d}^{\alpha_{1}}\end{pmatrix}.\]
Then there holds \(J_{\tilde{\varphi}}(z)=I+H_{\alpha_{1}}(z)\), where the matrix \(H_{\alpha_{1}}\in\mathbb{R}^{d\times d}\) has only one nonzero entry

\[[H_{\alpha_{1}}]_{d,j_{0}}=(y_{d}^{\alpha_{1}}-z_{d}^{\alpha_{1}})\frac{ \mathrm{d}h_{\alpha_{1}}(z_{j_{0}})}{\mathrm{d}z_{j_{0}}},\quad\text{with }\left|\frac{\mathrm{d}h_{\alpha_{1}}(z_{j_{0}})}{ \mathrm{d}z_{j_{0}}}\right|\leq 6\Delta_{j_{0}}^{-1}=\frac{6N}{\epsilon}.\]

Consequently,

\[\|J_{\tilde{\varphi}}(z)\|_{2} \leq 1+\|H_{\alpha_{1}}(z)\|_{2}\leq 1+6N\epsilon^{-1}|y_{d}^{ \alpha_{1}}-z_{d}^{\alpha_{1}}|,\] \[\max_{z\in\mathbb{R}^{d}}\|J_{\tilde{\varphi}^{-1}}(z)\|_{2} =\max_{z\in\mathbb{R}^{d}}\|(J_{\tilde{\varphi}}(z))^{-1}\|_{2}= \max_{z\in\mathbb{R}^{d}}\|I-H_{\alpha_{1}}(z)\|_{2}\leq 1+6N\epsilon^{-1}|y_{d}^ {\alpha_{1}}-z_{d}^{\alpha_{1}}|.\]

This completes the proof of the lemma. 

**Remark 2.6**.: _Similar to Remark 2.4, the forward process of the mapping \(\tilde{\varphi}_{\alpha}\) in Lemma 2.6 amounts to a coupling-based INN with six weight layers and two added layers, and \(d\) neurons on each layer._

Last we construct an INN \(\tilde{F}_{\mathrm{nn}}\) to approximate the map \(F\). This is the main result of this section.

**Theorem 2.1**.: _Given a finite collection of evaluations \(\{(x^{\alpha},y^{\alpha})\}_{\alpha\in[n-1]^{d}}\), for any tolerance \(\epsilon>0\), there exists a bi-Lipschitz continuous invertible mapping \(\tilde{F}_{\mathrm{nn}}\in\mathcal{A}_{\mathcal{F}}\), such that \(\|\tilde{F}_{\mathrm{nn}}(x^{\alpha})-y^{\alpha}\|_{2}<\epsilon\) for any \(\alpha\in[n-1]^{d}\), and the Lipschitz constants of \(\tilde{F}_{\mathrm{nn}}\) and \(\tilde{F}_{\mathrm{nn}}^{-1}\) are bounded by_

\[L_{\tilde{F}_{\mathrm{nn}}} \leq\frac{n}{n-1}\big{(}1+6Nc\big{)}(1+2N\epsilon)\big{(}1+6N \epsilon^{-1}c\big{)},\] \[L_{\tilde{F}_{\mathrm{nn}}^{-1}} \leq\frac{n}{n-1}\big{(}1+6Nc\big{)}(1+2N\epsilon)\big{(}1+6N \epsilon^{-1}c\big{)},\]

_with \(c:=\max_{i=1,\cdots N}\|y^{\alpha_{i}}-x^{\alpha_{i}}\|_{2}+n^{-1}\)._

Proof.: By Lemma 2.2, there exists an INN \(\eta\in\mathcal{A}_{\mathcal{F}}\) such that \(\Delta_{d}\big{(}\eta(x^{\alpha})\big{)}=N^{-1}\) for all \(\alpha\in[n-1]^{d}\) and

\[\max_{x\in\mathbb{R}^{d}}\|J_{\eta}(x)\|_{2}\leq\frac{n}{n-1}\quad\text{and} \quad\max_{x\in\mathbb{R}^{d}}\|J_{\eta^{-1}}(x)\|_{2}\leq\frac{n}{n-1}.\]

Then, by Lemma 2.3, we can construct a composite INN \(\varphi^{N}:=\varphi_{\alpha_{N}}\circ\cdots\circ\varphi_{\alpha_{1}}\), with \(\big{(}\varphi^{j}(\eta(x^{\alpha_{i}}))\big{)}^{\prime}=\big{(}y^{\alpha_{i} }\big{)}^{\prime}\) for any \(1\leq i\leq j\leq N\). Note that each \(\varphi_{\alpha_{i}}\), \(i=1,\cdots,N\), is the identity mapping in \((S_{\eta(x^{\alpha_{i}})}(\frac{1}{2N}))^{c}\), and \(\{S_{\eta(x^{\alpha_{i}})}(\frac{1}{2N})\}_{i=1}^{N}\) are disjoint from each other. Thus, we have

\[\max_{x\in\mathbb{R}^{d}}\|J_{\varphi^{N}}(x)\|_{2}\leq\max_{i=1,\cdots N,x\in \mathbb{R}^{d}}\|J_{\varphi_{\alpha_{i}}}(x)\|_{2}\leq 1+6N\max_{i=1,\cdots N}\|(y^{ \alpha_{i}}-x^{\alpha_{i}})^{\prime}\|_{2}.\]

Similarly, there holds

\[\max_{x\in\mathbb{R}^{d}}\|J_{(\varphi^{N})^{-1}}(x)\|_{2}\leq 1+6N\max_{i=1, \cdots N}\|(y^{\alpha_{i}}-x^{\alpha_{i}})^{\prime}\|_{2}.\]

Then by Lemma 2.5, we can construct \(\tilde{\eta}\in\mathcal{A}_{\mathcal{F}}\) such that

\[\Delta_{j_{0}}\big{(}\tilde{\eta}\circ\varphi^{N}\circ\eta(x^{\alpha})\big{)} \geq N^{-1}\epsilon,\quad\max_{x\in\mathbb{R}^{d}}\|J_{\tilde{\eta}}(x)\|_{2} \leq 1+2N\epsilon\quad\text{and}\quad\max_{x\in\mathbb{R}^{d}}\|J_{\tilde{ \eta}^{-1}}(x)\|_{2}\leq 1+2N\epsilon.\]

Then repeating the argument with Lemma 2.6, there exists an INN \(\tilde{\varphi}^{N}:=\tilde{\varphi}_{\alpha_{N}}\circ\cdots\circ\tilde{ \varphi}_{\alpha_{1}}\in\mathcal{A}_{\mathcal{F}}\) such that \((\tilde{\varphi}^{N}\circ\tilde{\eta}\circ\varphi^{N}\circ\eta(x^{\alpha_{i}})) )_{d}=(\tilde{\varphi}^{j}(z^{\alpha_{i}}))_{d}=y_{d}^{\alpha_{i}}\) for any \(1\leq i\leq j\leq N\) and

\[\max_{z\in\mathbb{R}^{d}}\|J_{\tilde{\varphi}}(z)\|_{2}\leq 1+6N\epsilon^{-1}\max_{i=1,\cdots N}|y_{d}^{\alpha_{i}}-z_{d}^{ \alpha_{i}}|,\] \[\max_{z\in\mathbb{R}^{d}}\|J_{\tilde{\varphi}^{-1}}(z)\|_{2}\leq 1+6N\epsilon^{-1}\max_{i=1,\cdots N}|y_{d}^{\alpha_{i}}-x_{d}^{ \alpha_{i}}|.\]

By (2.9), we can bound \(|y_{d}^{\alpha_{i}}-z_{d}^{\alpha_{i}}|\) by

\[|y_{d}^{\alpha_{i}}-z_{d}^{\alpha_{i}}|\leq|y_{d}^{\alpha_{i}}-x_{d}^{\alpha_{i}}| +|z_{d}^{\alpha_{i}}-x_{d}^{\alpha_{i}}|\leq|y_{d}^{\alpha_{i}}-x_{d}^{\alpha_{i}}| +\sum_{j=1}^{d-1}n^{-j}x_{j}^{\alpha_{i}}\leq|y_{d}^{\alpha_{i}}-x_{d}^{\alpha_{ i}}|+n^{-1}.\]
Finally, let \(\tilde{F}_{\mathrm{nn}}=\tilde{\varphi}^{N}\circ\tilde{\eta}\circ\varphi^{N} \circ\eta\in\mathcal{A}_{\mathcal{F}}\). This is the desired mapping. With \(c:=\max_{i=1,\cdots N}\|y^{\alpha_{i}}-x^{\alpha_{i}}\|_{2}+n^{-1}\), the preceding estimates together imply

\[L_{\tilde{F}_{\mathrm{nn}}}= \max_{z\in\mathbb{R}^{d}}\|J_{\tilde{F}_{\mathrm{nn}}}(z)\|\leq \frac{n}{n-1}\big{(}1+6Nc\big{)}(1+2N\epsilon)\big{(}1+6N\epsilon^{-1}c\big{)},\] \[L_{\tilde{F}_{\mathrm{nn}}^{-1}}= \max_{z\in\mathbb{R}^{d}}\|J_{\tilde{F}_{\mathrm{nn}}^{-1}}(z)\| \leq\frac{n}{n-1}\big{(}1+6Nc\big{)}(1+2N\epsilon)\big{(}1+6N\epsilon^{-1}c \big{)}.\]

This completes the proof of the theorem. 

**Remark 2.7**.: _By the construction in Theorem 2.1 and the analysis in Remarks 2.3-2.6, the INN \(\tilde{F}_{\mathrm{nn}}\) is a coupling-based INN with \(1+16N\) weight layers and \(1+6N\) added layers, and \(d\) neurons on each layer._

### Error estimation

Now we bound the error of the INN approximation \(F_{\mathrm{nn}}:=\tilde{F}_{\mathrm{nn}}\circ H_{\mathrm{nn}}^{r}\), where \(\tilde{F}_{\mathrm{nn}}\) (cf. Theorem 2.1) has \(1+16N\) weight layers, \(1+6N\) added layers and each layer having \(d\) neurons, and \(H_{\mathrm{nn}}^{r}(x):=(h_{\mathrm{nn}}^{r}(x_{1}),\cdots,h_{\mathrm{nn}}^{r} (x_{d}))\) (cf. Section 2.1.1) with a two-layer NN \(h_{\mathrm{nn}}^{r}\) of \(2(n+1)\) neurons.

**Theorem 2.2**.: _For any \(\epsilon>0\), define the bi-Lipschitz INNs \(F_{\mathrm{nn}}:=\tilde{F}_{\mathrm{nn}}\circ H_{\mathrm{nn}}^{r},F_{\mathrm{ nn}}^{-1}:=H_{\mathrm{nn}}^{r}\circ\tilde{F}_{\mathrm{nn}}^{-1}\in\mathcal{A}_{ \mathcal{F}}\), where the INNs \(\tilde{F}_{\mathrm{nn}}\), \(\tilde{F}_{\mathrm{nn}}^{-1}\) and \(H_{\mathrm{nn}}^{r}\) are defined in Theorem 2.1. Then the following error estimates hold_

\[\|F_{\mathrm{nn}}-F\|_{\mathcal{L}^{2}(K)}^{2}\leq 2\big{(}(L_{\tilde{F}_{\mathrm{nn}}}+L_{F})(1-r)\sqrt{d}n^{-1}+ \epsilon\big{)}^{2}+2L_{F}^{2}(2r-1)^{2}dn^{-2}\] \[+2\big{(}\tfrac{1}{2}(L_{\tilde{F}_{\mathrm{nn}}}+L_{F})\sqrt{d} n^{-1}+\epsilon\big{)}^{2}(1-r^{d}),\] \[\|F_{\mathrm{nn}}^{-1}-F^{-1}\|_{\mathcal{L}^{2}(F(K))}^{2}\leq 2L_{F}^{d}\Big{[}\Big{(}2+\big{(}L_{\tilde{F}_{\mathrm{nn}}^{-1}}+ L_{F^{-1}}\big{)}L_{F}(1-r)\Big{)}\sqrt{d}n^{-1}+L_{F^{-1}}\epsilon\Big{]}^{2}\] \[+2L_{F}^{d}(2r-1)^{2}dn^{-2}+8L_{F}^{d}(1-r^{d}).\]

Proof.: By the triangle inequality, we can decompose the total error \(\|F_{\mathrm{nn}}-F\|_{\mathcal{L}^{2}(K)}^{2}\) between the INN \(F_{\mathrm{nn}}\) and the bi-Lipschitz map \(F\) into two parts

\[\|F_{\mathrm{nn}}-F\|_{\mathcal{L}^{2}(K)}^{2}\leq 2\|\tilde{F}_{\mathrm{nn}}\circ H_{\mathrm{nn}}^{r}-F\circ H^{r} \|_{\mathcal{L}^{2}(K)}^{2}+2\|F\circ H^{r}-F\|_{\mathcal{L}^{2}(K)}^{2}\] \[= 2\|\tilde{F}_{\mathrm{nn}}\circ H^{r}-F\circ H^{r}\|_{\mathcal{L} ^{2}(K)}^{2}+2\|F\circ H^{r}-F\|_{\mathcal{L}^{2}(K)}^{2},\]

since \(H^{r}=H_{\mathrm{nn}}^{r}\), cf. the construction in Section 2.1.1. Then with the union of the (disjoint) hypercubes

\[K^{r}:=\bigcup_{\alpha\in[n-1]^{d}}[x_{1}^{\alpha},x_{1}^{\alpha}+\tfrac{r}{n} ]\times\cdots\times[x_{d}^{\alpha},x_{d}^{\alpha}+\tfrac{r}{n}],\]

and using the argument of Lemma 2.1, since \(|K\setminus K^{r}|=1-r^{d}\), we deduce

\[\|\tilde{F}_{\mathrm{nn}}\circ H^{r}-F\circ H^{r}\|_{\mathcal{L} ^{2}(K)}^{2}= \|\tilde{F}_{\mathrm{nn}}\circ H^{r}-F\circ H^{r}\|_{\mathcal{L} ^{2}(K^{r})}^{2}+\|\tilde{F}_{\mathrm{nn}}\circ H^{r}-F\circ H^{r}\|_{ \mathcal{L}^{2}(K\setminus K^{r})}^{2}\] \[\leq |K^{r}\|\|\tilde{F}_{\mathrm{nn}}\circ H^{r}-F\circ H^{r}\|_{ \mathcal{L}^{\infty}(K^{r})}^{2}+|K\setminus K^{r}|\|\tilde{F}_{\mathrm{nn}}- F\|_{\mathcal{L}^{\infty}(K)}^{2}\] \[\leq \|\tilde{F}_{\mathrm{nn}}\circ H^{r}-F\circ H^{r}\|_{\mathcal{L} ^{\infty}(K^{r})}^{2}+(1-r^{d})\|\tilde{F}_{\mathrm{nn}}-F\|_{\mathcal{L}^{ \infty}(K)}^{2}. \tag{2.15}\]

For any index \(\alpha\in[n-1]^{d}\), let \(\xi^{\alpha}:=\tilde{F}_{\mathrm{nn}}(x^{\alpha})-F(x^{\alpha})\), i.e., the approximation error at the point \(x^{\alpha}\). Then there holds \(\|\xi^{\alpha}\|_{2}\leq\epsilon\) by the construction of \(\tilde{F}_{\mathrm{nn}}\), cf. Lemmas 2.3 and 2.5. Note that

\[\|\tilde{F}_{\mathrm{nn}}\circ H^{r}-F\circ H^{r}\|_{\mathcal{L} ^{\infty}(K^{r})}\] \[= \|\tilde{F}_{\mathrm{nn}}\circ(H^{r}-\lim_{r\to 1^{-}}H^{r})+( \tilde{F}_{\mathrm{nn}}-F)\circ\lim_{r\to 1^{-}}H^{r}+F\circ(\lim_{r\to 1^{-}}H^{r}-H^{r})\|_{ \mathcal{L}^{\infty}(K^{r})}\] \[\leq \|\tilde{F}_{\mathrm{nn}}\circ(H^{r}-\lim_{r\to 1^{-}}H^{r})\|_{ \mathcal{L}^{\infty}(K^{r})}+\|(\tilde{F}_{\mathrm{nn}}-F)\circ\lim_{r\to 1^{-}}H^{r}\|_{ \mathcal{L}^{\infty}(K^{r})}\] \[+\|F\circ(\lim_{r\to 1^{-}}H^{r}-H^{r})\|_{\mathcal{L}^{\infty}(K^{r})}\]
\[\|H^{r}\circ F^{-1}-F^{-1}\|_{\mathcal{L}^{2}(F(K))}\leq L_{F}^{d}(2r-1)^{2}dn^{-2}.\]

Combining the preceding estimates completes the proof of the theorem. 

The above estimates are derived in the \(\mathcal{L}^{2}\) space, and similar results can be derived for the \(\mathcal{L}^{p}\) spaces, \(1\leq p<\infty\). The next result provides quantitative estimates for the INN approximations \(F_{\mathrm{nn}}\) and \(F_{\mathrm{nn}}^{-1}\).
**Corollary 2.1**.: _By choosing \(\epsilon=c_{\epsilon}n^{-1}\) with some \(c_{\epsilon}>0\) and_

\[r\geq c_{\mathrm{nn}}^{r}:=\max\big{(}\big{(}1-(L_{\tilde{F}_{\mathrm{nn}}}+L_{F })^{-2}\big{)}^{\frac{1}{d}},1-\big{(}L_{\tilde{F}_{\mathrm{nn}}^{-1}}+L_{F^{-1 }}\big{)}^{-1},(1-n^{-2})^{\frac{1}{d}}\big{)},\]

_the following error estimates hold_

\[\|F_{\mathrm{nn}}-F\|_{\mathcal{L}^{2}(K)}^{2} \leq 2\big{[}(3+L_{F}^{2})d+3c_{\epsilon}^{2}\big{]}n^{-2},\] \[\|F_{\mathrm{nn}}^{-1}-F^{-1}\|_{\mathcal{L}^{2}(F(K))}^{2} \leq 2L_{F}^{d}\big{[}\big{(}2(2+L_{F})^{2}+1\big{)}d+2L_{F^{-1}}c_{ \epsilon}+6\big{]}n^{-2}.\]

Proof.: By the choice of \(\epsilon\) and \(r\), we have \(L_{F}(1-r)\leq(L_{\tilde{F}_{\mathrm{nn}}}+L_{F})(1-r)\leq(L_{\tilde{F}_{ \mathrm{nn}}}+L_{F})(1-r^{d})^{\frac{1}{2}}\leq 1\), \(\big{(}L_{\tilde{F}_{\mathrm{nn}}^{-1}}+L_{F^{-1}}\big{)}(1-r)\leq 1\) and \(1-r^{d}\leq n^{-2}\). Then it follows from Theorem 2.1 that

\[\|F_{\mathrm{nn}}-F\|_{\mathcal{L}^{2}(K)}^{2}\leq 2\big{(}(L_{\tilde{F}_{\mathrm{nn}}}+L_{F})(1-r)\sqrt{d}+c_{ \epsilon}\big{)}^{2}n^{-2}+2L_{F}^{2}(2r-1)^{2}dn^{-2}\] \[+2\big{(}\tfrac{1}{2}(L_{\tilde{F}_{\mathrm{nn}}}+L_{F})\sqrt{d} +c_{\epsilon}\big{)}^{2}(1-r^{d})n^{-2}\] \[\leq 2\big{[}\big{(}\sqrt{d}+c_{\epsilon}\big{)}^{2}+L_{F}^{2}d+ \big{(}\tfrac{1}{2}\sqrt{d}+c_{\epsilon}n^{-1}\big{)}^{2}\big{]}n^{-2}\] \[\leq 2\big{[}2d+2c_{\epsilon}^{2}+L_{F}^{2}d+\tfrac{1}{2}d+2c_{ \epsilon}^{2}n^{-2}\big{]}n^{-2}\leq 2\big{[}(3+L_{F}^{2})d+3c_{\epsilon}^{2}\big{]}n^{-2}.\]

Similarly, we have

\[\|F_{\mathrm{nn}}^{-1}-F^{-1}\|_{\mathcal{L}^{2}(F(K))}^{2}\leq 2L_{F}^{d}\big{(}\big{(}2+(L_{\tilde{F}_{\mathrm{nn}}^{-1}}+L_{F ^{-1}})L_{F}(1-r)\big{)}\sqrt{d}+L_{F^{-1}}c_{\epsilon}\big{)}^{2}n^{-2}\] \[+2L_{F}^{d}(2r-1)^{2}dn^{-2}+8L_{F}^{d}(1-r^{d})\] \[\leq 2L_{F}^{d}\big{[}\big{(}(2+L_{F})\sqrt{d}+L_{F^{-1}}c_{\epsilon }\big{)}^{2}+d+4\big{]}n^{-2}\] \[\leq 2L_{F}^{d}\big{[}\big{(}2(2+L_{F})^{2}+1\big{)}d+2L_{F^{-1}c_{ \epsilon}}+6\big{]}n^{-2}.\]

This completes the proof of the corollary. 

**Remark 2.8**.: _The approximation error decays at a rate \(O(N^{-\frac{1}{d}})\), and thus it suffers from the usual curse of dimensionality. This result agrees with the fact that for Lipschitz maps, standard fully connected NNs also suffer from this curse similarly [15, Section 8.7] (up to a logarithmic factor). However, it is still unclear whether this is also the lower bound for coupling based INNs. In order to overcome the curse, additional structural information, e.g., manifold assumption or low complexity, is needed [15, Section 8.10]._

## 3 INN operator approximation

Now we consider the approximation of the following nonlinear operator equation using INNs:

\[F^{\dagger}(x)=y, \tag{3.1}\]

where \(F^{\dagger}:\mathcal{D}(F^{\dagger})\to\mathcal{I}(F^{\dagger})\) is a bi-Lipschitz continuous mapping with its domain \(\mathcal{D}(F^{\dagger})\subset\mathcal{X}\) and image space \(\mathcal{I}(F^{\dagger})\subset\mathcal{Y}\), and \(\mathcal{X}\) and \(\mathcal{Y}\) are infinite-dimensional separable Hilbert spaces with inner products \(\langle\cdot,\cdot\rangle\) and norms \(\|\cdot\|\), respectively. The model (3.1) can describe many important scientific problems. We aim at constructing a DNN that simultaneously approximates the forward process \(F^{\dagger}:x\to y\) and the inverse process \((F^{\dagger})^{-1}:y\to x\) by training the DNN on a finite collection of paired observations \(\{(x^{n},y^{n})\}_{n=1}^{N}\) (i.e., the training data). Throughout, we assume that the evaluation points \(x^{n}\)s are drawn independent and identically distributed (i.i.d.) with respect to an (unknown) probability measure \(\mu\) supported on \(\mathcal{X}\), and \(y^{n}\)s are i.i.d. with respect to the push-forward measure \(F^{\dagger}_{\sharp}\mu\).
### Principal component analysis

We first recall model reduction by principal component analysis (PCA). Following the construction in [10], we project any \(x\in\mathcal{X}\) and \(y\in\mathcal{Y}\) into finite-dimensional spaces by principal component analysis (PCA) on \(\mathcal{X}\) and \(\mathcal{Y}\), respectively. We approximate the identity mappings \(I_{\mathcal{X}}:\mathcal{X}\rightarrow\mathcal{X}\) and \(I_{\mathcal{Y}}:\mathcal{Y}\rightarrow\mathcal{Y}\) by the composition of two maps, known as the encoder and decoder, respectively, denoted by \(G_{\mathcal{U}}:\mathcal{U}\rightarrow\mathbb{R}^{d_{\mathcal{U}}}\) and \(G^{*}_{\mathcal{U}}:\mathbb{R}^{d_{\mathcal{U}}}\rightarrow\mathcal{U}\) where \(\mathcal{U}=\mathcal{X}\) or \(\mathcal{Y}\), which have a finite-dimensional range and domain, such that \(G^{*}_{\mathcal{U}}\circ G_{\mathcal{U}}\approx I_{\mathcal{U}}\). Specifically, let \(\nu\) denote a probability measure supported on the space \(\mathcal{U}\), and we make the assumption of a finite fourth moment: \(\mathbb{E}_{u\sim\nu}[\|u\|^{4}]<\infty\). We denote by \(\{u^{i}\}_{i=1}^{N}\) a finite collection of \(N\) samples drawn independent and identically distributed (i.i.d.) from the probability measure \(\nu\) that will be used as the training data on which PCA is based. We consider the empirical, non-centered covariance operator \(C^{\mathcal{U}}_{N}:=\frac{1}{N}\sum_{i=1}^{N}u^{i}\otimes u^{i}\), where \(\otimes\) denotes the outer product. \(C^{\mathcal{U}}_{N}\) is a non-negative, self-adjoint, trace-class operator on the Hilbert space \(\mathcal{U}\), of rank at most \(N\) with eigenvectors \(\{\phi^{N,i}_{\mathcal{U}}\}_{i=1}^{N}\) and its corresponding eigenvalues \(\lambda^{\mathcal{U}}_{N,1}\geq\lambda^{\mathcal{U}}_{N,2}\geq\cdots\geq \lambda^{\mathcal{U}}_{N,N}\geq 0\) in a nonincreasing order. We define that, with some truncation level \(0<d_{\mathcal{U}}\leq N\),

\[(\text{PCA encoder}) G_{\mathcal{U}}(u):= \big{(}\langle u,\phi^{N,1}_{\mathcal{U}}\rangle,\cdots,\langle u, \phi^{N,d_{\mathcal{U}}}_{\mathcal{U}}\rangle\big{)}^{t}\in\mathbb{R}^{d_{ \mathcal{U}}},\quad\forall u\in\mathcal{U},\] \[(\text{PCA decoder}) G^{*}_{\mathcal{U}}(u^{*}):= \sum_{i=1}^{d_{\mathcal{U}}}u^{*}_{i}\phi^{N,i}_{\mathcal{U}}, \quad\forall u^{*}=(u^{*}_{1},\cdots,u^{*}_{d_{\mathcal{U}}})^{t}\in\mathbb{R} ^{d_{\mathcal{U}}}.\]

Below we denote the associated projection operator via PCA by \(\mathcal{R}^{\mathcal{U}}_{N}\), which clearly depends on the random samples \(\{u^{i}\}_{i=1}^{N}\) used to compute the empirical covariance \(C^{\mathcal{U}}_{N}\).

With the PCA on the spaces \(\mathcal{X}\) and \(\mathcal{Y}\), we can now describe the reduced model. Specifically, we transform the evaluations \(\{(x^{\dagger,i},y^{\dagger,i})\}_{i=1}^{N}\) to \(\{x^{i}:=G_{\mathcal{X}}(x^{\dagger,i}),y_{i}:=G_{\mathcal{Y}}(y^{\dagger,i}) \}_{i=1}^{N}\), and finally return them back to \(\{\tilde{x}^{i}:=G^{*}_{\mathcal{X}}(x^{i}),\tilde{y}^{i}:=G^{*}_{\mathcal{Y}} (y^{i})\}_{i=1}^{N}\). Throughout, we take \(\{(\tilde{x}^{i},\tilde{y}^{i})\}_{i=1}^{N}\subset\{(x^{\dagger,i},y^{\dagger,i})\}_{i=1}^{N}\) which implies \(F^{\dagger}(\tilde{x})=\tilde{y}\), and define the following mapping

(3.2)

Then for any \(x^{i}\in G_{\mathcal{X}}(\mathcal{D}(F^{\dagger}))\), \(i=1,\ldots,N\), there holds \(\tilde{x}^{i}=G^{*}_{\mathcal{X}}(x^{i})\) and \(y^{i}=G_{\mathcal{Y}}(y^{\dagger,i})=G_{\mathcal{Y}}(\tilde{y}^{i})\). Since \(F^{\dagger}(\tilde{x}^{i})=\tilde{y}^{i}\), we have \(F(x^{i})=G_{\mathcal{Y}}\circ F^{\dagger}\circ G^{*}_{\mathcal{X}}(x^{i})=G_{ \mathcal{Y}}\circ F^{\dagger}(\tilde{x}^{i})=G_{\mathcal{Y}}(\tilde{y}^{i})=y^ {i}\). Further, we have the inverse mapping \(F^{-1}\) of \(F\) given by

\[F^{-1}= G_{\mathcal{X}}\circ(F^{\dagger})^{-1}\circ G^{*}_{\mathcal{Y}}:G_{ \mathcal{Y}}\big{(}\mathcal{I}(F^{\dagger})\big{)}\to G_{\mathcal{X}} \big{(}\mathcal{D}(F^{\dagger})\big{)},\]

since for any \(x\in G_{\mathcal{X}}\big{(}\mathcal{D}(F^{\dagger})\big{)}\), there holds \(F^{-1}\circ F(x)=G_{\mathcal{X}}\circ(F^{\dagger})^{-1}\circ G^{*}_{\mathcal{Y} }\circ G_{\mathcal{Y}}\circ F^{\dagger}\circ G^{*}_{\mathcal{X}}(x)=x\).

### Proposed approach for approximating bi-Lipschitz maps

Now we construct an approximation of the bi-Lipschitz map \(F^{\dagger}:\mathcal{X}\rightarrow\mathcal{Y}\) that approximates both \(F^{\dagger}\) and its inverse \((F^{\dagger})^{-1}\) simultaneously. It proceeds in two steps:

1. approximate the true mapping \(F^{\dagger}\) by the reduced one \(G^{*}_{\mathcal{Y}}\circ F\circ G_{\mathcal{X}}\) with the PCA encoder \(G_{\mathcal{X}}\) and decoder \(G^{*}_{\mathcal{Y}}\), with \(F\) given in (3.2), and also approximate the inverse mapping \((F^{\dagger})^{-1}\) by \(G^{*}_{\mathcal{X}}\circ F^{-1}\circ G_{\mathcal{Y}}\) with the PCA encoder \(G_{\mathcal{Y}}\) and the decoder \(G^{*}_{\mathcal{X}}\).
2. approximate the forward map \(G^{*}_{\mathcal{Y}}\circ F\circ G_{\mathcal{X}}\) and its inverse \(G^{*}_{\mathcal{X}}\circ F^{-1}\circ G_{\mathcal{Y}}\) (with \(F\) and \(F^{-1}\) defined on \(\mathbb{R}^{d}\)) by \(G^{*}_{\mathcal{Y}}\circ F_{\text{nn}}\circ G_{\mathcal{X}}\) and \(G^{*}_{\mathcal{X}}\circ F^{-1}_{\text{nn}}\circ G_{\mathcal{Y}}\) in Section 2, respectively, where \(F_{\text{nn}}\) is a coupling based INN.

The diagram in Fig. 1 illustrates the complete procedure of constructing and inferring the INN approximation. Given the training dataset \(\{x^{i},y^{i}\}_{i=1}^{N}\), we first generate the pairs of encoder and decoder on the space \(\mathcal{X}\) and \(\mathcal{Y}\), i.e., \((G_{\mathcal{X}},G^{*}_{\mathcal{X}})\) and \((G_{\mathcal{Y}},G^{*}_{\mathcal{Y}})\), and then learn the coupling-based NN \(F_{\text{nn}}\) by minimizing the following empirical loss

\[\mathcal{L}(\theta)=\sum_{i=1}^{N}\left(\frac{c_{0}}{2}\|G^{*}_{\mathcal{X}} \circ F^{-1}_{\text{nn}}\circ G_{\mathcal{Y}}(y^{i})-x^{i}\|_{\mathcal{X}}^{2}+ \frac{1}{2}\|G^{*}_{\mathcal{Y}}\circ F_{\text{nn}}\circ G_{\mathcal{X}}(x^{i})-y^ {i}\|_{\mathcal{Y}}^{2}\right),\]
where \(c_{0}>0\) is a scalar that controls the balance between the two terms in the loss, and \(\theta\) denotes the vector of trainable parameters in the INN \(F_{\mathrm{nn}}\) approximating the reduced map. The optimization problem is typically minimized by gradient type methods, e.g., Adam [30].

The obtained approximation involves several different sources of errors: the truncation / approximation error due to the use of PCA encoder / decoder, the approximation error due to the INN approximation to the bi-Lipschitz map, generalization error due to using a finite-number of samples (instead of the underlying distribution), and the optimization error due to the optimizer (which fails to find a global minimizer of the empirical loss \(\mathcal{L}(\theta)\), since the loss landscape is highly complex). Our analysis below focuses on the first two components. The key tool is to analyze the approximation error of the INN approximation in Section 2.

### Error analysis of INN operator approximation

Now we give the error analysis of the approximation for both forward and inverse processes, which involve both PCA model reduction and INN approximation of bi-Lipschitz map on finite-dimensional spaces. We prove that there exists an INN that achieves small approximation error in Theorem 3.2.

#### 3.3.1 Error estimate of PCA model reduction

Now we discuss the error due to the PCA. There are two sources of errors: one is due to the truncation, and the other due to random sampling (in approximating the covariance operator by an empirical one). We denote the covariance operator of the infinite-dimensional input data space \(\mathcal{X}\) and output data space \(\mathcal{Y}\) by \(C^{\mathcal{U}}:=\mathbb{E}_{u\sim\nu}[u\otimes u]\), with \(\mathcal{U}=\mathcal{X}\) or \(\mathcal{Y}\), respectively. By the i.i.d. assumption on the samples \(\{u^{i}\}_{i=1}^{N}\), we have \(C^{\mathcal{U}}=\mathbb{E}_{\{u^{i}\}\sim\nu}[C^{\mathcal{U}}_{N}]\), where \(\mathbb{E}[\cdot]\) denotes taking expectation with respect to the samples from the probability measure \(\nu\). We also define the eigenvectors of \(C^{\mathcal{U}}\) by \(\{\phi_{\mathcal{U}}^{i}\}_{i=1}^{\infty}\) and its corresponding nonnegative eigenvalues \(\lambda_{1}^{\mathcal{U}}\geq\lambda_{2}^{\mathcal{U}}\geq\cdots\geq\cdots\) in a nonincreasing order. Then the PCA encoder and decoder for the infinite-dimensional space \(\mathcal{U}\) are defined respectively by

\[\begin{split}\text{(PCA encoder)}&\quad G_{\mathcal{U}}(u ):=&\big{(}\langle u,\phi_{\mathcal{U}}^{1}\rangle,\cdots, \langle u,\phi_{\mathcal{U}}^{d_{\mathcal{U}}}\rangle\big{)}^{t}\in\mathbb{R} ^{d_{\mathcal{U}}},\quad\forall u\in\mathcal{U},\\ \text{(PCA decoder)}&\quad G_{\mathcal{U}}^{*}(u^{*}) :=&\sum_{i=1}^{d_{\mathcal{U}}}u_{i}^{*}\phi_{i}^{ \mathcal{U}},\quad\forall u^{*}=(u_{1}^{*},\cdots,u_{d_{\mathcal{U}}}^{*})^{ t}\in\mathbb{R}^{d_{\mathcal{U}}}.\end{split}\]

We define \(d:=d_{\mathcal{X}}=d_{\mathcal{Y}}\geq 2\) by choosing suitable truncation indices \(d_{\mathcal{X}}\) and \(d_{\mathcal{Y}}\). Let \(\tilde{u}=\mathcal{R}^{\mathcal{U}}_{N}u^{\dagger}=G^{*}_{\mathcal{U}}\circ G _{\mathcal{U}}u^{\dagger}\) (\(\mathcal{R}^{\mathcal{U}}_{N}\) is the orthogonal projection into the PCA subspace).

Figure 1: A schematic diagram for the proposed simultaneous approximations to \(F^{\dagger}\) and \((F^{\dagger})^{-1}\) based on the INN \(F_{\mathrm{nn}}\) and its inverse \(F_{\mathrm{nn}}^{-1}\).
Now, we bound the mean-squared errors \(\mathbb{E}_{\{x^{\dagger,i}\}\sim\nu}\big{[}\mathbb{E}_{x^{\dagger}\sim\nu}[\|x^{ \dagger}-\tilde{x}\|^{2}]\big{]}\) and \(\mathbb{E}_{\{y^{\dagger,i}\}\sim F_{\sharp}^{\dagger}\nu}\big{[}\mathbb{E}_{y ^{\dagger}\sim F_{\sharp}^{\dagger}\nu}[\|y^{\dagger}-\tilde{y}\|^{2}]\big{]}\) caused by computing the PCA encoder \(G_{\mathcal{U}}\) and decoder \(G_{\mathcal{U}}^{*}\) from the finite data \(\{(x^{\dagger,i},y^{\dagger,i})\}_{i=1}^{N}\) and truncating the order of the reduced model, on the spaces \(\mathcal{X}\) and \(\mathcal{Y}\), separately.

**Lemma 3.1**.: _Given a finite collection of evaluations \(\{u^{\dagger,i}\}_{i=1}^{N}\), which are i.i.d. with respect to a probability measure \(\nu\) supported on \(\mathcal{U}\), then with \(c_{\nu}=\Big{(}\mathbb{E}_{u\sim\nu}[\mathrm{tr}\big{(}(u\otimes u-C^{ \mathcal{U}})^{t}(u\otimes u-C^{\mathcal{U}})\big{)}]\Big{)}^{\frac{1}{2}}\), there holds_

\[\mathbb{E}_{\{u^{\dagger,i}\}\sim\nu}\big{[}\mathbb{E}_{u^{\dagger }\sim\nu}[\|u^{\dagger}-\tilde{u}\|^{2}]\big{]}\leq c_{\nu}\sqrt{d}N^{-\frac{ 1}{2}}+\sum_{j=d+1}^{\infty}\lambda_{j}^{\mathcal{U}}.\]

Proof.: The whole proof is shown in [10, Theorem 3.4] that

\[\mathbb{E}_{\{u^{\dagger,i}\}\sim\nu}\big{[}\mathbb{E}_{u^{\dagger }\sim\nu}[\|u^{\dagger}-\tilde{u}\|^{2}]\big{]}\leq c_{\nu}\sqrt{d}N^{-\frac{ 1}{2}}+\mathbb{E}_{u^{\dagger}\sim\nu}[\|u^{\dagger}-G_{\mathcal{U}}^{*}\circ G _{\mathcal{U}}(u^{\dagger})\|^{2}]=c_{\nu}\sqrt{d}N^{-\frac{1}{2}}+\sum_{j=d+1} ^{\infty}\lambda_{j}^{\mathcal{U}}.\]

This completes the proof of the lemma. 

Now we bound the error between the approximate mapping \(G_{\mathcal{Y}}^{*}\circ F\circ G_{\mathcal{X}}\) and the exact one \(F^{\dagger}\), and the analogue for the inverse \((\mathcal{F}^{\dagger})^{-1}\). In the statement, \(L_{F}\) denotes the Lipschitz constant of the map \(F\), and the constants \(c_{\mu}\) and \(c_{F^{\dagger}\mu}\) are defined in Lemma 3.1.

**Theorem 3.1**.: _For the operator \(F\) defined in (3.2), the following error estimates hold_

\[\mathbb{E}_{\{x^{\dagger,i}\}\sim\mu}\big{[}\mathbb{E}_{x^{\dagger }\sim\mu}[\|G_{\mathcal{Y}}^{*}\circ F\circ G_{\mathcal{X}}(x^{\dagger})-F^{ \dagger}(x^{\dagger})\|^{2}]\big{]}\] \[\leq 2(L_{F^{\dagger}}^{2}c_{\mu}+c_{F_{\sharp}^{\dagger}\mu})\sqrt{d }N^{-\frac{1}{2}}+2\Big{(}L_{F^{\dagger}}^{2}\sum_{j=d+1}^{\infty}\lambda_{j} ^{\mathcal{X}}+\sum_{j=d+1}^{\infty}\lambda_{j}^{\mathcal{Y}}\Big{)},\] \[\mathbb{E}_{\{y^{\dagger}_{i}\}\sim F_{\sharp}^{\dagger}\mu} \big{[}\mathbb{E}_{y^{\dagger}\sim F_{\sharp}^{\dagger}\mu}[\|G_{\mathcal{X}} ^{*}\circ F^{-1}\circ G_{\mathcal{Y}}(y^{\dagger})-(F^{\dagger})^{-1}(y^{ \dagger})\|^{2}]\big{]}\] \[\leq 2(L_{F^{\dagger-1}}^{2}c_{F_{\sharp}^{\dagger}\mu}+c_{\mu})\sqrt {d}N^{-\frac{1}{2}}+2\Big{(}L_{F^{\dagger-1}}^{2}\sum_{j=d+1}^{\infty}\lambda _{j}^{\mathcal{Y}}+\sum_{j=d+1}^{\infty}\lambda_{j}^{\mathcal{X}}\Big{)}.\]

Proof.: By the definition of the operator \(F\), the triangle inequality and Cauchy-Schwarz inequality, with \(\tilde{x}=G_{\mathcal{X}}^{*}\circ G_{\mathcal{X}}(x^{\dagger})\) and \(\tilde{y}=G_{\mathcal{Y}}^{*}\circ G_{\mathcal{Y}}\circ y^{\dagger}\), we have

\[\|G_{\mathcal{Y}}^{*}\circ F\circ G_{\mathcal{X}}(x^{\dagger})-F^{ \dagger}(x^{\dagger})\|^{2}=\|G_{\mathcal{Y}}^{*}\circ G_{\mathcal{Y}}\circ F ^{\dagger}\circ G_{\mathcal{X}}^{*}\circ G_{\mathcal{X}}(x^{\dagger})-F^{ \dagger}(x^{\dagger})\|^{2}\] \[\leq 2\|G_{\mathcal{Y}}^{*}\circ G_{\mathcal{Y}}\circ F^{\dagger}( \tilde{x})-G_{\mathcal{Y}}^{*}\circ G_{\mathcal{Y}}\circ F^{\dagger}(x^{ \dagger})\|^{2}+2\|G_{\mathcal{Y}}^{*}\circ G_{\mathcal{Y}}\circ y^{\dagger}-y ^{\dagger}\|^{2}\] \[\leq 2L_{G_{\mathcal{Y}}^{*}\circ G_{\mathcal{Y}}\circ F^{\dagger}}^{2} \|\tilde{x}-x^{\dagger}\|^{2}+2\|\tilde{y}-y^{\dagger}\|^{2}\leq 2L_{F^{ \dagger}}^{2}\|\tilde{x}-x^{\dagger}\|^{2}+2\|\tilde{y}-y^{\dagger}\|^{2},\]

since \(G_{\mathcal{Y}}^{*}\circ G_{\mathcal{Y}}\) is an orthogonal projection. Using the error estimate in Lemma 3.1, we derive

\[\mathbb{E}_{\{x^{\dagger,i}\}\sim\mu}\big{[}\mathbb{E}_{x^{\dagger }\sim\mu}[\|G_{\mathcal{Y}}^{*}\circ F\circ G_{\mathcal{X}}(x^{\dagger})-F^{ \dagger}(x^{\dagger})\|^{2}]\big{]}\] \[\leq 2L_{F^{\dagger}}^{2}\Big{(}c_{\mu}\sqrt{d}N^{-\frac{1}{2}}+\sum_ {j=d+1}^{\infty}\lambda_{j}^{\mathcal{X}}\Big{)}+2\Big{(}c_{F_{\sharp}^{\dagger} \mu}\sqrt{d}N^{-\frac{1}{2}}+\sum_{j=d+1}^{\infty}\lambda_{j}^{\mathcal{Y}}\Big{)}.\]

Similarly, there holds

\[\|G_{\mathcal{X}}^{*}\circ F^{-1}\circ G_{\mathcal{Y}}(y^{\dagger} )-F^{\dagger-1}(y^{\dagger})\|^{2}=\|G_{\mathcal{X}}^{*}\circ G_{\mathcal{X}} \circ F^{\dagger-1}\circ G_{\mathcal{Y}}^{*}\circ G_{\mathcal{Y}}(y^{\dagger})-(F ^{\dagger})^{-1}(y^{\dagger})\|^{2}\] \[\leq 2L_{G_{\mathcal{X}}^{*}\circ G_{\mathcal{X}}\circ F^{\dagger-1}}^{2} \|\tilde{y}-y^{\dagger}\|^{2}+2\|\tilde{x}-x^{\dagger}\|^{2}\leq 2L_{F^{\dagger-1}}^{2}\| \tilde{y}-y^{\dagger}\|^{2}+2\|\tilde{x}-x^{\dagger}\|^{2},\]

and hence

\[\mathbb{E}_{\{y^{\dagger,i}\}\sim F_{\sharp}^{\dagger}\mu}\big{[} \mathbb{E}_{y^{\dagger}\sim F_{\sharp}^{\dagger}\mu}[\|G_{\mathcal{X}}^{*} \circ F^{-1}\circ G_{\mathcal{Y}}(y^{\dagger})-F^{\dagger-1}(y^{\dagger})\|^{2}] \big{]}\]
\[\leq 2L^{2}_{F^{\dagger-1}}\Big{(}c_{F^{\dagger}_{\sharp}\mu}\sqrt{d}N^{- \frac{1}{2}}+\sum_{j=d+1}^{\infty}\lambda_{j}^{\mathcal{Y}}\Big{)}+2\Big{(}c_{ \mu}\sqrt{d}N^{-\frac{1}{2}}+\sum_{j=d+1}^{\infty}\lambda_{j}^{\mathcal{X}} \Big{)},\]

which directly completes the proof of the theorem. 

**Remark 3.1**.: _In the discussion, we have ignored the possible presence of data noise in \(y^{\dagger,i}\). Should there be any noise \(\{\xi_{i}\}_{i=1}^{N}\) in the observation, it will impact the accuracy of the reduced model. The mean squared errors indicate that the bound depends on the spectral decay of the covariance operators \(C^{\mathcal{X}}\) and \(C^{\mathcal{Y}}\) via the remainders \(\sum_{j=d+1}^{\infty}\lambda_{j}^{\mathcal{X}}\) and \(\sum_{j=d+1}^{\infty}\lambda_{j}^{\mathcal{Y}}\) and the number of samples \(N\) (at a rate \(N^{-\frac{1}{2}}\)). In practice, the covariance operator \(C^{\mathcal{X}}\) is often represented by an integral operator, and then its spectral decay can be characterized by the smoothness of the associated kernel [20]._

#### 3.3.2 Error estimate of INN operator approximation

Now we derive the full error estimate of the INN approximation of the forward and inverse mappings. The following theorem is the main result of this section.

**Theorem 3.2**.: _Let \(F_{\mathrm{nn}}:=\tilde{F}_{\mathrm{nn}}\circ H^{r}_{\mathrm{nn}}\) and \(F_{\mathrm{nn}}^{-1}=H^{r}_{\mathrm{nn}}\circ\tilde{F}_{\mathrm{nn}}^{-1}\) when \(\epsilon=c_{\epsilon}n^{-1}\) with some \(c_{\epsilon}>0\) and \(r\geq c^{r}_{\mathrm{nn}}\), where \(H^{r}_{\mathrm{nn}}(x):=(h^{r}_{\mathrm{nn}}(x_{1}),\cdots,h^{r}_{\mathrm{nn}} (x_{d}))\) with a two-layer NN \(h^{r}_{\mathrm{nn}}\) of \(2(N^{\frac{1}{d}}+1)\) neurons, and \(\tilde{F}_{\mathrm{nn}}\) has \(1+16N\) weight layers, \(1+6N\) added layers and \(N^{\frac{1}{2}}\) neurons. Then the expected errors of the whole systems are bounded by_

\[\mathbb{E}_{\{x^{\dagger}_{i}\}\sim\mu}\big{[}\mathbb{E}_{x^{ \dagger}\sim\mu}[\|G^{*}_{\mathcal{Y}}\circ F_{\mathrm{nn}}\circ G_{\mathcal{ X}}(x^{\dagger})-F^{\dagger}(x^{\dagger})\|^{2}]\big{]}\] \[\leq 2c_{\mathrm{nn}}N^{-\frac{2}{2}}+4(L^{2}_{F^{\dagger}}c_{\mu}+c_ {F^{\dagger}_{\sharp}\mu})\sqrt{d}N^{-\frac{1}{2}}+4(L^{2}_{F^{\dagger}}\sum_ {j=d+1}^{\infty}\lambda_{j}^{\mathcal{X}}+\sum_{j=d+1}^{\infty}\lambda_{j}^{ \mathcal{Y}}),\] \[\mathbb{E}_{\{y^{\dagger}_{i}\}\sim F^{\dagger}_{\sharp}\mu} \big{[}\mathbb{E}_{y^{\dagger}\sim F^{\dagger}_{\sharp}\mu}[\|G^{*}_{ \mathcal{X}}\circ F_{\mathrm{nn}}^{-1}\circ G_{\mathcal{Y}}(y^{\dagger})-(F^{ \dagger})^{-1}(y^{\dagger})\|^{2}]\big{]}\] \[\leq 2c^{r}_{\mathrm{nn}}N^{-\frac{2}{2}}+4(L^{2}_{F^{\dagger-1}}c_{F ^{\dagger}_{\sharp}\mu}+c_{\mu})\sqrt{d}N^{-\frac{1}{2}}+4\Big{(}L^{2}_{F^{ \dagger-1}}\sum_{j=d+1}^{\infty}\lambda_{j}^{\mathcal{Y}}+\sum_{j=d+1}^{\infty }\lambda_{j}^{\mathcal{X}}\Big{)},\]

_with \(c_{\mathrm{nn}}=2\big{[}(3+L^{2}_{F})d+3c_{\epsilon}^{2}\big{]}\) and \(c^{r}_{\mathrm{nn}}=2L^{d}_{F}[\big{(}2(2+L_{F})^{2}+1\big{)}d+2L_{F^{-1}}c_{ \epsilon}+6]\)._

Proof.: By the triangle inequality, we have

\[\|G^{*}_{\mathcal{Y}}\circ F_{\mathrm{nn}}\circ G_{\mathcal{X}}(x ^{\dagger})-F^{\dagger}(x^{\dagger})\|^{2}\leq 2\|G^{*}_{\mathcal{Y}}\circ F_{\mathrm{nn}}\circ G_{ \mathcal{X}}(x^{\dagger})-G^{*}_{\mathcal{Y}}\circ F\circ G_{\mathcal{X}}(x^{ \dagger})\|^{2}\] \[+2\|G^{*}_{\mathcal{Y}}\circ F\circ G_{\mathcal{X}}(x^{\dagger})-F ^{\dagger}(x^{\dagger})\|^{2}:=\mathrm{I}+\mathrm{II}.\]

It suffices to bound the two terms separately. For the term \(\mathrm{I}\), we have

\[\mathbb{E}_{\{x^{\dagger}_{i}\}\sim\mu}\big{[}\mathbb{E}_{x^{ \dagger}\sim\mu}[\mathrm{I}]\big{]}\leq 2\mathbb{E}_{\{x_{i}\}\sim G_{X^{\sharp}\mu}}\big{[}\mathbb{E}_{x \sim G_{X^{\sharp}\mu}}[\|G^{*}_{\mathcal{Y}}\circ F_{\mathrm{nn}}(x)-G^{*}_{ \mathcal{Y}}\circ F(x)\|^{2}]\big{]}\] \[\leq 2\mathbb{E}_{\{x_{i}\}\sim G_{X^{\sharp}\mu}}\big{[}\|G^{*}_{ \mathcal{Y}}\|^{2}\|F_{\mathrm{nn}}-F\|^{2}_{\mathcal{L}^{2}(K)}\big{]}\leq 2\|F_{ \mathrm{nn}}-F\|^{2}_{\mathcal{L}^{2}(K)}.\]

In view of Remark 2.1, there holds

Meanwhile, by Theorem 3.1, we have

\[\mathbb{E}_{\{x^{\dagger}_{i}\}\sim\mu}\big{[}\mathbb{E}_{x^{ \dagger}\sim\mu}[\Pi]\big{]}\leq 4(L^{2}_{F^{\dagger}}c_{\mu}+c_{F^{\dagger}_{\sharp}\mu})\sqrt{d}N^{-\frac{1}{ 2}}+4\Big{(}L^{2}_{F^{\dagger}}\sum_{j=d+1}^{\infty}\lambda_{j}^{\mathcal{X}}+ \sum_{j=d+1}^{\infty}\lambda_{j}^{\mathcal{Y}}\Big{)}.\]

Combining the last two estimates yields the first assertion. Similarly, there holds

\[\|G^{*}_{\mathcal{X}}\circ F_{\mathrm{nn}}^{-1}\circ G_{\mathcal{Y}}(y^{ \dagger})-(F^{\dagger})^{-1}(y^{\dagger})\|^{2}\leq 2\|G^{*}_{\mathcal{X}}\circ F_{\mathrm{nn}}^{-1}\circ G_{ \mathcal{Y}}(y^{\dagger})-F^{-1}(y^{\dagger})\|^{2}\]
\[+2\|G_{\mathcal{X}}^{*}\circ F^{-1}\circ G_{\mathcal{Y}}(y^{\dagger})-(F^{ \dagger})^{-1}(y^{\dagger})\|^{2}=:\mathrm{III}+\mathrm{IV}.\]

By Remark 2.1 and Theorem 3.1, we derive

\[\mathbb{E}_{\{y^{\dagger}_{i}\}\sim F^{\dagger}_{1}\mu}\big{[} \mathbb{E}_{y^{\dagger}\sim F^{\dagger}_{1}\mu}[\mathrm{III}]\big{]} \leq 2\|F_{\mathrm{nn}}^{-1}-F^{-1}\|_{\mathcal{L}^{2}(F(K))}^{2} \leq 2c^{\prime}_{\mathrm{nn}}N^{-\frac{2}{3}},\] \[\mathbb{E}_{\{y^{\dagger}_{i}\}\sim F^{\dagger}_{1}\mu}\big{[} \mathbb{E}_{y^{\dagger}\sim F^{\dagger}_{1}\mu}[\mathrm{IV}]\big{]} \leq 4\big{(}L^{2}_{F^{\dagger-1}}c_{F^{\dagger}_{1}\mu}+c_{\mu} \big{)}\sqrt{d}N^{-\frac{1}{2}}+4\Big{(}L^{2}_{F^{\dagger-1}}\sum_{j=d+1}^{ \infty}\lambda_{j}^{\mathcal{Y}}+\sum_{j=d+1}^{\infty}\lambda_{j}^{\mathcal{X }}\Big{)}.\]

Combining the preceding estimates completes the proof of the theorem. 

**Remark 3.2**.: _Theorem 3.2 gives only the existence of an INN that approximates the forward and inverse maps simultaneously to a certain tolerance. It does not address the important issue of realizing the approximation in practice via optimizing the empirical loss \(\mathcal{L}(\theta)\), for which one has to analyze also the statistical error and optimization error. The former would indicate how many samples are needed in order to achieve the tolerance, and the latter is notoriously challenging due to the complex landscape of the associated optimization problem. We leave the study of these important errors to future works._

## 4 Numerical experiments

In this section, we present preliminary numerical results to showcase the feasibility of the proposed INN based approach to simultaneously approximate the forward and inverse processes of a bi-Lipschitz map on infinite-dimensional spaces. Consider the following second-order elliptic PDE on a smooth bounded domain \(D\) (with a boundary \(\partial D\)):

\[\begin{cases}-\nabla\cdot(u(x)\nabla y(x))=f,&\text{in }D,\\ y=0,&\text{on }\partial D,\end{cases} \tag{4.1}\]

where \(u\in L^{\infty}(\Omega)\) is the diffusion coefficient in the PDE, and \(f\in L^{2}(D)\) is the given source term. The diffusion coefficient \(u\) is assumed to satisfy the standard ellipticity and boundedness assumption so that for any given \(u\), there exists a unique solution \(y\equiv y(u)\in H^{1}_{0}(\Omega)\). The map \(F\) of interest is defined by \(F:u\mapsto y\). Note that this setting has been extensively studied in the literature [14], and it often serves as a model problem for parameter identification of partial differential equations.

In the numerical experiments below, we take a unit square \(D=(0,1)^{2}\), and \(\mathrm{fix}\;f\equiv 1\). To represent the diffusion coefficient \(u\), we employ truncated Karhunen-Loeve expansion [19, 5]. We simulate \(M=10000\) samples of the coefficient \(u\), generated by the following Fourier representation

\[u(x,\xi)=2+\sum_{i,j=1}^{20}\frac{\xi_{i,j}}{i^{3}+j^{3}}\cos(i\pi x_{1})\cos (j\pi x_{2}),\quad x=[x_{1},x_{2}]\in D,\]

where \(\xi=[\xi_{1,1},\xi_{1,2},\cdots,\xi_{1,20},\xi_{2,1},\cdots,\xi_{20,20}]\in \mathbb{R}^{20\times 20}\) and each \(\xi_{i,j}\) follows the standard normal distribution \(N(0,1)\) (for each entry independently). For each realization \(u\), we discretize the boundary value problem (4.1) with a mesh size \(h=\frac{1}{50}\) using the standard piecewise linear Galerkin finite element method, implemented in FEniCS [37], to obtain the corresponding output \(y_{h}\in\mathbb{R}^{51\times 51}\). For the ease of exposition, we apply a one-to-one transformation \(T\) on \(u\), given by \(T(u)=\xi\in\mathbb{R}^{20\times 20}\), to transform the infinite-dimensional inputs to a finite-dimensional Euclidean space without additional discretization. Then we reduce the dimensionalities of both inputs \(\xi\) and outputs \(y_{h}\) by PCA (using the implementation in scikit-learn, available from [https://github.com/scikit-learn/scikit-learn/releases/tag/1.2.2](https://github.com/scikit-learn/scikit-learn/releases/tag/1.2.2)) with a reduced dimensionality \(d=10\) (which capture more than \(99\%\) of the spectral energy) and obtain the truncated inputs \(\hat{u}=[\hat{u}_{1},\cdots,\hat{u}_{10}]\in\mathbb{R}^{10}\) and truncated outputs \(\hat{y}=[\hat{y}_{1},\cdots,\hat{y}_{10}]\in\mathbb{R}^{10}\). We also record the \(10\) largest singular values of the empirical covariance to form normalized, in terms of the \(\ell^{1}\)-norm, weight vectors \(w_{u}\in\mathbb{R}^{10}\) and \(w_{y}\in\mathbb{R}^{10}\). This normalization amounts to working in suitable Sobolev space scales so that the map is nearly bi-Lipschitz, in accordance with the known stability estimates [26].
In our experiments, we use an identical NN architecture with 3 blocks. One block of the forward process and inverse process is sketched in Fig. 2, which are commonly used in an affine coupling layer [17, 2]. In the figure, we define the maps \(g_{o}\) (odd part), \(g_{e}\) (even part) and \(g_{c}\) (combination) as

\[g_{o}(u) =[u_{1},u_{3},\cdots,u_{9}],\quad\forall u\in\mathbb{R}^{10},\] \[g_{e}(u) =[u_{2},u_{4},\cdots,u_{10}],\quad\forall u\in\mathbb{R}^{10},\] \[g_{c}(u^{\prime},u^{\prime\prime}) =[u^{\prime}_{1},u^{\prime\prime}_{1},u^{\prime}_{2},u^{\prime \prime}_{2},\cdots,u^{\prime}_{5},u^{\prime\prime}_{5}],\quad\forall u^{\prime },u^{\prime\prime}\in\mathbb{R}^{5}.\]

Let \(\mathcal{T}_{r}\) be a training dataset of size \(N_{t}=|\mathcal{T}_{r}|\) (taken to be 100, 500, and 1000), and \(\mathcal{T}_{e}\) be the test data set of size \(|\mathcal{T}_{e}|=10000\). We train the INN for the forward and inverse processes as constructed above simultaneously, denoted by \(\Phi_{\text{nn}}\) and \(\Phi_{\text{nn}}^{-1}\), respectively, on the data in the training set \(\mathcal{T}_{r}\). We define \(\hat{y}_{\text{nn}}=\Phi_{\text{nn}}(\hat{u})\) and \(\hat{u}_{\text{nn}}=\Phi_{\text{nn}}^{-1}(\hat{y})\). The loss function \(\ell\) that we employ for the training is defined by

\[\ell(\hat{u},\hat{y})=\frac{c_{0}}{2}\sum_{(\hat{u},\hat{y})\in \mathcal{T}_{r}}\|(\hat{u}-\hat{u}_{\text{nn}})\odot w_{u}\|^{2}+\frac{1}{2} \sum_{(\hat{u},\hat{y})\in\mathcal{T}_{r}}\|(\hat{y}-\hat{y}_{\text{nn}}) \odot w_{y}\|^{2}, \tag{4.2}\]

where \(c_{0}\) is a penalty parameter to be specified according to the relative errors of forward and inverse processes, and \(\odot\) denotes the Hadamard product between two vectors. In the experiment, we fix \(c_{0}\) at \(0.001\) for problem (4.1), since the magnitude of the loss of the inverse process \(\sum_{(\hat{u},\hat{y})\in\mathcal{T}_{r}}\|(\hat{u}-\hat{u}_{\text{nn}}) \odot w_{u}\|^{2}\) is more than 100 times that of the forward process \(\sum_{(\hat{u},\hat{y})\in\mathcal{T}_{r}}\|(\hat{y}-\hat{y}_{\text{nn}}) \odot w_{y}\|^{2}\), and the Lipschitz constant of the inverse process is much bigger than that of the forward process. The optimization problem (4.2) is minimized by the Adam algorithm [30] with a learning rate \(r_{0}=0.001\).

First, we compare our results to that obtained by two different 5-layer fully-connected NN (FNN) for forward and inverse processes of problem (4.1), respectively. Likewise we run Adam [30] for at most

Figure 2: The forward affine coupling block with input \(\hat{u}\) and output \(\hat{u}^{\prime}\) (top) and the inverse affine coupling block with input \(\hat{y}\) and output \(\hat{y}^{\prime}\) (bottom). Both \(\phi_{W}\) and \(\phi_{b}\) are four-layer full connected NNs with ReLU.
500000 steps to minimize the resulting optimization problem. The numerical results are given in Table 1, where \(e_{g}\) and \(e_{a}\) denote respectively the smallest relative generalization error (on the test dataset \(\mathcal{T}_{e}\)) and the relative approximation error (of \(N_{t}\) training data points) of the NN that trained on the training dataset \(\mathcal{T}_{r}\) along the iteration trajectory. The numerical results show that the proposed INN approach can achieve both approximation and generalization errors comparable with that by the FNN for all sizes of training dataset for the forward process; and can also achieve the generalization accuracy largely comparable with that by the FNN for the inverse process, albeit FNN tends to do a slightly better job. Note that a carefully tuned \(c_{0}\) is needed for different problem settings. In practice, we only have access to a very limited amount of training data, for which both INN and FNN can achieve acceptable accuracy. However, INN allows training forward and inverse processes simultaneously using one single NN.

Now, we examine more closely the convergence behaviour of the Adam iterates for the training with 100, 500 and 1000 training data, see Fig. 3. We observe that both relative approximation and generalization errors for the forward and inverse processes decay rapidly at first several iterations and then reach a steady state. This behavior is commonly observed for many problems. Additionally it is observed that the generalization error \(e_{g}\) exhibits a semi-convergence phenomenon, especially for the inverse process: it first decreases steadily up to a finite number of iterations, and then starts to increase (occasionally in a dramatic manner) as the iteration further proceeds. Hence, suitable regularization might be necessary to overcome the phenomenon, e.g., early stopping or suitable explicit regularization. In contrast, the approximation error remains fairly stable throughout, indicating a stable optimization process. Compared with the inverse process, the forward process tends to be more stable numerically, since the forward map \(F\) enjoys a far more favorable stability estimate and the inverse map \(F^{-1}\) enjoys only conditional (Holder type) stability estimates (see [26] and references therein for details).

Finally, the results in Figs. 4 and 5 indicate that the training data size \(N_{t}=500\) is sufficient for problem (4.1) to achieve reasonable accuracy. A smaller training set yields lower approximation errors, but performs worse on the generalization; Larger training datasets lead to comparable performance in terms of both approximation and generalization. This observation agrees well with the intuition that within the supervised learning paradigm, often one needs a lot of training data in order to achieve good performance of the trained NNs.
