Abstract
The estimation of spatial signatures and spatial frequencies is crucial for several practical applications such as radar, sonar, and wireless communications. In this paper, we propose two generalized iterative estimation algorithms to the case in which a multidimensional (D) sensor array is used at the receiver. The first tensorbased algorithm is an D blind spatial signature estimator that operates in scenarios where the source’s covariance matrix is nondiagonal and unknown. The second tensorbased algorithm is formulated for the case in which the sources are uncorrelated and exploits the dualsymmetry of the covariance tensor. Additionally, a new tensorbased formulation is proposed for an shaped array configuration. Simulation results show that our proposed schemes outperform the stateoftheart matrixbased and tensorbased techniques.
1. Introduction
High resolution parameter estimation plays a fundamental role in array signal processing and has practical applications in radar, sonar, mobile communications, and seismology. In light of this, several techniques have been developed to increase the accuracy of the estimated parameters, from which we may cite the classical Multiple Signal Classification (MUSIC) [1] and Estimation of Signal Parameters via Rotational Invariance Technique (ESPRIT) [2]. However, their performance can be further improved by exploiting the multidimensional structure of the data by means of tensor modeling, which can include several signal dimensions such as space, time, frequency, and polarization. Tensor decompositions have been successfully employed in array signal processing for parameters estimation since they provide better identifiability conditions when compared to conventional matrixbased methods. Another advantage of tensorbased methods is the socalled “tensor gain” which manifests itself with more precise parameter estimates due to the good noise rejection capability of tensorbased signal processing, as shown in [3–6].
In regards to tensorbased methods for blind spatial signatures estimation, the Parallel Factor (PARAFAC) analysis decomposition [7] is widely applied due to its welldefined conditions for uniqueness [8]. As seen in [9], an iterative technique for PARAFAC decomposition such as Trilinear Alternating Least Squares (TALS) can be applied to estimate the directions of arrival of the sources. Closedform solutions such as the Standard Tensor ESPRIT (STE) [10] and ClosedForm PARAFAC [11] are also appealing, since these exploit the multidimensional structure in a noniterative fashion. Recently in [12], an iterative algorithm was proposed in a manner similar to Independent Component Analysis (ICA) based on the Orthogonal Procrustes Problem (OPP) and KhatriRao factorization [13] for a PARAFAC decomposition with dualsymmetry. This solution exploits the dualsymmetry property of the data tensor and can be applied in covariancebased array signal processing techniques. The method proposed in [14] is based on the Tucker decomposition [15] of a fourthorder covariance tensor and was elaborated for arrays with arbitrary structures, where a priori knowledge about the geometry of the sensor array is not required. However, a limitation of the method in [14] is the necessity of transmitting the same sequence of symbols in different time blocks, which results in a loss of spectral efficiency. The proposed solution is an algorithm for a multidimensional (D) sensor array in which the different dimensions of the array are exploited, thus dismissing the need to transmit a repeated sequence as in [14], as will be detailed later.
In this paper, two tensorbased approaches to the estimation of spatial signatures are presented. By using the signals received on a D sensor array, covariance tensors are calculated and solutions for correlated and uncorrelated sources are presented, respectively. For the former scenario, in which the source’s covariance structure is nondiagonal and unknown, the covariance tensor of the received data is formulated as a Tucker decomposition of order . Such a formulation yields a generalized Tucker model based D sensor array processing that deals with arbitrary source covariance structures. By assuming uncorrelated sources, we then show that the problem boils down to a PARAFAC decomposition, from which a method that exploits the dualsymmetry property of the covariance tensor is derived by considering the ideas rooted in [12]. For both Tucker and PARAFAC based models, the blind estimation of the spatial signatures is achieved by means of an alternating least squares (ALS) algorithm. The contributions of this paper are twofold: (i) we propose a covariancebased generalization of the Tucker decomposition for the blind spatial signature estimation problem with D sensor arrays and (ii) we establish a link between dualsymmetry decompositions and techniques based on covariancebased array signal processing for parameter estimation. The performance of the proposed algorithms is evaluated by Monte Carlo simulations, corroborating their gains over competing stateoftheart matrixbased and tensorbased techniques.
The rest of this paper is organized as follows: Section 2 briefly introduces tensor operations and decompositions. The signal model for an D sensor array is then presented in Section 3. In Section 4 a novel covariancebased tensor model for the received data is formulated and our blind spatial signature estimation algorithms are formulated. In Section 5 an approach for shaped sensor arrays is proposed. The computational complexity of the proposed methods is analyzed in Section 6. In Section 7, the advantages and disadvantages of the proposed methods are discussed. Simulation results are provided in Section 8, and the conclusions are drawn in Section 9.
Notation. Scalar values are represented by lowercase letters , vectors by bold lowercase letters , matrices by bold uppercase letters , and tensors by calligraphic letters . The symbols , , , and represent the transpose, conjugate transpose, pseudoinverse, and complex conjugate operations, respectively. operator generates a diagonal matrix from a vector . The th row of is represented by , while its th column is represented by . operator converts into a vector , while converts into a matrix. stands for a diagonal matrix constructed from the th row of . stands for the Frobenius norm of a matrix or tensor. “” operator stands for the vector outer product. The Kronecker product is represented by . The KhatriRao product between the matrices and , represented by , is defined as
2. Tensor Preliminaries
In the following, we briefly introduce for convenience the basics on operations involving tensors and tensor decompositions, which refer to [16, 17]. Firstly, we present the Tucker decomposition. Then, the PARAFAC decomposition is introduced and issues involving uniqueness are briefly discussed for both cases, which will be useful later. Then, we introduce the dualsymmetry property for these decompositions. The basic material presented in this section is exploited in later sections in the context of our blind spatial signature estimation problem.
2.1. Basic Tensor Operations
Let denote an th order tensor, th entry of which is denoted by . The fibers are the higherorder analogues of matrix rows and columns. The mode fibers of are vectors of size defined by fixing every index but . The mode unfolding operation, denoted by , stands for the process of reordering the elements of into a matrix by arranging its mode fibers to be the columns of the resulting matrix. The mode product between and a matrix along of the th mode, denoted by , is a tensor of size , obtained by taking the inner product between each mode fiber and the rows of the matrix ; that is, [16, 17]
2.2. Tucker Decomposition
The Tucker decomposition [15] represents a tensor of order as a multilinear transformation of a core tensor by factor matrices along each mode . In scalar form, the th order Tucker decomposition is given bywhere is th entry of the th mode factor matrix and is th entry of the core tensor . Using mode product notation, the Tucker decomposition can be written aswhich admits the following factorization in terms of the factor matrices and core tensor:
In general, the Tucker decomposition is not unique; that is, there are infinite solutions for , and that yield the same reconstructed version of the data tensor . However, in special cases where several elements of the core tensor are constrained to be equal to zero, that is, if the core tensor has some sparsity, the number of solutions may be finite, and the associated factor matrices and core tensor become unique up to trivial permutations and scaling ambiguities [18]. The Tucker based methods presented in Sections 4.2 and 5.1 belong to a special category where unique solutions exist.
2.3. PARAFAC Decomposition
The PARAFAC decomposition [7] expresses a tensor as a sum of rankone tensors; that is, where is the number of factors, also known as the rank of the decomposition, and is defined as the minimum number of rankone tensors that yield exactly.
The th order PARAFAC decomposition (6) can be seen as a special case of the Tucker decomposition (4) with a core tensor and for . The elements of the th order identity tensor are equal to one when all indices are equal and zero elsewhere. Using the mode product notation, the PARAFAC decomposition can be written aswhile the mode unfolding of can be expressed as
The th order PARAFAC decomposition is unique up to permutation and scaling ambiguities affecting the columns of factors matrices , if the following sufficient condition is satisfied [19]:where denotes the Kruskalrank of , defined as the maximum value such that any subset of columns is linearly independent [20].
Throughout this work, special attention is given to dualsymmetric tensors. The PARAFAC decomposition of a given tensor of order is said to have dualsymmetry if defined as follows:Note that this definition also applies to Tucker decomposition by simply replacing the identity tensor by an arbitrary core tensor of order .
3. Signal Model
We consider snapshots originating from the superposition of farfield narrowband signal sources sampled by a dimensional sensor array of size , where is the size of the th array dimension, . The matrix collects the samples received by the sensor array, which can be factored as [10] where (i) is the spatial signature matrix of the D array for and ;(ii) is the spatial signature matrix of the th dimension;(iii) is the array response in the th dimension to the th planar wavefront () which is function of the spatial frequency ;(iv) is the matrix containing the signal transmitted by the sources;(v) is the additive white Gaussian noise (assumed uncorrelated to the source signals).From (11), the sample covariance matrix of the signals received at the sensor array is given by where is the sample covariance matrix of the source signals and is the noise variance.
4. TensorBased Methods for Blind Spatial Signature Estimation
In this section, we propose two iterative algorithms to solve the blind spatial signature estimation problem in D sensor arrays. Initially, a novel multidimensional structure is formulated from the covariance matrix of the received data. Then, an alternating least squares (ALS) based algorithm for a Tucker decomposition of order is proposed. Finally, we derive a link between the method in [12] and a covariancebased blind spatial signature estimation problem.
4.1. Novel Covariance Tensor
With the intention of exploiting the multidimensional structure of the received signal, the noiseless sample covariance matrix (12), given by , is interpreted as a multimode unfolding of the noiseless covariance tensor of order , defined aswhere is the source covariance tensor, which has dimensions, each of size . Note that this tensor is dualsymmetric; that is, the factor matrix related to ()th dimension is equal to , and . The th frontal slice of is a diagonal matrix whose main diagonal is given by the th column of the covariance matrix . For instance, considering for the sake of notation, the following expression satisfies the relationship previously cited:where the matrix denotes the th frontal slice of the covariance tensor obtained by fixing its last two modes. The tensor follows a dualsymmetric Tucker decomposition of order with factor matrices and , and core tensor .
Considering the case in which the sources are uncorrelated and have unitary variance, we can rewrite (13) aswhere is the identity tensor of order in which each dimension has size . In this case, the covariance tensor follows a dualsymmetric PARAFAC decomposition of order .
In general, the Tucker decomposition does not impose restrictions on the core tensor structure, which makes this model more flexible. In the context of this paper, this characteristic reflects an arbitrary and unknown structure for the source’s covariance which can also be estimated from (13). In contrast, the PARAFAC decomposition (15) denotes a particular case of the Tucker decomposition when the sources’ signals are uncorrelated and the source covariance matrix is perfectly known (i.e., diagonal). In practice, this may not hold.
4.2. ALSTucker Algorithm
Our goal is to blindly estimate the spatial signature matrices and which refer to the different dimensions of the sensor array from the covariance tensor . For the sake of simplicity, from this point on, we consider . In matrixbased notation, the Tucker decomposition (13) allows the following factorization in terms of its factor matrices and core tensor in accordance with (5):wherewhile and , , denote the mode unfolding of the covariance tensor and the core tensor , respectively.
From the matrix unfoldings of , an ALS based algorithm is formulated to estimate the desired factor matrices. An estimate of the spatial signature matrix , associated with the th dimension of the covariance tensor, is obtained by solving the following least squares (LS) problem: whose analytic solution is given by
As discussed in Section 4.1, the Tucker decomposition does not impose restrictions on the structure of the core tensor and its estimation becomes necessary. Let be an unknown matrix of arbitrary structure. The following LS problem is formulated from the vectorization of the sample covariance matrix : from which an estimate of can be obtained as where .
Since (19) and (21) are nonlinear functions of the parameters to be estimated, the blind spatial signature estimation problem can be solved using a classical ALS iterative solution [21, 22]. The basic idea of the algorithm is to estimate one factor matrix at each step while the others remain fixed at the values obtained in previous steps. This procedure is repeated until convergence. The proposed generalized ALSTucker algorithm for D sensor arrays is summarized in Algorithm 1.

