In the previous chapters, we have viewed certain neural networks as discretisations of continuous systems. In the following, we will return to the classical view of discrete networks, which are additionally endowed with a certain structure. \emph{Invertible neural networks}, i.e. bijective neural networks, are such an example. They have been an emerging topic in deep learning over the last few years. They are naturally connected to neural networks that are inspired by ordinary differential equations, as flows of ODEs are themselves invertible.
Two of the main applications of invertible neural networks lie in memory-efficient backpropagation as well as generative modeling with density estimation. For simplicity, in the following we will use the abbreviation $f_k=f^k(\cdot, \theta^k)$ when appropriate. 
Much like in general, invertible networks are typically parametrised as a function composition $\Psi=f_K \circ \dots \circ f_1$,  with the additional constraint that each layer $f_k$ is a bijective function. The inverse can then simply be computed via $\Psi^{-1}={f_1}^{-1} \circ \dots \circ {f_K}^{-1}$. A main restriction this imposes on the architecture is the fact that for each layer $f^k : \FeatureSpace^k \times \GenParamSet^k \to \mathcal{X}^{k+1}$, one has $\dim ({\mathcal{X}^k})=\dim ({\mathcal{X}^{k+1}})$, where the $\mathcal{X}^k$ are the feature spaces. As a consequence, the dimension of the input space determines the dimensions of the feature spaces in a (fully) invertible neural network. 


\subsection{Types of invertible layers}
Invertible networks and layers can roughly be divided into two categories: Those that are algebraically invertible and those that are inverted with a numerical approximation scheme.

\subsubsection{Coupling layers} \label{sec:coupling}
Coupling layers \cite{dinh2014nice} work by splitting a layer's input into two parts and applying a suitable transformation that is easily invertible for one of the two parts. Mathematically, for an $M$-dimensional input, the index set $I=\{1,\dots,M\}$ is partitioned into two sets. In convolutional networks, the partition is often done channel-wise, i.e. the channels are split into two sets. A different type of partitioning is \emph{invertible downsampling} (also known as a \emph{masking} or \emph{squeeze} operation \cite{dinh2014nice}). 
With $x = (x_{1},x_{2}) \in \mathbb{R}^{M_1} \times \mathbb{R}^{M_2}$ and $y = (y_1,y_2)  \in \mathbb{R}^{M_1} \times \mathbb{R}^{M_2}$, a coupling layer $g: x\mapsto y$ is defined via the mapping
\begin{equation}
    \begin{aligned}
        y_{1} &= x_{1} \\
        y_{2} &= \gamma (x_{2}, f(x_{1})), \\ 
    \end{aligned} \label{eq:forward_coupling}
\end{equation}
where $f: \R^{M_1} \to  \R^{M_2}$ and $\gamma: \R^{M_2} \times  \R^{M_1} \to \R^{M_2}$ is invertible in the first argument. Then $g^{-1}: y \mapsto x$ is given by 
\begin{equation}
    \begin{aligned}
        x_{1} &= y_{1} \\
        x_{2} &= \gamma^{-1}(y_{2}, f(y_{1})), \\ 
    \end{aligned} \label{eq:inverse_coupling}
\end{equation}
where $\gamma^{-1}$ is the inverse of $\gamma$ in its first argument (for fixed second argument). The \emph{coupling law} $\gamma$ can be as simple as $\gamma: (a,b) \mapsto a+b$, with which $f$ is called an \emph{additive coupling layer} \cite{dinh2014nice}, such that $\gamma^{-1}: (c,b) \mapsto c-b$. Another commonly-used class of coupling layers are \emph{affine coupling layers} \cite{dinh2014nice, dinh2016density}. More complex, but in theory more expressive coupling laws can be constructed from strictly monotonic splines \cite{durkan2019neural}. Note that \eqref{eq:forward_coupling} and \eqref{eq:inverse_coupling} are shear mappings (Section \ref{sec:shears}), which in the case of ODEs are used to construct e.g. symplecticity-preserving numerical solutions.

There is in principle no restriction on the function $f$. This allows for the utilisation of arbitrarily expressive sub-networks $f$ (e.g. a sequence of convolutional layers with non-linear activation functions), without rendering $g$ (algebraically) non-invertible. The numerical stability of this inversion may still pose an issue. Stability guarantees (both for the forward and the inverse mapping) can for instance be controlled via the Lipschitz constant of $f$ and $\gamma$, as shown in \cite{behrmann2020on}. This mirrors stability considerations of ODEs (see Section \ref{sec:structure_preserving_ODE_formulations}), where guarantees can be formulated as conditions on the Lipschitz constants. Note that while $f$ has to map to $\R^{M_2}$, the individual layers which comprise $f$ may change the dimensionality throughout this sub-network -- only the final layer has to transform the data to the required space $\R^{M_2}$. Since coupling layers perform learnable computations only on a part of the input, the partitioning should change throughout the network, e.g. by switching the roles of $x_1$ and $x_2$.

