Combining blind source extraction with joint approximate diagonalization: Thin algorithms for ICA

COMBINING BLIND SOURCE EXTRACTION WITH JOINT APPROXIMATE DIAGONALIZATION: THIN ALGORITHMS FOR ICA Sergio Cruces Signal Processing Group Ingenieros, Camino Descubrimientos, 41092-Seville, Spain. http://viento.us.es/˜sergio ser[email protected] Andrzej Cichocki Brain Science Institute RIKEN, 2-1 Hirosawa, Wako-shi, Saitama, 351-0198 Japan. http://www.bsp.brain.riken.go.jp/ [email protected] ABSTRACT In this paper a multivariate contrast function is proposed for the blind signal extraction of a subset of the independent components from a linear mixture. This contrast combines the robustness of the joint approximate diagonalization techniques with the flexibility of the methods for blind signal extraction. Its maximization leads to hierarchical and simultaneous ICA extraction algorithms which are respectively based on the thin QR and thin SVD factorizations. The interesting similarities and differences with other existing contrasts and algorithms are commented. 1. INTRODUCTION In the last decades powerful criteria and algorithms has been developed to solve the problem of the analysis of the independent components in a linear mixture [1, 2, 3]. In the history of this problem one may distinguish to different approaches: the first one is usually named blind signal separation (BSS) and consists in the simultaneous estimation of all the latent independent components from the mixture; the second one is known as blind signal extraction (BSE) and consists in the estimation of only a subset of the independent components. The BSE problem can be regarded as more general and flexible than BSS because of the following two reasons: 1) BSE includes BSS as the particular the case where one is interested in all the independent components, 2) BSE has computational advantages over BSS if only a small subset of the independent components is of interest. These computational savings can be important in applications like MEG and EEG, where the decomposition of the observations considers a large number of possible independent components but only a few of them are of interest. This work has been supported by the CICYT project of the Government of Spain (Grant TIC2001-0751-C04-04). The original criteria for blind signal extraction were based in the hierarchical the recovery of the independent components [4, 5] one by one alternating extraction with deflation. These criteria have been extended to allow the simultaneous extraction of an arbitrary number of independent components [6, 7]. Another desired property of any BSS/BSE criterion is its robustness, understood in the sense of its ability to give accurate estimates from the available data. Popular techniques for blind signal separation [8, 9, 10] are robust in the previous sense, having criteria based on the joint approximate diagonalization of several cumulant slices. However, up to our knowledge, no equivalent approach for blind signal extraction has been obtained yet. Thus, the purpose of this contribution is to present algorithms that perform the simultaneously extraction a subset of independent components from the mixture by combining the information of several cumulant slices of the observations process. 2. SIGNAL MODEL As shown in figure 1, the chosen signal model for the observations x(t) = [x1(t),··· , xM(t)]Tobey the following equation x(t) = As(t) + n(t)(1) where s(t) = [s1(t),··· , sN(t)]Tis the signal vector process of Nindependent components, n(t)the noise vector process, and Ais a M×Nmixing matrix. We consider the following assumptions: A1 The components of s(t)are mutually independent, locally stationary and normalized to zero mean and unit variance. A2 The noise vector process n(t)is locally stationary, Gaussian, white (Rn(t2, t1) = δ(t2−t1)E[n(t)(n(t))H]) and with a known correlation matrix Rn(t, t) = 463 4th International Symposium on Independent Component Analysis and Blind Signal Separation (ICA2003), April 2003, Nara, Japan Axy U N + n MNP W sz G Fig. 1. Signal model for the blind extraction of Psources. E[n(t)(n(t))H]or one that can be accurately estimated from the observations, perhaps using factor analysis or any robust prewhitening technique. A3 The mixing matrix Ais full-column rank. A4 For a given subset of Pindependent components {s1(t),...,sP(t)}, that one wishes to extract, there exist a set Ω = {τk= (t(k) 1,...,t(k) q), k = 1,...,r : τk∈Rqif q > 2, τk∈R2\ {(t, t),∀t∈R}if q= 2}and a permutation σof the indices 1,...,P that sorts the following statistic of the components ψΩ(sj) = X τ∈Ω wτ|Cum(sj(t1),··· , sj(tq))|2 in such a way that these inequalities hold true ψΩ(sσ(1))>···> ψΩ(sσ(P))> ψΩ(sP+1)≥ · · · ≥ ψΩ(sN)(2) From A1 one obtains AAH=Rx(t, t)−Rn(t, t). From A2 and A3, using principal component analysis one can project the observations onto the signal subspace to reduce their dimensionality from Mto Nand also sphere the resulting signals. The N×Mprewhitenning system W= (AAH)−1/2gives the vector of preprocessed observations z(t) = Wx(t)(3) In order to extract Pindependent components (1≤P≤ N) from the mixture we use a P×Nmatrix Uwhich is semi-unitary (UUH=IP). This matrix multiplies the preprocessed observations to give the outputs or estimated sources y(t) = [y1(t),··· , yP(t)]Tas y(t) = Uz(t) = Gs(t) + UWn(t)(4) where G=UWA denotes the global transfer function from the independent components to the outputs. Assumption A4 guarantees that there exist and order relation among the sources that is maximized by those we want to extract. 3. EXTRACTION OF A SINGLE SOURCE We will first analyze the case of P= 1, i.e., the extraction of a single source. Then Uand Gare row vectors both of unit 2-norm and there is a single output y(t) = Uz(t). We propose to estimate the desired independent component by jointly maximizing a weighted square sum of cumulants of fixed1order q≥2, determined by the tuples τ= (t1,··· , tq)contained the set Ω. A contrast function that achieves this objective is given by ψΩ(y) = X τ∈Ω wτ|Cum (y(t1),··· , y(tq))|2 subject to kUk2= 1.(5) were wτare positive weighting terms. The chosen notation t1,...,tqin (5) is due to the fact that from A1 and A2 the observation process may be large term non-stationary, case for which the contrast can exploit the same cumulant slices at different times (using segments of quasi-stationary data). The problem with this approach is in the difficulty of the optimization of (5), which is highly non-linear with respect to U. The following theorem, whose proof is sketched in the Appendix, shows how to circumvent this difficulty by proposing a similar contrast function to (5) but whose dependence with respect to each of the extracting system candidates is quadratic, and thus, much more easy to optimize using algebraic methods. Theorem 1 Under the assumptions A1-A4, particularized for the extraction of only one independent component s1, there exist a set Ωfor which ψΩ(s1)> ψΩ(sj)∀j= 2,...,N. (6) Considering a set of qcandidates for the extracting system {U[1],...,U[q]}and the set of their respective outputs ¯y= {y[1],...,y[q]}, the following multivariate function ψΩ(¯y) = X τ∈Ω wτ  Cum(y[1](t1),··· , y[q](tq))   2 subject to kU[m]k2= 1, m = 1,...,q.(7) where wτ>0, is a contrast function whose global maximum leads to the extraction of the desired source, i.e., at this extreme point y[1](t) = ···=y[q](t) = s1(t). We should note that this contrast admits a least squares matching interpretation associated with the rank one approximation of cumulant tensors [11]. A good method to maximize the proposed contrast ψΩ(¯y)is to optimize it cyclically with respect to each one of the elements U[m], m = 1,...,q, while keeping fixed the others. In the following, the superindex ()[k]will continue denoting the k-th variable (cyclic notation) while the superindex ()(k)will indicate the variable taken value at the k-th iteration (sequential notation). Then, note that at the (k)th iteration the [(k 1The results of the paper also apply for an arbitrary combination of cumulants with different orders but, due to the somewhat more cumbersome notation it needs, this extension will be provided elsewhere. 464 mod q)+1] variable will be optimized. Due to the invariant property of ψΩ(¯y)with respect to permutations in its arguments, the cyclic maximization of the contrast is equivalent to the sequential maximization of the function φΩ(U(k)) = X τ∈Ω wτ|Cum(y(k)(t1), y(k−1)(t2),··· ...,y(k−q+1)(tq))|2 =U(k)M(k−1)U(k)H(8) with respect to the extraction system U(k)through iterations. Note that M(k−1) is a constant matrix (as long as U(k−1),··· ,U(k−q+1) are kept fixed) given by M(k−1) =X τ∈Ω wτc(k−1) zy(τ)c(k−1) zy(τ)H (9) c(k−1) zy(τ) = Cum(z(t1), y(k−1)(t2),··· , y(k−q+1)(tq)) Since the vector U(k)Hwhich maximizes φ(U(k))is the eigenvector associated to the dominant eigenvalue of M(k−1), we can move the extracting vector U(k−1) towards the solution with one or more iterations of any of the standard eigenpair finding methods. Starting from the previous solution, if one considers to use Literations of the power method to approximate the dominant eigenvector (in practice L= 1 works well), the following extraction algorithm is obtained U(0) =U(k−1) FOR l= 1 : L U(l)=Pτ∈Ωwτd(l−1) y(τ)c(k−1) zy(τ)H    Pτ∈Ωwτd(l−1) y(τ)c(k−1) zy(τ)H   2 (10) END U(k)=U(L) where d(l−1) y(τ) = U(l−1)c(k−1) zy (τ). 3.1. Convergence analysis The iterative optimization of the function φΩ(·)with respect to U(k)through iterations will result in a monotonous ascent sequence φΩ(U(0))≤φΩ(U(1))≤... ≤φΩ(U(k)) that maximizes φΩ(·). However, this property by itself does not guarantee the convergence to an extraction solution because deceptive local maxima of φΩ(U(k))might exist. The following theorem (whose proof is sketched in appendix B) shows that this is not the case. Theorem 2 Under assumptions A1-A4, the only local maxima of φΩ(U(k))correspond with solutions that extract one of the independent components of the mixture. 4. EXTRACTION OF SEVERAL SOURCES In the case of the extraction of Pindependent components (1≤P≤N) the P×Nmatrix U(k)is semi-unitary with rows U(k) i:, i = 1,...,P. A result proved in [7] states that any non-negative contrast ψΩ(y[1],...,y[q])designed for the extraction of a single source, that satisfies (2) and (20), can be used to construct a contrast function for the extraction of the independent components s1(t),··· , sP(t). Particularized to our case, this leads to the sequential optimization, of ΦΩ(U(k)) = P X i=1 φΩ(U(k) i:)s.t. U(k)(U(k))H=IP.(11) In table 1 we consider two choices for the optimization of the previous function which lead to the thin ICA algorithms (TICA). 4.1. Hierarchical extraction A first option (see step 5b, 1st choice) is to hierarchically maximize (11) with respect to the rows U(k) i:, i = 1,...,P in such a way that each ith row satisfies the constraints U(k) i:(U(k) j:)H=δij ∀j≤i, i.e., the first rows are less constrained than the last ones. Using Householder reflections the new update is expressed in a very compact and simple form, it is U(k)=QHwhere Qis a tall semi-unitary matrix of dimension N×Pwhich results from the thin QR decomposition2of the weighted statistic C(k−1) zy defined in step 5a of table 1. 4.2. Simultaneous extraction A second option (see step 5b, 2nd choice) is to simultaneously maximize ΦΩ(U(k))with respect to all the rows of U(k). After defining the Hermitian matrix of multipliers Λ, the gradient of the Lagrangian function associated to (11) is ∇LΩ(U(k)) = ∇ΦΩ(U(k))−ΛU(k)(12) Noting that ∇ΦΩ(U(k))weakly depends with U(k), only through a set of the diagonal terms, and approximating these diagonal terms by their current estimates D(k−1) y(τ)we obtain ∇ΦΩ(U(k))≈(C(k−1) zy )H. The solution of the equations (C(k−1) zy )H=ΛU(k)(13) U(k)(U(k))H=IP(14) ΛH=Λ(15) 2The thin QR decomposition and the thin Singular Value Decomposition both have, for P << N, a computational complexity of O(NP 2) flops. An efficient implementation of them can be found under the MatLab commands qr(·,0) and svd(·,0). 465 Table 1. Summary of the thin ICA algorithms. 1. Set P≤Nthe number of independent components to extract from x(t). 2. Prewhitening z(t) = Wx(t) 3. Initialization U(0) =IP×N;y(0)(t) = U(0)z(t); k= 1; 4. Estimate ∀τ∈Ωthe matrices C(k−1) zy (τ) = [c(k−1) zy1(τ),...,c(k−1) zyq(τ)] where c(k−1) zyi(τ)is defined in (18) and initialize U(0) =U(k−1) 5. FOR l= 1 : L (a) Compute the diagonal matrices D(l−1) y(τ) = diag U(l−1)C(k−1) zy (τ) and the weighted sum: C(l−1) zy =X τ∈Ω wτC(k−1) zy (τ)(D(l−1) y(τ))∗ (b) Choice 1. OR 2.                  1. Hierarchical approach: [Q,R] = qr(C(l−1) zy ,0) U(l)=QH 2. Simultaneous approach: [Q,∆P×P,V] = svd(C(l−1) zy ,0) U(l)=Vsign(∆P×P)QH END 6. Update U(k)=U(L)and estimate Pindependent components: y(k)(t) = U(k)z(t) 7. IF Convergence STOP ELSE k=k+ 1; RETURN TO 4 that maximizes ΦΩ(U(k))is [Q,∆P×P,V] = svd(C(k−1) zy ,0) U(k)=Vsign(∆P×P)QH(16) where svd(·,0) is the MatLab command for the thin singular value decomposition. The validity of the approximation can be further improved at iteration k, by considering Lsubiterations (l= 1,...,L) of a zigzag procedure which updates the estimates D(l−1) y(τ)of the diagonal terms as the solution changes. This is done in step 5. of table 1. 4.3. Projection onto the symmetric subspace From theorem 1 one a priori knows that the solutions U[m], m= 1,...,q which extract the independent components belong to the symmetric subspace U[1] =...=U[q]. This information can be exploited to improve the convergence of TICA algorithms just adding, at the end of each iteration, a projection step of U(k),...,U(k−q+1) onto this subspace U(k),U(k−1),...,U(k−q+1) →U(k),...,U(k)(17) This is obtained using definition (b) instead of (a) c(k−1) zyi(τ) =            a) Without projection: Cum(z(t1), y(k−1) i(t2),··· , y(k−q+1) i(tq)) b) With projection: Cum(z(t1), y(k−1) i(t2),··· , y(k−1) i(tq)) (18) 5. SIMULATIONS In this section we illustrate how the algorithm (in similarly with [9]) can be applied to obtain accurate estimates from reduced set of observations. In our example an array of 20 sensors registers 250 snapshots of the observations. These are a random instantaneous mixture of 10 independent signals, in presence of white additive Gaussian noise, and with a maximum signal to noise ratio of 15dB. The desired independent components are the three correlated signals that can be obtained, after normalization, from the filtering of three binary processes by the corresponding systems F1(z−1) = (1+0.4z−1+0.9z−2+.5z−3)−1,F2(z−1) = (1+0.6z−1− 0.3z−2)−1and F3(z−1) = (1 −0.7z−1)−1. The other seven independent components are samples of temporally i.i.d. uniform processes. We chose second order statistics q= 2 because for short data records like this they are usually the most reliable, and we set Ω = {(t1, t1−1),(t1, t1− 2),··· ,(t1, t1−7)}because these seven pairs guarantee that the considered independent components can be ordered according to (2). We run the algorithm in one hundred random experiments (with different random matrices, sources and noise samples). In each experiment we applied the simultaneous TICA algorithm, with P= 3 and L= 2, which extracted in all the cases the desired subset sources. As can be observed in figure 2 the convergence of the TICA algorithm is quite fast, requiring between 3and 6iterations. Unfortunately, no proper comparison can be given with other algorithms because SOBI and JADE cannot extract a subset of signals while Fast-ICA did not perform well for such a small data set. 6. DISCUSSION Recently, we notice that the contrast function proposed in theorem 1 admits a least squares cumulant matching inter466 pretation associated with the rank one approximation of a set of cumulant tensors. A closely related contrast and an alternating Least Squares technique similar to that we use were previously proposed in [11] to solve the Blind Source Separation problem. The results of this paper bring a complementary insight for the simultaneous extraction problem and lead to the proposal of the thin ICA algorithms. In sequel, we comment some of the interesting links of these algorithms with other existing approaches: •For P= 1,q= 4 and Ω = {(t, t, t, t)}the TICA algorithms with projection extract a single independent component by maximizing the modulo of the kurtosis of the output. In this case, both TICA algorithms particularize to the fixed point algorithm (Fast-ICA with cubic non-linearity) proposed in [5]. For arbitrary P≤Nthe thin ICA implementation based on the thin QR decomposition is equivalent to the hierarchical application of the fixed point algorithm with deflation. Similarly, the thin ICA implementation based on the SVD reduces to the fixed point algorithm with symmetric orthogonalization when P= N(BSS), and provides a novel extension of it for P < N (BSE). •For P=Nand q= 2 (alternatively q= 4) the TICA algorithms extract all the independent components using second order statistics (fourth order statistics) and, when it is possible, using also any nonstationarity of the independent components. The criterion (7) is equivalent to that of the SOBI [9] (JADE [8]) algorithm based on the joint approximate diagonalization of a certain set of cumulant slices, although, the implementation differs. •For 1< P < N and arbitrary qthe TICA algorithms extract Pindependent components using the joint optimization criteria (11). In this case, none of the previously cited algorithms can solve this problem: fixed point algorithms (Fast-ICA) does not perform a joint optimization of several statistics, while SOBI and JADE implementations are not suitable for extraction because the extended Jacobi plane rotations they use are not the most adequate technique for the estimation of a subset of eigenvectors. 7. CONCLUSIONS We have proposed a multivariate contrast function for the extraction of a subset of desired independent components from a linear mixture. This contrast function jointly optimizes several statistics of the same order and have no spurious maxima. We have suggested the thin ICA algorithms for the maximization of the contrast because they combine, 0 5 10 15 20 25 10−4 10−2 100 ITERATIONS PERFORMANCE INDEX 0 1 2 3 4 5 6 7 8 9 10 11 0 0.5 1 1.5 SAMPLE EXPERIMENT COEFFICIENTS OF |G| Fig. 2. Upper fig.: performance index Pindex(G) = (PN)−1PP i=1 kGi:k2 2/kGi:k2 ∞−1versus iterations. Continuous line is the median curve of convergence for 100 experiments, the dashed lines denote the 5th and 95th percentiles. Lower figure presents the coefficients of the different rows of a 3×10 matrix Gfor one sample experiment. at the same time, several of the advantages of other powerful techniques like Fast-ICA, JADE and SOBI. A. PROOF OF THEOREM 1 The proof starts observing that under assumptions A1-A4 the contrast (7) is theoretically unaffected by the additive noise. Then, using Cauchy-Schwarz’s inequality one obtains that each square cumulant within the summation is upper bounded by   Cum(y[1](t1),··· , y[q](t−τq))   2 ≤ N X j=1 |G[p] 1j|2|Cum(sj(t1),··· , sj(tq))|2 ·  N X j=1 Y m6=p |G[m] 1j|2  Since the global transfer vectors are normalized N X j=1 q Y m6=p |G[m] 1j|2≤ q Y m6=p N X j=1 |G[m] 1j|2= 1 (19) therefore,   Cum(y[1](t1),··· , y[q](tq))   2≤ N X j=1 |G[p] 1j|2 ·|Cum(sj(t1),··· , sj(tq)|2 467 Substituting these terms in (7) results in ψΩ(¯y)≤ N X j=1 |G[p] 1j|2ψ(sj)p= 1,...,q. (20) But from the ordering ψΩ(s1)> ψΩ(sj)∀j= 2,...,N in A4, and the constraint kG[p]k2= 1 we finally obtain ψΩ(y[1],...,y[q])≤ψΩ(s1),(21) which means that the upper bound coincides with the extraction of the desired independent component. Noting that the equality between both sides of (19) only holds when the row vectors are equal G[1] =··· =G[q]and aligned with one of the axis ejT, j = 1,··· , N, one can conclude that the global maximum of ψΩ(¯y)is only attained at the extraction of the desired source, i.e., when U[1] =··· =U[q]= (WAe1)Hwhere e1= (1,0,...,0)T. B. PROOF OF THEOREM 2 Any critical point U0of φΩ(·)has associated a global extraction system G0=U0WA. Following [12] we define the set of indices for which the elements of G0are nonzero I={m:G0 1m6= 0, m = 1,...,N},(22) From the proof of theorem 1 one can observe that any solution which extracts one of the independent components (Ihas cardinality one) is a local maximum of the contrast. Thus, the deceptive local maximum of φΩ(·), if they exit, should correspond with G0having at least two nonzero elements, however, these kind of points cannot be a maximum of the function φ(·)because we can always find local perturbation of them for which the function increases. For a local perturbation α∈RNwhere: kαksufficient small, PN j=1 αj= 0,αj= 0 ⇔j6∈ I, and such that |G1j|2=|G0 1j|2+αj, j = 1,...,N;(23) the function φ(·), at the perturbed point U, is written as φ(U) = (1 + γ(α))φ(U0) + β(α) + o(kαk2) where γ(α) = q(q−2) 4X j∈I      αj G0 1j      2 (24) β(α) = X τ∈Ω wτ      q 2X j∈I αj (G0 1j)q |G0 1j|2csj(τ)      2 (25) Since γ(α)≥0and β(α)>0for all q≥2we have that φ(U)> φ(U0),(26) and we conclude that any solution for which Ihas cardinality greater than 1 cannot be a local maximum of φΩ(·). C. REFERENCES [1] P. Comon, “Independent component analysis, a new concept?,” Signal Proc., vol. 3(36), pp. 287–314, 1994. [2] J. F. Cardoso, “Blind signal separation: Statistical principles,” Proceedings of the IEEE, vol. 86(10), pp. 2009– 2025, 1998. [3] A. Cichocki, S.-i. Amari, Adaptive Blind Signal and Image Processing, John Wiley & Sons, 2002. [4] N. Delfosse and P. Loubaton, “Adaptive blind separation of independent sources: A deflation approach,” Signal Processing, 1995, vol. 45, pp. 59–83. [5] A. Hyvarinen and E. Oja, “A fast fixed-point algorithm for independent component analysis,” Neural Computation, 1997, vol. 9, pp. 1483–1492. [6] S. Amari, “Natural gradient learning for overand under-complete bases in ICA,” Neural Computation, vol. 11, pp. 1875–1883, 1999. [7] S. Cruces, A. Cichocki, S-i. Amari, “On a new blind signal extraction algorithm: different criteria and stability analysis”, IEEE Signal Proc. Letters, vol. 9(8), pp. 233–236, 2002. [8] J.-F. Cardoso and A. Solumiac, “Blind beamforming for non Gaussian signals,” IEE Proceedings-F, vol. 140(6), pp. 362–370, 1993. [9] A. Belouchrani, K. Abel-Meraim, J.-F. Cardoso, and E. Moulines, “A blind source separation technique using second-order statistics,” IEEE Trans. on Signal Processing, vol. 45(2), pp. 434–444, 1997. [10] L. De Lathauwer, B. De-Moor, and J. Vandewalle, “Independent component analysis and (simultaneous) third-order tensor diagonalization”, IEEE Trans. on Signal Processing, vol. 49(10), pp. 2262-2271, 2001. [11] L. De Lathauwer, P. Comon, B. De-Moor, and J. Vandewalle, “Higher-order power method - application in independent component analysis,” in International Symposium on Nonlinear Theory and Applications NOLTA, Dec. 1995. [12] O. Shalvi and E. Weinstein, “New criteria for blind deconvolution,” IEEE Trans. on Information Theory, vol. 36(2), pp. 312–321, 1990. 468