Heavy-Tailed Analogues of the Covariance Matrix for ICA

TL;DR

This paper introduces HTICA, a practical and efficient algorithm for heavy-tailed ICA that outperforms existing methods, and analyzes how heavy tails influence various ICA algorithms through theoretical and experimental studies.

Contribution

The paper presents HTICA, a new practical algorithm for heavy-tailed ICA based on the centroid body, and provides a comprehensive analysis of heavy tails' impact on ICA algorithms.

Findings

HTICA outperforms existing algorithms in heavy-tailed regimes.

Some algorithms using covariance or higher moments can solve ICA with infinite second moments.

Explicit analytic representations of the centroid body improve practical efficiency.

Abstract

Independent Component Analysis (ICA) is the problem of learning a square matrix $A$, given samples of $X=AS$, where $S$ is a random vector with independent coordinates. Most existing algorithms are provably efficient only when each $S_i$ has finite and moderately valued fourth moment. However, there are practical applications where this assumption need not be true, such as speech and finance. Algorithms have been proposed for heavy-tailed ICA, but they are not practical, using random walks and the full power of the ellipsoid algorithm multiple times. The main contributions of this paper are: (1) A practical algorithm for heavy-tailed ICA that we call HTICA. We provide theoretical guarantees and show that it outperforms other algorithms in some heavy-tailed regimes, both on real and synthetic data. Like the current state-of-the-art, the new algorithm is based on the centroid body (a first moment analogue of the covariance matrix). Unlike the state-of-the-art, our algorithm is practically efficient. To achieve this, we use explicit analytic representations of the centroid body, which bypasses the use of the ellipsoid method and random walks. (2) We study how heavy tails affect different ICA algorithms, including HTICA. Somewhat surprisingly, we show that some algorithms that use the covariance matrix or higher moments can successfully solve a range of ICA instances with infinite second moment. We study this theoretically and experimentally, with both synthetic and real-world heavy-tailed data.

Figures

Equations

$$K_{X_R}=E\exp (-{\Vert X\Vert }^2/R^2)$$

$${\rho }_{X_R}(x)={\rho }_X(x)\exp (-{\Vert x\Vert }^2/R^2)/K_{X_R}.$$

$$P(|L_{ij}|>{ϵ}_1)\le \frac{8M^2}{{ϵ}_1^2{N}^{\gamma /3}}$$

$$P({\Vert L-D^{-1}\Vert }_2<n{ϵ}_1)$$

$$\ge P(\forall i\ne j,|L_{ij}|\le {ϵ}_1)$$

$$\ge 1-\frac{8n^2{M}^2}{{ϵ}_1^2{N}^{\gamma /3}}$$

$${\Vert D\Vert }_2=\frac{1}{{\sigma }_{\min }(D^{-1})}=\frac{1}{{\min }_{i\in [n]}{L}_{ii}}$$

$$P\left(\sum_{k\in[N]}\mathbbm{1}({s_{i}^{(k)}}^{2}\geq t^{2})\leq\frac{Np}{2}\right)\leq\exp\bigg{(}-\frac{Np}{8}\bigg{)}$$

$$N_{2}:=\max\bigg{\{}\bigg{(}\frac{24n^{2}M^{2}}{\epsilon_{1}^{2}\delta}\bigg{)}^{3/\gamma},\bigg{(}\frac{8M^{2}}{\epsilon_{1}}\bigg{)}^{\frac{1}{2}+\frac{1}{\gamma}}\bigg{\}}$$

$${\Vert L^{-1}-D\Vert }_2\le 2{\Vert D\Vert }_2{\Vert L-D^{-1}\Vert }_2\le 2\frac{4}{p^2{t}^4}\frac{t^4{p}^2ϵ}{8}=ϵ$$

$$\begin{split}P\left(\frac{1}{d_{i}}>N^{5}\right)&=P(L_{ii}>N^{4})=P\bigg{(}\sum_{j}^{N}{s_{i}^{(j)}}^{2}>N^{5}\bigg{)}\\ &\leq NP(S_{i}^{2}>N^{4})\leq NP(\lvert S_{i}\rvert>N^{2})\\ &\leq N\frac{\operatorname{\mathbb{E}}\lvert S_{i}\rvert}{N^{2}}=\frac{1}{N}\end{split}$$

$$\frac{\Vert {\tilde{C}}^{-1}-C^{-1}\Vert }{\Vert C^{-1}\Vert }\le \frac{\Vert C^{-1}E\Vert }{1-\Vert C^{-1}E\Vert }.$$

$${\Vert {\tilde{C}}^{-1}-C^{-1}\Vert }_2\le 2{\Vert C^{-1}\Vert }_2^2{\Vert E\Vert }_2.$$

$$h_{\Gamma Y}(a)=E(|a^{T}Y|)=(1/N)\sum _{i=1}^N|a^T{X}_i|.$$

$$N\ge {\left(\frac{16M{n}^4}{({ϵ}^{\prime }{)}^2}{\delta }^{\prime }\right)}^{\frac{1}{2}+\frac{3}{\gamma }}.$$

$$\Gamma X=\{\operatorname{\mathbb{E}}\bigl{(}\lambda(X)X\bigr{)}\mathrel{:}\text{$\lambda:\mathbb{R}^{n}\to[-1,1]$ is measurable}\}.$$

$$h_K(\theta )$$

$$=\underset{\lambda }{\sup }E(\lambda (X)⟨X,\theta ⟩)$$

$$=E({\lambda }^{\ast }(X)⟨X,\theta ⟩)$$

$$=E(|⟨X,\theta ⟩|).$$

$$\frac{1}{N}\sum _{i=1}^N{\lambda }_i{x}^{(i)}=q$$

$$-1\le {\lambda }_i\le 1\forall i.$$

$${\lambda }^{\ast }=\max \lambda$$

$$\frac{1}{N}\sum _{i=1}^N{\lambda }_i{x}^{(i)}=\lambda q$$

$$-1\le {\lambda }_i\le 1\forall i$$

$$h_K(y)$$

$$=\underset{{\lambda }_i\in \{-1,1\}}{\max }\frac{1}{N}\sum _{i=1}^N{\lambda }_i{x}^{(i)}\cdot y$$

$$=\frac{1}{N}\sum _{i=1}^N|x^{(i)}\cdot y|$$

$$=E(|X\cdot y|).$$

$$N\ge \max \left\{{\left(\frac{8Mn{s}_M^{1+\gamma }}{{ϵ}^2\delta }\right)}^{3/\gamma },{\left(\frac{8Mn{s}_M^{1+\gamma }}{ϵ}\right)}^{\frac{1}{2}+\frac{1}{\gamma }}\right\},$$

$$P\left(|\sum _{i=1}^N|a^T{x}^{(i)}|-E(|a^{T}X|)|>ϵ\right)\le \delta .$$

$$\Vert z-y\Vert \le {ϵ}_1.$$

