{'fraction_non_alphanumeric': 0.051383083213412825, 'fraction_numerical': 0.011489526942141209, 'mean_word_length': 3.9560763004471613, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 50, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
ICE-BeeM: Identifiable Conditional Energy-Based Deep Models

Despite the growing popularity of energy-based models, their identifiability properties are not well-understood. In this paper we establish sufficient conditions under which a large family of conditional energy-based models is identifiable in function space, up to a simple transformation. Our results build on recent developments in the theory of nonlinear ICA, showing that the latent representations in certain families of deep latent-variable models are identifiable. We extend these results to a very broad family of conditional energy-based models. In this family, the energy function is simply the dot-product between two feature extractors, one for the dependent variable, and one for the conditioning variable. We show that under mild conditions, the features are unique up to scaling and permutation. Second, we propose the framework of independently modulated component analysis (IMCA), a new form of nonlinear ICA where the indepencence assumption is relaxed. Importantly, we show that our energy-based model can be used for the estimation of the components: the features learned are a simple and often trivial transformation of the latents.


INTRODUCTION
A central question in unsupervised deep learning is how to learn nonlinear representations that are a perfect reconstruction of the true latent variable behind the data. This allows us to learn factors that are semantically meaningful, interpretable and useful for downstream tasks. However, this is a very difficult task: by definition, we never observe the latent variables; the only information directly available to us is given by the observed variables.
Nonetheless, a lot of effort went into the related goal of disentanglement (Burgess et al., 2018;Higgins et al., 2017;Esmaeili et al., 2018;Chen et al., 2018), showing promising empirical proof. These methods relied on introducing hyperparameters in the ELBO objective, and train a variational autoencoder (VAE) (Kingma and Welling, 2013) to learn disentangled representation. These approaches are somehow heuristic, and lacked theoretical identifiability guarantees. Later on, these models were proven to be non-identifiable, both empirically (Locatello et al., 2018) and formally (Khemakhem et al., 2020). In fact, this is not a novel problem, and dates to the early days of the theory of nonlinear Independent Component Analysis (ICA) (Hyvärinen and Pajunen, 1999).
Learning the true representations is only possible when the representation is identifiable: conversely, if multiple choices of representations map to the same density over the observations, then the representations are indistinguishable and therefore unidentifiable. The first sufficient conditions for identifiable models in deep neural networks were introduced by recent work in the theory of nonlinear ICA (Hyvärinen andMorioka, 2016, 2017;Hyvärinen et al., 2019). Identifiability was made possible by introducing an additionally observed random variable that controls the independence of the latent variables, which are then transformed into observations through a nonlinear and invertible mapping. The goal in nonlinear ICA is then to estimate the inverse transformation from observations to latents. The nonlinearity of the transformation introduces many degrees of freedom that render the unconditional nonlinear ICA problem ill-defined. Broadly speaking, the conditioning variable reduces the space of possible transformations just enough so that it's possible to estimate the true inverse transformation.
Most recently, Khemakhem et al. (2020) explicited a connection between nonlinear ICA and the framework of variational autoencoders. They showed how one can use the theory of nonlinear ICA to train identifiable VAEs, but also how to use VAEs to estimate the nonlinear ICA model. It was shortly followed by work by Sorrenson et al. (2020), where a similar connection was made to flow-based models (Rezende and Mohamed, 2015).
We extend this trend to a broad family of unnormalized conditional energy-based models (EBM), using insight from the nonlinear ICA theory. A first important contribution of this paper is to provide strong identifiability properties for a large family of conditional EBMs, while making less assumptions than previous work. Conditional EBMs define a conditional density to be the exponential of the negative energy function that depends on both the dependent and conditioning variables, which are both observed. The energy function is defined in two steps: we learn two feature extractors, parameterized by neural networks, one for each of the observed variables (dependent and conditioning); then, we set the energy function to be the dot-product of the learned features.
EBMs offer unparalleled flexibility, mainly because they do not require the modeled densities to be normalized nor easy to sample from. In fact, the energy model we suggest will have universal approximation capabilities. Furthermore, and more importantly, we can show that both feature extractors are identifiable under very mild conditions: they learn representations that are unique up to a linear transformation. In addition, by slightly altering the architecture of the model, we prove the linear transformation is essentially a permutation. Effectively, this makes our family of models very flexible and adaptable to practical problems. The family of models for which we prove identifiability is by far the least constrained of all identifiable deep models in the literature so far. We call this model Identifiable Conditional Energy-Based deep Models, or ICE-BeeM for short.
Our second contribution is to introduce the Independently Modulated Component Analysis (IMCA) framework: a deep latent variable model where the latents are nonindependent, with an arbitrary global dependency structure. To achieve identifiability, we assume that the latent densities are modulated in a factorizable way when conditioned by an auxiliary variable (such as time index, history, or another data source). This framework is a generalization of nonlinear ICA, and a thorough analysis of its identifiability is covered here. We show how our ICE-BeeM can estimate this generative model, thus connecting both the generative and non-generative views.
As a further, rather different application of our results, we show how identifiability of ICE-BeeM can be leveraged for transfer learning. In fact, we believe that the identifiability results are generally important for principled application of EBMs, whether for the purposes of disentanglement or otherwise.

IDENTIFIABLE CONDITIONAL ENERGY-BASED DEEP MODELS 2.1 MODEL DEFINITION
We introduce in this section a parametric model for conditional densities. Formally, we collect a dataset of observations of tuples (x, y), where x ∈ X ⊂ R dx is the main variable of interest, also called the dependent variable, and y ∈ Y ⊂ R dy is an auxiliary variable also called the conditioning variable. We model the conditional density of x given y using a family of conditional energy-based models (EBM).
To do so, first consider two feature extractors f θ (x) ∈ R dz and g θ (y) ∈ R dz , which we parameterize by neural networks, and θ is the vector of weights and biases. To alleviate notations, we will drop θ when it's clear which quantities we refer to. These feature extractors are used to define the conditional energy function: The parameter θ lives in the space Θ which is defined such that the normalizing constant Z is finite: Our family of conditional energy-based models has the form: p θ (x|y) = exp(−f θ (x) T g θ (y)) Z(y; θ) As we will see later, this choice of energy function is not restrictive, as our model has powerful theoretical guarantees: universal approximation capabilities and strong identifiability properties. There exists a multitude of methods we can use to estimate such model, like for instance Flow Contrastive Estimation (Gao et al., 2019) and Denoising Score Matching (Vincent, 2011).

IDENTIFIABILITY
As stated earlier, we want our model to learn meaningful representations of the dependent and conditioning variables. In particular, when learning two different models of the family (2) from the same dataset, we want the learned features to be very similar.
This similarity between representations is better expressed as equivalence relations on the parameters θ of the network, which would characterize the form of identifiability we will end up with for our energy model. This notion of identifiability up to equivalence class was introduced by Khemakhem et al. (2020) to address the fact that there typically exist many choices of neural network parameters θ that map to the same point in function-space. In our case, it is given by the following definition. Definition 1. Let ∼ f w and ∼ g w be equivalence relations on Θ defined as: where A and B are (d z × d z )-matrices of rank at least min(d z , d x ) 1 , and c and e are vectors.
Two parameters are thus considered equivalent if they parameterize two feature extractors that are equal up to a linear transformation. The subscript w in both equivalence relations stands for weak, as we will strengthen this relation in the next section.
Identifiability of the learned representations was the key motivation behind the design of the energy function. Loosely speaking, if the conditioning variable conditions the distribution of the dependent variable strong enough, then the learned features will be similarly constrained because they interact closely through the dot-product. This is stated rigorously by the following theorem, which we prove in Appendix A.1. Theorem 1 (Identifiable conditional EBMs). Assume that for any choice of parameter θ: 1 if A is not full rank, then the definition (3) isn't an equivalence relation because it is not necessarily symmetric; we can fix this by requiring that there exists A1, A2 such that f θ = A1f θ + c1 and f θ = A2f θ + c2, but we present the simpler version here for clarity. 2 its rank is equal to its smaller dimension. then p θ (x|y) = p θ (x|y) =⇒ θ ∼ g w θ where ∼ g w is defined in (4). Finally, if all assumptions 1-4 hold, then the matrices A in (3) and B in (4) have full rank equal to d z .
This initial form of identifiability requires very little assumptions on the feature extractors f and g. In fact, the assumptions of this theorem are easy to satisfy in practice, and we will see how in the following section. Note that this result also covers the case where the number of features is larger than the number of observed variables. As far as we know, this is the first identifiability result that extends to overcomplete 3 representations in the nonlinear setting.

IDENTIFIABLE ARCHITECTURES
We explore in this section how feasible the assumptions of Theorem 1 are. We will only discuss assumptions 1 and 2 since assumptions 3 and 4 are very similar.
We suppose that each of the networks f and g are parameterized as multi-layer perceptrons (MLP). The case of other neural network architectures like convolutional networks is left for future work. More specifically, consider an MLP with L layers, where each layer consists of a linear mapping with weight matrix W l ∈ R d l ×d l−1 and bias b l , followed by an activation function. Consider the following architecture: (a.) All activation functions are LeakyReLUs.
(b.) All weight matrices W l are full rank.
(c.) The row dimension of the weight matrices are either monotonically increasing or decreasing: This architecture is sufficient to satisfy the assumptions of Theorem 1, as we demonstrate in detail below.
Assumption 1 This assumption requires that the mapping be differentiable and have a full rank Jacobian. The following proposition explains how this can be achieved by an MLP in practice, the proof of which is presented in Appendix A.2.
Consider an MLP f whose architecture satisfies assumptions (a.), (b.) and (c.), then f has a full rank Jacobian almost everywhere. If in addition, In other words, f satisfies Assumption 1.
In practice, and in addition to this, we need to pay attention to the way the weights change due to the optimization algorithm. One way to ensure that the MLP remains full rank is by performing spectral normalization on the weights, so that we ensure the lowest singular value is positive, or by encouraging orthogonality of the weight matrices via a penalty in the objective function for instance.
Assumption 2 This assumption requires that the conditioning feature extractor, g, has an image that is rich enough. Intuitively, this lifts the amount of flexibility the main feature extractor f would need to have if g was otherwise very simple. It implies that the search for f will be naturally restricted to a smaller space, for which we can prove identifiability. This is formalized in the proposition below, the proof of which is provided in Appendix A.3.
Linear feature extractors The particular case of linear feature extractors is quite interesting. If d z ≤ d y and the feature extractor g satisfies the assumptions of Proposition 1, then assumption 2 is trivially satisfied.
On the other hand, if d z > d y , then assumption 2 can't hold when the network is linear. This signals that it is important to use deep nonlinear networks to parameterize the feature extractors, at least in the overcomplete case.

UNIVERSAL APPROXIMATION CAPABILITY
With such a potentially overcomplete network, we can further achieve universal approximation of the data distribution. It might initially seem that this is an impossible endeavor given the somehow restricted form of the energy function (1). It turns out that, if we also consider the dimension d z of f and g as an additional parameter, then we can find a tuple (θ, d z ) ∈ Θ × N such that p θ,dz (x|y) ≈ p(x|y) arbitrarily well. This is summarized by the following theorem, which we prove in Appendix A.6.
Theorem 2. Let p(x|y) be a conditional probability density. Assume that X and Y are compact Hausdorff spaces, and that p(x|y) > 0 almost surely ∀(x, y) ∈ X × Y. Then for each ε > 0, there exists (θ, d z ) ∈ Θ × N, where d z is the dimension of the feature extractor, such that This means that our model is capable of approximating any conditional distribution that is positive on its compact support arbitrarily well. In practice, the optimal dimension d z of the feature extractors can be estimated using cross-validation for instance. It is possible that to achieve a near perfect approximation, we require a value of d z that is larger than the dimension of the input. This is why it is crucial that our identifiability result from Theorem 1 covers the overcomplete case as well, and highlights the importance of our contribution in comparison to previous identifiable deep models.

EXTENSIONS FOR STRONGER IDENTIFIABILITY
The identifiability result given in section 2.2 stipulates that two different parameters that define the same conditional density map to feature extractors that are linearly related. We wish to strengthen this result akin to what was done by Khemakhem et al. (2020, Theorems 2,3). This however required a slightly modified model architecture.
We will discuss and prove results for the feature extractor f only, but as we saw in the previous section, we can have the same results for g by a simple transposition of the assumptions.
The stronger identifiability results are summarized by the following equivalence relation.
Definition 2. Let ∼ s be the equivalence relation on Θ defined as: where σ is a permutation of [[1, n]], a i is a non-zero scalar and b i is a scalar.
We propose two different extensions to our energy model's architecture which will both allow for this stronger form of identifiability, ∼ s . The first extension is for non-negative feature extractors. The second extension is based on augmenting the output of the feature extractor by additional activation functions, and use the concatenated output to define the energy function. Importantly, we will suppose that the output dimension d z is smaller than the input dimension d x .

NON-NEGATIVE FEATURE EXTRACTOR
Suppose in this section that the feature extractor f (x) ∈ R dz is non-negative for all x ∈ X ⊂ R dx , and that d z ≤ d x . The non-negativity will induce constraints on the matrix A defining the equivalence relation ∼ f w : loosely speaking, if A induces a rotation in space, then it will violate the non-negativity constraint, since the only rotation that maps the positive portion of the plan to itself is the identity. Thus A can only be a scaled permutation. This is made rigorous by the following theorem: Theorem 3. Assume the assumptions of Theorem 1 hold. Further assume that d z ≤ d x , and that for any choice of parameter θ, the feature extractor f θ is surjective onto R dx + , and each component f i,θ is non-negative and reaches zero. Then In practice, we can either square or compute the absolute value of the features f (x), or apply a ReLU to the final layer of the neural network parameterizing f , which effectively produces non-negative features. The non-negativity assumption is not really required, and could be replaced by existence of any lower bound, but it is used here for simplicity. We can also just assume that the minimum is reached in the limit, and the conclusion would still hold. A more detailed discussion of this, along with a proof of Theorem 3, are provided in Appendix A.4.

AUGMENTED FEATURE EXTRACTOR
In this section, we consider an uncostrained feature extractor f (x) ∈ R dz , which we augment by its square, effectively resulting in the 2d z -dimensional feature extractorf which we combine with a 2d z -dimensional feature extractorg(y) ∈ R 2dz for the conditioning variable y, to define an augmented energy functioñ The energy function is slightly more constrained, but this results in much stronger identifiability results, as we will see below. The advantage of this approach is that it doesn't require the feature extractors to be positive (or, more generally, to have a strict lower or upper bound). The drawback, however, is that the effective size of the feature extractor is 2d z , and this will be reflected in the assumptions of the following theorem, which we prove in Appendix A.5.
Theorem 4. Assume that d z ≤ d x , and that for any choice of parameter θ : 1. The feature extractor f θ is differentiable and surjective, and its Jacobian J f θ is full rank.
Theorems 3 and 4 are fundamental as they prove very strong identifiability results for a conditional deep energybased model. As far as we know, our results require the least amount of assumptions in recent theoretical work for functional identifiability of deep learning models (Khemakhem et al., 2020;Sorrenson et al., 2020). Most importantly, we do not make any assumption on the distribution of the learned features.
This, however, doesn't mean that our model can't be used to estimate latent variables in a latent variable models. The first application that comes to mind is nonlinear ICA.
In the next section, we develop a novel generative model that generalizes ICA to non-independent latent variables, and show how we can use our energy model to estimate them.

INDEPENDENTLY MODULATED COMPONENT ANALYSIS
In this section, we show how the energy-based model above relates to a generative, latent variable model. Until recently, results relating to identifiability of latent variable models were constrained to linear models (e.g., as in linear ICA), as it was acknowledged that the flexibility of nonlinear mappings could yield arbitrary latent variables (Hyvärinen and Pajunen, 1999). However, it is now understood that nonlinear latent variable models may also be identifiable provided some additional auxiliary variables (Khemakhem et al., 2020;Hyvärinen et al., 2019). The purpose of this auxiliary variable is to introduce additional constraints on the distribution of the latent variables, which are typically required to be conditionally independent given the auxiliary variable. More precisely, the latent variable z has a density of the form p(z|y) = i p i (z i |y).
However, the need for (conditional) independence in order to obtain identifiability results may sometimes be seen as a limitation. This is particularly the case in the context of learning disentangled representations, where independence is seen as too severe a restriction. We would like to relax the assumption of independence while maintaining identifiability. This was achieved before in the linear case (Monti and Hyvärinen, 2018;Hyvärinen and Hurri, 2004), and we next show how it may be achieved in the nonlinear setting.

MODEL DEFINITION
Assume we observe a random variable x ∈ R dx as a result of a nonlinear transformation h of a latent variable z ∈ R dz whose distribution is conditioned on an auxiliary variable y ∈ R dy that is also observed: We begin by considering the scenario where d x = d z = d and consider dimensionality reduction as a special case.
The main modelling assumption we make on the latent variable is that its density has the form: The exponential term factorizes across components: this means that the sufficient statistic is composed of d functions that are each a function of only one component z i of latent variable z. On the other hand, the base measure µ(z) is unconstrained, and not necessarily factorial.
Equations (7) and (8) together define a nonparametric model with parameters (h, T, λ, µ). For the special case µ(z) = i µ i (z i ), the distribution of z factorizes across dimensions, and the components z i are independent. Then the generative model gives rise to a nonlinear ICA model, and it was studied to a great depth by Khemakhem et al. (2020). In particular, the authors identified sufficient conditions for a strong identifiability of the model that allows for recovery of the latent variables up to permutation and nonlinear scaling.
We propose to generalise this model by allowing for an arbitrary base measure µ(z), i.e. the components of the latent variable are no longer independent, as µ doesn't necessarily factorise across dimensions. However, it is crucial that the components of the latent variable are independently modulated given the auxiliary variable y, and that through the term exp( i T i (z i ) T λ i (y)).
We call this new framework Independently Modulated Component Analysis (IMCA). In the next section we show that the strong identifiability guarantees developed for nonlinear ICA can be extended to IMCA, yielding a more general and more flexible principled framework for representation learning and disentanglement. A more thorough exposition of the theoretical properties of the IMCA framework are presented in Appendix C.

THEORETICAL ANALYSIS
We can heuristically understand the connection between the ICE-BeeM model and IMCA as follows. Under the IMCA model, the likelihood can be computed as Now, assume our energy-based model has learned the log probability density function exactly. If we equate, purely heuristically, the individual terms in (8) and (2), we have which means the f i may recover the latent variables up to the nonlinear functions which are the sufficient statistics T i , but we emphasize this is a heuristic result.
Next, we will rigorously prove the ability of the energybased model to estimate the IMCA model, as shown by the following theorems, under both of the weak and strong identifiability assumptions. To be able to prove strong identifiability results, we augment our feature extractor by output activation functions ). This is inspired by Section 3.2. We now state the main result of this section. All proofs for this section are presented in Appendix B. Theorem 5. Assume: (i) The observed data follows the exponential IMCA model of equations (7)-(8).
(ii) The mixing function h is a D 2 -diffeomorphism 4 .
(iii) The sufficient statistics T i are twice differentiable, and the functions T ij ∈ T i are linearly independent on any subset of X of measure greater than zero 5 . Furthermore, they all satisfy dim(T l ) ≥ 2, ∀l or dim(T l ) = 1 and T l is non-monotonic ∀l.
(iv) There exist k + 1 distinct points y 0 , . . . , y k such that the matrix We use a consistent estimator to fit the model (2) to the conditional density p(x|y), where we assume the feature extractor f (x) to be a D 2 -diffeomorphism and d-dimensional, and the pointwise nonlinearitiy H to be differentiable and k-dimensional, and the dimensions of its vector-valued components H l to be chosen from (dim(T 1 ), . . . , dim(T d )) without replacement.
Then, in the limit of infinite data, we can recover the latent variables up to an block permutation linear transformation and point-wise nonlinearities, that is where γ is a permutation of [[1, d]] such that dim(H i ) = dim(T γ(i) ) and A i a square invertible matrix.
Dimensionality reduction In practice, it is a natural desire to have the feature extractor reduce the dimension of the data, as it is usually very large. This has been achieved in nonlinear ICA before (Khemakhem et al., 2020;Hyvärinen and Morioka, 2016). It turns out that we can also incorporate dimensionality reduction in IMCA and its estimation by ICE-BeeM, under some assumptions.
Suppose that only n out of d components of the latent variable are modulated by the auxiliary variable y. In other words, we assume that the parameters λ n+1:d (y) are constant, and we can write its density as The term e d i=n+1 Ti(zi) T λi is simply absorbed into µ(z). In this case, the feature extractor f (x) ∈ R n is still capable of recovering the n modulated components. This is summarised by the following theorem.
Theorem 6. Assume the assumptions of Theorem 5 hold. Further assume: (vi) Only n < d components of the latent variable are modulated, and its density has the form (10).
(vii) The feature extractor f has the form f (x) = (f 1 (x), f 2 (x)) where f 1 (x) ∈ R n , and the auxiliary feature extractor g has the form g(y) = (g 1 (y), g 2 ) where g 1 (y) ∈ R n and g 2 is constant.
Then f 1 recovers the n modulated latent components as per equation (9).