In this approach the factor matrices are treated as independent variables; that is, the dualsymmetry property of the covariance tensor is not exploited. In this case, a final estimate of the spatial signature matrix associated with the th dimension of the array is given by
4.3. ALSProKRaft Algorithm
In this section, a link is established between the ALSProKRaft algorithm proposed initially in [12] and blind spatial signature estimation in array signal processing. The main idea behind this algorithm is to exploit the dualsymmetry property of the PARAFAC decomposition described in (15). Next, the ALSProKRaft algorithm is formulated in the context of this work. A more detailed description of the method can be found in the original reference.
The multimode unfolding of the PARAFAC decomposition in (15) can be rewritten aswhich assumes the following factorization:where can be obtained from the singular value decomposition of given by obeying the following structure:where is formed by the first columns of and is a diagonal matrix which contains the dominant singular values of . The matrix represents a unitary rotation factor.
Equation (26) describes an orthogonal Procrustes problem (OPP) [23], in which is a transformation matrix that maps to such that . An efficient estimate for is obtained minimizing the Frobenius norm of the residual error:This problem can be solved using a change of basis from the singular value decomposition of the matrixwhich leads to the following solution [23]:
From (26) and (29), an ALSbased iterative algorithm is formulated to estimate the spatial signature matrices from the PARAFAC decomposition (15). Firstly, individual estimates for each factor matrix , , are obtained by applying the multidimensional LS KhatriRao factorization (LSKRF) algorithm on the composite spatial signature matrix . Then, the matrix is obtained from (29). For more details and access to the pseudocode of this algorithm, we refer the interested reader to [12]. The ALSProKRaft algorithm for blind spatial signature estimation in D sensor arrays is summarized in Algorithm 2.

