Generalized generating function with tucker decomposition and alternating least squares for underdetermined blind identification

Generating function (GF) has been used in blind identification for real-valued signals. In this paper, the definition of GF is first generalized for complex-valued random variables in order to exploit the statistical information carried on complex signals in a more effective way. Then an algebraic structure is proposed to identify the mixing matrix from underdetermined mixtures using the generalized generating function (GGF). Two methods, namely GGF-ALS and GGF-TALS, are developed for this purpose. In the GGF-ALS method, the mixing matrix is estimated by the decomposition of the tensor constructed from the Hessian matrices of the GGF of the observations, using an alternating least squares (ALS) algorithm. The GGF-TALS method is an improved version of the GGF-ALS algorithm based on Tucker decomposition. More specifically, the original tensor, as formed in GGF-ALS, is first converted to a lower-rank core tensor using the Tucker decomposition, where the factors are obtained by the left singular-value decomposition of the original tensor’s mode-3 matrix. Then the mixing matrix is estimated by decomposing the core tensor with the ALS algorithm. Simulation results show that (a) the proposed GGF-ALS and GGF-TALS approaches have almost the same performance in terms of the relative errors, whereas the GGF-TALS has much lower computational complexity, and (b) the proposed GGF algorithms have superior performance to the latest GF-based baseline approaches.


Introduction
Blind identification (BI) of linear mixtures has recently attracted intensive research interest in many fields of signal processing including blind source separation (BSS).This work is devoted to BI of underdetermined mixtures with complex sources.Underdetermined mixtures are commonly encountered in many practical applications, such as in the radio communication context, where the reception of more sources than sensors becomes increasingly possible with the growth in reception bandwidth.In these applications, one often has to also deal with complex sources.One reason is that the communication signals are usually complex-valued such as the quadrature-amplitude modulation (QAM) signal, and minimum-shift keying (MSK) signal.Another reason is that frequency domain methods are often used for blind separation or identification from convolutive mixtures due to its computational efficiency [1,2], while the objective functions used in the frequency domain are usually defined on complex-valued variables.
A large number of methods for BI of underdetermined mixtures start from the assumption that the sources are sparse by nature (i.e., in its own domain such as the time domain) or could be made sparse in another domain (e.g., a transform domain).A predefined transform such as short-time Fourier transform (STFT) or a learned transform using, e.g., simultaneous codeword optimization (SimCO), is usually applied to sparsify the data [3,4] if the signal by nature is not sparse.Due to the sparsity of the sources, the scatter plot typically shows high signal values in the directions of the mixing vectors, which can be localized by using some clustering techniques [5,6].It should be noted that although some signals such as speech signals have some degree of sparsity in one domain or another, many other signals such as the majority of communication signals do not possess such a property.Hence, it is necessary to develop BI methods for the underdetermined mixtures that do not impose any sparsity constraint on the sources.
To this aim, many methods for BI of underdetermined mixtures turn to the use of various decomposition methods based on different data structures such as correlation [7,8] and higher-order cumulant [9][10][11][12][13][14] matrices.The main idea of these algorithms is to construct a tensor based on the cumulants of the observations and then to estimate the mixing matrix by the decomposition of such a tensor.This is notably the case for second-order blind identification of underdetermined mixtures (SOBIUM) [7], fourth-order blind identification of underdetermined mixtures (FOBIUM) [9], fourth-order-only blind identification (FOOBI) [10], FOOBI-2 [10], and blind identification of mixtures of sources using Redundancies in the daTa Hexacovariance matrix (BIRTH) [11,12] algorithms, which use second-order statistics tensors and fourth-and sixth-order cumulant tensors, respectively.A family of the methods named blind identification of over-complete mixtures of sources (BIOME) is proposed in [13], based on the even-order cumulants of the observations.However, all the methods proposed in [7][8][9][10][11][12][13][14] exploit only the statistical information contained in the data measured by second-order or higher-order statistics.
In order to exploit statistical information more effectively, a family of BI approaches was proposed in [15][16][17][18] by exploiting the statistical information with the characteristic function (CAF) or generating function (GF).In these works, the authors showed that the mixing matrix can be estimated up to trivial scaling and permutation indeterminacies by decomposing the tensor composed of partial derivatives of the GF.It is worth mentioning that the algorithms in [15][16][17] have been only applied to BI problems involving real-valued sources.In [18], the CAF approach was extended to the case of mixtures of complex-valued sources, which often occurs in digital communications.However, extra effort is required to obtain the correct real and imaginary combination of the mixing matrix since the real and imaginary parts of the mixing matrix are treated separately, leading to an increased computational cost due to the increased dimension of the matrix that needs to be processed.In this paper, we propose the Generalized Generating Function (GGF) to exploit the statistical information carried on the complex random variable.We show that the proposed GGF can exploit the statistical information carried on complex random variables in a more effective way than the GF presented in [18] due to the algebraic structure adopted by GGF (as detailed in Algebraic structure based on generalized generating function).Furthermore, a simple method for the mixing matrix estimation is derived based on tensor decomposition where the tensor is composed of the Hessian matrices of the GGF of the observations.
The remainder of this paper is organized as follows.In Problem formulation the BI problem is formulated and relevant assumptions are presented.In Algebraic structure based on generalized generating function we firstly generalize the definition of the GF for complex-valued random variables and then derive the corresponding core equation for BI.In Blind identification based on tensor decomposition, the GGF-ALS and GGF-TALS approaches are developed for the estimation of the mixing matrix.In the GGF-ALS algorithm, the mixing matrix is estimated by directly decomposing the tensor, constructed from the Hessian matrices of the second GGF of the observations, using the alternating least squares (ALS) algorithm.In the GGF-TALS algorithm, the Tucker decomposition is firstly applied to convert the original tensor to a lowerorder core tensor, then the mixing matrix is obtained by decomposing the core tensor with the ALS algorithm.Furthermore, the factors of Tucker decomposition are obtained by the left singular vectors of the original tensor's mode-3 matrix.Computer simulations are used to illustrate the performance of the proposed GGF approaches in Simulations and analysis.Finally, the paper is concluded in Conclusions.

