Independent component analysis

The data are represented by the observed random vector ${\displaystyle {\boldsymbol {x}}=(x_{1},\ldots ,x_{m})^{T}}$ and the hidden components as the random vector ${\displaystyle {\boldsymbol {s}}=(s_{1},\ldots ,s_{n})^{T}.}$ The task is to transform the observed data ${\displaystyle {\boldsymbol {x}},}$ using a linear static transformation ${\displaystyle {\boldsymbol {W}}}$ as ${\displaystyle {\boldsymbol {s}}={\boldsymbol {W}}{\boldsymbol {x}},}$ into a vector of maximally independent components ${\displaystyle {\boldsymbol {s}}}$ measured by some function ${\displaystyle F(s_{1},\ldots ,s_{n})}$ of independence.

The components ${\displaystyle x_{i}}$ of the observed random vector ${\displaystyle {\boldsymbol {x}}=(x_{1},\ldots ,x_{m})^{T}}$ are generated as a sum of the independent components ${\displaystyle s_{k}}$, ${\displaystyle k=1,\ldots ,n}$:

${\displaystyle x_{i}=a_{i,1}s_{1}+\cdots +a_{i,k}s_{k}+\cdots +a_{i,n}s_{n}}$

weighted by the mixing weights ${\displaystyle a_{i,k}}$.

The same generative model can be written in vector form as ${\displaystyle {\boldsymbol {x}}=\sum _{k=1}^{n}s_{k}{\boldsymbol {a}}_{k}}$, where the observed random vector ${\displaystyle {\boldsymbol {x}}}$ is represented by the basis vectors ${\displaystyle {\boldsymbol {a}}_{k}=({\boldsymbol {a}}_{1,k},\ldots ,{\boldsymbol {a}}_{m,k})^{T}}$. The basis vectors ${\displaystyle {\boldsymbol {a}}_{k}}$ form the columns of the mixing matrix ${\displaystyle {\boldsymbol {A}}=({\boldsymbol {a}}_{1},\ldots ,{\boldsymbol {a}}_{n})}$ and the generative formula can be written as ${\displaystyle {\boldsymbol {x}}={\boldsymbol {A}}{\boldsymbol {s}}}$, where ${\displaystyle {\boldsymbol {s}}=(s_{1},\ldots ,s_{n})^{T}}$.

Given the model and realizations (samples) ${\displaystyle {\boldsymbol {x}}_{1},\ldots ,{\boldsymbol {x}}_{N}}$ of the random vector ${\displaystyle {\boldsymbol {x}}}$, the task is to estimate both the mixing matrix ${\displaystyle {\boldsymbol {A}}}$ and the sources ${\displaystyle {\boldsymbol {s}}}$. This is done by adaptively calculating the ${\displaystyle {\boldsymbol {w}}}$ vectors and setting up a cost function which either maximizes the non-gaussianity of the calculated ${\displaystyle s_{k}={\boldsymbol {w}}^{T}{\boldsymbol {x}}}$ or minimizes the mutual information. In some cases, a priori knowledge of the probability distributions of the sources can be used in the cost function.

The original sources ${\displaystyle {\boldsymbol {s}}}$ can be recovered by multiplying the observed signals ${\displaystyle {\boldsymbol {x}}}$ with the inverse of the mixing matrix ${\displaystyle {\boldsymbol {W}}={\boldsymbol {A}}^{-1}}$, also known as the unmixing matrix. Here it is assumed that the mixing matrix is square (${\displaystyle n=m}$). If the number of basis vectors is greater than the dimensionality of the observed vectors, ${\displaystyle n>m}$, the task is overcomplete but is still solvable with the pseudo inverse.

With the added assumption of zero-mean and uncorrelated Gaussian noise ${\displaystyle n\sim N(0,\operatorname {diag} (\Sigma ))}$, the ICA model takes the form ${\displaystyle {\boldsymbol {x}}={\boldsymbol {A}}{\boldsymbol {s}}+n}$.

The mixing of the sources does not need to be linear. Using a nonlinear mixing function ${\displaystyle f(\cdot |\theta )}$ with parameters ${\displaystyle \theta }$ the nonlinear ICA model is ${\displaystyle x=f(s|\theta )+n}$.

The independent components are identifiable up to a permutation and scaling of the sources.[6] This identifiability requires that:

• At most one of the sources ${\displaystyle s_{k}}$ is Gaussian,
• The number of observed mixtures, ${\displaystyle m}$, must be at least as large as the number of estimated components ${\displaystyle n}$: ${\displaystyle m\geq n}$. It is equivalent to say that the mixing matrix ${\displaystyle {\boldsymbol {A}}}$ must be of full rank for its inverse to exist.

Signal mixtures tend to have Gaussian probability density functions, and source signals tend to have non-Gaussian probability density functions. Each source signal can be extracted from a set of signal mixtures by taking the inner product of a weight vector and those signal mixtures where this inner product provides an orthogonal projection of the signal mixtures. The remaining challenge is finding such a weight vector. One type of method for doing so is projection pursuit.[10][11]

Projection pursuit seeks one projection at a time such that the extracted signal is as non-Gaussian as possible. This contrasts with ICA, which typically extracts M signals simultaneously from M signal mixtures, which requires estimating a M × M unmixing matrix. One practical advantage of projection pursuit over ICA is that fewer than M signals can be extracted if required, where each source signal is extracted from M signal mixtures using an M-element weight vector.

We can use kurtosis to recover the multiple source signal by finding the correct weight vectors with the use of projection pursuit.

The kurtosis of the probability density function of a signal, for a finite sample, is computed as

${\displaystyle K={\frac {\operatorname {E} [(\mathbf {y} -\mathbf {\overline {y}} )^{4}]}{(\operatorname {E} [(\mathbf {y} -\mathbf {\overline {y}} )^{2}])^{2}}}-3}$

where ${\displaystyle \mathbf {\overline {y}} }$ is the sample mean of ${\displaystyle \mathbf {y} }$, the extracted signals. The constant 3 ensures that Gaussian signals have zero kurtosis, Super-Gaussian signals have positive kurtosis, and Sub-Gaussian signals have negative kurtosis. The denominator is the variance of ${\displaystyle \mathbf {y} }$, and ensures that the measured kurtosis takes account of signal variance. The goal of projection pursuit is to maximize the kurtosis, and make the extracted signal as non-normal as possible.

Using kurtosis as a measure of non-normality, we can now examine how the kurtosis of a signal ${\displaystyle \mathbf {y} =\mathbf {w} ^{T}\mathbf {x} }$ extracted from a set of M mixtures ${\displaystyle \mathbf {x} =(x_{1},x_{2},\ldots ,x_{M})^{T}}$ varies as the weight vector ${\displaystyle \mathbf {w} }$ is rotated around the origin. Given our assumption that each source signal ${\displaystyle \mathbf {s} }$ is super-gaussian we would expect:

1. the kurtosis of the extracted signal ${\displaystyle \mathbf {y} }$ to be maximal precisely when ${\displaystyle \mathbf {y} =\mathbf {s} }$.
2. the kurtosis of the extracted signal ${\displaystyle \mathbf {y} }$ to be maximal when ${\displaystyle \mathbf {w} }$ is orthogonal to the projected axes ${\displaystyle S_{1}}$ or ${\displaystyle S_{2}}$, because we know the optimal weight vector should be orthogonal to a transformed axis ${\displaystyle S_{1}}$ or ${\displaystyle S_{2}}$.

For multiple source mixture signals, we can use kurtosis and Gram-Schmidt Orthogonalization (GSO) to recover the signals. Given M signal mixtures in an M-dimensional space, GSO project these data points onto an (M-1)-dimensional space by using the weight vector. We can guarantee the independence of the extracted signals with the use of GSO.

In order to find the correct value of ${\displaystyle \mathbf {w} }$, we can use gradient descent method. We first of all whiten the data, and transform ${\displaystyle \mathbf {x} }$ into a new mixture ${\displaystyle \mathbf {z} }$, which has unit variance, and ${\displaystyle \mathbf {z} =(z_{1},z_{2},\ldots ,z_{M})^{T}}$. This process can be achieved by applying Singular value decomposition to ${\displaystyle \mathbf {x} }$,

${\displaystyle \mathbf {x} =\mathbf {U} \mathbf {D} \mathbf {V} ^{T}}$

Rescaling each vector ${\displaystyle U_{i}=U_{i}/\operatorname {E} (U_{i}^{2})}$, and let ${\displaystyle \mathbf {z} =\mathbf {U} }$. The signal extracted by a weighted vector ${\displaystyle \mathbf {w} }$ is ${\displaystyle \mathbf {y} =\mathbf {w} ^{T}\mathbf {z} }$. If the weight vector w has unit length, then the variance of y is also 1, that is ${\displaystyle \operatorname {E} [(\mathbf {w} ^{T}\mathbf {z} )^{2}]=1}$. The kurtosis can thus be written as:

${\displaystyle K={\frac {\operatorname {E} [\mathbf {y} ^{4}]}{(\operatorname {E} [\mathbf {y} ^{2}])^{2}}}-3=\operatorname {E} [(\mathbf {w} ^{T}\mathbf {z} )^{4}]-3.}$

The updating process for ${\displaystyle \mathbf {w} }$ is:

${\displaystyle \mathbf {w} _{new}=\mathbf {w} _{old}-\eta \operatorname {E} [\mathbf {z} (\mathbf {w} _{old}^{T}\mathbf {z} )^{3}].}$

where ${\displaystyle \eta }$ is a small constant to guarantee that ${\displaystyle \mathbf {w} }$ converges to the optimal solution. After each update, we normalize ${\displaystyle \mathbf {w} _{new}={\frac {\mathbf {w} _{new}}{|\mathbf {w} _{new}|}}}$, and set ${\displaystyle \mathbf {w} _{old}=\mathbf {w} _{new}}$, and repeat the updating process until convergence. We can also use another algorithm to update the weight vector ${\displaystyle \mathbf {w} }$.

Another approach is using negentropy[12][13] instead of kurtosis. Using negentropy is a more robust method than kurtosis, as kurtosis is very sensitive to outliers. The negentropy methods are based on an important property of Gaussian distribution: a Gaussian variable has the largest entropy among all continuous random variables of equal variance. This is also the reason why we want to find the most nongaussian variables. A simple proof can be found in Differential entropy.

${\displaystyle J(x)=S(y)-S(x)\,}$

y is a Gaussian random variable of the same covariance matrix as x

${\displaystyle S(x)=-\int p_{x}(u)\log p_{x}(u)du}$

An approximation for negentropy is

${\displaystyle J(x)={\frac {1}{12}}(E(x^{3}))^{2}+{\frac {1}{48}}(kurt(x))^{2}}$

A proof can be found in the original papers of Comon;[14][12] it has been reproduced in the book Independent Component Analysis by Aapo Hyvärinen, Juha Karhunen, and Erkki Oja[15] This approximation also suffers from the same problem as kurtosis (sensitivity to outliers). Other approaches have been developed.[16]