TRANSFER LEARNING BY EBM
As a practical application of our framework where identifiability is important, we consider metalearning, in particular multi-task and transfer learning. Assume we have N datasets, which could be, e.g., different subjects in biomedical settings, or different image datasets. This fits well with our framework, where y = 1, . . . , N is now the index of the dataset, or "task".
The key question in such a setting is how we can leverage all the observations to better model each single dataset, and especially transfer knowledge of existing model to a new dataset. To this end, we propose an intuitively appealing approach, where we approximate the log-pdf in y-th dataset p(x; y) by a linear combination of a learned "basis" functions f i,θ as where the g i (y) are scalar parameters as a function of y. This linear approximation is nothing else than a special case of ICE-BeeM, but here, we interpret such an approximation as a linear approximation in log-pdf space.
In fact, what we are doing is a kind of PCA in the set of probability distributions p(x; y). Such "probability space" PCA allows the models for the different datasets to learn from each other, as in the classical idea of denoising by projection onto the PCA subspace.
In transfer learning, we observe a new dataset, with distribution p(x; y new ) for y new = N + 1. Based on our decomposition, we approximate p(x; y new ) as in (11). This leads to a drastic simplification since we only need to estimate the k scalar parameters g i (y new ) for the new dataset, which is possible by any method for EBM estimation (as discussed in Section 6). In particular, it is very easy by score matching because this is an exponential family for fixed f i (Hyvärinen, 2007). In addition, the coefficients are likely to be sparse as well, which provides an additional penalty.
An important point here is that transfer learning is easier if we can uniquely identify the components f i in this decomposition. Reducing the transfer learning to estimation of the g i (y new ) clearly requires that we have estimated the true f i up to a linear transformation, which is the weaker form of identifiability in Theorem 1. Moreover, using a sparsity penalty is only meaningful if we have the true f i without any linear mixing, which requires the stronger identifiability in Theorems 3 and 4. In that case, rather surprisingly, the components of the factor-analytic model are identifiable due to constraints on the f i which are like factor loadings in classical terminology. Thus, we have an identifiable version of factor analysis in the space of probability densities, which may have wide applicability beyond transfer learning.
We mention in passing that another application where identifiability is known to be crucial is causal discovery (Monti et al., 2019). In concurrent work by Teshima et al. (2020), a related approach is proposed to domain adaptation.
It is important to note that the identifiability results presented above apply to conditional EBMs in general. As such, we may employ any of the wide variety of methods which have been proposed for the estimation of unnormalized EBMs. In this work we used two different options with good results for both: flow contrastive estimation (Gao et al., 2019) and denoising score matching (Vincent, 2011). Both methods can also be extended to the conditional case in a straightforward fashion.
Flow-contrastive estimation (FCE) can be seen as an extension of noise-contrastive estimation (Gutmann and Hyvärinen, 2012, NCE), which seeks to learn unnormalized EBMs by solving a surrogate classification task. The proposed classification task seeks to discriminate between the true data and some synthetic noise data based on the log-odds ratio of the EBM and the noise distribution. However, a limitation of NCE is the need to specify a noise distribution which can be sampled from and whose log-density can be evaluated pointwise but which also shares some of the empirical properties of the observed data. To address this concern Gao et al. (2019) propose to employ a flow model as the contrast noise distribution. FCE seeks to simultaneously learn both an unnormalized EBM as well as a flow model for the contrast noise in an alternating fashion. We naturally get a conditional version for FCE by learning a conditional EBM (Gao et al., 2019, eq. 12).
Score matching is another well-known method for learning unnormalized models (Hyvärinen, 2005). However, its computational implementation in deep networks is problematic, which is why Vincent (2011) proposed a stochastic approximation which can be interpreted as denoising the data, and which works efficiently in deep networks (Saremi et al., 2018;Song and Ermon, 2019). We extend the denoising score matching (DSM) objective to the conditional setting in a straightforward fashion, and the objective we minimize is where the expectation is computed with respect to observations x and y as well as noisy samples x ∼ N (x, σ 2 I).