Problem formulation
Considering the following linear mixture model where the stochastic vector z(t) ∈ C Q represents the observation signals, s(t) ∈ C P contains the unobserved source signals, and w(t) ∈ C Q denotes additive noise.From now on, the noise w(t) is simply ignored for convenience, except when running computer experiments.The unknown mixing matrix A ∈ C Q × P characterizes the way that the sources are acquired by the sensors.BI aims to estimate the mixing matrix from the observations based on the assumption that the source signals are statistically independent.The mixing matrix obtained may in turn be used to estimate the original source signals from the observations.In addition, we make the following assumptions: (i).The mixing matrix A is of full (row) rank.
(ii).The number P of sources is known.
(iii).The number of sensors is smaller than the number of sources, i.e., Q < P.

Core equation based on generalized generating function
For a real stochastic vector x ∈ R Q , the GF ϕ x (u) obtained by dropping the term of the square root of (−1) in the exponent of a CAF is defined as where u ∈ R Q is an arbitrary vector referred to as a processing point [15], and E[ ] denotes an expectation operator.Nevertheless, both the observation vector z and the mixing matrix A discussed in this paper belong to the complex field.Hence, a definition of GF for complex variables is required.One such definition has been presented in [18] as where R ⋅ ð Þ and I ⋅ ð Þ denote taking the real and imaginary parts from their arguments (i.e., complex-valued vectors) to form a real-valued vector of the same dimension.It is actually defined by assimilating C to R 2 .Thus the GF of a complex variable in ( 3) is defined as a function of the real and imaginary parts.In this paper, we generalize the definition of GF for real stochastic vector in (2) to the following complex form Note that the statistical information exploited by GF/ GGF is related to the number of processing points.Theoretically, a complete statistical description of the probability density function requires the evaluation of the GF/GGF at all (infinitely many) possible processing points.However, this often becomes computationally infeasible.In practice, such statistical information is obtained approximately by the evaluation of GF/GGF at a finite number of processing points.Hence, in comparison with the GF presented in [18], the GGF defined in (4) can exploit the statistical information carried on the complex variables more effectively when the number of the processing points stays the same, thanks to the incorporation of the imaginary part of the exponent to the function.Furthermore, as compared with the use of the GF in (3), using the GGF in (4) offers a simpler way for the estimation of the mixing matrix due to the exploitation of an elegant algebraic structure.Now, replacing z by its model and neglecting the noise contribution yield Defining φ z (u) = log ψ z (u), which is often referred to as the 'second' GGF, and using the source independence property, the second GGF of the observations can be rewritten as Consequently, by calculating the derivative of the conjugate gradient of φ z (u) with respect to u (more details can be found in Appendix 1), we can obtain the following core equation for the Hessian matrix ψ z (u), with where (•) * denotes the conjugate operator.It is necessary to point out that ψ s (A H u) is a diagonal matrix (more details can be found in Appendix 2).

