Year: 2019

Abstract:
Guarantees in deep learning are hard to achieve due to the interplay of flexible modeling schemes and complex tasks. Invertible neural networks (INNs), however, provide several mathematical guarantees by design, such as the ability to approximate non-linear diffeomorphisms. One less studied advantage of INNs is that they enable the design of bi-Lipschitz functions. This property has been used implicitly by various works to design generative models, memory-saving gradient computation, regularize classifiers, and solve inverse problems. In this work, we study Lipschitz constants of invertible architectures in order to investigate guarantees on stability of their inverse and forward mapping. Our analysis reveals that commonly-used INN building blocks can easily become noninvertible, leading to questionable “exact” log likelihood computations and training difficulties. We make use of numerical analysis tools to diagnose non-invertibility in practice. Finally, based on our theoretical analysis, we show how to guarantee numerical invertibility for one of the most common INN architectures.

heading: 1 INTRODUCTION
text: Invertible neural networks (INNs) have become a standard building block in the deep learning toolkit. Invertibility is useful for training generative models with exact likelihoods (Dinh et al., 2014; 2017; Kingma & Dhariwal, 2018; Kingma et al., 2016; Behrmann et al., 2019; Chen et al., 2019), increasing posterior flexibility in VAEs (Rezende & Mohamed, 2015; Tomczak & Welling, 2016; Papamakarios et al., 2017), learning transition operators in MCMC samplers (Song et al., 2017; Levy et al., 2017), computing memory-efficient gradients (Gomez et al., 2017; Donahue & Simonyan, 2019), allowing for bi-directional training (Grover et al., 2018), solving inverse problems (Ardizzone et al., 2019) and analysing adversarial robustness (Jacobsen et al., 2019).
The application space of INNs is rapidly growing and many approaches for constructing invertible architectures have been proposed. A common way to construct invertible networks is to use triangular coupling layers (Dinh et al., 2014; 2017; Kingma & Dhariwal, 2018), where dimension partitioning is interleaved with ResNet-type computation. Another approach is to use various forms of masked convolutions, generalizing the dimension partitioning approach of coupling layers (Song et al., 2019; Hoogeboom et al., 2019). To avoid dimension partitioning altogether, multiple approaches based on efficiently estimating the log-determinant of the Jacobian, necessary for applying the change of variable formula, have been proposed to allow for free-form Jacobian structure (Grathwohl et al., 2019; Behrmann et al., 2019; Chen et al., 2019).
From a mathematical perspective, invertible architectures enable several unique guarantees like:
• Enabling flexible approximation of non-linear diffeomorphisms (Rezende & Mohamed, 2015; Dinh et al., 2017; Kingma & Dhariwal, 2018; Chen et al., 2019) • Memory-saving gradient computation (Gomez et al., 2017; Donahue & Simonyan, 2019) • Fast analytical invertibility (Dinh et al., 2014) • Guaranteed preservation of mutual information and exact access to invariants of deep
networks (Jacobsen et al., 2018; 2019).
Despite the increased interest in invertible neural networks, little attention has been paid to guarantees on their numerical invertibility. Specifically, this means analyzing their ability to learn bi-Lipschitz
neural networks, i.e. Lipschitz continuous neural networks with a bound on the Lipschitz constant of the forward and inverse mapping.
While the stability analysis of neural networks has received significant attention e.g. due to adversarial examples (Szegedy et al., 2013), the focus here is only on bounding Lipschitz constants of the forward mapping. However, bounding the Lipschitz constant of the inverse mapping is of major interest, e.g. when reconstructing inputs from noisy or imprecise features. In fact, analytical invertibility as provided by some invertible architectures does not necessarily imply numerical invertibility in practice.
In this paper, we first discuss the relevance of controlling the bi-Lipschitz bounds of invertible networks. Afterwards we analyze Lipschitz bounds of commonly used invertible neural network building blocks. Our contributions are:
• We argue for forward and inverse stability analysis as a unified viewpoint on invertible network (non-)invertibility. To this end, we derive Lipschitz bounds of commonly-used invertible building blocks for their forward and inverse maps. • We numerically monitor and detect (non-)invertibility for different practical tasks such as
classification and generative modeling. • We show how this overlooked issue with non-invertibility can lead to questionable claims
when computing exact likelihoods with the change-of-variable formula. • Finally, we study spectral normalization as a stabilizer for one of the most commonly-used
family of INN architectures, namely additive coupling blocks.