\subsubsection{Invertible layers through iterative schemes} \label{sec:invertible_by_numerical_schemes}
Another class of invertible layers are those that are invertible with an iterative scheme. Invertible residual networks \cite{invertibleresnet} are a special case of the commonly-used residual networks~\cite{He2016resnet}, which allow for the inversion with a simple fixed-point iteration. A residual layer can be framed as a function $g: \mathbb{R}^M \to \mathbb{R}^M$ with
\begin{equation}
    \begin{aligned}
        g(x) = x + f(x),
    \end{aligned} \label{eq:invertible_resnet_formula}
\end{equation}
where $f: \mathbb{R}^M \to \mathbb{R}^M$ is a sub-network. For fixed $x$, let $y:=x+f(x)$, such that $x= y - f(x)$. A fixed point $z^\ast$ of the function $\Phi_y: z \mapsto y-f(z)$ is thus in the preimage of $y$ under $g$, i.e. $g(z^\ast)=y$. Note that 
\begin{equation}
    \begin{aligned}
        \| \Phi_y(a) - \Phi_y(b) \| = \| f(b) - f(a) \| \leq \text{Lip}(f) \cdot \| b-a \|,
    \end{aligned}
\end{equation}
where $\text{Lip}(f)$ is the Lipschitz-constant of $g$. According to Banach's fixed-point theorem, if $\text{Lip}(f)<1$, then $\Phi_y$ is guaranteed to have a \emph{unique} fixed point, which is the inverse of $y$ under $g$, i.e. $z^\ast=g^{-1}(y)$. This inverse can be approximated to arbitrary precision via the iteration
\begin{equation}
    x^{i} = y - f(x^{i-1})
\end{equation}
for any initial value $x^0$. While most common layer types (such as dense or convolutional layers equipped with activation functions with bounded derivative) are Lipschitz, the condition $\text{Lip}(f)<1$ is not necessarily met. In \cite{invertibleresnet}, the required Lipschitz constraint is enforced by \emph{spectral normalisation}:   
Let the linear layer $A_\theta$ depend linearly on its parameters $\theta$ (e.g. convolutional layers without non-linear activation functions and biases; these are linear layers that linearly depend on their kernel). Then $A_{\tilde{\theta}}$ parametrised by 
\begin{equation}\bar{\theta}:= c \cdot \theta / \| A_\theta \|_2 \label{eq:spectral_normalisation}\end{equation}
has Lipschitz constant $\text{Lip} (A_{\bar{\theta}}) = c$. As a consequence, the layer $f(x)=\Activation (A_{\bar{\theta}}\cdot x + b)$ has Lipschitz constant $\text{Lip}(f) \leq c$ for any activation function $\Activation$ with $\text{Lip}(\Activation) \leq 1$, where $b$ is a bias vector. Thus, by updating the layers' parameters via spectral normalisation according to \eqref{eq:spectral_normalisation} after each gradient step, the required Lipschitz constraint for invertibility is guaranteed, if one chooses $c<1$. In practice, the spectral norm $\| A_\theta \|_2$ is calculated with the power method. 
To save computations, the authors in \cite{invertibleresnet} only perform a single power iteration, but re-use the estimation of the leading singular vector from the previous training step as the initial guess. While the power method (with a finite number of iterations) technically only provides a lower bound to $\| A_\theta \|_2$, the authors in \cite{invertibleresnet} find the Lipschitz constraint to still be met in practice. \\
In the case of invertible residual networks, the connection to neural networks as numerical solutions to ODEs is particularly strong. As noted in Section \ref{sec:intro_resnets_ode}, ResNets can generally be viewed as Euler discretisations of an ODE, if one views the activations, weights and biases as observations of time-dependent variables. Fittingly, in \cite{invertibleresnet} the authors originally motivate the inversion of such a ResNets by looking at the dynamics of the associated ODE backwards in time.

\subsubsection{Linear Invertible Components}\label{linearInvertible}
Above, two general approaches for constructing invertible, nonlinear layers were presented. In the following, we list a few \emph{linear}, invertible layers, which are used to increase the expressivity of invertible networks.

