Blind Source Separation Based on Joint Diagonalization in R: The Packages JADE and BSSasymp

,


Introduction
The blind source separation (BSS) problem is, in its most simple form, the following: Assume that observations x 1 , . . ., x n are p-variate vectors whose components are linear combinations of the components of p-variate unobservable zero mean vectors z 1 , . . ., z n .If we consider p-variate vectors x and z as row vectors (to be consistent with the programming language R), the BSS model can be written as where A is an unknown full rank p × p mixing matrix and µ is a p-variate location vector.The goal is then to estimate an unmixing matrix, W = A −1 , based on the n×p data matrix X = [x ¦ 1 , . . ., x ¦ n ] ¦ , such that z i = (x i − µ)W ¦ , i = 1, . . ., n.Notice that BSS can also be applied in cases where the dimension of x is larger than that of z by applying a dimension reduction method at rst stage.In this paper we, however, restrict to the case where A is a square matrix.
The unmixing matrix W cannot be estimated without further assumptions on the model.There are three major BSS models which dier in their assumptions made upon z: In the independent component analysis (ICA), which is the most popular BSS approach, it is assumed that the components of z are mutually independent and at most one of them is Gaussian.ICA applies best to cases where also z 1 , . . ., z n are independent and identically distributed (iid).The two other main BSS models, the second order source separation (SOS) model and the second order nonstationary source separation (NSS) model, utilize temporal or spatial dependence within each component.In the SOS model, the components are assumed to be uncorrelated weakly (second-order) stationary time series with dierent time dependence structures.The NSS model diers from the SOS model in that the variances of the time series components are allowed to be nonstationary.All these three models will be dened in detail later in this paper.
None of the three models has a unique solution.This can be seen by choosing any p × p matrix C from the set C = {C : each row and column of C has exactly one non-zero element}. (2) Then C is invertible, A * = AC −1 is of full rank, the components of z * = zC ¦ are uncorrelated (and independent in ICA) and the model can be rewritten as x = z * A * ¦ .Thus, the order, signs and scales of the source components cannot be determined.This means that, for any given unmixing matrix W , also W * = CW with C ∈ C is a solution.
As the scales of the latent components are not identiable, one may simply assume that COV(z) = I p .Let then Σ = COV(x) = AA ¦ denote the covariance matrix of x, and further let Σ −1/2 be the symmetric matrix satisfying Σ −1/2 Σ −1/2 = Σ −1 .Then, for the standardized random variable x st = (x − µ)Σ −1/2 , we have that z = x st U ¦ for some orthogonal U [Miettinen et al., 2015, Theorem 1].Thus, the search for the unmixing matrix W can be separated into nding the whitening (standardization) matrix Σ −1/2 and the rotation matrix U .The unmixing matrix is then given by W = U Σ −1/2 .In this paper, we describe the R package JADE which oers several BSS methods covering all three major BSS models.In all of these methods, the whitening step is performed using the regular covariance matrix whereas the rotation matrix U is found via joint diagonalization.The concepts of simultaneous and approximate joint diagonalization are recalled in Section 2, and several ICA, SOS and NSS methods based on diagonalization are described in Sections 3, 4 and 5, respectively.As performance indices are widely used to compare dierent BSS algorithms, we dene some popular indices in Section 6.We also introduce the R package BSSasymp which includes functions for computing the asymptotic covariance matrices and their data-based estimates for most of the BSS estimators in the package JADE.Section 7 describes the R packages JADE and BSSasymp, and in Section 8 we illustrate the use of these packages via simulated and real data examples.
2 Simultaneous and approximate joint diagonalization 2. 1 Simultaneous diagonalization of two symmetric matrices Let S 1 and S 2 be two symmetric p × p matrices.If S 1 positive denite, then there is a nonsingular p × p matrix W and a diagonal p × p matrix D such that If the diagonal values of D are distinct, the matrix W is unique up to a permutation and sign changes of the rows.Notice that the requirement that either S 1 or S 2 is positive denite is not necessary; there are more general results on simultaneous diagonalization of two symmetric matrices, see for example Golub and Van Loan [2002].However, for our purposes the assumption on positive deniteness is not too strong.
The simultaneous diagonalizer can be solved as follows.First solve the eigenvalue/eigenvector problem and dene the inverse of the square root of S 1 as Next solve the eigenvalue/eigenvector problem The simultaneous diagonalizer is then

