Nonparametric Evaluation of Noisy ICA Solutions

Authors: syamantak  kumar
purnamrita  sarkar
peter  bickel
derek  bean
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 19:22:31

Abstract: Independent Component Analysis (ICA) was introduced in the 1980's as a model for Blind Source Separation (BSS), which refers to the process of recovering the sources underlying a mixture of signals, with little knowledge about the source signals or the mixing process. While there are many sophisticated algorithms for estimation, different methods have different shortcomings. In this paper, we develop a nonparametric score to adaptively pick the right algorithm for ICA with arbitrary Gaussian noise. The novelty of this score stems from the fact that it just assumes a finite second moment of the data and uses the characteristic function to evaluate the quality of the estimated mixing matrix without any knowledge of the parameters of the noise distribution. In addition, we propose some new contrast functions and algorithms that enjoy the same fast computability as existing algorithms like FASTICA and JADE but work in domains where the former may fail. While these also may have weaknesses, our proposed diagnostic, as shown by our simulations, can remedy them. Finally, we propose a theoretical framework to analyze the local and global convergence properties of our algorithms.

Main paper: 
 introduction: 
 Independent Component Analysis (ICA) was introduced in the 1980's as a model for Blind Source Separation (BSS) [13, 12] , which refers to the process of recovering the sources underlying a mixture of signals, with little knowledge about the source signals or the mixing process. It has become a powerful alternative to the Gaussian model, leading to PCA in factor analysis. A good account of the history, many results, and the role of dimensionality in noiseless ICA can be found in [28, 5] . In one of its first forms, ICA is based on observing n independent samples from a k-dimensional population:EQUATIONwhere B ∈ R d×k (d ≥ k) is an unknown mixing matrix, z z z ∈ R k is a vector of statistically independent mean zero "sources" or "factors" and g g g ∼ N (0 0 0, Σ) is a d dimensional mean zero Gaussian noise vector with covariance matrix Σ, independent of z z z. We denote diag (a a a) := cov (z z z) and S := cov (x x x) = B diag (a a a) B T + Σ. The goal of the problem is to estimate B.ICA was initially developed in the context of g g g = 0, now called noiseless ICA. In that form, it spawned an enormous literature [17, 39, 40, 1, 8, 27, 35, 37, 23, 11, 6, 42] and at least two wellestablished algorithms FASTICA [29] and JADE [9] which are widely used. The general model with nonzero g g g and unknown noise covariance matrix Σ, known as noisy ICA, has developed more slowly, with its own literature [4, 50, 49, 30, 31, 41] .In the following sections, we will show that various ICA and noisy ICA methods have distinct shortcomings. Some struggle with heavy-tailed distributions and outliers, while others require approximations of entropy-based objectives, which have their own challenges (see eg. [32] ).Although methods for noiseless cases may sometimes work in noisy settings [49] , they are not always reliable (e.g., see Figure 1 in [49] and Theorem 2 in [11] ). Despite the plethora of algorithms for noisy and noiseless ICA, the literature has largely been missing a diagnostic to decide which algorithm to pick for a given dataset. This is our primary goal.The paper is organized as follows. Section 2 contains background and notation. Section 3.1 contains the independence score, its properties, and a Meta algorithm 1 that utilizes this score. Section 3.3 introduces new contrast functions based on the natural logarithm of the characteristic function and moment generating function of x x x T u u u. Section 4 presents the global and local convergence results for noisy ICA. Section 5 demonstrates the empirical performance of our new contrast functions as well as the Meta algorithm applied to a variety of methods for inferring B. 
