Comput Optim Appl DOI 10.1007/s10589-007-9031-2 Algorithms with high order convergence speed for blind source extraction Pando Georgiev · Panos Pardalos · Andrzej Cichocki © Springer Science+Business Media, LLC 2007 Abstract A rigorous convergence analysis for the fixed point ICA algorithm of Hyvärinen and Oja is provided and a generalization of it involving cumulants of an arbitrary order is presented. We consider a specific optimization problem OP(p), p > 3, integer, arising from a Blind Source Extraction problem (BSE) and prove that every local maximum of OP(p) is a solution of (BSE) in sense that it extracts one source signal from a linear mixture of unknown statistically independent signals. An algorithm for solving OP(p) is constructed, which has a rate of convergence p − 1. Keywords Independent component analysis · Cumulants · Fixed point algorithm 1 Introduction The problem of blind source extraction BSE, which we shall consider is formulated as follows: we can observe sensor signals (random variables) x(k) = Research partially supported by NSF and NIH grants. P. Georgiev () ECECS Department, University of Cincinnati, ML 0030, Cincinnati, OH 45221-0030, USA e-mail: pgeorgie@ececs.uc.edu P. Georgiev Sofia University “St. Kl. Ohridski”, 5 J. Bourchier Blvd., 1164, Sofia, Bulgaria e-mail: pandogeorgiev@gmail.com P. Pardalos Center for Applied Optimization, University of Florida, Florida, USA e-mail: pardalos@cao.ise.ufl.edu A. Cichocki The Institute for Physical and Chemical Research (RIKEN), Brain Science Institute, Wako-shi, Saitama 351-0198, Japan e-mail: cia@bsp.brain.riken.jp P. Georgiev et al. [x1 (k), . . . , xm (k)]T which are described as x(k) = As(k), k = 1, 2, . . . , (1) where A is n × n non-singular unknown mixing matrix, s(k) = [s1 (k), . . . , sn (k)]T is a vector of unknown zero mean source signals. Our objective is to estimate the source signals sequentially one-by-one or simultaneously assuming that they are statistically independent. The literature about BSE problem (resp. independent component analysis, ICA) is huge (see for instance [2, 3, 5] and references therein). In this article we prove rigorously the convergence of the algorithm of Hyvarinen and Oja [6] under suitable conditions and generalize it to the case of cumulants of arbitrary order (grater than two). It is an example of an algorithm with high order convergence speed (bigger than cubic). Define the function ϕ : Rn → R by ϕp (w) = cump (wT x) where cump means the self-cumulant of order p: cump (s) = cumulant(s, . . . , s )    p (see [8] for definition and properties of the cumulants). We should mention the following property of the cumulants, which is used essentially: if si , i = 1, . . . , n, are statistically independent, then  n  n   p ci cump (si ). (2) cump ci si = i=1 i=1 Consider the maximization problems OP(p) maximize |ϕp (w)| under constraint w = 1 and DP(p) maximize |ψp (c)| under constraint c = 1 where ψp (c) = cump ( ni=1 ci si ). Denoting y = wT x, c = AT w we have T y=c s= n  ci si and ϕp (w) = ψp (AT w). (3) i=1 Without loss of generality we may assume that the matrix A is orthogonal (assuming that we have made the well know preprocessing called “prewhitening”, see [3] or [5]). Algorithms with high order convergence speed for blind source It is easy to see (using (3) and the orthogonality of A) that the problems DP(p) and OP(p) are equivalent in sense that w0 is a solution of OP(p) if and only if c0 = AT w0 is a solution of DP(p). We shall call this property “duality” between the original (observable) problem OP(p) and the dual (hidden) problem DP(p). A very useful observation is the following: if a vector c contains only one nonzero component, say ci0 = ±1, then the vector w = Ac gives extraction (say y(k)) of the source with index i0 , since y(k) := wT x(k) = cT AT x(k) = cT s(k) = si0 (k) ∀k = 1, 2, . . . . (4) We shall show that the solutions c of DP(p) have exactly one nonzero element, and by the above mentioned duality, we can obtain the vectors w = Ac as solutions of the original problem OP(p), and by (4) we achieve extraction. One interesting property of the optimization problem OP(p) is that it has exactly n solutions (up to sign) which are orthonormal and any of which gives extraction of one source signal. The fixed point algorithm finds one of its solutions. 2 Extraction via maximization of the absolute value of the cumulants We note that the idea of maximizing of cum4 (wT x) in order to extract one source from a linear mixture is already considered in [4]. We need the following lemma, which is generalization (with a simple proof) of a lemma in [4] (the case considered there is p = 4). Lemma 1 Consider the optimization problem: minimize (maximize) n  p ki v i i=1 subject to v = c > 0 where p > 2 and v = (v1 , . . . , vn ). Denote I + = {i ∈ {1, . . . , n}: ki > 0}, I − = {i ∈ {1, . . . , n}: ki < 0} and ei = (0, . . . , 0, 1, 0, . . . , 0), (1 is at the ith place). (a) If p is even, then the points of local minimum are exactly the vectors m± i = ±cei , ± − i ∈ I and the points of local maximum are exactly the vectors Mj = ±cej , j ∈ I +. (b) If p is odd, then the points of local minimum are exactly the vectors m− i = cei , + − + i ∈ I and mi = −cei , i ∈ I and the points of local maximum are exactly the − + − vectors M+ j = cej , j ∈ I and Mj = −cej , j ∈ I . P. Georgiev et al. Proof (a) Applying the Lagrange multipliers theorem for a point of a local optimum v = (v 1 , . . . , v m ), we write p−1 ki pv i − 2λv i = 0, i = 1, . . . , m, (5) where λ is a Lagrange multiplier. Multiplying (5) by v i and summing, we obtain pfopt = 2λc2 , where fopt means the value of f at the local optimum. Hence λ= p fopt . 2c2 (6) From (5) we obtain p−2 v i ki pv i − p fopt = 0, c2 or fopt ± ki c 2 whence v i is either 0, 1 p−2 . (7) Case 1. Assume that ki0 < 0 for some index i0 and v is a local minimum. Then obviously floc min < 0. According to the second order optimality condition [1], a point x0 is a local minimum if hT L′′ (x0 )h > 0 ∀h ∈ K(x0 ) = {h: hT x0 = 0}, h = 0, where L(x) = n  p ki xi − λ x2 − c2 i=1 is the Lagrange function. In our case, by (6) and (7) we obtain hT L′′ (v)h = n  p−2 p(p − 1)ki v i − 2λ h2i i=1     p 2 hi − h2i = 2 floc min (p − 2) c i∈I (8) i ∈I / where I is the set of those indexes i, for which v i is different from 0. We shall check the second order sufficient condition for a local minimum for the points m± i0 . We have K(m± i0 ) = {h: hi0 = 0}, Algorithms with high order convergence speed for blind source therefore, for h ∈ K(m± i0 ), h = 0 we have hT L′′ (m± i0 )h > 0, since hi0 = 0 and floc min < 0, i.e. the second order sufficient condition is satisfied and m± is a local minimum. By (4) it follows that for any vector v with at least two nonzeros the quadratic form (8) can take positive and negative values for different values of h, i.e. the necessary condition for a local minimum is not satisfied for such a vector. Case 2. Assume that kj > 0 and v is a local maximum. We apply Case 1 to the function −f and finish the proof. The assertion (b) is proved analogously.  In the following theorem we show that any local solution of OP(p) gives extraction of one source, so we don’t care about what type of solution we obtain: local or global. We characterize all solutions of the Lagrange equation for OP(p). Recall that the Lagrange equation of OP(p) is ϕp′ (w) sign(ϕp (w)) + 2λw = 0. (9) Theorem 2 Suppose that the matrix A in (1) is orthogonal. Then: (a) all local solutions of maximization problem OP(p) are exactly ±Aci , i = 1, . . . , n; (b) Any solution w of the Lagrange equation for OP(p) is characterized by the equations p−1 ci ψ ′ (c) = pcump (si )ci ∀i = 1, . . . , n, (10) where c = AT w. (c) We can extract the original source signals up to sign and permutation: yi (k) = wTi x(k) = ±sri (k), ri ∈ {1, . . . , n}. Proof (a) By property (2) and by Lemma 1 we have that DP(p) has exactly n solutions (up to sign), which are the vectors ±ei having ±1 in ith place and 0 in the other places. Now the conclusion of (a) is evident, since the problems OP(p) and DP(p) are equivalent (see the duality in the introduction). (b) Let w be a solution of (9). Since ϕp (w) = ψp (AT w), using the property of the cumulants (2) and the chain rule for differentiating the composite function ψp (AT w) with respect to w, we obtain: p−1 ϕp′ (w) = Aψp′ (c) = pA c1 p−1 cump (s1 ), . . . , cn cump (sn ) T . (11) Hence, using (9) we have Ac = w = ± since A is orthogonal. ϕp′ (w) ϕp′ (w) =± Aψp′ (c) Aψp′ (c) =± Aψ ′ (c) , ψ ′ (c) (12) P. Georgiev et al. Multiplying (12) by AT we obtain c=± ψ ′ (c) , ψ ′ (c) (13) or, componentwise, p−1 ci = ± pcump (si )ci ψ ′ (c) ∀i = 1, . . . , n, (14) which is (10). Analogously we obtain that any solution of (10) is a solution of (9).  The assertion (c) is already proved by (4). 3 A generalization of the fixed point algorithm Consider the following algorithm: w(l) = ϕp′ (w(l − 1)) ϕp′ (w(l − 1)) , l = 1, 2, . . . , (15) which is a generalization of the fixed point algorithm of Hyvärinen and Oja. The name is derived by the Lagrange equation (9), since (15) tries to find a solution of it iteratively, and this solution is a fixed point of the operator defined by the right-hand side of (15). Theorem 1 Assume that si are statistically independent, zero mean signals and the mixing matrix A is orthogonal. Let p ≥ 3 be a given integer number, cump (si ) = 0, i = 1, . . . , n and let  1  I (c) = arg max ci cump (si ) p−2 . 1≤i≤n Denote by W0 the set of all elements w ∈ Rn such that w = 1, the set I (AT w) contains only one element, say i(w), and ci(w) = 0. Then: (a) The complement of W0 has a measure zero. (b) If w(0) ∈ W0 then lim yl (k) = ±si0 (k) l→∞ ∀k = 1, 2, . . . , where yl (k) = w(l)T x(k) and i0 = i(w(0)). (c) The rate of convergence in (b) is of order p − 1. Proof (a) It is easy to see that the complement K of W0 is a finite union of proper subspaces of Rn , therefore this union has a measure zero. Indeed, K can be represented as n  K= Wi,j , j =i,j =1 Algorithms with high order convergence speed for blind source where Wi,j is a subspace such that I (AT w) contains at least the indices i and j . (b), (c) Using (11) with w = w(l) and c(l) = Aw(l) we obtain (like (12)): Ac(l) = w(l) = ϕp′ (w(l − 1)) ϕp′ (w(l − 1)) Aψp′ (c(l − 1)) = = Aψ ′ (c(l − 1)) , (16) ψ ′ (c(l − 1)) ∀l = 1, 2, . . . , (17) Aψp′ (c(l − 1)) since A is orthogonal. Multiplying (16) by AT we obtain c(l) = ψ ′ (c(l − 1)) ψ ′ (c(l − 1)) or, componentwise, ci (l) = σi p|cump (si )|ci (l − 1)p−1 ψ ′ (c(l − 1)) (18) for every i = 1, . . . , n, l = 1, 2, . . . , where σi = sign(cump (si )). From (18) we obtain  1  ci (l)cump (si ) p−2 1 σi [ci (l − 1)|cump (si )| p−2 ]p−1 . = (c1 (l − 1)p−1 cump (s1 ), . . . , cn (l − 1)p−1 cump (sn ))T  (19) From (19) it follows by induction that if the initial conditions satisfy ci (0)cump (si ) = 0 for some i ∈ {1, . . . , n}, then the denominator in (19) can not be zero for every l ≥ 0, so w(l) in (16) is well defined. Hence, if w(0) ∈ W0 , we obtain  1 1 p−1 σi ci (l − 1) |cump (si )| p−2 ci (l) |cump (si )| p−2 . (20) = ci0 (l) |cump (si0 )| σi0 ci0 (l − 1) |cump (si0 )| Now by (20) one proves easily by induction that for every l ≥ 0 and every i = 1, . . . , n ci (l) |cump (si )| ci0 (l) |cump (si0 )| 1 p−2 = σi σi0 l  ci (0) |cump (si )| ci0 (0) |cump (si0 )| 1 p−2 (p−1)l . (21) From (21) it follows that ci (l) → 0 when l → ∞ for i = i0 and |ci0 (l)| → 1 (since c(l) = 1), as the speed of the convergence is p − 1. Thus we proved: lim yl (k) = lim wT (l)x(k) = lim cT (l)s(k) l→∞ l→∞ = ±si0 (k) l→∞ ∀k = 1, 2, . . . . (22) Similar algorithm can be used for cumulants involving time delays, like in [7]. This will allow us to extract signals with zero cumulants.  P. Georgiev et al. 4 Examples Below we consider some examples, for concrete values of p. 1. p = 3. Then   ϕ3 (w) = cum3 (wT x) = E (wT x)3 and 2. p = 4. Then   ϕ3′ (w) = 3E (wT x)2 x .     ϕ4 (w) = cum4 (wT x) = E (wT x)4 − 3 E (wT x)2 and 2       ϕ4′ (w) = 4E (wT x)3 x − 12E (wT x)2 E (xxT ) w. We note that if the standard prewhitening is made (i.e. E{xxT } = In , A is orthogonal), the algorithm (15) recovers the fixed-point algorithm of Hyvarinen and Oja, i.e. w(l + 1) = E{(w(l)T x)3 x} − 3w . E{(w(l)T x)3 x} − 3w 3. p = 5. Then       ϕ5 (w) = E (wT x)5 − 10E (wT x)3 E (wT x)2 and         ϕ5′ (w) = 5E (wT x)4 x − 30E (wT x)2 x E (wT x)2 − 20E (wT x)3 E{xxT }w. If, in addition E{xxT } = In , then       ϕ5′ (w) = 5E (wT x)4 x − 30E (wT x)2 x − 20E (wT x)3 w. 4. p = 6. Then       ϕ6 (w) = E (wT x)6 − 15E (wT x)2 E (wT x)4    3  2 − 10 E (wT x)3 + 30 E (wT x)2 ,     ϕ6′ (w) = 6E (wT x)5 x − 30E{xxT }wE (wT x)4         − 60E (wT x)2 E (wT x)3 x − 60E (wT x)3 E (wT x)2 x   2 + 180 E (wT x)2 E{xxT }w. If, in addition, E{xxT } = I, then       ϕ6′ (w) = 6E (wT x)5 x − 30wE (wT x)4 − 60E (wT x)3 x     − 60E (wT x)3 E (wT x)2 x + 180w. Algorithms with high order convergence speed for blind source 5 Some remarks for practical implementation We have implemented the algorithm (15) in the case of six order cumulants, using a deflation procedure. We generated 50 sparse source signals and run the algorithm after random mixture. The number of iterations for extracting all signals one by one using six order cumulants was 313, and the number of iterations using fourth order cumulants was 339. But the computational time is bigger in the case of six order cumulants. This can be explained by the computational time needed for calculation of the six order cumulants, since they have more complex structure. We should mention also that the successful extraction of components depends of their statistical independence. In case when they are not so independent, the extraction is problematic and depends on the initial condition. References 1. Alekseev, V.M., Tihomirov, V.M., Fomin, S.V.: Optimal Control. Nauka, Moscow (1979) (in Russian) 2. Amari, S., Cichocki, A.: Adaptive blind signal processing—neural network approaches. Proc. IEEE 86(10), 2016–2048 (1998) (invited paper) 3. Cichocki, A., Amari, S.: Adaptive Blind Signal and Image Processing. Wiley, Chichester (2002) 4. Delfosse, N., Loubaton, P.: Adaptive blind separation of independent sources: a deflation approach. Signal Process. 45, 59–83 (1995) 5. Hyvarinen, A., Karhunen, J., Oja, E.: Independent Component Analysis. Wiley, New York (2001) 6. Hyvarinen, A., Oja, E.: A fast fixed-point algorithm for independent component analysis. Neural Comput. 9, 1483–1492 (1997) 7. Hyvarinen, A.: Blind source separation by nonstationarity of variance: a cumulant based approach. IEEE Trans. Neural Netw. 12(6), 1471–1474 (2001) 8. Shiryaev, A.N.: Probability. Graduate Texts in Mathematics. Springer, New York (1996) 