Approximate joint diagonalization
Exact diagonalization of a set of symmetric p × p matrices S 1 , . . ., S K , K > 2 is only possible if all matrices commute.As shown later in Sections 3, 4 and 5, in BSS this is, however, not the case for nite data and we need to perform approximate joint diagonalization, that is, we try to make W S K W ¦ as diagonal as possible.In practice, we have to choose a measure of diagonality M , a function that maps a set of p × p matrices to [0, ∞), and seek W that minimizes Usually the measure of diagonality is chosen to be where o(V ) has the same o-diagonal elements as V , and the diagonal elements are zero.In common principal component analysis for positive denite matrices, Flury [1984] used the measure Obviously the sum of squares criterion is minimized by the trivial solution W = 0.The most popular method to avoid this solution is to diagonalize one of the matrices, then transform the rest K −1 matrices, and approximately diagonalize them requiring the diagonalizer to be orthogonal.To be more specic, suppose that S 1 is a positive denite p × p matrix.Then nd ) ¦ , k = 2, . . ., K. Notice that in classical BSS methods, matrix S 1 is usually the covariance matrix, and the transformation is called whitening.Now if we measure the diagonality using the sum of squares of the o-diagonal elements, the approximate joint diagonalization problem is equivalent to nding an orthogonal p × p matrix U that minimizes Since the sum of squares remains the same when multiplied by an orthogonal matrix, we may equivalently maximize the sum of squares of the diagonal elements Several algorithms for orthogonal approximate joint diagonalization have been suggested, and in the following we describe two algorithms which are given in the R package JADE.For examples of nonorthogonal approaches, see R package jointDiag and references therein as well as Yeredor [2002].
The rjd algorithm uses Given's (or Jacobi) rotations to transform the set of matrices to a more diagonal form two rows and two columns at a time [Clarkson, 1988].Givens rotation matrix is given by In rjd algorithm the initial value for the orthogonal matrix U is I p .First, the value of θ is computed using the elements (S * k ) 11 , (S * k ) 12 and (S * k ) 22 , k = 2, . . ., K, and matrices U, S * 2 , . . ., S * K are then updated by Similarly all pairs i < j are gone through.When θ = 0, the Givens rotation matrix is identity and no more rotation is done.Hence, the convergence has been reached when θ is small for all pairs i < j.Based on vast simulation studies it seems that the solution of the rjd algorithm always maximizes the diagonality criterion (3).
In the deation based joint diagonalization (djd) algorithm the rows of the joint diagonalizer are found one by one [Nordhausen et al., 2012].Following the notations above, assume that S * 2 , . . ., S * K , K g 2, are the symmetric p × p matrices to be jointly diagonalized by an orthogonal matrix, and write the criterion (3) as where u j is the jth row of U .The sum (4) can then be approximately maximized by solving successively for each j = 1, . . ., p − 1, u j that maximizes under the constraint u r u ¦ j = δ rj , r = 1, . . ., j − 1. Recall that δ rj = 1 as r = j and zero otherwise.
The djd algorithm in the R package JADE is based on gradients, and to avoid stopping to local maxima, the initial value for each row is chosen from a set of random vectors so that criterion ( 5) is maximized in that set.The djd function also has an option to choose the initial values to be the eigenvectors of the rst matrix S * 2 which makes the function faster, but does not guarantee that the local maximum is reached.Recall that even if the algorithm nds the global maximum in every step, the solution only approximately maximizes the criterion (4).
In the djd function also criteria of the form can be used instead of (5), and if all matrices are positive denite, also The joint diagonalization plays an important role is BSS.In the next sections, we recall the three major BSS models, and corresponding separation methods which are based on the joint diagonalization.
All these mehods are included in the R package JADE.

Independent Component Analysis
The independent component model assumes that the source vector z in model ( 1) has mutually independent components.Based on this assumption, the mixing matrix A in ( 1) is not well-dened, therefore some extra assumptions are usually made.Common assumptions on the source variable z in the IC model are (IC1) the source components are mutually independent, (IC2) E(z) = 0 and E(z ¦ z) = I p , (IC3) at most one of the components is gaussian, and (IC4) each source component is independent and identically distributed, Assumption (IC2) xes the variances of the components, and thus the scales of the rows of A. Assumption (IC3) is needed as, for multiple normal components, the independence and uncorrelatedness are equivalent.Thus, any orthogonal transformation of normal components preserves the independence.
Classical ICA methods are often based on maximizing the non-Gaussianity of the components.This approach is motivated by the central limit theorem which, roughly speaking, says that the sum of random variables is more Gaussian than the summands.Several dierent methods to perform ICA are proposed in the literature.For general overviews, see for example Hyvärinen et al. [2001], Comon and Jutten [2010], Oja and Nordhausen [2012], Yu et al. [2014].
In the following, we review two classical ICA methods, FOBI and JADE, which utilize joint diagonalization when estimating the unmixing matrix.As the FOBI method is a special case of ICA methods based on two scatter matrices with so-called independence property [Oja et al., 2006], we will rst recall some related denitions.