EXPERIMENTS
We explore the importance of identifiability and the applicability of ICE-BeeM in a series of experiments. First, we quantify the performance of ICE-BeeM on nonlinear unmixing tasks under ICA and IMCA models. Second, we present a transfer learning example on MNIST.

SIMULATIONS
Nonlinear ICA experiments In order to demonstrate the capabilities of the proposed conditional energy based framework we run a series of simulations comparing ICE-BeeM to previous nonlinear ICA methods such as iVAE (Khemakhem et al., 2020) and TCL (Hyvärinen and Morioka, 2016). We generate 5-dimensional synthetic datasets following the nonlinear ICA model which is a special case of equation (7) where the base measure, µ(z), is factorial. As such, latent variables are conditionally independent given segment labels. The conditioning variable y is defined to be the segment index. Following Hyvärinen and Morioka (2016), the z are generated according to isotropic Gaussian distributions with distinct variances λ(y) determined by the segment index. A randomly initialized neural network with varying number of layers, L ∈ {2, 4}, was employed to generate the nonlinear mixing function h. Leaky ReLU was employed as the activation function in order to ensure the network was invertible. Throughout experiments we fix the number of segment to be 8 and vary the number of observations per segment.
For each nonlinear ICA method, an MLP of depth L was employed to recover latent variables. In the case of ICE-BeeM, f θ was parameterized by an MLP with L layers whilst g θ was defined to be a distinct vector per segment. FCE was employed to estimate network parameters. To evaluate the performance of the method, we compute the mean correlation coefficient (MCC) between the true latent variables and the recovered latents estimated by all three methods. Results are provided in Figure [1a], where we note that ICE-BeeM performs competitively with respect to both iVAE and TCL. We note that as the depth of the mixing network, L, increases the performance of all methods decreases.
Nonlinear IMCA experiments Second, we perform the same experiment but on data generated from an IMCA model where the base measure µ(z) is not factorial. We set it to be a Gaussian term with a fixed but non-diagonal covariance matrix. More specifically, we randomly generate an invertible and symmetric matrix Σ 0 ∈ R d×d , such that µ(z) ∝ e −0.5z T Σ −1 0 z . As before, we define λ(y) to be the distinct conditional variances. The precision matrix of each segment is now equal to Σ(y) −1 = Σ −1 0 + diag(λ(y)) −1 , meaning the latent variables are no longer conditionally independent. As before, FCE was employed to estimate parameters of ICE-BeeMmodels. Results are provided in Figure [1b] where ICE-BeeM outperforms alternative nonlinear ICA methods. This is because such other methods implicitly assume latent variables are conditionally independent and are therefore misspecified.

TRANSFER LEARNING ON MNIST
We also present an application of ICE-BeeM to transfer learning as described in Section 5. An ICE-BeeM model was trained on digits 0-7 using the conditional DSM objective. We follow Song and Ermon (2019) and employ a U-net architecture for f θ and allow g θ to be parameterized by a vector for each segment. After training, we fix f θ and learn new g θ (y new ) for the unseen classes. We compare against a baseline where both f θ and g θ are estimated using new classes only (i.e., there is no transfer learning).
Results are presented in Figure [1c] where we vary the sample size over new classes and report the conditional DSM objective. We note that for small sample sizes, the use of a pretrained f θ network improves performance, demonstrating effective transfer learning. The top panel of Figure [2] shows unseen samples (taken across all possible classes) which are assigned high confidence of belonging to the "new" class 8 after transfer learning, indicating that the ICE-BeeM model has learnt a reasonable distribution over unseen classes. By comparison the case where no transfer learning is employed (bottom panel), incorrectly assigns high confidences to other digits.