Note that when compared with conventional ALSbased PARAFAC solutions [21] formulated from the unfolding matrices (8), the ALSProKRaft algorithm becomes preferred since only half of the factors matrices needs to be estimated by exploiting the dualsymmetry property of the covariance tensor. This generally leads to a fast convergence rate of the algorithm.
5. Spatial Signature Estimation in Shaped Sensor Arrays
In this section, the blind spatial signature estimation problem is formulated for shaped array configuration. Considering that the receiver array is divided into smaller subarrays, the Tucker decomposition of a fourthorder tensor is formulated from the sample crosscorrelation matrix of the data received by the different subarrays. From this multidimensional structure the proposed generalized ALSTucker algorithm previously presented in Section 4.2 can be used to estimate the source’s spatial signatures.
5.1. CrossCorrelation Tensor for Shaped Sensor Arrays
In this approach, we consider an shaped sensor array equipped with sensors positioned in the  plane, as illustrated in Figure 1. Each linear array contains and sensors equally spaced at distances and , respectively. We consider that each linear array is divided into and smaller subarrays, respectively. Each subarray contains and sensors. The signal received at the th subarray, for , is given byand the signal received at the th subarray, , is given bywhere (i) is the spatial signature matrix of the first subarray (or reference subarray) for the th array dimension, ;(ii) and are diagonal matrices whose main diagonal is given by the th and th row of the matrices and , respectively. The rows of and capture the delays suffered by the signals impinging the th and th subarrays with respect to the reference subarray, which are defined based on the following spatial frequencies:where and are the azimuth and elevation angles of the th source, respectively.
From (30) and (31), let us introduce the following extended data matrices:or, more compactly,
In contrast with (12), in this approach we shall work with the following sample crosscorrelation matrix:
As mentioned in Section 4.1, we can see that the noiseless term in (35) denotes a multimode unfolding of the following crosscorrelation tensor of size :where . In this modeling, the tensor follows a fourthorder Tucker decomposition. By analogy with (4), we deduce the following correspondences:
The spatial signatures of the sources can be estimated from (36) by using the proposed generalized ALSTucker algorithm. Note that, in this case, the ALSTucker algorithm is simplified to a fourthorder tensor input.
For all the previously proposed algorithms, the final estimates for the spatial signature matrices are obtained when the convergence is declared. A usually adopted criterion for convergence is defined as , where denotes the residual error of the th iteration, defined aswhere is a noisy version of , is an additive complexvalued white Gaussian noise tensor, and is the covariance tensor reconstructed from the estimated factor matrices and core tensor. Since the ALSProKRaft algorithm exploits the dualsymmetry property of the data tensor the procedure in (22) is not necessary.
5.2. Estimation of the Spatial Frequencies
After the estimation of the spatial signatures matrices , , the final step is to estimate the spatial frequencies of the sources , . The final estimates can be computed from the average over the values obtained in each row of as follows:
6. Computational Complexity
In the following, we discuss the computational complexity of the iterative ALSTucker and ALSProKRaft algorithms. For the sake of simplicity, the computational complexity of the proposed methods is described in terms of the computational cost of the matrix SVD. For a matrix of size the number of floatingpoint operations associated with the SVD computation is [24]. The computational complexity of one Tucker iteration refers to the cost associated with the SVD used to calculate the matrix pseudoinverses in the least squares problems (18) and (20). The overall computational complexity per iteration of the ALSTucker algorithm equals the complexity of matrix SVDs associated with each estimated factor matrix according to (19) plus the complexity of one additional matrix SVD associated with the estimated core tensor according to (21). The overall computational cost per iteration of the ALSProKRaft algorithm equals the complexity of matrix SVDs associated with the application of the multidimensional LSKRF algorithm in (26) plus the complexity of one additional matrix SVD associated with the update of the unknown unitary rotation factor matrix according to (29).
7. Advantages and Disadvantages of the Proposed Methods
In this section, we discuss the advantages and disadvantages of the proposed methods to blind spatial signatures estimation in D sensor arrays. As previously stated in Section 4.3, the ALSProKRaft algorithm works on the assumption that the sample covariance matrix of the sources is perfectly known and diagonal. However, this is only true in the asymptotic case when a sufficiently large number of snapshots is assumed (i.e., ), as well as when the source signals are perfectly uncorrelated. In practice, this assumption is not guaranteed. On the other hand, the ALSTucker algorithm previously formulated in Section 4.2 naturally captures any structure for the sources covariance into the core tensor . Therefore, the assumption of uncorrelated source signals is not necessary for the ALSTucker algorithm, making it able to operate in scenarios where the source covariance structure is unknown and arbitrary (nondiagonal). Such scenarios occur, for instance, when the sample covariance is computed from a limited number of snapshots.
In contrast, the ALSTucker algorithm does not exploit the dualsymmetry property of the data covariance tensor and all factor matrices need to be estimated as independent variables. However, in the ALSProKRaft algorithm, only half of the factor matrices are estimated by exploiting the dualsymmetry property of the covariance tensor. Therefore, the ALSProKRaft algorithm is more computationally attractive than the ALSTucker algorithm. When compared with the stateoftheart matrixbased algorithms such as MUSIC, ESPRIT, and Propagator Method, the proposed tensorbased algorithms have the advantage of fully exploiting the multidimensional nature of the received signal in less specific scenarios, which leads to more accurate estimates. For instance, the ESPRIT algorithm was formulated for sensor arrays that obey the shift invariance property. On the other hand, the MUSIC algorithm has high computational complexity due to the search of parameters in the spatial spectrum.
8. Simulation Results
In the following, simulation results and performance evaluations of the ALSTucker and ALSProKRaft algorithms for D sensor arrays are presented. This section is divided into two parts. Firstly, the results related to Section 4 are presented and discussed. Then, the same is done for the shaped array approach presented in Section 5. Results are obtained from an average of 1000 independent Monte Carlo runs. In the first part of this section, we consider a uniform rectangular array (URA) positioned on the  plane. The th wavefront has direction of arrival , where and are elevation and azimuth angles, respectively.
In Figures 2 and 3, the performance is measured in terms of the root mean square error (RMSE) of the estimated spatial frequencies in terms of Signal to Noise Ratio (SNR). The relations between directions of arrival and spatial frequencies are given by (32), where denotes the distance between sensors in the th array dimension, which is assumed here equal to . In Figure 2, we consider Hadamard sequences for the sources, while in Figure 3 we consider BPSK modulated sources. In both cases, we have sensors (i.e., and equal to 8) and the sample covariance matrix of the received data (12) is calculated from a reduced number of samples. The total RMSE is defined as
From Figure 2, it can be seen that both ALSTucker and ALSProKRaft algorithms have similar performances in terms of RMSE, when Hadamard sequences are considered. In contrast, in Figure 3, when BPSK symbols are considered, a floor is exhibited by the ALSProKraft algorithm for high SNR values. This behavior occurs due to the modeling errors in the core tensor of the PARAFAC decomposition, which in turn arises due to the nonorthogonality of the source signals, resulting in a nondiagonal sample covariance matrix of the sources in this case. Note that in the ALSTucker algorithm the covariance matrix of the sources possess an arbitrary and unknown structure which discards possible constraints to the source signals. This difference makes the Tucker decomposition approach more attractive in those practical scenarios in which source uncorrelatedness is not guaranteed. When compared to matrixbased standard ESPRIT (SE) [2], matrixbased MUSIC algorithm to planar array configuration [25], and tensorbased standard ESPRIT (STE) [10], the proposed algorithms have improved accuracy in all the considered scenarios.
Figure 4 shows the convergence performance of the iterative algorithms. In this experiment, the median values of the normalized estimation error are plotted in terms of the number of iterations for different SNR. It is noteworthy that the ALSProKRaft algorithm has a faster convergence compared to the ALSTucker algorithm. This behavior is expected since ALSProKRaft exploits the dualsymmetry property of the data tensor, which results in estimating half as many factor matrices compared to the ALSTucker approach.
In the second part of this section, we consider a shaped configuration array. In Figure 5, we set sensors (i.e., and equal to 7) and samples. Each uniform linear array is divided into and subarrays, respectively. In this experiment, the performance of the proposed ALSTucker algorithm is compared to the stateoftheart matrixbased methods, namely, Propagator Method (PM) [26], MUSIC [27], and ESPRIT [28], all of them originally formulated for shaped arrays. Note that the ALSTucker algorithm presents an improved performance over its competitors, with more evidenced gains in the lowtomedium SNR range. For high SNR values, the performance of the MUSIC method comes close to that of our proposal. However, the ALSTucker algorithm dispenses any estimation procedure via bidimensional peak search as occurs with MUSIC, being the former more computationally attractive.
Figure 6 shows the performance of the ALSTucker by assuming different number of sensors. In this experiment, we consider samples. We can observe a better performance in terms of RMSE when the number sensors of the shaped array is increased. This is valid for all the simulated SNR values.
In Figure 7, we analyze the influence of the number of samples on the performance of the ALSTucker algorithm. This experiment considers the same parameters as the experiment of Figure 5, except the SNR value that is assumed fixed at 20 dB and the number of samples that varies between 50 and 3000. First, we can see that the performance of the algorithms improves by increasing the number of samples collected by sensor array, as expected. However, similar to Figure 5 the proposed ALSTucker algorithm outperforms the stateoftheart PM, ESPRIT, and MUSIC methods.
9. Conclusion
In this paper, two tensorbased approaches based on the Tucker and PARAFAC decompositions have been formulated to solve the blind spatial signatures estimation problem in multidimensional sensor arrays. First, we have proposed a covariancebased generalization of the Tucker decomposition for D sensor arrays. Then, a link between the ALSProKRaft algorithm and covariancebased array signal processing for blind spatial signatures estimation has been established. As another contribution, we have formulated a crosscorrelationbased fourthorder Tucker decomposition which makes the proposed ALSTucker algorithm applicable in scenarios composed by shaped array configurations. The two proposed tensor methods differ in the structure assumed for the source covariance. It is worth pointing out that, in realistic scenarios, when the received covariance matrix is calculated from a reduced number of samples, or snapshots, the ALSTucker algorithm becomes preferable since it operates with an arbitrary and unknown structure for the covariance of the source signals. In contrast, when the sources can be assumed to be uncorrelated, we can achieve improved performance by exploiting the dualsymmetry property of the covariance tensor, which makes the ALSProKRaft algorithm preferable since it provides good estimation accuracy with a smaller number of ALS iterations. Finally, when compared with other stateoftheart matrixbased and tensorbased techniques, the proposed tensorbased iterative algorithms have shown their effectiveness with remarkable gains in terms of estimation error.
Competing Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors would like to thank the Fundação Cearense de Apoio ao Desenvolvimento Científico e Tecnológico (FUNCAP), the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) under the PVE Grant no. 88881.030392/201301, the productivity Grant no. 303905/20140, the postdoctoral scholarship abroad (PDE) no. 207644/20152, and the Program Science without Borders, Aerospace Technology supported by CNPq and CAPES for the postdoctoral scolarship abroad (PDE) no. 207644/20152.