Our next topic is Independent Components Analysis (ICA). Similar to PCA, this will find a new basis in which to represent our data. However, the goal is very different.

As a motivating example, consider the "cocktail party problem." Here, $d$$d$dd speakers are speaking simultaneously at a party, and any microphone placed in the room records only an overlapping combination of the $d$$d$dd speakers' voices. But lets say we have $d$$d$dd different microphones placed in the room, and because each microphone is a different distance from each of the speakers, it records a different combination of the speakers' voices. Using these microphone recordings, can we separate out the original $d$$d$dd speakers' speech signals?

To formalize this problem, we imagine that there is some data $s\in {\mathbb{R}}^d$$s\in {\mathbb{R}}^d$s inR^(d)s \in \mathbb{R}^{d} that is generated via $d$$d$dd independent sources. What we observe is

where $A$$A$AA is an unknown square matrix called the mixing matrix. Repeated observations gives us a dataset $\{x^{(i)};i=1,\dots ,n\}$$\left\{x^{(i)};i=1,\dots ,n\right\}${x^((i));i=1,dots,n}\left\{x^{(i)} ; i=1, \ldots, n\right\}, and our goal is to recover the sources $s^{(i)}$$s^{(i)}$s^((i))s^{(i)} that had generated our data $(x^{(i)}=A{s}^{(i)})$$\left(x^{(i)}=A{s}^{(i)}\right)$(x^((i))=As^((i)))\left(x^{(i)}=A s^{(i)}\right).

In our cocktail party problem, $s^{(i)}$$s^{(i)}$s^((i))s^{(i)} is an $d$$d$dd-dimensional vector, and $s_j^{(i)}$$s_j^{(i)}$s_(j)^((i))s_{j}^{(i)} is the sound that speaker $j$$j$jj was uttering at time $i$$i$ii. Also, $x^{(i)}$$x^{(i)}$x^((i))x^{(i)} in an $d$$d$dd-dimensional vector, and $x_j^{(i)}$$x_j^{(i)}$x_(j)^((i))x_{j}^{(i)} is the acoustic reading recorded by microphone $j$$j$jj at time $i$$i$ii.

Let $W=A^{-1}$$W=A^{-1}$W=A^(-1)W=A^{-1} be the unmixing matrix. Our goal is to find $W$$W$WW, so that given our microphone recordings $x^{(i)}$$x^{(i)}$x^((i))x^{(i)}, we can recover the sources by computing $s^{(i)}=W{x}^{(i)}$$s^{(i)}=W{x}^{(i)}$s^((i))=Wx^((i))s^{(i)}=W x^{(i)}. For notational convenience, we also let $w_i^T$$w_i^T$w_(i)^(T)w_{i}^{T} denote the $i$$i$ii-th row of $W$$W$WW, so that

Thus, $w_i\in {\mathbb{R}}^d$$w_i\in {\mathbb{R}}^d$w_(i)inR^(d)w_{i} \in \mathbb{R}^{d}, and the $j$$j$jj-th source can be recovered as $s_j^{(i)}=w_j^T{x}^{(i)}$$s_j^{(i)}=w_j^T{x}^{(i)}$s_(j)^((i))=w_(j)^(T)x^((i))s_{j}^{(i)}=w_{j}^{T} x^{(i)}.

To what degree can $W=A^{-1}$$W=A^{-1}$W=A^(-1)W=A^{-1} be recovered? If we have no prior knowledge about the sources and the mixing matrix, it is easy to see that there are some inherent ambiguities in $A$$A$AA that are impossible to recover, given only the $x^{(i)}$$x^{(i)}$x^((i))x^{(i)}'s.

Specifically, let $P$$P$PP be any $d$$d$dd-by-$d$$d$dd permutation matrix. This means that each row and each column of $P$$P$PP has exactly one "1." Here are some examples of permutation matrices:

If $z$$z$zz is a vector, then $Pz$$Pz$PzPz is another vector that contains a permuted version of $z$$z$zz's coordinates. Given only the $x^{(i)}$$x^{(i)}$x^((i))x^{(i)}'s, there will be no way to distinguish between $W$$W$WW and $PW$$PW$PWP W. Specifically, the permutation of the original sources is ambiguous, which should be no surprise. Fortunately, this does not matter for most applications.