heading: 2 BACKGROUND AND MOTIVATION
text: Invertible neural networks are bijective functions with a parametrized forward mappingFθ : Rd → Rd with Fθ : x 7→ z, where θ ∈ Rp defines the parameter vector. Additionally, they define an inverse mapping F−1θ : Rd → Rd with F −1 θ : z 7→ x. This inverse can be given in closed-form (analytical inverse, e.g. Dinh et al. (2017); Kingma & Dhariwal (2018)) or approximated numerically (numerical inverse, e.g. Behrmann et al. (2019); Song et al. (2019)).
Before we discuss building blocks of invertible networks, we provide some background and motivation for studying forward and inverse stability.
Definition 1 (Lipschitz and bi-Lipschitz continuity). A function F : (Rd1 , ‖ · ‖) → (Rd2 , ‖ · ‖) is called Lipschitz continuous if there exists a constant L =: Lip(F ) such that
‖F (x1)− F (x2)‖ ≤ L‖x1 − x2‖, ∀x1, x2 ∈ Rd1 .
If an inverse F−1 : (Rd2 , ‖ · ‖)→ (Rd1 , ‖ · ‖) and a constant L∗ =: Lip(F−1) exists such that
‖F−1(y1)− F−1(y2)‖ ≤ L∗‖y1 − y2‖, ∀y1, y2 ∈ Rd2 ,
then F is called bi-Lipschitz continuous.
Remark 2. We focus on invertible functions F : (Rd, ‖ · ‖2)→ (Rd, ‖ · ‖2), i.e. functions where the domain and co-domain are of the same dimensionality d and the norm is given by the euclidian norm.
Lemma 3. (Rademacher (Federer, 1969, Theorem 3.1.6)) If F : Rd → Rd is a locally Lipschitz continuous function (i.e. functions whose restriction to a neighborhood around any point is Lipschitz), then F is differentiable almost everywhere. Moreover, if F is Lipschitz continuous, then
Lip(F ) = sup x∈Rd
‖JF (x)‖2,
where JF (x) is the Jacobian matrix of F at x and ‖JF (x)‖2 denotes its spectral norm.
Lipschitz bounds on the forward mapping are of crucial importance in several areas, including in adversarial example research (Szegedy et al., 2013), to avoid exploding gradients, or the training of Wasserstein GANs (Anil et al., 2019). The stability of the inverse, however, can have a similar impact. For instance, having a Lipschitz bound on the inverse may avoid vanishing gradients during training.
Given that deep-learning computations are carried out with limited precision, imprecision is always introduced in both the forward and backward passes, i.e., zδ = F (x) + δ and x̂δ = F−1(zδ). Instability in either pass will aggravate this problem, and essentially make the invertible network numerically non-invertible. To summarize, this problem occurs in the following situations:
• Numerical reconstruction of x, where features zδ are inexact due to limited precision (e.g. when computations are executed in single precision as common on modern hardware). • Reconstruction based on imprecise measurements from physical devices (e.g. when using
invertible networks for inverse problems (Ardizzone et al., 2019)). • Numerical re-computation of intermediate activations of the neural network to allow for
memory-efficient backpropagation (Gomez et al., 2017).
Furthermore, some computations are performed via numerical approximation, which in turn adds another source of imprecision that might be aggravated via instability. Examples include:
• Numerical forward computation, as in Neural ODEs (Chen et al., 2018) (numerical solver is used to approximate dynamic of ODE). • Numerical inverse computation, e.g. via fixed-point iterations as in i-ResNets (Behrmann
et al., 2019) or MintNet (Song et al., 2019) or via ODE-solvers for the backward dynamics as in Neural ODEs (Chen et al., 2018).
As an example of why bi-Lipschitz continuity is critical for numerical stability in invertible functions, let’s consider the simple mappings F1(x) = log(x), F−11 (z) = exp(z), and F2(x) = x, F −1 2 (z) = z. Though both functions tend to infinity when x → ∞ , F1 is much less stable. Consider the introduction of numerical imprecision as zδ = F1(x)+ δ where δ denotes the introduced imprecision. Then this imprecision is magnified in the inverse pass as:
||F−11 (z)− F −1 1 (z δ)||22 ≈ ||δ ∂F−11 (z δ)
∂zδ ||22 = ||δ exp
( zδ ) ||22. (1)
A similar example can be constructed for both the forward and backward passes, which speaks to the importance of bi-Lipschitz continuity. For an additional discussion on the connection of Lipschitz constants and numerical errors, we refer to Appendix B.

heading: 3 STABILITY OF INVERTIBLE NEURAL NETWORKS

heading: 3.1 LIPSCHITZ BOUNDS FOR BUILDING BLOCKS OF INVERTIBLE NETWORKS
text: Research on invertible networks has produced a large variety of architectural building blocks. Yet, the focus of prior work was on obtaining flexible architectures while maintaining invertibility guarantees. Here, we build on the work in (Behrmann et al., 2019), where bi-Lipschitz bounds were proven for invertible ResNets, by deriving Lipschitz bounds on the forward and inverse mapping of common building blocks. Together with an overview of common invertible building blocks, we provide our main results in Table 1. We chose these particular model classes in order to cover both coupling-based approaches and free-form approaches like Neural ODE (Chen et al., 2018) and i-ResNets (Behrmann et al., 2019). The derivations of the bounds are given in Appendix A. Note that the bounds provide the worst-case stability and serve mainly as a guideline for future designs of invertible building blocks.