$$\Vert w-q\Vert \le \Vert z-y\Vert \le {ϵ}_1.$$

$$r:=\left(1-\frac{1}{\alpha }\right)\frac{{\sigma }_{\min }(A)}{2\sqrt{n}}$$

$$N\ge {\left(\frac{8M{n}^2}{r^2\delta }\right)}^{\frac{1}{2}+\frac{1}{\gamma }}.$$

$$r$$

Full text

(1) A practical algorithm for heavy-tailed ICA that we call HTICA. We provide theoretical guarantees and show that it outperforms other algorithms in some heavy-tailed regimes, both on real and synthetic data. Like the current state-of-the-art, the new algorithm is based on the centroid body (a first moment analogue of the covariance matrix). Unlike the state-of-the-art, our algorithm is practically efficient. To achieve this, we use explicit analytic representations of the centroid body, which bypasses the use of the ellipsoid method and random walks.

(2) We study how heavy tails affect different ICA algorithms, including HTICA. Somewhat surprisingly, we show that some algorithms that use the covariance matrix or higher moments can successfully solve a range of ICA instances with infinite second moment. We study this theoretically and experimentally, with both synthetic and real-world heavy-tailed data.

1 Introduction

Independent component analysis (ICA) is a computational and statistical technique with applications in areas ranging from signal processing to machine learning and more. Formally, if $S$ is an $n$-dimensional random vector with independent coordinates and $A\in\mathbb{R}^{n\times n}$ is invertible, then the ICA problem is to estimate $A$ given access to i.i.d. samples of the mixed signals $X=AS$. We say that $X$ is generated by an ICA model $X=AS$. The recovery of $A$ (the mixing matrix) is possible only up to scaling and permutation of the columns. Moreover, for the recovery to be possible, the distributions of the random variables $S_{i}$ must not be Gaussian (except possibly one of them). Since its inception in the eighties (see [CJ10] for historical remarks), ICA has been thoroughly studied and a vast literature exists (e.g. [HKO01, CJ10]). The theory is well-developed and practical algorithms—e.g., FastICA [Hyv99], JADE [CS93]—are now available along with implementations, e.g. [CAS*+*]. However, to our knowledge, rigorous complexity analyses of these assume that the fourth moment of each component is finite: $\mathbb{E}(S_{i}^{4})<\infty$. If at least one of the independent components does not satisfy this assumption we will say that the input is in the heavy-tailed regime. Many ICA algorithms first preprocess the data to convert the given ICA model into another one where the mixing matrix $A$ has orthogonal columns; this step is often called whitening. We will instead call it orthogonalization, as this describes more precisely the desired outcome. Traditional whitening is a second order method that may not make sense in the heavy-tailed regime. In this regime, it is not clear how the existing algorithms would perform, because they depend on empirical estimation of various statistics of the data such as the covariance matrix or the fourth cumulant tensor, which diverge in general for heavy-tailed data. For example, for the covariance matrix in the mean-0 case this is done by taking the empirical average $(1/N)\sum_{i=1}^{N}x(i)x(i)^{T}$ where the $\{x(i)\}$ are i.i.d. samples of $X$. ICA in the heavy-tailed regime is of considerable interest, directly (e.g., [Kid01b, Kid01a, SYM01, CB04, CB05, SAML*+*05, WKZ09, JEK01, CS07, BC99]) and indirectly (e.g., [BG10, GTG09, WOH02]) and has applications in speech and finance. We also mention an informal connection with robust statistics: Algorithms solving heavy-tailed ICA might work by focusing on samples in a small (but high probability) region to get reliable statistics about the data and avoid the instability of the tail. Thus, if the data has outliers, the outliers are less likely to affect such an algorithm.

Recent theoretical work [AGNR15] proposed a polynomial time algorithm for ICA that works in the regime where each component $S_{i}$ has finite $(1+\gamma)$-moment for $\gamma>0$. This algorithm follows the two phases of several ICA algorithms: (i) Orthogonalize the independent components. The purpose of this step is to apply an affine transformation to the samples from $X$ so that the resulting samples correspond to an ICA model where the unknown matrix $A$ has orthogonal columns. (ii) Learn the matrix with orthogonal columns. Each of these two phases required new techniques: (1) Orthogonalization via uniform distribution in the centroid body. The input is assumed to be samples from an ICA model $X=AS$ where each $S_{i}$ is symmetrically distributed (w.l.o.g, see Sec. 2) and has at least $(1+\gamma)$-moments. The goal is to construct an orthogonalization matrix $B$ so that $BA$ has orthogonal columns. In [AGNR15], the inverse of the square root of the covariance matrix of the uniform distribution in the centroid body is one such matrix. (2) Gaussian damping. The previous step allows one to assume that the mixing matrix $A$ is orthogonal. The modified second step is: If $X$ has density $\rho_{X}(t)$ for $t\in\mathbb{R}^{n}$, then the algorithm constructs another ICA model $X_{R}=AS_{R}$ where $X_{R}$ has pdf proportional to $\rho_{X}(t)\exp({-{\lVert t\rVert}_{2}^{2}/R^{2}})$, where $R>0$ is a parameter chosen by the algorithm. This explains the term Gaussian damping. This achieves two goals: (1) All moments of $X_{R}$ and $S_{R}$ are finite. (2) The product structure of is retained. This follows from two facts: $A$ has orthogonal columns, and the Gaussian has independent components in any orthonormal basis. Because of these properties, the model can be solved by traditional ICA algorithms.

The algorithm in [AGNR15] is theoretically efficient but impractical. Their orthogonalization uses the ellipsoid algorithm for linear programming, which is not practical. It is not clear how to replace their use of the ellipsoid algorithm by practical linear programming tools, as their algorithm only has oracle access to a sort of dual and not an explicit linear program. Moreover, their orthogonalization technique uses samples uniformly distributed in the centroid body, generated by a random walk. This is computationally efficient in theory but, to the best of our knowledge, only efficient in practice for moderately low dimension.

Our contributions. Our contributions are experimental and theoretical. We provide a new and practical ICA algorithm, HTICA, building upon the previous theoretical work in [AGNR15]. HTICA works as follows: (1) Compute an orthogonalization matrix $B$. (2) Pre-multiply samples by $B$ to get an orthogonal model. (3) Damp the data, run an existing ICA algorithm. For step (1), we propose two theoretically sound and practically efficient ways below, orthogonalization via centroid body scaling and orthogonalization via covariance. Our algorithm is simpler and more efficient, but needs a more technical justification than the method in [AGNR15]. We demonstrate the effectiveness of HTICA on both synthetic and real-world data.