Further, there is no way to recover the correct scaling of the $w_i$$w_i$w_(i)w_{i}'s. For instance, if $A$$A$AA were replaced with $2A$$2A$2A2 A, and every $s^{(i)}$$s^{(i)}$s^((i))s^{(i)} were replaced with $(0.5)s^{(i)}$$(0.5)s^{(i)}$(0.5)s^((i))(0.5) s^{(i)}, then our observed $x^{(i)}=2A\cdot (0.5)s^{(i)}$$x^{(i)}=2A\cdot (0.5)s^{(i)}$x^((i))=2A*(0.5)s^((i))x^{(i)}=2 A \cdot(0.5) s^{(i)} would still be the same. More broadly, if a single column of $A$$A$AA were scaled by a factor of $\alpha$$\alpha$alpha\alpha, and the corresponding source were scaled by a factor of $1/\alpha$$1/\alpha$1//alpha1 / \alpha, then there is again no way to determine that this had happened given only the $x^{(i)}$$x^{(i)}$x^((i))x^{(i)}'s. Thus, we cannot recover the "correct" scaling of the sources. However, for the applications that we are concerned with – including the cocktail party problem – this ambiguity also does not matter. Specifically, scaling a speaker's speech signal $s_j^{(i)}$$s_j^{(i)}$s_(j)^((i))s_{j}^{(i)} by some positive factor $\alpha$$\alpha$alpha\alpha affects only the volume of that speaker's speech. Also, sign changes do not matter, and $s_j^{(i)}$$s_j^{(i)}$s_(j)^((i))s_{j}^{(i)} and $-s_j^{(i)}$$-s_j^{(i)}$-s_(j)^((i))-s_{j}^{(i)} sound identical when played on a speaker. Thus, if the $w_i$$w_i$w_(i)w_{i} found by an algorithm is scaled by any non-zero real number, the corresponding recovered source $s_i=w_i^{T}x$$s_i=w_i^{T}x$s_(i)=w_(i)^(T)xs_{i}=w_{i}^{T} x will be scaled by the same factor; but this usually does not matter. (These comments also apply to ICA for the brain/MEG data that we talked about in class.)

Are these the only sources of ambiguity in ICA? It turns out that they are, so long as the sources $s_i$$s_i$s_(i)s_{i} are non-Gaussian. To see what the difficulty is with Gaussian data, consider an example in which $n=2$$n=2$n=2n=2, and $s\sim \mathcal{N}(0,I)$$s\sim \mathcal{N}(0,I)$s∼N(0,I)s \sim \mathcal{N}(0, I). Here, $I$$I$II is the $2\mathrm{x}2$$2\mathrm{x}2$2x22 \mathrm{x} 2 identity matrix. Note that the contours of the density of the standard normal distribution $\mathcal{N}(0,I)$$\mathcal{N}(0,I)$N(0,I)\mathcal{N}(0, I) are circles centered on the origin, and the density is rotationally symmetric.

Now, suppose we observe some $x=As$$x=As$x=Asx=A s, where $A$$A$AA is our mixing matrix. Then, the distribution of $x$$x$xx will be Gaussian, $x\sim \mathcal{N}(0,A{A}^T)$$x\sim \mathcal{N}\left(0,A{A}^T\right)$x∼N(0,AA^(T))x \sim \mathcal{N}\left(0, A A^{T}\right), since