Coupling layers (Section \ref{sec:coupling}) work by dimension splitting, which in the context of convolutional neural networks usually consists in splitting the channels into two groups, which are processed independently from one another. This is in contrast to standard, multi-channel convolutions, where each output channel depends on all input channels. In order to overcome this limitation, the channels can subsequently be linearly combined (via an endomorphism). If the matrix of these coefficients is invertible, the automorphism created this way can be integrated as a linear layer into the invertible network framework. These invertible matrices can be parametrised in different ways:
One approach \cite{kingma2018glow} is to directly parametrise the desired matrix via a LU factorisation (with additional fixed permutation), which can then be inverted. While not discussed in the original publication, the diagonal entries themselves can be e.g. parametrised to be larger than some $\varepsilon>0$ to enforce stable invertibility. This is known as \emph{invertible} $1\times 1$ \emph{convolution}, because this 'channel mixing' can equivalently be realised as a convolution with (in 2D) 1-by-1 filters. An extension to this idea to \emph{convolutional} layers (with larger filters) is presented in \cite{pmlr-v97-hoogeboom19a}. 
A computationally particularly cheaply (and stably) invertible class of matrices are orthogonal matrices. In e.g. \cite{NIPS2019_8336}, these are parametrised as a sequence of Householder transformations. Other possible parametrisations of (special) orthogonal matrices include exponentials of skew-symmetric matrices, Cayley transforms or sequences of Givens rotations. For a discussion about the numerical optimisation of such classes of parameter matrices compare also the forthcoming Section \ref{sec:SPT} and in particular \eqref{eq:param1}.

Another type of linear, invertible layer are invertible down- and upsampling operations for image data. While classical down- and upsampling methods (such as bilinear interpolation or nearest-neighbour methods) from image processing are inherently non-invertible (as they change the dimensionality of the input), it is possible to construct \emph{invertible} up- and downsampling methods. Since the dimensionality of their output must be the same as the dimensionality of their input, decreasing (or increasing) the spatial dimensionality of an image \emph{must} be accompanied by a suitable increase (respectively decrease) in the number of channels. One such transformation for invertible downsampling is the \emph{pixel shuffle} transform \cite{dinh2014nice}, which subsamples the pixels and reorders them into new channels. The inverse of this transformation is an invertible upsampling operation. In \cite{etmann2020iunets}, this is generalised to learnable invertible up- and downsampling operations, which contain the (inverse) pixel shuffle as a special case. The general idea is to construct \emph{orthogonal} strided convolutional operators, where the kernel size matches the stride $s \in \mathbb{N}^d$, with spatial dimensionality $d$. It can be shown that by a specific reordering of an orthogonal $\tau$-by-$\tau$ matrix (where $\tau=s_1 \cdots s_d$) into a filter kernel, the resulting convolution is orthogonal. This implies that the inverse is simply given by the corresponding transposed convolution. The authors propose to parametrise the required orthogonal matrices via exponentiating skew-symmetric matrices. As the matrix exponential is a surjective mapping from the Lie algebra of skew-symmetric matrices to the special orthogonal matrices, \emph{any} special orthogonal matrix can be parametrised this way.

\subsection{Applications} 
\subsubsection{Memory-efficient backpropagation}\label{sec:memory_efficient_backprop} One possible area of application for invertible neural networks is memory-efficient backpropagation \cite{NIPS2017_6816}. Let $\HiddenFeature^{k+1}=\Layer(\HiddenFeature^{k},{\theta^k})$ be the output of a neural network's layer with nonlinear mapping $\Layer: \FeatureSpace^k \times \GenParamSet^k \rightarrow \FeatureSpace^{k+1}$, 
where $\theta^k$ are the layer's parameters. Let further $\Loss$ be the loss of the network. Then 
$$\nabla_{\theta^k} \Loss = \left( \partial_{\theta^k} \HiddenFeature^{k+1} \right)^\ast \nabla_{\HiddenFeature^{k+1}} \Loss = \left( \partial_{\theta^k} \Layer(\HiddenFeature^{k},\theta^k) \right)^\ast \nabla_{\HiddenFeature^{k+1}} \Loss$$
provides the weight-gradient necessary for the training of the network. For the calculation of this gradient, one needs both the gradient of the loss with respect to the output node (i.e. $\nabla_{\HiddenFeature^{k+1}} \Loss$), as well as the ability to calculate the derivative $\partial_{\theta^k} \Layer(\HiddenFeature^k, {\theta^k})$. Unless $\Layer$ is linear, the calculation of the derivative requires access to $\HiddenFeature^k$, meaning that $\HiddenFeature^k$ needs to be \emph{stored in memory}. This typically represents the bulk of the memory demand in training neural networks. If, however $\Layer$ is \emph{invertible}, one can simply calculate $\HiddenFeature^{k}$ from $\HiddenFeature^{k+1}$ via
$$\HiddenFeature^{k}=\Layer^{-1}(\HiddenFeature^{k+1},{\theta^k})$$ for the additional computational cost of calculating the inverse (where $f^{-1}$ is the inverse of $\Layer$ in its first argument). Thus, instead of storing activations in memory during the forward pass, intermediate activations can simply be successively reconstructed from the output layer's activation. The memory requirement of training an invertible network using this scheme is thus independent of the number of invertible layers.

