A Pseudo-Euclidean Iteration for Optimal Recovery in Noisy ICA

A Pseudo-Euclidean Iteration for Optimal Recovery in Noisy ICA

James Voss111Department of Computer Science and Engineering, the Ohio State University
vossj@cse.ohio-state.edu
   Mikhail Belkinfootnotemark:
mbelkin@cse.ohio-state.edu
   Luis Rademacherfootnotemark:
lrademac@cse.ohio-state.edu
Abstract

Independent Component Analysis (ICA) is a popular model for blind signal separation. The ICA model assumes that a number of independent source signals are linearly mixed to form the observed signals. We propose a new algorithm, PEGI (for pseudo-Euclidean Gradient Iteration), for provable model recovery for ICA with Gaussian noise. The main technical innovation of the algorithm is to use a fixed point iteration in a pseudo-Euclidean (indefinite “inner product”) space. The use of this indefinite “inner product” resolves technical issues common to several existing algorithms for noisy ICA. This leads to an algorithm which is conceptually simple, efficient and accurate in testing.

Our second contribution is combining PEGI with the analysis of objectives for optimal recovery in the noisy ICA model. It has been observed that the direct approach of demixing with the inverse of the mixing matrix is suboptimal for signal recovery in terms of the natural Signal to Interference plus Noise Ratio (SINR) criterion. There have been several partial solutions proposed in the ICA literature. It turns out that any solution to the mixing matrix reconstruction problem can be used to construct an SINR-optimal ICA demixing, despite the fact that SINR itself cannot be computed from data. That allows us to obtain a practical and provably SINR-optimal recovery method for ICA with arbitrary Gaussian noise.

1 Introduction

Independent Component Analysis refers to a class of methods aiming at recovering statistically independent signals by observing their unknown linear combination. There is an extensive literature on this and a number of related problems [Comon2010].

Specifically, in the ICA model, we observe -dimensional realizations of a latent variable model where denotes the th column of the mixing matrix and is the unseen latent random vector of “signals”. It is assumed that are independent and non-Gaussian. The source signals and entries of may be either real- or complex-valued. For simplicity, we will assume throughout that has zero mean, as this may be achieved in practice by centering the observed data.

Many ICA algorithms use the preprocessing “whitening” step whose goal is to orthogonalize the independent components. In the noiseless, case this is commonly done by computing the square root of the covariance matrix of . Consider now the noisy ICA model with additive -mean noise independent of . It turns out that the introduction of noise makes accurate recovery of the signals significantly more involved. Specifically, whitening using the covariance matrix does not work in the noisy ICA model as the covariance matrix combines both signal and noise. For the case when the noise is Gaussian, matrices constructed from higher order statistics (specifically, cumulants) can be used instead of the covariance matrix. However, these matrices are not in general positive definite and thus the square root cannot always be extracted. This limits the applicability of several previous methods, such as [AroraGMS12, 4203062, albera2004blind]. The GI-ICA algorithm proposed in [NIPS2013_5134] addresses this issue by using a complicated quasi-orthogonalization step followed by an iterative method.

In this paper (section 2), we develop a simple and practical one-step algorithm, PEGI (for pseudo-Euclidean Gradient Iteration) for provably recovering (up to the unavoidable ambiguities of the model) in the case when the noise is Gaussian (with an arbitrary, unknown covariance matrix). The main technical innovation of our approach is to formulate the recovery problem as a fixed point method in an indefinite “inner product” (pseudo-Euclidean) space.

The second contribution of the paper is combining PEGI with the analysis of objectives for optimal recovery in the noisy ICA model. In most applications of ICA (e.g., speech separation [makino2007blind], MEG/EEG artifact removal [vigario2000independent] and others) one cares about recovering the signals . This is known as the source recovery problem. This is typically done by first recovering the matrix (up to an appropriate scaling of the column directions).

At first, source recovery and recovering the mixing matrix appear to be essentially equivalent. In the noiseless ICA model, if in invertible222 can be replaced with in the discussion below for over-determined ICA. then recovers the sources. On the other hand, in the noisy model, the exact recovery of the latent sources becomes impossible even if is known exactly. Part of the “noise” can be incorporated into the “signal” preserving the form of the model. Even worse, neither nor are defined uniquely as there is an inherent ambiguity in the setting. There could be many equivalent decompositions of the observed signal as (see the discussion in section 3).