background and overview: 
 This section contains a general overview of the main concepts of ICA.Notation: We denote vectors using bold font as v v v ∈ R d . The dataset is represented as x x x (i) i∈[n] for x x x (i) ∈ R d . Scalar random variables and matrices are represented using upper case letters, respectively. I k denotes the kdimensional identity matrix. Entries of vector v v v (matrix M ) are represented as v i (M ij ). M (:, i) (M (i, :)) represents the i th column (row) of matrix M . M † represents the Moore-Penrose inverse of M . ∥.∥ represents the Euclidean l 2 norm for vectors and operator norm for matrices. ∥.∥ F similarly represents the Frobenius norm for matrices. E [.] ( E [.]) denotes the expectation (empirical average); x x x ⊥ ⊥ y y y denotes statistical independence. ∥.∥ ψ2 represents the subgaussian Orlicz norm (see Def. 2.5.6 in [46] ). X n = O P (a n ) implies that Xn an is stochastically bounded i.e, ∀ϵ > 0, there exists a finite M such that sup n P Xn an > M ≤ ϵ (See [45] ).Identifiability: One key issue in Eq 1 is the identifiability of A := B -11 , which holds up to permutation and scaling of the coordinates of z z z if, at most, one component of z z z is Gaussian (see [14, 43] ). If d = k and Σ are unknown, it is possible, as we shall see later, to identify the canonical versions of B and Σ. For d > k and unknown Σ, it is possible (see [15] ) to reduce to the d = k case. In this work, we assume d = k. For classical factor models, when z z z is Gaussian, identifiability can only be identified up to rotation and scale, leading to non-interpretability. This suggests fitting strategies focusing on recovering B -1 through nongaussianity as well as independence. In the noisy model (Eq 1), an additional difficulty is that recovery of the source signals becomes impossible even if B is known exactly. This is because, the ICA models x = B(z + r) + g and x = Bz + (Br + g) are indistinguishable (see [48] pages 2-3). This is resolved in [48] via maximizing the Signal-to-Interference-plus-Noise ratio (SINR). In classical work on ICA, a large class of algorithms (see [37] ) choose to optimize a measure of non-gaussianity, referred to as contrast functions.Contrast functions: Methods for estimating B typically optimize an empirical contrast function. Formally, we define a contrast function f (u u u|P ) by f : B k × P → R where B k is the unit ball in R k and P is essentially the class of mean zero distributions on R d . More concretely, for suitable choices of g, f (u u u|P ) ≡ f (u u u|x x x) = E x x x∼P [g(u u u T x x x)] for any random variable x x x ∼ P . The most popular example of a contrast function is the kurtosis, which is defined as f (u u u|P ) = E x x x∼P (u u u T x x x) 4 E (u u u T x x x) 2 2 -3.The aim is to maximize an empirical estimate, f (u u u|P ). Notable contrast functions are the negentropy, mutual information (used by INFOMAX [7] ), the tanh function, etc. (See [37] for further details).Fitting strategies: There are two broad categories of the existing fitting strategies. (I) Find A such that the components of Ax x x are both as nongaussian and as independent as possible. See, for example, the multivariate cumulant-based method JADE [9] and characteristic function-based methods [21, 11] . The latter we shall call the PFICA method. (II) Find successively the j th row of A, denoted by A(j, :) ∈ R d , j = 1, . . . , k such that A(j, :) T x x x are independent and each as nongaussian as possible, that is, estimate A(1, :), project, and then apply the method again on the residuals successively. The chief representative of these methods is FastICA [27] based on univariate kurtosis.From noiseless to noisy ICA: In the noiseless case one can first prewhiten the dataset using the empirical variance-covariance matrix, Σ, of x x x (i) , using Σ-1/2 . Therefore, WLOG, B can be assumed to be orthogonal. Searching over orthogonal matrices not only simplifies algorithm design for fitting 3.1 Independence score While there are many tests for independence [22, 34, 42, 6] , we choose the Characteristic functionbased one because it is easy to compute and has been widely used in normality testing [19, 20, 25] and ICA [11, 21, 44] . Let us start with noiseless ICA to build intuition. Say we have estimated the inverse of the mixing matrix, i.e. B -1 using a matrix F . Ideally, if F = B -1 , then Fx x x = z z z, where z z z is a vector of independent random variables. Kac's theorem [33] tells us that E exp(it T z z z) = k j=1 E [exp(it j z j )] if and only if z i are independent. In [21] , a novel characteristic functionbased objective (CHFICA), is further analyzed and studied by [11] (PFICA). Here one minimizes|E exp(it T Fx x x) - k j=1 E [exp(it j (Fx x x) j )] | over F .We refer to this objective as the uncorrected independence score. We propose to adapt the CHFICA objective using estimable parameters, S, to the noisy ICA setting. We will minimize:EQUATION)where S = E[x x xx x x T ] and can be estimated using the sample covariance matrix. Hence, we do not require knowledge of any model parameters. We refer to this score as the (corrected) independence score. The second terms in JOINT and PRODUCT (Eq 2) corrects the original score by canceling out the additional terms resulting from the Gaussian noise using the covariance matrix S of the data (estimated via the sample covariance). See [21] for a related score requiring knowledge of the unknown Gaussian covariance matrix Σ. We are now ready to present our first theoretical result for consistency of ∆(t t t, F |P ).Theorem 1. If F ∈ R k×k is invertible and the joint and marginal characteristic functions of all independent components,{z i } i∈[k] , are twice-differentiable, then ∀t t t ∈ R k , ∆(t t t, F |P ) = 0 iff F = DB -1where D is a permutation of an arbitrary diagonal matrix.Unfortunately, F is not uniquely defined if Σ is unknown as we have noted in Section 1. The proof of Theorem 1 is provided in the Appendix Section A.1. When using this score in practice, S is replaced with the sample covariance matrix of the data, and the expectations are replaced by the empirical estimates of the characteristic function. The convergence of this empirical score, ∆(t, F | P ), to the population score, ∆(t, F |P ), is given by the following theorem.Theorem 2. Let F := {F ∈ R k×k : ∥F ∥≤ 1}, x x x ∼ subgaussian(σ) and C k := max(1, k log (n) Tr(S)). Then, we have supF ∈F |E t t t∈N (0,I k ) ∆(t t t, F |P ) -E t t t∈N (0,I k ) ∆(t t t, F | P )|= O P   k 2 ∥S∥max(k, σ 4 ∥S∥) log 2 (nC k ) n  This bound shows that uniformly over F , the empirical average E t t t∈N (0,I k ) ∆(t t t, F | P ), is close to the population score. This guarantees that as long as the difference between the population scores of the two candidate algorithms is not too small, the meta-algorithm can pick up the better score. Remark 1 (Subgaussianity assumption). The subgaussianity assumption in Theorem 2 simply ensures the concentration of the sample covariance matrix since it is used in the score (see Theorem A.3, line 825). It can be relaxed if the concentration in the operator norm is ensured. Appendix Section A.2.1 contains experiments with more super-Gaussian source signals.The proof of this result is deferred to the Appendix section A.1.1. It is important to note that ∆(t t t, F | P ) is not easy to optimize. As we show in Section 3.2, objective functions that are computationally more amenable to optimize for ICA, e.g., cumulants, satisfy some properties (see Assumption 1). The independence score does not satisfy the first property (Assumption 1 (a)). In the noiseless case, whitening reduces the problem to B being orthogonal. This facilitates efficient gradient descent in the Stiefel manifold of orthogonal matrices [11] . However, in the noisy setting, the prewhitening matrix contains the unknown noise covariance Σ, making optimization hard. Furthermore, as noted in Remark 2, in practice, it is better to use E t t t∼N (0,I k ) ∆(t t t, F | P ) [11, 21, 25] instead of using a fixed vector, t t t, to evaluate ∆(t t t, F | P ). Any iterative optimization method requires repeatedly estimating this expectation at each gradient step. Therefore, instead of using this score directly as the optimization objective, we use it to choose between different contrast functions after extracting the demixing matrix with each. This process is referred to as the Meta algorithm (Algorithm 1).Algorithm 1 Meta-algorithm for choosing best candidate algorithm.Input: Algorithm list L, Data X ∈ R n×k for j in range[1, size (L)] do B j ← L j (X) {Extract mixing matrix B j using j th candidate algorithm} Remark 2 (Averaging over t t t). Algorithm 1 averages the independence score over t t t ∼ N (0, I k ). While Eq 2 defines the independence score for one direction t t t, in practice, however, there may be some directions such that a non-gaussian signal has a small score in that direction. Hence, it is desirable to average over t t t, following the convention in [11, 21, 25] .δ j ← E t t t∼N (0,I k ) ∆ t, B -1 j | P end for i * ← arg min j∈[size(L)] [δ j ] return B i *To gain an intuition of the independence score, we conduct an experiment where we use the dataset mentioned in Section 5 (Figure 2b ) and compute the independence score for B ′ = ϵB + (1 -ϵ) I, ϵ ∈ (0.5, 1). As we increase ϵ, B ′ approaches B, and hence the Amari error d B ′ ,B (see Eq 8) decreases. Figure 1 shows the independence score versus Amari error, indicating that the independence score accurately predicts solution quality even without knowing the true B.In this section, we provide the proof of Theorem 4. Recall the ICA model from Eq 1,x x x = Bz z z + g g g, diag (a a a) := cov (z z z) , Σ := cov (g g g) = B diag (a a a) B T + ΣWe provide the proof for the CGF-based contrast function. The proof for the CHF-based contrast function follows similarly. From the Proof of Theorem A.1 we have,∇f (u u u|P ) = E exp u u u T Bz z z Bz z z E [exp (u u u T Bz z z)] -B diag (a a a) B T u u u, ∇ 2 f (u u u|P ) = E exp u u u T Bz z z Bz z zz z z T B T E [exp (u u u T Bz z z)] - E exp u u u T Bz z z Bz z z E exp u u u T Bz z z z T B T E [exp (u u u T Bz z z)] 2 -B diag (a a a) B T = B             E exp u u u T Bz z z z z zz z z T E [exp (u u u T Bz z z)] - E exp u u u T Bz z z z z z E exp u u u T Bz z z z z z T E [exp (u u u T Bz z z)] 2 Covariance of z z z ∼ P (z z z) exp u u u T Bz z z E [exp (u u u T Bz z z)] -diag (a a a)             B TThe new probability densityP (z z z) exp u u u T Bz z z E [exp (u u u T Bz z z)]is an exponential tilt of the original pdf, and since {z i } d i=1 are independent, and the new tilted density also factorizes over the z i 's, therefore, the covariance under this tilted density is also diagonal.Let C := BD 0 B T . We denote the i th column of B as B .i . Define functionsr i (u u u) := ln E exp u u u T B .i z i -Var (z i ) u u u T B .i B T .i u u u 2 g i (u u u) := ∇r i (u u u) Then f (u u u|P ) = d i=1 r i (u u u), ∇f (u u u|P ) = d i=1 g i (u u u). For fixed point iteration, u u u k+1 = ∇f (C -1 u u u k |P ) ∥∇f (C -1 u u u k |P )∥ .Consider the function ∇f C -1 u u u|P . Let e i denote the i th axis-aligned unit basis vector of R n . Since B is full-rank, we can denote α α α = B -1 u u u. We have,∇f C -1 u u u|P = d i=1 g i C -1 u u u = d i=1 E exp u u u T C -1 T B .i z i B .i z i E exp u u u T (C -1 ) T B .i z i -Var (z i ) B .i B T .i C -1 u u u = d i=1 E exp u u u T B T -1 D -1 0 e e e i z i B .i z i E exp u u u T (B T ) -1 D -1 0 e e e i z i -Var (z i ) B .i e e e T i D -1 0 B -1 u = d i=1   E exp αi (D0) ii z i z i E exp αi (D0) ii z i -Var (z i ) α i (D 0 ) ii   B .iFor t ∈ R, define the function q i (t) : R → R as -q i (t) := E exp t (D0) ii z i z i E exp t (D0) ii z i -Var (z i ) t (D 0 ) ii Note that q i (0) = E [z i ] = 0. For t = (t 1 , t 2 , • • • t d ) ∈ R d , let q (t) := [q 1 (t 1 ) , q 2 (t 2 ) • • • q d (t d )] T ∈ R d . Then, ∇f C -1 u|P = Bq (α α α) Therefore, if u u u t+1 = Bα α α t+1 , then Bα α α t+1 = Bq (α α α t ) ∥Bq (α α α t )∥ =⇒ α α α t+1 = q (α α α t ) ∥Bq (α α α t )∥ =⇒ ∀i ∈ [d], (α α α t+1 ) i = q i ((α α α t ) i ) ∥Bq (α α α t )∥At the fixed point, α α α * = e e e 1 ∥Be e e 1 ∥ and Bq q q (α α α * ) ∥Bq q q (α α α * )∥ 2= Be e e 1 . Therefore,∀i ∈ [2, d], (α α α t+1 ) i -α * i = (α α α t+1 ) i = q i ((α α α t ) i ) ∥Bq (α α α t )∥ (A.26) for i = 1, (α α α t+1 ) 1 -α * 1 = q 1 ((α α α t ) i ) ∥Bq (α α α t )∥ - 1 ∥Be e e 1 ∥ (A.27)Note the smoothness assumptions on q i (.) mentioned in Theorem 4, ∀i ∈ [d],1. sup t∈[-∥B -1 ∥2,∥B -1 ∥2] |q i (t) |≤ c 1 2. sup t∈[-∥B -1 ∥2,∥B -1 ∥2] |q ′ i (t) |≤ c 2 3. sup t∈[-∥B -1 ∥2,∥B -1 ∥2] |q ′′ i (t) |≤ c 3Since ∀k, ∥u u u k ∥= 1, therefore, ∀k, ∥α α α t ∥≤ ∥B -1 ∥ 2 . We seek to prove that the sequence {α α α t } n k=1 converges to α α α * .Question: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes]Justification: In the Experiments (Section 5, there are some noise/sample-size regimes where other algorithms may perform slightly better than our proposed contrast functions, but we show that even then our diagnostic can pick out the best algorithm.Guidelines:• If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable. • Depending on the contribution, reproducibility can be accomplished in various ways.For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed. • While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. , with an open-source dataset or instructions for how to construct the dataset). (d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility.In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results. 
desired properties of contrast functions: 
 The following properties are often useful for designing provable optimization algorithms for ICA.Assumption 1 (Properties of contrast functions). Let f be a contrast function defined for a mean zero random variable X. Then,(a) f (u|X + Y ) = f (u|X) + f (u|Y ) for X ⊥ ⊥ Y (b) f (0|X) = 0, f ′ (0|X) = 0, f ′′ (0|X) = 0 (c) f (u|G) = 0, ∀u for G ∼ N (0, 1) (d) WLOG, f (u|X) = f (-u|X) (symmetry)Properties (a) and (c) ensure that f (u|G + X) = f (u|X) for non-gaussian X and independent non-gaussian G, which means that the additive independent Gaussian noise does not change f . Property (c) ensures that |f (u|G)| is minimal for a Gaussian. Property (d) holds without loss of generality because one can always symmetrize the distribution. 2In Section 2, we provided a discussion of the individual shortcomings of different contrast functions for existing contrast functions. Before we introduce new contrast functions in this section, we revisit the algorithmic issues posed by the added Gaussian noise with unknown Σ in Eq 1.Prewhitening the data is challenging for noisy ICA because E[x x xx x x T ] includes the unknown Gaussian covariance matrix Σ. [4] show that it is enough to quasi orthogonalize the data, i.e., multiply it by a matrix, which makes the data have a diagonal covariance matrix (not necessarily the identity). Subsequent work [50, 49] uses cumulants and the Hessian of f (u u u) to construct a matrix C := BDB T where D is a diagonal matrix, and then use this matrix to achieve quasi-orthogonalization. In general, D may not be positive semidefinite (e.g. when components of z z z have kurtosis of different signs). To remedy this, in [49] , the authors propose computing the C matrix using the Hessian of f (u u u) at some fixed unit vector and then perform an elegant power method in the pseudo-Euclidean space:EQUATIONA pseudo-Euclidean space is a generalization of the Euclidean space used by [48] . Here the produce between vectors u u u, v v v is given as u ⊤ Av, Algorithm 2 Pseudo-Euclidean Power Iteration for ICA (Algorithm 2 in [48] ) B is the recovered matrix for the ICA model in Eq 1. Ã is a running estimate of B † .Input: C ∈ R k×k , ∇f Ã ← 0, B ← 0 for j in range[1, k] do Draw u uniformly at random from S k-1 while Convergence (up to sign) do u ← B Ãu u ← ∇f C † u / ∇f C † u 2 end while Bj ← u, Ãj ← C † Bj / C † Bj ⊤ Bj ⊤ end for return B, ÃIn this section, we present two new contrast functions based on the logarithm of the characteristic function (CHF) and the cumulant generating function (CGF). These do not depend on a particular cumulant, and the first only requires a finite second moment, which makes it suitable for heavytailed source signals. We will use f (u u u|P ) or f (u u u|x x x) where x x x ∼ P interchangeably to represent the population contrast function. In constructing both of the following contrast functions, we use the fact that, like the cumulants, the CGF and CHF-based contrast functions satisfy Assumption 1 (a). To satisfy Assumption 1 (c), we subtract out the part resulting from the Gaussian noise, which leads to the additional terms involving u u u ⊤ Su u u.CHF-based contrast function: Recall that S denotes the covariance matrix of x x x. We maximize the absolute value of following:EQUATIONThe intuition is that this is exactly zero for zero mean Gaussian data and maximizing this leads to extracting non-gaussian signal from the data. It also satisfies Assumption 1 a), b) and c).The cumulant generating function also has similar properties (see [53] for the noiseless case). In this case, we maximize the absolute value of following:EQUATIONLike CHF, this vanishes iff x x x is mean zero Gaussian, and satisfies Assumption 1 a), b) and c). However, it cannot be expected to behave well in the heavy-tailed case. In the Appendix (Section A.1, Theorem A.1), we show that the Hessian of both functions obtained from Eq 4 and Eq 5 evaluated at some u u u is, in fact, of the form BDB T where D is some diagonal matrix.Putting all the cases together in Eq A. 45 , 3) Laplace(0, 0.05) (mean 0 and standard deviation 0.05), 4) Exponential(5), and 5) and 6) a Student's t-distribution with 3 and 5 degrees of freedom, respectively. The Meta algorithm closely follows the best candidate algorithm even in the presence of many super-Gaussian signals. In this section, we provide additional experiments for image-demixing using ICA. We mix images using flattening and linearly mixing them with a 4 × 4 matrix B (i.i.d entries ∼ N (0, 1)) and Wishart noise (ρ = 0.001). Demixing is performed using the SINR-optimal demixing matrix (see Section 2) and the results are shown in Figure A .1 along with their corresponding independence scores. The CHF-based method recovers the original sources well, upto sign. The Kurtosis and CGF-based method fails to recover the second source. This is consistent with their higher independence score. The Meta algorithm selects CHF from candidates CHF, CGF, Kurtosis, FastICA, and JADE. 
