Blind Source Separation Using Renyi’s Mutual Information	
	
	Kenneth E. Hild II*	Deniz Erdogmus
	Computational NeuroEngineering Lab	Computational NeuroEngineering Lab
	University of Florida	University of Florida
	Gainesville, FL 32611	Gainesville, FL 32611
	k.hild@ieee.org	deniz@grove.ufl.edu
	
	José Príncipe		
	Computational NeuroEngineering Lab
	University of Florida
	Gainesville, FL 32611		
	principe@cnel.ufl.edu
	
Abstract
A blind source separation algorithm is proposed that is based on the minimization of Renyi’s mutual information by means of non-parametric pdf estimation. The two-stage process consists of spatial whitening and a series of Givens rotations, which allows the removal of the joint entropy so that the resulting cost function is the sum of marginal entropies. This formulation avoids the problems encountered by previous MMI (minimum mutual information) methods; namely, the problems of pdf inaccuracy due to truncation of series expansion and the sparsity of data samples needed to estimate joint pdfs in high-dimensional spaces. The performance is compared to an H-J network and Bell and Sejnowski’s ME method, which indicate that the new method is able to find a better solution in far fewer iterations.
1	Introduction
Information theoretic criteria have been applied to the problem of blind source separation using two main methods. They are the ME (maximum entropy) and MMI (minimum mutual information) methods. ME maximizes the joint output entropy and MMI seeks to minimize the mutual information (MI) between the outputs. The mutual information between the N outputs, y, can be written in terms of Shannon’s entropy as,

*www.cnel.ufl.edu.
      			        (1)
which is the sum of the marginal entropies minus the joint entropy, where y is a (N x 1) vector and yi is the ith output. Note that throughout this paper, vectors and matrices are denoted with bold letters and scalars with normal font-weight. Equation 1 represents the canonical contrast function for blind separation in that it “expresses everything required for separation and nothing else” [1]. The difficulty in using MMI is in the estimation of the marginal entropies which are usually derived by approximating the output probability density functions (pdf) with truncated polynomial expansions [2]-[4]. Such approximations, however, are not valid when any of the pdfs in question is far from Gaussianity [2], and because it is an approximation, it is no longer, in general, guaranteed to converge to the same solution as MMI. In order for ME to converge to the MMI solution, on the other hand, component-wise non-linearities are required on the outputs that are equal to the original source cumulative distribution functions (cdf) [3]. Since these are unknown, they must either be approximated during training or assumed to be of a given form. In the first case, the complexity of the algorithm is greatly increased and in the latter case, the solution is no longer guaranteed to be optimal.
A different approach to MMI was taken by Xu et al in [5]. In their work, Renyi’s entropy was used in place of Shannon’s and the MI was estimated by using the Cauchy-Schwartz inequality between the joint pdf and the product of the marginals. Renyi’s entropy, which is given by,
     		        (2)
where fY(y) is the pdf of the output, is a family of entropy definitions that include Shannon’s definition in the limit as a approaches 1. Whereas the MMI methods listed above require an explicit estimate of the (marginal) pdfs, the pdf estimations in [5] are implicit. This is a result of using Renyi’s entropy with a = 2 and Parzen windowing with Gaussian kernels to estimate the pdfs. In this case, the infinite-limit integral in Equation 2 disappears. The resulting equation is simply a double summation over all possible pairs of L points in the output space, where L is a parameter under the control of the user. One very nice feature is that no approximation or truncation is involved for either the joint or marginal entropies (aside from the implicit pdf estimation). One of the drawbacks of the algorithm, however, is that the joint entropy must be estimated. This measurement must be done in an N-dimensional space using L points, and non-parametric pdf estimation, such as Parzen windows, is ill-posed in high-dimensional spaces due to the sparseness of points [5].
As can be seen from the above, the performance of MMI methods suffer from the truncations involved in the estimation of marginal entropies or in the non-parametric estimation of the N-dimensional joint entropy. Likewise, the performance of ME suffers when the non-linearities used are not equal to the source cdfs. However, if Renyi’s entropy formulation were used to estimate only the marginal pdfs, then all of these problems are avoided. Indeed, such an algorithm exists, and is outlined in the following sections.
2	Cost function for Renyi’s mutual information using rotations
Renyi’s version of mutual information is equal to [6],
.		        (3)
Notice that this is not equal to the sum of (Renyi’s) marginal entropies minus the (Renyi’s) joint entropy, which is given as,
.	        (4)
However, both Equations 3 and 4 are non-negative [6] and both evaluate to zero when and only when the joint pdf can be written as a product of the marginals, which occurs identically when the sources have been successfully separated. From this, we know that both Equations 3 and 4 produce the desired solution when they are minimized. Therefore, the goal for source separation using Renyi’s MI can once again be stated as minimizing the sum of (Renyi’s) marginal entropies minus the (Renyi’s) joint entropy. An explanation has already  been given why it is disadvantageous to use Renyi's formulation for the joint entropy, but the joint entropy can be avoided altogether in the following manner. A two-step algorithm is used where the first stage performs spatial whitening and the second stage performs an N-dimensional rotation. Notice that,

	        (5)