Orthogonalization via centroid body scaling. We propose a more practical orthogonalization matrix than the one from [AGNR15] (orthogonalization via the uniform distribution in the centroid body, mentioned before). First, consider the centroid body of random vector $X$, denoted $\Gamma X$ (this is really a function of the distribution of $X$; formal definition in Sec. 2). For intuition, it is helpful to think of the centroid body as an ellipsoid whose axes are aligned with the independent components of $X$. The centroid body is in general not an ellipsoid, but it has certain symmetries aligned with the independent components. Let random vector $Y$ be a scaling of $X$ along every ray so that points at infinity are mapped to the boundary of $\Gamma X$, the origin is mapped to itself and the scaling interpolates smoothly. One such scaling is obtained in the following way: It is helpful to consider how far a point is in its ray with respect to the boundary of $\Gamma X$. This is given by the Minkoswki functional of $\Gamma X$, denoted $p:\mathbb{R}^{n}\to\mathbb{R}$, which maps the boundary of $\Gamma X$ to 1 and interpolates linearly along every ray. We can then achieve the desired scaling by first mapping a given point to the boundary point on its ray (the mapping $x\mapsto x/p(x)$) and then using the function $\tanh$, which maps $[0,\infty)$ to $[0,1]$ with $\tanh(0)=0$ and $\lim_{x\to\infty}\tanh(x)=1$ to determine the final scale along the ray, namely, $\tanh p(x)$. More formally, our scaling is the following: Let $Y$ be $\frac{\tanh p(X)}{p(X)}X$. We show in Sec. 4.1 that $B=\operatorname{Cov}(Y)^{-1/2}$ is an orthogonalization matrix when $\operatorname{Cov}(Y)$ is invertible. In order to make this practical, one needs a practical estimator of the Minkowski functional of $\Gamma X$ from a sample of $X$. In Sec. 4.1 and 5, we describe such an algorithm and provide a theoretical justification, including finite sample estimates. The proposed algorithm is much simpler and practical than the one described in [AGNR15]. In particular, it avoids the use of the ellipsoid algorithm by the use of a closed-form linear programming representation of the centroid body (Prop. 10, Lemma 11) and new approximation guarantees between the empirical (sample estimate) and true centroid body of a heavy-tailed distribution. In Sec. 4.1, we discuss our practical implementation and show results where orthogonalization via centroid body scaling produces results with smaller error.

Orthogonalization via covariance. Previously, (e.g., in [CB04]), the empirical covariance matrix was used for whitening in the heavy-tailed regime and, surprisingly, worked well in some situations. Unfortunately, the understanding of this was quite limited . We give a theoretical explanation for this phenomenon in a fairly general heavy-tailed regime: Covariance-based orthogonalization works well when each component $S_{i}$ has finite $(1+\gamma)$-moment, where $\gamma>0$. We also study this algorithm in experimental settings. As we will see, while orthogonalization via covariance improves over previous algorithms, in general orthogonalization via centroid body has better performance because it has better numerical stability; but there are some situations where orthogonalization via covariance matrix is better.

Empirical Study. We perform experiments on both synthetic and real data to see the effect of heavy-tails on ICA.

In the synthetic data setting, we generate samples from a fixed heavy-tailed distribution and study how well the algorithm can recover a random mixing matrix (Sec. 3).

To study the algorithm with real data, we use recordings of human speech provided by [Don09]. This involves a room with different arrangements of microphones, and six humans speaking independently. The speakers are recorded individually, so we can artificially mix them and have access to a ground truth. We study the statistical properties of the data, observing that it does indeed behave as if the underlying processes are heavy-tailed. The performance of our algorithm shows improvement over using FastICA on its own.

2 Preliminaries

Heavy-tailed distributions arise in a wide range of applications (e.g., [Nol15]). They are characterized by the slow decay of their tails. Examples of heavy-tailed distributions include the Pareto and log-normal distributions.

We denote the pdf of random variable $Z$ by $\rho_{Z}$. We will assume that our distributions are symmetric, that is $\rho(x)=\rho(-x)$ for $x\in\mathbb{R}$. As observed in [AGNR15], this is without loss of generality for our purposes. This follows from the fact that if $X=AS$ is an ICA model, and if we let $X^{\prime}=AS^{\prime}$ be an i.i.d. copy of the same model, then $X-X^{\prime}=A(S-S^{\prime})$ is an ICA model with components of $S-S^{\prime}$ having symmetric pdfs. One further needs to check that if the components of $S$ are away from Gaussians then the same holds for $S-S^{\prime}$; see [AGNR15]. We formulate our algorithms for the symmetric case; the general case immediately reduces to the symmetric case.

For $K\subseteq\mathbb{R}^{n}$, $K_{\epsilon}$ denotes the set of points that are at distance at most $\epsilon$ from $K$. The set $K_{-\epsilon}$ is all points for which an $\epsilon$-ball around them is still contained in $K$. The ${n}$-dimensional $\ell_{p}$ ball is denoted as $B_{p}^{n}$.

An important related family of distributions is that of stable distributions (e.g., [Nol15]). In general, the density of a stable distribution has no closed form, but is fully defined by four real-valued parameters. Some stable distributions do admit a closed form, such as the Cauchy and Gaussian distributions. For us the most important parameter is $\alpha\in(0,2]$, known as the stability parameter; we will think of the other three parameters as being fixed to constants.

We use the notation $\operatorname{poly}(\cdot)$ to indicate a function which is asymptotically upper bounded by a polynomial expression of the given variables.

If $\alpha=2$, the distribution is Gaussian (the only non-heavy-tailed stable distribution), and if $\alpha=1$, it is the Cauchy distribution.

Definition 1** (Centroid body).**

Let $X\in\mathbb{R}^{n}$ be a random vector with finite first moment, that is, for all $u\in\mathbb{R}^{n}$ we have $\operatorname{\mathbb{E}}(\lvert\langle{u},{X}\rangle\rvert)<\infty$. The centroid body of $X$ is the compact convex set, denoted $\Gamma X$, whose support function is $h_{\Gamma X}(u)=\operatorname{\mathbb{E}}(\lvert\langle{u},{X}\rangle\rvert)$. For a probability measure $\mathbb{P}$, we define $\Gamma\mathbb{P}$, the centroid body of $\mathbb{P}$, as the centroid body of any random vector distributed according to $\mathbb{P}$.

Note that for the centroid body to be well-defined, the mean of the data must be finite. This excludes, for instance, the Cauchy distribution from consideration in the present work.

3 HTICA and experiments

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 $f_{\eta}(x)$ proportional to ${(\lvert x\rvert+1.5)^{-\eta}}$, which is symmetric, and for $\eta>1$, $f_{\eta}$ is the density of a distribution which has finite $k<\eta-1$ moment. The signal $S$ is generated with each $S_{i}$ independently distributed from $f_{\eta_{i}}$. The mixing matrix $A\in\mathbb{R}^{{n}\times{n}}$ is generated with each coordinate i.i.d. $\mathcal{N}(0,1)$, columns normalized to unit length. To compare the quality of recovery, the columns of the estimated mixing matrix, $\tilde{A}$ are permuted to align with the closest matching column of $A$, 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 [ACY*+*96]; the results have similar behavior and are not presented here.