heading: 3.2 CONTROLLING STABILITY OF BUILDING BLOCKS
text: As shown in Table 1, there are many factors that influence the stability of INNs. Of particular importance are the Lipschitz constants Lip(g) of the sub-network g for i-ResNets (Behrmann et al., 2019) and affine coupling blocks (Dinh et al., 2014), and Lip(s), Lip(t) for additive coupling blocks (Dinh et al., 2017). Whereas computing the Lipschitz constants of neural networks is NP-hard (Virmaux & Scaman, 2018), there is a simple data-independent upper bound:
Lip(g) ≤ L∏ i=1 ‖Ai‖2, for g(x) = AL ◦ φ ◦AL−1 ◦ · · · ◦A2 ◦ φ ◦A1, (2)
where Ai are linear layers, ‖ · ‖2 is the spectral norm and φ a contractive activation function (Lip(φ) ≤ 1). The above bound was used by (Behrmann et al., 2019) in conjunction with spectral normalization (Miyato et al., 2018; Gouk et al., 2018) to ensure a contractive residual block g. In particular, this employs a normalization via:
Ã = κ A
σ̂1 , with σ̂1 ≈ σ1 = ‖A‖2 (approx. via power-method),
where κ > 0 is a coefficient that sets the approximate upper bound on the spectral norm of each linear layer Ai. Thus, by setting an appropriate coefficient κ depending on the targeted Lipschitz bound of the building block, this approach enables one to control both forward and inverse stability. Note that the above discussion can be generalized to other `p-norms, see (Chen et al., 2019).
However, this is not sufficient when using affine coupling blocks because their bound on the Lipschitz constant holds only locally. In particular, it depends on the regions of the inputs x to the coupling block. While inputs to the first layer are usually bounded by the nature of the data, obtaining bounds for intermediate activations is less straightforward. One interesting avenue for future work could be local regularizers like gradient penalties (Gulrajani et al., 2017), where spectral normalization could be used post-hoc to certify stability.
Lastly, we use ActNorm (Kingma & Dhariwal, 2018) in several architectures and avoid small diagonal terms which would yield large Lipschitz constants in the inverse (see Table 1) by adding a positive constant. Further stabilization could be achieved via bounding the scaling.

heading: 4 NUMERICAL EXPERIMENTS
text: In this section, we study numerical invertibility for several objectives and architecture settings. This section is structured by task:
1. Classification: we show that INN classifiers can become non-invertible on CIFAR-10 and discuss consequences for memory-efficient backpropagation.
2. Density estimation: we analyze the numerical invertibility of SOTA trained density models. 3. Generative modeling: we study the stability of adversarially trained INN generators,
discuss the consequences for likelihood evaluation, and stabilize an additive-coupling based INN generator using spectral normalization.
4. Decorrelation: we perform an in-depth study of the effect of different architecture settings on a simple task, where both stable and unstable solutions are possible. Furthermore, we show that spectral normalization is effective at stabilizing additive-coupling based flows.
We use the following measures to diagnose the numerical instability of invertible models:
• Reconstruction error. We measure the `2-distance between the input x and its reconstruction, i.e. ||x(i) − F−1θ (Fθ(x(i)))||2. • Conditioning of the Jacobian and max/min singular values. For forward stability, we
are interested in the behavior of the Jacobian JF (x), while for inverse stability we are interested in the Jacobian of the inverse mapping JF−1(x). We compute the singular values of the Jacobians using the SVD, which allows us to compute its condition number. 1
While the reconstruction error allows us to quantitatively monitor non-invertibility even before reconstruction artifacts are perceptible, the linear approximation JF (x) and its singular values provide insights into unstable directions of the forward (very large singular values of JF (x)) and inverse (very small singular values of JF (x)) mapping. Both measures were also used in (Jacobsen et al., 2018), where an ill-conditioned inverse was observed.
4.1 CLASSIFICATION WITH INVERTIBLE MODELS
In this section, we show that when training an INN for classification, there exist both stable and unstable solutions. We compare a stable model that uses additive coupling and an unstable model that uses affine coupling and ActNorm both inside and between the blocks—for additional experimental details see Appendix C. These models achieve similar test accuracies of 90.2% and 90.5%, respectively (Figure 6, Appendix C); however, we note that the goal of this experiment is not to achieve SOTA accuracy on CIFAR-10. Rather, we aim to show that models trained for classification with reasonable accuracy, can vary greatly with respect to stability. To observe the differences in stability, we plot the reconstruction results in Figure 2.
An important use-case of INNs is to enable memoryefficient training, by re-computing activations in the backward pass rather than storing them in memory during the forward pass (Gomez et al., 2017). This approach enables e.g. large-scale generative modeling (Donahue & Simonyan, 2019) and scaling segmentation networks to high-resolution medical images (Brügger et al., 2019).
Re-computing the activations, however, relies on a numerically precise inverse mapping. To better understand the effect of numerical errors, we perform an analysis similar to Gomez et al. (2017): we track the angle between the true and memory-saving gradients during training (Figure 1). As expected from the reconstruction results, we observe that the affine model yields gradients very different from the true gradient; in fact, after approximately 20 epochs of training, the memory-saving gradients of the affine model contain numerically infinite or nan values. Thus, it would not be possible to train the affine model successfully using memory-saving gradients.
1 Using the SVD is feasible for CIFAR-10, but becomes prohibitively expensive for ImageNet (with images of size 3× 256× 256); for such larger images, one could instead use the Lanczos algorithm (Lanczos, 1950) to find the largest and smallest singular values, to compute the condition number.

heading: 4.2 ANALYZING STATE-OF-THE-ART DENSITY MODELS
text: In this section, we analyze the invertibility of trained density models. In particular, we expose non-invertibility in models that otherwise appear stable by optimizing in input space to find examples that are poorly reconstructed by the model. Here, we take a trained Glow model (Kingma & Dhariwal, 2018)2 and optimize the input using Projected Gradient Descent (PGD) (Madry et al., 2018). Our goal is to find a point x′ in the domain of the invertible model such that the reconstruction F−1(F (x′)) differs from x′. In particular, we start with a datapoint x and use PGD to find a perturbed example x′ that has high reconstruction error via Eq. 3 (additional details in Appendix E).
argmax ||x′−x||∞≤
||x′ − F−1(F (x′))||2. (3)
As shown in Figure 3, this attack is effective for finding examples that are perceptually identical to test examples, yet induce large reconstruction errors. This attack can be understood as a worstcase invertibility diagnosis; however, we note that unsuccessful attacks can be due to algorithmic
2We used the PyTorch implementation from https://github.com/y0ast/Glow-PyTorch. This trained model achieves a likelihood of 3.39 bits-per-dimension.
issues and thus do not necessarily imply stable invertible models. Also, it is not always clear how to get gradients that are able to exploit numerical instabilities, leaving room for improvement via gradient-free methods (e.g., the boundary attack from Brendel et al. (2018)).
We present in Appendix E additional experiments using the PGD attack on the additive Glow model from (Kingma & Dhariwal, 2018) trained on CelebA, where we show that the model becomes non-invertible outside the valid input range of images. Furthermore, we analyze a trained residual flow model (Chen et al., 2019) on CIFAR-10, where we cannot find such dramatic non-invertible inputs, which is to be expected given the stability bounds of i-ResNet blocks (see Table 1).