Now, let $R$$R$RR be an arbitrary orthogonal (less formally, a rotation/reflection) matrix, so that $R{R}^T=R^{T}R=I$$R{R}^T=R^{T}R=I$RR^(T)=R^(T)R=IR R^{T}=R^{T} R=I, and let $A^{\mathrm{\prime }}=AR$$A^{\mathrm{\prime }}=AR$A^(')=ARA^{\prime}=A R. Then if the data had been mixed according to $A^{\mathrm{\prime }}$$A^{\mathrm{\prime }}$A^(')A^{\prime} instead of $A$$A$AA, we would have instead observed $x^{\mathrm{\prime }}=A^{\mathrm{\prime }}s$$x^{\mathrm{\prime }}=A^{\mathrm{\prime }}s$x^(')=A^(')sx^{\prime}=A^{\prime} s. The distribution of $x^{\mathrm{\prime }}$$x^{\mathrm{\prime }}$x^(')x^{\prime} is also Gaussian, $x^{\mathrm{\prime }}\sim \mathcal{N}(0,A{A}^T)$$x^{\mathrm{\prime }}\sim \mathcal{N}\left(0,A{A}^T\right)$x^(')∼N(0,AA^(T))x^{\prime} \sim \mathcal{N}\left(0, A A^{T}\right), since ${\mathbb{E}}_{s\sim \mathcal{N}(0,I)}\left[x^{\mathrm{\prime }}{\left(x^{\mathrm{\prime }}\right)}^T\right]=\mathbb{E}\left[A^{\mathrm{\prime }}s{s}^T{\left(A^{\mathrm{\prime }}\right)}^T\right]=\mathbb{E}[ARs{s}^T(AR{)}^T]=AR{R}^T{A}^T=A{A}^T.$${\mathbb{E}}_{s\sim \mathcal{N}(0,I)}\left[x^{\mathrm{\prime }}{\left(x^{\mathrm{\prime }}\right)}^T\right]=\mathbb{E}\left[A^{\mathrm{\prime }}s{s}^T{\left(A^{\mathrm{\prime }}\right)}^T\right]=\mathbb{E}\left[ARs{s}^T(AR{)}^T\right]=AR{R}^T{A}^T=A{A}^T.$E_(s∼N(0,I))[x^(')(x^('))^(T)]=E[A^(')ss^(T)(A^('))^(T)]=E[ARss^(T)(AR)^(T)]=ARR^(T)A^(T)=AA^(T).\mathbb{E}_{s \sim \mathcal{N}(0, I)}\left[x^{\prime}\left(x^{\prime}\right)^{T}\right]=\mathbb{E}\left[A^{\prime} s s^{T}\left(A^{\prime}\right)^{T}\right]=\mathbb{E}\left[A R s s^{T}(A R)^{T}\right]=A R R^{T} A^{T}=A A^{T} . Hence, whether the mixing matrix is $A$$A$AA or $A^{\mathrm{\prime }}$$A^{\mathrm{\prime }}$A^(')A^{\prime}, we would observe data from a $\mathcal{N}(0,A{A}^T)$$\mathcal{N}\left(0,A{A}^T\right)$N(0,AA^(T))\mathcal{N}\left(0, A A^{T}\right) distribution. Thus, there is no way to tell if the sources were mixed using $A$$A$AA and $A^{\mathrm{\prime }}$$A^{\mathrm{\prime }}$A^(')A^{\prime}. There is an arbitrary rotational component in the mixing matrix that cannot be determined from the data, and we cannot recover the original sources.

Our argument above was based on the fact that the multivariate standard normal distribution is rotationally symmetric. Despite the bleak picture that this paints for ICA on Gaussian data, it turns out that, so long as the data is not Gaussian, it is possible, given enough data, to recover the $d$$d$dd independent sources.

Before moving on to derive the ICA algorithm proper, we first digress briefly to talk about the effect of linear transformations on densities.

Suppose a random variable $s$$s$ss is drawn according to some density $p_s(s)$$p_s(s)$p_(s)(s)p_{s}(s). For simplicity, assume for now that $s\in \mathbb{R}$$s\in \mathbb{R}$s inRs \in \mathbb{R} is a real number. Now, let the random variable $x$$x$xx be defined according to $x=As$$x=As$x=Asx=A s (here, $x\in \mathbb{R},A\in \mathbb{R}$$x\in \mathbb{R},A\in \mathbb{R}$x inR,A inRx \in \mathbb{R}, A \in \mathbb{R}). Let $p_x$$p_x$p_(x)p_{x} be the density of $x$$x$xx. What is $p_x$$p_x$p_(x)p_{x}?

Let $W=A^{-1}$$W=A^{-1}$W=A^(-1)W=A^{-1}. To calculate the "probability" of a particular value of $x$$x$xx, it is tempting to compute $s=Wx$$s=Wx$s=Wxs=W x, then then evaluate $p_s$$p_s$p_(s)p_{s} at that point, and conclude that "$p_x(x)=p_s(Wx)$$p_x(x)=p_s(Wx)$p_(x)(x)=p_(s)(Wx)p_{x}(x)=p_{s}(W x)." However, this is incorrect. For example, let $s\sim$$s\sim$s∼s \sim Uniform $[0,1]$$[0,1]$[0,1][0,1], so $p_s(s)=1\{0\le s\le 1\}$$p_s(s)=1\{0\le s\le 1\}$p_(s)(s)=1{0 <= s <= 1}p_{s}(s)=1\{0 \leq s \leq 1\}. Now, let $A=2$$A=2$A=2A=2, so $x=2s$$x=2s$x=2sx=2 s. Clearly, $x$$x$xx is distributed uniformly in the interval $[0,2]$$[0,2]$[0,2][0,2]. Thus, its density is given by $p_x(x)=(0.5)1\{0\le x\le 2\}$$p_x(x)=(0.5)1\{0\le x\le 2\}$p_(x)(x)=(0.5)1{0 <= x <= 2}p_{x}(x)=(0.5) 1\{0 \leq x \leq 2\}. This does not equal $p_s(Wx)$$p_s(Wx)$p_(s)(Wx)p_{s}(W x), where $W=0.5=A^{-1}$$W=0.5=A^{-1}$W=0.5=A^(-1)W=0.5=A^{-1}. Instead, the correct formula is $p_x(x)=p_s(Wx)|W|$$p_x(x)=p_s(Wx)|W|$p_(x)(x)=p_(s)(Wx)|W|p_{x}(x)=p_{s}(W x)|W|.