We consider recovered signals of the form for a choice of demixing matrix . Signal recovery is considered optimal if the coordinates of maximize Signal to Interference plus Noise Ratio (SINR) within any fixed model . Note that the value of SINR depends on the decomposition of the observed data into “noise” and “signal”: .

Surprisingly, the SINR optimal demixing matrix does not depend on the decomposition of data into signal plus noise. As such, SINR optimal ICA recovery is well defined given access to data despite the inherent ambiguity in the model. Further, it will be seen that the SINR optimal demixing can be constructed from and the directions of the columns of (which are also well-defined across signal/noise decompositions).

Our SINR-optimal demixing approach combined with the PEGI algorithm provides a complete SINR-optimal recovery algorithm in the ICA model with arbitrary Gaussian noise. We note that the ICA papers of which we are aware that discuss optimal demixing do not observe that SINR optimal demixing is invariant to the choice of signal/noise decomposition. Instead, they propose more limited strategies for improving the demixing quality within a fixed ICA model. For instance, joho2000overdetermined show how SINR-optimal demixing can be approximated with extra sensors when assuming a white additive noise, and koldovsky2007asymptotic discuss how to achieve asymptotically low bias ICA demixing assuming white noise within a fixed ICA model. However, the invariance of the SINR-optimal demixing matrix appears in the array sensor systems literature [Chevalier199927].

Finally, in section 4, we demonstrate experimentally that our proposed algorithm for ICA outperforms existing practical algorithms at the task of noisy signal recovery, including those specifically designed for beamforming, when given sufficiently many samples. Moreover, most existing practical algorithms for noisy source recovery have a bias and cannot recover the optimal demixing matrix even with infinite samples. We also show that PEGI requires significantly fewer samples than GI-ICA [NIPS2013_5134] to perform ICA accurately.

1.1 The Indeterminacies of ICA

Notation: We use to denote the entry-wise complex conjugate of a matrix , to denote its transpose, and to denote its conjugate transpose.

Before proceeding with our results, we discuss the somewhat subtle issue of indeterminacies in ICA. These ambiguities arise from the fact that the observed may have multiple decompositions into ICA models and .

Noise-free ICA has two natural indeterminacies. For any nonzero constant , the contribution of the th component to the model can equivalently be obtained by replacing with and with the rescaled signal . To lessen this scaling indeterminacy, we use the convention333Alternatively, one may place the scaling information in the signals by setting for each . that throughout this paper. As such, each source (or equivalently each ) is defined up to a choice of sign (a unit modulus factor in the complex case). In addition, there is an ambiguity in the order of the latent signals. For any permutation of (where ), the ICA models and are indistinguishable. In the noise free setting, is said to be recovered if we recover each column of up to a choice of sign (or up to a unit modulus factor in the complex case) and an unknown permutation. As the sources are only defined up to the same indeterminacies, inverting the recovered matrix to obtain a demixing matrix works for signal recovery.

In the noisy ICA setting, there is an additional indeterminacy in the definition of the sources. Consider to be -mean axis-aligned Gaussian random vector. Then, the noisy ICA model in which is considered part of the latent source signal , and the model in which is part of the noise are indistinguishable. In particular, the latent source and its covariance are ill-defined. Due to this extra indeterminacy, the lengths of the columns of no longer have a fully defined meaning even when we assume . In the noisy setting, is said to be recovered if we obtain the columns of up to non-zero scalar multiplicative factors and an arbitrary permutation.

The last indeterminacy is the most troubling as it suggests that the power of each source signal is itself ill-defined in the noisy setting. Despite this indeterminacy, it is possible to perform an SINR-optimal demixing without additional assumptions about what portion of the signal is source and what portion is noise. In section 3, we will see that SINR-optimal source recovery takes on a simple form: Given any solution which recovers up to the inherent ambiguities of noisy ICA, then is an SINR-optimal demixing matrix.

1.2 Related Work and Contributions