3.1 Heavy-tailed ICA when $A$ 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 [AGNR15], Gaussian damping is a preprocessing technique that converts data from an ICA model $X=AS$, where $A$ is unitary (columns are orthogonal with unit $l_{2}$-norm) to data from a related ICA model $X_{R}=AS_{R}$, where $R>0$ is a parameter to be chosen. The independent components of $S_{R}$ have finite moments of all orders and so the existing algorithms can estimate $A$.

Using samples of $X$, we construct the damped random variable $X_{R}$, with pdf $\rho_{X_{R}}(x)\propto\rho_{X}(x)\exp({-{\lVert x\rVert}^{2}/R^{2}})$. To normalize the right hand side, we can estimate

[TABLE]

so that

[TABLE]

If $x$ is a realization of $X_{R}$, then $s=A^{-1}x$ is a realization of the random variable $S_{R}$ and we have that $S_{R}$ has pdf $\rho_{S_{R}}(s)=\rho_{X_{R}}(x)$. To generate samples from this distribution, we use rejection sampling on samples from $\rho_{X}$. When performing the damping, we binary search over $R$ so that about 25% of the samples are rejected. For more details about the technical requirements for choosing $R$, see [AGNR15].

Figure 1 shows that, when $A$ is already a perfectly orthogonal matrix, but where $S$ 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 $\log\cosh$ (“FastICA - tanh”), (3) JADE, and (4) Second Order Joint Diagonalization as in, e.g., [Car89] .

3.2 Experiments on synthetic heavy-tailed data