${\displaystyle J(y)=k_{1}(E(G_{1}(y)))^{2}+k_{2}(E(G_{2}(y))-E(G_{2}(v))^{2}}$

A choice of ${\displaystyle G_{1}}$ and ${\displaystyle G_{2}}$ are

${\displaystyle G_{1}={\frac {1}{a_{1}}}\log(\cosh(a_{1}u))}$ and ${\displaystyle G_{2}=-\exp(-{\frac {u^{2}}{2}})}$

Infomax ICA[17] is essentially a multivariate, parallel version of projection pursuit. Whereas projection pursuit extracts a series of signals one at a time from a set of M signal mixtures, ICA extracts M signals in parallel. This tends to make ICA more robust than projection pursuit.[18]

The projection pursuit method uses Gram-Schmidt orthogonalization to ensure the independence of the extracted signal, while ICA use infomax and maximum likelihood estimate to ensure the independence of the extracted signal. The Non-Normality of the extracted signal is achieved by assigning an appropriate model, or prior, for the signal.

The process of ICA based on infomax in short is: given a set of signal mixtures ${\displaystyle \mathbf {x} }$ and a set of identical independent model cumulative distribution functions(cdfs) ${\displaystyle g}$, we seek the unmixing matrix ${\displaystyle \mathbf {W} }$ which maximizes the joint entropy of the signals ${\displaystyle \mathbf {Y} =g(\mathbf {y} )}$, where ${\displaystyle \mathbf {y} =\mathbf {Wx} }$ are the signals extracted by ${\displaystyle \mathbf {W} }$. Given the optimal ${\displaystyle \mathbf {W} }$, the signals ${\displaystyle \mathbf {Y} }$ have maximum entropy and are therefore independent, which ensures that the extracted signals ${\displaystyle \mathbf {y} =g^{-1}(\mathbf {Y} )}$ are also independent. ${\displaystyle g}$ is an invertible function, and is the signal model. Note that if the source signal model probability density function ${\displaystyle p_{s}}$ matches the probability density function of the extracted signal ${\displaystyle p_{\mathbf {y} }}$, then maximizing the joint entropy of ${\displaystyle Y}$ also maximizes the amount of mutual information between ${\displaystyle \mathbf {x} }$ and ${\displaystyle \mathbf {Y} }$. For this reason, using entropy to extract independent signals is known as infomax.

Consider the entropy of the vector variable ${\displaystyle \mathbf {Y} =g(\mathbf {y} )}$, where ${\displaystyle \mathbf {y} =\mathbf {Wx} }$ is the set of signals extracted by the unmixing matrix ${\displaystyle \mathbf {W} }$. For a finite set of values sampled from a distribution with pdf ${\displaystyle p_{\mathbf {y} }}$, the entropy of ${\displaystyle \mathbf {Y} }$ can be estimated as:

${\displaystyle H(\mathbf {Y} )=-{\frac {1}{N}}\sum _{t=1}^{N}\ln p_{\mathbf {Y} }(\mathbf {Y} ^{t})}$

The joint pdf ${\displaystyle p_{\mathbf {Y} }}$ can be shown to be related to the joint pdf ${\displaystyle p_{\mathbf {y} }}$ of the extracted signals by the multivariate form:

${\displaystyle p_{\mathbf {Y} }(Y)={\frac {p_{\mathbf {y} }(\mathbf {y} )}{|{\frac {\partial \mathbf {Y} }{\partial \mathbf {y} }}|}}}$

where ${\displaystyle \mathbf {J} ={\frac {\partial \mathbf {Y} }{\partial \mathbf {y} }}}$ is the Jacobian matrix. We have ${\displaystyle |\mathbf {J} |=g'(\mathbf {y} )}$, and ${\displaystyle g'}$ is the pdf assumed for source signals ${\displaystyle g'=p_{s}}$, therefore,

${\displaystyle p_{\mathbf {Y} }(Y)={\frac {p_{\mathbf {y} }(\mathbf {y} )}{|{\frac {\partial \mathbf {Y} }{\partial \mathbf {y} }}|}}={\frac {p_{\mathbf {y} }(\mathbf {y} )}{p_{\mathbf {s} }(\mathbf {y} )}}}$

therefore,

${\displaystyle H(\mathbf {Y} )=-{\frac {1}{N}}\sum _{t=1}^{N}\ln {\frac {p_{\mathbf {y} }(\mathbf {y} )}{p_{\mathbf {s} }(\mathbf {y} )}}}$

We know that when ${\displaystyle p_{\mathbf {y} }=p_{s}}$, ${\displaystyle p_{\mathbf {Y} }}$ is of uniform distribution, and ${\displaystyle H({\mathbf {Y} })}$ is maximized. Since

${\displaystyle p_{\mathbf {y} }(\mathbf {y} )={\frac {p_{\mathbf {x} }(\mathbf {x} )}{|{\frac {\partial \mathbf {y} }{\partial \mathbf {x} }}|}}={\frac {p_{\mathbf {x} }(\mathbf {x} )}{|\mathbf {W} |}}}$

where ${\displaystyle |\mathbf {W} |}$ is the absolute value of the determinant of the unmixing matrix ${\displaystyle \mathbf {W} }$. Therefore,

${\displaystyle H(\mathbf {Y} )=-{\frac {1}{N}}\sum _{t=1}^{N}\ln {\frac {p_{\mathbf {x} }(\mathbf {x} ^{t})}{|\mathbf {W} |p_{\mathbf {s} }(\mathbf {y} ^{t})}}}$

so,

${\displaystyle H(\mathbf {Y} )={\frac {1}{N}}\sum _{t=1}^{N}\ln p_{\mathbf {s} }(\mathbf {y} ^{t})+\ln |\mathbf {W} |+H(\mathbf {x} )}$

since ${\displaystyle H(\mathbf {x} )=-{\frac {1}{N}}\sum _{t=1}^{N}\ln p_{\mathbf {x} }(\mathbf {x} ^{t})}$, and maximizing ${\displaystyle \mathbf {W} }$ does not affect ${\displaystyle H_{\mathbf {x} }}$, so we can maximize the function

${\displaystyle h(\mathbf {Y} )={\frac {1}{N}}\sum _{t=1}^{N}\ln p_{\mathbf {s} }(\mathbf {y} ^{t})+\ln |\mathbf {W} |}$

to achieve the independence of extracted signal.

If there are M marginal pdfs of the model joint pdf ${\displaystyle p_{\mathbf {s} }}$ are independent and use the commonly super-gaussian model pdf for the source signals ${\displaystyle p_{\mathbf {s} }=(1-\tanh(\mathbf {s} )^{2})}$, then we have

${\displaystyle h(\mathbf {Y} )={\frac {1}{N}}\sum _{i=1}^{M}\sum _{t=1}^{N}\ln(1-\tanh(\mathbf {w_{i}^{T}x^{t}} )^{2})+\ln |\mathbf {W} |}$

In the sum, given an observed signal mixture ${\displaystyle \mathbf {x} }$, the corresponding set of extracted signals ${\displaystyle \mathbf {y} }$ and source signal model ${\displaystyle p_{\mathbf {s} }=g'}$, we can find the optimal unmixing matrix ${\displaystyle \mathbf {W} }$, and make the extracted signals independent and non-gaussian. Like the projection pursuit situation, we can use gradient descent method to find the optimal solution of the unmixing matrix.

Maximum likelihood estimation (MLE) is a standard statistical tool for finding parameter values (e.g. the unmixing matrix ${\displaystyle \mathbf {W} }$) that provide the best fit of some data (e.g., the extracted signals ${\displaystyle y}$) to a given a model (e.g., the assumed joint probability density function (pdf) ${\displaystyle p_{s}}$ of source signals).[18]

The ML "model" includes a specification of a pdf, which in this case is the pdf ${\displaystyle p_{s}}$ of the unknown source signals ${\displaystyle s}$. Using ML ICA, the objective is to find an unmixing matrix that yields extracted signals ${\displaystyle y=\mathbf {W} x}$ with a joint pdf as similar as possible to the joint pdf ${\displaystyle p_{s}}$ of the unknown source signals ${\displaystyle s}$.

MLE is thus based on the assumption that if the model pdf ${\displaystyle p_{s}}$ and the model parameters ${\displaystyle \mathbf {A} }$ are correct then a high probability should be obtained for the data ${\displaystyle x}$ that were actually observed. Conversely, if ${\displaystyle \mathbf {A} }$ is far from the correct parameter values then a low probability of the observed data would be expected.

Using MLE, we call the probability of the observed data for a given set of model parameter values (e.g., a pdf ${\displaystyle p_{s}}$ and a matrix ${\displaystyle \mathbf {A} }$) the likelihood of the model parameter values given the observed data.

We define a likelihood function ${\displaystyle \mathbf {L(W)} }$ of ${\displaystyle \mathbf {W} }$:

${\displaystyle \mathbf {L(W)} =p_{s}(\mathbf {W} x)|\det \mathbf {W} |.}$

This equals to the probability density at ${\displaystyle x}$, since ${\displaystyle s=\mathbf {W} x}$.

Thus, if we wish to find a ${\displaystyle \mathbf {W} }$ that is most likely to have generated the observed mixtures ${\displaystyle x}$ from the unknown source signals ${\displaystyle s}$ with pdf ${\displaystyle p_{s}}$ then we need only find that ${\displaystyle \mathbf {W} }$ which maximizes the likelihood ${\displaystyle \mathbf {L(W)} }$. The unmixing matrix that maximizes equation is known as the MLE of the optimal unmixing matrix.

It is common practice to use the log likelihood, because this is easier to evaluate. As the logarithm is a monotonic function, the ${\displaystyle \mathbf {W} }$ that maximizes the function ${\displaystyle \mathbf {L(W)} }$ also maximizes its logarithm ${\displaystyle \ln \mathbf {L(W)} }$. This allows us to take the logarithm of equation above, which yields the log likelihood function

${\displaystyle \ln \mathbf {L(W)} =\sum _{i}\sum _{t}\ln p_{s}(w_{i}^{T}x_{t})+N\ln |\det \mathbf {W} |}$

If we substitute a commonly used high-Kurtosis model pdf for the source signals ${\displaystyle p_{s}=(1-\tanh(s)^{2})}$ then we have

${\displaystyle \ln \mathbf {L(W)} ={1 \over N}\sum _{i}^{M}\sum _{t}^{N}\ln(1-\tanh(w_{i}^{T}x_{t})^{2})+\ln |\det \mathbf {W} |}$

This matrix ${\displaystyle \mathbf {W} }$ that maximizes this function is the maximum likelihood estimation.