where R is a rotation matrix, y = RTx, and det R is the determinant of R, which in this case is 1. This shows that Renyi’s joint entropy is invariant to rotations. Therefore, minimizing Renyi’s MI is equivalent to minimizing the sum of Renyi’s marginal entropies at the output of a rotation matrix whose inputs have been spatially whitened. Hence, the optimization problem reduces to finding,
    			        (6)
This architecture is still completely general and the cost function has the advantage we are seeking. Comon also used the two-step procedure just described, with the difference that the problem was formulated in terms of (Shannon’s) negentropies and the approximations of the negentropy utilized truncated pdf expansions [4]. 
3	System description

The overall system block diagram for the blind source separation using Renyi’s  MMI is given in Figure 1 for a two-input system. The equations can be given in
Figure 1: System block diagram (N = 2).
matrix notation by constructing an (N x M) source, whitened observation, and output matrices, s, x, and y, respectively, where the number of sources and observations equal N and each signal has M samples. These equations are,
				        (7)
				        (8)
where W = FL-1/2, F is the matrix of eigenvectors of the observations (HTs), L is the (corresponding) eigenvalue matrix, E[xxT] = IN due to the spatial whitening, and IN is the (N x N) identity matrix. The rotation matrix, R, is constructed from the product of N(N-1)/2 Givens rotation matrices in the following manner,
				        (9)
where Rij is an (N x N) identity matrix, IN, with four elements modified. The elements located at IN(i,i) and IN(j,j) are changed to cosqij and the elements at IN(i,j) and IN(j,i) are changed to sinqij and –sinqij , respectively, where IN(i,j) is the element of IN located at the ith row and jth column. This method has one rotation angle for each Rij, whose purpose is to perform a rotation in the output space in the plane specified by the ith and jth (orthogonal) basis vectors. For this method, it does not matter the order in which the Rij are multiplied or whether pre-multiplication or post-multiplication is used as long as consistency is maintained.
The gradient of R with respect to qij is also needed for the algorithm development. It is constructed in a manner similar to Rij. In this case, however, the products are taken of N(N-1)/2 matrices Pkl , where Pkl = Rkl for all indices kl ¹ ij, because these are constants as far as qij is concerned. For index ij, however, Pij equals the ÑRij, which is constructed from an (N x N) matrix of zeros, ZN, by replacing ZN(i,i), ZN(i,j), ZN(j,i), and ZN(j,j) with –sinqij, cosqij, -cosqij, -sinqij, respectively. The construction of ÑijR is represented with the following equation,
				      (10)
4	Stochastic gradient descent algorithm
When Parzen windowing is used with a Gaussian kernel, the estimate for Renyi’s quadratic (a = 2) marginal entropy becomes [5],
	      (11)
where G(x,2s2) is a Gaussian pdf, and yi(j) is the jth sample of output yi. Substituting this into Equation 6 and taking the derivative with respect to qij produces,
     			      (12)

where (ÑijR)i is the ith column of ÑijR and x(m) is the (N x 1) vector of the whitened observations. The overall update equation for stochastic gradient descent of the cost function given by Equation 6 is then,
    ,			      (13)
where q(n) and Dq(n) are (N(N-1)/2 x 1) vectors of angles and h is the step size.
5	Results
Figure 2A shows the performance surface for speech data consisting of two male speakers that is instantaneously mixed and where L = 100. There are several important items about this plot. It shows that the solution is periodic with a period of p/2, not 2p. This is expected since there are a total of four combinations of sign changes and permutations of two signals, resulting in a period of 2p/4 = p/2. As can be seen from the figure, the minimum value is found at  –0.5, which is the correct location. Also notice that there are no local minima in this example (within each period of p/2), although this result is not proven analytically. This is unlike the performance surface for ME, which can be expected to have many local maxima due to the non-linearities used.
Figure 2B shows the convergence of the single angle parameter when N = 2 and the mixing matrix imposes a p/4 rotation. The experiment was run ten times, each time by selecting a random L = 100 points from the (whitened) observations. Only ten or fifteen iterations are needed for all but one of the solutions to converge, and eight out of ten found the correct solution with a high degree of accuracy.
In order to compare this algorithm with other methods that use a different criteria, a performance criteria based on the product of the mixing and demixing matrices is used. This criteria, which will be referred to as signal-to-distortion ratio, SDR, is given as,
     		      (14)