We now present the results of HTICA using different orthogonalization techniques: (1) Orthogonalization via covariance (Section 4.2 (2) Orthogonalization via the centroid body (Section 4.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 (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. A partial explanation is that, given the results in Figure 1, we know that equal exponents favor FastICA without damping and orthogonalization (identity in Figure 2). 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 $R$ is sensitive to scaling. A scaled-up distribution will be somewhat hurt because fewer samples will survive damping.

3.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. [Kid00]) that speech data can be modeled by $\alpha$-stable distributions. For an $\alpha$-stable distribution, with $\alpha\in(0,2)$, only the moments of order less than $\alpha$ will be finite. We present here some results on a data set of human speech according to the standard cocktail party model, from [Don09].

The physical setup of the experiments (the human speakers and microphones) is shown in Figure 3.

To estimate whether the data is heavy-tailed, as in [Kid00], we estimate parameter $\alpha$ of a best-fit $\alpha$-stable distribution. This estimate is in Figure 4 for one of the data sets collected. We can see that the estimated $\alpha$ is clearly in the heavy-tailed regime for some signals.

Using data from [Don09], we perform the same experiment as in Section 3.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 4 demonstrates that HTICA (orthogonalizing with centroid body scaling, Section 4.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.

4 New approach to orthogonalization and a new analysis of empirical covariance

As noted above, the technique in [AGNR15], 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 [AGNR15]). 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 3 that the covariance of this subset of samples is enough to orthogonalize the mixing matrix $A$. 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 $A$. 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 4 containing the singular value and condition number of the mixing matrix $BA$ in the approximately orthogonal ICA model.

4.1 Orthogonalization via centroid body scaling

In [AGNR15], 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 $p(x)$, the Minkowski functional of the centroid body. The complete procedure is stated in Subroutine 1.

We now explain how to estimate the Minkowski functional. The Minkowski functional was informally described in Section 1. The Minkowski functional of $\Gamma X$ is formally defined by $p(x):=\inf\{t>0\mathrel{:}x\in t\Gamma X\}$. Our estimation of $p(x)$ is based on an explicit linear program (LP) (10) that gives the Minkowski functional of the centroid body of a finite sample of $X$ exactly and then arguing that a sample estimate is close to the actual value for $\Gamma X$. For clarity of exposition, we only analyze formally a special case of LP (10) that decides membership in the centroid body of a finite sample of $X$ (LP (9)) and approximate membership in $\Gamma X$. This analysis is in Section 5. Accuracy guarantees for the approximation of the Minkowski functional follow from this analysis.

Lemma 2** ([AGNR15]).**

Let $U$ be a family of $n$-dimensional product distributions. Let $\bar{U}$ be the closure of $U$ under invertible linear transformations. Let $Q(\mathbb{P})$ be an $n$-dimensional distribution defined as a function of $\mathbb{P}\in\bar{U}$. Assume that $U$ and $Q$ satisfy:

For all $\mathbb{P}\in U$, $Q(\mathbb{P})$ is absolutely symmetric. 2. 2.

$Q$* is linear equivariant (that is, for any invertible linear transformation $T$ we have $Q(T\mathbb{P})=TQ(\mathbb{P})$).* 3. 3.

For any $\mathbb{P}\in\bar{U}$, $\operatorname{Cov}(Q(\mathbb{P}))$ is positive definite.

Then for any symmetric ICA model $X=AS$ with $\mathbb{P}_{S}\in U$ we have $\operatorname{Cov}(Q(\mathbb{P}_{X}))^{-1/2}$ is an orthogonalizer of $X$.

Lemma 3**.**

*Let $X$ be a random vector drawn from an ICA model $X=AS$ such that for all $i$ we have $\operatorname{\mathbb{E}}\lvert S_{i}\rvert=1$ and $S_{i}$ is symmetrically distributed. Let $Y=\frac{\tanh p(X)}{p(X)}X$ where $p(X)$ is the Minkoswki functional of $\Gamma X$.

Then $\operatorname{Cov}(Y)^{-1/2}$ is an orthogonalizer of $X$.*

Proof.

We will be applying Lemma 2. Let $U$ denote the set of absolutely symmetric product distributions $\mathbb{P}_{W}$ over $\mathbb{R}^{n}$ such that $\operatorname{\mathbb{E}}\lvert W_{i}\rvert=1$ for all $i$. For $\mathbb{P}_{V}\in\bar{U}$, let $Q(\mathbb{P}_{V})$ be equal to the distribution obtained by scaling $V$ as described earlier, that is, distribution of $\alpha V$, where $\alpha=\frac{\tanh p(V)}{p(V)}$, $p(V)$ is the Minkoswki functional of $\Gamma\mathbb{P}_{V}$.

For all $\mathbb{P}_{W}\in U$, $W_{i}$ is symmetric and $\operatorname{\mathbb{E}}\lvert W_{i}\rvert=1$ which implies that $\alpha W$, that is, $Q(\mathbb{P}_{W})$ is absolutely symmetric. Let $\mathbb{P}_{V}\in\bar{U}$. Then $Q(\mathbb{P}_{V})$ is equal to the distribution of $\alpha V$. For any invertible linear transformation $T$ and measurable set $\mathcal{M}$, we have $Q(T\mathbb{P}_{V})(\mathcal{M})=Q(\mathbb{P}_{TV})(\mathcal{M})=\mathbb{P}_{\alpha TV}(\mathcal{M})=\mathbb{P}_{\alpha V}(T^{-1}\mathcal{M})=TQ(\mathbb{P}_{V})(\mathcal{M})$.

Thus $Q$ is linear equivariant. Let $\mathbb{P}\in\bar{U}$. Then there exist $A$ and $\mathbb{P}_{W}\in U$ such that $\mathbb{P}=A\mathbb{P}_{W}$. We get $\operatorname{Cov}(Q(\mathbb{P}))=\operatorname{Cov}(AQ(\mathbb{P}_{W}))$. Let $W_{\alpha}=\alpha W$. Thus, $\operatorname{Cov}(AQ(\mathbb{P}_{W}))=A\operatorname{\mathbb{E}}(W_{\alpha}{W}_{\alpha}^{T})A^{T}$ where $\operatorname{\mathbb{E}}(W_{\alpha}{W}_{\alpha}^{T})$ is a diagonal matrix with elements $\operatorname{\mathbb{E}}(\alpha^{2}W_{i}^{2})$ which are non-zero because we assume $\operatorname{\mathbb{E}}\lvert W_{i}\rvert=1$. This implies that $\operatorname{Cov}(Q(\mathbb{P}))$ is positive definite and thus by Lemma 2, $\operatorname{Cov}(Y)^{-1/2}$ is an orthogonalizer of $X$. ∎

4.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 [HKO01]: Given samples from the ICA model $X=AS$, we compute the empirical covariance matrix $\tilde{\Sigma}$ which tends to the true covariance matrix as we take more samples and set $B=\tilde{\Sigma}^{-1/2}$. Then one can show that $BA$ is a rotation matrix, and thus by pre-multiplying the data by $B$ we obtain an ICA model $Y=BX=(BA)S$, where the mixing matrix $BA$ 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., [CB04]), the empirical covariance matrix was used for whitening in the heavy-tailed regime with good empirical performance; [CB04] 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 $B:=\tilde{\Sigma}^{-1/2}$) “works” under considerably more general conditions, namely if $(1+\gamma)$-moment exists for $\gamma>0$ for each independent component $S_{i}$. By “works” we mean that instead of whitening the data (that is $BA$ 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 $BA$ is close to a matrix whose columns are orthogonal. In other words, $(BA)^{T}(BA)$ is close to a diagonal matrix (in a sense made precise in Theorem 5).

Let $X$ be a real-valued symmetric random variable such that $\mathbb{E}(\lvert X\rvert^{1+\gamma})\leq M$ for some $M>1$ and $0<\gamma<1$. The following lemma from [AGNR15] says that the empirical average of the absolute value of $X$ converges to the expectation of $|X|$. The proof, which we omit, follows an argument similar to the proof of the Chebyshev’s inequality. Let $\tilde{\mathbb{E}}_{N}[\lvert X\rvert]$ be the empirical average obtained from $N$ independent samples $X^{(1)},\ldots,X^{(N)}$, i.e., $(\lvert X^{(1)}\rvert+\dotsb+\lvert X^{(N)}\rvert)/N$.

Lemma 4**.**

Let $\epsilon\in(0,1)$. With the notation above, for $N\geq\left(\frac{8M}{\epsilon}\right)^{\frac{1}{2}+\frac{1}{\gamma}}$, we have $\Pr[\lvert\tilde{\mathbb{E}}_{N}[\lvert X\rvert]-\mathbb{E}[\lvert X\rvert]\rvert>\epsilon]\leq\frac{8M}{\epsilon^{2}N^{\gamma/3}}$.

Theorem 5** (Orthogonalization via covariance matrix).**

Let $X$ be given by ICA model $X=AS$. Assume that there exist $t,p,M>0$ and $\gamma\in(0,1)$ such that for all $i$ we have

(a) $\operatorname{\mathbb{E}}(\lvert S_{i}\rvert^{1+\gamma})\leq M<\infty$,

(b) (normalization) $\operatorname{\mathbb{E}}\lvert S_{i}\rvert=1$, and

*(c) $\Pr(\lvert S_{i}\rvert\geq t)\geq p$. Let $x^{(1)},\dotsc,x^{(N)}$ be i.i.d. samples according to $X$. Let $\tilde{\Sigma}=(1/N)\sum_{k=1}^{N}x^{(k)}{x^{(k)}}^{T}$ and $B=\tilde{\Sigma}^{-1/2}$. Then for any $\epsilon,\delta\in(0,1)$, ${\lVert(BA)^{T}BA-D\rVert}_{2}\leq\epsilon$ for a diagonal matrix $D$ with diagonal entries $d_{1},\dotsc,d_{n}$ satisfying $0<d_{i},1/d_{i}\leq\max\{2/pt^{2},N^{4}\}$ for all $i$ with probability $1-\delta$ when $N\geq\operatorname{poly}(n,M,1/p,1/t,1/\epsilon,1/\delta)$. *

Proof idea. For $i\neq j$ we have $\operatorname{\mathbb{E}}(S_{i}S_{j})=0$ (due to our symmetry assumption on $S$) and $\operatorname{\mathbb{E}}(\lvert S_{i}S_{j}\rvert)=\operatorname{\mathbb{E}}(\lvert S_{i}\rvert)\operatorname{\mathbb{E}}(\lvert S_{j}\rvert)<\infty$. We have $(BA)^{T}BA=L^{-1}$, where $L=(1/N)\sum_{k=1}^{N}s^{(k)}{s^{(k)}}^{T}$. The off-diagonal entries of $L$ converge to [math]: We have $L_{i,j}=\mathbb{E}S_{i}S_{j}=(\mathbb{E}S_{i})(\mathbb{E}S_{j})$. Now by our assumption that $(1+\gamma)$-moments exist, Lemma 4 is applicable and implies that empirical average $\tilde{\mathbb{E}}S_{i}$ tends to the true average $\mathbb{E}S_{i}$ 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 $X_{i}$ and hence $S_{i}$ have been centered). The diagonal entries of $L$ 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 $L$ the diagonal entries highly dominate the off-diagonal entries, then the same is true of $L^{-1}$.

$\square$

Proof.

We have $(BA)^{T}BA=L^{-1}$, where $L=(1/N)\sum_{k=1}^{N}s^{(k)}{s^{(k)}}^{T}$. By assumption, $\operatorname{\mathbb{E}}L_{ij}=0$ for $i\neq j$. Note that $\operatorname{\mathbb{E}}\lvert s_{i}s_{j}\rvert^{1+\gamma}\leq M^{2}$ and so by Lemma 4, for $i\neq j$,

[TABLE]

when $N\geq(\frac{8M^{2}}{\epsilon_{1}})^{\frac{1}{2}+\frac{1}{\gamma}}$.