3.1
Scatter Matrix and Independence Property Let F x denote the cdf of a p-variate random vector x.A matrix valued functional S(F x ) is called a scatter matrix if it is positive denite, symmetric and ane equivariant in the sense that S(F Ax+b ) = AS(F x )A ¦ for all x, full rank matrices p × p matrices A and all p-variate vectors b.Oja et al. [2006] noticed that the simultaneous diagonalization of any two scatter matrices with the independence property yields the ICA solution.The issue was further studied in Nordhausen et al. [2008a].A scatter matrix S(F x ) with the independence property is dened to be a diagonal matrix for all x with independent components.An example of a scatter matrix with the independence property is the covariance matrix, but what comes to most scatter matrices, they do not possess the independence property (for more details, see Nordhausen and Tyler [2015]).However, it was noticed in Oja et al. [2006] that if the components of x are independent and symmetric, then S(F x ) is diagonal for any scatter matrix.Thus a symmetrized version of a scatter matrix S sym (F x ) = S(F x1−x2 ), where x 1 and x 2 are independent copies of x, always has the independence property, and can be used to solve the ICA problem.
The ane equivariance of the matrices, which are used in the simultaneous diagonalization and approximate joint diagonalization methods, imply the ane equivariance of the unmixing matrix estimator.More precisely, if the unmixing matrices W and W * correspond to x and x * = xB ¦ , respectively, then xW ¦ = x * W * ¦ (up to sign changes of the components) for all p × p full rank matrices B. This is a desirable property of an unmixing matrix estimator as it means that the separation result does not depend on the mixing procedure.It is easy to see that the ane equivariance also holds even if S 2 , . . ., S K , K g 2, are only orthogonal equivariant.

FOBI
One of the rst ICA methods, FOBI (fourth order blind identication) introduced by Cardoso [1989], uses simultaneous diagonalization of the covariance matrix and the matrix based on the fourth moments, respectively.Notice that both S 1 and S 2 are scatter matrices with the independence property.The unmixing matrix is the simultaneous diagonalizer W satisfying where the diagonal elements of D are the eigenvalues of S 2 (F z ) given by E[z 4 i ] + p − 1, i = 1, . . ., p. Thus, for a unique solution, FOBI requires that the independent components have dierent kurtosis values.The statistical properties of FOBI are studied in Ilmonen et al. [2010a] and Miettinen et al. [2015].

JADE
The JADE (joint approximate diagonalization of eigenmatrices) algorithm [Cardoso and Souloumiac, 1993] can be seen as a generalization of FOBI since both of them utilize fourth moments.For a recent comparison of these two methods, see Miettinen et al. [2015].Contrary to FOBI, the kurtosis values do not have to be distinct in JADE.The improvement is gained by increasing the number of matrices to be diagonalized as follows.Dene, for any p × p matrix M , the fourth order cumulant matrix as where x st is a standardized variable.Notice that C(I p ) is the matrix based on the fourth moments used in FOBI.Write then E ij = e ¦ i e j , i, j = 1, . . ., p, where e i is a p-vector with the ith element one and others zero.In JADE (after the whitening) the matrices C(E ij ), i, j = 1, . . ., p are approximately jointly diagonalized by an orthogonal matrix.The rotation matrix U thus maximizes the approximate joint diagonalization criterion JADE is ane equivariant even though the matrices C(E ij ), i, j = 1, . . ., p, are not orthogonal equivariant.If the eighth moments of the independent components are nite, then the vectorized JADE unmixing matrix estimate has a limiting multivariate normal distribution.For the asymptotic covariance matrix and a detailed discussion about JADE, see Miettinen et al. [2015].
The JADE estimate jointly diagonalizes p 2 matrices.Hence its computational load grows quickly with the number of components.Miettinen et al. [2013] suggested a quite similar, but faster method, called k-JADE which is computationally much simpler.The k-JADE method whitens the data using FOBI and then jointly diagonalizes The value k f p can be chosen by the user and corresponds basically to the guess of the largest multiplicity of identical kurtosis values of the sources.If k is larger or equal to the largest multiplicity, then k-JADE and JADE seem to be asymptotically equivalent.