Independent Component Analysis is probably the most used model for Blind Signal Separation. It has seen numerous applications and has generated a vast literature, including in the noisy and underdetermined settings. We refer the reader to the books [Hyvarinen2001, Comon2010] for a broad overview of the subject.

It was observed early on by cardoso1991super that ICA algorithms based soley on higher order cumulant statistics are invariant to additive Gaussian noise. This observation has allowed the creation of many algorithms for recovering the ICA mixing matrix in the noisy and often underdetermined settings. Despite the significant work on noisy ICA algorithms, they remain less efficient, more specialized, or less practical than the most popular noise free ICA algorithms.

Research on cumulant-based noisy ICA can largely be split into several lines of work which we only highlight here. Some algorithms such as FOOBI [cardoso1991super] and BIOME [albera2004blind] directly use the tensor structure of higher order cumulants. In another line of work, de1996independent and yeredor2002non have suggested algorithms which jointly diagonalize cumulant matrices in a manner reminiscent of the noise-free JADE algorithm [Cardoso1993]. In addition, Yeredor2000 and DBLP:journals/corr/GoyalVX13 have proposed ICA algorithms based on random directional derivatives of the second characteristic function.

Each line of work has its advantages and disadvantages. The joint diagonalization algorithms and the tensor based algorithms tend to be practical in the sense that they use redundant cumulant information in order to achieve more accurate results. However, they have a higher memory complexity than popular noise free ICA algorithms such as FastICA [Hyvarinen2000]. Moreover, the tensor methods (FOOBI and BIOME) require the latent source signals to have positive order (, a predetermined fixed integer) cumulants as they rely on taking a matrix square root. Finally, the methods based on random directional derivatives of the second characteristic function rely heavily upon randomness in a manner not required by the most popular noise free ICA algorithms.

We continue a line of research started by AroraGMS12 and NIPS2013_5134 on fully determined noisy ICA which addresses some of these practical issues by using a deflationary approach reminiscent of FastICA. Their algorithms thus have lower memory complexity and are more scalable to high dimensional data than the joint diagonalization and tensor methods. However, both works require a preprocessing step (quasi-orthogonalization) to orthogonalize the latent signals which is based on taking a matrix square root. AroraGMS12 require each latent signal to have positive fourth cumulant in order to carry out this preprocessing step. In contrast, NIPS2013_5134 are able to perform quasi-orthogonalization with source signals of mixed sign fourth cumulants; but their quase-orthogonalization step is more complicated and can run into numerical issues under sampling error. We demonstrate that quasi-orthogonalization is unnecessary. We introduce the PEGI algorithm to work within a (not necessarily positive definite) inner product space instead. Experimentally, this leads to improved demixing performance. In addition, we handle the case of complex signals.

Finally, another line of work attempts to perform SINR-optimal source recovery in the noisy ICA setting. It was noted by koldovsky2006methods that for noisy ICA, traditional ICA algorithms such as FastICA and JADE actually outperform algorithms which first recover in the noisy setting and then use the resulting approximation of to perform demixing. It was further observed that is not the optimal demixing matrix for source recovery. Later, koldovsky2007blind proposed an algorithm based on FastICA which performs a low SINR-bias beamforming.

2 Pseudo-Euclidean Gradient Iteration ICA

In this section, we introduce the PEGI algorithm for recovering in the “fully determined” noisy ICA setting where . PEGI relies on the idea of Gradient Iteration introduced NIPS2013_5134. However, unlike GI-ICA NIPS2013_5134, PEGI does not require the source signals to be orthogonalized. As such, PEGI does not require the complicated quasi-orthogonalization preprocessing step of GI-ICA which can be inaccurate to compute in practice. We sketch the Gradient Iteration algorithm in Section 2.1, and then introduce PEGI in Section 2.2. For simplicity, we limit this discussion to the case of real-valued signals. We show how to construct PEGI for complex-valued signals in Appendix A.

In this section we assume a noisy ICA model such that is arbitrary Gaussian and independent of . We also assume that , that is known, and that the columns of are linearly independent.

2.1 Gradient Iteration with Orthogonality