heading: 4.3 GENERATIVE MODELING WITH INVERTIBLE MODELS
text: Generative models based on invertible networks have been predominantly trained using maximum likelihood estimation (MLE). Another viable approach is to train them adversarially (ADV), as done in Flow-GAN (Danihelka et al., 2017; Grover et al., 2018). Flow-GAN is appealing as it can result in a generator capable of producing high-quality samples (as in GANs), while also giving access to exact density estimates, which GANs lack. Prior work (Danihelka et al., 2017; Grover et al., 2018) has compared these two techniques for training flows (MLE vs ADV); the main conclusion of these studies was that training with MLE yields good likelihoods but relatively poor samples, while training with a GAN loss yields good samples but likelihoods orders of magnitude worse than MLE training.
Here, we analyze in depth the effect of Flow-GAN training on numerical stability. We use networks with repeated additive coupling layers, and ActNorm between blocks. We examine two architectures, both with 3 levels, i.e., ‘squeeze’ between levels (additional details in Appendix D):
• Stable: having a depth of 4 (i.e., 4 blocks per level) and spectral normalization applied to all the convolution layers, see section 3.2 for details. • Unstable: having a depth of 16, ActNorm within coupling layer, and no spectral normaliza-
tion applied to the convolutional layers.
Table 2 shows that models trained with MLE objective can achieve good bits-per-dimenstion (BPD), but models trained with ADV can achieve better sample quality as measured by the Frechet Inception Distance (FID), a common measure of sample quality (Heusel et al., 2017; Lucic et al., 2018). This justifies why considering training INN with objective functions other than MLE is desired.
Broken Flow-GAN. In Figure 4, we show that an INN trained only with adversarial loss can become non-invertible, depending on the architecture. We perform forward and inverse passes repeatedly on the same mini-batch. The unstable model shows visible reconstruction errors quickly. Table 3 shows that the unstable model has BPD orders of magnitude larger than the stable model.
Flow-GAN “Likelihood”. Typically likelihood is computed by a change of variables, which assumes invertibility. When a network is numerically non-invertible, the assumption breaks, and the computed value becomes some numerical approximation to the density. The model used in Figure 4 consists of additive coupling blocks, ActNorm, and squeezing operations. All these operations have data independent log-determinant Jacobian. Thus, obtaining a numerical value via the change of variable formula is straightforward. However, it is unclear what this numerical value represents, likely it cannot be trusted as true likelihood due to the lack of invertibility. In this case, we advocate for ensuring invertibility using the remedies discussed here to make sure BPD values are trustworthy. In terms of sample quality, however, we currently observe a tradeoff since the stable model yields higher FID scores than the unstable model (Table 2, includes also MLE-models for comparison). Yet, we believe that proper tuning e.g. of spectral normalization could remove this tradeoff.
Lastly, one can adjust for the likelihood by adjusting the prior (see Appendix I). In summary, in this section we point out that Flow-GAN can become non-invertible, in which case the computed likelihood cannot be taken as ground truth likelihood (Grover et al., 2018). In sum, INN trained with MLE are stable, but the sample quality is worse than those trained with ADV. Yet, training with ADV loss might make INN non-invertible, which defeats the purpose of using INN in the first place. Hence, when training with alternative objectives, numerical stability is a crucial property to consider.

heading: 4.4 DECORRELATION TASK
text: As the last part of our empirical study, we use a simple decorrelation task to benchmark the stability of invertible models. In particular, we compute the correlation matrix C ∈ Rd×d via
Cθj,l = 1
N N∑ i=1
( Fθ(x (i))j − µ̂j ) ( Fθ(x (i))l − µ̂l )
σ̂j σ̂l , (4)
where µ̂j is the estimated mean over output samples Fθ(x(i)) and σ̂j the estimated standard deviation. Then, we optimize the parameters θ to minimize the off-diagonal correlation, i.e.
min θ ‖Cθ − diag(Cθ)‖F , (5)
where ‖ · ‖F is the Frobenius-norm 3. Decorrelation objectives have been used in (Cogswell et al., 2015) to reduce overfitting and in (Cheung et al., 2014) to disentangle hidden activations.
This objective serves as a good task for our purposes for two reasons:
1. Decorrelation is a simpler objective than optimizing outputs z = Fθ(x) to follow a factorized Gaussian as in Normalizing Flows (Rezende & Mohamed, 2015). Furthermore, it will show that changing the objective to less standard tasks can lead to larger instabilities compared to using INNs on more common tasks such as density estimation.
2. Decorrelation allows multiple solutions using invertible mappings, where both stable and unstable transforms are equally valid for the given objective. See Appendix F for a motivation based on a simple 2D toy example.
In summary, the decorrelation objective offers an environment to study which INN components steer the mapping towards stable or unstable solutions, that are equally plausible for the given task.
In our experiments shown in Figure 5, we focus on coupling-based models like Glow (Kingma & Dhariwal, 2018), which are analytically invertible and thus allow to a simpler analysis compared to models relying on numerical inversion like i-ResNets (Behrmann et al., 2019). We evaluate the effects of different architectural choices on numerical stability, including additive vs. affine coupling layers, ActNorm, and architecture depth. For ActNorm we study two settings: 1) between coupling blocks, 2) inside blocks, i.e. as part of the function g in additive blocks or s and t in affine blocks. Details on the architectures/ training schemes and extended results are provided in Appendix H.
3We provide example PyTorch code for the decorrelation objective in Appendix G.