Second Order Source Separation
In second order source separation (SOS) model, the source vectors compose a p-variate time series z = (z t ) t=0,±1,±2,... that satises (SOS1) E(z t ) = 0 and E(z ¦ t z t ) = I p , and Above assumptions imply that the components of z are weakly stationary and uncorrelated time series.
In the following we will shortly describe two classical (SOS) methods, which yield ane equivariant unmixing matrix estimates.
The AMUSE (Algorithm for Multiple Unknown Signals Extraction) [Tong et al., 1990] algorithm uses the method of simultaneous diagonalization of two matrices.In AMUSE, the matrices to be diagonalized are the covariance matrix and the autocovariance matrix with chosen lag τ , that is, The unmixing matrix W τ then satises The requirement for distinct eigenvalues implies that the autocorrelations with the chosen lag need to be unequal for the source components.Notice that, as the population quantity S τ (F x ) is symmetric, the estimate Ŵτ is obtained by diagonalizing the sample covariance matrix and the symmetrized autocovariance matrix with lag τ .The sample autocovariance matrix with lag τ is given by and the symmetrized autocovariance matrix with lag τ , respectively.
It has been shown that the choice of the lag has a great impact on the performance of the AMUSE estimate [Miettinen et al., 2012].However, without any preliminary knowledge of the uncorrelated components it is dicult to choose the best lag for the problem at hand.Cichocki and Amari [2002] simply recommend to start with τ = 1, and check the diagonal elements of the estimate D. If there are two almost equal values, another value for τ should be chosen.Belouchrani et al. [1997] provide a natural approximate joint diagonalization method for SOS model.In SOBI (Second Order Blind Identication) the data is whitened using the covariance matrix S 0 (F x ) = COV(x).The K matrices for rotation are then autocovariance matrices with distinct lags τ 1 , . . ., τ K , that is, S τ1 (F x ), . . ., S τ K (F x ).The use of dierent joint diagonalization methods yields estimates with dierent properties.For details about deation-based algorithm (djd) in SOBI see Miettinen et al. [2014a], and for details about SOBI using the rjd see Miettinen et al. [2016].General agreement seems to be that in most cases, the use of rjd in SOBI is preferrable.
The problem of choosing the set of lags τ 1 , . . ., τ K for SOBI is not as important as the choice of lag τ for AMUSE.Among the signal processing community, K = 12 and τ k = k for k = 1, . . ., K, are conventional choices.Miettinen et al. [2014a] argue that, when the deation-based joint diagonalization is applied, one should rather take too many matrices than too few.The same suggestion applies to SOBI using the rjd.If the time series are linear processes, the asymptotic results in Miettinen et al. [2016] provide tools to choose the set of lags, see also Example 2 in Section 8.

Nonstationary Source Separation
The SOS model assumptions are sometimes considered to be too strong.The NSS model is a more general framework for cases where the data are ordered observations.In addition to the basic BSS model ( 1) assumptions, the following assumptions on the source components are made: ) is positive denite and diagonal for all t, (NSS3) E(z ¦ t z t+τ ) is diagonal for all t and τ .
Hence the source components are uncorrelated and they have a constant mean.However, the variances are allowed to change over time.Notice that this denition diers from the block-nonstationary model, where the time series can be divided into intervals so that the SOS model holds for each interval.NSS-SD, NSS-JD and NSS-TD-JD are algorithms that take into account the nonstationarity of the variances.For the description of the algorithms dene where T is a nite time interval and τ ∈ {0, 1, . ..}.
The NSS-SD unmixing matrix simultaneously diagonalizes S T1,0 (F x ) and S T2,0 (F x ), where T 1 , T 2 ¢ [1, n] are separate time intervals.T 1 and T 2 should be chosen so that S T1,0 (F x ) and S T2,0 (F x ) are as dierent as possible.
Corresponding approximate joint diagonalization method is called NSS-JD.The data is whitened using the covariance matrix S [1,n],0 (F x ) computed from all the observations.After whitening, the K covariance matrices S T1,0 (F x ), . . ., S T K ,0 (F x ) related to time intervals T 1 , . . ., T K are diagonalized with an orthogonal matrix.
Both NSS-SD and NSS-JD algorithms ignore the possible time dependence.Assume that the full time series can be divided into K time intervals T 1 , . . ., T K so that, in each interval, the SOS model holds approximately.Then the autocovariance matrices within the intervals make sense, and the NSS-TD-JD algorithm is applicable.Again, the covariance matrix S [1,n],0 (F x ) whitens the data.Now the matrices to be jointly diagonalized are S Ti,τj (F x ), i = 1, . . ., K, j = 1, . . ., L. When selecting the intervals one should take into account the lengths of the intervals so that the random eect is not too large when the covariances and the autocovariances are computed.A basic rule in signal processing community is to have 12 intervals if the data is large enough, or K < 12 intervals such that each interval contains at least 100 observations.Notice that NSS-SD and NSS-JD (as well as AMUSE and SOBI) are special cases of NSS-TD-JD.Naturally, NSS-SD, NSS-JD and NSS-TD-JD are all ane equivariant.
For further details on NSS methods see for example Choi and Cichocki [2000a,b], Choi et al. [2001], Nordhausen [2014].Notice that asymptotic results are not yet available for any of these NSS methods.