CONCLUSION
Unsupervised feature learning is one of the fundamental challenges in machine learning. Despite empirical success, a limitation of current methods is their lack of strong theoretical grounding. Arguably, the gold standard in terms of theory is identifiability, which implies that latent variables can be uniquely recovered from data.
We proposed a new identifiable conditional energy-based deep model, or ICE-BeeM for short, for unsupervised representation learning. This is probably the first energybased model to benefit from rigorous identifiability results. Crucially, the model benefits from the tremendous flexibility and generality of EBMs. We even prove a universal approximation capability for the model.
We further prove a fundamental connection between EBMs and latent variable models, showing that ICE-BeeM is able to estimate nonlinear ICA, as a special case. In fact, it can even estimate a generalized version where the components do not need to be independent: they only need to be independently modulated by another variable such as a time index, history or noisy labels.
Identifiability is fundamental for meaningful and principled disentanglement; it is necessary to make any interpretation of the features meaningful; it is also crucial in such applications as causal discovery and transfer learning. The present results go further than any identifiability results hitherto and extend them to the EBM framework. We believe this paves the way for many new applications of EBMs, by giving them a theoretically sound basis. We divide the Appendix into 4 main sections: in section A, we prove the identifiability of ICE-BeeM and its universal approximation capability; in section B, we show how ICE-BeeM estimates IMCA; in section C, we provide a thorough theoretical analysis of the IMCA framework and draw parallels to the identifiability results in nonlinear ICA; finally, in section D, we discuss the estimation algorithms we used with ICE-BeeM and how they can be extended to the conditional setting.

A IDENTIFIABILITY OF THE CONDITIONAL ENERGY-BASED MODEL
Recall the form of our conditional energy model We present in this section the proofs for the different forms of identifiability that is guaranteed for the feature extractors f and g. We will focus on the proofs for the feature extractor f , as the proofs for the feature extractor g are very similar. For the rest of the Appendix, we will denote by d = d x , m = d y and n = d z .

A.1 PROOF OF THEOREM 1
We start by proving the main theoretical result of this paper, which applies to all dimensions of the feature extractor.
Theorem 1 (Identifiable conditional EBMs). Assume: 1. The feature extractor f is differentiable, and its Jacobian J f is full rank.
A is a (d z × d z )-matrix of rank at least min(d z , d x ).
If, instead or in addition, we assume that: 3. The feature extractor g is differentiable, and its Jacobian J g is full rank.
4. There exist n + 1 points x 0 , . . . , x n such that the matrix B is a (d z × d z )-matrix of rank at least min(d z , d x ).
Finally, if d z ≥ max(d x , d y ) and all assumptions 1-4 hold, then the matrices A and B have full rank (equal to d z ).
Proof. We will only prove this theorem for the feature extractor f . The proof for g is very similar. Suppose assumptions 1 and 2 hold.
Consider two parameters θ andθ such that p θ (x|y) = pθ(x|y) Then, by applying the logarithm to both sides, we get: Consider the points y 0 , . . . , y n provided by assumption 2 for g θ . We plug each of these points in (18) to obtain n + 1 such equations. We subtract the first equation for y 0 from the remaining n equations, and write the resulting eauqtions in matrix form: where R = (. . . , g θ (y l ) − g θ (y 0 ), . . . ),R = (. . . , gθ(y l ) − gθ(y 0 ), . . . ), and b = (. . . , log Z(y l ;θ) Z(y l ;θ) − log Z(y 0 ;θ) Z(y 0 ;θ) , . . . ). Since R is invertible (by assumption 2), we multiply by its inverse from the left to get: where A = R −1R and c = R −1 b. Now since f θ is differentiable and its Jacobian is full rank (assumption 1), by differentiating the last equation we deduce that rank(A) ≥ min(n, d), which in turn proves that θ ∼ f wθ . Finally, suppose that in addition, assumptions 4 holds. Then there exists x 0 , . . . x n such that Q := (. . . , f θ (x i ) − f θ (x 0 ), . . . ). Plugging these n + 1 points into eauqtion (20), and subtracting the first equation for x 0 from the remaining n equations, we get Since Q is an n × n invertible matrix, we conclude that A is also invertible, which concludes the proof.
Proof under weaker assumptions Assumption 2 of full rank Jacobian can be weakened without changing the conclusion of Theorem 1. In fact, this assumption is only used right after equation (20) to prove that the matrix A has a rank that is at least equal to min(n, d). Suppose instead that 1.' There exists a point x 0 ∈ R d where the Jacobian J f θ of f θ exists and is invertible Then by computing the differential of equation (20) at x 0 (assuming that J fθ (x 0 ) exists), we can make the same conclusion on the rank of A.
A.2 PROOF OF PROPOSITION 1 Proposition 1. Consider an MLP with L layers, where each layer consists of a linear mapping with weight matrix W l ∈ R d l ×d l−1 and bias b l , followed by an activation function. Assume a. All activation functions are LeakyReLUs.
b. All weight matrices W l are full rank.
c. The row dimension of the weight matrices are either monotonically increasing or decreasing: Then the MLP has a full rank Jacobian almost everywhere. If in addition, d L ≤ d 0 , then the MLP is surjective.
Proof. Denote by x the input to the MLP, and by x l the output of layer l, x l = W l x l−1 + b l (23) with h in equation (24) is applied to each element of its input, and α ∈ (0, 1).
Denote by v l ∈ R d l the vector whose elements are which is undefined if x l k = 0, and by V l = diag(v l ). Note that V l is a function of its input, and thus of x, but we keep this implicit for simplicity. Using these notations, and the fact that h is piece-wise linear, we can write, x where v l k is the k-th row of V l . For each x / ∈ N , we have that V l is full rank, and, using Lemma 2, V l is also a full rank matrix.
While it is true that b l k and v l k are functions of x, yet they only take a finite number of values. Thus, the set is included in the union over all the values taken by b j k and v j k up to layer l. For each of these values, the set becomes a dot product between a row of V j which is independent of the input x, and is nonzero because V j is full rank; such set has measure zero in R d . Thus, N is included in a finite union of sets of measure zero, which implies that it also has measure zero.
Now, for all x / ∈ N , ∂x L ∂x exists, and can be computed using the chain rule: which is full rank. Thus, the MLP has a full rank Jacobian almost everywhere.
The surjectivity is easy to prove since h is surjective and so is x l as a function of x l−1 if d l−1 ≥ d l and rank(W l ) = d l . Proof. The proof in the case m ≥ n ≥ p can be found in (Arbel et al., 2018, Lemma 10), but we provide a proof here for completeness, and for the other case m ≤ n ≤ p.
Let R n * := R n \ {0}, and λ min (A) the smallest eigenvalue of A. Recall that for a matrix A ∈ R n×m , with m ≥ n, Thus, if the null space of N is non trivial, then σ min (N) = 0, and the inequality is satisfied. Otherwise, we have Similarly, if the null space of M T is non trivial, then σ min (M T ) = σ min (M) = 0, and the inequality holds. Otherwise, which concludes the proof.
Lemma 2. Consider a finite sequence of matrices (M i ) 1≤i≤p , with M i ∈ R ni−1×ni . If M i is full rank for all i ∈ [[1, p]], and either n 0 ≤ n 1 ≤ . . . ≤ n p or n 0 ≥ n 1 ≥ . . . ≥ n p , then the product M 1 M 2 . . . M p is also full rank.
Proof. If two matrices M 1 and M 2 with ordered dimensions are full rank, then σ min (M 1 ) > 0 and σ min (M 2 ) > 0. According to Lemma 1, this implies that σ min (M 1 M 2 ) > 0, and that M 1 M 2 is full rank. The proof for p ≥ 3 is done by induction on p.

A.3 PROOF OF PROPOSITION 2
Proposition 2. Consider an MLP g with L layers, where each layer consists of a linear mapping with weight matrix W l ∈ R d l ×d l−1 and bias b l , followed by an activation function. Assume a. All activation functions are LeakyReLUs.
b. All weight matrices W l are full rank.
c. The row dimension of the weight matrices are either monotonically increasing or decreasing: d. The intercepts b l,i ∈ b l are distinct.
Proof. Let y 0 be an arbitrary point in R d0 . Without loss of generality, suppose that g(y 0 ) = 0. This is because y → g(y) − g(y 0 ) is still an MLP that satisfies all the assumptions above. If for any choice of points y 1 to y d L , the matrix R defined above isn't invertible, then this means that g(R d0 ) is necessarily included in a subspace of R d L of dimension at most d L − 1. In other words, this would imply that the functions g 1 , . . . , g d L are not linearly independent. We will consider two cases, and prove that this is not possible.
First, if d L ≤ d 0 , then by the conclusion of Proposition 1, g is surjective. Thus, its image can't be included in a subspace of R d L of dimension at most d L − 1.
Second, if d L > d 0 , then using Lemmas 5 and 6, we can prove that g 1 , . . . , g d L are linearly independent.
Thus, if such points do not exist, we always run into a contradiction, and we can conclude thatthere exist d L + 1 points y 0 , . . . , y d L such that the matrix R = g(y 1 ) − g(y 0 ), . . . , g(y d L ) − g(y 0 ) is invertible.
Proof under weaker conditions Note that the proof argument used for the overcomplete case can be used for the undercomplete as well. This same argument can be proved for ReLU as the nonlinearity instead of LeakyReLU. We chose to give the proof for, and suggest to use the latter because it is needed for Proposition 1.
Lemma 3. Let A be an n × n invertible matrix. Denote by a n the n-th row of A. Then the matrix B ∈ R n+1,n+1 such that is invertible for any choice of γ 1 , . . . , γ n−1 , and for λ = 1.
Proof. Denote by b i the i-th row of B. Let α 1 , . . . , α n+1 such that Then in particular, by looking at the first n lines of this vectorial equation, we have that n−1 i=1 α i a i + (α n + α n+1 )a n = 0. Since A is invertible, its rows are linearly independent, and thus α n = −α n+1 and α i = 0, ∀i < n. Plugging this back into equation (33), and looking closely at the last equation, we have that (1 − λ)α n = 0, and we conclude that α n+1 = α n = 0 (because λ = 1), and that B is invertible.