The gradient iteration relies on the properties of cumulants. We will focus on the fourth cumulant, though similar constructions may be given using other even order cumulants of higher order. For a zero-mean random variable , the fourth order cumulant may be defined as  [see Comon2010, Chapter 5, Section 1.2]. Higher order cumulants have nice algebraic properties which make them useful for ICA. In particular, has the following properties: (1) (Independence) If and are independent, then . (2) (Homogeneity) If is a scalar, then . (3) (Vanishing Gaussians) If is normally distributed then .

We consider the following function defined on the unit sphere: . Expanding using the above properties we obtain:

Taking derivatives we obtain:

(1)
(2)

where is a diagonal matrix with entries .

NIPS2013_5134 introduced GI-ICA as a fixed point algorithm under the assumption that the columns of are orthogonal but not necessarily unit vectors. The main idea is that the update is a form of a generalized power iteration. From equation (1), each may be considered as a direction in a hidden orthogonal basis of the space. During each iteration, the coordinate of is raised to the rd power and multiplied by a constant. Treating this iteration as a fixed point update, it was shown that given a random starting point, this iterative procedure converges rapidly to one of the columns of (up to a choice of sign). The rate of convergence is cubic.

However, the GI-ICA algorithm requires a somewhat complicated preprocessing step called quasi-orthogonalization to linearly transform the data to make columns of orthogonal. Quasi-orthogonalization makes use of evaluations of Hessians of the fourth cumulant function to construct a matrix of the form where has all positive diagonal entries—a task which is complicated by the possibility that the latent signals may have fourth order cumulants of differing signs—and requires taking the matrix square root of a positive definite matrix of this form. However, the algorithm used for constructing under sampling error is not always positive definite in practice, which can make the preprocessing step fail. We will show how our PEGI algorithm makes quasi-orthogonalization unnecessary, in particular, resolving this issue.

2.2 Gradient Iteration in a Pseudo-Euclidean Space

We now show that the gradient iteration can be performed using in a pseudo-Euclidean space in which the columns of are orthogonal. The natural candidate for the “inner product space” would be to use defined as . Clearly, gives the desired orthogonality property. However, there are two issues with this “inner product space”: First, it is only an inner product space when is invertible. This turns out not to be a major issue, and we move forward largely ignoring this point. The second issue is more fundamental: We only have access to in the noise free setting where . In the noisy setting, we have access to matrices of the form from equation (2) instead.

Inputs: Unit vector , ,
repeat
     
     
until Convergence (up to sign)
return
Algorithm 1 Recovers a column of up to a scaling factor if is generically chosen.

consider a pseudo-Euclidean inner product defined as follows: Let where is a diagonal matrix with non-zero diagonal entries, and define by . When contains negative entries, this is not a proper inner product since is not positive definite. In particular, may be negative. Nevertheless, when , gives that the columns of are orthogonal in this space.

We define functions by such that for any , then is the expansion of in its basis. Continuing from equation (1), for any we see is the gradient iteration recast in the space. Expanding in its basis, we obtain

(3)

which is a power iteration in the unseen coordinate system. As no assumptions are made upon the values, the scalings which were not present in eq. (1) cause no issues. Using this update, we obtain Alg. 1, a fixed point method for recovering a single column of up to an unknown scaling.

Before proceeding, we should clarify the notion of fixed point convergence in Algorithm 1. We say that the sequence converges to up to sign if there exists a sequence such that each and as . We have the following convergence guarantee.

Theorem 1.

If is chosen uniformly at random from , then with probability 1, there exists such that the sequence defined as in Algorithm 1 converges to up to sign. Further, the rate of convergence is cubic.

Due to limited space, we omit the proof of Theorem 1. It is similar to the proof of [NIPS2013_5134, Theorem 4].

In practice, we test near convergence by checking if we are still making significant progress. In particular, for some predefined , if there exists a sign value such that , then we declare convergence achieved and return the result. As there are only two choices for , this is easily checked, and we exit the loop if this condition is met.

Full ICA Recovery Via the Pseudo-Euclidean GI-Update. We are able to recover a single column of up to its unknown scale. However, for full recovery of , we would like (given recovered columns ) to be able to recover a column such that on demand.