heading: 5 RELATED WORK
text: Invertibility and stability of deep networks. The inversion from activations in standard neural networks to inputs has been studied in various works, e.g. via optimization in input space (Mahendran & Vedaldi, 2014). Linking invertibility and inverse stability for relu-networks was e.g. done in Behrmann et al. (2018). However, few works study the stability of INNs: Gomez et al. (2017) study the numerical errors in the gradient computation when using their memory-efficient backpropagation variant. Similarly to our empirical analysis, (Jacobsen et al., 2018) computed the SVD of the Jacobian of a trained i-RevNet and observed an ill-conditioned Jacobian. Lastly, the i-ResNet architecture (Behrmann et al., 2019) yields bi-Lipschitz bounds by design.
On the other hand, the stability of neural networks has been of major interest due to the problem of exploding and vanishing gradients, and more recently due to adversarial examples (Szegedy et al., 2013) and training of Wasserstein GANs (Arjovsky et al., 2017). See e.g. (Anil et al., 2019) for a promising approach to learn flexible Lipschitz neural networks.
Invertible building blocks. Besides the invertible building blocks we studied in Table 1, several other approaches were proposed. Most prominently, autogressive models like MAF (Papamakarios et al., 2017) or IAF (Kingma et al., 2016) provide invertible models that are not studied in our analysis. Furthermore, several newer coupling layers that require numerical inversion have been introduced (Jaini et al., 2019; Durkan et al., 2019). Besides the coupling-based approaches, multiple approaches (Chen et al., 2018; Behrmann et al., 2019; Chen et al., 2019; Song et al., 2019) use numerical inversion schemes, where the interplay of numerical errors due to stability and errors due to the numerical approximation of the inverse adds another dimension to the study of invertibility.
Fixed-Point arithmetic and limited precision. Maclaurin et al. (2015); MacKay et al. (2018) implement invertible computation using fixed-point numbers, with specially-designed schemes to store information that is “lost” when bits are shifted due to multiplication/division, enabling exact invertibility at the cost of additional memory usage. As Gomez et al. (2017) point out, this approach allows exact numerical inversion when using additive coupling blocks independent of stability issues. However, our stability analysis aims for a broadly applicable methodology beyond the special case of additive coupling. Lastly, there may be connections to deep learning using limited precision, see e.g. (Gupta et al., 2015), which could provide more insights into our observed numerical errors.

heading: 6 CONCLUSION
text: Numerical instability is an important concern for the practical application of invertible models. If for instance analytical invertibility does not carry through to the numerical computation due to instabilities or numerical errors, the consequences can be arbitrarily severe. As shown in our experiments, this can impact memory-efficient backpropagation (Gomez et al., 2017) and thus significantly reduces the usability of invertible networks if not handled appropriately. Flow-GAN illustrates another application where non-invertibility poses a serious threat, as instabilities can strongly influence or even break likelihood-computation.
In this paper, we shed light on the underlying causes of instability by deriving Lipschitz bounds on many of the atomic building blocks commonly used to construct INNs. From a practical standpoint, we used diagnostics to measure stability and provided an empirical framework to benchmark stability. Further, we have shown how to guarantee stability for one of the most common INN architectures. We hope that this will inspire future work to view numerical stability as a crucial axis in the design of new building blocks and architectures for invertible neural networks.

