Using basic Tensorlab features for ICA¶

Instantaneous mixtures¶

An instantaneous mixture \(\vec{x}(t) \in \R^{Q}\) of \(R\) independent components \(\vec{s}(t) \in \R^{R}\) can be written as

in which \(\mat{M}\) contains the mixing coefficients and temporally i.i.d. noise is captured by the additive noise term \(\vec{n}(t)\). We consider two auxiliary signals that are identically uniformly distributed over \([-1,1]\). Eeach signal has \(10^4\) samples. The two actual source signals are generated by passing these auxiliary signals through finite impulse response filters with coefficients \([-1,1,1]\) and \([1,1,1]\), respectively:

N = 10000; R = 2;
Saux = rand(N,R)*2-1;
S(:,1) = filter([-1 1 1],1,Saux(:,1));
S(:,2) = filter([1 1 1],1,Saux(:,2));

Figure 69 shows the source histograms. The samples are confined to the interval \([-3,3]\). The mixing matrix is taken equal to

Independent identically distributed Gaussian noise is added with a signal-to-noise-ratio of 10 dB:

M = [1/sqrt(2) 2/sqrt(5); -1/sqrt(2) 1/sqrt(5)];
X = S*M.';
SNR = 10;
Xnoisy = noisy(X,SNR);

Figure 69: Histogram of the source samples.

Second-order statistics¶

Consider a tensor \(\ten{C}^{(2)}_x\) of which the slices are lagged covariance matrices. This tensor can be generated with the scov command:

lags = [0 1 2];
C2x = scov(Xnoisy,lags);

Each observed lagged covariance matrix equals the corresponding lagged covariance matrix of the source signals, multiplied by the mixing matrix in the first and second mode. As the lagged covariance matrices of the source signals are approximately diagonal, the tensor is approximately low-rank, as illustrated in Figure 70.

Figure 70: Tensorizing the data by stacking lagged covariance matrices.

The mixing matrix can be recovered by decomposing \(\ten{C}^{(2)}_x\) in a minimal number of rank-1 terms. We use the basic cpd command without exploiting the symmetry; hence, two estimates of the mixing matrix are obtained:

Uest = cpd(C2x,R);
Mest_mode1 = Uest{1};
Mest_mode2 = Uest{2};

Since the mixing matrix is square non-singular, the source signals can subsequently be estimated as

The cpderr command can be used to compute the relative error on the estimates of the mixing matrix, and to fix the indeterminacies for visualization purposes:

[err_scov_mode1,P1,D1,Mestscaled] = cpderr(M,Mest_mode1)
[err_scov_mode2,P2,D2] = cpderr(M,Mest_mode2)
Srec = Xnoisy*pinv(Mestscaled.');

The joint distributions of the signals are visualized in Figure 71. The joint distribution of the source signals is confined to the square \([-3,3]\times [-3,3]\), and the components vary independently along the horizontal and vertical axis according to their distribution function (cf. Figure 69). The edges of the figure in the middle illustrate that the observed signals cannot take values independent of one another. As a matter of fact, the observed joint distribution is equal to the joint distribution of the source signals, linearly transformed by the mixing matrix \(\mat{M}\). The axes correspond to the columns of the identity matrix (left), the mixing matrix \(\mat{M}\) (middle), and the matrix \(\tilde{\mat{M}}^{-1}\cdot\mat{M}\) (right), with \(\tilde{\mat{M}}\) the estimate of \(\mat{M}\). The figure on the right shows that the estimated source signals are close to independent.

Figure 71: Visualization of the joint distributions of the source signals (left), the observed signals (middle) and the estimated source signals (right). The axes correspond to the columns of the identity matrix (left), the mixing matrix \(\mat{M}\) (middle), and the matrix \(\tilde{\mat{M}}^{-1}\cdot\mat{M}\) (right), with \(\tilde{\mat{M}}\) the estimate of \(\mat{M}\).

Fourth-order statistics¶

Consider the fourth-order cumulant \(\ten{C}^{(4)}_x\) of the observed signals, which can be computed using the cum4 command:

C4x = cum4(Xnoisy);

The fourth-order cumulant of the sources is approximately diagonal because of the mutual statistical independence; hence, the tensor \(\ten{C}^{(4)}_x\) can be decomposed in rank-1 terms as follows:

with \(\ten{C}^{(4)}_s = \text{diag}_4(c_{s_1}^{(4)},c_{s_2}^{(4)},\ldots,c_{s_R}^{(4)})\). All four factor matrices are equal to the mixing matrix. We only use the first factor matrix as an estimate:

Uest = cpd(C4x,R);
Mest_mode1 = Uest{1};
[err_cum,P,D] = cpderr(M,Mest_mode1)

Second-order statistics and missing values¶

If infinitely many samples are available, the i.i.d. Gaussian noise only affects the values on the diagonal of the covariance matrix for lag zero. These entries can be omitted from the tensor \(\ten{C}^{(2)}_x\) as follows:

C2xi = C2x;
for i = 1:size(C2x,1)
  C2xi(i,i,1) = NaN;
end

The cpd command can be used to decompose the incomplete tensor C2xi:

Uest = cpd(C2xi,R);
Mest_mode1 = Uest{1};
[err_scovi,P,D] = cpderr(M,Mest_mode1)

The demo Using advanced Tensorlab features for ICA elaborates on noise modeling using advanced Tensorlab features.