The idea behind the simultaneous recovery of all columns of is two-fold. First, instead of just finding columns of using Algorithm 1, we simultaneously find rows of . Then, using the recovered columns of and rows of , we project onto the orthogonal complement of the recovered columns of within the pseudo-inner product space.

Recovering rows of .

Suppose we have access to a column (which may be achieved using Algorithm 1). Let denote the th row of . Then, we note that recovers up to an arbitrary, unknown constant . However, the constant may be recovered by noting that . As such, we may estimate as .

1:Inputs: ,
2:,
3:for  to  do
4:     Draw uniformly at random from .
5:     repeat
6:         
7:         .
8:     until Convergence (up to sign)
9:     
10:     
11:end for
12:return ,
Algorithm 2 Full ICA matrix recovery algorithm. Returns two matrices: (1) is the recovered mixing matrix for the noisy ICA model , and (2) is a running estimate of .

Enforcing Orthogonality During the GI Update. Given access to a vector (where is the projection onto the orthogonal complements of the range of ), some recovered columns , and corresponding rows of , we may zero out the components of corresponding to the recovered columns of . Letting , then . In particular, is orthogonal (in the space) to the previously recovered columns of . This allows the non-orthogonal gradient iteration algorithm to recover a new column of .

Using these ideas, we obtain Algorithm 2, which is the PEGI algorithm for recovery of the mixing matrix in noisy ICA up to the inherent ambiguities of the problem. Within this Algorithm, step 6 enforces orthogonality with previously found columns of , guaranteeing that convergence to a new column of .
Practical Construction of . In our implementation, we set , as it can be shown from equation (2) that with . This deterministically guarantees that each latent signal has a significant contribution to .

3 SINR Optimal Recovery in Noisy ICA

In this section, we demonstrate how to perform SINR optimal ICA within the noisy ICA framework given access to an algorithm (such as PEGI) to recover the directions of the columns of . To this end, we first discuss the SINR optimal demixing solution within any decomposition of the ICA model into signal and noise as . We then demonstrate that the SINR optimal demixing matrix is actually the same across all possible model decompositions, and that it can be recovered. The results in this section hold in greater generality than in section 2. They hold even if (the underdetermined setting) and even if the additive noise is non-Gaussian.

Consider an demixing matrix, and define the resulting approximation to . It will also be convenient to estimate the source signal one coordinate at a time: Given a row vector , we define . If (the th row of ), then is our estimate to the th latent signal . Within a specific ICA model , signal to intereference-plus-noise ratio (SINR) is defined by the following equation:

(4)

is the variance of the contribution of th source divided by the variance of the noise and interference contributions within the signal.

Given access to the mixing matrix , we define . Since , this may be rewritten as . Here, may be estimated from data, but due to the ambiguities of the noisy ICA model, (and specifically its column norms) cannot be estimated from data.

koldovsky2006methods observed that when is a white Gaussian noise, jointly maximizes for each , i.e., takes on its maximal value at . Below in Proposition 2, we generalize this result to include arbitrary non-spherical, potentially non-Gaussian noise.

It is interesting to note that even after the data is whitened, i.e. , the optimal SINR solution is different from the optimal solution in the noiseless case unless is an orthogonal matrix, i.e. . This is generally not the case, even if is white Gaussian noise.

Proposition 2.

For each , is a maximizer of .

The proof of Proposition 2 is deferred to appendix B.

Since is scale invariant, Proposition 2 implies that any matrix of the form where is a diagonal scaling matrix (with non-zero diagonal entries) is an -optimal demixing matrix. More formally, we have the following result.

Theorem 3.

Let be an matrix containing the columns of up to scale and an arbitrary permutation. That is, there exists a permutation of and non-zero constants such that for each . Then, is a maximizer of .

By Theorem 3, given access to a matrix which recovers the directions of the columns of , then is the SINR-optimal demixing matrix. For ICA in the presence of Gaussian noise, the directions of the columns of are well defined simply from , that is, the directions of the columns of do not depend on the decomposition of into signal and noise (see the discussion in section 1.1 on ICA indeterminacies). The problem of SINR optimal demixing is thus well defined for ICA in the presence of Gaussian noise, and the SINR optimal demixing matrix can be estimated from data without any additional assumptions on the magnitude of the noise in the data.