heading: A DERIVATIONS OF LIPSCHITZ BOUNDS
text: The bounds for invertible ResNets are taken from (Behrmann et al., 2019). For Neural ODEs (Chen et al., 2018), one needs to consider a Lipschitz constant Lip(F ) that holds for all t ∈ [0, T ], i.e.
‖F (t, x1)− F (t, x2)‖2 ≤ Lip(F )‖x1 − x2‖2, for all t ∈ [0, T ].
Then, the claimed bound is a standard result, see e.g. (Ascher, 2008, Theorem 2.3). Note that the inverse is given by dy(t)dt = −F (y(t), t), hence the same bound holds. In the subsequent subsections, we derive the bounds for coupling layers.
A.1 DERIVATION OF LIPSCHITZ BOUND FOR ADDITIVE COUPLING LAYERS
Consider an additive coupling block defined as
F (x)I1 = xI1 F (x)I2 = xI2 + g(xI1),
where I1, I2 is a disjoint partition of indices {1, ..., d} of the same cardinality, i.e. |I1| = |I2| = d2 . Further, xI1 , xI2 correpsonds to the corresponding dimension of x ∈ Rd and g : R d 2 → R d2 . By Lemma 3, it is
Lip(F ) = sup x∈Rd
‖JF (x)‖2.
Thus, in order to obtain a bound on the Lipschitz constant, it is helpful to look into the structure of the Jacobian. If the partitions I1 and I2 correspond to the first and last d2 indices, the Jacobian has a lower-block structure with an identity diagonal, i.e.
JF (x) =
( I 0
Jg(x) I
) .
By using this structure, we can derive the following upper bound:
Lip(F )2 = sup x∈Rd ‖JF (x)‖22
= sup x∈Rd sup ‖x∗‖2=1
‖JF (x)x∗‖22
= sup x∈Rd sup ‖x∗‖2=1
‖(JF (x)x∗)I1‖22 + ‖(JF (x)x∗)I2‖22
= sup x∈Rd sup ‖x∗‖2=1
‖x∗I1‖ 2 2 + ‖x∗I2 + Jg(x)x ∗ I1‖ 2 2
≤ sup x∈Rd sup ‖x∗‖2=1 ‖x∗I1‖ 2 2 + (‖x∗I2‖2 + ‖Jg(x)x ∗ I1‖2) 2 (6)
= sup x∈Rd sup ‖x∗‖2=1
‖x∗I1‖ 2 2 + ‖x∗I2‖ 2 2 + 2‖x∗I2‖2‖Jg(x)x ∗ I1‖2 + ‖Jg(x)x ∗ I1‖ 2 2
= sup x∈Rd sup ‖x∗‖2=1
‖x∗‖22 + 2‖x∗I2‖2‖Jg(x)x ∗ I1‖2 + ‖Jg(x)x ∗ I1‖ 2 2
= sup x∈Rd sup ‖x∗‖2=1
1 + 2‖x∗I2‖2‖Jg(x)x ∗ I1‖2 + ‖Jg(x)x ∗ I1‖ 2 2
= sup x∈Rd sup ‖x∗‖2=1
1 + 2‖Jg(x)x∗I1‖2 + ‖Jg(x)x ∗ I1‖ 2 2
= sup x∈Rd sup ‖x∗‖2=1
( 1 + ‖Jg(x)x∗I1‖2 )2 = sup x∈Rd (1 + ‖Jg(x)‖2)2
⇒ Lip(F ) ≤ 1 + Lip(g).
Furthermore, the inverse of F can be obtained via the simple algebraic transformation (y := F (x))
F−1(y)I1 = yI1 F−1(y)I2 = yI2 − g(yI1).
Since the only difference to the forward mapping is the minus sign, the Lipschitz bound for the inverse is the same as for the forward mapping.
A.2 DERIVATION OF LIPSCHITZ BOUND FOR AFFINE COUPLING LAYERS
Since the structure of the forward and inverse mapping for affine coupling layers has some differences, we split the derivation of the Lipschitz bounds into two sections. First, we start with the forward mapping and then reuse several steps for the bounds on the inverse mapping.
A.2.1 DERIVATION FOR THE FORWARD MAPPING
Consider an affine coupling block defined as
F (x)I1 = xI1 F (x)I2 = xI2 g(s(xI1)) + t(xI1),
where g(·) 6= 0 for all XI2 and I1, I2 as before. The Jacobian for this operation has the structure
JF (x) =
( I 0
DI(xI2)Dg′(xI1)Js(xI1) + Jt(xI1) Dg(s(xI1))
) ,
where D are following diagonal matrices DI(xI2) = diag ( (xI2)1, . . . , (xI2)|I2| ) ,
Dg′(xI1) = diag ( g′(s(xI2)1, . . . , g ′(s(xI2)|I2|) )
Dg(s(xI1)) = diag ( g(s(xI2)1, . . . , g(s(xI2)|I2|) ) .
Denote
M(x) := DI(xI2)Dg′(xI1)Js(xI1) + Jt(xI1).
By using an analogous derivation as in equation 6 (up to the inequality sign), we get
Lip(F )2 ≤ sup x∈Rd sup ‖x∗‖2=1 ‖x∗I1‖ 2 2 +
( ‖Dg(s(xI1))x∗I2‖2 + ‖M(x)x ∗ I1‖2 )2 = sup x∈Rd max i∈[|I1|] (1, Dg(s(xI1)i)) 2 + 2 max i∈[|I1|] (Dg(s(xI1)i))‖M(x)‖2 + ‖M(x)‖22
≤ sup x∈Rd max i∈[|I1|] (1, Dg(s(xI1)i)) 2 + 2 max i∈[|I1|] (1, Dg(s(xI1)i))‖M(x)‖2 + ‖M(x)‖22
= sup x∈Rd ( max i∈[|I1|] (1, Dg(s(xI1)i)) + ‖M(x)‖2 )2
⇐⇒ Lip(F ) ≤ max i∈[|I1|] (1, Dg(s(xI1)i)) + sup x∈Rd ‖M(x)‖2.
Next, we will look into the structure of M(x) to derive a more precise bound. Since inputs x are assumed to be bounded as x ∈ [a, b]d, it holds
‖DI(xI2)‖2 ≤ max(|a|, |b|).
Furthermore, let the derivative g′ of the element-wise function g be globally bounded by c, i.e. supx∈R g ′(x) ≤ cg′ . Then, it is
‖Dg′(xI1)‖2 ≤ cg′ .
In a similar manner as in section A.1, the spectral norm of the Jacobian of the scale-function s and translation-function t can be bounded by their Lipschitz constant, i.e.
‖Js(xI1)‖2 ≤ Lip(s) ‖Jt(xI1)‖2 ≤ Lip(t).
By using above bounds, we obtain
sup x∈Rd
‖M(x)‖22 ≤ max(|a|, |b|) · c · Lip(s) + Lip(t).
If we further assume, that the elementwise-function g is globally upper bounded by cg and we insert above bounds, we obtain
Lip(F ) ≤ max(1, cg) + max(|a|, |b|) · cg′ · Lip(s) + Lip(t).
A.2.2 DERIVATION FOR THE INVERSE MAPPING
For the affine coupling block from section A.2.1, the inverse is defined as
F−1(y)I1 = yI1 F−1(y)I2 = (yI2 − t(xI1)) g(s(yI1)),
where g(·) 6= 0 for all XI2 , I1, I2 as before and denotes elementwise division. The Jacobian for this operation has the structure
JF (x) =
( I 0
M∗(y) D 1 g (s(xI1))
) ,
where D 1 g (s(xI1)) denotes a diagonal matrix, as before. Furthermore, M∗ is defined as
M∗(y) = DI(yI2)D( 1g ) ′(s(yI1))Js(yI1)−D( 1g )′(s(yI1))Js(yI1)DI(t(yI1))−D 1g (s(yI1))Jt(yI1),
where D( 1g ) ′(s(xI1)) also denotes a diagonal matrix. Using analogous arguments as in section A.2.1, we obtain the bound
Lip(F−1) ≤ max i∈[|I1|] (1, D 1 g (s(xI1)i)) + sup x∈Rd ‖M∗(x)‖2.
Hence, we need to further bound the spectral norm of M∗. First, assume that 1g , the derivative ( 1 g )′ and translation t is globally upper bounded by c 1
g , c( 1g ) ′ and ct respectively. Furthermore consider
bounded inputs y ∈ [a∗, b∗]d. Then we obtain the bound
sup x∈Rd
‖M∗(x)‖22 ≤ max(|a∗|, |b∗|) · c( 1g )′ · Lip(s) + c( 1g )′ · Lip(s) · ct + c 1g · Lip(t).
Hence, we can bound the Lipschitz constant of the inverse of an affine block as
Lip(F−1) ≤ max i∈[|I1|] (1, c 1 g ) + max(|a∗|, |b∗|) · c( 1g )′ · Lip(s) + c( 1g )′ · Lip(s) · ct + c 1g · Lip(t).

heading: B NUMERICAL ERRORS AND LIPSCHITZ CONSTANTS
text: In a general setting, connecting numerical errors e.g. due to floating point operations to Lipschitz constants of the underlying mapping in a quantitative manner is not straightforward. For example, numerical errors due to limited precision occurs when summing to floating point numbers. As discussed in (Gomez et al., 2017), this occurs in additive coupling layers and is one source of numerical errors we observe in our experiments.
To formalize the connection to the Lipschitz constant, consider the following two mappings:
F (x) = z, (analytical exact computation) Fδ(x) = z + δ =: zδ, (floating point inexact computation)
In order to bound the error in the reconstruction due to the imprecision in the forward mapping, let xδ1 = F −1(zδ). Now consider
‖x− xδ1‖2 ≤ Lip(F−1)‖z − zδ‖2 = Lip(F−1)‖δ‖2,
where the Lipschitz constant of the inverse is used to bound the influence of the numerical error in the forward mapping. However, similarly to the forward mapping, the inverse mapping can also be imprecise. Thus, we introduce
F−1δ (zδ) = xδ1 + δ2 := xδ2
to formalize the numerical error in the inverse mapping. Hence, we obtain the bound
‖x− (xδ1 + δ2)‖2 ≤ ‖x− xδ1‖2 + ‖δ2‖2 ≤ Lip(F−1)‖z − zδ‖2 + ‖δ2‖2 = Lip(F−1)‖δ‖2 + ‖δ2‖2,
where the numerical errors of the mapping are denoted via δ (forward) and δ2 (inverse). While obtaining quantitative values for δ and δ2 for a model as complex as deep neural networks is hard, above formalization still provides insights into a potential role of the inverse stability when reconstructing inputs.

heading: C ADDITIONAL DETAILS FOR CLASSIFICATION EXPERIMENTS
text: Both affine and additive coupling-based models have the same architecture, that consists of 3 levels, 16 blocks per level, and 128 hidden channels. Each level consists of a sequence of residual blocks that operate on the same dimensionality. Between levels, the input is spatially downsampled by 2× in both width and height, while the number of channels is increased by 4×. Each residual block consists of a chain of 3× 3, 1× 1, 3× 3 convolutions, with ReLU activations in between. We trained on CIFAR-10 for 200 epochs, using SGD with Nesterov momentum 0.9 and weight decay 5e-4, with initial learning rate 0.01, decayed by a factor of 5 at epochs 60, 120, and 160. We found that
using a smaller initial learning rate of 0.01 was important for training the INN classifier, as opposed to the standard initial learning rate of 0.1 used for ResNets (Zagoruyko & Komodakis, 2016). We used standard data augmentation (random cropping and horizontal flipping).

heading: D EXPERIMENTAL DETAILS FOR FLOW-GAN
text: For the experiments in Section 4.3, the ADV models are trained with standard binary cross entropy loss and the hyperparameters in Table 4, whereas the MLE models are trained with learning rate 1e-3.

heading: E EXTENDED RESULTS FOR CRAFTED, NON-INVERTIBLE INPUTS
text: PGD Setup. To find non-invertible inputs for Glow and Residual Flows, we used PGD (Eq. 3) with = 0.1 and step size 0.01. For the Glow model in Section 4.2, we consistently found inputs with severe reconstruction errors (as shown in Figure 3) in fewer than 10 PGD iterations. For the Residual Flow model analyzed in this section, we ran 200 iterations of PGD. In each iteration, the pixel values of the perturbed image were clipped to the valid input range that the respective model was trained on ([−0.5, 0.5] for Glow and [0, 1] for the Residual Flow).
Crafted Inputs for Residual Flows. We also applied the PGD attack from Section 4.2 to a Residual Flow (Chen et al., 2019) pre-trained on CIFAR-10 (Figure 7).4 We find that, while there are visible differences between the crafted input x′ and its reconstruction x̂′, the reconstruction errors are less severe for the Residual Flow compared to Glow (analyzed in Section 4.2).
Instability Outside the Range of Training Inputs. Here, we applied the PGD attack from Section 4.2 to a Glow model pre-trained on Celeb-A (Liu et al., 2015). 5 This model was trained on images normalized to the range [−0.5, 0.5]. While the PGD attack was not successful at finding adversarial inputs in [−0.5, 0.5], it succeeded when the range was increased to [−0.7, 0.7] (by using = 0.2 and not clipping the perturbed inputs to [−0.5, 0.5]), yielding the example shown in Figure 8. Thus, we found that invertible models can become numerically non-invertible on out-of-distribution data.
4We used the pre-trained model from https://github.com/rtqichen/residual-flows. 5We used the model from https://github.com/openai/glow.

heading: F MOTIVATING THE DECORRELATION TASK
text: To motivate the decorrelation task as simple toy environment for stability of invertible models consider the following task (and its visualization in Figure 9):
• Consider input data x that is distributed via a standard normal distribution, i.e. x ∼ N (0, I) (Figure 9 (left)). • Assume the data is transformed by a rotation matrix R and a diagonal matrix D, i.e. y = RDx (Figure 9 (2. from left)). • Goal: decorrelate transformed data y using an invertible mapping.
Since correlation is independent of scale (scaling by standard deviation of the data), at least the two solutions A2 = D−1RT (Figure 9 (2. from right)) and A1 = RT (Figure 9 (right)) are equally valid for the given decorrelation task. However, the conditioning of the mappings A1 and A2 can be largely different if the scaling matrix D has a high condition number. Hence, this task both offers a stable solution, namely A1, and a (potentially) unstable solution A2.
To conclude, decorrelation can allow multiple solutions with different stability. Hence, decorrelation is a natural simple task to study which solution the INN picks. Furthermore, guiding the network to a stable solution is a justified strategy for this task and it is not expected to harm performance.

heading: G DECORRELATION EXAMPLE CODE
text: Here we provide an example implementation of the decorrelation objective used in Section 4.4, that minimizes the norm of the off-diagonal entries in the correlation matrix.
Listing 1: Example PyTorch code to implement the decorrelation loss used in our experiments. z = model(img) z_flat = z.view(z.size(0), -1) z_flat = z_flat - z_flat.mean(dim=0) # Subtract mean z_flat = z_flat / (z_flat.std(dim=0) + 1e-8) # Standardize correlation = (torch.mm(z_flat.t(), z_flat) / (z_flat.size(0)-1)) loss = torch.norm(correlation - torch.diag(torch.diagonal(correlation))) optimizer.zero_grad() loss.backward() optimizer.step()

heading: H EXTENDED RESULTS FOR DECORRELATION
text: Decorrelation Experiment Details Here we provide additional details on the model architectures and training schemes we used in our numerical experiments.
For all the decorrelation experiments, we used a 3-level model with blocks of depth 16. We used Adam (Kingma & Ba, 2015) with fixed learning rate 1e-4 and no weight decay, and trained on mini-batches of size 64. The CIFAR-10 images were normalized to the range [-0.5, 0.5], and were dequantized with uniform noise in [0, 1e-6].
Effect of Model Depth. Furthermore, we investigated the effect of network depth on stability since its an additional influence factor besides the selection of each invertible building block. Starting with a 3-level additive model with ActNorm, we vary the depth of the blocks between {4, 16, 32} and train with the decorrelation objective. The quantitative reconstruction errors and condition numbers of the Jacobians are shown in Figure 10. As expected, deeper architectures become unstable faster than shallow ones.
Loss Plots. Here we show that all the model variants investigated in the decorrelation experiments achieve their objective, i.e., the loss decreased enough for the correlation matrices to be diagonal.
Evolution of Condition Numbers, Max & Min Singular Values for Decorrelation. Here we plot the condition numbers, maximum and minimum singular values during training for the decorrelation task. We include plots for all settings discussed in the main paper as well as the appendix.

heading: I REFITTING PRIOR IN FLOW-GAN
text: In general, if the model is not optimized with forward KL, DKL(Pdata||Pθ), as in the case when optimizing with maximum likelihood, we cannot be sure F (x), x ∼ Pdata is best fitted with a standard Normal. Hence, a reasonable strategy, without changing anything in the learned network, is to refit the prior parameters. Here we simply optimize for maximum likelihood (as typically done for flow models) while only fitting a diagonal variance in prior. This can be interpreted as increasing the entropy of our model. In Kingma & Dhariwal (2018), they observed the opposite phenomenon that a model trained with maximum likelihood generates better samples after decreasing the entropy in the prior. See Figure 21 for samples after refitting the prior.

Review summaries:
The paper claims that for invertible neural networks, mathematical guarantees on invertibility is not enough, and we also require numerical invertibility. To this end, the lipschitz constants/condition numbers of Jacobians of both the forward and inverse maps of invertible NNs based on coupling layers are examined mathematically and experimentally. The paper also displays cases that expose non-invertibility in these architectures via gradient-based construction of adversarial inputs, as well as a decorrelation benchmark task, and show that spectral normalization can be a remedy for stabilizing these flows.

This paper analyses the numerical invertibility of analytically invertible neural networks (INN). The numerical invertibility depends on the Lipschitz constant of the respective transformation. The paper provides Lipschitz bounds on the components of building blocks for certain INN architectures, which would guarantee numerical stability. Furthermore, this paper shows empirically, that the numerical invertibility can indeed be a problem in practice.