Estimating ψ z (u)
In this subsection, we discuss how to consistently estimate the Hessian matrix ψ z (u).Under the ergodicity assumption, the mean value of a random variable can be estimated by a time average.Hence, we can estimate the GGF of the observation vector as The conjugate gradient ς z (u) = ∂ ψ z (u)/∂ u* which is a Q × 1 vector can be estimated by Based on the above analysis, the Hessian matrix ψ z (u) of the second GGF φ z (u) can be obtained as

Blind identification based on tensor decomposition
In this section, the GGF-ALS and GGF-TALS algorithms are developed for the estimation of the mixing matrix.
In the GGF-ALS algorithm, the mixing matrix is estimated by decomposing the tensor which is formed of the Hessian matrices of the GGF of the observations.An improved version of the GGF-ALS algorithm, i.e., the GGF-TALS algorithm, is also developed, where the Tucker decomposition is firstly employed to convert the original tensor as used in the GGF-ALS algorithm to a lower-rank core tensor, and the mixing matrix is then estimated by decomposing the core tensor with the ALS algorithm.

The GGF-ALS algorithm
Evaluating the Hessian matrices ψ z (u) of GGF at a series of processing points u 1 , u 2 ,..., u K and using Equation ( 7), one can obtain the following joint diagonalization (JD) problem The problem we need to address is to estimate the mixing matrix A based on the set {ψ z (u 1 ), ⋯, ψ z (u K )}.For the determined/overdetermined case, it is obvious that JD methods, such as the AC-DC method [19], can be used to estimate the mixing matrix.However, this method does not work when Q < P i.e., in the underdetermined case.As shown in [7], the JD problem ( 9) can be seen as a particular case of the parallel factor (PARAFAC) decomposition, also known as canonical decomposition (CAND), of the third-order tensor M ∈ C QÂQÂK built by stacking the K matrices ψ x (u k ) along the third mode.Specifically, the tensor M ∈ C QÂQÂK is built by stacking ψ z (u 1 ), ψ z (u 2 ),..., ψ z (u K ) as follows: which we write as where º denotes the tensor outer product, and a l and d l represent the lth column of A and D , respectively.In this way, the mixing matrix A can be estimated by solving the following problem.Given the third-order tensor M ∈ C QÂQÂK , we can compute its CAND with P components of the rank-one tensors that best approximates M , i.e., min where ‖‖ F is the Frobenius norm.Several algorithms exist for the computation of tensor decomposition.The standard way for computing the tensor decomposition is by using an ' ALS' algorithm [20].Several improved versions, such as the enhanced line search (ELS) [21] and extrapolating search direction (ESD) [22], are proposed to accelerate the rate of convergence of the ALS.Hence, the ALS is chosen here to compute the CAND.
To a large extent, the practical importance of tensor decomposition stems from its uniqueness properties.It is clear that the tensor decomposition can only be unique up to a permutation of the rank-1 terms and scaling of the factors of the rank-1 terms.Therefore, we consider the tensor decomposition (11) as essentially unique if any other matrix pair A ' and D ' , that satisfies ( 11) is related to A and D via with Δ 1 , Δ 2 ∈ C P×P being diagonal matrices, satisfying Δ 1 Δ 1 Δ 2 = I, and P ∈ R P×P being a permutation matrix.
The k-rank.The Kruskal rank or k-rank of a matrix A, denoted by κ A , is the maximal number λ such that any set of λ columns of A is linearly independent.
Theorem 1.The tensor decomposition of (11) is essentially unique if [23] We call a property generic when it holds with probability one.Generically, the mixing matrix is of full rank and of full k-rank when the parameters it involves are drawn from continuous probability densities.Hence, in practice, κ A = min(Q, P) and κ D = min(K, P).
In summary, we come to the following conclusion: when Q ≥ P, P ≥ 2, then the generic essential uniqueness is guaranteed for K ≥ 2; when Q < P and if K ≥ P, then the generic essential uniqueness is guaranteed for P ≤ 2Q − 2, if K < P, then the generic essential uniqueness is guaranteed for