Finally, we note that in the noise-free case, the SINR-optimal source recovery simplifies to be .

Corollary 4.

Suppose that is a noise free (possibly underdetermined) ICA model. Suppose that contains the columns of up to scale and permutation, i.e., there exists diagonal matrix with non-zero entries and a permutation matrix such that . Then is an SINR-optimal demixing matrix.

Proof.

By Theorem 3, is an SINR-optimal demixing matrix. Expanding, we obtain: . ∎

Corollary 4 is consistent with known beamforming results. In particular, it is known that is an optimal (in terms of minimum mean squared error) beamforming matrix for underdetermined ICA [van1988beamforming, section 3B].

4 Experimental Results

(a) Accuracy under additive Gaussian noise.
(b) Bias under additive Gaussian noise.
Figure 1: SINR performance comparison of ICA algorithms.

We now compare the proposed PEGI algorithm with several existing ICA algorithms. In addition to qorth+GI-ICA (that is, GI-ICA with the quasi-orthogonalization preprocessing step), we use the following algorithmic baselines:
JADE [Cardoso1993] is a popular fourth cumulant based ICA algorithm designed for the noise free setting. We use the implementation of JADE.
FastICA [Hyvarinen2000] is a popular ICA algorithm designed for the noise free setting based on a deflationary approach of recovering one component at a time. We use the implementation of FastICA.
1FICA [koldovsky2007asymptotic, koldovsky2007blind] is a variation of FastICA with the tanh contrast function designed to have low bias for performing SINR-optimal beamforming in the presence of Gaussian noise.
Ainv is the oracle demixing algorithm which uses as the demixing matrix.
SINR-opt is the oracle demixing algorithm which demixes using to achieve an SINR-optimal demixing.

We compare these algorithms on simulated data with . We constructed mixing matrices with condition number 3 via a reverse singular value decomposition (). The matrices and were random orthogonal matrices, and was chosen to have 1 as its minimum and 3 as its maximum singular values, with the intermediate singular values chosen uniformly at random. We drew data from a noisy ICA model where was chosen to be malaligned with . We set where is a constant defining the noise power. It can be shown that is the ratio of the maximum directional noise variance to the maximum directional signal variance. We generated 100 matrices for our experiments with 100 corresponding ICA data sets for each sample size and noise power. When reporting results, we apply each algorithm to each of the 100 data sets for the corresponding sample size and noise power and we report the mean performance. The source distributions used in our ICA experiments were the Laplace and Bernoulli distribution with parameters 0.05 and 0.5 respectively, the -distribution with 3 and 5 degrees of freedom respectively, the exponential distribution, and the uniform distribution. Each distribution was normalized to have unit variance, and the distributions were each used twice to create 14-dimensional data. We compare the algorithms using either SINR or the SINR loss from the optimal demixing matrix (defined by SINR Loss = [Optimal SINR Achieved SINR]).

In Figure 1, we compare our proprosed ICA algorithm with various ICA algorithms for signal recovery. In the PEGI-+SINR algorithm, we use PEGI- to estimate , and then perform demixing using the resulting estimate of , the formula for SINR-optimal demixing. It is apparent that when given sufficient samples, PEGI-+SINR provides the best SINR demixing. JADE, FastICA-tanh, and 1FICA each have a bias in the presence of additive Gaussian noise which keeps them from being SINR-optimal even when given many samples.

Figure 2: Accuracy comparison of PEGI using pseudo-inner product spaces and GI-ICA using quasi-orthogonalization.

n Figure 1, we compare algorithms at various sample sizes. The PEGI-+SINR algorithm relies more heavily on accurate estimates of fourth order statistics than JADE, and the FastICA-tanh and 1FICA algorithms do not require the estimation of fourth order statistics. For this reason, PEGI-+SINR requires more samples than the other algorithms in order to be run accurately. However, once sufficient samples are taken, PEGI-+SINR outperforms the other algorithms including 1FICA which is designed to have low SINR bias.