Lemma 4. Consider n affine functions
Then there exist n non-empty regions H 1 , . . . , H n of R n defined by the signs of the functions f i (for instance, H = {x ∈ R n |∀i, f i (x) > 0}) such that the matrix S n ∈ R n×n defined as S n i,j = sign x∈Hi (f j (x)) is invertible.
Proof. We will prove this Lemma by induction on n the number of functions f i . Denote by The sign of f i changes if we cross the hyperplan V i .
First, suppose that n = 2. If a 1 ∝ a 2 , then the hyperplans V 1 and V 2 are not parallel, and they divide R n into 4 regions. This implies that the regions are not empty. The same reasoning can be used if a 1 ∝ a 2 , but b 1 = b 2 which concludes this step. Second, suppose that there exists n regions H 1 , . . . H n such that the the matrix S n is invertible. Consider the affine function f n+1 = a T n+1 x + b n+1 . The hyperplan V n+1 = {x ∈ R d |f n+1 (x) = 0} intersects at least one of the regions H 1 , . . . H n . This is because for all i ≤ n, either a n+1 ∝ a i , or b n+1 = b i . Suppose without loss of generality that this region is H n . Denote byH n = {x ∈ R n |x ∈ H n , f n+1 (x) < 0} ⊂ H n . Now consider the matrixS n such that S n n,j = sign x∈Hn (f j (x)) andS n i,j = S n i,j . BecauseH n ⊂ H n , we have that sign x∈Hn (f j (x)) = sign x∈Hn (f j (x)) and thusS n = S n , which implies thatS n is also invertible. Now define H n+1 = {x ∈ R n |x ∈ H n , f n+1 (x) > 0} ⊂ H n . Again, the inclusion implies that sign x∈Hn (f j (x)) = sign x∈Hn (f j (x)). Finally, consider the regions H 1 , . . . , H n−1 ,H n , H n+1 , and the matrix S n+1 defined on those regions. Then where u i = sign x∈Hi f n+1 (x) and s n n is the n-th line of S n . According to Lemma 3, S n+1 is invertible, which achieves the proof.
Lemma 5. Let h denote a LeakyReLU activation function with slope λ ∈ (0, 1). Consider n > d scalar functions such that the matrix A ∈ R n×d whose rows are the a i is full column rank, and such that b 1 , . . . , b n are distinct. Then the functions g 1 , . . . , g n are linearly independent. Proof.
First, since A is full column rank and b 1 , . . . , b n are distinct, the assumptions of Lemma 4 are met for the function f 1 , . . . , f n , and we conclude that there exists n regions H 1 , . . . , H n such that S n = sign x∈Hi (f j (x)) i,j is invertible.
Define the matrixS where we replace all entries of S n by λ if they are equal to −1. ThenS is invertible (in fact, to see this, consider the proof of the previous lemma with the slightly unconventional choice of sign function sign(x) = λ if x < 0). Now consider α 1 , . . . , α n such that Let k ∈ [[1, n]], and evaluate this equation at x ∈ H k . After taking the gradient with respect to x, we get Denote bys k the k-th line of the matrixS, and define e l = (α 1 a 1,l , . . . , α n a n,l ) ∈ R n . We can write the l-th line of equation (36) as:s T k e l = 0 Collating these equations for a fixed l and k ∈ [[1, n]], we get which implies that e l = 0 because S is invertible. In particular, e l = 0 for all l ∈ [[1, d]] implies that A T α = 0. Since A is full column rank by assumption, we deduce that α = 0, and we conclude that the functions g 1 , . . . , g n are linearly independent.
Lemma 6. Let h denote a LeakyReLU activation function with slope λ ∈ (0, 1), and t = (t 1 , . . . , t n ) be a vectorvalued function defined on R d , such that t 1 , . . . , t n are linearly independent piece-wise affine functions. Consider m > n scalar functions g i : x ∈ R d → h(a T i t(x) + b i ), such that the matrix A ∈ R n×d whose rows are the a i is full column rank, and such that b 1 , . . . , b n are distinct. Then the functions g 1 , . . . , g m are linearly independent.
Proof. In Lemma 5, we've shown that applying LeakyReLU to n affine functions makes them linearly independent. This is only possible if their slopes form a full column rank matrix, and their intercepts are all indistinct. What this condition ensures is that their signs change at sufficiently varied locations and not all simultaneously so that we the output of LeakyReLU is not everywhere the same. If now we aplly LeakyReLU to m piece-wise stationary functions that satisfy the assumptions above, then it is more likely to find these regions in space which would make them linearly independent. Then, we can follow the reasoning of Lemmas 5 and 4, making the small necessary adjustment to account for piece-wise affine functions, which would prove this Lemma.

A.4 PROOF OF THEOREM 3
We will prove here a more general version where we assume that each component f i of the feature extractor f has a global minimum that is reached, instead of being necessarily non-negative.
Theorem 3. Assume the assumptions of Theorem 1 hold. Further assume that n ≤ d, and that each f i has a global minimum that is reached at least in the limit, and the feature extractor f = (f 1 , . . . , f n ) is surjective onto the set that is defined by the lower bounds of the f i . Then where ∼ s is defined as follows: where σ is a permutation of [[1, n]], a i is a non zero scalar and b i is a scalar.
Proof. Consider two different parameters θ andθ such that: To simplify notations, denote by f = f θ andf = fθ. We start the proof from the conclusion of Theorem 1, since its assumptions hold: where A is an invertible n × n matrix and c a constant vector. Without loss of generality, we can suppose that f i has an infimum equal to zero, simply by subtracting inf f i , and including in c, and similarly forf . We will also suppose that the infimums are reached, as the next argument would hold if we change exact minima by limits.
Now since f ≥ 0 and is surjective, then there exists . We know thatf ≥ 0 and is surjective, and so h is also surjective, and its image includes R n + . Let I = (e 1 , . . . , e n ) be the matrix of canonical basis vectors, or positive scalar multiples of the canonical basis vectors e i . These must be mapped to the non-negative quadrant, so AI must be non-negative, which implies that A must be non-negative.
Denote by B = A −1 . B is also non-negtaive for the same reasons described above. Denote the rows of A by a i and the columns of B by b j . We have by definition of inverse: where if i = j then δ ij = 1, else δ ij = 0. Now, assume there is a row a k which has at least two non-zero entries. By the property above, d − 1 of the vectors b j must have zero dot-product with that vector. By non-negativity of B and A, those d − 1 vectors must have zeros in the at least two indices corresponding to the non-zeros of a k . But that means they can only span a d − 2-dimensional subspace, and all the b j together can only span a d − 1-dimensional subspace. This is in contradiction of the invertibility of B. Thus, each a i can have only one non-zero entry, which, together with the invertibility of A, proves it is a scaled permutation matrix.

A.5 PROOF OF THEOREM 4
Theorem 4. Assume that n ≤ d, and that: 1. The feature extractor f is differentiable and surjective, and its Jacobian J f is full rank.
2. There exist 2n + 1 points y 0 , . . . , y 2n such that the matrix R = g(y 1 ) −g(y 0 ), . . . ,g(y 2n ) −g(y 0 ) of size 2n × 2n is invertible. Then Proof. Similarly to the proof of Theorem 3, we pass the features f i through the nonlinear function H i (f i ) = (f i , f 2 i ) which produces the augmented featuresf introduced in section 3.2.
Consider two different parameters θ andθ such that: Since we have similar assumptions to Theorem 1, we will skip the first part of the proof and make the same conclusion, where the equivalence up to linear transformation here applies to H(f θ ) and H (fθ): where A is a 2n × 2n matrix of rank at least n because J f and J H are full rank (A is not necessarily invertible yet, but this will be proven later) and c a constant vector. By replacing H by its expression, we get: where each A (i) is an n × n matrix, and c = (α, β). To simplify notations, denote by h = fθ. We will also drop reference to θ andθ. The first n lines in the previous equation are: and the last n lines are: Fix an index i in equations (47) and (48). To alleviate notations and reduce the number of subscripts and superscripts, we introduce a j = A (1) ij , α = α i and β = β i . This proof is done in 5 steps. Note that the surjectivity assumption is key for the rest of the proof, and it requires that we set the dimension of the feature extractor to be lower than the dimension of the observations. By equating equations (48) and (47) after squaring, we get, using our new notations: Step 1 First, since h is surjective, there exists a point where it is equal to zero. Evaluating equation (49) at this point shows that β = α 2 .
Step 2 Second, the left hand side of equation (49) has terms raised to the power 4. These terms grow to infinity much faster than the rest of the terms of the rhs and the lhs. It is thus equal to zero. More rigorously, consider the vectors e l (y) = (0, . . . , y, . . . , 0) ∈ R n where the only non zero entry is y at the l-th position. Each of these vectors has a pre-image by h (since it is surjective), which we denote by x l (y). By evaluating equation (49) at each of these points, we get (a l y + b l y 2 + α) 2 = c l y + d l y 2 + β Devide both sides of this equation by y 4 , then take the limit y → ∞. The right hand side will converge to 0, while the left hand side will converge to b l , which shows that b l = 0. By doing this process for all l ∈ [[1, n]], we can show that b = 0.
Step 3 So far, we've shown that (49) becomes, after expanding the square in the lhs, and writing j a j h j (x) = a T h(x): Let's again consider the vectors e l (y) from earlier, and their preimages x l (y). By evaluating (51) at the points x l (y), we get a 2 l y 2 + 2αa l y + α 2 = c l y + d l y 2 + α 2 (52) Divide both sides by y, and take the limit y → 0. The lhs converges to 2αa l , while the rhs converges to c l . Since this is valid for all l ∈ [[1, n]], we conclude that c = 2αa. It also follows that d = a.
Step 4 Injecting this back into equation (51), and writing j d j h 2 , we are left with: By applying the trace operator to both sides of this equation, and rearranging terms, we get which is of the form trace(C T B(x)) = 0. This is a dot product on the space S n of n × n symmetric matrices (both C and B(x) are symmetric!), which is a vector space of dimension n(n+1)