Now let $D:=\operatorname{diag}(L^{-1}_{11},L^{-1}_{22},\dots,L^{-1}_{{n}{n}})$. Then when $\lvert L_{ij}\rvert<\epsilon_{1}$ for all $i\neq j$, we have ${\lVert L-D^{-1}\rVert}_{2}\leq{\lVert L-D^{-1}\rVert}_{F}\leq n\epsilon_{1}$. The union bound then implies

[TABLE]

when $N\geq(\frac{8M^{2}}{\epsilon_{1}})^{\frac{1}{2}+\frac{1}{\gamma}}$.

Next, we aim to bound ${\lVert D\rVert}_{2}$ which can be done by writing

[TABLE]

where $L_{ii}=(1/N)\sum_{k=1}^{N}{s_{i}^{(k)}}^{2}$. Consider the random variable $\mathbbm{1}(s_{i}^{2}\geq t^{2})$. We can calculate $\operatorname{\mathbb{E}}\sum_{j}\mathbbm{1}({s_{i}^{(j)}}^{2}\geq t^{2})\geq Np$ and use a Chernoff bound to see

[TABLE]

and when $\sum_{k\in[N]}\mathbbm{1}({s_{i}^{(k)}}^{2}\geq t^{2})\geq\frac{Np}{2}$, we have $L_{ii}\geq t^{2}p/2$. Then with probability at least $1-n\exp(-Np/8)$, all entries of $D^{-1}$ are at least $t^{2}p/2$. Using this, if $N\geq N_{1}:=(8/p)\ln(3n/\delta)$ then ${\lVert D\rVert}_{2}\leq 2/pt^{2}$ with probability at least $1-\delta/3$.

Similarly, suppose that ${\lVert D\rVert}_{2}\leq 2/pt^{2}$ and choose $\epsilon_{1}=\min\{\frac{t^{4}p^{2}}{4n}\cdot\frac{\epsilon}{2},\frac{1}{pt^{2}}\}$ and

[TABLE]

so that when $N\geq N_{2}$, we have ${\lVert L-D^{-1}\rVert}_{2}\leq 1/(2{\lVert D\rVert}_{2})$ and ${\lVert L-D^{-1}\rVert}_{2}\leq t^{4}p^{2}\epsilon/8$ with probability at least $1-\delta/3$. Invoking (7), when $N\geq\max\{N_{1},N_{2}\}$, we have

[TABLE]

with probability at least $1-2\delta/3$.

Finally, we upper bound $1/d_{i}$ for a fixed $i$ by using Markov’s inequality:

[TABLE]

so that $1/d_{i}\leq N^{4}$ for all $i$ with probability at least $1-\delta/3$ when $N\geq N_{3}:=n/3\delta$. Therefore, when $N\geq\max\{N_{1},N_{2},N_{3}\}$, we have ${\lVert L^{-1}-D\rVert}_{2}\leq\epsilon$, $d_{i}\leq 2/pt^{2}$, and $1/d_{i}\leq N^{4}$ for all $i$ with overall probability at least $1-\delta$. ∎

We used the following technical result.

Lemma 6**.**

Let ${\lVert\cdot\rVert}$ be a matrix norm such that ${\lVert AB\rVert}\leq{\lVert A\rVert}{\lVert B\rVert}$. Let matrices $C,E\in\mathbb{R}^{n\times n}$ be such that ${\lVert C^{-1}E\rVert}_{2}\leq 1$, and let $\tilde{C}=C+E$. Then

[TABLE]

This implies that if ${\lVert E\rVert}_{2}={\lVert\tilde{C}-C\rVert}_{2}\leq 1/(2{\lVert C^{-1}\rVert}_{2})$, then

[TABLE]

In Theorem 5, 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 4, where the condition number of “centroid” is much smaller than the condition number of “covariance.”

5 Membership oracle for the centroid body, without polarity

We will see now how to implement an $\epsilon$-weak membership oracle for $\Gamma X$ directly, without using polarity. We start with an informal description of the algorithm and its correctness.

The algorithm implementing the oracle (Subroutine 2) is the following: Let $q\in\mathbb{R}^{n}$ be a query point. Let $X_{1},\dotsc,X_{N}$ be a sample of random vector $X$. Given the sample, let $Y$ be uniformly distributed in $\{X_{1},\dotsc,X_{N}\}$. Output YES if $q\in\Gamma Y$, else output NO.

Idea of the correctness of the algorithm: If $q$ is not in $(\Gamma X)_{\epsilon}$, then there is a hyperplane separating $q$ from $(\Gamma X)_{\epsilon}$. Let $\{x\mathrel{:}a^{T}x=b\}$ be the hyperplane, satisfying ${\lVert a\rVert}=1$, $a^{T}q>b$ and $a^{T}x\leq b$ for every $x\in(\Gamma X)_{\epsilon}$. Thus, we have $h_{(\Gamma X)_{\epsilon}}(a)\leq b$ and $h_{\Gamma X}(a)\leq b-\epsilon$. We have

[TABLE]

By Lemma 4, $(1/N)\sum_{i=1}^{N}\lvert a^{T}X_{i}\rvert$ is within $\epsilon$ of $\operatorname{\mathbb{E}}\lvert a^{T}X\rvert=h_{\Gamma X}(a)\leq b-\epsilon$ when $N$ is large enough with probability at least $1-\delta$ over the sample $X_{1},\dotsc,X_{N}$. In particular, $h_{\Gamma Y}(a)\leq b$, which implies $q\notin\Gamma Y$ and the algorithm outputs NO, with probability at least $1-\delta$.

If $q$ is in $(\Gamma X)_{-\epsilon}$, let $y=q+\epsilon\hat{q}\in\Gamma X$. We will prove the following claim:

Informal claim (Lemma 13): For $p\in\Gamma X$, for large enough $N$ and with probability at least $1-\delta$ there is $z\in\Gamma Y$ so that ${\lVert z-p\rVert}\leq\epsilon/10$.

This claim applied to $p=y$ to get $z$, convexity of $\Gamma Y$ and the fact that $\Gamma Y$ contains $B\simeq\sigma_{\min}(A)B_{2}^{n}$ (Lemma 9) imply that $q\in\operatorname{conv}(B\cup\{z\})\subseteq\Gamma Y$ and the algorithm outputs YES.

We will prove the claim now. Let $p\in\Gamma X$. By the dual characterization of the centroid body (Proposition 10), there exists a function $\lambda:\mathbb{R}^{n}\to\mathbb{R}$ such that $p=\operatorname{\mathbb{E}}(\lambda(X)X)$ with $-1\leq\lambda\leq 1$. Let $z=\frac{1}{N}\sum_{i=1}^{N}\lambda(X_{i})X_{i}.$ We have $\operatorname{\mathbb{E}}_{X_{i}}(\lambda(X_{i})X_{i})=p$ and $\operatorname{\mathbb{E}}_{X_{i}}(\lvert\lambda(X_{i})X_{i}\rvert^{1+\gamma})\leq\operatorname{\mathbb{E}}_{X_{i}}(\lvert X_{i}\rvert^{1+\gamma})\leq M$. By Lemma 4 and a union bound over every coordinate we get $\mathbb{P}({\lVert p-z\rVert}\geq\epsilon)\leq\delta$ for $N$ large enough.