In order to avoid clutter, we did not include qorth+GI-ICA-+SINR (the SINR optimal demixing estimate constructed using qorth+GI-ICA- to estimate ) in the figures 1 and 1. It is also assymptotically unbiased in estimating the directions of the columns of , and similar conclusions could be drawn using qorth+GI-ICA- in place of PEGI-. However, in Figure 2, we see that PEGI-+SINR requires fewer samples than qorth+GI-ICA-+SINR to achieve good performance. This is particularly highlighted in the medium sample regime.

On the Performance of Traditional ICA Algorithms for Noisy ICA.

An interesting observation [first made in koldovsky2006methods] is that the popular noise free ICA algorithms JADE and FastICA perform reasonably well in the noisy setting. In Figures 1 and 1, they significantly outperform demixing using for source recovery. It turns out that this may be explained by a shared preprocessing step. Both JADE and FastICA rely on a whitening preprocessing step in which the data are linearly transformed to have identity covariance. It can be shown in the noise free setting that after whitening, the mixing matrix is a rotation matrix. These algorithms proceed by recovering an orthogonal matrix to approximate the true mixing matrix . Demixing is performed using . Since the data is white (has identity covariance), then the demixing matrix is an estimate of the SINR-optimal demixing matrix. Nevertheless, the traditional ICA algorithms give a biased estimate of under additive Gaussian noise.

References

Appendix A PEGI for Complex Signals

In Section 2, we showed how to perform gradient iteration ICA within a pseudo-Euclidean inner product space. In this appendix, we show how this PEGI algorithm can be extended to include complex valued signals. For clarity, we repeat the entire PEGI algorithmic construction from Section 2 with the necessary modifications to handle the complex setting.

Throughout this appendix, we assume a noisy ICA model where is an arbitrary Gaussian noise independent of . We also assume that , that is known, and that the columns of are linearly dependent.

a.1 Fourth Cumulants of Complex Variables

The gradient iteration relies on the properties of cumulants. We will focus on the fourth cumulant, though similar constructions may be given using other even order cumulants of higher order. We will use two versions of the fourth cumulant which capture slightly different fourth order information. For a zero-mean random variable , they may be defined as and . For real random variables, these two definitions are equivalent, and they come from two different conjugation schemes when constructing the fourth order cumulant [see Comon2010, Chapter 5, Section 1.2]. However, in general, only is guaranteed to be real valued. The higher order cumulants have nice algebraic properties which make them useful for ICA:

  1. (Independence) If and are independent random variables, then and .

  2. (Homogeneity) If is a scalar, then and .

  3. (Vanishing Gaussians) If is normally distributed then and .

In this appendix, we consider a noisy ICA model where is a -mean (possibly complex) Gaussian and independent of . We consider the following functions defined on the unit sphere: and . Then, expanding using the above properties we obtain:

Using similar reasoning, it can be seen that .

It turns out that some slightly non-standard notions of derivatives are most useful in constructing the gradient iteration in the complex setting. We use real derivatives for the gradient and we use the complex Hessian. In particular, expanding , we use the gradient operator . We make use of the operators and to define . Applying this version of the Hessian is different than using real derivatives as in the gradient operation.

Taking derivatives, we obtain:

(5)
(6)

where is a diagonal matrix with entries .

a.2 Gradient Iteration in a Pseudo-Euclidean Space

We now demonstrate that the gradient iteration can be performed using a generalized notion of an inner product space in which the columns of are orthogonal. The natural candidate for the “inner product space” would be to use defined as . Clearly, gives the desired orthogonality property. However, there are two issues with this “inner product space”: First, it is only an inner product space when is non-singular (invertible). This turns out not to be a major issue, and we will move forward largely ignoring this point. The second issue is more fundamental: We only have access to the matrix in the noise free setting where . In the noisy setting, we have access to matrices of the form from equation (A.1) instead.

We consider a pseudo-Euclidean inner product defined as follows: Let where is a diagonal matrix with non-zero diagonal entries, and define by . When contains negative entries, this is not a proper inner product since is not positive definite. In particular, may be negative. Nevertheless, when , gives that the columns of are orthogonal in this space.