\subsubsection{Invertible networks as sub-networks}
In practice, many classical applications of neural networks such as classification, regression or segmentation do not typically map between vector spaces of the same dimensionality. Hence, a bijective function that maps between these two spaces does not exist, which is why a fully invertible neural network typically cannot solve the desired task. It is however possible to use an invertible neural networks as a \emph{subnetwork} in a neural network. For example, if $\Psi: \FeatureSpace \to \FeatureSpace$ is an invertible network and $\eta: \FeatureSpace \to \LabelSpace$ is a suitable non-invertible layer, the combined network $\eta \circ \Psi: \FeatureSpace \to \LabelSpace$ can for instance be used for a classification task of predicting labels in $\LabelSpace$ from features $\FeatureSpace$. As demonstrated in \cite{jacobsen2018irevnet}, such a neural network can have competitive performance to a comparably-sized residual network \cite{He2016resnet} on ImageNet \cite{deng2009imagenet}. Adding an output layer which transforms between spaces mirrors the approach for discretised ODEs, see Section \ref{sec:intro_resnets_ode}.

The use of invertible networks as sub-networks for example allows for the utilisation of the memory efficient backpropagation for these sub-networks (cf. Section \ref{sec:memory_efficient_backprop}), while the respective gradients of non-invertible sub-networks can be calculated conventionally. 

\subsubsection{Density estimation and generative modeling}
Aside from generative adversarial networks (GANs) and variational autoencoders (VAEs), \emph{normalising flows} are another class of machine learning models that can be used to artificially generate data. While e.g. GANs are able to generate realistic-looking images \cite{karras2019style}, they lack the ability to estimate the likelihood of data under the generative model at hand. Likewise, VAEs can only estimate a lower bound of the likelihood -- the variational lower bound. This is in contrast to normalising flows, which are trained via maximum likelihood estimation. 

Like most generative models, normalising flows generate data from a simple base distribution (usually Gaussian) via a learned transformation. Let the random variable $z$ have probability density function $q$, for which we will write $z \sim q(z)$. For any diffeomorphism $f$, it holds that $$x:=f^{-1}(z) \sim q(z) \cdot \left| \det (\partial_z f^{-1}(z)) \right|^{-1}$$
due to the change-of-variables theorem\footnote{Here it suffices to consider invertible functions $f$ that are only \emph{locally} diffeomorphic in $x$. This is of practical relevance, because many neural networks are only locally differentiable due to the use of piecewise differentiable activation functions such as ReLU.} \cite{bogachev2007measure}. This means that for the probability density of $x$ (denoted $p$), the log-likelihood of $x$ can be expressed as
\begin{equation}
    \log p(x) = \log q(f(x)) + \log \left| \det (\partial_x f(x))  \right|. \label{eq:log_likelihood_normalising}
\end{equation}
Let $f$ be parametrised by an invertible neural network, i.e. $f(x)=\Psi(x,\theta)$ for all $x \in \FeatureSpace$. Given training data $(\Feature_n)_{n=1}^{N}$, the minimiser of
\begin{align}
    \min_{\theta \in \Theta} \left\{ E(\theta) = -\frac1N \left( \sum_{n=1}^N -\log q(\Psi(x_n,\theta)) + \log \left| \det (\partial_{x_n} \Psi (x_n,\theta))  \right| \right)\right\} . \label{eq1:normalising_flows_minimisation}
\end{align}
is a maximum likelihood estimator of the training data under the neural network. Models trained this way are called \emph{normalising flows}, because they are usually trained to convert complicated data into normal distributions. In this framework, artificial data can be generated by sampling $f^{-1}(z)$, where $z\sim q(z)$. An example of learning the 'two half-moons' dataset this way is provided in Figure \ref{fig:normalising_flow}. This case highlights that continuously transforming a normal distribution into a distribution whose probability density function has disconnected support is not possible. However, the normalising flow instead assigns low probability density to areas outside of the support, which results in a distribution that resembles the true distribution.