Comput Optim Appl DOI 10.1007/s10589-007-9031-2 Algorithms with high order convergence speed for blind source extraction Pando Georgiev · Panos Pardalos · Andrzej Cichocki © Springer Science+Business Media, LLC 2007 Abstract A rigorous convergence analysis for the fixed point ICA algorithm of Hyvärinen and Oja is provided and a generalization of it involving cumulants of an arbitrary order is presented. We consider a specific optimization problem OP(p), p > 3, integer, arising from a Blind Source Extraction problem (BSE) and prove that every local maximum of OP(p) is a solution of (BSE) in sense that it extracts one source signal from a linear mixture of unknown statistically independent signals. An algorithm for solving OP(p) is constructed, which has a rate of convergence p − 1. Keywords Independent component analysis · Cumulants · Fixed point algorithm 1 Introduction The problem of blind source extraction BSE, which we shall consider is formulated as follows: we can observe sensor signals (random variables) x(k) = Research partially supported by NSF and NIH grants. P. Georgiev () ECECS Department, University of Cincinnati, ML 0030, Cincinnati, OH 45221-0030, USA e-mail: pgeorgie@ececs.uc.edu P. Georgiev Sofia University “St. Kl. Ohridski”, 5 J. Bourchier Blvd., 1164, Sofia, Bulgaria e-mail: pandogeorgiev@gmail.com P. Pardalos Center for Applied Optimization, University of Florida, Florida, USA e-mail: pardalos@cao.ise.ufl.edu A. Cichocki The Institute for Physical and Chemical Research (RIKEN), Brain Science Institute, Wako-shi, Saitama 351-0198, Japan e-mail: cia@bsp.brain.riken.jp P. Georgiev et al. [x1 (k), . . . , xm (k)]T which are described as x(k) = As(k), k = 1, 2, . . . , (1) where A is n × n non-singular unknown mixing matrix, s(k) = [s1 (k), . . . , sn (k)]T is a vector of unknown zero mean source signals. Our objective is to estimate the source signals sequentially one-by-one or simultaneously assuming that they are statistically independent. The literature about BSE problem (resp. independent component analysis, ICA) is huge (see for instance [2, 3, 5] and references therein). In this article we prove rigorously the convergence of the algorithm of Hyvarinen and Oja [6] under suitable conditions and generalize it to the case of cumulants of arbitrary order (grater than two). It is an example of an algorithm with high order convergence speed (bigger than cubic). Define the function ϕ : Rn → R by ϕp (w) = cump (wT x) where cump means the self-cumulant of order p: cump (s) = cumulant(s, . . . , s )    p (see [8] for definition and properties of the cumulants). We should mention the following property of the cumulants, which is used essentially: if si , i = 1, . . . , n, are statistically independent, then  n  n   p ci cump (si ). (2) cump ci si = i=1 i=1 Consider the maximization problems OP(p) maximize |ϕp (w)| under constraint w = 1 and DP(p) maximize |ψp (c)| under constraint c = 1 where ψp (c) = cump ( ni=1 ci si ). Denoting y = wT x, c = AT w we have T y=c s= n  ci si and ϕp (w) = ψp (AT w). (3) i=1 Without loss of generality we may assume that the matrix A is orthogonal (assuming that we have made the well know preprocessing called “prewhitening”, see [3] or [5]). Algorithms with high order convergence speed for blind source It is easy to see (using (3) and the orthogonality of A) that the problems DP(p) and OP(p) are equivalent in sense that w0 is a solution of OP(p) if and only if c0 = AT w0 is a solution of DP(p). We shall call this property “duality” between the original (observable) problem OP(p) and the dual (hidden) problem DP(p). A very useful observation is the following: if a vector c contains only one nonzero component, say ci0 = ±1, then the vector w = Ac gives extraction (say y(k)) of the source with index i0 , since y(k) := wT x(k) = cT AT x(k) = cT s(k) = si0 (k) ∀k = 1, 2, . . . . (4) We shall show that the solutions c of DP(p) have exactly one nonzero element, and by the above mentioned duality, we can obtain the vectors w = Ac as solutions of the original problem OP(p), and by (4) we achieve extraction. One interesting property of the optimization problem OP(p) is that it has exactly n solutions (up to sign) which are orthonormal and any of which gives extraction of one source signal. The fixed point algorithm finds one of its solutions. 2 Extraction via maximization of the absolute value of the cumulants We note that the idea of maximizing of cum4 (wT x) in order to extract one source from a linear mixture is already considered in [4]. We need the following lemma, which is generalization (with a simple proof) of a lemma in [4] (the case considered there is p = 4). Lemma 1 Consider the optimization problem: minimize (maximize) n  p ki v i i=1 subject to v = c > 0 where p > 2 and v = (v1 , . . . , vn ). Denote I + = {i ∈ {1, . . . , n}: ki > 0}, I − = {i ∈ {1, . . . , n}: ki < 0} and ei = (0, . . . , 0, 1, 0, . . . , 0), (1 is at the ith place). (a) If p is even, then the points of local minimum are exactly the vectors m± i = ±cei , ± − i ∈ I and the points of local maximum are exactly the vectors Mj = ±cej , j ∈ I +. (b) If p is odd, then the points of local minimum are exactly the vectors m− i = cei , + − + i ∈ I and mi = −cei , i ∈ I and the points of local maximum are exactly the − + − vectors M+ j = cej , j ∈ I and Mj = −cej , j ∈ I . P. Georgiev et al. Proof (a) Applying the Lagrange multipliers theorem for a point of a local optimum v = (v 1 , . . . , v m ), we write p−1 ki pv i − 2λv i = 0, i = 1, . . . , m, (5) where λ is a Lagrange multiplier. Multiplying (5) by v i and summing, we obtain pfopt = 2λc2 , where fopt means the value of f at the local optimum. Hence λ= p fopt . 2c2 (6) From (5) we obtain p−2 v i ki pv i − p fopt = 0, c2 or fopt ± ki c 2 whence v i is either 0, 1 p−2 . (7) Case 1. Assume that ki0 < 0 for some index i0 and v is a local minimum. Then obviously floc min < 0. According to the second order optimality condition [1], a point x0 is a local minimum if hT L′′ (x0 )h > 0 ∀h ∈ K(x0 ) = {h: hT x0 = 0}, h = 0, where L(x) = n  p ki xi − λ x2 − c2 i=1 is the Lagrange function. In our case, by (6) and (7) we obtain hT L′′ (v)h = n  p−2 p(p − 1)ki v i − 2λ h2i i=1     p 2 hi − h2i = 2 floc min (p − 2) c i∈I (8) i ∈I / where I is the set of those indexes i, for which v i is different from 0. We shall check the second order sufficient condition for a local minimum for the points m± i0 . We have K(m± i0 ) = {h: hi0 = 0}, Algorithms with high order convergence speed for blind source therefore, for h ∈ K(m± i0 ), h = 0 we have hT L′′ (m± i0 )h > 0, since hi0 = 0 and floc min < 0, i.e. the second order sufficient condition is satisfied and m± is a local minimum. By (4) it follows that for any vector v with at least two nonzeros the quadratic form (8) can take positive and negative values for different values of h, i.e. the necessary condition for a local minimum is not satisfied for such a vector. Case 2. Assume that kj > 0 and v is a local maximum. We apply Case 1 to the function −f and finish the proof. The assertion (b) is proved analogously.  In the following theorem we show that any local solution of OP(p) gives extraction of one source, so we don’t care about what type of solution we obtain: local or global. We characterize all solutions of the Lagrange equation for OP(p). Recall that the Lagrange equation of OP(p) is ϕp′ (w) sign(ϕp (w)) + 2λw = 0. (9) Theorem 2 Suppose that the matrix A in (1) is orthogonal. Then: (a) all local solutions of maximization problem OP(p) are exactly ±Aci , i = 1, . . . , n; (b) Any solution w of the Lagrange equation for OP(p) is characterized by the equations p−1 ci ψ ′ (c) = pcump (si )ci ∀i = 1, . . . , n, (10) where c = AT w. (c) We can extract the original source signals up to sign and permutation: yi (k) = wTi x(k) = ±sri (k), ri ∈ {1, . . . , n}. Proof (a) By property (2) and by Lemma 1 we have that DP(p) has exactly n solutions (up to sign), which are the vectors ±ei having ±1 in ith place and 0 in the other places. Now the conclusion of (a) is evident, since the problems OP(p) and DP(p) are equivalent (see the duality in the introduction). (b) Let w be a solution of (9). Since ϕp (w) = ψp (AT w), using the property of the cumulants (2) and the chain rule for differentiating the composite function ψp (AT w) with respect to w, we obtain: p−1 ϕp′ (w) = Aψp′ (c) = pA c1 p−1 cump (s1 ), . . . , cn cump (sn ) T . (11) Hence, using (9) we have Ac = w = ± since A is orthogonal. ϕp′ (w) ϕp′ (w) =± Aψp′ (c) Aψp′ (c) =± Aψ ′ (c) , ψ ′ (c) (12) P. Georgiev et al. Multiplying (12) by AT we obtain c=± ψ ′ (c) , ψ ′ (c) (13) or, componentwise, p−1 ci = ± pcump (si )ci ψ ′ (c) ∀i = 1, . . . , n, (14) which is (10). Analogously we obtain that any solution of (10) is a solution of (9).  The assertion (c) is already proved by (4). 3 A generalization of the fixed point algorithm Consider the following algorithm: w(l) = ϕp′ (w(l − 1)) ϕp′ (w(l − 1)) , l = 1, 2, . . . , (15) which is a generalization of the fixed point algorithm of Hyvärinen and Oja. The name is derived by the Lagrange equation (9), since (15) tries to find a solution of it iteratively, and this solution is a fixed point of the operator defined by the right-hand side of (15). Theorem 1 Assume that si are statistically independent, zero mean signals and the mixing matrix A is orthogonal. Let p ≥ 3 be a given integer number, cump (si ) = 0, i = 1, . . . , n and let  1  I (c) = arg max ci cump (si ) p−2 . 1≤i≤n Denote by W0 the set of all elements w ∈ Rn such that w = 1, the set I (AT w) contains only one element, say i(w), and ci(w) = 0. Then: (a) The complement of W0 has a measure zero. (b) If w(0) ∈ W0 then lim yl (k) = ±si0 (k) l→∞ ∀k = 1, 2, . . . , where yl (k) = w(l)T x(k) and i0 = i(w(0)). (c) The rate of convergence in (b) is of order p − 1. Proof (a) It is easy to see that the complement K of W0 is a finite union of proper subspaces of Rn , therefore this union has a measure zero. Indeed, K can be represented as n  K= Wi,j , j =i,j =1 Algorithms with high order convergence speed for blind source where Wi,j is a subspace such that I (AT w) contains at least the indices i and j . (b), (c) Using (11) with w = w(l) and c(l) = Aw(l) we obtain (like (12)): Ac(l) = w(l) = ϕp′ (w(l − 1)) ϕp′ (w(l − 1)) Aψp′ (c(l − 1)) = = Aψ ′ (c(l − 1)) , (16) ψ ′ (c(l − 1)) ∀l = 1, 2, . . . , (17) Aψp′ (c(l − 1)) since A is orthogonal. Multiplying (16) by AT we obtain c(l) = ψ ′ (c(l − 1)) ψ ′ (c(l − 1)) or, componentwise, ci (l) = σi p|cump (si )|ci (l − 1)p−1 ψ ′ (c(l − 1)) (18) for every i = 1, . . . , n, l = 1, 2, . . . , where σi = sign(cump (si )). From (18) we obtain  1  ci (l)cump (si ) p−2 1 σi [ci (l − 1)|cump (si )| p−2 ]p−1 . = (c1 (l − 1)p−1 cump (s1 ), . . . , cn (l − 1)p−1 cump (sn ))T  (19) From (19) it follows by induction that if the initial conditions satisfy ci (0)cump (si ) = 0 for some i ∈ {1, . . . , n}, then the denominator in (19) can not be zero for every l ≥ 0, so w(l) in (16) is well defined. Hence, if w(0) ∈ W0 , we obtain  1 1 p−1 σi ci (l − 1) |cump (si )| p−2 ci (l) |cump (si )| p−2 . (20) = ci0 (l) |cump (si0 )| σi0 ci0 (l − 1) |cump (si0 )| Now by (20) one proves easily by induction that for every l ≥ 0 and every i = 1, . . . , n ci (l) |cump (si )| ci0 (l) |cump (si0 )| 1 p−2 = σi σi0 l  ci (0) |cump (si )| ci0 (0) |cump (si0 )| 1 p−2 (p−1)l . (21) From (21) it follows that ci (l) → 0 when l → ∞ for i = i0 and |ci0 (l)| → 1 (since c(l) = 1), as the speed of the convergence is p − 1. Thus we proved: lim yl (k) = lim wT (l)x(k) = lim cT (l)s(k) l→∞ l→∞ = ±si0 (k) l→∞ ∀k = 1, 2, . . . . (22) Similar algorithm can be used for cumulants involving time delays, like in [7]. This will allow us to extract signals with zero cumulants.  P. Georgiev et al. 4 Examples Below we consider some examples, for concrete values of p. 1. p = 3. Then   ϕ3 (w) = cum3 (wT x) = E (wT x)3 and 2. p = 4. Then   ϕ3′ (w) = 3E (wT x)2 x .     ϕ4 (w) = cum4 (wT x) = E (wT x)4 − 3 E (wT x)2 and 2       ϕ4′ (w) = 4E (wT x)3 x − 12E (wT x)2 E (xxT ) w. We note that if the standard prewhitening is made (i.e. E{xxT } = In , A is orthogonal), the algorithm (15) recovers the fixed-point algorithm of Hyvarinen and Oja, i.e. w(l + 1) = E{(w(l)T x)3 x} − 3w . E{(w(l)T x)3 x} − 3w 3. p = 5. Then       ϕ5 (w) = E (wT x)5 − 10E (wT x)3 E (wT x)2 and         ϕ5′ (w) = 5E (wT x)4 x − 30E (wT x)2 x E (wT x)2 − 20E (wT x)3 E{xxT }w. If, in addition E{xxT } = In , then       ϕ5′ (w) = 5E (wT x)4 x − 30E (wT x)2 x − 20E (wT x)3 w. 4. p = 6. Then       ϕ6 (w) = E (wT x)6 − 15E (wT x)2 E (wT x)4    3  2 − 10 E (wT x)3 + 30 E (wT x)2 ,     ϕ6′ (w) = 6E (wT x)5 x − 30E{xxT }wE (wT x)4         − 60E (wT x)2 E (wT x)3 x − 60E (wT x)3 E (wT x)2 x   2 + 180 E (wT x)2 E{xxT }w. If, in addition, E{xxT } = I, then       ϕ6′ (w) = 6E (wT x)5 x − 30wE (wT x)4 − 60E (wT x)3 x     − 60E (wT x)3 E (wT x)2 x + 180w. Algorithms with high order convergence speed for blind source 5 Some remarks for practical implementation We have implemented the algorithm (15) in the case of six order cumulants, using a deflation procedure. We generated 50 sparse source signals and run the algorithm after random mixture. The number of iterations for extracting all signals one by one using six order cumulants was 313, and the number of iterations using fourth order cumulants was 339. But the computational time is bigger in the case of six order cumulants. This can be explained by the computational time needed for calculation of the six order cumulants, since they have more complex structure. We should mention also that the successful extraction of components depends of their statistical independence. In case when they are not so independent, the extraction is problematic and depends on the initial condition. References 1. Alekseev, V.M., Tihomirov, V.M., Fomin, S.V.: Optimal Control. Nauka, Moscow (1979) (in Russian) 2. Amari, S., Cichocki, A.: Adaptive blind signal processing—neural network approaches. Proc. IEEE 86(10), 2016–2048 (1998) (invited paper) 3. Cichocki, A., Amari, S.: Adaptive Blind Signal and Image Processing. Wiley, Chichester (2002) 4. Delfosse, N., Loubaton, P.: Adaptive blind separation of independent sources: a deflation approach. Signal Process. 45, 59–83 (1995) 5. Hyvarinen, A., Karhunen, J., Oja, E.: Independent Component Analysis. Wiley, New York (2001) 6. Hyvarinen, A., Oja, E.: A fast fixed-point algorithm for independent component analysis. Neural Comput. 9, 1483–1492 (1997) 7. Hyvarinen, A.: Blind source separation by nonstationarity of variance: a cumulant based approach. IEEE Trans. Neural Netw. 12(6), 1471–1474 (2001) 8. Shiryaev, A.N.: Probability. Graduate Texts in Mathematics. Springer, New York (1996)