global and local convergence: loss landscape of noisy ica: 
 Here we present sufficient conditions for global and local convergence of a broad class of contrast functions. This covers the cumulant-based contrast functions and our CHF and CGF-based proposals.The classical argument for global optima of kurtosis for noiseless ICA in [18] assumes WLOG that B is orthogonal. This reduces the optimization over u u u into one for v v v = Bu u u. Since B is orthogonal, ∥u u u∥= 1 translates into ∥v v v∥= 1. It may seem that this idea should extend to the noisy case due to property 1 a) and c) by optimizing over v v v = Bu u u over potentially non-orthogonal mixing matrices B.However, for non-orthogonal B, the norm constraint onv v v is no longer v v v T v v v = 1, rendering the original argument invalid.In what follows, we extend the original argument to the noisy ICA setting by introducing pseudo-euclidean constraints. Our framework includes a large family of contrast functions, including cumulants.To our knowledge, this is the first work to characterize the loss functions in the pseudo-euclidean space. In contrast, [50] provides global convergence of the cumulant-based methods by a convergence argument of the power method itself.Consider the contrast function f (u u u|P ) := E x x x∼P g x x x T u u u and recall the definition of the quasi orthogonalization matrix C = BDB T , where D is a diagonal matrix. For simplicity, WLOG let us assume that D is invertible. We now aim to find the optimization objective that leads to the pseudo-Euclidean update in Eq 3. For f C -1 u u u|P = E x x x∼P g u u u T C -1 x x x , consider the following:f C -1 u u u|P = i∈[k] E g B -1 u u u i z i /D ii = i∈[k] E [g (α i zi )] = i∈[k] f (α i /D ii |z i ) (6)where we define α α α := B -1 u u u and zi = z i /D ii . We now examine f (C -1 u u u) subject to the "pseudo" norm constraint u u u T C -1 u u u = 1. The key point is that for suitably defined f , one can construct a matrix C = BDB T from the data even when B is unknown. So our new objective is to optimizef C -1 u u u|P s.t. u u u T C -1 u u u = 1Using the Lagrange multiplier method, optimizing the above simply gives the power method update (Eq 3) u u u ∝ ∇f (C -1 u u u). Furthermore, optimizing Eq 6 leads to the following transformed objective:EQUATIONNow we provide our theorem about maximizing f (C -1 u u u|P ). Analogous results hold for minimization. We will denote the corresponding derivatives evaluated at t 0 as f ′ (t 0 |z i ) and f ′′ (t 0 |z i ). Theorem 3. Let C be a matrix of the form BDB T , whered i := D ii = f ′′ (u i |z i ) for some random u i . Let S + = {i : d i > 0}.Consider a contrast function f : R → R. Assume that for every non-Gaussian independent component, Assumption 1 holds and the third derivative, f ′′′ (u|X), does not change the sign in the half-lines [0, ∞), (∞, 0]. Then f C -1 u u u x x x with the constraint of ⟨u u u, u u u⟩ C -1 = 1 has local maxima at B -1 u u u = e e e i , i ∈ S + . All solutions satisfying e e e T i B -1 u u u ̸ = 0, ∀i ∈ [1, k] are minima, and all solutions with e e e T i B -1 u u u ̸ = 0 for i ∈ (1, k) are saddle points. These are the only fixed points.Note that the above result immediately implies that for C = -C, the same optimization in Theorem 3 will have maxima at all u u u such that B -1 u u u = e e e i , i ∈ S -. Corollary 3.1. 2k th cumulant-based contrast functions for k > 1, have maxima and minima at the directions of the columns of B, provided {z z z i } i∈[k] have a corresponding non-zero cumulant.Proof. This proof (see Appendix Section A.1) is immediate since cumulants satisfy (a), (c). For (c) and (d), note that f (u|Z) is simply u 2j κ 2j (Z), where κ 2j (Z) is the 2j th cumulant of Z.Theorem 3 can be easily extended to the case where C is not invertible. The condition on the third derivative, f ′′′ (u|X), may seem strong, and it is not hard to construct examples of random variables Z where this fails for suitable contrast functions (see [36] ). However, for such settings, it is hard to separate Z from a Gaussian. See the Appendix A.4.1 for more details.In this section, we establish a local convergence result for the power iteration-based method described in Eq 3. Define ∀t ∈ R, ∀i ∈ [d], the functions q i (t) := f ′ αi Dii z i . Then, without loss of generality, we have the following result for the first corner: [-∥B -1 ∥ 2 , ∥B -1 ∥ 2 ] and define c 1 , c 2 , c 3 > 0 such that ∀i ∈ [d], sup x∈B |qi(x)| c1 , |q ′ i (x)| c2 , |q ′′ i (x)| c3 ≤ 1. Define ϵ := ∥B∥ F ∥Be1∥ 2 max c3+c2∥Be1∥ |q1(α * 1 )| , c 2 2 |q1(α * 1 )| 2 , c3∥Be1∥ |q ′ 1 (α * 1 )| . Let ∥α α α 0 -α α α * ∥ 2 ≤ R, R ≤ max c2 c3 , 1 5ϵ . Then, we have ∀t ≥ 1, ∥α α α t+1 -α α α * ∥ ∥α α α t -α α α * ∥ ≤ 1 2 .Therefore, we can establish linear convergence without the third derivative constraint required for global convergence (see Theorem 3), by assuming some mild regularity conditions on the contrast functional and a sufficiently close initialization. The proof is included in Appendix Section A.1.2. Remark 3. Let α α α = B -1 u u u denote the demixed direction, u u u. Theorem 4 shows that if the initial u u u 0 is such that α α α 0 is close to the ground truth α α α * (which is WLOG taken to be e e e 1 ), then there is geometric convergence to α α α * . |q 1 (α * 1 )| and |q ′ 1 (α * 1 )| essentially quantify how non-gaussian the corresponding independent component is since they are zero for a mean zero Gaussian. When these are large quantities, ϵ is small, and the initialization radius ∥α α α 0 -α α α * ∥ 2 can be relatively larger. 