The GGF-TALS algorithm
In order to ensure the GGF-ALS algorithm has robust performance, a large value is often chosen for K in the tensor M ∈ C QÂQÂK .Nevertheless, this will lead to a heavy computational load.To reduce the computational complexity, the Tucker decomposition [23][24][25] is firstly applied to represent the tensor M as a lower-rank core tensor, whose size is much smaller than the tensor M. Then the ALS algorithm is used to perform the CAND of the core tensor.In this way, the computational complexity can be reduced dramatically.
Matricization.Matricization, also known as unfolding or flatting, is the process of turning an N-way tensor into a matrix.The mode-n matricization of a tensor F ∈ C I 1 ÂI 2 ⋯ÂI N is denoted by F (n) which arranges the mode-n fibers to be the columns of the resulting matrix, that is, the tensor element (i 1 , ⋯, i N ) is mapped to the matrix element (i n , j) where The n-rank.Let F be an Nth-order tensor of size In other words, if we let r n ¼ rank n F ð Þ for n = 1,…,N, then we can say that F is a rank r 1 , ⋯, r N tensor.
The n-mode product.The n-mode (matrix) product of a tensor F ∈ C I 1 ÂI 2 ⋯ÂI N with a matrix U JÂI n is denoted by F Â n U and is of size Tucker decomposition is a form of higher-order principal component analysis (PCA).It decomposes a tensor into a core tensor multiplied by a matrix along each mode as shown in Figure 1.Thus, in the three-way case where F ∈ C I 1 ÂI 2 ÂI 3 , we have where G∈C J 1 ÂJ 2 ÂJ 3 is the core tensor; 3 are the factor matrices (which are usually column unitary for real-valued data) and can be thought of as the principal components in each mode; and U ∈ C I 1 ÂI 2 ÂI 3 represents errors or noise.
With the help of Tucker decomposition, the tensor M , which is composed of the Hessian matrices ψ z (u k ) of the second GGF of the observations, can be compressed.Since the mixing matrix A is of full rank and the Hessian matrices ψ s (A H u k ) of the second GGF of the sources are diagonal, the mode-n (n = 1, 2, 3) matrices of tensor On the other hand, if we assume rank 3 M ð Þ ¼ L and L ≤ K, then M is a rank-(Q,Q,L) tensor.Thus, the Tucker decomposition of tensor M is the so-called Tucker1 decomposition [23] where T ∈ C QÂQÂL is the core tensor, I ∈ R Q × Q is an identity matrix, and G ∈ C K × L is a column-unitary matrix.This is equivalent to a standard two-dimensional PCA since where T (3) ∈ C L × QQ is the mode-3 matrix of core tensor T .It is obvious that (18) corresponds to the PCA of M (3) .Therefore, G ∈ C K × L consists of the L-leading left singular vectors of M (3) .Since K > L and G is a column-unitary matrix, it is straightforward to derive Therefore, the core tensor can be obtained by Because the first and second factors of the Tucker decomposition in (17) are identity matrices, the core tensor T ∈ C QÂQÂL is also a symmetric tensor [23], as (1 ) (2) (3) The Tucker decomposition of a three-way tensor F .For data compression, the size of the core tensor G is smaller than that of F (I 1 ≤ J 1 , I 2 ≤ J 2 , and I 3 ≤ J 3 ), and A (1) , A (2) , and A (3) are thin matrices.
for the tensor M ∈ C QÂQÂK and the CAND of the core tensor is as follows where a l and d l are the lth column of A and D, respectively.The CAND process can also be implemented by the ALS algorithm.Nevertheless, our goal is to estimate the mixing matrix A, whereas the CAND of core tensor T ∈ C QÂQÂL only conduces A. Hence, it is necessary to derive the mixing matrix A based on A. Since the first and second factors of Tucker decomposition (17) are identity matrices, it is straightforward to derive

Computational analysis
In this subsection, we aim at giving an insight into the numerical complexity of the proposed algorithm.For the GGF-ALS algorithm, the ALS algorithm is directly used to decompose the tensor M ∈ C QÂQÂK ; therefore, the computational complexity is O(3PKQ 2 + QKP 2 + Q 2 P 2 ) per iteration.For the GGF-TALS algorithm, the computational complexity is dominated by the Tucker decomposition of the original tensor M ∈ C QÂQÂK (realized by the singularvalue decomposition (SVD) of the mode-3 matrix M (3) ) and the decomposition of the core tensor T ∈ C QÂQÂL using the ALS algorithm.Therefore, the computational complexities for these two operations per iteration are O(Q 6 ) and O(3PLQ 2 + QLP 2 + Q 2 P 2 ), respectively.Table 1 shows the computational complexity for the GGF-ALS, and GGF-TALS methods.
It is worth mentioning that a large value is usually chosen for the number of processing points K in order to accumulate sufficient statistical information from data, while the rank of the core tensor order L is often chosen to be much smaller than K.For example, K = 100, L = 8 are chosen typically in our simulations as shown in the next section.The SVD is required to be calculated only once, whereas the ALS algorithm needs multiple iterations to achieve convergence, for example, 1,000 iterations as in our experiments.Moreover, the number of sources P and the number of sensors Q are usually small, for example, P = 4, Q = 3 in our simulations.For these reasons, it can be readily derived from Table 1 that the computational cost of the GGF-TALS algorithm is, in practice, much lower than that of the GGF-ALS algorithm (note that the complexity of SVD becomes negligible in this case).