BSS performance criteria
The performance of dierent BSS methods using real data is often dicult to evaluate since the true unmixing matrix is unknown.In simulations studies, however, the situation is dierent, and in the literature many performance indices have been suggested to measure the performance of dierent methods.For a recent overview see Nordhausen et al. [2011], for example.
The package JADE contains several performance indices but in the following we will only introduce two of them.Both performance indices are functions of the so-called gain matrix, Ĝ, which is a product of the unmixing matrix estimate Ŵ and the true mixing matrix, that is, Ĝ = Ŵ A. Since the order, the signs and the scales of the source components cannot be estimated, the gain matrix of an optimal estimate does not have to be identity, but equivalent to the identity in the sense that Ĝ = C for some C ∈ C, where C is given in (2).
The Amari error [Amari et al., 1996] is dened as where ĝij denotes the ijth element of Ĝ.The range of the Amari error values is [0, 1], and a small value corresponds to a good separation performance.The Amari error is not scale invariant.Therefore, when dierent algorithms are compered, the unmixing matrices should be scaled in such a way that the corresponding rows of dierent matrices are of equal length.
The minimum distance index [Ilmonen et al., 2010b] is dened as where ∥ • ∥ is the matrix (Frobenius) norm and C is dened in (2).Also the MD index is scaled to have a maximum value 1, and M D( Ĝ) = 0 if and only if Ĝ ∈ C. The MD index is ane invariant.The statistical properties of the index are thoroughly studied in Ilmonen et al. [2010b] and Ilmonen et al. [2012].
A feature that makes the minimum distance index especially attractive in simulation studies is that its value can be related to the asymptotic covariance matrix of an estimator Ŵ .If Ŵ → A −1 and √ n vec( Ŵ A − I p ) → N p 2 (0, Σ), which is for example the case for FOBI, JADE, AMUSE and SOBI, then the limiting distribution of n(p − 1)M D( Ĝ) 2 has an expected value where D p,p = p i=1 e ¦ i e i ¹ e ¦ i e i , with ¹ denoting the kronecker product and e i a p-vector with ith element one and others zero.
Notice that tr (I p 2 − D p,p )Σ(I p 2 − D p,p ) is the sum of the o-diagonal elements of Σ and therefore a natural measure of the variation of the unmixing matrix estimate Ŵ .We will make use of this relationship later in one of our examples.