2
. If we can show that the matrix C is orthognal to a basis of S n , then we can conclude that C = 0.
For this, let (e j ) 1≤j≤n be the Euclidean basis of R n , where each vector e j has one non-zero entry equal to 1 at index j, and let (E ij ) 1≤i≤n,1≤j≤n be the Euclidean basis of R n×n , where each matric E ij has only one non-zero entry equal to 1 at row i and column j. Now since h is surjective, there exists x j such that h(x j ) = e j , and h(x j )h(x j ) T = e j e T j = E jj . The n different x j give us our first n matrices we will use to construct a basis of S n . We now need to find n(n−1) 2 remaining basis matrices. For this, consider the sums (e j + e l ) 1≤j<l≤n , of which there is exactly n(n−1) 2 . Each of these sums of vectors have a primage x j,l by h, and h(x j,l )h(x j,l ) T = (e j + e l )(e j + e l ) T = E jj + E ll + (E il + E li ), which is a matrix in S n that is linearly independent of all E jj , and all other (e s + e t )(e s + e t ) T where (s, t) = (j, l) because they have non-zero entries at different rows and columns.
We have then found a total of n(n+1) 2 different vectors (x 1 , . . . , x n , x 1,2 , . . . , x n−1,n ) such that their images by hh T form a basis of S n . If we now evaluate equation (54) at each of these points, we find that the matrix aa T − diag(d) is orthogonal to a basis of S n , which implies that it is neceassirly equal to 0. This in turn implis that aa T is a diagonal matrix, and that a j a l = 0 for all j = l, which implies that at most one a j is non-zero.
Step 5 So far, we have proven that, among other things, A (2) i,j = 0 for all i, j. We now go back to equation (47), which we can write as: Both f and h are differentiable, and according to assumption 2, J f has rank n (it is full rank and n ≤ d). Thus, by differentiting the last equation, we conclude that A (1) has rank n, and is thus invertible.
ij . This is valid for all i ∈ [[1, n]]. Now since A (1) is invertible, the non-zero entry A (1) ij has to be in a different column for each row, otherwise some rows will be linearly dependent. Thus, there exists a permutation σ of [[1, n]], such that A (1) iσ(i) = 0, and we deduce that which concludes the proof.
From the second conclusion of step 3, we have that d = a, which implies that A (4) = A (1) , which in turn means that A is full rank.

A.6 PROOF OF THEOREM 2
Theorem 2. Let p(x|y) be a conditional probability density. Assume that X and Y are compact Hausdorff spaces, and that p(x|y) > 0 almost surely ∀(x, y) ∈ X × Y. Then for each ε > 0, there exists (θ, n) ∈ Θ × N, where n is the dimension of the feature extractor, such that sup x,y |p θ (x|y) − p(x|y)| < ε.
Proof. We consider here two cases.
Continuous auxiliary variable Recall the form of our model: By parameterizing each of f i , g i as neural networks, these functions can approximate continuous function on their respective domains arbitrarily well. According to Lemma 7, this implies that any continuous function on X × Y can be approximated arbitrarily well by a term of the form −f (x) T g(y).
Thus, any contiuous function can be approximated by log p θ (x|y) + log Z(y) for some θ, where Z(y) captures the difference in scale between the function in question and the normalized density p θ (x|y). We apply this result to log p(x|y): for any ε > 0, there exists (θ, n) ∈ Θ × N such that: Since p(x|y) > 0 a.s. on X × Y, log p(x|y) is finite and bounded. So is the term − n i=1 f i (x; θ)g i (y; θ). We can then use the fact that exp is Lipschitz on compacts to conclude for p(x|y), to conclude that: where K is the Lipschitz constant of exp, which concludes the proof.
Discrete auxiliary variable If y is discrete and Y is compact, then y only takes finitely many values. In this case, we do not need Lemma 7 for the proof. g(y) can simply be a lookup table, and we learn different approximations for each fixed value of y, since f has the universal approximation capability, which concludes the proof.
Denote by C(X ) (respectively C(Y) and C(X × Y)) the Banach algebra of continuous functions from X (respectively Y and X × Y) to R. For any subsets of functions F X ⊂ C(X ) and be the set of all linear combinations of products of functions from F X and F Y to R. The energy function defining our model belongs to this last set. Finally, universal approximation is expressed in terms of density: for instance, the set of functions F X have universal approximation of C(X ) if it is dense in it, i.e. for any function in C(X ), we can always find a limit of a sequence of functions of F X that converges to it. We mathematically express density by writing F X = C(X ).
Let F X (respectively F Y ) be the set of deep neural networks with input in X (respectively in Y). The universal approximation capability is summarised in the following Lemma.
Lemma 7 (Universal approximation capability). Suppose the following: (i) X and Y are compact Hausdorff spaces.
(ii) F X = C(X ) and F Y = C(Y) then F X ⊗ F Y = C(X × Y). All completions here are with respect to the infinity norm.
Proof. We prove this theorem in two steps: 1. We first prove that F X ⊗ F Y is dense in C(X ) ⊗ C(Y) using the hypotheses of Theorem 2.
2. we prove that C(X ) ⊗ C(Y) is dense in C(X × Y) Theorem 7.
Step 1 Let ε > 0. Let h ∈ C(X ) ⊗ C(Y). Then there exists k ∈ N and functions f i ∈ C(X ) and We then have Using this, we conclude that Step 2 We will use the Stone-Weirstrass theorem for this step. It is enough to show that: (iv) C(X ) ⊗ C(Y) separates points in X × Y (see Definition 3).
To prove (i), we use the fact that every finite product of compact spaces is a compact space, and every finite product of Hausdorff spaces is a Hausdorff space. Points (ii) and (iii) are easy to verify. To prove (iv), let (x, y) and (x , y ) be distinct points in X × Y. Assume that x = x (we proceed similarly if y = y ). Define the continuous function f ∈ C(X ) such that f (x) = 0 and f (x ) = 0. Then for g = 1 ∈ C(Y), we have f (x)g(y) = f (x) = 0 = f (x )g(y ).
All the conditions required to use the Stone-Weirestrass Theorem are verified, and we can conclude that Conclusion Combining the results of steps 1 and 2, we conclude that F X ⊗ F Y is dense in C(X × Y).
To prove this theorem, we first need to introduce some definitions and a useful theorem.
Definition 3. Let K be a compact Hausdorff space. Consider the Banach algebra C(K) equipped with the supremum norm f ∞ = sup t∈K |f (t)|. Then: 1. A ∈ C(K) is a unital sub-algebra if: (ii) for all f, g ∈ A and α, β ∈ R, we have αf + βg ∈ A and f g ∈ A.

A ⊂ C(K) separates points of
Theorem 7 (Stone-Weirstrass). Let K be a compact Hausdorff space, and A a unital sub-algebra of C(K) which separates points of K. Then A is dense in C(K).
Proof. A proof to this theorem can be found in many references, for instance Brosowski and Deutsch (1981).

B LATENT VARIABLE ESTIMATION IN GENERATIVE MODELS
Recall the generative model of IMCA: we observe a random variable x ∈ R d as a result of a nonlinear transformation h of a latent variable z ∈ R d whose distribution is conditioned on an auxiliary variable y that is also observed: We assume the latent variable in the IMCA model has a density of the form where µ is not necessarily factorial.
Further, we will suppose that the density p(z|y) belongs to the following subclass of the exponential families, introduced by Khemakhem et al. (2020): Definition 4 (Strongly exponential). We say that an exponential family distribution is strongly exponential if for any subset X of R the following is true: where Λ is the Lebesgue measure.
If we suppose that only n out of d components of the latent variable are modulated by the auxiliary variable y (equivalently, if we suppose that the parameters λ n+1:d (y) are constant), then we can write its density as The term e d i=n+1 Ti(zi) T λi is absorbed into µ(z). This last expression will be useful for dimensionality reduction.
To estimate the latent variables of the IMCA model, we fit an augmented version of our energy model where H(f (x)) = (H 1 (f 1 (x)), . . . , H d (f d (x))), and each H l is a (nonlinear) output activation. An example of such map is H l (x) = (x, x 2 ).
In this section, we present the proofs for the estimation of the Independently Modulated Component Analysis by an identifiable energy model. These proofs are based on similar ideas and techniques to previous proofs, but are different enough that we can't forgo them.