One of the main difficulties in designing normalising flows lies in the evaluation (respectively differentiation) of the determinant term. Note that for a neural network of the form \eqref{eq:neuralnet_general_form} it holds that $$\det (\partial_{z^0} \Psi (z^0,\theta))=\det (\partial_{z^{K-1}} z^K) \cdots \det  (\partial_{z^{0}} z^1),$$
such that the computation of the determinant for the whole network can be performed in a layer-by-layer fashion. Still, na\"{i}vely computing each determinant by the Leibniz formula yields prohibitive computation times. Depending on the types of layers used, the structure of the individial layers' Jacobian makes employing different determinant identities possible. If e.g. the Jacobian is lower triangular (as is the case with coupling layers \cite{dinh2014nice}), the determinant reduces to the product of the Jacobian diagonal entries. Invertible residual networks \cite{invertibleresnet} on the other hand induce no such structure. By using the Lipschitz constraint on $f$ in equation \eqref{eq:invertible_resnet_formula} and using the identity from \cite{withers2010log}, the authors show that in their case
$$\log \left| \det (\partial_{z^0} \Psi (z^0,\theta)) \right| = \text{trace}\left( \log \partial_{z^0} \Psi (z^0,\theta)\right)$$
holds. The matrix logarithm is then approximated as a truncated Taylor series, whereas for the evaluation of the trace, the Hutchinson estimator \cite{hutchinson1990stochastic} is employed. An extension to this concept is presented in \cite{chen2019residual}, where the truncation with a fixed number of steps is replaced by a 'Russian roulette' estimator, a type of Monte-Carlo estimator. This estimator introduces a stochastic truncation of the series, which results in an unbiased estimator of the log-likelihood \eqref{eq:log_likelihood_normalising}.\\
In the context of normalising flows, another connection to ODEs becomes apparent: As noted in \cite{chen2018neural}, a continuous version of the change-of-variables theorem may be used to formulate a continuous normalising flow, which in turn can be solved numerically via a discretisation scheme. A maximum likelihood estimator can in turn be obtained by training a Neural ODE \cite{chen2018neural}.

\begin{figure}
    \centering
    \includegraphics[width=.3\textwidth]{pics/moons_real.png}\hspace{2mm}\includegraphics[width=.3\textwidth]{pics/moons_fake.png}\hspace{2mm}\includegraphics[width=.3\textwidth]{pics/moons_generated.png}
    \caption{Left: 'Two half-moons' dataset. Middle: An approximation of the probability density function of the 'two half-moons' dataset by a normalising flow, generated by the code from \cite{chen2019residual}. Right: Points sampled from the approximated distribution.}
    \label{fig:normalising_flow}
\end{figure}

\subsection{Open problems}
\subsubsection{Fast inverses without coupling}
Coupling layers (Section \ref{sec:coupling}) allow both for quick forward and inverse computations. Their expressivity is however hindered by the fact that expressive, learned transformations are only applied to a part of their input, while the remaining portion of the input is not transformed. However, while this problem can be mitigated by using multiple coupling layers (with different partitions of the input), it may be desirable to construct invertible layers, where each output neuron depends on all input neurons (as is the case with 'conventional' fully-connected or convolutional layers). Invertible residual layers (Section \ref{sec:invertible_by_numerical_schemes}) on the other hand fulfill this criterion, but they rely on a numerical inversion via a fixed-point iteration, which requires multiple evaluations of the forward operation. Furthermore, in order to guarantee their invertibility, a Lipschitz constant needs to be controlled, which requires additional computation time. For practical purposes, it would be desirable to combine both methods' advantages in order to obtain both fast and expressive invertible layers.

\subsubsection{Stability guarantees}
As discussed in \cite{behrmann2020on}, controlling the Lipschitz constant of invertible layers is a conceptionally simple, but practically costly method of guaranteeing a certain stability of the inverse. Furthermore, as discussed above, they need to be controlled in order to guarantee the convergence of the iterative scheme for the inversion of residual layers. In the interest of computational efficiency, is there a computationally cheaper way to control the Lipschitz constant than via the proposed power method? Work in this direction has already been performed in \cite{chen2019residual}, as the authors experiment with Lipschitz constant corresponding to mixed norms. 

In a similar vein, are there alternative ways of guaranteeing stability other than by controlling the Lipschitz constant?
{"config": "arxiv", "file": "2006.03364/invertible.tex"}