Functionality of the packages
The package JADE is freely available from the Comprehensive R Archive Network at http://CRAN.R-project.org/package=JADEand comes under the GNU General Public Licence (GPL) 2.0 or higher licence.
The main functions of the package implement the blind source separation methods described in the previous sections.The function names are selfexplanatory being FOBI, JADE and k_JADE for ICA, AMUSE and SOBI for SOS and NSS.SD, NSS.JD and NSS.TD.JD for NSS.All functions usually take as an input either a numerical matrix X or as alternative a multivariate time series object of class ts.The functions have method appropriate arguments like for example which lags to choose for AMUSE and SOBI.
All functions return an object of the S3-class bss which contains at least an object W, which is the estimated unmixing matrix, and S containing the estimated (and centered) sources.Depending on the chosen function also other information is stored.The methods available for the class are print: which prints basically all information except the sources S. coef: which extracts the unmixing matrix W. plot: which produces a scatter plot matrix of S using pairs for ICA methods and a multivariate time series plot using plot.tsfor other BSS methods.
To extract the sources S the helper function bss.components can be used.
The functions which use joint approximate diagonalization of several matrices provide the user an option to choose the method for joint diagonalization from the list below.
rjd: for joint diagonalization using Givens rotations.frjd: which is basically the same as rjd, but has less options and is much faster as being implemented in C.
From our experience the function frjd, when appropriate, seems to obtain the best results.
In addition, the JADE package provides two other functions for joint diagonalization.The function FG is designed for diagonalization of real positive-denite matrices and cjd is the generalization of rjd to the case of complex matrices.For details about all functions for joint diagonalization see also their help pages.More functions for joint diagonalization are also available in the R package jointDiag [Gouy-Pailler, 2009].
To evaluate the performance of BSS methods using simulation studies, performance indices are needed.The package provides for this purpose the functions amari_error, ComonGAP, MD and SIR.Our personal favorite is the MD-function which implements the minimum distance index described in Section 6.
For further details on all the functions see their help pages and the references therein.
For ICA, many alternative methods are implemented in other R packages.Examples include fas-tICA [Marchini et al., 2013], fICA [Miettinen et al., 2014b], mlica2 [Teschendor, 2012], PearsonICA [Karvanen, 2009] and ProDenICA [Hastie and Tibshirani, 2010].None of these ICA methods uses joint diagonalization in estimation.Two slightly overlapping packages with JADE are ICS [Nordhausen et al., 2008b] which provides a generalization of the FOBI method, and ica [Helwig, 2014] which includes the JADE algorithm.In current practise JADE and fastICA (implemented for example in the packages fas-tICA and fICA) seem to be the most often used ICA methods.Other newer ICA methods, as for example ICA via product density estimation as provided in the package ProDenICA, are often computationally very intensive as the sample size is usually high in typical ICA applications.
To the best of our knowledge there are currently no other R packages for SOS or NSS available.Many methods included in the JADE package are also available in the MATLAB toolbox ICALAB [Cichocki et al., 2014] which accompanies the book of Cichocki and Amari [2002].A collection of links to JADE implementations for real and complex values in dierent languages like MATLAB, C and python as well as some joint diagonalization functions for MATLAB are available on J.-F.Cardoso's homepage http://perso.telecom-paristech.fr/~cardoso/guidesepsou.html.
The package BSSasymp is freely available from the Comprehensive R Archive Network at http: //CRAN.R-project.org/package=BSSasymp and comes under the GNU General Public Licence (GPL) 2.0 or higher licence.
There are two kinds of functions in the package.The rst set of functions compute the asymptotic covariance matrices of the vectorized mixing and unmixing matrix estimates under dierent BSS models.The others estimate the covariance matrices based on a data matrix.The package BSSasymp includes functions for several estimators implemented in package JADE.They are FOBI and JADE in the IC model and AMUSE, deation-based SOBI and regular SOBI in the SOS model.The asymtotic covariance matrices for FOBI and JADE estimates are computed using the results in Miettinen et al. [2015].For the limiting distributions of AMUSE and SOBI estimates, see Miettinen et al. [2012] and Miettinen et al. [2016], respectively.
Functions ASCOV_FOBI and ASCOV_JADE compute the theoretical values for covariance matrices.The argument sdf is the vector of source density functions standardized to have mean zero and variance equal to one, supp is a two column matrix, whose rows give the lower and the upper limits used in numerical integration for each source component and A is the mixing matrix.The corresponding functions ASCOV_SOBIdefl and ASCOV_SOBI in the SOS model take as input the matrix psi, which gives the MAcoecients of the source time series, the vector of integers taus for the lags, a matrix of fourth moments of the innovations Beta (default value is for gaussian innovations) and the mixing matrix A.
Functions ASCOV_FOBI_est, ASCOV_JADE_est, ASCOV_SOBIdefl_est and ASCOV_SOBI_est can be used for approximating the covariance matrices of the estimates.They are based on asymptotical results, and therefore the sample size should not be very small.Argument X can be either the observed data or estimated source components.When argument mixed is set TRUE, X is considered as observed data and the unmixing matrix is rst estimated using the method of interest.The estimation of the covariance matrix of the SOBI estimate is also based on the assumption that the time series are stationary linear processes.If the time series are gaussian, then the asymptotic variances and covariances depend only on the autocovariances of the source components.Argument M gives the number of autocovariances to be used in the approximation of the innite sums of autocovariances.Thus, M should be the largest lag for which any of the source time series has non-zero autocovariance.In the non-gaussian case the coecients of the linear processes need to be computed.In functions ASCOV_SOBIdefl_est and ASCOV_SOBI_est, ARMA-parameter estimation is used and arguments arp and maq x the order of ARMA series, respectively.There are also faster functions ASCOV_SOBIdefl_estN and ASCOV_SOBI_estN, which assume that the time series are gaussian and do not estimate the MA-coecients.The argument taus is to dene the lags of the SOBI estimate.
All functions for the theoretical asymptotic covariance matrices return lists with ve components.A and W are the mixing and unmixing matrices and COV_A and COV_W are the corresponding asymptotic covariance matrices.In simulations studies where the MD index is used as performance criterion, the sum of the variance of the o-diagonal values is of interest (recall Section 6 for details).This sum is returned as object EMD in the list.
The functions for the estimated asymptotic covariance matrices return similar lists as their theoretic counterparts excluding the component EMD.

Examples
In this section we provide four examples to demonstrate how to use the main functions in the packages JADE and BSSasymp.In Section 8.1 we show how dierent BSS methods can be compared using an articial data.Section 8.2 demonstrates how the package BSSasymp can help in selecting the best method for the source separation.In Section 8.3 we show how a typical simulation study for the comparison of BSS estimators can be performed using the packages JADE and BSSasymp, and nally, in Section 8.4 a real data example is considered.In these examples the dimension of the data is relatively small, but for example in Joyce et al. [2004] SOBI has been succesfully applied to analyze EEG data where the electrical activity of the brain is measured by 128 sensors on the scalp.As mentioned earlier, computation of the JADE estimate for such high-dimensional data is demanding because of the large number of matrices and the use of k-JADE is recommended then.
In the examples we use the option options(digits = 4) in R 3.2.1 [R Core Team, 2015] together with the packages JADE 1.9-93, BSSasymp 1.0-2 and tuneR 1.2.1 [Ligges et al., 2014] for the output.Random seeds (when applicable) are provided for reproducibility of examples.