experiments: 
 In this section, we provide experiments to compare the fixed-point algorithms based on the characteristic function (CHF), the cumulant generating function (CGF) (Eqs 4 and 5) with the kurtosis-based algorithm (PEGI-κ 4 [49] ). We also compare against noiseless ICA 3 algorithms -FastICA, JADE, and PFICA. These six algorithms are included as candidates for the Meta algorithm (see Algorithm 1). We also present the Uncorrected Meta algorithm (Unc. Meta) to denote the Meta algorithm with the uncorrected independence score proposed in CHFICA [21] . Table 1 and Figure 2 provide experiments on simulated data, and Figure 3 provides an application for image demixing [16] . We provide additional experiments on denoising MNIST images in the Appendix Section A.2.Experimental setup: Similar to [49] , the mixing matrix B is constructed as B = U ΛV T , where U, V are k dimensional random orthonormal matrices, and Λ is a diagonal matrix with Λ ii ∈ [1, 3] . The covariance matrix Σ of the noise g g g follows the Wishart distribution and is chosen to be ρ k RR T , where k is the number of sources, and R is a random Gaussian matrix. Higher values of the noise power ρ can make the noisy ICA problem harder. Keeping B fixed, we report the median of 100 random runs of data generated from a given distribution (different for different experiments). The quasi-orthogonalization matrix for CHF and CGF is initialized as B BT using the mixing matrix, B, estimated via PFICA. The performance of CHF and CGF is based on a single random initialization of the vector in the power method (see Algorithm 2). Our experiments were performed on a Macbook Pro M2 2022 CPU with 8 GB RAM.Error Metrics: Due to the inherent ambiguity in signal recovery discussed in Section 2, we measure error using the accuracy of estimating B. We report the widely used Amari error [2, 24, 11, 10] for our results. For estimated demixing matrix B, define W = B -1 B, after normalizing the rows of B -1 and B -1 . Then we measure d B,B as -EQUATIONVariance reduction using Meta: We first show that the independence score can also be used to pick the best solution from many random initializations. In Figure 2 Effect of varying kurtosis: For our second experiment (see Table 1 ), we use k = 5 independent components, sample size n = 10 5 and noise power ρ = 0.2 from a Bernoulli(p) distribution, where we vary p from 0.001 to 0. κ 4 := (1 -6p (1 -p))/(p (1 -p)).Observe that the Meta algorithm (shaded in red) performs at par or better than the best candidate algorithms. FastICA did not converge for zero-kurtosis data.perform differently (best candidate algorithm is highlighted in bold font) for different p. In particular, characteristic function-based methods like PFICA and CHF perform poorly for small values of p. We attribute this to the fact that the characteristic function is close to one for small p. Kurtosis-based algorithms like PEGI, JADE, and FASTICA perform poorly for kurtosis close to zero. Furthermore, the Uncorrected Meta algorithm performs worse than the Meta algorithm since it shadows PFICA.Effect of varying noise power: For our next two experiments, we use k = 9, out of which 3 are from Uniform -√ 3, √ 3 , 3 are from Exponential(5) and 3 are from Bernoulli 1 2 -1 12and hence have zero kurtosis. In these experiments, we vary the sample size (Figure 2b ), n fixing ρ = 0.2, and the noise power (Figure 2a ) , ρ fixing n = 10 5 , respectively. Note that with most mixture distributions, it is easily possible to have low or zero kurtosis. We include such signals in our data to highlight some limitations of PEGI-κ 4 and show that Algorithm 1 can choose adaptively to obtain better results. Figure 2a shows that large noise power leads to worse performance for all methods. PEGI, JADE, and FastICA perform poorly consistently, because of some zero kurtosis components. CGF suffers because of heavy-tailed components. However, CHF and PFICA perform well consistently. The Meta algorithm mostly picks the best algorithm except for a few points, where the difference between the two leading algorithms is small. The Uncorrected Meta algorithm (which uses the independence score without the additional correction term introduced for noise) performs identically to PFICA.Effect of varying sample size: Figure 2b shows that the noiseless ICA methods have good performance at smaller sample sizes. However, they suffer a bias in performance compared to their noiseless counterparts as the sample size increases. The Meta algorithm can consistently pick the best option amongst the candidates, irrespective of the distribution, leading to significant improvements in performance. What is interesting is that up to a sample size of 10 4 , PFICA dominates other algorithms and Meta performs like PFICA. However, after that CHF dominates, and one can see that Meta starts to have a similar trend as CHF. We also see that the Uncorrected Meta algorithm (which uses the independence score without the additional correction term introduced for noise) has a near identical error as PFICA and has a bias for large n.Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material?Answer: [Yes] Justification: We make the code for our experiments available.Guidelines:• The answer NA means that paper does not include experiments requiring code.• Please see the NeurIPS code and data submission guidelines (https://nips.cc/ public/guides/CodeSubmissionPolicy) for more details. • While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for a new open-source benchmark). • The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https: //nips.cc/public/guides/CodeSubmissionPolicy) for more details. • The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc. • The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only a subset of experiments are reproducible, they should state which ones are omitted from the script and why. • At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable). • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
conclusion: 
 ICA is a classical problem that aims to extract independent non-Gaussian components. The vast literature on both noiseless and noisy ICA introduces many different inference methods based on a variety of contrast functions for separating the non-Gaussian signal from the Gaussian noise. Each has its own set of shortcomings. We aim to identify the best method for a given dataset in a data-driven fashion. In this paper, we propose a nonparametric score, which is used to evaluate the quality of the solution of any inference method for the noisy ICA model. Using this score, we design a Meta algorithm, which chooses among a set of candidate solutions of the demixing matrix. We also provide new contrast functions and a computationally efficient optimization framework. While they also have shortcomings, we show that our diagnostic can remedy them. We provide uniform convergence properties of our score and theoretical results for local and global convergence of our methods. Simulated and real-world experiments show that our Meta-algorithm matches the accuracy of the best candidate across various distributional settings.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort. 12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected? Answer: [Yes] Justification: MATLAB implementations (under the GNU General Public License) can be found at -FastICA and JADE. The code for PFICA was provided on request by the authors. Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 13. New Assets Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets? Answer: [NA] Justification: The paper does not release new assets to require documentation or licensing. Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: The Appendix is organized as follows -• Section A.1 proves Theorems 1, 2, A.1, 3 and 4 • Section A.2 provides additional experiments for noisy ICA • Section A.3 provides the algorithm to compute independence scores via a sequential procedure • Section A.4 explores the third derivative constraint in Theorem 3 and provides examples where it holds • Section A.5 contains surface plots of the loss landscape for noisy ICA using the CHF-based contrast functionsProof of Theorem 1. We will give an intuitive argument for the easy direction for k = 2. We will show that if F = DB -1 for some permutation of a diagonal matrix D, then after projecting using that matrix, the data is of the form z z z + g g g ′ for some Gaussian vector g g g ′ . Let k = 2. Let the entries of this vector be z 1 + g ′ 1 and z 2 + g ′ 2 , where z i , i ∈ {1, 2} are mean zero independent non-Gaussian random variables and g ′ i , i ∈ {1, 2} are mean zero possibly dependent Gaussian variables. The Gaussian variables are independent of the non-Gaussian random variables. Let t i , i ∈ {1, 2} be arbitrary but fixed real numbers. Let t t t = (t 1 , t 2 ) and let Σ g ′ denote the covariance matrix of the Gaussian g ′ . Assume for simplicity var(z i ) = 1 for i ∈ {1, 2}. Denote by Λ := cov(Fx x x) = F SF T = I + Σ g ′ .The JOINT part of the score is given by:EQUATIONThe PRODUCT part of the score is given byEQUATIONIt is not hard to see that JOINT equals PRODUCT. The same argument generalizes to k > 2. We next provide proof for the harder direction here. Suppose that ∆ F (t t t|P ) = 0, i.e, E exp(it t t T Fx x x) exp -t t t T diag F SF T t t t 2 - k j=1 E [exp(it j (Fx x x) j )] exp - t t t T F SF T t t t 2 = 0 (A.9)We then prove that F must be of the form DB -1 for D being a permutation of a diagonal matrix. Then, taking the logarithm, we haveln E exp(it t t T Fx x x) - k j=1 ln (E [exp(it j (Fx x x) j )]) = 1 2 t t t T diag F SF T -F SF T t t t = 1 2 t t t T diag F ΣF T -F ΣF T t t tUnder the ICA model we have, x x x = Bz z z + g g g. Therefore, from Eq A.9, using the definition of the characteristic function of a Gaussian random variable and noting that cov (g g g) = Σ, we haveln E exp(it t t T Fx x x) = ln E exp(it t t T F Bz z z) + ln E exp(it t t T Fg g g) = ln E exp(it t t T F Bz z z) - 1 2 t t t T F ΣF T t t tand similarly, To simplify the notation, let ψ j (t) := E e itzj , f j ≡ log ψ j , M := F B. The above is then equivalent tok j=1 ln (E [exp(it j (Fx x x) j )]) = k j=1 ln (E [exp(it j (F Bz z z) j )]) - 1 2 t t t T diag F ΣF T t2 k=1 M 1k M 2k f ′′ k (M 1k t 1 + M 2k t 2 ) ≡ const. (A.11)Note that M is invertible since F, B are invertible. Therefore, let s ≡ M T t, then2 k=1 M 1k M 2k f ′′ k (s k ) ≡ const.which implies that eitherM 1k M 2k = 0 or f ′′ k (s k ) ≡ const, for k = 1, 2.For any k ∈ {1, 2}, since z k is non-Gaussian, its logarithmic characteristic function, i.e. f k cannot be a quadratic function, soM 1k M 2k = 0which implies that each column of M has a zero. Since M is invertible, thus M is either a diagonal matrix or a permutation of a diagonal matrix. Now to consider k > 2, we just need to apply the above k = 2 argument to pairwise entries of {t 1 , • • • , t k } with other entries fixed. Hence proved.Proof of Theorem 3. We start with considering the following constrained optimization problem -sup u u u T C -1 u u u=1 f C -1 u u u|PBy an application of lagrange multipliers, for the optima u u u opt , we haveνC -1 u u u opt = C -1 ∇f (C -1 ) T u u u opt |P , u u u T opt C -1 u u u opt = 1 with multiplier ν ∈ R.This leads to a fixed-point iteration very similar to the one in [49] , which is the motivation for considering such an optimization framework. We now introduce some notations to simplify the problem.Let u u u = Bα α α, z ′ i := zi di and a ′ i := var(z ′ i ) = ai d 2 i, where a i := var(z i ). Then,sup u u u T C -1 u u u=1 f C -1 u u u|P = sup α α α T D -1 α α α=1 f D -1 α α α|z z zWe will use f (u u u|P ) or f (u u u|x x x) where x x x ∼ P interchangeably. By Assumption 1-(a), we see thatf D -1 α α α|z z z = k i=1 f αi di |z i .For simplicity of notation, we will useh i (α i /d i ) = f (α i /d i |z i ), since the functional form can be different for different random variables z i . Therefore, we seek to characterize the fixed points of the following optimization problem -sup α α α k i=1 h i (α i /d i ) (A.12) s.t k i=1 α 2 i d i = 1, d i ̸ = 0 ∀ i ∈ [k]where α α α, d ∈ R k . We have essentially moved from optimizing in the ⟨., .⟩ C -1 space to the ⟨., .⟩ D -1 space. We find stationary points of the LagrangianL(α α α, λ) := k i=1 h i α i d i -λ k i=1 α 2 i d i -1 (A.13)The components of the gradient of L(α α α, λ) w.r.t α α α are given as∂ ∂α j L(α α α, λ) = 1 d j h ′ j α j d j -λ α j d j ,At the fixed point (α α α, λ) we have,∀j ∈ [k], h ′ j α j d j -λα j = 0By Assumption 1-(b), for α i = 0, the above is automatically zero. However, for {j : α j ̸ = 0}, we have,λ = h ′ j αj dj α j (A.14)So, for all {j : α j ̸ = 0}, we must have the same value and sign of 1 αj h ′ j αj dj . By assumption in the theorem statement, h ′′′ i (x) does not change sign on the half-lines x > 0 and x < 0. Along with Assumption 1-(b), this implies that∀ x ∈ [0, ∞), sgn(h i (x)) = sgn(h ′ i (x)) = sgn(h ′′ i (x)) = sgn(h ′′′ i (x)) = κ 1 (A.15) ∀ x ∈ (-∞, 0], sgn(h i (x)) = -sgn(h ′ i (x)) = sgn(h ′′ i (x)) = -sgn(h ′′′ i (x)) = κ 2 (A.16) where κ 1 , κ 2 ∈ {-1, 1} are constants. Furthermore, we note that Assumption 1-(d) ensures that h i (x) is a symmetric function. Therefore, ∀ x ∈ R, sgn(h i (x)) = sgn(h ′′ i (x)) = κ, κ ∈ {-1, 1}. Then, since d i = h ′′ i (u i ), sgn (d i ) = sgn (h i (x)) , ∀x ∈ R (A.17)Now, using A.14,sgn(λ) = sgn h ′ j α j d j × sgn(α j ) = sgn α j d j × sgn h j α j d j× sgn(α j ), using Eq A.15 and A.16= sgn (d j ) × sgn h j α j d j = 1, using Eq A.17(A.18)Keeping this in mind, we now compute the Hessian, H ∈ R k×k , of the lagrangian, L(α α α, λ) at the fixed point, (α α α, λ). Recall that, we have h ′ i (0) = 0 and h ′′ i (0) = 0. Thus, for {i :α i = 0}, H ii = -λ/d i (A.19)This implies that sgn (d i H ii ) = sgn (λ) = -1 for {i : α i = 0}, using Eq A.18.For {i : α i ̸ = 0}, we haveH ij = ∂ 2 ∂α i ∂α j L(α α α, λ) α α α = 1(i = j)   h ′′ i αi di d 2 i - λ d i   = 1(i = j)   h ′′ i αi di d 2 i - h ′ i αi di α i d i   , for α i ̸ = 0, using A.14 = 1(i = j) 1 d i α i α i d i h ′′ i α i d i -h ′ i α i d i , for α i ̸ = 0We consider the pseudo inner product space ⟨., .⟩ D -1 for optimizing α α α. Furthermore, since we are in this pseudo-inner product space, we have⟨v v v, Hv v v⟩ D -1 = v v vD -1 Hv v vSo we will consider the positive definite-ness of the matrix H := D -1 H to characterize the fixed points. Recall that for a differentiable convex function f , we have ∀ x, y ∈ dom(f) ⊆ R n f (y) ≥ f (x) + ∇f (x) T (y -x)Therefore, for {i :α i ̸ = 0}, we have sgn (d i H ii ) = sgn (d i ) × sgn (d i α i ) × sgn α i d i h ′′ i α i d i -h ′ i α i d i = sgn (α i ) × sgn h ′′′ i α i d i , using convexity/concavity of h ′ (.) = sgn (α i ) × sgn h ′ i α i d i, using A.15 and A.16= sgn (α i ) × sgn (λ) × sgn (α i ) , using A.14 = sgn (λ) = 1, using Eq A.17 (A.20)Let S := {i : α i ̸ = 0}. We then have the following cases -1. Assume that only one α i is nonzero, then ∀v v v orthogonal to this direction, ⟨v v v, Hv v v⟩ D -1 < 0 by Eq A. 19 . Thus this gives a local maxima.2. Assume more than one α i are nonzero, but |S|< k. Then we have for i ̸ ∈ S, Hjj < 0 using A. 19 . For i ∈ S, Hii > 0 from Eq A. 20 . Hence these are saddle points.3. Assume we have S = [k], i.e. ∀i, α i ̸ = 0. In this case, H ≻ 0. So, we have a local minima.This completes our proof.Theorem A.1. Consider the data generated from the noisy ICA model (Eq 1). If f (u u u|P ) is defined as in Eq 4 or Eq 5, we have ∇ 2 f (u u u|P ) = BD u B T , for some diagonal matrix D u , which can be different for the differerent contrast functions. Remark 4. The above theorem is useful because it shows that the Hessian of the contrast functions based on the CHF (eq 4 or CGF (eq 5) is of the form BDB T , where D is some diagonal matrix. This matrix can be precomputed at some vector u u u and used as the matrix C in the power iteration update (see eq 3) in Algorithm 2 of [49] .Proof of Theorem A.1. First, we show that the CHF-based contrast function (see Eq 4), the Hessian, is of the correct form.f (u u u|P ) = log E exp(iu u u T Bz z z) + log E exp(-iu u u T Bz z z) + u u u T B diag(a a a)B T u u u,where diag(a a a) is the covariance matrix of z z z. For simplicity of notation, define ϕ(z z z; u u u):= E exp(iu u u T Bz z z) . ∇f (u u u|P ) = iB      E exp(iu u u T Bz z z)z z z ϕ(z z z; u u u) µ(u u u;z z z) -i E exp(-iu u u T z z z)z z z ϕ(z z z; -u u u)      + 2B diag(a a a)B T u u u Define H(u u u; z z z) := E exp(iu u u T Bz z z)z z zz z z T ϕ(z z z; u u u) -µ(u u u; z z z)µ(u u u; z z z) T .Taking a second derivative, we have:∇ 2 f (u u u|P ) = B (-H(u u u; z z z) -H(z z z; -u u u) + 2 diag(a a a)) B TAll we have to show at this point is that H(u u u; z z z) is diagonal. We will evaluate the k, ℓ entry. Let B i denote the i th column of B. Let Y \k,ℓ := u u u T B j̸ =k,ℓ z j . Using independence of the components of z z z, we have, for k ̸ = ℓ,E exp(iu u u T B k z k + iu u u T B ℓ z ℓ + Y \kℓ )z k z ℓ ϕ(z z z; u u u) = E exp(iu u u T B k z k )z k E exp(iu u u T B ℓ z ℓ )z ℓ E [exp(iu u u T B k z k )] E [exp(iu u u T B ℓ z ℓ )] (A.21)Now we evaluate:e e e T k µ(u u u; z z z) :=E exp(iu u u T Bz z z)z k ϕ(z z z; u u u) = E exp(iu u u T B k z k )z k E [exp(iu u u T B k z k )] (A.22)Using Eqs A.21 and A.22, we see that indeed H(u u u; z z z) is diagonal. Now, we prove the statement about the CGF-based contrast function.We have,EQUATIONCovariance of z z z ∼ P(z z z).exp u u u T Bz z zE [exp (u u u T Bz z z)] -diag (a a a)             B TThe new probability density P (z z z) . exp(u u u T Bz z z) E[exp(u u u T Bz z z)] is an exponential tilt of the original pdf, and since {z i } d i=1 are independent, and the new tilted density also factorizes over the z i 's, therefore, the covariance under this tilted density is also diagonal.In this section, we provide the proof for Theorem 2. First, we state some preliminary results about the uniform convergence of smooth function classes which will be useful for our proofs.Let D := sup θ, θ∈T ρ X θ, θ denote the diameter of set T, and let N X (δ; T) denote the δ-covering number of T in the ρ X metric. Then, we have the following standard result: Proposition A.2. (Proposition 5.17 from [51] ) Let {X θ , θ ∈ T} be a zero-mean subgaussian process with respect to the metric ρ X . Then for any δ ∈ [0, D] such that N X (δ; T) ≥ 10, we have E sup θ, θ∈T X θ -X θ ≤ 2E sup γ,γ ′ ∈T;ρ X (γ,γ ′ )≤δ X θ -X θ + 4 D 2 log (N X (δ; T))sup f ∈C 1 n n i=1 f X (i) -E [f (X)] ≤ 2 1 n E X,ϵ sup f ∈C n i=1 ϵ i f X (i) + δwith probability at least Proof. Define E = {t t t : ∥t t t∥≤ √ 2k log n}. We will use the fact that P (E c ) ≤ 2/n (Example 2.12, [51] ). Also note that by construction, ∆(t t t, F |P ) ≤ 1. Observe that:sup F ∈F |E t t t∈N (0,I k ) ∆(t t t, F |P ) -E t t t∈N (0,I k ) ∆(t t t, F | P )| ≤ sup F ∈F E t t t∈N (0,I k ) |∆(t t t, F |P ) -∆(t t t, F | P )| ≤ E t t t∈N (0,I k ) sup F ∈F |∆(t t t, F |P ) -∆(t t t, F | P )| ≤ E t t t∈N (0,I k ) sup F ∈F ∆(t t t, F |P ) -∆(t t t, F | P ) E + E t t t∈N (0,I k ) sup F ∈F |∆(t t t, F |P ) -∆(t t t, F | P ) |E c P (E c ) ≤ E t t t∈N (0,I k ) sup F ∈F ∆(t t t, F |P ) -∆(t t t, F | P ) E + 2/n = O P   k 2 ∥S∥max(k, σ 4 ∥S∥) log 2 (nC k ) n   + 2/n, using Theorem A.5The second inequality follows from Jensen's inequality, the fourth follows from E[sup(.)] ≤ sup(E[.]).Theorem A.5. Let F := {F ∈ R k×k : ∥F ∥≤ 1}. Assume that x x x ∼ subgaussian(σ). We have: We next have that, with probability at leastsup F ∈F |∆(t1 -1/n, diag F (S -S)F T 2 ≤ F (S -S)F T 2 ≤ C σ 2 ∥S∥ k log n nThus with probability at least 1 -1/n, using the inequality |1 -e x |≤ 2|x|, ∀x ∈ [-1, 1], Eq A.25 leads to:ϕ(t t t, F | P ) exp -t t t T diag(F SF T )t t t -exp -t t t T diag(F SF T )t t t ≤ 2C∥t t t∥ 2 σ 2 ∥S∥ k log n nNote that the matrices diag F (S -S)F T and F SF T are positive semi-definite and t t t is unit-norm. Observe that, using Lemmas A.6 and A. Proof. Let ϕ(t t t, F ; x x x (i) ) = exp(it t t T Fx x x (i) ). Next, we note that ∥F ∥ 2 ≤ 1 imply F T t t t ≤ ∥t t t∥. Therefore,sup F ∈F 1 n n j=1 exp(it t t T Fx x x (j) ) -E x x x exp(it t t T Fx x x) ≤ sup u u u∈R k ,∥u u u∥≤1 1 n n j=1 exp(iu u u T x x x (j) ) -E x x x exp(iu u u T x x x) Let f X (u u u) = i ϵ i ξ(u u u; x x x (i)). We will argue for unit vectors u u u, and later scale them by ∥t t t∥. DefineB k := u u u ∈ R k , ∥u u u∥ ≤ 1 , ξ(u u u; x x x) := exp(iu u u T x x x) and let d(u u u, u u u ′ ) 2 := i (ξ(u u u; x x x (i) ) -ξ(u u u ′ ; x x x (i) )) 2Then we have,∥∇ u u u exp(iu u u T x x x)∥ 2 = ∥x x x∥|-sin(u u u T x x x) + i cos(u u u T x x x)| ≤ ∥x x x∥Then, ξ(u u u; x x x (i) ) -ξ(t t t; x x x (i) ) ≤ min ∥∇ u u u ′′ ξ(u u u; x x x (i) )∥ 2 ∥u u u -u u u ′ ∥ , 2EQUATIONLet S := i x x x (i) x x x (i) T /n and τ n := n Tr S . We next have,D 2 := max u u u,u u u ′ d(u u u, u u u ′ ) 2 ≤ min 4n, max u u u,u u u ′ i ∥x x x (i) ∥ 2 2 ∥u u u -u u u ′ ∥ 2 ≤ 4 min {n, τ n } Next, we bound the covering number N (δ; B k , d u u u ). Note that N (δ; B k , d u u u ) ≤ N (δ/ √ τ n ; B k , ∥.∥ 2 )since,d (u u u, u u u ′ ) 2 = i (ξ(u u u; x x x (i) ) -ξ(u u u ′ ; x x x (i) )) 2 ≤ i ∥x x x (i) ∥ 2 2 ∥u u u -u u u ′ ∥ 2 = i Tr x x x (i) x x x (i) T ∥u u u -u u u ′ ∥ 2Using Cauchy-Schwarz, we have:|f X (u u u) -f X (u u u ′ )|≤ n i=1 ϵ i |ξ(u u u; x x x (i) ) -ξ(u u u; x x x (i) )|≤ nd(u u u, u u u ′ ) 2Therefore, we haveN (δ; B k , d u u u ) ≤ 3 √ τn δ k. Therefore, using Proposition A.2, Now we recall that ∥F T t t t∥≤ ∥t t t∥, since all vectors t t t appear as t t t T F T x x x (i) , this simply leads to scaling τ and ∥S∥ by ∥t t t∥ 2 . Dividing both sides by n, we have E x x x,ϵ sup F ∈F |f X (u u u) -f X (u u u ′ )| ≤ 2E x x x sup d(u u u,u u u ′ )≤δ |f X (u u u) -f X (u u u ′ )| + 4E x x x D log N (δ; B k , d u u u ) ≤ 2E x x x inf δ>0 δ √ n + 8 min {n, τ n } 2k log 3 √ τ n δ = O E x x x (k min(n, τ n ) log (max(τ n , n)) ,1 n E x x x,ϵ sup F ∈F |f X (t t t) -f X (t t t ′ )| ≤ k n O min(1,sup F ∈F ψ(t t t, F | P ) -ψ(t t t, F |P ) ≤ O P ∥t t t∥ k 2 ∥S∥log n n Proof. Let B k := u u u ∈ R k , ∥u u u∥ ≤ 1 Now define ψ j (t t t, F ; x x x) = E[exp(it j F T j x x x)]. Also define f (j) X (F ) = 1 n i ϵ i ψ j (t t t, F ; x x x (i)). Thus, using the same argument as in Lemma A.6 for t t t ← ∥t t t∥e e e j and x x x (i) ← t j x x x (i) , supF ∈F |ψ j (t t t, F | P ) -ψ j (t t t, F |P )|= O P |t j | k∥S∥log n nFinally, we see that:sup F ∈F ψ(t t t, F | P ) -ψ(t t t, F |P ) ≤ j O P |t j | k∥S∥log n n = O P k∥S∥log n n   k j=1 |t j |   = ∥t t t∥O P k 2 ∥S∥log n nThe above is true because, for|a i |, |b i |≤ 1, i = 1, . . . , k, k i=1 a i - k i=1 b i = k-1 j=0 i≤j b i (a j+1 -b j+1 ) k i=j+2 a i ≤ k-1 j=0 |a j+1 -b j+1 | A.1.Consider the function g i (y) := q i (y i ) ∥Bq (y y y) ∥ for y y y ∈ R d . We start by computing the gradient, ∇ y g i (y).Lemma A.8. Let g i (y) := q i (y i ) ∥Bq (y y y) ∥ , then we have[∇ y g i (y)] j =            1 ∥Bq (y y y) ∥ ∂q i (y i ) ∂y i - q i (y i ) ∥Bq (y y y) ∥ 2 ∂∥Bq (y y y) ∥ ∂q i (y i )∂q i (y i ) ∂y i , for j = i q ′ j (y j ) ∥Bq (y y y) ∥ q i (y i ) e e e T j B T Bq (y y y) ∥Bq (y y y) ∥ 2 , for j ̸ = iProof. The derivative w.r.t y i is given as -∂g i (y y y) ∂y i = 1 ∥Bq (y y y) ∥ ∂q i (y i ) ∂y i - q i (y i ) ∥Bq (y y y) ∥ 2 ∂∥Bq (y y y) ∥ ∂y i = 1 ∥Bq (y y y) ∥ ∂q i (y i ) ∂y i - q i (y i ) ∥Bq (y y y) ∥ 2 ∂∥Bq (y y y) ∥ ∂q i (y i ) ∂q i (y i ) ∂y i∇ y y y ∥Ay y y∥ 2 = 1 ∥Ay y y∥ A T Ay y y Therefore,∂g i (y y y) ∂y i = 1 ∥Bq (y y y) ∥ ∂q i (y i ) ∂y i - q i (y i ) ∥Bq (y y y) ∥ 2 1 ∥Bq (y y y) ∥e e e T i B T Bq (y y y) Lemma A.11. ∀ y y y ∈ R d satisfying (A.34), we have∂q i (y i ) ∂y i = q ′ i (y i ) ∥Bq (y y y) ∥ 1 - q i (y i ) e e1. For i = 1, 1 -5ϵ∥y y y -α α α * ∥≤ q 1 (α * 1 ) v (y y y|1) ≤ 1 + 5ϵ∥y y y -α α α * ∥ 2. For i ̸ = 1, |q 1 (α * 1 ) v (y y y; i)| ≤ (1 + 5ϵ∥y y y -α α α * ∥)∥Be e e i ∥ ∥Be e e 1 ∥Proof. Define θ := ϵ∥y y y -α α α * ∥ for convenience of notation. We have,e e e T i B T Bq (y y y) ∥Bq (y y y)∥ 2 = q 1 (α * 1 )e e e T i B T Be e e 1 ∥Bq (y y y)∥ 2 +e e e T i B T Bδ (y y y) ∥Bq (y y y) ∥ 2 , using definition of δ (y y y) (A.40)Therefore, if i = 1, then using Lemma A.10,1 (1 + θ) 2 + q 1 (α * 1 )e e e T 1 B T Bδ (y y y) ∥Bq (y y y) ∥ 2 ≤ q 1 (α * 1 ) e e e T 1 B T Bq (y y y) ∥Bq (y y y)∥ 2 ≤ 1 (1 -θ) 2 + q 1 (α * 1 )e e e T 1 B T Bδ (y y y) ∥Bq (y y y) ∥ 2 (A.41)Using the following inequalities -EQUATION) 2 ≤ 1 + 5θ, θ ≤ 1 5and observing that using Eq A.34, and Lemma A.10,q 1 (α * 1 ) e e e T 1 B T Bδ (y y y) ∥Bq (y y y) ∥ 2 ≤ |q 1 (α * 1 )| ∥Be e e 1 ∥ (1 -θ) 2 |q 1 (α * 1 )| 2 ∥Be e e 1 ∥ 2 ∥B∥ ∥δ (y y y)∥ ≤ 1 (1 -θ) 2 c 2 ∥B∥ |q 1 (α * 1 )| ∥Be e e 1 ∥ ∥y y y -α α α * ∥ ≤ ϵ (1 -θ) 2 ∥y y y -α α α * ∥Therefore we have from Eq A.41,1 -5θ ≤ q 1 (α * 1 ) v (y y y; i) ≤ 1 + 5θ (A.42)where we used θ := ϵ∥y y y -α α α * ∥. For the case of i ̸ = 1, using Lemma A.10 we have We now operate under the assumption that Eq A.34 holds for y y y = α α α t and inductively show that it holds for y y y = α α α t+1 as well. Recall the function g i (y y y) := q i (y i ) ∥Bq (y y y) ∥. By applying the mean-value theorem for g i (.) for the points α α α t and α α α * , we have from Eq A.26 and A.27 -|(α α α t+1 ) i -α * i | = |g i (α α α t ) -g i (α α α * )| = ∇g i (β β β i ) T (α α α t -α α α * ) for β β β i := (1 -λ i ) α α α t + λ i α α α * , λ i ∈ (0, 1) ≤ ∥∇g i (β β β i ) ∥∥α α α t -α α α * ∥ (A.44)Note the induction hypothesis assumes that Eq A.34 is true for x = α α α t . Since ∀i, λ i ∈ (0, 1), therefore Eq A.34 holds for all β β β i as well. Squaring and adding Eq A.44 for i ∈ [d] and taking a square-root, we have∥α α α t+1 -α α α * ∥ ≤ ∥α α α t -α α α * ∥ k i=1 ∥∇g i (β β β i ) ∥ 2 ≤ ∥α α α t -α α α * ∥ k i=1 k j=1   ∂g i (y y y) ∂y j β β β i   2 (A.45)Let us consider the expression G ij := ∂g i (y y y)∂y j β β β i , ∀i, j ∈ [k].For the purpose of the subsequent analysis, we define θ := ϵ ∥α α α t -α α α * ∥. We divide the analysis into the following cases -Case 1 : i = 1, j = 1Using Lemma A.8,|G 11 | = q ′ 1 ((β β β 1 ) 1 ) ∥Bw (β β β 1 ) ∥ 1 - q 1 ((β β β 1 ) 1 )e e e T 1 B T Bq q q (β β β 1 ) ∥Bq q q (β β β 1 ) ∥ 2 From Lemma A.9,|q ′ 1 ((β β β 1 ) 1 )| ≤ |q ′ 1 (α * 1 )| + c 3 |(β β β 1 ) 1 -α * 1 | = |q ′ 1 (α * 1 )| 1 + c 3 |(β β β 1 ) 1 -α * 1 | |q ′ 1 (α * 1 )| ≤ |q ′ 1 (α * 1 )| 1 + c 3 ∥α α α t -α α α * ∥ |q ′ 1 (α * 1 )| ≤ |q ′ 1 (α * 1 )| (1 + θ) and, q 1 ((β β β 1 ) 1 ) ≤ |q 1 (α * 1 ) |+c 2 |(β β β 1 ) 1 -α * 1 | = |q 1 (α * 1 ) | 1 + c 2 |(β β β 1 ) 1 -α * 1 | |q 1 (α * 1 ) ≤ |q 1 (α * 1 ) | 1 + c 2 ∥α α α t -α α α * ∥ |q 1 (α * 1 ) | ≤ |q 1 (α * 1 ) |(1 + θ) From Lemma A.10, ∥Bq q q (β β β 1 ) ∥ ≥ (1 -θ) |q 1 (α * 1 )| ∥Be e e 1 ∥ From Lemma A.11 and θ ≤ 1 10 , -6θ ≤ 1 - q 1 ((β β β 1 ) 1 ) e e e T 1 B T Bw (β β β 1 ) ∥Bw (β β β 1 ) ∥ 2 ≤ 6θTherefore,|G 11 | ≤ 6 1 + θ 1 -θ |q ′ 1 (α * 1 )| θ ∥Be e e 1 ∥|q 1 (α * 1 )| ≤ 7.5ϵ |q ′ 1 (α * 1 )| ∥Be e e 1 ∥|q 1 (α * 1 )| ∥α α α t -α α α * ∥ , since θ ≤ 1 10Case 2 : i = 1, j ̸ = 1Using Lemma A.8,|G 1j | = q ′ j (β β β 1 ) j ∥Bq q q (β β β 1 ) ∥ q 1 ((β β β 1 ) 1 )e e e T j B T Bq q q (β β β 1 ) ∥Bq q q (β β β 1 ) ∥ 2 From Lemma A.9,q ′ j (β β β 1 ) j ≤ c 3 (β β β 1 ) j -α * j ≤ c 3 (α α α t ) j -α * j , since α * j = 0From Lemma A.10,∥Bq q q (β β β 1 ) ∥≥ (1 -θ) |q 1 (α * 1 )| ∥Be e e 1 ∥ From Lemma A.9,|q 1 ((β β β 1 ) 1 )| ≤ |q 1 (α * 1 )| + c 2 |(β β β 1 ) 1 -α * 1 | = |q 1 (α * 1 )| 1 + c 2 |(β β β 1 ) 1 -α * 1 | |q 1 (α * 1 )| ≤ |q 1 (α * 1 )| (1 + θ) From Lemma A.11,e e e T j B T Bq q q (β β β 1 ) ∥Bq q q (β β β 1 )∥ 2 ≤ (1 + 5θ) |q 1 (α * 1 ) |∥Be e e j ∥ ∥Be e e 1 ∥ Therefore,|G 1j | ≤ c 3 1 + θ 1 -θ (1 + 5θ) ∥Be e e j ∥ |q 1 (α * 1 ) |∥Be e e 1 ∥ 2 (α α α t ) j -α * j ≤ 2c 3 ∥Be e e j ∥ |q 1 (α * 1 ) |∥Be e e 1 ∥ 2 (α α α t ) j -α * j Case 3 : i ̸ = 1, j = 1Using Lemma A.8,|G i1 | = q ′ 1 ((β β β i ) 1 ) ∥Bq q q (β β β i ) ∥ q i ((β β β i ) i ) e e e T 1 B T Bq q q (β β β i ) ∥Bq q q (β β β i ) ∥ 2From Lemma A.9,|q ′ 1 ((β β β i ) 1 )| ≤ |q ′ 1 (α * 1 )| + c 3 |(β β β i ) 1 -α * 1 | = |q ′ 1 (α * 1 )| 1 + c 3 |(β β β i ) 1 -α * 1 | |q ′ 1 (α * 1 )| ≤ |q ′ 1 (α * 1 )| (1 + θ) From Lemma A.10, ∥Bq q q (β β β i ) ∥≥ (1 -θ) |q 1 (α * 1 )| ∥Be e e 1 ∥ From Lemma A.9,|q i ((β β β i ) i )| ≤ c 3 2 ((β β β i ) i -α * i ) 2 ≤ c 3 2 ((α α α t ) i -α * i ) 2 From Lemma A.11, e e e T 1 B T Bq q q (β β β i ) ∥Bq q q (β β β i ) ∥ 2 ≤ 1 + 5θ |q 1 (α * 1 )| Therefore, |G i1 | ≤ c 3 2 (1 + θ) (1 + 5θ) 1 -θ |q ′ 1 (α * 1 )| |q 1 (α * 1 )| 2 ∥Be e e 1 ∥ ((α α α t ) i -α * i ) 2 ≤ c 3 |q ′ 1 (α * 1 )| |q 1 (α * 1 )| 2 ∥Be e e 1 ∥ ((α α α t ) i -α * i ) 2 ≤ c 3 R |q ′ 1 (α * 1 )| |q 1 (α * 1 )| 2 ∥Be e e 1 ∥ |(α α α t ) i -α * i | ≤ c 2 |q ′ 1 (α * 1 )| |q 1 (α * 1 )| 2 ∥Be e e 1 ∥ |(α α α t ) i -α * i | Case 4 : i ̸ = 1, j ̸ = 1, i = jUsing Lemma A.8,|G ii | = q ′ i ((β β β i ) i ) ∥Bq q q (β β β i ) ∥ 1 - q i ((β β β i ) i ) e e e T i B T Bq q q (β β β i ) ∥Bq q q (β β β i ) ∥ 2From Lemma A.9,|q ′ i ((β β β i ) i )| ≤ c 3 (β β β i ) i -α * j ≤ c 3 |(α α α t ) i -α * i | From Lemma A.10, ∥Bq q q (β β β i ) ∥≥ (1 -θ) |q 1 (α * 1 )| ∥Be e e 1 ∥ From Lemma A.9,|q i ((β β β i ) i )| ≤ c 3 2 ((β β β i ) i -α * i ) 2 ≤ c 3 2 ((α α α t ) i -α * i ) 2 ≤ c 3 R ∥α α α t -α α α * ∥ ≤ c 2 ∥α α α t -α α α * ∥ From Lemma A.11,e e e T i B T Bq q q (β β β i )∥Bq q q (β β β i ) ∥ 2 ≤ (1 + 5θ)∥Be e e i ∥ |q 1 (α * 1 ) |∥Be e e 1 ∥ Then,q i ((β β β i ) i ) e e e T i B T Bq q q (β β β i ) ∥Bq q q (β β β i ) ∥ 2 ≤ (1 + 5θ) ∥α α α t -α α α * ∥ ∥Be e e i ∥ ∥Be e e 1 ∥ c 2 |q 1 (α * 1 ) | ≤ (1 + 5θ) ∥α α α t -α α α * ∥ ∥B∥ F ∥Be e e 1 ∥ c 2 |q 1 (α * 1 ) | ≤ (1 + 5θ) ϵ ∥α α α t -α α α * ∥ ≤ 2θTherefore,|G ii | ≤ c 3 |(α α α t ) i -α * i | (1 -θ) |q 1 (α * 1 )| ∥Be e e 1 ∥ ≤ 2c 3 |q 1 (α * 1 )| ∥Be e e 1 ∥ |(α α α t ) i -α * i | Case 5 : i ̸ = 1, j ̸ = 1, i ̸ = jUsing Lemma A.8,|G ij | = q ′ j (β β β i ) j ∥Bq q q ((β β β i )) ∥ q i ((β β β i ) i) e e e T j B T Bq q q ((β β β i )) ∥Bq q q ((ββ β i )) ∥ 2From Lemma A.9,q ′ j (β β β i ) j ≤ c 3 (β β β i ) j -α * j ≤ c 3 (α α α t ) j -α * jFrom Lemma A.10,∥Bq q q (β β β i ) ∥≥ (1 -θ) |q 1 (α * 1 )| ∥Be e e 1 ∥ From Lemma A.9,|q i ((β β β i ) i )| ≤ c 3 2 ((β β β i ) i -α * i ) 2 ≤ c 3 2 ((α α α t ) i -α * i ) 2 ≤ c 3 R 2 |(α α α t ) i -α * i | From Lemma A.11,e e e T j B T Bq q q (β β β i ) ∥Bq q q (β β β i )∥ 2 ≤ (1 + 5θ) |q 1 (α * 1 ) |∥Be e e j ∥ ∥Be e e 1 ∥ Therefore,|G ij | ≤ c 2 3 R 2 1 + 5θ 1 -θ ∥Be e e j ∥ |q 1 (α * 1 )| 2 ∥Be e e 1 ∥ 2 (α α α t ) j -α * j |(α α α t ) i -α * i | ≤ c 2 3 R ∥Be e e j ∥ |q 1 (α * 1 )| 2 ∥Be e e 1 ∥ 2 (α α α t ) j -α * j |(α α α t ) i -α * i | A.1.2.In this experiment, we use the ICA-based denoising technique proposed in [38] to compare candidate Noisy ICA algorithms and show that the Meta algorithm can pick the best-denoised image based on the independence score proposed in our work.We use the noisy MNIST dataset and further add entrywise Gaussian noise with variance proportional to sin 2 (c 1 π(i + j))(c 1 = 50) to the (i, j) th pixel. Training images are flattened to create a 784dimensional vector, and PCA is performed to reduce dimensionality to 25, on which subsequently ICA is performed. The original and denoised images along with their independence score are shown in Figure A.2. We note that CHF-based denoising provides qualitatively better results, while CGF-based denoising provides the worst results. This is consistent with their corresponding independence scores.In this section, we provide an algorithm for the computation of the independence score proposed in the manuscript, but in a sequential manner. The detailed algorithm is provided as Algorithm 3. We assume that we have access to a matrix of the form C = BDB T .The power method (see Eq 3) essentially extracts one column of B (up to scaling) at every step. [49] provides an elegant way to use the pseudo-Euclidean space to successively project the data onto columns orthogonal to the ones extracted. This enables one to extract each column of the mixing matrix. For completeness, we present the steps of this projection. After extracting the i th (f) Demixed Images using FastICA (∆(t t t, F |P ) = 7.1 × 10 -3 ).(g) Demixed Images using JADE (∆(t t t, F |P ) = 6.8 × 10 -3 ). where g g g = Fg g g, where F is some k × k matrix. So we have min(ℓ, k) vectors of the form z i B i + g g g i , and when ℓ < k an additional vector which contains all independent random variables z j , j > ℓ along with a mean-zero Gaussian vector. Then we can project each vector to a scalar using unit random vectors and check for independence using the Independence Score defined in 2. When ℓ = k, then all we need are the j = 1, . . . , k projections on the k directions identified via ICA in Eq A. 47 .As an example, we conduct an experiment (Figure A.3) , where we fix a mixing matrix B using the same generating mechanism in Section 5. Now we create vectors q which interpolate between B(:, 1) and a fixed arbitrary vector orthogonal to B(:, 1). As the interpolation changes, we plot the independence score for direction q. The dataset is the 9-dimensional dataset, which has independent components from many different distributions (see Section 5). This plot clearly shows that there is a clear negative correlation between the score and the dot product of a vector with the column B(:, 1). To be concrete, when q has a small angle with B(:, 1), the independence score is small, and the error bars are also very small. However, as the dot product decreases, the score grows and the error bars become larger.Algorithm 3 Independence Score after extracting ℓ columns. U and V are running estimates of B and B -1 upto ℓ columns and rows respectively. Input X ∈ R n×k , U ∈ R k×ℓ , V ∈ R ℓ×k , Number of random projections M k 0 ← min(ℓ + 1, k) for j in range[1, ℓ] do Y j ← X (U (:, j)V (j, :)) T end for if ℓ < k then Y ℓ+1 ← X I -V T U T end if for j in range[1, M ] do t t t ← random unit vector in R k for a in range[1, k 0 ] do W (:, j) ← Y j t t t end for Let α α α ∈ R 1,k0 represent a row of W S ← cov(W ) γ ← i,j Sij β β β ← i Sii s(j) = | Ê exp( j iα j ) exp(-β β β) -Consider the probability distribution function (see [36] ) f (x) = ke -x 2 2 (a + b cos(cπx)) for constants a, b, c > 0. For f (x) to be a valid pdf, we require that f (x) ≥ 0, and∞ -∞ f (x) dx = 1. Therefore, we have, ∞ -∞ ke -x 2 2 (a + b cos(cπx)) dx = 1 (A.48)Noting standard integration results, we have that,∞ -∞ e -x 2 2 dx = √ 2π and ∞ -∞ e -x 2 2 cos(cπx) dx = √ 2πe -c 2 π 2 2 (A.49) Therefore, k = 1 √ 2π(a+be -c 2 π 2 2 ). Some algebraic manipulations yield the following:M GF f (x) (t) = 1 (a + be -c 2 π 2 2 ) e t 2 2 a + be -c 2 π 2 2 cos(cπt) Therefore, ln(M GF f (x) (t)) = -ln(a + be -c 2 π 2 2 ) + t 2 2 + ln(a + be -c 2 π 2 2 cos(cπt)).We now compute the mean, µ, and variance σ 2 .The mean, µ can be given as :µ = ∞ -∞ f (x)x dx = 0 since f (x) is an even function.The variance, σ 2 can therefore be given as :σ 2 = k a ∞ -∞ e -x 2 2 x 2 dx + b ∞ -∞ e -x 2 2 x 2 cos(cπx) dx Now, we note that ∞ -∞ e -x 2 2 x 2 dx = √ 2π and ∞ -∞ e -x 2 2 x 2 cos(cπx) dx = e -c 2 π 2 2 √ 2π(1 -c 2 π 2 ) Therefore, σ 2 = k a √ 2π + be -c 2 π 2 2 √ 2π(1 -c 2 π 2 ) = 1 - bc 2 π 2 e -c 2 π 2 2 a + be -c 2 π 2 2The symmetrized CGF can therefore be written as sym CGF (t) := ln(M GF f (x) (t)) + ln(M GF f (x) (-t)) = -2 ln(a + be -c 2 π 2 2 ) + t 2 + 2 ln(a + be -c 2 π 2 2 cos(cπt)) Therefore, g(t) = sym CGF (t) -(1 -bc 2 π 2 e -c 2 π 2 2 a+be -c 2 π 2 2)t 2 can be written as :g(t) = -2 ln(a + be -c 2 π 2 2 ) + 2 ln(a + be -c 2 π 2 2 cos(cπt)) + bc 2 π 2 e -c 2 π 2 2 a + be -c 2 π 2 2 t 2Finally, g ′ (t) can be written as :g ′ (t) = d dt 2 ln(a + be -c 2 π 2 2 cos(cπt)) + bc 2 π 2 e -c 2 π 2 2 a + be -c 2 π 2 2 t 2 = - 2bcπe -c 2 π 2 2 sin(cπt) a + be -c 2 π 2 2 cos(cπt) + 2bc 2 π 2 e -c 2 π 2 2 a + be -c 2 π 2 2 tg ′′ (t) can be written as :g ′′ (t) = -2bc 2 π 2 e -c 2 π 2 2 a cos(cπt) + be -c 2 π 2 2 (a + be -c 2 π 2 2 cos(cπt)) 2 + 2bc 2 π 2 e -c 2 π 2 2 a + be -c 2 π 2 2and g ′′′ (t) can be evaluated as:g ′′′ (t) = 2bc 3 π 3 e -c 2 π 2 2 sin(cπt) a 2 -2b 2 e -c 2 π 2 -abe -c 2 π 2 2 cos(cπt) (a + be -c 2 π 2 2 cos(cπt)) 3 Lets set a = 2, b = -1, c = 4 π . Then, we have that the pdf, f (x) = ke -x 2 2 (2 -cos(4x)) = ke -x 2 2 (1 + 2 sin 2 (2x)).The corresponding functions are :EQUATION2 -e -8 cos(4t) Closeness to a Gaussian MGF: It can be seen there that g ′′′ (t) changes sign in the half-line. However, the key is to note that the MGF of this distribution is very close to that of a Gaussian and can be made arbitrarily small by varying the parameters a and b. This renders the CGF-based contrast function ineffectual. This leads us to believe that Assumption 1(d) is not far from being necessary to ensure global optimality of the corresponding objective function. Guidelines:g(t) := sym CGF (t) -var(X)t 2 = -2 ln(2 -e -8 ) + 2 ln(2 -e -8 cos(4t)) - 16e -8 2 -e -8 t 2 g ′ (t) = 8e -8 sin(4t) 2 -e -8 cos(4t) - 32e -8• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper.If B is not invertible, this can be replaced by the Moore-Penrose inverse.Note that yi := xi -x ⌊n/2⌋+i , i = 1 • • • ⌊n/2⌋ has a symmetric distribution, and also remains a noisy version of a linear combination of source signals with the same mixing matrix.MATLAB implementations (under the GNU General Public License) can be found at -FastICA and JADE. The code for PFICA was provided on request by the authors. 


Decision: This paper presents a modification of the noiseless characteristic function-based ICA to the noisy ICA without requiring the knowledge of the Gaussian noise distribution parameters. Two computationally efficient contrast functions are introduced in the paper. Reviewers claimed that the origin of these new contrast functions was not clear, but during the reviewer-author discussion period it has been clarified up to a certain point. One of serious concerns is that the paper seems to a collection of not clearly-developed ideas, rather than exploring coherently and clearly a single idea. I  believe this is partially true. During the AC-reviewer discussion period, there was a mix of support and concern. One reviewer with negative comments agreed that he/she would not mind accepting the paper even though he/she wanted to keep his/her rating as a minority. Consider all the comments raised by reviewers, I think that paper has a contribution in the direction of noisy ICA. I would like to ask the authors to revise the paper carefully, improving the clarity in the final version.

Reviews: review: paper_summary: paper_summary: This work proposes a modification of the original CHFICA characteristic function to the noisy ICA case without requiring knowledge of the noise distribution parameters.
The modified independence score is then used to select the best out of multiple ICA methods based on the score assigned to their solutions, in a per dataset basis. In practice, an average score is evaluated over random directions t, which renders the score computationally inefficient for direct Minimization.
The empirical estimate of the new independence score is shown to converge to the population score if the data covariance is concentrated.
Empirical evidence is shown that the new empirical score correlates well with the well-established Amari error.

This work also introduces two computationally efficient contrast functions to be Maximized (CHF-based and CGF-based) for BSS, both of which do not require higher moments and, thus, work well on data with near-zero kurtosis, like Bernoulli(p) with low p.
Theoretical analyses of local and global convergence are also included for a large class of smooth contrast functions that meet certain desirable properties (Assumption 1). Properties (a) and (c) are critical as they yield contrast functions that are not influenced by additive independent Gaussian noise.
The CHF contrast only requires a finite second moment and is suitable for heavy-tailed sources.
The CGF contrast is not appropriate for heavy-tailed sources.
Global convergence: The Hessian of both contrasts is shown to have the form C = BDB^T. The maxima of the Hessian-adjusted contrasts is shown to occur at B^-1 u = e_i for positive elements of D. Fast sequential power method optimization in a pseudo-Euclidean space defined by C^-1 was shown by [48]. Before optimization, C need only be computed once for some random vector u and just reused during the power iteration updates. A requirement of no sign change of the third derivative of the contrast over the half-lines is required to sufficiently distinguish between Gaussian and near-0 kurtosis source distributions.
Local convergence: Linear convergence is attainable without the third derivative requirement. Geometric convergence is attained at some low epsilon based on how non-gaussian the source is.
Contrary to prior work on cumulant-based methods, convergence of contrasts meeting Assumption 1 was established based on the characterization of the contrasts in the pseudo-Euclidean space and not on the convergence of the power method.

Experimental results support Meta's ability to select the optimal solution over various algorithms, showing variance reduction in the estimates, uniform performance across high and near-0 kurtosis.
Evaluation of the CHF and CGF contrasts show that CHF outperforms other methods at higher noise, but only at sample size n > 20000.
The Meta algorithm picks up the slack of CHF by selecting the competing PFICA solution when n < 20000.
strength_weakness: strength_weakness: Quality:
 - The code is not friendly to readers as there is no 1-to-1 correspondence with the notation in the paper. Suggestion: improve documentation and tidy up the codebase (especially comments and function arguments) for readability.

Clarity:
 - Despite the well-written manuscript and thorough methodology and proofs, the paper contains many typos (both in text and proofs) and inconsistencies in notation (especially between Appendix and main manuscript) that limit the clarity. Besides fixing these typos, explicitly include some of the "obvious" steps currently omitted in the derivation; that would greatly improve readability. Also add some high-level intuitions as noted below.
 - Clarify the origin of new contrast functions and how they came about. It is unclear what motivated or led to these otherwise "semi-arbitrary" functions.
 - Define O_P(.)
 - It appears that the new independence score does not meet some of the properties in Assumption 1. Please, clarify if that is indeed the case and, if so, which properties specifically.
 - Appendix A.3: Clarify that Algo 2 is NOT specific to contrasts based on power iteration only, but any sequential method. It is suggested to make one separate sub-section covering how the projection is done and another for the Algo description itself. Also clarify what the "Total" (non-sequential) Algo would look like, for comparison. 
 - Lines 119-122: t.T should be bold face.
 - Eq. 2: missing last t in the PRODUCT term.
 - Section 3.2: Change "Properties of Contrast Functions 1." to "Assumption 1. (Properties of Contrast Functions)"
 - Line 175: remove "them"
 - Line 182: "note" --> "not"
 - Footnote 2 (page 5): add "and remains a linear combination of z"
 - Theorem 4, line 254: Assumption 1(d) is NOT about the third derivative.
 - Line 271: Define the value of k in this experiment. I see 11 mentioned in passing in line 786...
 - Line 277: Fix contradicting statements: "CHF and CGF are initialized using the B estimated via PFICA" and "CHF and CGF is based on a single random initialization". 
 - Line 283: "also used" --> "also be used"
 - Table 1 caption: Clarify that since median is being reported, Meta is better than the best method “on average”, but on any individual experiment it is identical to the best method.
 - Line 299: refs to 2b and 2c are in flipped order.
 - Line 315: In the noise power experiment, it is stated that "the difference between the two leading algorithms is small." Therefore, the same acknowledgement should be indicated in line 315 after "CHF dominates, and one can see that Meta starts following CHF"
 - Line 456: Assumption 1(d) is NOT about the 3rd derivative constraint...
 - Line 464: Missing ' in g.
 - Line 472: For consistent notation, it should be Delta(t,F | P) = 0, not Delta_F(t) = 0 
 - Lines 474-477: there is an extra transpose symbol at the last t (several times).
 - Line 510: (e) --> (d)
 - Notation and conditions in lines 524-529 is not consistent with Theorem 3 in the main manuscript.
 - Line 777: w_i is not defined, g_i should be ~g_i.
 - Lines 784-789: the r defined here has no relation to the r in Algo 2. Please pick another letter.
 - Algo 2:
	 - Clarify why it is necessary to evaluate Y_j per component and why it would not be equivalent to just compute W = X V^T U^T r
	 - for j in range[1, M] : M is not defined, but I think it would be l + 1? Also, the letter j is used for many different things inside the same loop. Change to " for p in range[1, l + 1] ", and s(j) --> s(p).
	 - " W(:,j) ← Y_j r " --> " W(:,a) ← Y_a r "
	 - I suspect r is equivalent to t in Eq. 2, if so change it to t.
 - Appendix A.4: Verify that all occurrences of Assumption 1 (d) are correct (i.e., referring to symmetry not the third derivative)
questions: questions: Questions were included along with the weaknesses.
limitations: limitations: Yes, the authors addressed the limitations of the work.
paper_summary: paper_summary: The presented paper focuses on the problem of noisy ICA which remains a significant challenge in classical machine learning.
The authors introduces a nonparametric independent score extending the work in [21] to evaluate the estimation of the demixing matrices without requiring any prior knowledge of underlying noise distribution parameters. 
Further. the authors propose some new contrast functions and provides a very interesting discussion on the convergence for the presented contrast functions of noisy ICA which can be also applied for cumulant based contrast functions.
strength_weakness: strength_weakness: Two main weaknesses stood out to me in the paper:
1.  The proposed study heavily rely on Gaussian noise assumptions. The performance of suggested score might not hold in cases with different noise characteristics. Similarly, the proposed contrast functions could have limited scope of applicability in real life scenarios..
2. The meta algorithm's effectiveness is contingent on the pool of candidate algorithms it selects from.
questions: questions: In equation 2 in the exponent of the product term \mathbf{t} seems to be missing?
limitations: limitations: The methods focus primarily on Gaussian noise, which might limit their applicability in scenarios where noise distributions do not conform to this assumption.
The global convergence analysis heavily depends on the assumption that the noise is Gaussian. This dependency raises questions about the generalizability of the convergence results to other types of noise.
paper_summary: paper_summary: This paper proposes a nonparametric score to adaptively pick the best noisy ICA algorithm from a set of candidates. This “independence score” is based on the characteristic function-based objective (CHF) introduced by Eriksson&Koivunen in 2003. 

In practice, this independence score evaluates the inverse mixing matrix obtained by an ICA algorithm without requiring accessing to the true sources. 

In addition, the paper proposes some new contrast functions and algorithms and present simulation results showing the effectiveness of the proposed independence score.
strength_weakness: strength_weakness: •	The proposed approach focused only in the recovery of the proper unmixing matrix (B^{-1}) while in noisy ICA, the final goal is to recover the clean sources. There is no discussion about how to recover the sources after the proper unmixing matrix is obtained.

•	The experimental part focused on synthetically mixed signals only. The paper would be greatly benefited from the inclusion of some real-world source separation problem.

•	Theorem 2 assumes Subgaussian mixtures, which ensures the concentration of the sample covariance matrix. This assumtion is rather strong and it would be good to analyze what happen if it is not met.
questions: questions: Regarding assumption of Theorem 2 (sub-Gaussianity). I understand that assumption is needed to guaranty the concentration of the norm operator, which is a rather strong assumption. Do you think is possible to relax that assumption and keep the method working? Have you make any experiments with super-Gaussian sources, for example?
limitations: limitations: Yes, the experimental results section discuss some limitations of the proposed algorithm in highly noisy scenarios.
paper_summary: paper_summary: The authors consider the noisy ICA problem. They first propose a numerical objective function that can be used as a guide to assess the quality of an existing ICA algorithm. This function is not suitable for optimization and therefore, the authors propose to use it in a meta-algorithm, where the purpose is to select the best ICA solutions out of several, where each solution is possibly produced by a separate algorithm. Then, they propose new contrast functions that are specifically suited for the noisy ICA problem. Finally, they study critical points of contrast functions so as to validate their use in noisy ICA problems. They provide numerical experiments to back up their claims.
strength_weakness: strength_weakness: The paper touches on too many things and is not very readable. For instance, I don't see the connection between the meta algorithm and the new contrast functions. Similarly, I would have welcomed more detail on the development of the meta algorithm, and the choice of randomization in its computation. 
Occasionally, something is introduced out of nowhere and it's hard to understand the motivation. For instance, I don't see the development for the new contrast functions. 
Overall, I'd have preferred a paper that coherently and clearly explores a single idea. Currently, it reads like a collection of not clearly-developed ideas.
questions: questions: - line 123, "We propose to adapt the CHFICA objective using estimable parameters to the noisy ICA setting" : I don't understand this sentence. What are the "estimable parameters" here? Is it $S$?
- eqn 2 : What's the motivation for the additional terms (the second factors in both "JOINT" and "PRODUCT")?
- eqn 2 : In the last term, is a $t$ missing?
- Thm 2 : How does knowing this bound help? Please discuss how this theorem is useful in the text.
- line 175, "...in this section them, ..." : typo?
- line 214 : What is a pseudo-Euclidean space?
limitations: limitations: N/A

Comments: title: Clarification on how the different ideas are tied together
Post author rebuttal comments
Thanks
comment: Thank you for your questions and for giving us an opportunity to further clarify. We address them below:

**Regarding $O_P$ and $O$**: The $X_n=O_P(a_n)$ implies that $X_n/a_n$ is stochastically bounded (or more technically uniformly tight), i.e. $\forall \epsilon>0$, there exists a finite $M$ such that $\sup_n P(|X_n/a_n|>M)\leq \epsilon$ (See [A]). We will clarify this.

**Regarding Independence Score and Property (a)**: We will do so.

**Regarding dimension of dataset**: It should be the 9-dimensional dataset. We will fix this typographical error.

**Regarding random initialization of CGF and CHF**: Both are correct. Recall that, for the power iteration-based algorithm of Voss et. al (2015) we need both a quasi-orthogonalization matrix which is an estimator of a matrix of the form $BDB^T$, and a unit vector $u$. The quasi-orthogonalization matrices for CHF and CGF use the $B$ matrix estimated via PFICA. But the initial vector $u$ for the power iteration is a random unit vector. We will clarify this further in the revised manuscript.

**Regarding the origin of new contrast functions**: We started from the fact that, like the cumulants, the cumulant generating function (and similarly the CHF-based counterpart) satisfies Assumption 1a). In order to satisfy 1c), we wish to subtract out the part resulting from the Gaussian noise, which led to the additional terms involving $u^{\top}Su$. We will add this explanation to better address your remark.

References :
[A] Vaart, A. W. van der. 1998. Asymptotic Statistics. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press.
Thank you for your response.
I have some follow up questions:

1. How is O_P different from just O?

2. Include a note that "The new independence score does not satisfy property (a) of Assumption 1"

3. Clarify which one is the 11-dimensional dataset mentioned in line 786... (line 297 says k=9).

4. Line 277: So which is the correct statement: "CHF and CGF are initialized using the B estimated via PFICA" or "CHF and CGF is based on a single random initialization"?

5. Origin of new contrast functions: How do you arrive at Eq. (4) from the CHF? Likewise, how do you arrive at Eq. (4) from the CGF? Although the conditions in Assumption 1 are satisfied by the new contrast functions, the conditions do not specify what f() should be. Currently the train of thought/rationale/derivation from CHF -> f() is missing (likewise from CGF -> f() ). Assumption 1 is insufficient to clarify this point.
Dear Reviewer uqTT,

We are grateful for your response and for giving us another opportunity to explain why we believe that the different ideas or components in this paper really belong together. To show the efficacy of our meta-algorithm, we need to have a candidate pool of contrast functions whose weaknesses complement each other. For example, kurtosis is a widely used cumulant-based method for noisy ICA that has been previously shown to suffer in the presence of heavy-tailed distributions or distributions with small kurtosis. So, we also need to have contrast functions that are not cumulant-based and are less sensitive to heavy-tailed source signals. While there are many such contrast functions for noiseless ICA, there are not many with provable guarantees for noisy ICA. This is why we designed our CHF-based contrast function and developed a theoretical framework for analysis under the noisy ICA model. We explained this between lines 87-110. We will be happy to add a longer discussion about this at the beginning of the paper to better address your remarks.
I thank the authors for their responses. They do help clarify the specific points I raised. I'd be willing to slightly increase my score. However, my main criticism, that the paper contains a number of disparate ideas which may as well be split, still stands -- to be fair, it may not be easy to address that with a minor revision.
Dear reviewer c27m, 

Thank you - we are grateful for your review and response.

We just wanted to clarify that the first weakness (re: recovery of source signals) is not exactly a weakness because it is a property of the noisy ICA model. We addressed this in lines 60-65 of our paper, and we are happy to add a more detailed discussion. Similarly, the third point (re: sub-gaussianity assumption in theorem 2) can be replaced with a much weaker condition on finite moments (as pointed out in our rebuttal).
I appreciate the authors' responses to my comments. However, despite the weaknesses not being serious, they still persist, so I will maintain my original score.