More generally, if $s$$s$ss is a vector-valued distribution with density $p_s$$p_s$p_(s)p_{s}, and $x=As$$x=As$x=Asx=A s for a square, invertible matrix $A$$A$AA, then the density of $x$$x$xx is given by

where $W=A^{-1}$$W=A^{-1}$W=A^(-1)W=A^{-1}.

Remark. If you're seen the result that $A$$A$AA maps $[0,1{]}^d$$[0,1{]}^d$[0,1]^(d)[0,1]^{d} to a set of volume $|A|$$|A|$|A||A|, then here's another way to remember the formula for $p_x$$p_x$p_(x)p_{x} given above, that also generalizes our previous 1-dimensional example. Specifically, let $A\in {\mathbb{R}}^{d\times d}$$A\in {\mathbb{R}}^{d\times d}$A inR^(d xx d)A \in \mathbb{R}^{d \times d} be given, and let $W=A^{-1}$$W=A^{-1}$W=A^(-1)W=A^{-1} as usual. Also let $C_1=[0,1{]}^d$$C_1=[0,1{]}^d$C_(1)=[0,1]^(d)C_{1}=[0,1]^{d} be the $d$$d$dd-dimensional hypercube, and define $C_2=\{As:s\in C_1\}\subseteq {\mathbb{R}}^d$$C_2=\left\{As:s\in C_1\right\}\subseteq {\mathbb{R}}^d$C_(2)={As:s inC_(1)}subeR^(d)C_{2}=\left\{A s: s \in C_{1}\right\} \subseteq \mathbb{R}^{d} to be the image of $C_1$$C_1$C_(1)C_{1} under the mapping given by $A$$A$AA. Then it is a standard result in linear algebra (and, indeed, one of the ways of defining determinants) that the volume of $C_2$$C_2$C_(2)C_{2} is given by $|A|$$|A|$|A||A|. Now, suppose $s$$s$ss is uniformly distributed in $[0,1{]}^d$$[0,1{]}^d$[0,1]^(d)[0,1]^{d}, so its density is $p_s(s)=1\{s\in C_1\}$$p_s(s)=1\left\{s\in C_1\right\}$p_(s)(s)=1{s inC_(1)}p_{s}(s)=1\left\{s \in C_{1}\right\}. Then clearly $x$$x$xx will be uniformly distributed in $C_2$$C_2$C_(2)C_{2}. Its density is therefore found to be $p_x(x)=1\{x\in C_2\}/\mathrm{vol}\left(C_2\right)$$p_x(x)=1\left\{x\in C_2\right\}/\mathrm{vol}\left(C_2\right)$p_(x)(x)=1{x inC_(2)}//vol(C_(2))p_{x}(x)=1\left\{x \in C_{2}\right\} / \operatorname{vol}\left(C_{2}\right) (since it must integrate over $C_2$$C_2$C_(2)C_{2} to 1$)$$)$)). But using the fact that the determinant of the inverse of a matrix is just the inverse of the determinant, we have $1/\mathrm{vol}\left(C_2\right)=1/|A|=\left|A^{-1}\right|=|W|$$1/\mathrm{vol}\left(C_2\right)=1/|A|=\left|A^{-1}\right|=|W|$1//vol(C_(2))=1//|A|=|A^(-1)|=|W|1 /\operatorname{vol}\left(C_{2}\right)=1 /|A|=\left|A^{-1}\right|=|W|. Thus, $p_x(x)=1\{x\in C_2\}|W|=1\{Wx\in C_1\}|W|=p_s(Wx)|W|$$p_x(x)=1\left\{x\in C_2\right\}|W|=1\{Wx\in \left.C_1\right\}|W|=p_s(Wx)|W|$p_(x)(x)=1{x inC_(2)}|W|=1{Wx in{:C_(1)}|W|=p_(s)(Wx)|W|p_{x}(x)=1\left\{x \in C_{2}\right\}|W|=1\{W x \in\left.C_{1}\right\}|W|=p_{s}(W x)|W|.

We are now ready to derive an ICA algorithm. We describe an algorithm by Bell and Sejnowski, and we give an interpretation of their algorithm as a method for maximum likelihood estimation. (This is different from their original interpretation involving a complicated idea called the infomax principal which is no longer necessary given the modern understanding of ICA.)