Example 1
A classical example of the application of BSS is the so-called cocktail party problem.To separate the voices of p speakers, we need p microphones in dierent parts of the room.The microphones record the mixtures of all p speakers and the goal is then to recover the individual speeches from the mixtures.To illustrate the problem the JADE package contains in its subfolder datafiles three audio les which are often used in BSS 1 .For demonstration purpose we mix the audio les and try to recover the original sounds again.The cocktail party problem data can be created using the packages To use the les in R we use functions from the package tuneR and then delete again the downloaded les.

R> play(x1) R> play(x2) R> play(x3) R> play(x4)
To demonstrate the use of BSS methods, assume now that we have observed the mixture of unknown source signals plotted in Figure 2. The aim is then to estimate the original sound signals based on this observed data.The question is then, which method to use.Based on Figure 2, the data are neither iid nor second order stationary.Nevertheless, we rst apply JADE, SOBI and NSSTDJD with their default settings: R> jade <-JADE(X) R> sobi <-SOBI(Xt) R> nsstdjd <-NSS.TD.JD(Xt)  All three objects are then of class bss and for demonstration purposes we look at the output of the call to SOBI.

R> MD(coef(jade
MD indices show that NSSTDJD performs best and that JADE is the worst method here.This result is in agreement with how well the data meets the assumptions of each method.The SOBI with the second set of lags is better than the default SOBI.In Section 8.2 we show how the package BSSasymp can be used to select a good set of lags.
To play the sounds recovered by NSSTDJD, one can use the function bss.components to extract the estimated sources and convert them back to audio.

8.3
Example 3 In simulation studies usually several estimators are compared and it is of interest to study which of the estimators performs best under the given model and also how fast the estimators converge to their limiting distributions.In the following we will perform a simulation study similar to that of Miettinen et al. [2016] and compare the performances of FOBI, JADE and 1-JADE using the package BSSasymp.
Consider the ICA model where the three source component distributions are exponential, uniform and normal distributions, all of them centered and scaled to have unit variances.Due to the ane equivariance of the estimators, the choice of the mixing matrix does not aect the performances, and we can choose A = I 3 for simplicity.
We rst create a function ICAsim which generates the data and then computes the MD indices using the unmixing matrices estimated with the three ICA methods.The arguments in ICAsim are a vector of dierent sample sizes (ns) and the number of repetitions (repet).The function then returns a data frame with the variables N, fobi, jade and kjade, which includes the used sample size and the obtained MD index value for each run and for the three dierent methods.For each of the sample sizes, 250,500,1000,2000,4000,8000,16000 and 32000, we then generate 2000 repetitions.Notice that this simulation will take a while.R> set.seed(123)R> N <-2^(( -2):5) * 1000 R> MDs <-ICAsim(ns = N, repet = 2000) Besides the nite sample performances of dierent methods, we are interested in seeing how quickly the estimators converge to their limiting distributions.The relationship between the minimum distance index and the asymptotic covariance matrix of the unmixing matrix estimate was described in Section 6.To compute (6) we rst compute the asymptotic covariance matrices of the unmixing matrix estimates Ŵ .Since all three independent components in the model have nite eighth moments, all three estimates have a limiting multivariate normal distribution [Ilmonen et al., 2010a, Miettinen et al., 2015].The functions ASCOV_FOBI and ASCOV_JADE compute the asymptotic covariance matrices of the corresponding unmixing matrix estimates Ŵ and the mixing matrix estimates Ŵ −1 .As arguments, one needs the source density functions standardized so that the expected value is zero and the variance equals to one, and the support of each density function.The default value for the mixing matrix is the identity matrix.

R> f1 <-function
Let us next look at the simulation results concerning the FOBI method in more detail.First notice that the rows of the FOBI unmixing matrices are ordered according to the kurtosis values of resulting independent components.Since the source distributions f1, f2 and f3 are not ordered accordingly, the unmixing matrix fobi$W is dierent from the identity matrix.
To make use of the relationship between the minimum distance index and the asymptotic covariance matrices, we need to extract the asymptotic variances of the o-diagonal elements of such Ŵ A that converges to I 3 .In fact, these variances are the second, third, fourth, fth, seventh and ninth diagonal element of fobi$COV_W, but there is also an object fobi$EMD, which directly gives the sum of the variances as given in (6).
R> fobi$EMD [1] 40.45 The corresponding value for JADE can be obtained as follows.
R> jade$EMD [1] 23.03 Based on these results we can conclude that for this ICA model, the theoretically best separation method is JADE.The value n(p − 1)M D( Ĝ) 2 for JADE should converge to 23.03 and that for FOBI to 40.45.Since all three components have the dierent kurtosis values, 1-JADE is expected to have the same limiting behavior as JADE.
To compare the theoretical values to their nite sample counterparts, we next compute the average values of n(p − 1)M D( Ĝ) 2 for each sample size and each estimator, and plot them together with their limiting expected values in Figure 3.

8.4
Example 4 So far we have considered examples where the true sources and the mixing matrix have been known.
In our last example we use a real data set which includes electrocardiography (ECG) recordings of a pregnant woman.ECG measures the electrical potential, generated by the heart muscle, from the body surface.The electrical activity produced by the heart beats of a fetus can then be detected by measuring the potential on the mother's skin.As the measured signals are mixtures of the fetus's and the mother's heart beats, the goal is to use the BSS method to separate these two heart beats as well as some possible artifacts from each other.In this context it is useful to know that a fetus's heart is supposed to beat faster than that of the mother.For a more detail discussion on the data and of the use of BSS in this context, see De Lathauwer et al. [1995].
In this ECG recording, eight sensors have been placed on the skin of the mother, the rst ve in the stomach area and the other three in the chest area.The data was obtained as foetal_ecg.datfrom http://homes.esat.kuleuven.be/~smc/daisy/daisydata.html 2 and is also provided in the supplementary les of the JSS paper Miettinen et al. [2017].In this ECG recording, eight sensors have been places on the skin of the mother, the rst ve in the stomach area and the other three in the chest area.We rst load the data assuming it is in the working directory and plot it in Figure 4. R> library("JADE") R> library("BSSasymp") R> dataset <-matrix(scan(paste0("foetal_ecg.dat")),2500, 9, byrow = TRUE) Read 22500 items R> X <-dataset[ , 2:9] R> plot.ts(X,nc = 1, main = "Data") Figure 4 shows that the mother's heartbeat is clearly the main element in all of the signals.The heart beat of the fetus is visible in some signals -most clearly in the rst one.
We next scale the components to have unit variances to make the elements of the unmixing matrix larger.Then the JADE estimate is computed and resulting components are plotted in Figure 5.
R> scale(X, center = FALSE, scale = apply(X, 2, sd)) R> jade <-JADE(X) R> plot.ts(bss.components(jade),nc = 1, main = "JADE solution") From Figure 5 it is seen that the rst three components are related to the mother's heartbeat and the fourth component is related to the fetus's heartbeat.Since we are interested in the fourth component, we pick up the corresponding coecients from the fourth row of the unmixing matrix estimate.For demonstration purposes, we also derive their standard errors in order to see how much uncertainty is included in the results.These would be useful for example when selecting the best BSS method in a case where estimation accuracy of only one component is of interest, as opposed to Example 2 where the whole unmixing matrix was considered.Furthermore, we can test, for example, whether the recordings from the mother's chest area contribute to the estimate of the fourth component (fetus's heartbeat), i.e., whether the last three elements of the fourth row of the unmixing are non-zero.Since the JADE estimate is asymptotically multivariate normal, we can compute the Wald test statistic related to the null hypothesis H 0 : ((W ) 46 , (W ) 47 , (W ) 48 ) = (0, 0, 0).Notice that ascov$COV_W is the covariance matrix estimate of the vector built from the columns of the unmixing matrix estimate.Therefore we create the vector w and hypothesis matrix L accordingly.The sixth, seventh and eighth element of the fourth row of the 8 × 8 matrix are the 5 • 8 + 4 = 44th, 6 • 8 + 4 = 52nd and 7 • 8 + 4 = 60th elements of w, respectively.

Conclusions
In this paper we have introduced the R packages JADE and BSSasymp which contain several practical tools for blind source separation.
Package JADE provides methods for three common BSS models.The functions allow the user to perform blind source separation in cases where the source signals are (i) independent and identically distributed, (ii) weakly stationary time series, or (iii) time series with nonstationary variance.All BSS methods included in the package utilize either simultaneous diagonalization of two matrices or approximate joint diagonalization of several matrices.In order to make the package self-contained we have included in it several algorithms for joint diagonalization.Two of the algorithms, deation-based joint diagonalization and joint diagonalization using Givens rotations, are described in detail in this paper.
Package BSSasymp provides tools to compute the asymptotic covariance matrices as well as their data-based estimates for most of the BSS estimators included in the package JADE.The functions allow the user to study the uncertainty in the estimation either in simulation studies or in practical applications.
Notice that package BSSasymp is the rst R package so far to provide such variance estimation methods for practitioners.
We have provided four examples to introduce the functionality of the packages.The examples show in detail (i) how to compare dierent BSS methods using articial example (cocktail-party problem) or simulated data, (ii) how to select a best method for the problem at hand, and (iii) how to perform blind source separation with real data (ECG recording).

Figure 1 :
Figure 1: Original sound and noise signals.

2Figure 3 :
Figure 3: Simulation results based on 2000 repetitions.The dots give the average values of n(p − 1)M D( Ĝ) 2 for each sample size, and the horizontal lines are the expected values of the limiting of the limiting distributions of n(p − 1)M D( Ĝ) 2 for the FOBI method and the two JADE methods.

Figure 4 :
Figure 4: Electrocardiography recordings of a pregnant woman.