5.1 Formal Argument

Lemma 7** ([AGNR15]).**

Let $S=(S_{1},\dots,S_{n})\in\mathbb{R}^{n}$ be an absolutely symmetrically distributed random vector such that $\operatorname{\mathbb{E}}(\lvert S_{i}\rvert)=1$ for all $i$. Then $B_{1}^{{n}}\subseteq\Gamma S\subseteq[-1,1]^{n}$. Moreover, ${n}^{-1/2}B_{2}^{n}\subseteq(\Gamma S)^{\circ}\subseteq\sqrt{{n}}B_{2}^{n}$.

Lemma 8** ([AGNR15]).**

Let $X$ be a random vector on $\mathbb{R}^{n}$. Let $A:\mathbb{R}^{n}\to\mathbb{R}^{n}$ be an invertible linear transformation. Then $\Gamma(AX)=A(\Gamma X)$.

Lemma 9**.**

Let $S=(S_{1},\dots,S_{n})\in\mathbb{R}^{n}$ be an absolutely symmetrically distributed random vector such that $\operatorname{\mathbb{E}}(\lvert S_{i}\rvert)=1$ and $\operatorname{\mathbb{E}}(\lvert S_{i}\rvert^{1+\gamma})\leq M<\infty$ for all $i$. Let $S^{(i)},i=1,\dotsc,N$ be a sample of i.i.d. copies of $S$. Let $Y$ be a random vector, uniformly distributed in $S^{(1)},\dotsc,S^{(N)}$. Then $(1-\epsilon^{\prime})B_{1}^{{n}}\subseteq\Gamma Y$ whenever

[TABLE]

Proof.

From Lemma 7 we know $\pm e_{i}\in\Gamma S$. It is enough to apply Lemma 13 to $\pm e_{i}$ with $\epsilon=\epsilon^{\prime}/\sqrt{{n}}$ and $\delta=\delta^{\prime}/(2{n})$. This gives, for any $\theta\in S^{{n}-1}$, $h_{\Gamma Y}(\theta)\geq h_{\Gamma S}(\theta)-\epsilon\geq h_{B_{1}^{n}}(\theta)-\epsilon\geq(1-\sqrt{{n}}\epsilon)h_{B_{1}^{n}}(\theta)=(1-\epsilon^{\prime})h_{B_{1}^{n}}(\theta)$. In particular, $\Gamma Y\supseteq(1-\epsilon^{\prime})B_{1}^{n}$. ∎

Proposition 10** (Dual characterization of centroid body).**

Let $X$ be a ${n}$-dimensional random vector with finite first moment, that is, for all $u\in\mathbb{R}^{n}$ we have $\operatorname{\mathbb{E}}(\lvert\langle{u},{X}\rangle\rvert)<\infty$. Then

[TABLE]

Proof.

Let $K$ denote the rhs of the conclusion.We will show that $K$ is a non-empty, closed convex set and show that $h_{K}=h_{\Gamma X}$, which implies (8).