Figure 2: (A). Performance surface (lower plot) and derivative (upper plot) for N = 2 sources as a function of the single rotation parameter. (B). Adaptation of angle parameter versus iterations run 10 times with L = 100 different random samples.

where q = RTWTHT , qi is the ith row of q, and max(qi) is the maximum element in the row vector. The numerator expression is the power of the signal in the ith output and the denominator is the distortion power, which is found by subtracting the total power at the ith output minus the signal power. It is desired for this equation to be invariant to permutations and individual source gains, since these are indeterminants in the blind source separation problem. This is in fact the case when it is assumed that the sources have equal power. In case they do not, it is a simple matter to incorporate these factors into Equation 14.
Renyi’s MMI method is compared to two other methods, a H-J network [7] and Bell and Sejnowski’s ME method [8]. The H-J network requires two functions to be defined; these functions were chosen to be f1(x) = x and f2(x) = x3. Both the H-J network and the ME method were trained using a batch size of 200 samples. Likewise, the Renyi’s MMI method used L = 200 points. The three methods are tested with a linear, instantaneous mixture of 10 sources. The 10 sources include speech from five male speakers, four female speakers and one source that is music. The values of the (10 x 10) mixing matrix were chosen uniformly between the values of –1 and 1. Figure 3 shows the resulting plot of SDR as a function of number of iterations for the three algorithms.
It should be kept in mind that the amount of computation for Renyi’s MMI method is O(L2), whereas it is O(N) for the other two methods. However, the number of samples required for convergence is approximately two orders of magnitude less for Renyi’s MMI method. A larger step size was tried for the ME method in order to help speed convergence, but the resulting SDR after convergence was around 8 dB. Also notice that the final solution for Renyi’s is slightly better than the H-J network and is 10 dB superior to the ME method. One interesting note is that Renyi’s MMI method extracted all of its information from only 200 random samples. Although the iterations occurred over this fixed data set, if it were desired, new samples could be selected and used for each iteration. On the other hand, the other methods used the entire length of the data, 35000 samples, multiple times. This is particularly important when trying to track a changing environment. 
6	Conclusions
A blind source separation algorithm has been developed that is based on the minimization of Renyi’s mutual information. While this algorithm is more complex 

Figure 3.  SDR plots for N = 10. Notice that the x-axis is logarithmic.

than its predecessors, it has been shown to outperform the H-J network and the Bell and Sejnowski ME method in terms of final SDR and in its ability to extract information from a few data samples. Finally, it has been shown that the speed of convergence in terms of the number of iterations is also significantly improved.
Acknowledgments
This work was partially supported by NSF ECS #9900394.
References
[1] J. F. Cardoso, “Blind signal separation,” Proc. of the IEEE, Vol. 86, No. 10, Oct. 1998.
[2] A. Hyvarinen, “Survey on independent component analysis,” Neural Computing Surveys 2, pp. 94-128, 1999.
[3] H. H. Yang and S. Amari, "Adaptive online learning algorithms for blind separation: maximum entropy and minimum mutual information," Neural Computation, Vol. 9, No. 7, pp. 1457-1482, Oct. 1997.
[4]  P. Comon, "Independent component analysis, a new concept?," Signal Proc., Vol. 36, pp. 287-314, 1994.
[5] D. Xu, J.C. Principe, J. Fisher III, and H-C Wu, “A novel measure for independent component analysis (ICA),” ICASSP ’98.
[6] A. Rényi, Probability Theory, North-Holland Publishing Company, Amsterdam, 1970.
[7] C. Jutten and J. Herault, “Blind separation of sources, part I: An adaptive algorithm based on neuromimetic architecture,” Signal Proc. Vol. 24, 1991.
[8] A. J. Bell and T. J. Sejnowski, “An Information-Maximization Approach to Blind Separation and Blind Deconvolution,” Neural Computation, Vol. 7, No. 6, pp. 1129-1159, Nov. 1995.