Simulations and analysis
In this section, simulations are provided to illustrate the performance of the proposed GGF approaches for underdetermined mixtures of complex sources.The performance of the tested algorithms is evaluated and compared in terms of the relative error performance index (PI) versus the sample size and the signal-to-noise ratio (SNR) of the observations.Here the relative error PI is defined as [7] PI ¼ EfjjA−A ^jj=jjAjjg, in which the norm is the Frobenius norm and Â represents the optimally ordered and scaled estimate of the mixing matrix A.
First, we compare the performance of the GGF-ALS algorithm with that of the GGF-TALS algorithm and investigate the influence of the rank of the core tensor L on the performance of the GGF-TALS algorithm.To this end, the following two simulation experiments are conducted.In the first simulation, we evaluate the performance of the GGF-TALS algorithm with different core tensor rank L for a fixed number of processing points K, and compare it with GGF-ALS with the same number of processing points.In this simulation, K is chosen to be 100 for the GGF-ALS algorithm, both the real and imaginary parts of processing points are randomly drawn from [−1; 1], the SNR of the observations ranges from 0 to 25 dB, and the number of samples is 4,000.The core tensor rank for the GGF-TALS algorithm is chosen as L = 10,8,5, respectively.The threshold value described in (12) to stop the ALS algorithm is 10 −5 , and 100 Monte Carlo experiments are run.
The average performance of the GGF-TALS algorithm versus the core tensor rank for a fixed number of processing points is shown in Figure 2. We can see that the

Method
Operation Complexity

GGF-ALS ALS (one iteration)
performance curves of the GGF-TALS and that of the GGF-ALS almost coincide when the rank of the core tensor is larger than 8, whereas the performance of the GGF-TALS algorithm slightly deteriorates when the rank of the core tensor is less than 8.This indicates that the GGF-TALS algorithm maintains the identification accuracy offered by GGF-ALS despite the fact that the core tensors used by GGF-TALS have a much lower rank than that of the original tensor used by GGF-ALS.
In the second experiment, we investigate the performance of the GGF-TALS algorithm with different number of processing points K when the rank of the core tensor L is fixed.The simulation conditions are the same as those in the previous simulation except the number of processing points for the GGF-ALS algorithm which are 20, 40, and 100 respectively, and the rank of the core tensor for the GGF-TALS which is 8.The average performance of the GGF-TALS versus the number of processing points for the fixed number of the core tensor rank is shown in Figure 3.
We can see that the average performance of GGF-TALS is consistent with that of GGF-ALS when both algorithms use the same number of processing points.For instance, when the number of processing points is 100, the performance of GGF-TALS is close to the performance of GGF-ALS, so is for K = 40.Therefore, the GGF-TALS algorithm consistently offers similar performance to GGF-ALS when the number of processing points is varied.
Second, we compare the performance of the GGF-ALS and GGF-TALS algorithms with that of the latest GF-based BI algorithm.Here, the LEMACFAC-2 in [18] is chosen as the baseline algorithm.The number of processing points and the value of processing points for the GGF-ALS, GGF-TALS, and LEMACFAC-2 are the same.Specifically, the number of processing points is 100, both the real and imaginary parts of processing points are randomly drawn from [−1; 1], and the rank of the core tensor is 8.The threshold value described in (12) to stop the LM algorithm [26] exploited in LEMACFAC-2 method is also 10 −5 .
Figure 4 shows the PI of the tested algorithms as a function of the SNR when 4,000 samples are used.It can be seen from this figure that the performance of the GGF-ALS is consistent with that of the GGF-TALS, and both perform better than LEMACFAC-2.Figure 5 shows   where diag(•) represents an operator forming a diagonal matrix by assigning its arguments to the entries of the main diagonal.

Figure 2
Figure2Performance of GGF-TALS versus core tensor order for a fixed number of processing points.The results for GGF-ALS are also shown for comparison.

Figure 3
Figure3The performance of GGF-TALS versus the number of processing points.The results for GGF-ALS are also shown for comparison.

Figure 4 1 ;Figure 5
Figure 4  The performance of the tested algorithms versus the SNR of the observations.