By definition, $K$ is a non-empty bounded convex set. To see that it is closed, let $(y_{k})_{k}$ be a sequence in $K$ such that $y_{k}\to y\in\mathbb{R}^{n}$. Let $\lambda_{k}$ be the function associated to $y_{k}\in K$ according to the definition of $K$. Let $\mathbb{P}_{X}$ be the distribution of $X$. We have ${\lVert\lambda_{k}\rVert}_{L^{\infty}(\mathbb{P}_{X})}\leq 1$ and, passing to a subsequence $k_{j}$, $(\lambda_{k_{j}})$ converges to $\lambda\in L^{\infty}(\mathbb{P}_{X})$ in the weak-$*$ topology $\sigma(L^{\infty}(\mathbb{P}_{X}),L^{1}(\mathbb{P}_{X}))$, where $-1\leq\lambda\leq 1$. 555This is a standard argument, see [Bre11] for the background. Map $x\mapsto x_{i}$ is in $L^{1}(\mathbb{P}_{X})$. [Bre11, Theorem 4.13] gives that $L^{1}(\mathbb{P}_{X})$ is a separable Banach space. [Bre11, Theorem 3.16] (Banach-Alaoglu-Bourbaki) gives that the unit ball in $L^{\infty}(\mathbb{P}_{X})$ is compact in the weak-* topology. [Bre11, Theorem 3.28] gives that the unit ball in $L^{\infty}(\mathbb{P}_{X})$ is metrizable and therefore sequentially compact in the weak-* topology. Therefore, any bounded sequence in $L^{\infty}(\mathbb{P}_{X})$ has a convergent subsequence in the weak-* topology. This implies $\lim_{j}\operatorname{\mathbb{E}}(\lambda_{k_{j}}(X)X_{i})=\lim_{j}\int_{\mathbb{R}^{n}}\lambda_{k_{j}}(x)x_{i}\mathop{}\!\mathrm{d}\mathbb{P}_{X}(x)=\int_{\mathbb{R}^{n}}\lambda(x)x_{i}\mathop{}\!\mathrm{d}\mathbb{P}_{X}(x)=\operatorname{\mathbb{E}}(\lambda(X)X_{i})$. Thus, we have $y=\lim_{j}y_{k_{j}}=\lim_{j}\operatorname{\mathbb{E}}((\lambda_{k_{j}}(X)X)=\operatorname{\mathbb{E}}(\lambda(X)X)$ and $K$ is closed.

To conclude, we compute $h_{K}$ and see that it is the same as the definition of $h_{\Gamma X}$. In the following equations $\lambda$ ranges over functions such that $\lambda:\mathbb{R}^{n}\to\mathbb{R}$ is Borel-measurable and $-1\leq\lambda\leq 1$.

[TABLE]

∎

Lemma 11** (LP).**

Let $X$ be a random vector uniformly distributed in $\{x^{(i)}\}_{i=1}^{N}\subseteq\mathbb{R}^{n}$. Let $q\in\mathbb{R}^{n}$. Then:

$\Gamma X=\frac{1}{N}\sum_{i=1}^{N}[-x^{(i)},x^{(i)}]$. 2. 2.

Point $q~{}\in~{}\Gamma X$ iff there is a solution $\lambda~{}\in~{}\mathbb{R}^{N}$ to the following linear feasibility problem:

[TABLE] 3. 3.

Let $\lambda^{*}$ be the optimal value of (always feasible) linear program

[TABLE]

with $\lambda^{*}=\infty$ if the linear program is unbounded. Then the Minkowski functional of $\Gamma X$ at $q$ is $1/\lambda^{*}$.

Proof.

This is proven in [McM71]. It is also a special case of Proposition 10. We include an argument here for completeness. Let $K:=\frac{1}{N}\sum_{i=1}^{N}[-x^{(i)},x^{(i)}]$. We compute $h_{K}$ to see it is the same as $h_{\Gamma X}$ in the definition of $\Gamma X$ (Definition 1). As $K$ and $\Gamma X$ are non-empty compact convex sets, this implies $K=\Gamma X$. 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 12** (Correctness of Subroutine 2).**

Let $X=AS$ be given by an ICA model such that for all $i$ we have $\operatorname{\mathbb{E}}(\lvert S_{i}\rvert^{1+\gamma})\leq M<\infty$, $S_{i}$ is symmetrically distributed and normalized so that $\operatorname{\mathbb{E}}\lvert S_{i}\rvert=1$. Then, given a query point $q\in\mathbb{R}^{n}$, $\epsilon>0$, $\delta>0$, $s_{M}\geq\sigma_{\max}(A)$, and $s_{m}\leq\sigma_{\min}(A)$, Subroutine 2 is an $\epsilon$-weak membership oracle for $q$ and $\Gamma X$ with probability $1-\delta$ using time and sample complexity $\operatorname{poly}(n,M,1/s_{m},s_{M},1/\epsilon,1/\delta).$

Proof.

Let $Y$ be uniformly random in $(x^{(i)})_{i=1}^{N}$. There are two cases corresponding to the guarantees of the oracle:

• •

Case $q\notin(\Gamma X)_{\epsilon}$. Then there is a hyperplane separating $q$ from $(\Gamma X)_{\epsilon}$. Let $\{x\in\mathbb{R}^{n}\mathrel{:}a^{T}x=b\}$ be the separating hyperplane, parameterized so that $a\in\mathbb{R}^{n}$, $b\in\mathbb{R}$, ${\lVert a\rVert}=1$, $a^{T}q>b$ and $a^{T}x\leq b$ for every $x\in(\Gamma X)_{\epsilon}$. In this case $h_{(\Gamma X)_{\epsilon}}(a)\leq b$ and $h_{\Gamma X}(a)\leq b-\epsilon$. At the same time, $h_{\Gamma Y}(a)=\operatorname{\mathbb{E}}(\lvert a^{T}Y\rvert)=(1/N)\sum_{i=1}^{N}\lvert a^{T}x^{(i)}\rvert$.

We want to apply Lemma 4 to $a^{T}X$ to get that $h_{\Gamma Y}(a)=(1/N)\sum_{i=1}^{N}\lvert a^{T}x^{(i)}\rvert$ is within $\epsilon$ of $h_{\Gamma X}(a)=\operatorname{\mathbb{E}}(\lvert a^{T}X\rvert)$. For this we need a bound on the $(1+\gamma)$-moment of $a^{T}X$. We use the bound from [AGNR15, Equation (10)]: $\operatorname{\mathbb{E}}(\lvert a^{T}X\rvert^{1+\gamma})\leq M{n}s_{M}^{1+\gamma}$. Lemma 4 implies that for

[TABLE]

we have

[TABLE]

In particular, with probability at least $1-\delta$ we have $h_{\Gamma Y}(a)\leq b$, which implies $q\notin\Gamma Y$ and, by Lemma 11, Subroutine 2 outputs NO.

• •

Case $q\in(\Gamma X)_{-\epsilon}$. Let $y=q+\epsilon\hat{q}=q(1+\frac{\epsilon}{{\lVert q\rVert}})$. Let $\alpha=1+\frac{\epsilon}{{\lVert q\rVert}}$. Then $y\in\Gamma X$. Invoke Lemma 13 for i.i.d. sample $(x^{(i)})_{i=1}^{N}$ of $X$ with $p=y$ and $\epsilon$ equal to some $\epsilon_{1}>0$ to be fixed later to conclude $y\in(\Gamma Y)_{\epsilon_{1}}$. That is, there exist $z\in\Gamma Y$ such that

[TABLE]

Let $w=z/\alpha$. Given (12) and the relationships $y=\alpha q$ and $z=\alpha w$, we have

[TABLE]

From Lemma 9 with $\epsilon^{\prime}=1/2$ and equivariance of the centroid body (Lemma 8) we get $\Gamma Y\supseteq\frac{\sigma_{\min}(A)}{2\sqrt{{n}}}B_{2}^{n}$. This and convexity of $\Gamma Y$ imply $\operatorname{conv}\{\frac{\sigma_{\min}(A)}{2\sqrt{{n}}}B_{2}^{n}\cup\{z\}\}\subseteq\Gamma Y$. In particular, the ball around $w$ of radius

[TABLE]

is contained in $\Gamma Y$. The choice $\epsilon_{1}=r\geq$ and (13) imply $q\in\Gamma Y$ and Subroutine 2 outputs YES whenever

[TABLE]

To conclude, remember that $q\in(\Gamma X)_{-\epsilon}$. Therefore ${\lVert q\rVert}+\epsilon\leq\sqrt{{n}}\sigma_{\max}(A)$ (from Lemma 7 and equivariance of the centroid body, Lemma 8). This implies

[TABLE]

The claim follows. ∎

Lemma 13**.**

Let $X$ be a ${n}$-dimensional random vector such that for all coordinates $i$ we have $\operatorname{\mathbb{E}}(\lvert X_{i}\rvert^{1+\gamma})\leq M<\infty$. Let $p\in\Gamma X$. Let $(X^{(i)})_{i=1}^{N}$ be an i.i.d. sample of $X$. Let $Y$ be uniformly random in $(X^{(i)})_{i=1}^{N}$. Let $\epsilon>0$, $\delta>0$. If $N\geq\left(\frac{8M{n}^{2}}{\epsilon^{2}\delta}\right)^{\frac{1}{2}+\frac{3}{\gamma}}$, then, with probability at least $1-\delta$, $p\in(\Gamma Y)_{\epsilon}$.

Proof.

By Proposition 10, there exists a measurable function $\lambda:\mathbb{R}^{n}\to\mathbb{R}$, $-1\leq\lambda\leq 1$ such that $p=\operatorname{\mathbb{E}}(X\lambda(X))$. Let

[TABLE]

By Proposition 10, $z\in\Gamma Y$.

We have $\operatorname{\mathbb{E}}_{X^{(i)}}(X^{(i)}\lambda(X^{(i)}))=p$ and, for every coordinate $j$,

[TABLE]

By Lemma 4 and for any fixed coordinate $j$ we have, over the choice of $(X^{(i)})_{i=1}^{N}$,

[TABLE]

whenever $N\geq(8M\sqrt{n}/\epsilon)^{\frac{1}{2}+\frac{1}{\gamma}}$. A union bound over ${n}$ choices of $j$ gives:

[TABLE]

So $\mathbb{P}({\lVert p-z\rVert}\geq\epsilon)\leq\delta$ whenever

[TABLE]

and $N\geq(8M\sqrt{n}/\epsilon)^{\frac{1}{2}+\frac{1}{\gamma}}$. The claim follows. ∎