B.1 ASSUMPTIONS
We will decompose Theorem 5 into two sub-theorems, which will make the proof easier to understand, but also more adaptable into future work. For the sake of clarity, we will seperate its assumptions into smaller assumptions, and refer to them when needed in the proofs.
(iii) The sufficient statistics T i in (64) are differentiable, and the functions T ij ∈ T i are linearly independent on any subset of X of measure greater than zero.
(iv) There exist k + 1 distinct points y 0 , . . . , y k such that the matrix (v) We fit the model (67) to the conditional density p(x|y), where we assume the feature extractor f (x) to be differentiable, d-dimensional, and the pointwise nonlinearitiy H to be differentiable and k-dimensional, and the dimension of its vector-valued components H l to be chosen from (dim(T 1 ), . . . , dim(T d )) without replacement.
(vii) The mixing function h is a D 2 -diffeomorphisms.
(viii) The feature extractor f in (67) is a D 2 -diffeomorphism.
(vii') The mixing function h is a C 1 -diffeomorphism.
(viii') The feature extractor f in (67) is a C 1 -diffeomorphism, and the nonlinearities H l have a unique extremum.
(ix) Only n ≤ d components of the latent variable are modulated, and its density has the form (66).
The feature extractor f has the form f (x) = (f 1 (x), f 2 (x)) where f 1 (x) ∈ R n , and the auxiliary feature extractor g has the form g(y) = (g 1 (y), g 2 ) where g 1 (y) ∈ R n and g 2 is constant.

B.2 LEMMAS
We rely on the following Lemmas from Khemakhem et al. (2020), which we state below in the interest of completeness. Lemma 8. Consider an exponential family distribution with k ≥ 2 components. Then the components of the sufficient statistic T are linearly independent. Lemma 9. Consider a strongly exponential family distribution such that its sufficient statistic T is differentiable almost surely. Then T i = 0 almost everywhere on R for all 1 ≤ i ≤ k. Lemma 10. Consider a strongly exponential distribution of size k ≥ 2 with sufficient statistic T(x) = (T 1 (x), . . . , T k (x)). Further assume that T is differentiable almost everywhere. Then there exist k distinct values x 1 to x k such that (T (x 1 ), . . . , T (x k )) are linearly independent in R k . Lemma 11. Consider a strongly exponential distribution of size k ≥ 2 with sufficient statistic T. Further assume that T is twice differentiable almost everywhere. Then almost everywhere on R. Lemma 12. Consider n strongly exponential distributions of size k ≥ 2 with respective sufficient statistics T j = (T j,1 , . . . T j,k ), 1 ≤ j ≤ n. Further consider that the sufficient statistics are twice differentiable. Define the vectors e (j,i) ∈ R 2n , such that e (j,i) = 0, . . . , 0, T j,i , T j,i , 0, . . . , 0 , where the non-zero entries are at indices (2j, 2j + 1). Let x := (x 1 , . . . , x n ) ∈ R n . Then the matrix e(x) := (e (1,1) (x 1 ), . . . , e (1,k) (x 1 ), . . . e (n,1) (x n ), . . . , e (n,k) (x n )) of size (2n × nk) has rank 2n almost everywhere on R n .

B.3 PROOFS
As mentioned above, we decompose Theorem 5 into two smaller results, summarised in what follows by Theorems 5a and 5b. Theorem 5a. Assume assumptions (i)-(v) hold. Then, after convergence of our model p θ (x|y) to the true density p(x|y), we can recover the latent variables up to an invertible linear transformation and point-wise nonlinearities, i.e.
where A is an invertible matrix.
Proof. We fit our density model (67) to the conditional density p(x|y), setting the dimension of the feature extractor f to be equal to d, and the dimensions of the output nonlinearities H l chosen from (dim(T 1 ), . . . , dim(T d )), as per assumption (v): by doing the change of variable x = h(z), taking the log on both sides, we get: Let y 0 , . . . , y k be the points provided by assumption (iv) of the theorem, where k = i k i , and k i = dim(T i ). Define λ(y) = λ(y) − λ(y 0 ), Γ(y) = Γ(y) − Γ(y 0 ), g(y) = g(y) − g(y 0 ) and Z(y) = log Z(y) − log Z(y 0 ). We plug each of those y l in (72) to obtain k + 1 such equations. We subtract the first equation for y 0 from the remaining k equations to get for l = 1, . . . , k: The crucial point here is that the non factorial terms µ(g(x)) andμ(g(x)) cancel out when we take these differences. This is what allows us to generalize the identifiability results of nonlinear ICA to the context of IMCA.
Expressing (73) for all points y l in matrix form, we get: By assumption (iv), L is invertible, and thus we can write where c = L −T b and A = L −TLT .
To prove that A is invertible, we first take the gradient of equation (75) with respect to z. The Jacobian J T of T is a matrix of size k × d. Its columns are independent because each T i is only a function of z i , and thus the non-zero entries of each column are in different rows. This means that its rank is d . This is not enough to prove that A is invertible though. For that, we consider the functions T i for which k i > 1: for each of these functions, using Lemma 10, there exists points z We plug these vectors into equation (75)  . Now the matrix M is of size k × dk max , and it has exactly k independent columns by definition of the points z (j) . This means that M is of rank k, which in turn implies that rank(A) ≥ k. Since A is a k × k matrix, we conclude that A is invertible.
The theorem above shows a first step in identifiability which holds up to a linear transformation. This is similar to Hyvärinen et al. (2019), but here we allow for dependencies between components. We can further sharpen the result, in line with Khemakhem et al. (2020) even in this non-independent case as follows: Theorem 5b. Assume assumptions (i)-(v) hold. Further assume that either assumptions (vi)-(viii) or assumptions (vi')-(viii') hold. Then equation ( where γ is a permutation of [[1, d]] such that dim(H i ) = dim(T γ(i) ) and A i a square invertible matrix.
Proof. We prove this theorem seperately for both sets of assumptions.
The assumptions of Theorem 5a hold, and so we have for an invertible A ∈ R k×k . We will index A by four indices (i, l, a, b), where 1 ≤ i ≤ d, 1 ≤ l ≤ k i refer to the rows and 1 ≤ a ≤ d, 1 ≤ b ≤ k a to the columns.
Let dy (y). Each line of equation (78) can be written as: for i ≤ d, l ≤ k i . The first step is to show that v i (y) is a function of only one y ji , for all i ≤ d. by differentiating (79) with respect to y s , s ≤ d: and by differentiating (80) with respect to y t , s < t ≤ d: This equation is valid for all pairs (s, t), t > s. Define B a (y) : , M(y) := (B 1 (y), C 1 (y), . . . , B d (y), C d (y)), e (a,b) := (0, . . . , 0, T a,b , T a,b , 0, . . . , 0) ∈ R 2d , such that the non-zero entries are at indices (2a, 2a + 1) and e(y) := (e (1,1) (y 1 ), . . . , e (1,k1) (y 1 ), . . . , e (d,1) (y d ), . . . , e (d,k d ) (y d )) ∈ R 2d×k . Then by grouping equation (81) for all valid pairs (s, t) and pairs (i, l) and writing it in matrix form, we get: Now by Lemma 12, we know that e(y) has rank 2d almost surely on Z. Since A is invertible, it is full rank, and thus rank(e (y)A) = 2d almost surely on Z. It suffices then to multiply by its pseudo-inverse from the right to get In particular, C a (y) = 0 for all 1 ≤ a ≤ d. This means that the Jacobian of v at each y has at most one non-zero entry in each row. By invertibility and continuity of J v , we deduce that the location of the non-zero entries are fixed and do not change as a function of y. We deduce that there exists a permutation σ of [[1, d]] such that each of the v i (y) = v i (y σ(i) ), and the same would apply to v −1 . Without any loss of generality, we assume that σ is the identity.
This function is a pointwise function because H and v −1 are such functions. Plugging this back into equation (78) yields: The last equation is valid for every component: By differentiating both sides with respect to z s where s = i we get By Lemma 8, we get A i,l,s,b = 0 for all 1 ≤ b ≤ k. Since (86) is valid for all l and all s = i, we deduce that the matrix A has a block diagonal form: which achieves the proof.
One-dimensional sufficient statistics: assumptions (vi')-(viii') We now suppose that k i = 1, ∀i.  (79) for every component as: which is the same as equation (45) of Khemakhem et al. (2020). All the assumptions required to prove their theorem are met in our case, and the rest of their proof would simply apply here to prove that A is a permutation matrix.
Then f 1 recovers only the modulated latent components as per Theorem 5b.
Proof. The proof of Theorem 5a in this case is unchanged. Simply, we update the total dimension of matrix L here to k = n i=1 dim(T i ). when we evaluate equation (72) on these points y 0 , . . . , y k , the constant term g 2 and the non-modulated components cancel out, and we are left with the equatioñ We then use similar arguments to the proof of Theorem 5a to conclude that where A ∈ R n a square invertible matrix. At this point, we can make the same conclusion as Theorem 5a, while reducing the dimension of the latent space.
We now explain how we can extend Theorem 5b to the lower dimensional latent space case. Note that we still assume that f = (f 1 , f 2 ) is a diffeomorphism per assumptions (viii) and (viii'). We can then still define v = (f • h) −1 .
We consider now two cases like in the proof of Theorem 5b.
One-dimensional sufficient statistics Let D = A −1 and h i,a = D i,a H a (y a ) − D i,a c a . We can still write equation (90) like equation (88) as for all i ≤ n. The assumptions required for the proof are still met, despite reducing the dimension from d to n. This interesting fact is also used for the proof of Theorem 4 as well, which achieves this part of the proof.

Multi-dimensional sufficient statistics We rewrite equation (90)
for all i ≤ n, l ≤ k i . We proceed similarly to the proof of Theorem 5b, replacing all mentions of d by n and keeping all differentiations to indices t, s ≤ n, up to equation (83), after which we can conclude that v s i v t i = 0 for all i ≤ n, and all s, t ≤ n. This is not enought to conclude that each of the v i is only function of one y ji .
For that, we go back to equation (92) and differentiate it with respect to y s , s > n: which is valid for all i ≤ n, l ≤ k i . Since A is invertible, we can conclude that T a,b (v a (y))v s a (y) = 0 for all a ≤ n and s > n. Since we only consider strongly exponential distributions (assumption (iii)), and using proposition 9, we conclude that T a,b (v a (y)) = 0 almost everywhere, and that v s a (y) = 0, for all s > n. This, in addition to the fact that v s i v t i = 0 for all i ≤ n, and all s, t ≤ n allows us to conclude that the first n components of v are each only a function of one different y j because v is a diffeomorphism and its Jacobian is continuous. Finally, we can use this fact to deduce that A is a block permutation matrix, which achieves the proof.

C INDEPENDENTLY MODULATED COMPONENT ANALYSIS
As mentioned in section 4, linear latent variable models (Hyvärinen and Oja, 2000) and more recently nonlinear latent variable models may be identifiable provided some additional auxiliary variables (Khemakhem et al., 2020;Hyvärinen et al., 2019). The purpose of this auxiliary variable serves to introduce additional constraints over the distribution over latent variables, which are typically required to be conditionally independent given the auxiliary variable. This avenue of research has thus formalized the trade-off between expressivity of the mapping between latents to observations (from linear to nonlinear) and distributional assumptions over latent variables (from independent to conditionally independent given auxiliary variables).
We would like to relax the assumption of independence while maintaining identifiability, resulting in the framework of Independently Modulated Component Analysis (IMCA). In this section of the Appendix, we will give a detailed analysis of the IMCA model indepenedently of any estimation method, drawing parallels to the identifiability results of the nonlinear ICA model presented in Khemakhem et al. (2020).

C.1 DEFINITION OF THE GENERATIVE MODEL
Assume we observe a random variable x ∈ R d as a result of a nonlinear transformation h of a latent variable z ∈ R d whose distribution is conditioned on an auxiliary variable y that is also observed: The main modelisation assumption we make is on the latent variable distribution, given by the following definition, where u is a dummy variable.
Definition 5 (Exponentially factorial distributions). We say that a multivariate exponential family distribution is exponentially factorial if its density p(u) has the form We assume that the latent variable in the IMCA model has a conditional exponentially factorial distribution, where the parameters of the exponential family are a function of the auxiliary variable y: Equations (94) and (95) together define the nonparametric IMCA model with parameters (h, T, λ, µ). Most importantly, we allow for an arbitrary base measure µ(z), i.e. the components of the latent variable must no longer be independent, as µ doesn't necessarily factorise across dimensions. The crucial assumption is that the components of the latent variables are independently modulated given the auxiliary variable y, and that through the term exp( i T i (z i ) T λ i (y)).

C.2 IDENTIFIABILITY
The concenpt of identifiability is core to this work. As such, it is important to understand the different views one can have of this concept.
Often, the identifiability form we can prove for a model is equality of the parameters up to some indeterminacies. This can be understood as an equivalence relation between parameters. Identifiability in this context implies that the equivalence class of the ground truth parameter can be uniquely recovered from observations. This is relevant only if the definition of the equivalence class is sufficiently narrow and specific to be able to make meaningful conclusions.
One example of such equivalence relations can be found in linear ICA: the mixing matrix is uniquely recovered up to a scaled permutation. The permutation is irrelevant, and the scaling is circumvented by whitening the data. This is a good example of an equivalence class that doesn't restrict the practical utility of the ICA model.
An example of indeterminacy which is relevant to us here can be found in variational inference of latent variable models: two parameters are equivalent if they map to the same inference distribution (Khemakhem et al., 2020). This is the definition we will be using in this work. We will say that a generative model is identifiable if we can uniquely recover the latent variables, as given by the following definition.
Definition 6. Consider two different sets of parameters (h, T, λ, µ) and (h,T,λ,μ), defining two densities p and p . We say that the IMCA model is strongly identifiable if where γ is a permutation, A i is an invertible matrix, and b i a vector, ∀i ∈ [[1, d]].
We say that it is weakly identifiable if where A is an invertible matrix, and b a vector.

C.3 THEORETICAL ANALYSIS
In this section, we develop the theory of IMCA. We will give sufficient conditions that guarantee a strong identifiability of the latent components, and discuss a degenerate case where we only obtain a weaker form of identifiability.

C.3.1 Definitions
We will first define some sets of distributions which are subsets of the exponential family distribution. We will use u as a dummy variable, and introduce the definitions for the unconditional case. Note that all these definitions apply to the conditional case, when the parameters of the exponential family are a function of an auxiliary variable y. For completeness, we restate here Definition 4.
We say that that a multivariate distribution is strongly exponential if all its univariate marginals are.
In other words, the density of a strongly exponential distribution has almost surely the exponential component in its expression and can only be reduced to the base measure on a set of measure zero. This definition is very general, and is satisfied by all the usual exponential family distributions like the Gaussian, Laplace, Pareto, Chi-squared, Gamma, Beta, etc. We will only prove identifiability results for strongly exponential families. The non-strongly exponential case will be explored in future work.
There is a certain class of exponential families for which we can only prove a weak form of identifiability. Loosely speaking, this is because this class doesn't constrain the latent space enough.
Definition 8 (Quasi-location exponential distributions). We say that a univariate exponential family distribution with density p(u) = µ(u)e T(u) T θ−Γ(θ) is in the quasi-location family if: (i) dim(T) = 1 (ii) T is monotonic (either non-decreasing or non-increasing) We say that that a multivariate distribution is quasi-location exponential if all its univariate marginals are.
As a simple illustration, the Gaussian family with fixed variance is a quasi-location family, but with fixed mean it is not. This is because in the first case, the sufficient statistic is T (u) = u which is a monotonic scalar function, while in the second case it is T (u) = u 2 , a non-monotonic scalar function.

C.3.2 Identifiability of the General Case
As mentioned in section 4, the IMCA model described by equations (94) and (95) generalises previous nonlinear ICA models by relaxing the independence assumption required for the latent variables. We propose here to extend the identifiability theory of nonlinear ICA developed in Hyvärinen et al. (2019); Khemakhem et al. (2020) to this new framework.
We start by providing a weaker form of identifiability guarantee that applies to the general case, including quasi-location families.
(II) The mixing function h : R d → R d is invertible.
(III) The conditional latent distribution p(z|y) is strongly exponential (definition 7), and its sufficient statistic is differentiable.
This theorem extends the basic identifiability result of Khemakhem et al. (2020, Theorem 1). It is fundamental as it proves a general identifiability results without the restriction of having independent latent variables. This was previously not considered to be possible and could only be demonstrated in very specific circumstances and under very restrictive additional assumptions (e.g., Monti and Hyvärinen (2018) require both non-negativity and orthonormality of a mixing matrix in the linear case). In the nonlinear case, to prove Theorem 8, we still require that the latent variables are only dependent through the base measure, while still being independently modulated through the auxiliary variable y. This (and the necessity of having an auxiliary variable) is the price to pay for obtaining identifiability in a nonlinear setting.

C.3.3 Identifiability of the Non Quasi-location Family
The identifiability result of Theorem 8 is weak because of the presence of the linear transformation A in equation (97). It turns out that by excluding the quasi-location family (definition 8), we can remove this matrix and achieve a stronger form of identifiability. The main technical result of this paper is the following.
Theorem 9. Assume that the assumptions of Theorem 8 hold. Further assume one of the two following sets of assumptions: (V) The sufficient statistic in (95) is twice differentiable and dim(T l ) ≥ 2, ∀l.
(VI) The mixing function h is a D 2 -diffeomorphism 6 .
or 6 invertible, all second order cross-derivatives of the function and its inverse exist but aren't necessarily continuous (V)' dim(T l ) = 1 and T l is non-monotonic ∀l.
(VI)' The mixing function h is a C 1 -diffeomorphism 7 .
Then the IMCA model is strongly identifiable.
This form of identifiability mirrors the strongest results proven in the nonlinear ICA (Khemakhem et al., 2020, Theorems 2,3), without requiring that the latent components be independent. As far as we know, this is the first proof of the kind for nonlinear representation learning. We further note that this theorem generalizes even existing identifiability theory of the linear case. The mixed case where we have both cases where some sufficient statistics are of dimension greater than 2 and some are univariate and non-monotonic will be studied in future work.

C.4 ESTIMATION OF IMCA BY SELF-SUPERVISED LEARNING
A recent development in nonlinear ICA is given by Hyvärinen et al. (2019) where the authors assume they observe data x = h(z) following a noiseless conditional nonlinear ICA model p(z|y) = i p i (z i |y) For estimation, they rely on a self-supervised binary discrimination task based on randomization to learn the unmixing function. More specifically, from a dataset of observations and auxiliary variables pairs D = {x (i) , y (i) }, they construct a randomized dataset D * = {x (i) , y * } where y * is randomly drawn from the observed distribution of y. To distinguish between both datasets, a deep logistic regression is used. The last hidden layer of the neural network is a feature extractor whose purpose is to extract the relevant features which will allow to distinguish between the two datasets. Surprisingly, this estimation technique works for IMCA, and is summarized by the following theorem. Theorem 10. Self-supervised nonlinear ICA estimation algorithms presented in Hyvärinen and Morioka (2016) Consider two different sets of parameters (h, T, λ, µ) and (h,T,λ,μ), defining two conditional latent densities p(z|y) andp(z|y). Suppose that the density of the observations arising from these two different models are equal: p(x|y) =p(x|y) log p(g(x)|y) − log det J −1 h (x) = log p(g(x)|y) − log |det Jg(x)| (100) log µ(g(x)) + T(g(z)) T λ(y) − Γ(y) − log |det J g (x)| = logμ(g(x)) +T(g(z)) Tλ (y) −Γ(y) − log |det Jg(x)| Let y 0 , . . . , y k be the points provided by assumption (IV) of the theorem for T, where k = i k i , and k i = dim(T i ).
We plug each of those y l in (101) to obtain k + 1 such equations. Then, we subtract the first equation for y 0 from the remaining k equations to get for l = 1, . . . , k: T(z) T (λ(y l ) − λ(y 0 )) − G(y l ) =T(z) T (λ(y l ) −λ(y 0 )) −G(y l ) where we grouped terms that are only a function of y l in G andG.
Most importantly, both base measure terms disappear after taking the differences, which is the key enabler of identifiability in the IMCA framework.
The rest of the proof is similar to the proof of Khemakhem et al. (2020, Theorem 1). The only difference is that we don't restrict the sufficient statistics to have equal dimensions, and so we can't use the proof technique from Khemakhem et al. (2020, Theorem 1) without any modification. We present an alternative technique in the proof of Theorem 5, which we refer too for more details. We then conclude that T(h −1 (x)) = AT(h −1 (x)) + b which implies that the model is weakly identifiable. 7 invertible, all partial derivatives of the function and its inverse exist and are continuous