Geometric Methods for Robust Data Analysis in High Dimension
Joseph Anderson
TL;DR
This paper introduces the first provably efficient algorithms for Independent Component Analysis (ICA) with heavy-tailed data, utilizing the centroid body and convex geometry tools to improve robustness and applicability in high-dimensional data analysis.
Contribution
It presents novel algorithms that leverage the centroid body and convex geometry for robust ICA and Gaussian mixture model learning in high-dimensional, heavy-tailed data settings.
Findings
First efficient ICA algorithm for heavy-tailed data.
Effective transformation of samples from complex distributions.
Successful recovery of mixture parameters in high dimensions.
Abstract
Machine learning and data analysis now finds both scientific and industrial application in biology, chemistry, geology, medicine, and physics. These applications rely on large quantities of data gathered from automated sensors and user input. Furthermore, the dimensionality of many datasets is extreme: more details are being gathered about single user interactions or sensor readings. All of these applications encounter problems with a common theme: use observed data to make inferences about the world. Our work obtains the first provably efficient algorithms for Independent Component Analysis (ICA) in the presence of heavy-tailed data. The main tool in this result is the centroid body (a well-known topic in convex geometry), along with optimization and random walks for sampling from a convex body. This is the first algorithmic use of the centroid body and it is of independent theoretical…
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8
Figure 9
Figure 10
Figure 11
Figure 12
Figure 13
Figure 14
Figure 15
Figure 16
Figure 17
Figure 18
Figure 19
Figure 20
Figure 21
Figure 22
Figure 23
Figure 24Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsBlind Source Separation Techniques · Machine Learning and Algorithms · Sparse and Compressive Sensing Techniques
MethodsIndependent Component Analysis
\includeversion
vstd
\authordegrees
B.S, M.S \unitDepartment of Computer Science and Engineering
\coadvisorname
Luis Rademacher \coadvisornameAnastasios Sidiropolous \memberMikhail Belkin \memberFacundo Mémoli
Geometric Methods for Robust Data Analysis in High Dimension
Joseph Timothy Anderson
Abstract
Data-driven applications are growing. Machine learning and data analysis now finds both scientific and industrial application in biology, chemistry, geology, medicine, and physics. These applications rely on large quantities of data gathered from automated sensors and user input. Furthermore, the dimensionality of many datasets is extreme: more details are being gathered about single user interactions or sensor readings. All of these applications encounter problems with a common theme: use observed data to make inferences about the world. Our work obtains the first provably efficient algorithms for Independent Component Analysis (ICA) in the presence of heavy-tailed data. The main tool in this result is the centroid body (a well-known topic in convex geometry), along with optimization and random walks for sampling from a convex body. This is the first algorithmic use of the centroid body and it is of independent theoretical interest, since it effectively replaces the estimation of covariance from samples, and is more generally accessible.
We demonstrate that ICA is itself a powerful geometric primitive. That is, having access to an efficient algorithm for ICA enables us to efficiently solve other important problems in machine learning. The first such reduction is a solution to the open problem of efficiently learning the intersection of halfspaces in , posed in [43]. This reduction relies on a non-linear transformation of samples from such an intersection of halfspaces (i.e. a simplex) to samples which are approximately from a linearly transformed product distribution. Through this transformation of samples, which can be done efficiently, one can then use an ICA algorithm to recover the vertices of the intersection of halfspaces.
Finally, we again use ICA as an algorithmic primitive to construct an efficient solution to the widely-studied problem of learning the parameters of a Gaussian mixture model. Our algorithm again transforms samples from a Gaussian mixture model into samples which fit into the ICA model and, when processed by an ICA algorithm, result in recovery of the mixture parameters. Our algorithm is effective even when the number of Gaussians in the mixture grows with the ambient dimension, even polynomially in the dimension. In addition to the efficient parameter estimation, we also obtain a complexity lower bound for a low-dimension Gaussian mixture model.
\disscopyright
\dedication
*For my father,
who always believed that we can do better.*
Acknowledgements.
There are many who make this work possible. I thank my excellent faculty committee: Misha Belkin, Facundo Mémoli, Tasos Sidiropoulos, and Luis Rademacher. Much of this thesis was in collaboration with them, and I am in their debt. Special thanks go to Luis Rademacher for his years of advice, support, and teaching. The departments at Saint Vincent College and The Ohio State University have provided me with excellent environments for research and learning. Throughout my education I’ve had many helpful friends and collaborators; they each provided interesting and spirited discussion and insight. I am particularly grateful to Br. David Carlson for his instruction, advice, support, and friendship over the years. A large portion of my early adulthood was spent in his tutelage, and I am convinced nothing could have made it a better experience. I thank my family, especially my parents and sister, for their unfathomable support and encouragement. Finally, I thank my wife Sarah for her unfailing love and being an inspiring example of a smart, kind, and wonderful person. This work was supported by NSF grants CCF 1350870, and CCF 1422830.
{vita}
\dateitem
23 February 1990Born - Jeannette, PA
\dateitem
2012B.S. Computing & Information Science and Mathematics
Saint Vincent College
\dateitem
2016M.S. Computer Science & Engineering
The Ohio State University
\dateitem
2012-presentGraduate Teaching Associate
The Ohio State University.
{publist}
\researchpubs
\pubitem
J. Anderson, N. Goyal, A. Nandi, L. Rademacher. “Heavy-Tailed Analogues of the Covariance Matrix for ICA”. Thirty-First AAAI Conference on Artificial Intelligence (AAAI-17), Feb. 2017.
\pubitem
J. Anderson, N. Goyal, A. Nandi, L. Rademacher. “Heavy-Tailed Independent Component Analysis”. 56th Annual IEEE Symposium on Foundations of Computer Science (FOCS), IEEE 290-209, 2015.
\pubitem
J. Anderson, M. Belkin, N. Goyal, L. Rademacher, J. Voss. “The More The Merrier: The Blessing of Dimensionality for Learning Large Gaussian Mixtures”. Conference on Learning Theory (COLT), JMLR W&CP 35:1135-1164, 2014.
\pubitem
J. Anderson, N. Goyal, L. Rademacher. “Efficiently Learning Simplices”. Conference on Learning Theory (COLT), JMLR W&CP 30:1020-1045, 2013.
\pubitem
J. Anderson, M. Gundam, A. Joginipelly, D. Charalampidis. “FPGA implementation of graph cut based image thresholding”. Southeastern Symposium on System Theory (SSST), March. 2012.
{fieldsstudy}
\majorfield
Computer Science and Engineering
{studieslist}\studyitem
Theory of Computer ScienceLuis Rademacher \studyitemMachine LearningMikhail Belkin \studyitemMathematicsFacundo Mémoli
Contents
-
2.4 Orthogonalization via the uniform distribution in the centroid body
-
2.5.1 The fourth cumulant of Gaussian damping of heavy-tailed distributions
-
2.8.1 Heavy-tailed ICA when is orthogonal: Gaussian damping and experiments
-
3.3 Subroutine for finding the vertices of a rotated standard simplex
-
3.5 The local and global maxima of the 3rd moment of the standard simplex and the isotropic simplex
-
4.2 Learning GMM means using underdetermined ICA: The basic idea
List of Figures
- 2.1 The error of ICA with and without damping (solid lines and dashed lines, resp.), with unitary mixing matrix. The error is averaged over ten trials, in 3 and 10 dimensions where and , resp.
- 2.2 The difference between the errors of FastICA with and without damping in 2 dimensions, averaged over 40 trials. For a single cell, the parameters are given by the coordinates, . Red indicates that FastICA without damping does better than FastICA with damping, white indicates that the error difference is 0 and the blue indicates that FastICA with damping performs better than without damping. Black indicates that FastICA without damping failed (did not return two independent components).
- 2.3 The Frobenius error of the recovered mixing matrix with the ‘pow3’ and ‘tanh’ contrast functions, on 10-dimensional data, averaged over ten trials. The mixing matrix is random with unit norm columns, not orthogonal. In the left and middle figures, the distribution has while in the right figure, (see Section 2.8.2 for a discussion).
- 2.4 Microphone (blue circles) and human speaker (red “”) layouts for the “perimeterreg” (left) and “ceilreg” (right) voice separation data sets.
- 2.5 Stability parameter estimates for two voice recording data sets. Each line is one of the speech signals, we can see that the stability parameter is typically less than 2, and in the left data set, most of the signals have parameter close to 2, while in the right, the parameter is more commonly between 1.7 and 2.
- 2.6 (Left): Error of estimated mixing matrix on the “perimeterreg” data, averaged over ten trials.
- 2.7 1) Stability parameter estimates of each component in the “perimeterreg” data. Values below 2 are in the heavy-tailed regime. 2): Smallest singular value and condition number of the orthogonalization matrix computed via the centroid body and the covariance. The data was sampled with parameter .
Chapter 1 Introduction
The “curse of dimensionality” is a well-known problem encountered in the study of algorithms. This curse describes a phenomenon found in many fields of applied mathematics, for instance in numerical analysis, combinatorics, and computational geometry. In machine learning, many questions arise regarding how well existing algorithms scale as the dimension of the input data increases. These questions are well motivated in practice, as data includes more details about observed objects or events. The primary contribution of this research is to demonstrate that several important problems in machine learning are, in fact, efficiently solvable in high-dimension. This work outlines two main contributions toward this goal. First, we use the centroid body, from convex geometry, as a new algorithmic tool, giving the first algorithm for Independent Component Analysis which can tolerate heavy-tailed data. Second, we show that ICA itself can be used as an effective algorithmic primitive and that, with access to an efficient ICA algorithm, one can efficiently solve two other important open problems in machine learning: learning an intersection of halfspaces in , and learning the parameters of a Gaussian Mixture Model.
The fact that ICA can be used as a new algorithmic primitive, as we demonstrate, introduces new understanding of statistical data analysis. The second contribution in particular, where we show an efficient algorithm to learn the vertices of a simplex, has a key step in the reduction where one pre-processes the input in a carefully chosen – but simple – manner which results in data that will fit the standard ICA model. This pre-processing step is a random non-linear scaling inspired by the study of balls [16]. This scaling procedure is non-trivial and serves as an example of when one can exploit the geometry of a problem to gain statistical insight, and to frame the problem in a new light which has been thoroughly studied.
1.1 Organization of this Thesis
The rest of this thesis is organized as follows.
Chapter 2 will introduce the signal separation framework that will be of importance throughout this work: Independent Component Analysis (ICA). We will develop the general theory behind ICA and our method for improving the state-of-the art ICA method, extending the framework to be approachable when the input data has “heavy-tailed” properties. Our algorithm for ICA is theoretically sound, and comes with provable guarantees for polynomial time and sample complexity. We then present a more practical variation on our new ICA algorithm, and demonstrate its effectiveness on both synthetic and real-world heavy-tailed data. This chapter is based on [10] and [9].
Chapter 3 details an important first step in demonstrating that ICA can be used an effective algorithmic primitive. This chapter shows that ICA can be used to recover an arbitrary simplex with polynomial sample size. This chapter is based on work published in [11].
Chapter 4 presents a second reduction to ICA which recovers the parameters of a Gaussian Mixture Model efficiently in high dimension. Furthermore, we give a complexity lower bound for learning a Gaussian Mixture in low dimension which also yields a complexity lower bound for ICA itself in certain situations. This chapter is based on [8].
Chapter 2 Robust signal separation via convex geometry
The blind source separation problem is the general problem of recovering underlying “source signals” that have been mixed in some unknown way and are presented to an observer. Independent component analysis (ICA) is a popular model for blind source separation where the mixing is performed linearly. Formally, if is an -dimensional random vector from an unknown product distribution and is an invertible linear transformation, one is tasked with recovering the matrix and the signal , using only access to i.i.d. samples of the transformed signal, namely . Due to natural ambiguities, the recovery of is possible only up to the signs and permutations of the columns. Moreover, for the recovery to be possible the distributions of the random variables must not be a Gaussian distribution. ICA has applications in diverse areas such as neuroscience, signal processing, statistics, machine learning. There is vast literature on ICA; see, e.g., [31, 57, 33].
Since the formulation of the ICA model, a large number of algorithms have been devised employing a diverse set of techniques. Many of these existing algorithms break the problem into two phases: first, find a transformation which, when applied to the observed samples, gives a new distribution which is isotropic, i.e. a rotation of a (centered) product distribution; second, one typically uses an optimization procedure for a functional applied to the samples, such as the fourth directional moment, to recover the axes (or basis) of this product distribution.
To our knowledge, all known efficient algorithms for ICA with provable guarantees require higher moment assumptions such as finiteness of the fourth or higher moments for each component . Some of the most relevant works, e.g. algorithms of [39, 43], explicitly require the fourth moment to be finite. Algorithms in [111, 48], which make use of the characteristic function also seem to require at least the fourth moment to be finite: while the characteristic function exists for distributions without moments, the algorithms in these papers use the second or higher derivatives of the (second) characteristic function, and for this to be well-defined one needs the moments of that order to exist. Furthermore, certain anticoncentration properties of these derivatives are needed which require that fourth or higher moments exist.
Thus the following question arises: is ICA provably efficiently solvable when the moment condition is weakened so that, say, only the second moment exists, or even when no moments exist? By heavy-tailed ICA we mean the ICA problem with weak or no moment conditions (the precise moment conditions will be specified when needed).
Our focus in this chapter will be efficient algorithms for heavy-tailed ICA with provable guarantees and finite sample analysis. While we consider this problem to be interesting in its own right, it is also of interest in practice in a range of applications, e.g. [62, 64, 95, 28, 29, 91, 108]. The problem could also be interesting from the perspective of robust statistics because of the following informal connection: algorithms solving heavy-tailed ICA might work by focusing on samples in a small (but not low-probability) region in order to get reliable statistics about the data and ignore the long tail. Thus if the data for ICA is corrupted by outliers, the outliers are less likely to affect such an algorithm.
In this work, heavy-tailed distributions on the real line are those for which low order moments are not finite. Specifically, we will be interested in the case when the fourth or lower order moments are not finite as this is the case that is not covered by previous algorithms. We hasten to clarify that in some ICA literature the word heavy-tailed is used with a different and less standard meaning, namely distributions with positive kurtosis; this meaning will not be used in the present work.
Heavy-tailed distributions arise in a wide variety of contexts including signal processing and finance; see [78, 82] for an extensive bibliography. Some of the prominent examples of heavy-tailed distributions are the Pareto distribution with shape parameter which has moments of order less than , the Cauchy distributions, which has moments of order less than ; many more examples can be found on the Wikipedia page for heavy-tailed distributions. An abundant (and important in applications) supply of heavy-tailed distributions comes from stable distributions; see, e.g., [78]. There is also some theoretical work on learning mixtures of heavy-tailed distributions, e.g., [36, 27].
In several applied ICA models with heavy tails it is reasonable to assume that the distributions have finite first moment. In applications to finance (e.g., [30]), heavy tailed distributions are commonly used to model catastrophic but somewhat unlikely scenarios. A standard measures of risk in that literature, the so called conditional value at risk [87], is only finite when the first moment is finite. Therefore, it is reasonable to assume for some of our results that the distributions have finite first moment.
2.1 Main result
Our main result is an efficient algorithm that can recover the mixing matrix in the model when each has moments for a constant . The following theorem states more precisely the guarantees of our algorithm. The theorem below refers to the algorithm Fourier PCA [48] which solves ICA under the fourth moment assumption. The main reason to use this algorithm is that finite sample guarantees have been proved for it; we could have plugged in any other algorithm with such guarantee. The theorem below also refers to Gaussian damping, which is an algorithmic technique we introduce in this chapter and will be explained shortly.
Theorem 2.1.1** (Heavy-tailed ICA).**
Let be an ICA model such that the distribution of is absolutely continuous, for all we have and normalized so that , and the columns of have unit norm. Let be such that for each if has finite fourth moment then its fourth cumulant satisfies . Then, given , , , , Algorithm 1 combined with Gaussian damping and Fourier PCA outputs such that there are signs and a permutation satisfying with time and sample complexity and with probability at least . Here is a parameter of the distributions of the as described below. The degree of the polynomial is .
We note here that the assumption that has an absolutely continuous distribution is mostly for convenience in the analysis of Gaussian damping and not essential. In particular, it is not used in Algorithm 1.
Intuitively, in the theorem statement above measures how large a ball we need to restrict the distribution to, which has at least a constant (actually suffices) probability mass and, moreover, each when restricted to the interval has fourth cumulant at least . We show that all sufficiently large satisfy the above conditions and we can efficiently compute such an ; see the discussion after Theorem 2.1.2 (the restatement in Sec. 2.5). For standard heavy-tailed distributions, such as the Pareto distribution, behaves nicely. For example, consider the Pareto distribution with shape parameter and scale parameter , i.e. the distribution with density for and [math] otherwise. For this distribution it’s easily seen that suffices for the cumulant condition to be satisfied.
Theorem 2.1.1 requires that the -moment of the components be finite. However, if the matrix in the ICA model is unitary (i.e. , or in other words, is a rotation matrix) then we do not need any moment assumptions:
Theorem 2.1.2**.**
Let be an ICA model such that is unitary (i.e., ) and the distribution of is absolutely continuous. Let be such that for each if has finite fourth moment then . Then, given , Gaussian damping combined with Fourier PCA outputs such that there are signs and a permutation satisfying in time and sample complexity and with probability at least . Here is a parameter of the distributions of the as described above.
Idea of the algorithm.
Like many ICA algorithms, our algorithm has two phases: first orthogonalize the independent components (reduce to the pure rotation case), and then determine the rotation. In our heavy-tailed setting, each of these phases requires a novel approach and analysis in the heavy-tailed setting.
A standard orthogonalization algorithm is to put in isotropic position using the covariance matrix . This approach requires finite second moment of , which, in our setting, is not necessarily finite. Our orthogonalization algorithm (Section 2.4) only needs finite -absolute moment and that each is symmetrically distributed. (The symmetry condition is not needed for our ICA algorithm, as one can reduce the general case to the symmetric case, see Section 2.5.2). In order to understand the first absolute moment, it is helpful to look at certain convex bodies induced by the first and second moment. The directional second moment \mathbb{E}_{X}\bigl{(}(u^{T}X)^{2}\bigr{)} is a quadratic form in and its square root is the support function of a convex body, Legendre’s inertia ellipsoid, up to some scaling factor (see [73] for example). Similarly, one can show that the directional absolute first moment is the support function of a convex body, the centroid body. When the signals are symmetrically distributed, the centroid body of inherits these symmetries making it absolutely symmetric (see Section 2.2 for definitions) up to an affine transformation. In this case, a linear transformation that puts the centroid body in isotropic position also orthogonalizes the independent components (Lemma 10). In summary, the orthogonalization algorithm is the following: find a linear transformation that puts the centroid body of in isotropic position. One such matrix is given by the inverse of the square root of the covariance matrix of the uniform distribution in the centroid body. Then apply that transformation to to orthogonalize the independent components.
We now discuss how to determine the rotation (the second phase of our algorithm). The main idea is to reduce heavy-tailed case to a case where all moments exist and to use an existing ICA algorithm (from [48] in our case) to handle the resulting ICA instance. We use Gaussian damping to achieve such a reduction. By Gaussian damping we mean to multiply the density of the orthogonalized ICA model by a spherical Gaussian density.
We elaborate now on our contributions that make the algorithm possible.
Centroid body and orthogonalization.
The centroid body of a compact set was first defined in [81]. It is defined as the convex set whose support function equals the directional absolute first moment of the given compact set. We generalize the notion of centroid body to any probability measure having finite first moment (see Section 2.2 for the background on convexity and Section 2.2.5 for our formal definition of the centroid body for probability measures).
In order to put the centroid body in approximate isotropic position, we estimate its covariance matrix. For this, we use uniformly random samples from the centroid body. There are known methods to generate approximately random points from a convex body given by a membership oracle. We implement an efficient membership oracle for the centroid body of a probability measure with moments. The implementation works by first implementing a membership oracle for the polar of the centroid body via sampling and then using it via the ellipsoid method (see [50]) to construct a membership oracle for the centroid body. As far as we know this is the first use of the centroid body as an algorithmic tool.
An alternative approach to orthogonalization in ICA one might consider is to use the empirical covariance matrix of even when the distribution is heavy-tailed. A specific problem with this approach is that when the second moment does not exist, the diagonal entries would be very different and grow without bound. This problem gets worse when one collects more samples. This wide range of diagonal values makes the second phase of an ICA algorithm very unstable.
Linear equivariance and high symmetry.
A fundamental property of the centroid body, for our analysis, is that the centroid body is linearly equivariant, that is, if one applies an invertible linear transformation to a probability measure then the corresponding centroid body transforms in the same way (already observed in [81]). In a sense that we make precise (Lemma 9), high symmetry and linear equivariance of an object defined from a given probability measure are sufficient conditions to construct from such object a matrix that orthogonalizes the independent components of a given ICA model. This is another way to see the connection between the centroid body and Legendre’s ellipsoid of inertia for our purposes: Legendre’s ellipsoid of inertia of a distribution is linearly equivariant and has the required symmetries.
Gaussian damping.
Here we confine ourselves to the special case of ICA when the ICA matrix is unitary, that is . A natural idea to deal with heavy-tailed distributions is to truncate the distribution in far away regions and hope that the truncated distribution still gives us a way to extract information. In our setting, this could mean, for example, that we consider the random variable obtained from conditioned on the even that lies in the ball of radius centered at the origin. Instead of the ball we could restrict to other sets. Unfortunately, in general the resulting random variable does not come from an ICA model (i.e., does not have independent components in any basis). Nevertheless one may still be able to use this random variable for recovering . We do not know how to get an algorithmic handle on it even in the case of unitary . Intuitively, restricting to a set breaks the product structure of the distribution that is crucial for recovering the independent components.
We give a novel technique to solve heavy-tailed ICA for unitary . No moment assumptions on the components are needed for our technique. We call this technique Gaussian damping. Gaussian damping can also be thought of as restriction, but instead of being restriction to a set it is a “restriction to a spherical Gaussian distribution.” Let us explain. Suppose we have a distribution on with density . If we restrict this distribution to a set (which we assume to be nice: full-dimensional and without any measure theoretic issues) then the density of the restricted distribution is [math] outside and is proportional to for . One can also think of the density of the restriction as being proportional to the product of and the density of the uniform distribution on . In the same vein, by restriction to the Gaussian distribution with density proportional to we simply mean the distribution with density proportional to . In other words, the density of the restriction is obtained by multiplying the two densities. By Gaussian damping of a distribution we mean the distribution obtained by this operation.
Gaussian damping provides a tool to solve the ICA problem for unitary by virtue of the following properties: (1) The damped distribution has finite moments of all orders. This is an easy consequence of the fact that Gaussian density decreases super-polynomially. More precisely, one dimensional moment of order given by the integral is finite for all for any distribution. (2) Gaussian damping retains the product structure. Here we use the property of spherical Gaussians that it’s the (unique) class of spherically symmetric distributions with independent components, i.e., the density factors: (we are hiding a normalizing constant factor). Hence the damped density also factors when expressed in terms of the components of (again ignoring normalizing constant factors):
[TABLE]
Thus we have converted our heavy-tailed ICA model into another ICA model where and are obtained by Gaussian damping of and , resp. To this new model we can apply the existing ICA algorithms which require at least the fourth moment to exist. This allows us to estimate matrix (up to signs and permutations of the columns).
It remains to explain how we get access to the damped random variable . This is done by a simple rejection sampling procedure. Damping does not come free and one has to pay for it in terms of higher sample and computational complexity, but this increase in complexity is mild in the sense that the dependence on various parameters of the problem is still of similar nature as for the non-heavy-tailed case. A new parameter is introduced here which parameterizes the Gaussian distribution used for damping. We explained the intuitive meaning of after the statement of Theorem 2.1.1 and that it’s a well-behaved quantity for standard distributions.
Gaussian damping as contrast function.
Another way to view Gaussian damping is in terms of contrast functions, a general idea that in particular has been used fruitfully in the ICA literature. Briefly, given a function , for on the unit sphere in , we compute . Now the properties of the function as varies over the unit sphere, such as its local extrema, can help us infer properties of the underlying distribution. In particular, one can solve the ICA problem for appropriately chosen contrast function . In ICA, algorithms with provable guarantees use contrast functions such as moments or cumulants (e.g., [39, 43]). Many other contrast functions are also used.
Gaussian damping furnishes a novel class of contrast functions that also leads to provable guarantees. E.g., the function given by is in this class. We do not use the contrast function view in this work.
Previous work related to damping.
To our knowledge the damping technique, and more generally the idea of reweighting the data, is new in the context of ICA. But the general idea of reweighting is not new in other somewhat related contexts as we now discuss. In robust statistics (see, e.g., [55]), reweighting idea is used for outlier removal by giving less weight to far away data points. Apart from this high-level similarity we are not aware of any closer connections to our setting; in particular, the weights used are generally different.
Another related work is [22], on isotropic PCA, affine invariant clustering, and learning mixtures of Gaussians. This work uses Gaussian reweighting. However, again we are unaware of any more specific connection to our problem.
Finally, in [111, 48] a different reweighting, using a “Fourier weight” (here is a fixed vector and is a data point) is used in the computation of the covariance matrix. This covariance matrix is useful for solving the ICA problem. But as discussed before, results here do not seem to be amenable to our heavy-tailed setting.
2.2 Preliminaries
In this section, we review technical definitions and results which will be used in subsequent sections.
For a random vector , we denote the distribution function of as and, if it exists, the density is denoted .
For a real-value random variable , cumulants of are polynomials in the moments of . For , the th cumulant is denoted . Denoting . examples: . In general, cumulants can be defined as the coefficients in the logarithm of the moment generating function of :
[TABLE]
The first two cumulants are the same as the expectation and the variance, resp. Cumulants have the property that for two independent r.v.s we have (assuming that the first moments exist for both and ). Cumulants are order- homogeneous, i.e. if and is a r.v., then . The first two cumulants of the standard Gaussian distribution are the mean [math] and the variance , and all subsequent Guassian cumulants have value [math].
An -dimensional convex body is a compact convex subset of with non-empty interior. We say a convex body is absolutely symmetric if . Similarly, we say random variable (and its distribution) is absolutely symmetric if, for any choice of signs , has the same distribution as . We say that an -dimensional random vector is symmetric if has the same distribution as . Note that if is symmetric with independent components (mutually independent coordinates) then its components are also symmetric.
For a convex body , the Minkowski functional of is . The Minkowski sum of two sets is the set .
The singular values of a matrix will be ordered in the decreasing order: . By we mean .
We say that a matrix is unitary if , or in other words is a rotation matrix. (Normally for matrices with real-valued entries such matrices are called orthogonal matrices and the word unitary is reserved for their complex counterparts, however the word orthogonal matrix can lead to confusion in the present work.) For matrix , denote by the spectral norm and by the Frobenius norm. We will need the following inequality about the stability of matrix inversion (see for example [102, Chapter III, Theorem 2.5]).
Lemma 1**.**
Let be a matrix norm such that . Let matrices be such that , and let . Then
[TABLE]
This implies that if , then
[TABLE]
2.2.1 Heavy-Tailed distributions
In general, heavy-tailed distributions are those whose tail probabilities go to zero more slowly than an inverse exponential. Formally, one may write that the distribution of is heavy-tailed if diverges as approaches infinity. However, for our applications, we need only care about the existence of higher moments of the distribution. The primary statistical assumption made in the following work assumes only that for some . The addition of 1 in the above is primarily so that means are well-defined. Naturally, from this assumption, one can construct heavy-tailed distributions, e.g. with density proportional to , which will have first moment, up to the moment, but no higher.
2.2.2 Stable Distributions
The Stable Distributions are an important class of probability distributions characterized by one key property: the family itself is closed under addition and scaler multiplication. That is, if you have two random variables and with the same stable distribution, will also be a stable distribution with slightly different parameters. Two well-known members of this family are the Gaussian and Cauchy distribution, both of which will be important throughout this work.
Stable distributions are fully characterized by four parameters and can be written as where is the stable parameter ( for Gaussian and for Cauchy), is a skewness parameter, is scale, and is location. The PDF and CDF are not expressible analytically in general, but stable distributions are those for which the characteristic function can be written as
[TABLE]
where is the sign function and
[TABLE]
A useful property which will be used later, and gives a formalization of the closure property mentioned above is the following: if are iid with stable density and we let
[TABLE]
where , then will have density
[TABLE]
where
[TABLE]
2.2.3 Convex optimization
We need the following result: Given a membership oracle for a convex body one can implement efficiently a membership oracle for , the polar of . This follows from applications of the ellipsoid method from [50]. Specifically, we use the following facts: (1) a validity oracle for can be constructed from a membership oracle for [50, Theorem 4.3.2]; (2) a membership oracle for can be constructed from a validity oracle for [50, Theorem 4.4.1].
The definitions and theorems in this section all come (occasionally with slight rephrasing) from [50] except for the notion of -weak oracle. [The definitions below use rational numbers instead of real numbers. This is done in [50] as they work out in detail the important low level issues of how the numbers in the algorithm are represented as general real numbers cannot be directly handled by computers. These low-level details can also be worked out for the arguments here, but as is customary, we will not describe these and use real numbers for the sake of exposition.] In this section is a convex body. For , the distance of to is given by . Define and . Let denote the set of rational numbers.
Definition 1** ([50]).**
The -weak membership problem for is the following: Given a point and a rational number , either (i) assert that , or (ii) assert that . An -weak membership oracle for is an oracle that solves the weak membership problem for . For , an -weak membership oracle for acts as follows: Given a point , with probability at least it solves the -weak membership problem for , and otherwise its output can be arbitrary.
Definition 2** ([50]).**
The -weak validity problem for is the following: Given a vector , a rational number , and a rational number , either (i) assert that for all , or (ii) assert that for some . The notion of -weak validity oracle and -weak validity oracle can be defined similarly to Def. 1.
Definition 3** ([50, Section 2.1]).**
We say that an oracle algorithm is an oracle-polynomial time algorithm for a certain problem defined on a class of convex sets if the running time of the algorithm is bounded by a polynomial in the encoding length of and in the encoding length of the possibly existing further input, for every convex set in the given class.
The encoding length of a convex set will be specified below, depending on how the convex set is presented.
Theorem 2.2.1** (Theorem 4.3.2 in [50]).**
Let and . There exists an oracle-polynomial time algorithm that solves the weak validity problem for every convex body contained in the ball of radius and containing a ball of radius centered at given by a weak membership oracle. The encoding length of is plus the length of the binary encoding of , and .
We remark that the Theorem 4.3.2 as stated in [50] is stronger than the above statement in that it constructs a weak violation oracle (not defined here) which gives a weak validity oracle which suffices for us. The algorithm given by Theorem 4.3.2 makes a polynomial (in the encoding length of ) number of queries to the weak membership oracle.
Lemma 2** (Lemma 4.4.1 in [50]).**
There exists an oracle-polynomial time algorithm that solves the weak membership problem for , where is a convex body contained in the ball of radius and containing a ball of radius centered at [math] given by a weak validity oracle. The encoding length of is plus the length of the binary encoding of and .
Our algorithms and proofs will need more quantitative details from the proofs of the above theorem and lemma. These will be mentioned when we need them.
2.2.4 Algorithmic convexity
We state here a special case of a standard result in algorithmic convexity: There is an efficient algorithm to estimate the covariance matrix of the uniform distribution in a centrally symmetric convex body given by a weak membership oracle.
The result follows from the random walk-based algorithms to generate approximately uniformly random points from a convex body [41, 60, 70]. Most papers use access to a membership oracle for the given convex body. In this work we only have access to an -weak membership oracle. As discussed in [41, Section 6, Remark 2], essentially the same algorithm implements efficient sampling when given an -weak membership oracle. The problem of estimating the covariance matrix of a convex body was introduced in [60, Section 5.2]. That paper analyzes the estimation from random points. There has been a sequence of papers studying the sample complexity of this problem [21, 89, 45, 80, 15, 107, 2, 101].
Theorem 2.2.2**.**
Let be a centrally symmetric convex body given by a weak membership oracle so that . Let . Then there exists a randomized algorithm that, when given access to the weak membership of and inputs , it outputs a matrix such that with probability at least over the randomness of the algorithm, we have
[TABLE]
The running time of the algorithm is .
Note that (2.3) implies . The guarantee in (2.3) has the advantage of being invariant under linear transformations in the following sense: if one applies an invertible linear transformation to the underlying convex body, the covariance matrix and its estimate become and , respectively. These matrices satisfy
[TABLE]
This fact will be used later.
2.2.5 The centroid body
The main tool in our orthogonalization algorithm is the centroid body of a distribution, which we use as a first moment analogue of the covariance matrix. In convex geometry, the centroid body is a standard (see, e.g., [81, 73, 44]) convex body associated to (the uniform distribution on) a given convex body. Here we use a generalization of the definition from the case of the uniform distribution on a convex body to more general probability measures. Let be a random vector with finite first moment, that is, for all we have . Following [81], consider the function . Then it is easy to see that , is positively homogeneous, and is subadditive. Therefore, it is the support function of a compact convex set [81, Section 3], [93, Theorem 1.7.1], [44, Section 0.6]. This justifies the following definition:
Definition 4** (Centroid body).**
Let be a random vector with finite first moment, that is, for all we have . The centroid body of is the compact convex set, denoted , whose support function is . For a probability measure , we define , the centroid body of , as the centroid body of any random vector distributed according to .
The following lemma says that the centroid body is equivariant under linear transformations. It is a slight generalization of statements in [81] and [44, Theorem 9.1.3].
Lemma 3**.**
Let be a random vector on . Let be an invertible linear transformation. Then .
Proof.
. ∎
Lemma 4** ([44, Section 0.8], [93, Remark 1.7.7]).**
Let be a convex body with support function and such that the origin is in the interior of . Then is a convex body and has radial function for .
Lemma 4 implies that testing membership in is a one-dimensional problem if we have access to the support function : We can decide if a point is in by testing if instead of needing to check for all . In our application . We can estimate by taking the empirical average of where the are samples of . This leads to an approximate oracle for which will suffice for our application. The details are in Sec. 2.3.
2.3 Membership oracle for the centroid body
In this section we provide an efficient weak membership oracle (Subroutine 2) for the centroid body of the r.v. . This is done by first providing a weak membership oracle (Subroutine 1) for the polar body . We begin with a lemma that shows that under certain general conditions the centroid body is “well-rounded.” This property will prove useful in the membership tests.
Lemma 5**.**
Let be an absolutely symmetrically distributed random vector such that for all . Then . Moreover, .
Proof.
The support function of is (Def. 4). Then, for each canonical vector , . Thus, is contained in . Moreover, since , we get .
We claim now that each canonical vector is contained in . To see why, first note that since the support function is along each canonical direction, will touch the facets of the unit hypercube. Say, there is a point that touches facet associated to canonical vector . But the symmetry of the s implies that is absolutely symmetric, so that is also in the centroid body. Convexity implies that is in the centroid body. The same argument applied to all canonical vectors implies that they are all contained in the centroid body, and this with convexity implies that the centroid body contains . In particular, it contains . ∎
2.3.1 Mean estimation using moments
We will need to estimate the support function of the centroid body in various directions. To this end we need to estimate the first absolute moment of the projection to a direction. Our assumption that each component has finite -moment will allow us to do this with a reasonable small probability of error. This is done via the following Chebyshev-type inequality.
Let be a real-valued symmetric random variable such that for some and . Then we will prove that the empirical average of the expectation of converges to the expectation of .
Let be the empirical average obtained from independent samples , i.e., .
Lemma 6**.**
Let . With the notation above, for , we have
[TABLE]
Proof.
Let be a threshold whose precise value we will choose later. We have
[TABLE]
By the union bound,
[TABLE]
Define a new random variable by
[TABLE]
Using the symmetry of we have
[TABLE]
By the Chebyshev inequality and (2.5) we get
[TABLE]
Putting (2.4) and (2.6) together for we get
[TABLE]
Choosing
[TABLE]
the RHS of the previous equation becomes . The choice of is made to minimize the RHS; we ignore integrality issues. Pick so that . To estimate , note that
[TABLE]
Hence
[TABLE]
We want which is equivalent to . We set . Then, for putting together (2.7) and (2.9) gives
[TABLE]
Setting (so that ) and expressing the RHS of the last equation in terms of via (2.8) (and eliminating ), and using our assumptions , to get a simpler upper bound, we get
[TABLE]
Condition , when expressed in terms of via (2.8), becomes
[TABLE]
∎
2.3.2 Membership oracle for the polar of the centroid body
As mentioned before, our membership oracle for (Subroutine 1) is based on the fact that is the radial function of , and that is the directional absolute first moment of , which can be efficiently estimated by sampling.
Lemma 7** (Correctness of Subroutine 1).**
Let be a constant and be given by a symmetric ICA model such that for all we have and normalized so that . Let . Given , , Subroutine 1 is an -weak membership oracle for with using time and sample complexity . The degree of the polynomial is .
Proof.
Recall from Def. 1 that we need to show that, with probability at least , Subroutine 1 outputs TRUE when and FALSE when ; otherwise, the output can be either TRUE or FALSE arbitrarily.
Fix a point and let denote the direction of . The algorithm estimates the radial function of along , which is (see Lemma 4). In the following computation, we simplify the notation by using . It is enough to show that with probability at least the algorithm’s estimate, , of the radial function is within of the true value, .
For , i.i.d. copies of , the empirical estimator for is .
We want to apply Lemma 6 to . For this we need a bound on its -moment. The following simple bound is sufficient for our purposes: Let . Then for all . Then
[TABLE]
Lemma 6 implies that, for to be fixed later, and for
[TABLE]
we have .
Lemmas 5 and 3 give that for . It follows that . If and , then we have
[TABLE]
which in turn gives, when ,
[TABLE]
Plugging in the value of and, in turn of , into (2.11) gives that it suffices to take
[TABLE]
∎
2.3.3 Membership oracle for the centroid body
We now describe how the weak membership oracle for the centroid body is constructed using the weak membership oracle for , provided by Subroutine 1.
We will use the following notation: For a convex body , , such that , oracle is an -weak membership oracle for . Similarly, oracle is an -weak validity oracle. Lemma 5 along with the equivariance of (Lemma 3) gives . Then . Set and .
Detailed description of Subroutine 2.
There are two main steps:
Use Subroutine 1 to get an -weak membership oracle for . Theorem 4.3.2 of [50] (stated as Theorem 2.2.1 here) is used in Lemma 8 to get an algorithm to implement an -weak validity oracle running in oracle polynomial time; invokes a polynomial number of times, specifically (see proof of Lemma 8). The proof of Theorem 4.3.2 can be modified so that . 2. 2.
Lemma 4.4.1 of [50] (Lemma 2 here) gives an algorithm to construct an -weak membership oracle from . Lemma 4.4.1 in [50] shows calls once, with (where is the query point).
Lemma 8** (Correctness of Subroutine 2).**
Let be given by a symmetric ICA model such that for all we have and normalized so that . Then, given a query point , , , , and , Subroutine 2 is an -weak membership oracle for and with probability using time and sample complexity
Proof.
We first prove that (abbreviated to hereafter) works correctly. To this end we need to show that for any given input, acts as an -weak validity oracle with probability at least . Oracle makes queries to . If the answer to all these queries were correct then Theorem 4.3.2 from [50] would apply and would give that outputs an answer as expected. Since these queries are adaptive we cannot directly apply the union bound to say that the probability of all of them being correct is at least . However, a more careful bound allows us to do essentially that.
Let be the sequence of queries, where depends on the result of the previous queries. For , let be the event that the answer to query by Subroutine 1 is not correct according to the definition of the oracle it implements. These events are over the randomness of Subroutine 1 and event involves the randomness of , as the queries could be adaptively chosen. By the union bound, the probability that all answers are correct is at least . It is enough to show that . To see this, we can condition on the randomness associated to . That makes deterministic, and the probability of failure is now just the probability that Subroutine 1 fails. More precisely, , so that
[TABLE]
This proves that the first step works correctly. Correctness of the second step follows directly because the algorithm for construction of the oracle involves a single call to the input oracle as mentioned in Step 2 of the detailed description.
Finally, to prove that the running time of Subroutine 2 is as claimed the main thing to note is that, as mentioned in Step 1 of the detailed description, is polynomially small in and is polynomially small in and so is polynomially small in .
∎
2.4 Orthogonalization via the uniform distribution in the centroid body
The following lemma says that linear equivariance allows orthogonalization:
Lemma 9**.**
Let be a family of -dimensional product distributions. Let be the closure of under invertible linear transformations. Let be an -dimensional distribution defined as a function of . Assume that and satisfy:
For all , is absolutely symmetric. 2. 2.
* is linear equivariant (that is, for any invertible linear transformation we have ).* 3. 3.
For any , is positive definite.
Then for any symmetric ICA model with we have is an orthogonalizer of .
Proof.
Consider a symmetric ICA model with . Assumptions 1 and 3 imply is diagonal and positive definite. This with Assumption 2 gives . Let (the unique symmetric positive definite square root). We have for some unitary matrix (see [53, pg 406]).
Thus, has orthogonal columns, that is, it is an orthogonalizer for . ∎
The following lemma applies the previous lemma to the special case when the distribution is the uniform distribution on .
Lemma 10**.**
Let be a random vector drawn from a symmetric ICA model such that for all we have . Let Y be uniformly random in . Then is an orthogonalizer of .
Proof.
We will use Lemma 9. After a scaling of each , we can assume without loss of generality that . This will allow us to use Lemma 5. Let = { is an absolutely symmetric product distribution and , for all }. For , let be the uniform distribution on the centroid body of . For all , the symmetry of the ’s implies that , that is, , is absolutely symmetric. By the equivariance of (from Lemma 3) and Lemma 5 it follows that is linear equivariant. Let . Then there exist and such that . So we get . From Lemma 5 we know so that is a diagonal matrix with positive diagonal entries. This implies that is positive definite and thus by Lemma 9, is an orthogonalizer of . ∎
Theorem 2.4.1** (Correctness of algorithm 1).**
Let be given by a symmetric ICA model such that for all we have and normalized so that . Then, given , , , , Algorithm 1 outputs a matrix so that , for a diagonal matrix with diagonal entries satisfying . with probability at least using time and sample complexity.
Proof.
From Lemma 5 we know . Using the equivariance of (Lemma 3), we get . Thus, to satisfy the roundness condition of Theorem 2.2.2 we can take , .
Let be the estimate of computed by the algorithm. Let be the estimate of obtained from according to how covariance matrices transform under invertible linear transformations of the underlying random vector. As in the proof of Lemma 9, we have and for some unitary matrix . Thus, we have . It is natural then to set . Let be the diagonal entries of . We have, using Lemma 1,
[TABLE]
As in (2.2), we show that is small:
We first bound , the diagonal entries of . Let and . We find simple estimates of these quantities: We have and is the maximum variance of along coordinate axes. From Lemma 5 we know , so that and . Similarly, , where is the smallest diagonal entry of . In other words, it is the minimum variance of along coordinate axes. From Lemma 5 we know , so that for all and Lemma 11 below implies . That is, .
From the bounds on , Theorem 2.2.2 and the fact discussed after it, we have
[TABLE]
when .
This in (2.12) with Theorem 2.2.2 again gives
[TABLE]
The claim follows by setting . The sample and time complexity of the algorithm comes from the calls to Subroutine 2. The number of calls is given by Theorem 2.2.2. This leads to the complexity as claimed. ∎
Lemma 11**.**
Let be an absolutely symmetric convex body such that contains the segment (where is the first canonical vector). Let be uniformly random in . Then .
Proof.
Let be a diagonal linear transformation so that isotropic. Let be the first entry of . It is known that any -dimensional isotropic convex body is contained in the ball of radius [73, 99],[59, Theorem 4.1]. Note that contains the segment . This implies . Also, by isotropy we have, . The claim follows. ∎
2.5 Gaussian damping
In this section we give an efficient algorithm for the heavy-tailed ICA problem when the ICA matrix is a unitary matrix; no assumptions on the existence of moments of the will be required.
The basic idea behind our algorithm is simple and intuitive: using we construct another ICA model , where is a parameter which will be chosen later. The components of have light-tailed distributions; in particular, all moments exist. We show how to generate samples of efficiently using samples of . Using the new ICA model, the matrix can be estimated by applying existing ICA algorithms.
For a random variable we will denote its the probability density function by . The density of is obtained by multiplying the density of by a Gaussian damping factor. More precisely,
[TABLE]
Define
[TABLE]
then
[TABLE]
We will now find the density of . Note that if is a value of , and is the corresponding value of , then we have
[TABLE]
where we used that is a unitary matrix so that . Also, follows from the change of variable formula and the fact that . We also used crucially the fact that the Gaussian distribution is spherically-symmetric. We have now specified the new ICA model , and what remains is to show how to generate samples of .
Rejection sampling.
Given access to samples from we will use rejection sampling (see e.g. [Robert–Casella] ) to generate samples from .
Generate . 2. 2.
Generate . 3. 3.
If , output ; else, go to the first step.
The probability of outputting a sample with a single trial in the above algorithm is . Thus, the expected number of trials in the above algorithm for generating a sample is .
We now choose . There are two properties that we want sufficiently large so as to satisfy: (1) and where and are constants. Such a choice of exists and can be made efficiently; we outline this after the statement of Theorem 2.1.2. Thus, the expected number of trials in rejection sampling before generating a sample is bounded above by . The lower bound on will also be useful in bounding the moments of the , where is the random variable obtained by Gaussian damping of with parameter , that is to say
[TABLE]
Define
[TABLE]
and let be the product of over . By we denote the vector with its th element removed, then notice that
[TABLE]
We can express the densities of individual components of as follows:
[TABLE]
This allows us to derive bounds on the moments of :
[TABLE]
We now state Theorem 4.2 from [48] in a special case by setting parameters and in that theorem to for . The algorithm analyzed in Theorem 4.2 of [48] is called Fourier PCA.
Theorem 2.5.1**.**
[48]** Let be given by an ICA model where is unitary and the are mutually independent, for some positive constant , and . For any with probability at least , Fourier PCA will recover vectors such that there exist signs and a permutation satisfying
[TABLE]
using samples. The running time of the algorithm is also of the same form.
Combining the above theorem with Gaussian damping gives the following theorem. As previously noted, since we are doing rejection sampling in Gaussian damping, the expected number of trials to generate samples of is . One can similarly prove high probability guarantees for the number of trials needed to generate samples.
We remark that the choice of in the above theorem can be made algorithmically in an efficient way. Theorem 2.5.2 below shows that as we increase the cumulant goes to infinty. This shows that for any there exists so as to satisfy the condition of the above theorem, namely . We now briefly indicate how such an can be found efficiently (same sample and computational costs as in Theorem 2.1.2 above): For a given , we can certainly estimate from samples, i.e. by the empirical mean of samples . This allows us to search for so that is as large as we want. This also gives us an upper bound on the fourth moment via Eq. (2.15). To ensure that the fourth cumulants of all are large, note that for we have . We can estimate this quantity empirically, and minimize over on the unit sphere (the minimization can be done, e.g., using the algorithm in [43]). This would give an estimate of and allows us to search for an appropriate .
For the algorithm to be efficient, we also need . This is easily achieved as we can empirically estimate using the number of trials required in rejection sampling, and search for sufficiently large that makes the estimate sufficiently larger than .
2.5.1 The fourth cumulant of Gaussian damping of heavy-tailed distributions
It is clear that if r.v. is such that and , then as . However, it does not seem clear when we have as well. We will show that in this case we also get as .
We will confine our discussion to symmetric random variables for simplicity of exposition; for the purpose of our application of the theorem this is w.l.o.g. by the argument in Sec. 2.5.2.
Theorem 2.5.2**.**
Let be a symmetric real-valued random variable with . Then as .
Proof.
Fix a symmetric r.v. with ; as previously noted, if then the theorem is easily seen to be true. Since is symmetric and we will be interested in the fourth cumulant, we can restrict our attention to the positive part of . So in the following we will actually assume that is a positive random variable. Fix to be any large positive constant. Fix a small positive constant . Also fix another small positive constant . Then there exists such that
[TABLE]
Let . Recall that . Note that if (which we assume in the sequel), then
[TABLE]
Since as , by choosing sufficiently large we can ensure that
[TABLE]
Moreover, we choose to be sufficiently large so that . Then
[TABLE]
Summarizing the previous sequence of inequalities:
[TABLE]
Now
[TABLE]
Now note that for our choice of the parameters. Thus , and by our assumption with . ∎
2.5.2 Symmetrization
As usual we work with the ICA model . Suppose that we have an ICA algorithm that works when each of the component random variable is symmetric, i.e. its probability density function satisfies for all , with a polynomial dependence on the upper bound on the fourth moment of the and inverse polynomial dependence on the lower bound on the fourth cumulants of the . Then we show that we also have an algorithm without the symmetry assumption and with a similar dependence on and . We show that without loss of generality we may restrict our attention to symmetric densities, i.e. we can assume that each of the has density function satisfying . To this end, let be an independent copy of and set . Similarly, let . Clearly, the and have symmetric densities. The new random variables still satisfy the ICA model: . Moreover, the moments and cumulants of the behave similarly to those of the : For the fourth moment, assuming it exists, we have . The inequality above can easily be proved using the binomial expansion and Hölder’s inequality:
[TABLE]
The final inequality follows from the fact that for each term in the LHS, e.g. we have .
For the fourth cumulant, again assuming its existence, we have .
Thus if the fourth cumulant of is away from [math] then so is the fourth cumulant of .
2.6 Putting things together
In this section we combine the orthogonalization procedure (Algorithm 1 with performance guarantees in Theorem 2.4.1) with ICA for unitary via Gaussian damping to prove our main theorem, Theorem 2.1.1.
As noted in the introduction, intuitively, in the theorem statement above measures how large a ball we need to restrict the distribution to so that there is at least a constant (or if needed) probability mass in it and moreover each when restricted to the interval has the fourth cumulant at least . Formally, is such that , where can be chosen, and for simplicity, we will fix to . Moreover, satisfies that for all and , where is the Gaussian damping with parameter of .
In Sec. 2.5.2 we saw that the moments and cumulants of the non-symmetric random variable behave similarly to those of the symmetric random variable. So by the argument of Sec. 2.5.2 we assume that our ICA model is symmetric. Theorem 2.4.1 shows that Algorithm 1 gives us a new ICA model with the ICA matrix having approximately orthogonal columns. We will apply Gaussian damping to this new ICA model. In Theorem 2.4.1, it was convenient to use the normalization for all . But for the next step of Gaussian damping we will use a different normalization, namely, the columns of the ICA matrix have unit length. This will require us to rescale (in the analysis, not algorithmically) the components of appropriately as we now describe.
Algorithm 1 provides us with a matrix such that the columns of are approximately orthogonal: where is a diagonal matrix. Thus, we can rewrite our ICA model as , where . We rescale , the th column of , by multiplying it by . Denoting by the diagonal matrix with the th diagonal entry , the matrix obtained after the above rescaling of is and we have . We can again rewrite our ICA model as . Setting and we can rewrite our ICA model as . This is the model we will plug into the Gaussian damping procedure. Had Algorithm 1 provided us with perfect orthogonalizer (so that ) we would obtain a model where is unitary. We do get however that . To continue with a more standard ICA notation, from here on we will write for and for .
Applying Gaussian damping to gives us a new ICA model as we saw in Sec. 2.5. But the model we have access to is . We will apply Gaussian damping to it to get the r.v. . Formally, is defined starting with the model just as we defined starting with the model (recall that for a random variable , we denote its probability density function by ):
[TABLE]
where . The parameter has been chosen so that and for all . By the discussion after Theorem 2.1.2 (the restatement in Sec. 2.5), this choice of can be made efficiently. (The discussion there is in terms of the directional moments of , but note that the directional moments of also give us directional moments of . We omit further details.) But now since the matrix in our ICA model is only approximately unitary, after applying Gaussian damping the obtained random variable is not given by an ICA model (in particular, it may not have independent coordinates in any basis), although it is close to in a sense to be made precise soon.
Because of this, Theorem 2.5.1 is not directly usable for plugging in the samples of . To address this discrepancy we will need a robust version of Theorem 2.5.1 which also requires us to specify in a precise sense that and are close. To this end, we need some standard terminology from probability theory. The characteristic function of r.v. is defined to be , where . The cumulant generating function, also known as the second characteristic function, is defined by . The algorithm in [48] estimates the second derivative of and computes its eigendecomposition. (In [111] and [48], this second derivative is interpreted as a kind of covariance matrix of but with the twist that a certain “Fourier” weight is used in the expectation computation for the covariance matrix. We will not use this interpretation here.) Set , the Hessian matrix of . We can now state the robust version of Theorem 2.5.1.
Theorem 2.6.1**.**
Let be an -dimensional random vector given by an ICA model where is unitary and the are mutually independent, and for positive constants and . Also let . Suppose that we have another random variable that is close to in the following sense:
[TABLE]
for any with . Moreover, for . When Fourier PCA is given samples of it will recover vectors such that there exist signs and a permutation satisfying
[TABLE]
in samples and time complexity and with probability at least .
While this theorem is not stated in [48], it is easy to derive from their proof of Theorem 2.5.1; we now briefly sketch the proof of Theorem 2.6.1 indicating the changes one needs to make to the proof of Theorem 2.5.1 in [48].
Proof.
Ideally, for input model with unitary, algorithm Fourier PCA would proceed by diagonalizing . But it can only compute an approximation which is the empirical estimate for . For all with , it is shown that with high probability we have
[TABLE]
Then, a matrix perturbation argument is invoked to show that if the diagonalization procedure used in Fourier PCA is applied to instead of , one still recovers a good approximation of . This previous step uses a random chosen from a Gaussian distribution so that the eigenvalues of are sufficiently spaced apart for the eigenvectors to be recoverable (the assumptions on the distribution ensure that the requirement of is satisfied with high probability). The only property of used in this argument is (2.20). To prove Theorem 2.6.1, we show that the estimate is also good:
[TABLE]
where we ensured that by taking sufficiently many samples of to get a good estimate with probability at least ; as in [48], a standard concentration argument shows that samples suffice for this purpose. Thus the diagonalization procedure can be applied to . The upper bound of above translates into error in the final recovery guarantee, with the extra factor coming from the eigenvalue gaps of . ∎
To apply Theorem 2.6.1 to our situation, we need
[TABLE]
This will follow from the next lemma.
Note that
[TABLE]
(The gradient is a row vector.) Thus, to show (2.21) it suffices to show that each expression on the RHS of the previous equation is appropriately close:
Lemma 12**.**
Let , and let such that is unitary and . Let and be the random variables obtained by applying Gaussian damping to the ICA models and , resp. Then, for , we have
[TABLE]
Proof.
We will only prove the first inequality; proofs of the other two are very similar and will be omitted. In the second equality in the displayed equations below we use that . One way to see this is to think of the two integrals as expectations: .
[TABLE]
Now
[TABLE]
We have
[TABLE]
and so
[TABLE]
For the first summand in (2.24) note that
[TABLE]
This follows from the fact that for real we have
[TABLE]
So,
[TABLE]
where the last inequality used our assumption that .
We will next bound second summand in (2.24). Note that and . Since , we get . We will use that for which is satisfied by our assumption. Now the second summand in (2.24) can be bounded as follows.
[TABLE]
Combining our estimates gives
[TABLE]
Finally, to bound , note that
[TABLE]
This we just upper-bounded above by .
Thus we have the final estimate
[TABLE]
The proofs of the other two upper bounds in the lemma follow the same general pattern with slight changes. ∎
We are now ready to prove Theorem 2.1.1.
Proof of Theorem 2.1.1.
We continue with the context set after the statement of Theorem 2.1.1. The plan is to apply Theorem 2.6.1 to and . To this end we begin by showing that the premise of Theorem 2.6.1 is satisfied.
Theorem 2.4.1, with and to be specified, provides us with a matrix such that the columns of are approximately orthogonal: for some diagonal matrix . Now set , where . Theorem 2.4.1 implies that
[TABLE]
Then
[TABLE]
because by (2.26). For as above, there exists a unitary such that
[TABLE]
by Lemma 13 below.
By our choice of , the components of satisfy and (via (2.15) and our choice ). Hence . The latter bound via (2.27) gives for all .
Finally, Lemma 12 with (2.22) and simple estimates give
We are now ready to apply Theorem 2.6.1 with and . This gives that Fourier PCA produces output such that there are signs and permutation such that
[TABLE]
with sample and time complexity. Choose so that the RHS of (2.28) is .
The number of samples and time needed for orthogonalization is
. Substituting the value of the previous bound becomes . The probability of error, coming from the applications of Theorem 2.4.1 and 2.6.1 is at most . ∎
Lemma 13**.**
Let be such that . Then there exists a unitary matrix such that .
Proof.
This is related to a special case of the so-called orthogonal Procrustes problem [46, Section 12.4.1], where one looks for a unitary matrix that minimizes . A formula for an optimal is , where is the singular value decomposition of , with singular values . Although we do not need the fact that this minimizes , it is good for our purpose:
[TABLE]
By our assumption
[TABLE]
This implies . The claim follows. ∎
2.7 Improving orthogonalization
As noted above, the technique in [10], while being provably efficient and correct, suffers from practical implementation issues. Here we discuss two alternatives: orthogonalization by centroid body scaling and orthogonalization by using the empirical covariance. The former, orthogonalization via centroid body scaling, uses the samples already present in the algorithm rather than relying on a random walk to draw samples which are approximately uniform in the algorithm’s approximation of the centroid body (as is done in [10]). This removes the dependence on random walks and the ellipsoid algorithm; instead, we use samples that are distributed according to the original heavy-tailed distribution but non-linearly scaled to lie inside the centroid body. We prove in Lemma 14 that the covariance of this subset of samples is enough to orthogonalize the mixing matrix . Secondly, we prove that one can, in fact, “forget” that the data is heavy tailed and orthogonalize by using the empirical covariance of the data, even though it diverges, and that this is enough to orthogonalize the mixing matrix . However, as observed in experimental results, in general this has a downside compared to orthogonalization via centroid body in that it could cause numerical instability during the “second” phase of ICA as the data obtained is less well-conditioned. This is illustrated directly in the table in Figure 2.5 containing the singular value and condition number of the mixing matrix in the approximately orthogonal ICA model.
2.7.1 Orthogonalization via centroid body scaling
In [10], another orthogonalization procedure, namely orthogonalization via the uniform distribution in the centroid body is theoretically proven to work. Their procedure does not suffer from the numerical instabilities and composes well with the second phase of ICA algorithms. An impractical aspect of that procedure is that it needs samples from the uniform distribution in the centroid body.
We described orthogonalization via centroid body in Section 1, except for the estimation of , the Minkowski functional of the centroid body. The complete procedure is stated in Subroutine 3.
We now explain how to estimate the Minkowski functional. The Minkowski functional was informally described in Section 1. The Minkowski functional of is formally defined by . Our estimation of is based on an explicit linear program (LP) (2.36) that gives the Minkowski functional of the centroid body of a finite sample of exactly and then arguing that a sample estimate is close to the actual value for . For clarity of exposition, we only analyze formally a special case of LP (2.36) that decides membership in the centroid body of a finite sample of (LP (2.35)) and approximate membership in . This analysis is in Section 2.7.3. Accuracy guarantees for the approximation of the Minkowski functional follow from this analysis.
Lemma 14**.**
*Let be a random vector drawn from an ICA model such that for all we have and is symmetrically distributed. Let where is the Minkoswki functional of .
Then is an orthogonalizer of .*
Proof.
We will be applying Lemma 9. Let denote the set of absolutely symmetric product distributions over such that for all . For , let be equal to the distribution obtained by scaling as described earlier, that is, distribution of , where , is the Minkoswki functional of .
For all , is symmetric and which implies that , that is, is absolutely symmetric. Let . Then is equal to the distribution of . For any invertible linear transformation and measurable set , we have .
Thus is linear equivariant. Let . Then there exist and such that . We get . Let . Thus, where is a diagonal matrix with elements which are non-zero because we assume . This implies that is positive definite and thus by Lemma 9, is an orthogonalizer of . ∎
2.7.2 Orthogonalization via covariance
Here we show the somewhat surprising fact that orthogonalization of heavy-tailed signals is sometimes possible by using the “standard” approach: inverting the empirical covariance matrix. The advantage here, is that it is computationally very simple, specifically that having heavy-tailed data incurs very little computational penalty on the process of orthogonalization alone. It’s standard to use covariance matrix for whitening when the second moments of all independent components exist [57]: Given samples from the ICA model , we compute the empirical covariance matrix which tends to the true covariance matrix as we take more samples and set . Then one can show that is a rotation matrix, and thus by pre-multiplying the data by we obtain an ICA model , where the mixing matrix is a rotation matrix, and this model is then amenable to various algorithms. In the heavy-tailed regime where the second moment does not exist for some of the components, there is no true covariance matrix and the empirical covariance diverges as we take more samples. However, for any fixed number of samples one can still compute the empirical covariance matrix. In previous work (e.g., [28]), the empirical covariance matrix was used for whitening in the heavy-tailed regime with good empirical performance; [28] also provided some theoretical analysis to explain this surprising performance. However, their work (both experimental and theoretical) was limited to some very special cases (e.g., only one of the components is heavy-tailed, or there are only two components both with stable distributions without finite second moment).
We will show that the above procedure (namely pre-multiplying the data by ) “works” under considerably more general conditions, namely if -moment exists for for each independent component . By “works” we mean that instead of whitening the data (that is is rotation matrix) it does something slightly weaker but still just as good for the purpose of applying ICA algorithms in the next phase. It orthogonalizes the data, that is now is close to a matrix whose columns are orthogonal. In other words, is close to a diagonal matrix (in a sense made precise in Theorem 2.7.1).
Let be a real-valued symmetric random variable such that for some and . The following lemma from [10] says that the empirical average of the absolute value of converges to the expectation of . The proof, which we omit, follows an argument similar to the proof of the Chebyshev’s inequality. Let be the empirical average obtained from independent samples , i.e., .
Theorem 2.7.1** (Orthogonalization via covariance matrix).**
Let be given by ICA model . Assume that there exist and such that for all we have
(a) ,
(b) (normalization) , and
*(c) . Let be i.i.d. samples according to . Let and . Then for any , for a diagonal matrix with diagonal entries satisfying for all with probability when . *
Proof idea. For we have (due to our symmetry assumption on ) and . We have , where . The off-diagonal entries of converge to [math]: We have . Now by our assumption that -moments exist, Lemma 6 is applicable and implies that empirical average tends to the true average as we increase the number of samples. The true average is [math] because of our assumption of symmetry (alternatively, we could just assume that the and hence have been centered). The diagonal entries of are bounded away from [math]: This is clear when the second moment is finite, and follows easily by hypothesis (c) when it is not. Finally, one shows that if in the diagonal entries highly dominate the off-diagonal entries, then the same is true of .
Proof.
We have , where . By assumption, for . Note that and so by Lemma 6, for ,
[TABLE]
when .
Now let . Then when for all , we have . The union bound then implies
[TABLE]
when .
Next, we aim to bound which can be done by writing
[TABLE]
where . Consider the random variable . We can calculate and use a Chernoff bound to see
[TABLE]
and when , we have . Then with probability at least , all entries of are at least . Using this, if then with probability at least .
Similarly, suppose that and choose and
[TABLE]
so that when , we have and with probability at least . Invoking (2.2), when , we have
[TABLE]
with probability at least .
Finally, we upper bound for a fixed by using Markov’s inequality:
[TABLE]
so that for all with probability at least when . Therefore, when , we have , , and for all with overall probability at least . ∎
We used Lemma 1.
In Theorem 2.7.1, the diagonal entries are lower bounded, which avoids some degeneracy, but they could still grow quite large because of the heavy tails. This is a real drawback of orthogonalization via covariance. HTICA, using the more sophisticated orthogonalization via centroid body scaling does not have this problem. We can see this in the right table of Figure 2.5, where the condition number of “centroid” is much smaller than the condition number of “covariance.”
2.7.3 New Membership oracle for the centroid body
We will now describe and theoretically justify a new and practically efficient -weak membership oracle for , which is a black-box that can answer approximate membership queries in . More precisely:
Definition 5**.**
The -weak membership problem for is the following: Given a point and a rational number , either (i) assert that , or (ii) assert that . An -weak membership oracle for is an oracle that solves the weak membership problem for . For , an -weak membership oracle for acts as follows: Given a point , with probability at least it solves the -weak membership problem for , and otherwise its output can be arbitrary.
We start with an informal description of the algorithm and its correctness.
The algorithm implementing the oracle (Subroutine 4) is the following: Let be a query point. Let be a sample of random vector . Given the sample, let be uniformly distributed in . Output YES if , else output NO.
Idea of the correctness of the algorithm: If is not in , then there is a hyperplane separating from . Let be the hyperplane, satisfying , and for every . Thus, we have and . We have
[TABLE]
By Lemma 6, is within of when is large enough with probability at least over the sample . In particular, , which implies and the algorithm outputs NO, with probability at least .
If is in , let . We will prove the following claim:
Informal claim (Lemma 17): For , for large enough and with probability at least there is so that .
This claim applied to to get , convexity of and the fact that contains (Lemma 15) imply that and the algorithm outputs YES.
We will prove the claim now. Let . By the dual characterization of the centroid body (Proposition 1), there exists a function such that with . Let We have and . By Lemma 6 and a union bound over every coordinate we get for large enough.
Formal Argument
Lemma 15**.**
Let be an absolutely symmetrically distributed random vector such that and for all . Let be a sample of i.i.d. copies of . Let be a random vector, uniformly distributed in . Then whenever
[TABLE]
Proof.
From Lemma 5 we know . It is enough to apply Lemma 17 to with and . This gives, for any , . In particular, . ∎
Proposition 1** (Dual characterization of centroid body).**
Let be a -dimensional random vector with finite first moment, that is, for all we have . Then
[TABLE]
Proof.
Let denote the rhs of the conclusion.We will show that is a non-empty, closed convex set and show that , which implies (2.34).
By definition, is a non-empty bounded convex set. To see that it is closed, let be a sequence in such that . Let be the function associated to according to the definition of . Let be the distribution of . We have and, passing to a subsequence , converges to in the weak- topology , where . 111This is a standard argument, see [MR2759829] for the background. Map is in . [MR2759829, Theorem 4.13] gives that is a separable Banach space. [MR2759829, Theorem 3.16] (Banach-Alaoglu-Bourbaki) gives that the unit ball in is compact in the weak-* topology. [MR2759829, Theorem 3.28] gives that the unit ball in is metrizable and therefore sequentially compact in the weak-* topology. Therefore, any bounded sequence in has a convergent subsequence in the weak-* topology. This implies . Thus, we have and is closed.
To conclude, we compute and see that it is the same as the definition of . In the following equations ranges over functions such that is Borel-measurable and .
[TABLE]
∎
Lemma 16** (LP).**
Let be a random vector uniformly distributed in . Let . Then:
. 2. 2.
Point iff there is a solution to the following linear feasibility problem:
[TABLE] 3. 3.
Let be the optimal value of (always feasible) linear program
[TABLE]
with if the linear program is unbounded. Then the Minkowski functional of at is .
Proof.
This is proven in [72]. It is also a special case of Proposition 1. We include an argument here for completeness. Let . We compute to see it is the same as in the definition of (Definition 4). As and are non-empty compact convex sets, this implies . We have
[TABLE] 2. 2.
This follows immediately from part 1. 3. 3.
This follows from part 1 and the definition of Minkowski functional.∎
Proposition 2** (Correctness of Subroutine 4).**
Let be given by an ICA model such that for all we have , is symmetrically distributed and normalized so that . Then, given a query point , , , , and , Subroutine 4 is an -weak membership oracle for and with probability using time and sample complexity
Proof.
Let be uniformly random in . There are two cases corresponding to the guarantees of the oracle:
- •
Case . Then there is a hyperplane separating from . Let be the separating hyperplane, parameterized so that , , , and for every . In this case and . At the same time, .
We want to apply Lemma 6 to to get that is within of . For this we need a bound on the -moment of . We use the bound given in 2.11.
- •
Case . Let . Let . Then . Invoke Lemma 17 for i.i.d. sample of with and equal to some to be fixed later to conclude . That is, there exist such that
[TABLE]
Let . Given (2.37) and the relationships and , we have
[TABLE]
From Lemma 15 with and equivariance of the centroid body (Lemma 3) we get . This and convexity of imply . In particular, the ball around of radius
[TABLE]
is contained in . The choice and (2.38) imply and Subroutine 4 outputs YES whenever
[TABLE]
To conclude, remember that . Therefore (from Lemma 5 and equivariance of the centroid body, Lemma 3). This implies
[TABLE]
The claim follows. ∎
Lemma 17**.**
Let be a -dimensional random vector such that for all coordinates we have . Let . Let be an i.i.d. sample of . Let be uniformly random in . Let , . If , then, with probability at least , .
Proof.
By Proposition 1, there exists a measurable function , such that . Let
[TABLE]
By Proposition 1, .
We have and, for every coordinate ,
[TABLE]
By Lemma 6 and for any fixed coordinate we have, over the choice of ,
[TABLE]
whenever . A union bound over choices of gives:
[TABLE]
So whenever
[TABLE]
and . The claim follows. ∎
2.8 Empirical Study
In this section, we show experimentally that heavy-tailed data poses a significant challenge for current ICA algorithms, and compare them with HTICA in different settings. We observe some clear situations where heavy-tails seriously affect the standard ICA algorithms, and that these problems are frequently avoided by using the heavy-tailed ICA framework. In some cases, HTICA does not help much, but maintains the same performance of plain FastICA.
To generate the synthetic data, we create a simple heavy-tailed density function proportional to , which is symmetric, and for , is the density of a distribution which has finite moment. The signal is generated with each independently distributed from . The mixing matrix is generated with each coordinate i.i.d. , columns normalized to unit length. To compare the quality of recovery, the columns of the estimated mixing matrix, are permuted to align with the closest matching column of , via the Hungarian algorithm. We use the Frobenius norm to measure the error, but all experiments were also performed using the well-known Amari index [6]; the results have similar behavior and are not presented here.
2.8.1 Heavy-tailed ICA when is orthogonal: Gaussian damping and experiments
Focusing on the third step above, where the mixing matrix already has orthogonal columns, ICA algorithms already suffer dramatically from the presence of heavy-tailed data. As proposed in [10], Gaussian damping is a preprocessing technique that converts data from an ICA model , where is unitary (columns are orthogonal with unit -norm) to data from a related ICA model , where is a parameter to be chosen. The independent components of have finite moments of all orders and so the existing algorithms can estimate .
Using samples of , we construct the damped random variable , with pdf . To normalize the right hand side, we can estimate
[TABLE]
so that
[TABLE]
If is a realization of , then is a realization of the random variable and we have that has pdf . To generate samples from this distribution, we use rejection sampling on samples from . When performing the damping, we binary search over so that about 25% of the samples are rejected. For more details about the technical requirements for choosing , see [10].
Figure 2.1 shows that, when is already a perfectly orthogonal matrix, but where may have heavy-tailed coordinates, several standard ICA algorithms perform better after damping the data. In fact, without damping, some do not appear to converge to a correct solution. We compare ICA with and without damping in this case: (1) FastICA using the fourth cumulant (“FastICA - pow3”), (2) FastICA using (“FastICA - tanh”), (3) JADE, and (4) Second Order Joint Diagonalization as in, e.g., [26] .
2.8.2 Experiments on synthetic heavy-tailed data
We now present the results of HTICA using different orthogonalization techniques: (1) Orthogonalization via covariance (Section 2.7.2 (2) Orthogonalization via the centroid body (Section 2.7.1) (3) the ground truth, directly inverting the mixing matrix (oracle), and (4) No orthogonalization, and also no damping (for comparison with plain FastICA) (identity).
The “mixed” regime in the left and middle of Figure 2.3 (where some signals are not heavy-tailed) demonstrates a very dramatic contrast between different orthogonalization methods, even when only two heavy-tailed signals are present.
In the experiment with different methods of orthogonalization it was observed that when all exponents are the same or very close, orthogonalization via covariance performs better than orthogonalization via centroid and the true mixing matrix as seen in Figure 2.3. A partial explanation is that, given the results in Figure 2.1, we know that equal exponents favor FastICA without damping and orthogonalization (identity in Figure 2.3). The line showing the performance with no orthogonalization and no damping (“identity”) behaves somewhat erratically, most likely due the presence of the heavy-tailed samples. Additionally, damping and the choice of parameter is sensitive to scaling. A scaled-up distribution will be somewhat hurt because fewer samples will survive damping.
2.8.3 ICA on speech data
While the above study on synthetic data provides interesting situations where heavy-tails can cause problems for ICA, we provide some results here which use real-world data, specifically human speech. To study the performance of HTICA on voice data, we first examine whether the data is heavy-tailed. The motivation to use speech data comes from observations by the signal processing community (e.g. [63]) that speech data can be modeled by -stable distributions. For an -stable distribution, with , only the moments of order less than will be finite. We present here some results on a data set of human speech according to the standard cocktail party model, from [40].
The physical setup of the experiments (the human speakers and microphones) is shown in Figure 2.4.
To estimate whether the data is heavy-tailed, as in [63], we estimate parameter of a best-fit -stable distribution. This estimate is in Figure 2.5 for one of the data sets collected. We can see that the estimated is clearly in the heavy-tailed regime for some signals.
Using data from [40], we perform the same experiment as in Section 2.8.2: generate a random mixing matrix with unit length columns, mix the data, and try to recover the mixing matrix. Although the mixing is synthetic, the setting makes the resulting mixed signals same as real. Specifically, the experiment was conducted in a room with chairs, carpet, plasterboard walls, and windows on one side. There was natural noise including vents, computers, florescent lights, and traffic noise through the windows.
Figure 2.5 demonstrates that HTICA (orthogonalizing with centroid body scaling, Section 2.7.1) applied to speech data yields some noticeable improvement in the recovery of the mixing matrix, primarily in that it is less susceptible to data that causes FastICA to have large error “spikes.” Moreover, in many cases, running only FastICA on the mixed data failed to even recover all of the speech signals, while HTICA succeeded. In these cases, we had to re-start FastICA until it recovered all the signals.
Chapter 3 Learning simplices
We are given uniformly random samples from an unknown convex body in , how many samples are needed to approximately reconstruct the body? It seems intuitively clear, at least for , that if we are given sufficiently many such samples then we can reconstruct (or learn) the body with very little error. For general , it is known to require samples [47] (see also [65] for a similar lower bound in a different but related model of learning). This is an information-theoretic lower bound and no computational considerations are involved. As mentioned in [47], it turns out that if the body has few facets (e.g. polynomial in ), then polynomial in samples are sufficient for approximate reconstruction. This is an information-theoretic upper bound and no efficient algorithms (i.e., with running time poly) are known. (We remark that to our knowledge the same situation holds for polytopes with poly vertices.) Here, we study the reconstruction problem for the special case when the input bodies are restricted to be (full-dimensional) simplices. We show that in this case one can in fact learn the body efficiently. More precisely, the algorithm knows that the input body is a simplex but only up to an affine transformation, and the problem is to recover this affine transformation. This answers a question of [43, Section 6].
The problem of learning a simplex is also closely related to the well-studied problem of learning intersections of half-spaces. Suppose that the intersection of half-spaces in is bounded, and we are given poly uniformly random samples from it. Then our learning simplices result directly implies that we can learn the half-spaces. This also has the advantage of being a proper learning algorithm, meaning that the output of the algorithm is a set of half-spaces, unlike many of the previous algorithms.
Previous work.
Perhaps the first approach to learning simplices that comes to mind is to find a minimum volume simplex containing the samples. This can be shown to be a good approximation to the original simplex. (Such minimum volume estimators have been studied in machine learning literature, see e.g. [94] for the problem of estimating the support of a probability distribution. We are not aware of any technique that applies to our situation and provides theoretical guarantees.) However, the problem of finding a minimum volume simplex is in general NP-hard [79]. This hardness is not directly applicable for our problem because our input is a random sample and not a general point set. Nevertheless, we do not have an algorithm for directly finding a minimum volume simplex; instead we use ideas similar to those used in Independent Component Analysis (ICA). ICA studies the following problem: Given a sample from an affine transformation of a random vector with independently distributed coordinates, recover the affine transformation (up to some unavoidable ambiguities). [43] gave an efficient algorithm for this problem (with some restrictions on the allowed distributions, but also with some weaker requirements than full independence) along with most of the details of a rigorous analysis (a complete analysis of a special case can be found in [13]; see also [106] for a generalization of ICA to subspaces along with a rigorous analysis). The problem of learning parallelepipeds from uniformly random samples is a special case of this problem. [43] asked if one could learn other convex bodies, and in particular simplices, efficiently from uniformly random samples. [76] gave a simpler and rigorous algorithm and analysis for the case of learning parallelepipeds with similarities to the popular FastICA algorithm of [56]. The algorithm in [76] is a first order algorithm unlike Frieze et al.’s second order algorithm.
The algorithms in both [43, 76] make use of the fourth moment function of the probability distribution. Briefly, the fourth moment in direction is , where is the random variable distributed according to the input distribution. The moment function can be estimated from the samples. The independent components of the distribution correspond to local maxima or minima of the moment function, and can be approximately found by finding the local maxima/minima of the moment function estimated from the sample.
More information on ICA including historical remarks can be found in [58, 34]. Ideas similar to ICA have been used in statistics in the context of projection pursuit since the mid-seventies. It is not clear how to apply ICA to the simplex learning problem directly as there is no clear independence among the components. Let us note that [43] allow certain kinds of dependencies among the components, however this does not appear to be useful for learning simplices.
Learning intersections of half-spaces is a well-studied problem in learning theory. The problem of PAC-learning intersections of even two half-spaces is open, and there is evidence that it is hard at least for sufficiently large number of half-spaces: E.g., [67] prove that learning intersections of half-spaces in (for constant ) is hard under standard cryptographic assumptions (PAC-learning is possible, however, if one also has access to a membership oracle in addition to random samples [68]). Because of this, much effort has been expended on learning when the distribution of random samples is some simple distribution, see e.g. [66, 104, 105] and references therein. This line of work makes substantial progress towards the goal of learning intersections of half-spaces efficiently, however it falls short of being able to do this in time polynomial in and ; in particular, these algorithms do not seem to be able to learn simplices. The distribution of samples in these works is either the Gaussian distribution or the uniform distribution over a ball. [43] and [47] consider the uniform distribution over the intersection. Note that this requires that the intersection be bounded. Note also that one only gets positive samples in this case unlike other work on learning intersections of half-spaces. The problem of learning convex bodies can also be thought of as learning a distribution or density estimation problem for a special class of distributions.
[49] show how to reconstruct a polytope with vertices in , given its first moments in random directions. In our setting, where we have access to only a polynomial number of random samples, it’s not clear how to compute moments of such high orders to the accuracy required for the algorithm of [49] even for simplices.
A recent and parallel work of [7] is closely related to ours. They show that tensor decomposition methods can be applied to low-order moments of various latent variable models to estimate their parameters. The latent variable models considered by them include Gaussian mixture models, hidden Markov models and latent Dirichlet allocations. The tensor methods used by them and the local optima technique we use seem closely related. One could view our work, as well as theirs, as showing that the method of moments along with existing algorithmic techniques can be applied to certain unsupervised learning problems.
In Section 3.6, we give the theorems and proofs which allow the efficient reduction to ICA.
Our results
For clarity of the presentation, we use the following machine model for the running time: a random access machine that allows the following exact arithmetic operations over real numbers in constant time: addition, subtraction, multiplication, division and square root.
The estimation error is measured using total variation distance, denoted (see Section 3.1).
Theorem 3.0.1**.**
There is an algorithm (Algorithm 2 below) such that given access to random samples from a simplex , with probability at least over the sample and the randomness of the algorithm, it outputs vectors that are the vertices of a simplex so that . The algorithm runs in time polynomial in , and .
As mentioned earlier, our algorithm uses ideas from ICA. Our algorithm uses the third moment instead of the fourth moment used in certain versions of ICA. The third moment is not useful for learning symmetric bodies such as the cube as it is identically [math]. It is however useful for learning a simplex where it provides useful information, and is easier to handle than the fourth moment. One of the main contributions of our work is the understanding of the third moment of a simplex and the structure of local maxima. This is more involved than in previous work as the simplex has no obvious independence structure, and the moment polynomial one gets has no obvious structure unlike for ICA.
The probability of success of the algorithm can be “boosted” so that the dependence of the running time on is only linear in as follows: The following discussion uses the space of simplices with total variation distance as the underlying metric space. Let be the target distance. Take an algorithm that succeeds with probability and error parameter to be fixed later (such as Algorithm 2 with ). Run the algorithm times to get simplices. By a Chernoff-type argument, at least simplices are within of the input simplex with probability at least .
By sampling, we can estimate the distances between all pairs of simplices with additive error less than in time polynomial in and so that all estimates are correct with probability at least . For every output simplex, compute the number of output simplices within estimated distance . With probability at least both of the desirable events happen, and then necessarily there is at least one output simplex, call it , that has output simplices within estimated distance . Any such must be within of the input simplex. Thus, set .
While our algorithm for learning simplices uses techniques for ICA, we have to do substantial work to make those techniques work for the simplex problem. We also show a more direct connection between the problem of learning a simplex and ICA: a randomized reduction from the problem of learning a simplex to ICA. The connection is based on a known representation of the uniform measure on a simplex as a normalization of a vector having independent coordinates. Similar representations are known for the uniform measure in an -dimensional ball (denoted ) [16] and the cone measure on the boundary of an ball [92, 83, 98] (see Section 3.1 for the definition of the cone measure). These representations lead to a reduction from the problem of learning an affine transformation of an ball to ICA. These reductions show connections between estimation problems with no obvious independence structure and ICA. They also make possible the use of any off-the-shelf implementation of ICA. However, the results here do not supersede our result for learning simplices because to our knowledge no rigorous analysis is available for the ICA problem when the distributions are the ones in the above reductions.
Idea of the algorithm.
The new idea for the algorithm is that after putting the samples in a suitable position (see below), the third moment of the sample can be used to recover the simplex using a simple FastICA-like algorithm. We outline our algorithm next.
As any full-dimensional simplex can be mapped to any other full-dimensional simplex by an invertible affine transformation, it is enough to determine the translation and linear transformation that would take the given simplex to some canonical simplex. As is well-known for ICA-like problems (see, e.g., [43]), this transformation can be determined up to a rotation from the mean and the covariance matrix of the uniform distribution on the given simplex. The mean and the covariance matrix can be estimated efficiently from a sample. A convenient choice of an -dimensional simplex is the convex hull of the canonical vectors in . We denote this simplex and call it the standard simplex. So, the algorithm begins by picking an arbitrary invertible affine transformation that maps onto the hyperplane . We use a so that is an isotropic111See Section 3.1. simplex. In this case, the algorithm brings the sample set into isotropic position and embeds it in using . After applying these transformations we may assume (at the cost of small errors in the final result) that our sample set is obtained by sampling from an unknown rotation of the standard simplex that leaves the all-ones vector (denoted from now on) invariant (thus this rotation keeps the center of mass of the standard simplex fixed), and the problem is to recover this rotation.
To find the rotation, the algorithm will find the vertices of the rotated simplex approximately. This can be done efficiently because of the following characterization of the vertices: Project the vertices of the simplex onto the hyperplane through the origin orthogonal to and normalize the resulting vectors. Let denote this set of points. Consider the problem of maximizing the third moment of the uniform distribution in the simplex along unit vectors orthogonal to . Then is the complete set of local maxima and the complete set of global maxima (Theorem 3.5.1). A fixed point-like iteration (inspired by the analysis of FastICA [56] and of gradient descent in [76]) starting from a random point in the unit sphere finds a local maximum efficiently with high probability. By the analysis of the coupon collector’s problem, repetitions are highly likely to find all local maxima.
Idea of the analysis.
In the analysis, we first argue that after putting the sample in isotropic position and mapping it through , it is enough to analyze the algorithm in the case where the sample comes from a simplex that is close to a simplex that is the result of applying a rotation leaving invariant to the standard simplex. The closeness here depends on the accuracy of the sample covariance and mean as an estimate of the input simplex’s covariance matrix and mean. A sample of size guarantees ([3, Theorem 4.1], [100, Corollary 1.2]) that the covariance and mean are close enough so that the uniform distributions on and are close in total variation. We show that the subroutine that finds the vertices (Subroutine 5), succeeds with some probability when given a sample from . By definition of total variation distance, Subroutine 5 succeeds with almost as large probability when given a sample from (an argument already used in [76]). As an additional simplifying assumption, it is enough to analyze the algorithm (Algorithm 2) in the case where the input is isotropic, as the output distribution of the algorithm is equivariant with respect to affine invertible transformations as a function of the input distribution.
3.1 Preliminaries
An -simplex is the convex hull of points in that do not lie on an -dimensional affine hyperplane. It will be convenient to work with the standard -simplex living in defined as the convex hull of the canonical unit vectors ; that is
[TABLE]
The canonical simplex living in is given by
[TABLE]
Note that is the facet of opposite to the origin.
Let denote the -dimensional Euclidean ball.
The complete homogeneous symmetric polynomial of degree in variables , denoted , is the sum of all monomials of degree in the variables:
[TABLE]
Also define the -th power sum as
[TABLE]
For a vector , we define
[TABLE]
Vector denotes the all ones vector (the dimension of the vector will be clear from the context).
A random vector is isotropic if and . A compact set in is isotropic if a uniformly distributed random vector in it is isotropic. The inradius of an isotropic -simplex is , the circumradius is .
The total variation distance between two probability measures is for measurable . For two compact sets , we define the total variation distance as the total variation distance between the corresponding uniform distributions on each set. It can be expressed as
[TABLE]
This identity implies the following elementary estimate:
Lemma 18**.**
Let be two compact sets in . Let such that . Then .
Proof.
We have . Triangle inequality implies the desired inequality. ∎
Lemma 19**.**
Consider the coupon collector’s problem with coupons where every coupon occurs with probability at least . Let . Then with probability at least all coupons are collected after trials.
Proof.
The probability that a particular coupon is not collected after that many trials is at most
[TABLE]
The union bound over all coupons implies the claim. ∎
For a point , is the standard norm. The unit ball is defined by
[TABLE]
The Gamma distribution is denoted as and has density , with shape parameters . is the exponential distribution, denoted . The Gamma distribution also satisfies the following additivity property: If is distributed as and is distributed as , then is distributed as .
The cone measure on the surface of centrally symmetric convex body in [16, 92, 83, 98] is defined by
[TABLE]
It is easy to see that is uniform on for .
From [92] and [83] we have the following representation of the cone measure on :
Theorem 3.1.1**.**
Let be iid random variables with density proportional to . Then the random vector is independent of . Moreover, is distributed according to .
From [16], we also have the following variation, a representation of the uniform distribution in :
Theorem 3.1.2**.**
Let be iid random variables with density proportional to . Let be a random variable distributed as , independent of . Then the random vector
[TABLE]
is uniformly distributed in .
See e.g. [20, Section 20] for the change of variable formula in probability.
3.2 Computing the moments of a simplex
The -th moment over is the function
[TABLE]
In this section we present a formula for the moment over . Similar more general formulas appear in [69]. We will use the following result from [51] for :
[TABLE]
From the above we can easily derive a formula for integration over :
[TABLE]
Now
[TABLE]
The variant of Newton’s identities for the complete homogeneous symmetric polynomial gives the following relations which can also be verified easily by direct computation:
[TABLE]
[TABLE]
Divide the above integral by the volume of the standard simplex to get the moment:
[TABLE]
3.3 Subroutine for finding the vertices of a rotated standard simplex
In this section we solve the following simpler problem: Suppose we have samples from a rotated copy of the standard simplex, where the rotation is such that it leaves invariant. The problem is to approximately estimate the vertices of the rotated simplex from the samples.
We will analyze our algorithm in the coordinate system in which the input simplex is the standard simplex. This is only for convenience in the analysis and the algorithm itself does not know this coordinate system.
As we noted in the introduction, our algorithm is inspired by the algorithm of [76] for the related problem of learning hypercubes and also by the FastICA algorithm in [56]. New ideas are needed for our algorithm for learning simplices; in particular, our update rule is different. With the right update rule in hand the analysis turns out to be quite similar to the one in [76].
We want to find local maxima of the sample third moment. A natural approach to do this would be to use gradient descent or Newton’s method (this was done in [43]). Our algorithm, which only uses first order information, can be thought of as a fixed point algorithm leading to a particularly simple analysis and fast convergence. Before stating our algorithm we describe the update rule we use.
We will use the abbreviation . Then, from the expression for we get
[TABLE]
Solving for we get
[TABLE]
While the above expressions are in the coordinate system where the input simplex is the canonical simplex, the important point is that all terms in the last expression can be computed in any coordinate system that is obtained by a rotation leaving invariant. Thus, we can compute as well independently of what coordinate system we are working in. This immediately gives us the algorithm below. We denote by the sample third moment, i.e., for samples. This is a polynomial in , and the gradient is computed in the obvious way. Moreover, the gradient of the sample moment is clearly an unbiased estimator of the gradient of the moment; a bound on the deviation is given in the analysis (Lemma 20). For each evaluation of the gradient of the sample moment, we use a fresh sample.
It may seem a bit alarming that the fixed point-like iteration is squaring the coordinates of , leading to an extremely fast growth (see Equation 3.1 and Subroutine 1). But, as in other algorithms having quadratic convergence like certain versions of Newton’s method, the convergence is very fast and the number of iterations is small. We show below that it is , leading to a growth of that is polynomial in and . The boosting argument described in the introduction makes the final overall dependence in to be only linear in .
We state the following subroutine for instead of (thus it is learning a rotated copy of instead of ). This is for notational convenience so that we work with instead of .
Lemma 20**.**
Let be a constant, , and . Suppose that Subroutine 1 uses a sample of size for each evaluation of the gradient and runs for iterations. Then with probability at least Subroutine 1 outputs a vector within distance from a vertex of the input simplex. With respect of the process of picking a sample and running the algorithm, each vertex is equally likely to be the nearest.
Note that if we condition on the sample, different vertices are not equally likely over the randomness of the algorithm. That is, if we try to find all vertices running the algorithm multiple times on a fixed sample, different vertices will be found with different likelihoods.
Proof.
Our analysis has the same outline as that of [76]. This is because the iteration that we get is the same as that of [76] except that cubing is replaced by squaring (see below); however some details in our proof are different. In the proof below, several of the inequalities are quite loose and are so chosen to make the computations simpler.
We first prove the lemma assuming that the gradient computations are exact and then show how to handle samples. We will carry out the analysis in the coordinate system where the given simplex is the standard simplex. This is only for the purpose of the analysis, and this coordinate system is not known to the algorithm. Clearly, . It follows that,
[TABLE]
Now, since we choose randomly, with probability at least one of the coordinates of is greater than all the other coordinates in absolute value by a factor of at least , where . (A similar argument is made in [76] with different parameters. We briefly indicate the proof for our case: The probability that the event in question does not happen is less than the probability that there are two coordinates and such that their absolute values are within factor , i.e. . The probability that for given this event happens can be seen as the Gaussian area of the four sectors (corresponding to the four choices of signs of ) in the plane each with angle less than . By symmetry, the Gaussian volume of these sectors is . The probability that such a pair exists is less than .) Assuming this happens, then after iterations, the ratio between the largest coordinate (in absolute value) and the absolute value of any other coordinate is at least . Thus, one of the coordinates is very close to and others are very close to [math], and so is very close to a vertex of the input simplex.
Now we drop the assumption that the gradient is known exactly. For each evaluation of the gradient we use a fresh subset of samples of points. Here is chosen so that each evaluation of the gradient is within -distance from its true value with probability at least , where will be set at the end of the proof. An application of the Chernoff bound yields that we can take ; we omit the details. Thus all the evaluations of the gradient are within distance from their true values with probability at least .
We assumed that our starting vector has a coordinate greater than every other coordinate by a factor of in absolute value; let us assume without loss of generality that this is the first coordinate. Hence . When expressing in terms of the gradient, the gradient gets multiplied by (we are assuming ), keeping this in mind and letting we get for
[TABLE]
If , where will be determined later, then we get
[TABLE]
Else,
[TABLE]
where we used and in the last inequality.
We choose so that
[TABLE]
For this, or equivalently suffices.
For satisfying (3.3) we have . It then follows from (3.3) that the first coordinate continues to remain the largest in absolute value by a factor of at least after each iteration. Also, once we have , we have for all .
(3.3) gives that after iterations we have
[TABLE]
Now if is such that , we will be guaranteed that . This condition is satisfied if we have , or equivalently . Now using (3.3) it suffices to choose so that . Thus we can take .
Hence we get . It follows that for , the -distance from the vertex is at most for ; we omit easy details.
Now we set our parameters: and and satisfies all the constraints we imposed on . Choosing , we get that the procedure succeeds with probability at least . Now setting gives the overall probability of error , and the number of samples and iterations as claimed in the lemma. ∎
3.4 Learning simplices
In this section we give our algorithm for learning general simplices, which uses the subroutine from the previous section. The learning algorithm uses an affine map that maps some isotropic simplex to the standard simplex. We describe now a way of constructing such a map: Let be a matrix having as columns an orthonormal basis of in . To compute one such , one can start with the -by- matrix that has ones in the diagonal and first column, everything else is zero. Let be a QR-decomposition of . By definition we have that the first column of is parallel to and the rest of the columns span . Given this, let be the matrix formed by all columns of except the first. We have that the set is the set of vertices of a regular -simplex. Each vertex is at distance
[TABLE]
from the origin, while an isotropic simplex has vertices at distance from the origin. So an affine transformation that maps an isotropic simplex in to the standard simplex in is .
To simplify the analysis, we pick a new sample to find every vertex, as this makes every vertex equally likely to be found when given a sample from an isotropic simplex. (The core of the analysis is done for an isotropic simplex; this is enough as the algorithm’s first step is to find an affine transformation that puts the input simplex in approximately isotropic position. The fact that this approximation is close in total variation distance implies that it is enough to analyze the algorithm for the case of exact isotropic position, the analysis carries over to the approximate case with a small loss in the probability of success. See the proof below for the details.) A practical implementation may prefer to select one such sample outside of the for loop, and find all the vertices with just that sample—an analysis of this version would involve bounding the probability that each vertex is found (given the sample, over the choice of the starting point of gradient descent) and a variation of the coupon collector’s problem with coupons that are not equally likely.
Theorem 3.0.1.
As a function of the input simplex, the distribution of the output of the algorithm is equivariant under invertible affine transformations. Namely, if we apply an affine transformation to the input simplex, the distribution of the output is equally transformed.222To see this: the equivariance of the algorithm as a map between distributions is implied by the equivariance of the algorithm on any given input sample. Now, given the input sample, if we apply an affine transformation to it, this transformation is undone except possibly for a rotation by the step . A rotation may remain because of the ambiguity in the characterization of . But the steps of the algorithm that follow the definition of are equivariant under rotation, and the ambiguous rotation will be removed at the end when is applied again in the last step. The notion of error, total variation distance, is also invariant under invertible affine transformations. Therefore, it is enough to analyze the algorithm when the input simplex is in isotropic position. In this case (see Section 2.2) and we can set so that with probability at least (by an easy application of Chernoff’s bound), for some to be fixed later. Similarly, using results from [3, Theorem 4.1], a choice of implies that the empirical second moment matrix
[TABLE]
satisfies with probability at least . We have and this implies . Now, is an iid sample from a simplex . Simplex is close in total variation distance to some isotropic simplex333The isotropic simplex will typically be far from the (isotropic) input simplex, because of the ambiguity up to orthogonal transformations in the characterization of . . More precisely, Lemma 21 below shows that
[TABLE]
with probability at least .
Assume for a moment that are from . The analysis of Subroutine 5 (fixed point-like iteration) given in Lemma 20 would guarantee the following: Successive invocations to Subroutine 5 find approximations to vertices of within Euclidean distance for some to be determined later and . We ask for each invocation to succeed with probability at least with . Note that each vertex is equally likely to be found. The choice of is so that, if all invocations succeed (which happens with probability at least ), then the analysis of the coupon collector’s problem, Lemma 19, implies that we fail to find a vertex with probability at most . Overall, we find all vertices with probability at least .
But in reality samples are from , which is only close to . The estimate from (3.4) with appropriate gives
[TABLE]
which implies that the total variation distance between the joint distribution of all samples used in the loop and the joint distribution of actual samples from the isotropic simplex is at most , and this implies that the loop finds approximations to all vertices of when given samples from with probability at least . The points in are still within Euclidean distance of corresponding vertices of .
To conclude, we turn our estimate of distances between estimated and true vertices into a total variation estimate, and map it back to the input simplex. Let . As maps an isotropic simplex to a standard simplex, we have that is an isometry, and therefore the vertices of are within distance of the corresponding vertices of . Thus, the corresponding support functions are uniformly within
[TABLE]
of each other on the unit sphere. This and the fact that imply
[TABLE]
Thus, by Lemma 18, and this implies that the total variation distance between the uniform distributions on and the input simplex is at most . Over all random choices, this happens with probability at least . We set . ∎
Lemma 21**.**
Let be an -dimensional isotropic simplex. Let be an -by- positive definite matrix such that . Let be an -dimensional vector such that . Let be an -by- matrix such that . Let be the simplex . Then there exists an isotropic simplex such that .
Proof.
We use an argument along the lines of the orthogonal Procrustes problem (nearest orthogonal matrix to , already in [76, Proof of Theorem 4]): Let be the singular value decomposition of . Let be an orthogonal matrix (that approximates ). Let .
We have . Let , be the minimum and maximum singular values of , respectively. This implies:
[TABLE]
As , , is orthogonal and is convex, we have
[TABLE]
Also,
[TABLE]
This in (3.5) gives
[TABLE]
This and Lemma 18 imply
[TABLE]
The estimate on gives , . Thus
[TABLE]
∎
3.5 The local and global maxima of the 3rd moment of the standard simplex and the isotropic simplex
In this section we study the structure of the set of local maxima of the third moment as a function of the direction (which happens to be essentially as discussed in Section 3.2). This is not necessary for our algorithmic result, however it gives insight into the geometry of the third moment (the location of local maxima/minima and stationary points) and suggests that more direct optimization algorithms like gradient descent and Newton’s method will also work, although we will not prove that.
Theorem 3.5.1**.**
Let be an isotropic simplex. Let be random in . Let be the set of normalized vertices of . Then is a complete set of local maxima and a complete set of global maxima of given by .
Proof idea: Embed the simplex in . Show that the third moment is proportional to the complete homogeneous symmetric polynomial of degree 3, which for the relevant directions is proportional to the sum of cubes. To conclude, use first and second order optimality conditions to characterize the set of local maxima.
Proof.
Consider the standard simplex
[TABLE]
and identify it with via a linear map so that . Let be random in . Consider given by . Let be the equivalent feasible set for the embedded problem. We have for any and some constant independent of . To get the theorem, it is enough to show that the local maxima of in are precisely the normalized versions of the projections of the canonical vectors onto the hyperplane orthogonal to . According to Section 3.2, for we have
[TABLE]
Using a more convenient but equivalent constant, we want to enumerate the local maxima of the problem
[TABLE]
The Lagrangian function is
[TABLE]
The first order condition is , that is,
[TABLE]
Consider this system of equations on for any fixed , . Let , . The first order condition says , where is convex and is affine. That is, the s can take at most two different values. As our optimization problem (3.6) is symmetric under permutation of the coordinates, we conclude that, after putting the coordinates of a point in non-increasing order, if is a local maximum of (3.6), then must be of the form
[TABLE]
where and there are exactly s and s, for .
We will now study the second order necessary condition (SONC) to eliminate from the list of candidates all vectors with . It is easy to see that the surviving vectors are exactly the promised scaled projections of the canonical vectors. This vectors must all be local and global maxima: At least one of them must be a global maximum as we are maximizing a continuous function over a compact set and all of them have the same objective value so all of them are local and global maxima.
The SONC at asks for the Hessian of the Lagrangian to be negative semidefinite when restricted to the tangent space to the constraint set at [71, Section 11.5]. We compute the Hessian (recall that is the vector of the squared coordinates of ):
[TABLE]
[TABLE]
where is the -by- matrix having the entries of in the diagonal and 0 elsewhere.
A vector in the tangent space is any such that , , and definiteness of the Hessian is determined by the sign of for any such , where
[TABLE]
Suppose is a critical point with . To see that such a cannot be a local maximum, it is enough to show , as in that case we can take to make the second derivative of positive in the direction .
In terms of , the constraints of (3.6) are , , and this implies , . The inner product between the first order condition (3.7) and implies . It is convenient to consider the change of variable , as now candidate critical points are parameterized by certain discrete values of in . This gives , and
[TABLE]
This implies
[TABLE]
In , the function given by is convex and symmetric around , where it attains its global minimum value, , which is positive. ∎
3.6 Probabilistic results used
In this section we show the probabilistic results underlying the reductions from learning simplices and balls to ICA. The results are Theorems 3.6.1 and 3.6.2. They each show a simple non-linear rescaling of the respective uniform distributions that gives a distribution with independent components (Definition 6).
Theorem 3.6.1 below gives us, in a sense, a “reversal” of the representation of the cone measure on , seen in Theorem 3.1.1. Given any random point in the standard simplex, we can apply a simple non-linear scaling and recover a distribution with independent components.
Definition 6**.**
We say that a random vector has independent components if it is an affine transformation of a random vector having independent coordinates.
Theorem 3.6.1**.**
Let be a uniformly random vector in the -dimensional standard simplex . Let be a random scalar distributed as . Then the coordinates of are iid as .
Moreover, if is an invertible linear transformation, then the random vector has independent components.
Proof.
In the case where , Theorem 3.1.1 restricted to the positive orthant implies that for random vector , if each is an iid exponential random variable , then has the same (joint) distribution as . Given the measurable function , has the same distribution as . That is, and have the same distribution444See [52, Theorem 1.1] for a similar argument in this context..
For the second part, we know by linearity. By the previous argument the coordinates of are independent. This implies that has independent components. ∎
The next lemma complements the main result in [16], Theorem 1 (Theorem 3.1.2 elsewhere here). They show a representation of the uniform distribution in , but they do not state the independence that we need for our reduction to ICA.
Lemma 22**.**
Let . Let be iid random variables with density proportional to . Let be a nonnegative random variable with distribution and independent of . Then the random vector
[TABLE]
is independent of .
Proof idea. We aim to compute the joint density, showing that it is a product of individual densities. To avoid complication, we raise everything to the th power, which eliminates extensive use of the chain rule involved in the change of variables.
Proof.
It is enough to show the claim conditioning on the orthant in which falls, and by symmetry it is enough to prove it for the positive orthant. Let random variable . Since raising (strictly) positive numbers to the th power is injective, it suffices to show that the random vector
[TABLE]
is independent of the random vector .
First, let be the interior of the support of , that is and consider and where
[TABLE]
and
[TABLE]
The random vector has a density supported on and
[TABLE]
is one-to-one from onto . Let be the determinant of its Jacobian. This Jacobian is
[TABLE]
which, by adding each of the first rows to the last row, reduces to
[TABLE]
the determinant of which is trivially .
We have that is nonzero in . Thus, has density supported on given by
[TABLE]
It is easy to see555See for example [16, proof of Theorem 3]. that each has density and thus has density by the additivity of the Gamma distribution. We then compute the joint density
[TABLE]
Since is independent of ,
[TABLE]
where
[TABLE]
The result follows. ∎
With this in mind, we show now our analog of Theorem 3.6.1 for .
Theorem 3.6.2**.**
Let be a uniformly random vector in . Let be a random scalar distributed as . Then the coordinates of are iid, each with density proportional to . Moreover, if is an invertible linear transformation, then the random vector given by has independent components.
Proof.
Let where each is iid as . Also, let be an independent random variable distributed as . Let S=\bigl{(}\sum_{i=1}^{n}|G_{i}|^{p}+W\bigr{)}^{1/p}.
By Lemma 22 and Theorem 3.1.2 we know has the same joint distribution as . Then for the measurable function , we immediately have has the same distribution as and thus has the same distribution as .
For the second part, since is a scalar, we have . By the previous argument we have that the coordinates of are independent. Thus, has independent components. ∎
This result shows that one can obtain a vector with independent components from a sample in a linearly transformed ball. In Section 3.7 we show that they are related in such as way that one can recover the linear transformation from the independent components via ICA.
3.7 Reducing simplex learning to ICA
We now show randomized reductions from the following two natural statistical estimation problems to ICA:
Problem 1** (simplex).**
Given uniformly random points from an -dimensional simplex, estimate the simplex.
This is the same problem of learning a simplex as in the rest of the section, we just restate it here for clarity.
To simplify the presentation for the second problem, we ignore the estimation of the mean of an affinely transformed distribution. That is, we assume that the ball to be learned has only been linearly transformed.
Problem 2** (linearly transformed balls).**
Given uniformly random points from a linear transformation of the -ball, estimate the linear transformation.
These problems do not have an obvious independence structure. Nevertheless, known representations of the uniform measure in an ball and the cone measure (defined in Section 3.1) on the surface of an ball can be slightly extended to map a sample from those distributions into a sample with independent components by a non-linear scaling step. The use of a non-linear scaling step to turn a distribution into one having independent components has been done before [96, 97], but there it is applied after finding a transformation that makes the distribution axis-aligned. This alignment is attempted using ICA (or variations of PCA) on the original distribution [96, 97], without independent components, and therefore the use of ICA is somewhat heuristic. One of the contributions of our reduction is that the rescaling we apply is “blind”, namely, it can be applied to the original distribution. In fact, the distribution does not even need to be isotropic (“whitened”). The distribution resulting from the reduction has independent components and therefore the use of ICA on it is well justified.
The reductions are given in Algorithms 3 and 4. To state the reductions, we denote by the invocation of an ICA routine. It takes samples of a random vector , where the coordinates of are independent, and returns an approximation to a square matrix such that is isotropic and has independent coordinates. The theory of ICA [32, Theorem 11] implies that if is isotropic and at most one coordinate is distributed as a Gaussian, then such an exists and it satisfies , where is a permutation matrix and is a diagonal matrix with diagonal entries in . We thus need the following definition to state our reduction: Let . That is, the uniform distribution in is isotropic.
As we do not state a full analysis of any particular ICA routine, we do not state explicit approximation guarantees.
Algorithm 3 works as follows: Let be an -dimensional random vector with iid coordinates distributed as . Let be the matrix having columns for . Let be random according to the distribution that results from scaling in the algorithm. Theorem 3.6.1 implies that and have the same distribution. Also, is isotropic and and have the same distribution. Thus, the discussion about ICA earlier in this section gives that the only separating matrices are such that where is a permutation matrix and is a diagonal matrix with diagonal entries in . That is, . As the last row of is all ones, the sign change step in Algorithm 3 undoes the effect of and recovers the correct orientation.
Similarly, Algorithm 4 works as follows: Let be a random vector with iid coordinates, each with density proportional to . Let be random according to the distribution that results from scaling in the algorithm. Theorem 3.6.2 implies that and have the same distribution. Also, is isotropic and we have and have the same distribution. Thus, the discussion about ICA earlier in this section gives that the only separating matrices are such that where is a permutation matrix and is a diagonal matrix with diagonal entries in . That is, . The fact that is symmetric with respect to coordinate permutations and sign changes implies that and is the same as . When , the assumptions in the discussion above about ICA are satisfied and Algorithm 4 is correct. When , the distribution of the scaled sample is Gaussian and this introduces ambiguity with respect to rotations in the definition of , but this ambiguity is no problem as it is counteracted by the fact that the ball is symmetric with respect to rotations.
Chapter 4 Learning Gaussian Mixture Models
The question of recovering a probability distribution from a finite set of samples is one of the most fundamental questions of statistical inference. While classically such problems have been considered in low dimension, more recently inference in high dimension has drawn significant attention in statistics and computer science literature.
In particular, an active line of investigation in theoretical computer science has dealt with the question of learning a Gaussian Mixture Model in high dimension. This line of work was started in [37] where the first algorithm to recover parameters using a number of samples polynomial in the dimension was presented. The method relied on random projections to a low dimensional space and required certain separation conditions for the means of the Gaussians. Significant work was done in order to weaken the separation conditions and to generalize the result (see e.g., [38, 12, 103, 1, 42]). Much of this work has polynomial sample and time complexity but requires strong separation conditions on the Gaussian components. A completion of the attempts to weaken the separation conditions was achieved in [18] and [74], where it was shown that arbitrarily small separation was sufficient for learning a general mixture with a fixed number of components in polynomial time. Moreover, a one-dimensional example given in [74] showed that an exponential dependence on the number of components was unavoidable unless strong separation requirements were imposed. Thus the question of polynomial learnability appeared to be settled. It is worth noting that while quite different in many aspects, all of these papers used a general scheme similar to that in the original work [38] by reducing high-dimensional inference to a small number of low-dimensional problems through appropriate projections.
However, a surprising result was recently proved in [54]. The authors showed that a mixture of Gaussians in dimension could be learned using a polynomial number of samples, assuming a non-degeneracy condition on the configuration of the means. The result in [54] is inherently high-dimensional as that condition is never satisfied when the means belong to a lower-dimensional space. Thus the problem of learning a mixture gets progressively computationally easier as the dimension increases, a “blessing of dimensionality!” It is important to note that this was quite different from much of the previous work, which had primarily used projections to lower-dimension spaces.
Still, there remained a large gap between the worst case impossibility of efficiently learning more than a fixed number of Gaussians in low dimension and the situation when the number of components is equal to the dimension. Moreover, it was not completely clear whether the underlying problem was genuinely easier in high dimension or our algorithms in low dimension were suboptimal. The one-dimensional example in [74] cannot answer this question as it is a specific worst-case scenario, which can be potentially ruled out by some genericity condition.
We take a step to eliminate this gap by showing that even very large mixtures of Gaussians can be polynomially learned. More precisely, we show that a mixture of Gaussians with equal known covariance can be polynomially learned as long as is bounded from above by a polynomial of the dimension and a certain more complex non-degeneracy condition for the means is satisfied. We show that if is high enough, these non-degeneracy conditions are generic in the smoothed complexity sense. Thus for any fixed , generic Gaussians can be polynomially learned in dimension .
Further, we prove that no such condition can exist in low dimension. A measure of non-degeneracy must be monotone in the sense that adding Gaussian components must make the condition number worse. However, we show that for points uniformly sampled from there are (with high probability) two mixtures of unit Gaussians with means on non-intersecting subsets of these points, whose distance is and which are thus not polynomially identifiable. More generally, in dimension the distance becomes . That is, the conditioning improves as the dimension increases, which is consistent with our algorithmic results.
To summarize, our contributions are as follows:
(1) We show that for any , a mixture of Gaussians in dimension can be learned in time and number of samples polynomial in and a certain “condition number” . For sufficiently high dimension, this results in an algorithm polynomial from the smoothed analysis point of view (Theorem 4.0.1). To do that we provide smoothed analysis of the condition number using certain results from [90] and anti-concentration inequalities. The main technical ingredient of the algorithm is a new “Poissonization” technique to reduce Gaussian mixture estimation to a problem of recovering a linear map of a product distribution known as underdetermined Independent Component Analysis (ICA). We combine this with the recent work on efficient algorithms for underdetermined ICA from [48] to obtain the necessary bounds.
(2) We show that in low dimension polynomial identifiability fails in a certain generic sense (see Theorem 4.0.3). Thus the efficiency of our main algorithm is truly a consequence of the “blessing of dimensionality” and no comparable algorithm exists in low dimension. The analysis is based on results from approximation theory and Reproducing Kernel Hilbert Spaces.
Moreover, we combine the approximation theory results with the Poissonization-based technique to show how to embed difficult instances of low-dimensional Gaussian mixtures into the ICA setting, thus establishing exponential information-theoretic lower bounds for underdetermined ICA in low dimension. To the best of our knowledge, this is the first such result in the literature.
We discuss our main contributions more formally now. The notion of Khatri–Rao power of a matrix is defined in Section 4.1.
Theorem 4.0.1** (Learning a GMM with Known Identical Covariance).**
Suppose and let . Let be an -dimensional GMM, i.e. , , and . Let be the matrix whose column is . If there exists so that , then Algorithm 6 recovers each to within accuracy with probability . Its sample and time complexity are at most
[TABLE]
where , , r\geq\big{(}\max_{i}{\lVert\mu_{i}\rVert}+1)/(\min_{i}{\lVert\mu_{i}\rVert})\big{)}, are bounds provided to the algorithm, and .
We note that the requirement that is due to the invocation of Theorem 1.3 of [48]; it should not be difficult, however, to adapt the algorithm to use a method similar to that of [54] to handle the case where .
Given that the means have been estimated, the weights can be recovered using the tensor structure of higher order cumulants (see Section 4.1 for the definition of cumulants). This is shown in Appendix 4.5.
We show that is large in the smoothed analysis sense, namely, if we start with a base matrix and perturb each entry randomly to get , then is likely to be large:
Theorem 4.0.2**.**
For , let be an arbitrary matrix. Let be a randomly sampled matrix with each entry iid from , for . Then, for some absolute constant , {vstd}
[TABLE]
We point out the simultaneous and independent work of [19], where the authors prove learnability results related to our Theorems 4.0.1 and 4.0.2. We now provide a comparison. The results in [19], which are based on tensor decompositions, are stronger in that they can learn mixtures of axis-aligned Gaussians (with non-identical covariance matrices) without requiring to know the covariance matrices in advance. Their results hold under a smoothed analysis setting similar to ours. To learn a mixture of roughly Gaussians up to an accuracy of their algorithm has running time and sample complexity and succeeds with probability at least , where the means are perturbed by adding an -dimensional Gaussian from . On the one hand, the success probability of their algorithm is much better (as a function of , exponentially close to as opposed to polynomially close to , as in our result). On the other hand, this comes at a high price in terms of the running time and sample complexity: The polynomial above has degree exponential in , unlike the degree of our bound which is polynomial in . Thus, in this respect, the two results can be regarded as incomparable points on an error vs running time (and sample complexity) trade-off curve. Our result is based on a reduction from learning GMMs to ICA which could be of independent interest given that both problems are extensively studied in somewhat disjoint communities. Moreover, our analysis can be used in the reverse direction to obtain hardness results for ICA.
The technique of Poissonization is, to the best of our knowledge, new in the GMM setting, though we note that it has been previously applied in computational learning. For instance, in [ValiantV13], it has been used for the estimation of properties of discrete probability distributions such as the support size and the entropy.
Finally, in Section 4.4 we show that in low dimension the situation is very different from the high-dimensional generic efficiency given by Theorems 4.0.1 and 4.0.2: The problem is generically hard. More precisely, we show:
Theorem 4.0.3**.**
Let be a set of points uniformly sampled from . Then with high probability there exist two mixtures with equal number of unit Gaussians , centered on disjoint subsets of , such that, for some , {vstd}
[TABLE]
Here we would like to note that the assumption that is random is convenient as it provides a natural model for genericity, guarantees (with high probability) small fill (Section 4.4) and also ensures that the means of the Gaussian mixture components of and are not too close. In particular, it is not difficult to verify that with high probability any pair of the random means are at least separated. However the randomness assumption is not essential. In fact, the above theorem will hold for an arbitrary set of points with sufficiently small fill (see Section 4.4).
Combining the above lower bound with our reduction provides a similar lower bound for ICA; see a discussion on the connection with ICA below. Our lower bound gives an information-theoretic barrier. This is in contrast to conjectured computational barriers that arise in related settings based on the noisy parity problem (see [54] for pointers). The only previous information-theoretic lower bound for learning GMMs we are aware of is due to [74] and holds for two specially designed one-dimensional mixtures.
4.1 Preliminaries
The singular values of a matrix will be ordered in the decreasing order: . By we mean .
For a real-valued random variable , the cumulants of are certain polynomials in the moments of . For , the th cumulant is denoted . Denoting , we have, for example: and . In general, cumulants can be defined as certain coefficients of a Taylor expansion of the logarithm of the moment generating function of : . The first two cumulants are the same as the expectation and the variance, resp. Cumulants have the property that for two independent random variables we have (assuming that the first moments exist for both and ). Cumulants are degree- homogeneous, i.e. if and is a random variable, then . The third and higher cumulants of the Gaussian distribution are 0.
4.1.1 Gaussian Mixture Model.
For , define Gaussian random vectors with distribution where and . Let be an integer-valued random variable which takes on value with probability , henceforth called weights. (Hence .) Then, the random vector drawn as is said to be a Gaussian Mixture Model (GMM) . The sampling of can be interpreted as first picking one of the components according to the weights, and then sampling a Gaussian vector from component . We will be primarily interested in the mixture of identical Gaussians of known covariance. In particular, there exists known such that for each . Letting , and denoting by the random variable which takes on the th canonical vector with probability , we can write the GMM model as follows:
[TABLE]
In this formulation, acts as a selector of a Gaussian mean. Conditioning on , we have , which is consistent with the GMM model.
Given samples from the GMM, the goal is to recover the unknown parameters of the GMM, namely the means and the weights .
4.1.2 Underdetermined ICA.
In the basic formulation of ICA, the observed random variable is drawn according to the model , where is a latent random vector whose components are independent random variables, and is an unknown mixing matrix. The probability distributions of the are unknown except that they are not Gaussian. The ICA problem is to recover to the extent possible. The underdetermined ICA problem corresponds the case . We cannot hope to recover fully because if we flip the sign of the th column of , or scale this column by some nonzero factor, then the resulting mixing matrix with an appropriately scaled will again generate the same distribution on as before. There is an additional ambiguity that arises from not having an ordering on the coordinates : If is a permutation matrix, then gives a new random vector with independent reordered coordinates, gives a new mixing matrix with reordered columns, and provides the same samples as since is the inverse of . As is a permutation of the columns of , this ambiguity implies that we cannot recover the order of the columns of . However, it turns out that under certain genericity requirements, we can recover up to these necessary ambiguities, that is to say we can recover the directions (up to sign) of the columns of , even in the underdetermined setting.
It will be important for us to work with an ICA model where there is Gaussian noise in the data: , where is an additive Gaussian noise independent of , and the covariance of given by is in general unknown and not necessarily spherical. We will refer to this model as the noisy ICA model.
We define the flattening operation from a tensor to a vector in the natural way. Namely, when is a tensor, then where is a bijection with indices running from to . Roughly speaking, each index is being converted into a digit in a base number up to the final offset by 1. This is the same flattening that occurs to go from a tensor outer product of vectors to the Kronecker product of vectors.
The ICA algorithm from [48] to which we will be reducing learning a GMM relies on the shared tensor structure of the derivatives of the second characteristic function and the higher order multi-variate cumulants. This tensor structure motivates the following form of the Khatri-Rao product:
Definition 7**.**
Given matrices , a column-wise Khatri-Rao product is defined by , where is the th column of , is the th column of , denotes the Kronecker product and is flattening of the tensor into a vector. The related Khatri-Rao power is defined by ( times).
This form of the Khatri-Rao product arises when performing a change of coordinates under the ICA model using either higher order cumulants or higher order derivative tensors of the second characteristic function.
4.1.3 ICA Results.
Theorem 4.1.1 (Theorem 4.1.1, from [48]) allows us to recover up to the necessary ambiguities in the noisy ICA setting. The theorem establishes guarantees for an algorithm from [48] for noisy underdetermined ICA, UnderdeterminedICA. This algorithm takes as input a tensor order parameter , number of signals , access to samples according to the noisy underdetermined ICA model with unknown noise, accuracy parameter , confidence parameter , bounds on moments and cumulants and , a bound on the conditioning parameter , and a bound on the cumulant order . It returns approximations to the columns of up to sign and permutation.
ICA Results.
The following theorem (from [48]) allows us to recover up to the necessary ambiguities in the noisy ICA setting. The theorem establishes guarantees for an algorithm from [48] for noisy underdetermined ICA, UnderdeterminedICA. This algorithm takes as input a tensor order parameter , number of signals , access to samples according to the noisy underdetermined ICA model with unknown noise, accuracy parameter , confidence parameter , bounds on moments and cumulants and , a bound on the conditioning parameter , and moment order . It returns approximations to the columns of up to sign and permutation.
Theorem 4.1.1** ([48]).**
*Let a random vector be given by an underdetermined ICA model with unknown Gaussian noise where , with unit norm columns, and the covariance matrix are unknown. Let be such that . Let be such that for each , there is a satisfying and , and . Moreover, suppose that the noise also satisfies the same moment condition: for any unit vector (this is satisfied if we have where is the maximum eigenvalue of ). Then algorithm UnderdeterminedICA returns a set of -dimensional vectors so that for some permutation of and signs we have for all . Its sample and time complexity are . *
4.2 Learning GMM means using underdetermined ICA: The basic idea
In this section we give an informal outline of the proof of our main result, namely learning the means of the components in GMMs via reduction to the underdetermined ICA problem. Our reduction will be discussed in two parts. The first part gives the main idea of the reduction and will demonstrate how to recover the means up to their norms and signs, i.e. we will get . We will then present the reduction in full. It combines the basic reduction with some preprocessing of the data to recover the ’s themselves. The reduction relies on some well-known properties of the Poisson distribution stated in the lemma below; its proof can be found in Appendix 4.6.4.
Lemma 23**.**
Fix a positive integer , and let be such that . If and then for all and are mutually independent.
Basic Reduction: The main idea.
Recall the GMM from equation (4.1) is given by . Henceforth, we will set . We can write the GMM in the form , which is similar in form to the noisy ICA model, except that does not have independent coordinates. We now describe how a single sample of an approximate noisy ICA problem is generated.
The reduction involves two internal parameters and that we will set later. We generate a Poisson random variable , and we run the following experiment times: At the th step, generate sample from the GMM. Output the sum of the outcomes of these experiments: .
Let be the random variable denoting the number of times samples from the th component were chosen in the above experiment. Thus . Note that are not observable although we know their sum. By Lemma 23, each has distribution , and the random variables are mutually independent. Let .
For a non-negative integer , we define where the are iid according to . In this definition can be a r.v., in which case the are sampled independent of . Using to indicate that two random variables have the same distribution, then
[TABLE]
If there were no Gaussian noise in the GMM (i.e. if we were sampling from a discrete set of points) then the model in (4.2) becomes , which is the ICA model without noise, and so we could recover up to necessary ambiguities. However, the model fails to satisfy even the assumptions of the noisy ICA model, both because is not independent of and because is not distributed as a Gaussian random vector.
As the covariance of the additive Gaussian noise is known, we may add additional noise to the samples of to obtain a good approximation of the noisy ICA model. Parameter , the second parameter of the reduction, is chosen so that with high probability we have . Conditioning on the event we draw according to the rule
[TABLE]
where , , and are drawn independently conditioned on . Then, conditioned on , we have .
Note that we have only created an approximation to the ICA model. In particular, restricting can be accomplished using rejection sampling, but the coordinate random variables would no longer be independent. We have two models of interest: (1) , a noisy ICA model with no restriction on , and (2) the restricted model.
We are unable to produce samples from the first model, but it meets the assumptions of the noisy ICA problem. Pretending we have samples from model (1), we can apply Theorem 4.1.1 to recover the Gaussian means up to sign and scaling. On the other hand, we can produce samples from model (2), and depending on the choice of , the statistical distance between models (1) and (2) can be made arbitrarily close to zero. It will be demonstrated that given an appropriate choice of , running UnderdeterminedICA on samples from model (2) is equivalent to running UnderdeterminedICA on samples from model (1) with high probability, allowing for recovery of the Gaussian mean directions up to some error.
Full reduction.
To be able to recover the without sign or scaling ambiguities, we add an extra coordinate to the GMM as follows. The new means are with an additional coordinate whose value is for all , i.e. \mu_{i}^{\prime}:=\left(\begin{array}[]{c}\mu_{i}\\ 1\end{array}\right). Moreover, this coordinate has no noise. In other words, each Gaussian component now has an covariance matrix \Sigma^{\prime}:=\left(\begin{array}[]{cc}\Sigma&0\\ 0&0\end{array}\right). It is easy to construct samples from this new GMM given samples from the original: If the original samples were , then the new samples are where u^{\prime}_{i}:=\left(\begin{array}[]{c}u_{i}\\ 1\end{array}\right). The reduction proceeds similarly to the above on the new inputs.
Unlike before, we will define the ICA mixing matrix to be such that it has unit norm columns. The role of matrix in the basic reduction will now be played by . Since we are normalizing the columns of , we have to scale the ICA signal obtained in the basic reduction to compensate for this: Define . Thus, the ICA models obtained in the full reduction are:
[TABLE]
where we define \eta^{\prime}(\tau)=\left(\begin{array}[]{c}\eta(\tau)\\ 0\end{array}\right). As before, we have an ideal noisy ICA model (4.3) from which we cannot sample, and an approximate noisy ICA model (4.4) which can be made arbitrarly close to (4.3) in statistical distance by choosing appropriately. With appropriate application of Theorem 4.1.1 to these models, we can recover estimates (up to sign) of the columns of .
By construction, the last coordinate of now tells us both the signs and magnitudes of the : Let be the vector consisting of the first coordinates of , and let be the last coordinate of . Then
[TABLE]
with the sign indeterminancy canceling in the division.
4.3 Correctness of the Algorithm and Reduction
Subroutine 5 is used to capture the main process of the reduction: Let be the covariance matrix of the GMM, be an integer chosen as input, and a threshold value also computed elsewhere and provided as input. Let . If is larger , the subroutine returns a failure notice and the calling algorithm halts immediately. A requirement, then, should be that the threshold is chosen so that the chance of failure is very small; in our case, is chosen so that the chance of failure is half of the confidence parameter given to Algorithm 6. The algorithm then goes through the process described in the full reduction: sampling from the GMM, lifting the sample by appending a 1, then adding a lifted gaussian noise so that the total noise has distribution . The resulting sample is from the model given by (4.4).
Algorithm 6 works as follows: it takes as input the parameters of the GMM (covariance matrix, number of means), tensor order (as required by UnderdeterminedICA), error parameters, and bounds on certain properties of the weights and means. The algorithm then calculates various internal paremeters: a bound on directional covariances, error parameters to be split between the “Poissonization” process and the call to UnderdeterminedICA, the threshold parameter and Poisson parameter to be used in Subroutine 5, and values explicitly needed by the proof of UnderdeterminedICA in [48]. Other internal values are given by a constanct and a polynomial ; these are determined by the proof of Theorem 4.0.1. Essentially, is a constant so that one can cleanly compute a value of that will involve a polynomial, , of all the other parameters. The algorithm then calls UnderdeterminedICA, but instead of giving samples from the GMM, it allows access to Subroutine 5. It is then up to UnderdeterminedICA to generate samples as needed (bounded by the polynomial in Theorem 4.0.1). In the case that Subroutine 5 returns a failure, the entire algorithm process halts, and returns nothing. If no failure occurs, the matrix returned by UnderdeterminedICA will be the matrix of normalized means embedded in , and the algorithm de-normalizes, removes the last row, and then has approximations to the means of the GMM.
The bounds are used instead of actual values to allow flexibility – in the context under which the algorithm is invoked – on what the algorithm needs to succeed. However, the closer the bounds are to the actual values, the more efficient the algorithm will be.
We show that the reduction results in a model that allows one to use the ICA algorithm UnderdeterminedICA presented in [48] (see Section 2.2, Theorem 4.1.1). That is, we can choose the parameters, and , for the Poissonization and threshold that yield a model which will satisfy the assumptions of Theorem 4.1.1. This also involves translating the bounds on cumulants and higher moments of our Poisson signals and added noise to obtain feasible bounds for Theorem 4.1.1.
Sketch of correctness argument.
The proof of correctness of Algorithm 6 has two main parts. In the first part (section 4.3), we analyze the sample complexity of recovering the Gaussian means using UnderdeterminedICA when samples are taken from the ideal noisy ICA model (4.3). In the second part, we note that we do not have access to the ideal model (4.3), and that we can only sample from the approximate noisy ICA model (4.4) using the full reduction. Choosing appropriately, we use total variation distance to argue that with high probability, running UnderdeterminedICA with samples from the approximate noisy ICA model will produce equally valid results as running UnderdeterminedICA with samples from the ideal noisy ICA model. The total variation distance bound is explored in section 4.3.
These ideas are combined in section 4.3 to prove the correctness of Algorithm 6. One additional technicality arises from the implementation of Algorithm 6. Samples can be drawn from the noisy ICA model using rejection sampling on . In order to guarantee Algorithm 6 executes in polynomial time, when a sample of needs to be rejected, Algorithm 6 terminates in explicit failure. To complete the proof, we argue that with high probability, Algorithm 6 does not explicitly fail.
Error Analysis of the Ideal Noisy ICA Model
The proposed full reduction from Section 4.2 provides us with two models. The first is a noisy ICA model from which we cannot sample:
[TABLE]
The second is a model that fails to satisfy the assumption that has independent coordinates, but it is a model from which we can sample:
[TABLE]
Both models rely on the choice of two parameters, and . The dependence on is explicit in the models. The dependence on can be summarized in the unrestricted model as independently of each other, and .
The probability of choosing will be seen to be exponentially small in . For this reason, running UnderdeterminedICA with polynomially many samples from model (4.5) will with high probability be equivalent to running the ICA Algorithm with samples from model (4.6). This notion will be made precise later using total variation distance.
For the remainder of this subsection, we proceed as if samples are drawn from the ideal noisy ICA model (4.5). Thus, to recover the columns of , it suffices to run UnderdeterminedICA on samples of . Theorem 4.1.1 can be used for this analysis so long as we can obtain the necessary bounds on the cumulants of , moments of , and the moments of . We define and . Then, the cumulants of are bounded by the following lemma:
Lemma 24**.**
Given , for each . In particular, then .
Proof.
By construction, . By the homogeneity property of univariate cumulants,
[TABLE]
As , . The cumulants of the Poisson distribution are given in Lemma 37. It follows that . ∎
The bounds on the moments of for each can be computed using the following lemma:
Lemma 25**.**
For , we have .
Proof.
Let denote a random variable drawn from . It is known (see [86]) that
[TABLE]
where denotes Stirling number of the second kind. Using Lemma 36, it follows that
[TABLE]
Since where , it follows that . ∎
The absolute moments of Gaussian random variables are well known. For completeness, the bounds are provided in Lemma 38 of Appendix 4.6.6.
Defining ; vectors , , and similarly and for later; and choosing , we can now show a polynomial bound for the error in recovering the columns of using UnderdeterminedICA.
Theorem 4.3.1** (ICA specialized to the ideal case).**
Suppose that samples of are taken from the unrestricted ICA model (4.5) choosing parameter and a constant. Suppose that UnderdeterminedICA is run using these samples. Suppose . Fix and . Then with probability , when the number of samples is:
[TABLE]
the columns of are recovered within error up to their signs. That is, denoting the columns returned from UnderdeterminedICA by , there exists and a permutation of such that for each .
Proof.
Obtaining the sample bound is an exercise of rewriting the parameters associated with the model in a way which can be used by Theorem 4.1.1. In what follows, where new parameters are introduced without being described, they will correspond to parameters of the same name defined in and used by the statement of Theorem 4.1.1.
Parameter is fixed. We must choose and such that and is bounded away from [math]. It suffices to choose . By Lemma 24, for each . As , we have that for each , giving a somewhat more natural condition number. In the notation of Theorem 4.1.1, we have a constant
[TABLE]
such that for each .
Now we consider the upper bound on the absolute moments of both and on . As the Poisson distribution takes on non-negative values, it follows that takes on non-negative values. Thus, the moments and absolute moments of coincide. Using Lemma 25, we have that . Thus, for to bound the th moment of , it suffices that . Noting that
[TABLE]
it suffices that , giving a more natural condition number.
Now we bound the absolute moments of the Gaussian distribution. As , it follows that is odd. Given a unit vector , it follows from Lemma 38 that
[TABLE]
gives a clear upper bound for , and gives a clear upper bound to . As such, it suffices that in order to guarantee that . Using the obtained bounds for from the Poisson and Normal variables, it suffices that be taken such that
[TABLE]
to guarantee that bounds all required order absolute moments.
We can now apply Theorem 4.1.1, using the parameter values , from (4.8), and from (4.9). Then with probability ,
[TABLE]
samples suffice to recover up to sign the columns of within accuracy. More precisely, letting give the columns produced by UnderdeterminedICA, then there exists parameters such that captures the sign indeterminacy, and a permutation on such that for each .
The poly bound in (4.10) is equivalent to the poly bound in (4.7). ∎
Theorem 4.3.1 allows us to recover the columns of up to sign. However, what we really want to recover are the means of the original Gaussian mixture model, which are the columns of . Recalling the correspondence between and laid out in section 4.2, the Gaussian means which form the columns of are related to the columns of by the rule . Using this rule, we can construct estimate the Gaussian means from the estimates of the columns of . By propogating the errors from Theorem 4.3.1, we arrive at the following result:
Theorem 4.3.2** (Recovery of Gaussian means in Ideal Case).**
Suppose that UnderdeterminedICA is run using samples of from the ideal noisy ICA model (4.5) choosing parameters and a constant. Define such that . Suppose further that . Let be the returned estimates of the columns of (from model (4.5)) by UnderdeterminedICA. Let for each . Fix error parameters and . When at least
[TABLE]
samples are used, then with probability there exists a permutation of such that for each .
Proof.
Let (to be chosen later) give a desired bound on the errors of the columns of . Then, from Theorem 4.3.1, using
[TABLE]
samples suffices with probability to produce column estimates such that for an unknown permutation and signs , give -close estimates of the columns respectively of . In order to avoid notational clutter, we will assume without loss of generality that is the identity map, and hence that holds.
This proof proceeds in two steps. First, we replace the dependencies in (4.13) on parameters from the lifted GMM model generated by the full reduction with dependencies based on the GMM model we are trying to learn. Then, we propagate the error from recovering the columns to that of recovering .
Step 1: GMM Dependency Replacements.
In the following two claims, we consider alternative lower bounds for for recovering column estimators which are -close up to sign to the columns of . In particular, so long as we use at least as many samples of as in (4.13) when calling UnderdeterminedICA, then will be recovered with the desired precision with probability .
Claim**.**
*The dependence in (4.13) can be replaced by a
dependence.*
Proof of Claim..
By construction, \mu_{\max}^{\prime}=\left(\begin{array}[]{c}\mu_{\max}\\ 1\end{array}\right). By the triangle inequality,
[TABLE]
where is a polynomial of with coefficients bounded by . The maximal power of in is . It follows that . ∎
Claim**.**
The in (4.13) can be replaced by a dependence.
Proof of Claim.
First define to be the unnormalized version of . That is, . Then, implies . Thus, .
Next, we note that \underline{A}^{\prime}=\left(\begin{array}[]{c}A\\ \mathbf{1}\end{array}\right) where is an all ones row vector. It follows that the rows of are a strict subset of the rows of . Thus,
[TABLE]
Finally, and . It follows that . Chaining together inequalities yields:
[TABLE]
As , the triangle inequality implies . As we require the dependency of at least samples, it suffices to have the replacement dependency of samples. ∎
Thus, it is sufficient to call UnderdeterminedICA with
[TABLE]
samples to achieve the desired accuracy on the returned estimates of the columns of with probability .
Step 2: Error propagation.
What remains to be shown is that an appropriate choice of enforces by propagating the error.
Recall that A_{i}^{\prime}=\left(\begin{array}[]{c}\mu_{i}\\ 1\end{array}\right)\cdot{\lVert\left(\begin{array}[]{c}\mu_{i}\\ 1\end{array}\right)\rVert}^{-1}, making . Thus,
[TABLE]
We have that:
[TABLE]
which follows in part by applying (4.16) for the left summand and noting that is a unit vector for the right summand, giving the bound . Continuing with the restriction that ,
[TABLE]
Then, in order to guarantee that , it suffices to choose such that
[TABLE]
which occurs when
[TABLE]
As , the restriction holds automatically for the choice of in (4.17). The sample bound from (4.15) contains the dependency . Propagating the error gives a replacement dependency of
[TABLE]
as is non-negative. This propagated dependency is reflected in (4.12). ∎
Distance of the Sampled Model to the Ideal Model
An important part of the reduction is that the coordinates of are mutually independent. Without the threshold , this is true (c.f. Lemma 23). However, without the threshold, one cannot know how to add more noise so that the total noise on each sample is iid. We show that we can choose the threshold large enough that the samples still come from a distribution with arbitrarily small total variation distance to the one with truly independent coordinates.
Lemma 26**.**
Fix . Let for . Let , If , , and , then .
Proof.
By the Chernoff bound (See Theorem A.1.15 in [5]),
[TABLE]
For any , letting , we get
[TABLE]
To get , it suffices that . Note that
[TABLE]
If , then we have
[TABLE]
which then implies it suffices that
[TABLE]
which holds for , giving the desired result.
∎
Lemma 27**.**
Let , , and be iid with distribution . If then
[TABLE]
Proof.
By Lemma 26 implies for every . The union bound gives us the desired result. ∎
We have now that if we choose our threshold large enough, our samples can be statistically close (See Appendix 4.6.7) to ones that would come from the truly independent distribution. This claim is made formal now:
Lemma 28**.**
Fix . Let . Let be a poisson distribution with parameter and have corresponding density . Let be a discrete distribution with density when and 0 otherwise. Then .
Proof.
Since we are working with discrete distributions, we can write
[TABLE]
Then we can compute
[TABLE]
∎
Proof of Theorem 4.0.1
We now show that after the reduction is applied, we can use the ICA routine given in [48] to learn the GMM. Instead of requiring exact values of each parameter, we simply require a bound on each. The algorithm remains polynomial on those bounds, and hence polynomial on the true values.
Proof.
The algorithm is provided parameters: Covariance matrix , upper bound on tensor order parameter , access to samples from a mixture of identical, spherical Gaussians in with covariance , confidence parameter , accuracy parameter , upper bound , upper bound on the norm of the mixture means , lower bound so , and r\geq\big{(}\max_{i}{\lVert\mu_{i}\rVert}+1)/(\min_{i}{\lVert\mu_{i}\rVert})\big{)}.
The algorithm then needs to fix the number of samples , sampling threshold , Poisson parameter , and two new errors and so that . For simplicity, we will take . Then fix for . Recall that is the matrix whose th column is . Let be the matrix whose th column is .
Step 1
Assume that after drawing samples from Subroutine 5, the signals are mutually independent (as in the “ideal” model given by (4.5)) and the mean matrix satisfies . Then by Theorem 4.3.2, with probability of error , the call to UnderdeterminedICA in Algorithm 6 recovers the columns of to within and up to a permutation using samples where
[TABLE]
where is the bound on promised by Theorem 4.3.2 and is all its arguments except the dependence in . So then we have that with at least samples in this “ideal” case, we can recover approximations to the true means in up to a permutation and within distance.
Step 2
We need to show that after getting samples from the reduction, the resulting distribution is still close in total variation to the independent one. We will choose a new . Let . Given , Lemma 28 shows that for , with probability , .
Take iid random variables with distribution and density . Let be a distribution with density function . Let be iid random variables with distribution . Denote the joint distribution of the ’s by and the joint distribution of the ’s as . By the union bound and the fact that total variation distance satisfies the triangle inequality,
[TABLE]
Then for our choice of , by Lemma 26 and Lemma 28, we have
[TABLE]
By the same union bound argument, the probability that the algorithm fails (when ) is at most , since it has to draw samples. So with high probability, the algorithm does not fail; otherwise, it still does not take more than polynomial time, and will terminate instead of returning a false result.
Step 3
We know that is at least a polynomial which can be written in terms of the dependence on as . This means there will be a power of which dominates all of the factors in , and in particular, will be for some . It then suffices to choose so that , where
[TABLE]
Then, with the proper choice of (to be specified shortly), from step 2 we have
[TABLE]
Since it suffices to choose so that
[TABLE]
Finally, we claim that
[TABLE]
is enough for the desired bound on the sample size. Observe that .
An useful fact is that for general , satisfies . This captures the essence of our situation nicely. Letting play the role of , play the role of and play the role of , to satisfy (4.19), it suffices that
[TABLE]
We can see that and by construction. But we also get which implies . Thus for our choice of , which also preserves the requirement in Step 2, there is a corresponding set of choices for , where the required sample size remains polynomial as
[TABLE]
where we used the bound . ∎
4.4 Smoothed Analysis
In this section, we prove that in high enough dimension, a Gaussian Mixture will satisfy the non-degeneracy conditions discussed above. We start with a base matrix and add a perturbation matrix with each entry coming iid from for some . (We restrict the discussion to the second power for simplicity; extension to higher power is straightforward.) As in [48], it will be convenient to work with the multilinear part of the Khatri–Rao product: For a column vector define , a subvector of , given by for . Then for a matrix we have .
Theorem 4.4.1**.**
With the above notation, for any base matrix with dimensions as above, we have, for some absolute constant , {vstd}
[TABLE]
Theorem 4.0.2 follows by noting that .
Proof.
In the following, for a vector space (over the reals) denotes the distance between vector and subspace ; more precisely, . We will use a lower bound on , found in Appendix 4.6.2.
With probability , the columns of the matrix are linearly independent. This can be proved along the lines of a similar result in [48]. Fix and let be a unit vector orthogonal to the subspace spanned by the columns of other than column . Vector is well-defined with probability . Then the distance of the ’th column from the span of the rest of the columns is given by
[TABLE]
Now note that this is a quadratic polynomial in the random variables . We will apply the anticoncentration inequality of \processifversionvstdCarbery–Wright [23] to this polynomial to conclude that the distance between the ’th column of and the span of the rest of the columns is unlikely to be very small (see Appendix 4.6.3 for the precise result).
Using , the variance of our polynomial in (4.21) becomes
[TABLE]
In our application, random variables for are not standard Gaussians but are iid Gaussian with variance , and our polynomial does not have unit variance. After adjusting for these differences using the estimate on the variance of above, Lemma 32 gives . Therefore, by the union bound over the choice of .
Now choosing , Lemma 31 gives . ∎
We note that while the above discussion is restricted to Gaussian perturbation, the same technique would work for a much larger class of perturbations. To this end, we would require a version of the Carbery-Wright anticoncentration inequality which is applicable in more general situations. We omit such generalizations here.
The curse of low dimensionality for Gaussian mixtures
In this section we prove that for small there is a large class of superpolynomially close mixtures in with fixed variance. This goes beyond the specific example of exponential closeness given in [74] as we demonstrate that such mixtures are ubiquitous as long as there is no lower bound on the separation between the components.
Specifically, let be a cube . We will show that for any two sets of points and in , with fill (we say that has fill , if there is a point of within distance of any point of ), there exist two mixtures with means on non-overlapping subsets of , which are (nearly) exponentially close in in the norm. Note that the fill of a sample from the uniform distribution on the cube can be bounded (with high probability) by .
We start by defining some of the key objects. Let be the unit Gaussian kernel. Let be the integral operator corresponding to the convolution with a unit Gaussian: . Let be any subset of points in . Let be the kernel matrix corresponding to , . It is known to be positive definite. The interpolant is defined as , where the coefficients are chosen so that . It is easy to see that such interpolant exists and is unique, obtained by solving a linear system involving .
We will need some properties of the Reproducing Kernel Hilbert Space corresponding to the kernel . In particular, we need the bound and the reproducing property, .
Lemma 29**.**
Let be any positive function with norm supported on and let . If fill is , there exists such that
[TABLE]
Proof.
From [85], Corollary 5.1 (taking ) we have that for some and sufficiently small
[TABLE]
and
[TABLE]
Note that the norm is on while we need to control the norm on .
That, however, does not allow us to directly control the function behavior outside of the interval. To do that we need a bound on the RKHS norm of . We first observe for any , . Thus, from the reproducing property of RKHS, . Using properties of RKHS with respect to the operator (see, e.g., Proposition 10.28 of [109])
[TABLE]
Thus
[TABLE]
∎
Theorem 4.4.2**.**
Let and be any two subsets of with fill . Then there exist two Gaussian mixtures and (with positive coefficients summing to one, but not necessarily the same number of components), which are centered on two non-intersecting subsets of and such that
[TABLE]
for some constant .
Proof.
To simplify the notation we assume that . The general case follows verbatim, except that the interval of integration, , and its complement need to be replaced by the sphere of radius and its complement respectively.
Let and be the interpolants as above. We see that . and are both linear combinations of Gaussians possibly with negative coefficients and so is . By collecting positive and negative coefficients we write , where, and are mixtures with positive coefficients only.
Put , , where and are non-overlapping subsets of . Now we need to ensure that the coefficients can be normalized to sum to .
Let , . By integrating over the interval , and since is strictly positive on the interval, it is easy to see that (and by the same token ), where is some universal constant. We have
[TABLE]
[TABLE]
Noticing that the first summand is bounded by and the integral in the second summand is smaller than for sufficiently small, it follows immediately, that for some and sufficiently small.
Hence, we have
[TABLE]
Collecting exponential inequalities completes the proof.
∎
Theorem 4.0.3.
For convenience we will use a set of points instead of . Clearly it does not affect the exponential rate. By a simple covering set argument (cutting the cube into cubes with size ) and basic probability (the coupon collector’s problem), we see that the fill of points is at most with probability . Hence, given points, we have . We see that with a smaller probability (but still close to for large ), we can sample points times and still have the same fill on each group of .
Pairing the sets of points into pairs of sets arbitrarily and applying Theorem 4.4.2 (to points) we obtain pairs of exponentially close mixtures with at most components each. If one of the pairs has the same number of components, we are done. If not, by the pigeon-hole principle for at least two pairs of mixtures and the differences of the number of components (an integer number between [math] and ) must coincide. Assume without loss of generality that has no more components than and has no more components than .Taking and completes the proof. ∎
4.5 Recovery of Gaussian Weights
Multivariate cumulant tensors and their properties.
Our technique for the recovery of the Gaussian weights relies on the tensor properties of multivariate cumulants that have been used in the ICA literature.
Given a random vector , the moment generating function of is defined as . The cumulant generating function is the logarithm of the moment generating function: .
Similarly to the univariate case, multivariate cumulants are defined using the coefficients of the Taylor expansion of the cumulant generating function. We use both and to denote the order- cross cumulant between the random variables . Then, the cross-cumulants are given by the formula \kappa_{Y}^{j_{1},\dots,j_{\ell}}=\frac{\partial}{\partial t_{j_{1}}}\cdots\frac{\partial}{\partial t_{j_{\ell}}}g_{\rule{0.0pt}{4.52083pt}Y}(t)\big{|}_{t=0}. When unindexed, will denote the full order- tensor containing all cross-cumulants, with the order of the tensor being made clear by context. In the special case where , we obtain the order- univariate cumulant ( repeated times) previously defined. We will use some well known properties of multivariate cumulants, found in Appendix 4.6.1.
The most theoretically justified ICA algorithms have relied on the tensor structure of multivariate cumulants, including the early, popular practical algorithm JADE [25]. In the fully determined ICA setting in which the number source signals does not exceed the ambient dimension, the papers [14] and [17] demonstrate that ICA with additive Gaussian noise can be solved in polynomial time and using polynomial samples. The tensor structure of the cumulants was (to the best of our knowledge) first exploited in [24] and later in [4] to solve underdetermined ICA. Finally, [48] provides an algorithm with rigorous polynomial time and sampling bounds for underdetermined ICA in the presence of Gaussian noise.
Weight recovery (main idea).
Under the basic ICA reduction (see section 4.2) using the Poisson distribution with parameter , we have that is observed such that and . As has already been recovered, what remains to be recovered are the weights . These can be recovered using the tensor structure of higher order cumulants. The critical relationship is captured by the following Lemma:
Lemma 30**.**
Suppose that gives a noisy ICA model. When is of order , then .
Proof.
It is easily seen that the Gaussian component has no effect on the cumulant:
[TABLE]
Then, we expand :
[TABLE]
But, by independence, whenever fails to hold. Thus,
[TABLE]
Flattening yields: . ∎
In particular, we have that with the probability of sampling from the th Gaussian. Given knowledge of and the cumulants of the Poisson distribution, we can recover the Gaussian weights.
Theorem 4.5.1**.**
Suppose that is the unrestricted noisy ICA model from the basic reduction (see section 4.2). Let be such that has linearly independent columns, and let be its Moore-Penrose pseudoinverse. Let be of order . Then is the vector of mixing weights of the Gaussian mixture model.
Proof.
From Lemma 37, . Then Lemma 30 implies that . Multiplying on the left by gives the result. ∎
4.6 Addendum
4.6.1 Properties of Cumulants
The following properties of multivariate cumulants are well known and are largely inherited from the definition of the cumulant generating function:
- •
(Symmetry) Let give a permutation of indices. Then, .
- •
(Multilinearity of coordinate random variables) Given constants , then
[TABLE]
Also, given a scalar random variable , then
[TABLE]
with symmetry implying the additive multilinear property for all other coordinates.
- •
(Independence) If there exists such that and are independent random variables, then the cross-cumulant . Combined with multilinearity, it follows that when there are two independent random vectors and , then .
- •
(Vanishing Gaussians) When , then for the Gaussian random variable , .
4.6.2 Rudelson-Vershynin subspace bound
Lemma 31** (\processifversionvstdRudelson–Vershynin [90]).**
If has columns , then denoting , we have
[TABLE]
where as usual .
4.6.3 Carbery-Wright anticoncentration
The version of the anticoncentration inequality we use is explicitly given in [75] which in turn follows immediately from [23]:
Lemma 32** ([75]).**
Let be a multilinear polynomial of degree . Suppose that when for all . Then there exists an absolute constant such that for and ,
[TABLE]
4.6.4 Lemmas on the Poisson Distribution
The following lemmas are well-known; see, e.g., [35]. We provide proofs for completeness.
Lemma 33**.**
If and then .
Proof.
[TABLE]
∎
Lemma 34**.**
Fix a positive integer , and let be such that . If and then for all and are mutually independent.
Proof.
The first part of the lemma (i.e., for all ) follows from Lemma 33. For the second part, let’s prove it for the binomial case (); the general case is similar.
[TABLE]
∎
4.6.5 Bounds on Stirling Numbers of the Second Kind
The following bound comes from [84] Theorem 3.
Lemma 35**.**
If and are integers, then .
From this, we can derive a somewhat looser bound on the Stirling numbers of the second kind which does not depend on :
Lemma 36**.**
If such that , then .
Proof.
The Stirling number of the second kind gives a count of the number of ways of splitting a set of labeled objects into unlabeled subsets. In the case where , then As , it is clear that for these choices of and , . By the restriction , when , then giving that . As such, the only remaining cases to consider are when and , the cases where Lemma 35 applies.
When and , then
[TABLE]
which is slightly stronger than the desired upper bound. ∎
4.6.6 Values of Higher Order Statistics
In this appendix, we gather together some of the explicit values for higher order statistics of the Poisson and Normal distributions required for the analysis of our reduction from learning a Gaussian Mixture Model to learning an ICA model from samples.
Lemma 37** (Cumulants of the Poisson distribution).**
Let . Then, for every positive integer .
Proof.
The moment generating function of the Poisson distribution is given by . The cumulant generating function is thus . The th derivative is given by .
By definition, . ∎
Lemma 38** (Absolute moments of the Gaussian distribution).**
For a Gaussian random variable , the absolute moments of are given by:
[TABLE]
The case that is even in Lemma 38 is well known, and can be found for instance in [61, Section 3.4]. For general , it is known (see [110]) that
[TABLE]
When is odd, is an integer, allowing the Gamma function to simplify to a factorial: . This gives the case where is odd in Lemma 38.
4.6.7 Total Variation Distance
Total variation is a type of statistical distance metric between probability distributions. In words, the total variation between two measures is the largest difference between the measures on a single event. Clearly, this distance is bounded above by 1.
For probability measures and on a sample space with sigma-algebra , the total variation is denoted and defined as:
[TABLE]
Equivalently, when and are distribution functions having densities and , respectively,
[TABLE]
where is an arbitrary positive measure for which and are absolutely continuous.
More specifically, when and are discrete distributions with known densities, we can write
[TABLE]
where we choose that simply assigns unit measure to each atom of (in this case, absolute continuity is trivial since only when is empty and thus must also be 0). For more discussion, one can see Definition 15.3 in [77] and Sect. 11.6 in [88].
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] D. Achlioptas and F. Mc Sherry. On spectral learning of mixture of distributions. In The 18th Annual Conference on Learning Theory , 2005.
- 2[2] Radoslaw Adamczak, Alexander Litvak, Alain Pajor, and Nicole Tomczak-Jaegermann. Quantitative estimates of the convergence of the empirical covariance matrix in logconcave ensembles. J. Amer. Math. Soc. , 233:535–561, 2011.
- 3[3] Radosław Adamczak, Alexander E. Litvak, Alain Pajor, and Nicole Tomczak-Jaegermann. Quantitative estimates of the convergence of the empirical covariance matrix in log-concave ensembles. J. Amer. Math. Soc. , 23(2):535–561, 2010.
- 4[4] Laurent Albera, Anne Ferréol, Pierre Comon, and Pascal Chevalier. Blind identification of overcomplete mixtures of sources (biome). Linear algebra and its applications , 391:3–30, 2004.
- 5[5] Noga Alon and Joel H Spencer. The probabilistic method . Wiley, 2004.
- 6[6] Shun Amari, Andrzej Cichocki, Howard H Yang, et al. A new learning algorithm for blind signal separation. Advances in neural information processing systems , pages 757–763, 1996.
- 7[7] Anima Anandkumar, Rong Ge, Daniel Hsu, Sham M. Kakade, and Matus Telgarsky. Tensor decompositions for learning latent variable models. Co RR , abs/1210.7559, 2012.
- 8[8] Joseph Anderson, Mikhail Belkin, Navin Goyal, Luis Rademacher, and James Voss. The more, the merrier: the blessing of dimensionality for learning large gaussian mixtures. In Proceedings of The 27th Conference on Learning Theory , pages 1135–1164, 2014.