We define functions by such that for any , then is the expansion of in its basis. Continuing from equation (5), for any we see

is the gradient iteration recast in the space. Expanding in its basis, we obtain

(7)

which is a power iteration in the unseen coordinate system. As no assumptions are made upon the values, the scalings which were not present in equation (5) cause no issues. Using this update, we obtain Algorithm 3, a fixed point method for recovering a single column of up to an unknown scaling.

Inputs: (A unit vector), ,
repeat
     
     
until Convergence (up to a unit modulus factor)
return
Algorithm 3 Recovers a column of up to an unknown scaling factor when is generically chosen.

Before proceeding, we should clarify the notion of fixed point convergence in Algorithm 3. We say that the sequence converges to up to a unit modulus factor if there exists a sequence of constants such that each and as . We have the following convergence guarantee.

Theorem 5.

If is chosen uniformly at random from . Then with probability 1, there exists such that the sequence defined as in Algorithm 3 converges to a up to a unit modulus factor. Further, the rate of convergence is cubic.

Due to space limitations, we omit the proof of Theorem 5. However, its proof is very similar that of an analogous result for the GI-ICA algorithm [NIPS2013_5134, Theorem 4].

In practice, we test near convergence by testing if we are still making significant progress. In particular, for some predefined , if there exists a unit modulus constant such that , then we declare convergence achieved and return the result. We may determine using the following fact.

Fact 6.

Suppose that and are non-orthogonal unit modulus vectors. The expression is minimized by the choice of .

Letting , we exit the loop if .

a.3 Full ICA Recovery Via the Pseudo-Euclidean GI-Update

We are able to recover a single column of in noisy ICA. However, for full matrix recovery, we would like (given recovered columns ) to be able to recover a column such that on demand.

The main idea behind the simultaneous recovery of all columns of is two-fold. First, instead of just finding columns of using Algorithm 3, we simultaneously find rows of . Then, using the recovered columns of and rows of , we may project onto the orthogonal complement of the recovered columns of within the pseudo-Euclidean inner product space.

Recovering rows of .

Suppose we have access to a column (which may be achieved using Algorithm 3). Let denote the th row of . Then, we note that recovers up to an arbitrary, unknown constant . However, the constant may be recovered by noting that . As such, we may estimate as .

1:Inputs: ,
2:,
3:for  to  do
4:     Draw uniformly at random from .
5:     repeat
6:         
7:         .
8:     until Convergence (up to a unit modulus factor)
9:     
10:     
11:end for
12:return ,
Algorithm 4 Full ICA matrix recovery algorithm. Estimates and returns two matrices: (1) is the recovered mixing matrix for the noisy ICA model , and (2) is a running estimate of .

Enforcing Orthogonality During the GI Update.

Given access to , some recovered columns , and corresponding rows of , we may zero out the components of corresponding to the recovered columns of . Letting , then . In particular, is orthogonal (in the space) to the previously recovered columns of . This allows us to modify the non-orthogonal gradient iteration algorithm to recover a new column of .

Using these ideas, we obtain the Algorithm 4 for recovery of the ICA mixing matrix. Within this Algorithm, step 6 enforces orthogonality with previously found columns of , guaranteeing that convergence is to a new column of .

Practical Construction of

We suggest the choice of , as it can be shown from equation (A.1) that with . This deterministically guarantees that each latent signal has a significant contribution to .

Appendix B Proof of Proposition 2

Proof.

This proof is based on the connection between two notions of optimality, minimum mean squared error and SINR. The mean squared error of the recovered signal from th latent signal is defined as . It has been shown [joho2000overdetermined, equation 39] that jointly minimizes the mean squared errors of the recovered signals. In particular, if , then is a minimizer of .

We will first show that finding a matrix which minimizes the mean squared error has the side effect of maximizing the magnitude of the Pearson correlations for each , where . We will then demonstrate that if is a maximizer of , then is a maximizer of . These two facts imply the desired result. We will use the convention that is 0 if .

We fix a . We have:

Letting , we obtain

(8)

Further, with equality if and only if is real and non-negative. As such, all global minima of are contained in the set , and we may restrict our investigation to this set.

We define a function such that under the change of variable and