We suppose that the distribution of each source $s_j$$s_j$s_(j)s_{j} is given by a density $p_s$$p_s$p_(s)p_{s}, and that the joint distribution of the sources $s$$s$ss is given by

Note that by modeling the joint distribution as a product of marginals, we capture the assumption that the sources are independent. Using our formulas from the previous section, this implies the following density on $x=As=$$x=As=$x=As=x=A s= $W^{-1}s$$W^{-1}s$W^(-1)sW^{-1} s:

All that remains is to specify a density for the individual sources $p_s$$p_s$p_(s)p_{s}.

Recall that, given a real-valued random variable $z$$z$zz, its cumulative distribution function (cdf) $F$$F$FF is defined by $F\left(z_0\right)=P(z\le z_0)={\int }_{-\mathrm{\infty }}^{z_0}{p}_z(z)dz$$F\left(z_0\right)=P\left(z\le z_0\right)={\int }_{-\mathrm{\infty }}^{z_0} p_z(z)dz$F(z_(0))=P(z <= z_(0))=int_(-oo)^(z_(0))p_(z)(z)dzF\left(z_{0}\right)=P\left(z \leq z_{0}\right)=\int_{-\infty}^{z_{0}} p_{z}(z) d z and the density is the derivative of the cdf: $p_z(z)=F^{\mathrm{\prime }}(z)$$p_z(z)=F^{\mathrm{\prime }}(z)$p_(z)(z)=F^(')(z)p_{z}(z)=F^{\prime}(z).

Thus, to specify a density for the $s_i$$s_i$s_(i)s_{i}'s, all we need to do is to specify some cdf for it. A cdf has to be a monotonic function that increases from zero to one. Following our previous discussion, we cannot choose the Gaussian cdf, as ICA doesn't work on Gaussian data. What we'll choose instead as a reasonable "default" cdf that slowly increases from 0 to 1, is the sigmoid function $g(s)=1/(1+e^{-s})$$g(s)=1/\left(1+e^{-s}\right)$g(s)=1//(1+e^(-s))g(s)=1 /\left(1+e^{-s}\right). Hence, $p_s(s)=g^{\mathrm{\prime }}(s)$$p_s(s)=g^{\mathrm{\prime }}(s)$p_(s)(s)=g^(')(s)p_{s}(s)=g^{\prime}(s).^[1]

The square matrix $W$$W$WW is the parameter in our model. Given a training set $\{x^{(i)};i=1,\dots ,n\}$$\left\{x^{(i)};i=1,\dots ,n\right\}${x^((i));i=1,dots,n}\left\{x^{(i)} ; i=1, \ldots, n\right\}, the log likelihood is given by

We would like to maximize this in terms $W$$W$WW. By taking derivatives and using the fact (from the first set of notes) that ${\mathrm{\nabla }}_W|W|=|W|{\left(W^{-1}\right)}^T$${\mathrm{\nabla }}_W|W|=|W|{\left(W^{-1}\right)}^T$grad_(W)|W|=|W|(W^(-1))^(T)\nabla_{W}|W|=|W|\left(W^{-1}\right)^{T}, we easily derive a stochastic gradient ascent learning rule. For a training example $x^{(i)}$$x^{(i)}$x^((i))x^{(i)}, the update rule is:

where $\alpha$$\alpha$alpha\alpha is the learning rate.

After the algorithm converges, we then compute $s^{(i)}=W{x}^{(i)}$$s^{(i)}=W{x}^{(i)}$s^((i))=Wx^((i))s^{(i)}=W x^{(i)} to recover the original sources.

Remark. When writing down the likelihood of the data, we implicitly assumed that the $x^{(i)}$$x^{(i)}$x^((i))x^{(i)}'s were independent of each other (for different values of $i$$i$ii; note this issue is different from whether the different coordinates of $x^{(i)}$$x^{(i)}$x^((i))x^{(i)} are independent), so that the likelihood of the training set was given by $\prod _{i}p(x^{(i)};W)$$\prod _i p\left(x^{(i)};W\right)$prod_(i)p(x^((i));W)\prod_{i} p\left(x^{(i)} ; W\right). This assumption is clearly incorrect for speech data and other time series where the $x^{(i)}$$x^{(i)}$x^((i))x^{(i)}'s are dependent, but it can be shown that having correlated training examples will not hurt the performance of the algorithm if we have sufficient data. However, for problems where successive training examples are correlated, when implementing stochastic gradient ascent, it sometimes helps accelerate convergence if we visit training examples in a randomly permuted order. (I.e., run stochastic gradient ascent on a randomly shuffled copy of the training set.)