# Tensor-Based Methods for Blind Spatial Signature Estimation in Multidimensional Sensor Arrays

Authors: Paulo R. B. Gomes; André L. F. de Almeida; João Paulo C. L. da Costa; João C. M. Mota; Daniel Valle de Lima; Giovanni Del Galdo
Journal: International Journal of Antennas and Propagation
 (2017)

---

## 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 (R-D) sensor array is used at the receiver. The first tensor-based algorithm is an R-D blind spatial signature estimator that operates in scenarios where the source’s covariance matrix is nondiagonal and unknown. The second tensor-based algorithm is formulated for the case in which the sources are uncorrelated and exploits the dual-symmetry of the covariance tensor. Additionally, a new tensor-based formulation is proposed for an L-shaped array configuration. Simulation results show that our proposed schemes outperform the state-of-the-art matrix-based and tensor-based techniques.


---

## Body

## 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 classicalMultiple Signal Classification (MUSIC) [1] andEstimation 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 matrix-based methods. Another advantage of tensor-based methods is the so-called “tensor gain” which manifests itself with more precise parameter estimates due to the good noise rejection capability of tensor-based signal processing, as shown in [3–6].In regards to tensor-based methods for blind spatial signatures estimation, theParallel Factor(PARAFAC) analysis decomposition [7] is widely applied due to its well-defined conditions for uniqueness [8]. As seen in [9], an iterative technique for PARAFAC decomposition such asTrilinear Alternating Least Squares (TALS) can be applied to estimate the directions of arrival of the sources. Closed-form solutions such as the Standard Tensor ESPRIT (STE) [10] and Closed-Form 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 Khatri-Rao factorization [13] for a PARAFAC decomposition with dual-symmetry. This solution exploits the dual-symmetry property of the data tensor and can be applied in covariance-based array signal processing techniques. The method proposed in [14] is based on the Tucker decomposition [15] of a fourth-order 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 (R-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 tensor-based approaches to the estimation of spatial signatures are presented. By using the signals received on aR-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 2R. Such a formulation yields a generalized Tucker model based R-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 dual-symmetry 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 covariance-based generalization of the Tucker decomposition for the blind spatial signature estimation problem with R-D sensor arrays and (ii) we establish a link between dual-symmetry decompositions and techniques based on covariance-based array signal processing for parameter estimation. The performance of the proposed algorithms is evaluated by Monte Carlo simulations, corroborating their gains over competing state-of-the-art matrix-based and tensor-based techniques.The rest of this paper is organized as follows: Section2 briefly introduces tensor operations and decompositions. The signal model for an R-D sensor array is then presented in Section 3. In Section 4 a novel covariance-based tensor model for the received data is formulated and our blind spatial signature estimation algorithms are formulated. In Section 5 an approach for L-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 a, vectors by bold lowercase letters a, matrices by bold uppercase letters A, and tensors by calligraphic letters A. The symbols T, H, †, and ∗ represent the transpose, conjugate transpose, pseudoinverse, and complex conjugate operations, respectively. diaga operator generates a diagonal matrix from a vector a. The ith row of A∈CI×R is represented by Ai,:∈C1×R, while its rth column is represented by A:,r∈CI×1. vecA operator converts A into a vector a∈CIR×1, while unvecI×Ra converts a into a I×R matrix. DiA stands for a diagonal matrix constructed from the ith row of A. ·F stands for the Frobenius norm of a matrix or tensor. “∘” operator stands for the vector outer product. The Kronecker product is represented by ⊗. The Khatri-Rao product between the matrices A∈CI×R and B∈CJ×R, represented by ⋄, is defined as(1)A⋄B=A:,1⊗B:,1,…,A:,R⊗B:,R.


## 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 dual-symmetry 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

LetX∈CI1×⋯×IN denote an Nth order tensor, i1,…,iNth entry of which is denoted by xi1,…,iN. The fibers are the higher-order analogues of matrix rows and columns. The n-mode fibers of X are vectors of size In defined by fixing every index but in. The n-mode unfolding operation, denoted by Xn, stands for the process of reordering the elements of X into a matrix by arranging its n-mode fibers to be the columns of the resulting matrix. The n-mode product between X and a matrix A∈CJ×In along of the nth mode, denoted by X×nA, is a tensor of size I1×⋯×In-1×J×In+1×⋯×IN, obtained by taking the inner product between each n-mode fiber and the rows of the matrix A; that is, [16, 17](2)Y=X×nA⟺Yn=AXn.


### 2.2. Tucker Decomposition

The Tucker decomposition [15] represents a tensor X∈CI1×⋯×IN of order N as a multilinear transformation of a core tensor G∈CR1×⋯×RN by factor matrices An=a1n,a2n,…,aRnn∈CIn×Rn along each mode n=1,2,…,N. In scalar form, the Nth order Tucker decomposition is given by(3)xi1,i2,…,iN=∑r1=1R1∑r2=1R2⋯∑rN=1RNgr1,r2,…,rNai1,r11ai2,r22⋯aiN,rNN,where ain,rnn is in,rnth entry of the nth mode factor matrix An∈CIn×Rn and gr1,r2,…,rN is r1,…,rNth entry of the core tensor G∈CR1×⋯×RN. Using n-mode product notation, the Tucker decomposition can be written as(4)X=G×1A1×2A2⋯×NAN,which admits the following factorization in terms of the factor matrices and core tensor:(5)Xn=AnGnAN⊗⋯An+1⊗An-1⊗⋯A1T.In general, the Tucker decomposition is not unique; that is, there are infinite solutions forAn,n=1,…,N, and G that yield the same reconstructed version of the data tensor X. 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 X∈CI1×⋯×IN as a sum of R rank-one tensors; that is,(6)X=∑r=1Rar1∘ar2∘⋯∘arN, where R is the number of factors, also known as therank of the decomposition, and is defined as the minimum number of rank-one tensors that yield X exactly.TheNth order PARAFAC decomposition (6) can be seen as a special case of the Tucker decomposition (4) with a core tensor G=IN,R and Rn=R for n=1,…,N. The elements of the Nth order identity tensor IN,R are equal to one when all indices are equal and zero elsewhere. Using the n-mode product notation, the PARAFAC decomposition can be written as(7)X=IN,R×1A1×2A2⋯×NAN,while the n-mode unfolding of X can be expressed as(8)Xn=AnAN⋄⋯An+1⋄An-1⋄⋯A1T.TheNth order PARAFAC decomposition is unique up to permutation and scaling ambiguities affecting the columns of factors matrices An∈CIn×R,n=1,…,N, if the following sufficient condition is satisfied [19]:(9)∑n=1NκAn≥2R+N-1,where κAn denotes the Kruskal-rank of An, defined as the maximum value κ such that any subset of κ columns is linearly independent [20].Throughout this work, special attention is given to dual-symmetric tensors. The PARAFAC decomposition of a given tensorX∈CI1×⋯×I2N of order 2N is said to have dual-symmetry if defined as follows:(10)X=I2N,R×1A1×2A2⋯×NAN×N+1A1∗×N+2A2∗⋯×2NAN∗.Note that this definition also applies to Tucker decomposition by simply replacing the identity tensor I2N,R by an arbitrary core tensor G of order 2N.


## 2.1. Basic Tensor Operations

LetX∈CI1×⋯×IN denote an Nth order tensor, i1,…,iNth entry of which is denoted by xi1,…,iN. The fibers are the higher-order analogues of matrix rows and columns. The n-mode fibers of X are vectors of size In defined by fixing every index but in. The n-mode unfolding operation, denoted by Xn, stands for the process of reordering the elements of X into a matrix by arranging its n-mode fibers to be the columns of the resulting matrix. The n-mode product between X and a matrix A∈CJ×In along of the nth mode, denoted by X×nA, is a tensor of size I1×⋯×In-1×J×In+1×⋯×IN, obtained by taking the inner product between each n-mode fiber and the rows of the matrix A; that is, [16, 17](2)Y=X×nA⟺Yn=AXn.


## 2.2. Tucker Decomposition

The Tucker decomposition [15] represents a tensor X∈CI1×⋯×IN of order N as a multilinear transformation of a core tensor G∈CR1×⋯×RN by factor matrices An=a1n,a2n,…,aRnn∈CIn×Rn along each mode n=1,2,…,N. In scalar form, the Nth order Tucker decomposition is given by(3)xi1,i2,…,iN=∑r1=1R1∑r2=1R2⋯∑rN=1RNgr1,r2,…,rNai1,r11ai2,r22⋯aiN,rNN,where ain,rnn is in,rnth entry of the nth mode factor matrix An∈CIn×Rn and gr1,r2,…,rN is r1,…,rNth entry of the core tensor G∈CR1×⋯×RN. Using n-mode product notation, the Tucker decomposition can be written as(4)X=G×1A1×2A2⋯×NAN,which admits the following factorization in terms of the factor matrices and core tensor:(5)Xn=AnGnAN⊗⋯An+1⊗An-1⊗⋯A1T.In general, the Tucker decomposition is not unique; that is, there are infinite solutions forAn,n=1,…,N, and G that yield the same reconstructed version of the data tensor X. 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 X∈CI1×⋯×IN as a sum of R rank-one tensors; that is,(6)X=∑r=1Rar1∘ar2∘⋯∘arN, where R is the number of factors, also known as therank of the decomposition, and is defined as the minimum number of rank-one tensors that yield X exactly.TheNth order PARAFAC decomposition (6) can be seen as a special case of the Tucker decomposition (4) with a core tensor G=IN,R and Rn=R for n=1,…,N. The elements of the Nth order identity tensor IN,R are equal to one when all indices are equal and zero elsewhere. Using the n-mode product notation, the PARAFAC decomposition can be written as(7)X=IN,R×1A1×2A2⋯×NAN,while the n-mode unfolding of X can be expressed as(8)Xn=AnAN⋄⋯An+1⋄An-1⋄⋯A1T.TheNth order PARAFAC decomposition is unique up to permutation and scaling ambiguities affecting the columns of factors matrices An∈CIn×R,n=1,…,N, if the following sufficient condition is satisfied [19]:(9)∑n=1NκAn≥2R+N-1,where κAn denotes the Kruskal-rank of An, defined as the maximum value κ such that any subset of κ columns is linearly independent [20].Throughout this work, special attention is given to dual-symmetric tensors. The PARAFAC decomposition of a given tensorX∈CI1×⋯×I2N of order 2N is said to have dual-symmetry if defined as follows:(10)X=I2N,R×1A1×2A2⋯×NAN×N+1A1∗×N+2A2∗⋯×2NAN∗.Note that this definition also applies to Tucker decomposition by simply replacing the identity tensor I2N,R by an arbitrary core tensor G of order 2N.


## 3. Signal Model

We considerK snapshots originating from the superposition of M far-field narrowband signal sources sampled by a R-dimensional sensor array of size N1×N2×⋯×NR, where Nr is the size of the rth array dimension, r=1,…,R. The matrix X∈CN×K collects the samples received by the sensor array, which can be factored as [10](11)X=A1⋄A2⋯⋄ARS+V, where (i)
A=A1⋄A2⋯⋄AR∈CN×M is the spatial signature matrix of the R-D array for r=1,…,R and N=∏r=1R‍Nr;(ii)
Ar=a1r,…,aMr∈CNr×M is the spatial signature matrix of the rth dimension;(iii)
amr=1ej·μmrej·2μmr⋯ej·Nr-1μmrT∈CNr×1 is the array response in the rth dimension to the mth planar wavefront (m=1,…,M) which is function of the spatial frequency μmr;(iv)
S=s1,…,sK∈CM×K is the matrix containing the signal transmitted by the sources;(v)
V=v1,…,vK∈CN×K is the additive white Gaussian noise (assumed uncorrelated to the source signals).From (11), the sample covariance matrix R^∈CN×N of the signals received at the sensor array is given by (12)R^≜1KXXH≈ARsAH+σv2I≈A1⋄⋯⋄ARRsA1⋄⋯⋄ARH+σv2I,where Rs=1/KSSH is the sample covariance matrix of the source signals and σv2 is the noise variance.


## 4. Tensor-Based Methods for Blind Spatial Signature Estimation

In this section, we propose two iterative algorithms to solve the blind spatial signature estimation problem inR-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 2N is proposed. Finally, we derive a link between the method in [12] and a covariance-based 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 Ro=R^-σv2I∈CN×N, is interpreted as a multimode unfolding of the noiseless covariance tensor Ro∈CN1×N2×⋯×NR×N1×N2×⋯×NR of order 2N, defined as(13)Ro=Rs×1A1×2A2⋯×RAR×R+1A1∗×R+2A2∗⋯×2RAR∗,where Rs is the source covariance tensor, which has 2R dimensions, each of size M. Note that this tensor is dual-symmetric; that is, the factor matrix related to (R+r)th dimension is equal to Ar∗, and Nr=NR+rr=1,…,R. The mth frontal slice of Rs is a diagonal matrix whose main diagonal is given by the mth column of the covariance matrix Rs. For instance, considering R=2 for the sake of notation, the following expression satisfies the relationship previously cited:(14)Rs:,:,m,m=diagRs:,mm=1,…,M,where the matrix Rs:,:,m,m∈CM×M denotes the mth frontal slice of the covariance tensor Rs obtained by fixing its last two modes. The tensor Ro follows a dual-symmetric Tucker decomposition of order 2R with factor matrices Ar and Ar∗,r=1,…,R, and core tensor Rs.Considering the case in which the sources are uncorrelated and have unitary variance, we can rewrite (13) as(15)Ro=I2R,M×1A1×2A2⋯×RAR×R+1A1∗×R+2A2∗⋯×2RAR∗,where I2R,M is the identity tensor of order 2R in which each dimension has size M. In this case, the covariance tensor Ro follows a dual-symmetric PARAFAC decomposition of order 2R.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 covarianceRs 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. ALS-Tucker Algorithm

Our goal is to blindly estimate the spatial signature matricesAr and Ar∗,r=1,…,R which refer to the different dimensions of the sensor array from the covariance tensor Ro. For the sake of simplicity, from this point on, we consider A(R+r)=Ar∗. In matrix-based notation, the Tucker decomposition (13) allows the following factorization in terms of its factor matrices and core tensor in accordance with (5):(16)Ror=ArΔr,where(17)Δr=RsrA2R⊗⋯⊗Ar+1⊗Ar-1⊗⋯⊗A1T,while Ror and Rsr,  r=1,…,2R, denote the n-mode unfolding of the covariance tensor Ro and the core tensor Rs, respectively.From the matrix unfoldings ofRo, an ALS based algorithm is formulated to estimate the desired factor matrices. An estimate of the spatial signature matrix A^rr=1,…,2R, associated with the rth dimension of the covariance tensor, is obtained by solving the following least squares (LS) problem:(18)A^r=argminArRor-ArΔrF2, whose analytic solution is given by(19)A^r=RorΔr†.As discussed in Section4.1, the Tucker decomposition does not impose restrictions on the structure of the core tensor and its estimation becomes necessary. Let Rs be an unknown matrix of arbitrary structure. The following LS problem is formulated from the vectorization of the sample covariance matrix R^:(20)vecR^s=argminRsvecR^-A∗⊗AvecRsF2, from which an estimate of R^s can be obtained as(21)vecR^s=A∗⊗A†vecR^, where A=A1⋄A2⋯⋄AR∈CN×M.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 ALS-Tucker algorithm for R-D sensor arrays is summarized in Algorithm 1.Algorithm 1:Summary of the ALS-Tucker algorithm.
(1) Initialize A^(r)∈CNr×M for r=2,…,2R and the core tensorR^s randomly. (2) for  r=1,…,2R  do According to (19), obtain an estimate for the matrix A^(r): A^(r)=R(r)Δ(r)†. Note. The matrix Δr, described in (17), is updated by fixingA^r calculated previously. end (3) According to (21), obtain an estimate for R^s: vecR^s=A∗⊗A†vecR^; R^s=unvecM×MvecR^s. (4)  Using (14), reconstruct the core tensor R^s from R^s. (5)  Repeat Steps (2)–(4) until convergence.In this approach the factor matrices are treated as independent variables; that is, the dual-symmetry property of the covariance tensor is not exploited. In this case, a final estimate of the spatial signature matrix associated with therth dimension of the array is given by(22)A^finalr=A^r+A^R+r∗2.


### 4.3. ALS-ProKRaft Algorithm

In this section, a link is established between the ALS-ProKRaft algorithm proposed initially in [12] and blind spatial signature estimation in array signal processing. The main idea behind this algorithm is to exploit the dual-symmetry property of the PARAFAC decomposition described in (15). Next, the ALS-ProKRaft 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 as(23)Rmm=A1⋄⋯⋄ARA1⋄⋯⋄ARH,which assumes the following factorization:(24)Rmm=Rmm1/2·Rmm1/2H,where Rmm1/2∈CN×M can be obtained from the singular value decomposition of Rmm given by(25)Rmm=U·Σ·VH, obeying the following structure:(26)Rmm1/2=UM·ΣM·TH=A1⋄⋯⋄AR,where U[M]∈CN×M is formed by the first M columns of U and ΣM∈CM×M is a diagonal matrix which contains the M dominant singular values of Rmm. The matrix T represents a unitary rotation factor.Equation (26) describes an orthogonal Procrustes problem (OPP) [23], in which T is a transformation matrix that maps UM·ΣM to A1⋄⋯⋄AR such that UM·Σ[M]·TH=A1⋄⋯⋄AR. An efficient estimate for T is obtained minimizing the Frobenius norm of the residual error:(27)argminTUM·ΣM·TH-A1⋄⋯⋄ARF.This problem can be solved using a change of basis from the singular value decomposition of the matrix(28)A1⋄⋯⋄ARH·UM·ΣM=UP·ΣP·VPH,which leads to the following solution [23]:(29)T^=UP·VPH.From (26) and (29), an ALS-based iterative algorithm is formulated to estimate the spatial signature matrices from the PARAFAC decomposition (15). Firstly, individual estimates for each factor matrix A^r,  r=1,…,R, are obtained by applying the multidimensional LS Khatri-Rao factorization (LS-KRF) algorithm on the composite spatial signature matrix A^=A^1⋄A^2⋯⋄A^R. Then, the matrix T^ is obtained from (29). For more details and access to the pseudocode of this algorithm, we refer the interested reader to [12]. The ALS-ProKRaft algorithm for blind spatial signature estimation in R-D sensor arrays is summarized in Algorithm 2.Algorithm 2:Summary of the ALS-ProKRaft algorithm.
i = 0; Initialize T^i=0=IM. (1)  According to (25), obtain UM and ΣM from the SVD of the multimode unfolding matrixRmm. (2)  i=i+1; (3)  According to (26), obtain estimates for A^(i)(r) for r=1,…,R by applying the multidimensional LS-KRF algorithm on UM⋅ΣM⋅T^i-1H. (4)  According to (29), compute the SVD for the matrix A^i1⋄⋯⋄A^iRH⋅UM⋅ΣM=UP⋅ΣP⋅VPH and obtainT^(i)=UP⋅VPH. (5)  Repeat Steps (2)–(4) until convergence.Note that when compared with conventional ALS-based PARAFAC solutions [21] formulated from the unfolding matrices (8), the ALS-ProKRaft algorithm becomes preferred since only half of the factors matrices needs to be estimated by exploiting the dual-symmetry property of the covariance tensor. This generally leads to a fast convergence rate of the algorithm.


## 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 Ro=R^-σv2I∈CN×N, is interpreted as a multimode unfolding of the noiseless covariance tensor Ro∈CN1×N2×⋯×NR×N1×N2×⋯×NR of order 2N, defined as(13)Ro=Rs×1A1×2A2⋯×RAR×R+1A1∗×R+2A2∗⋯×2RAR∗,where Rs is the source covariance tensor, which has 2R dimensions, each of size M. Note that this tensor is dual-symmetric; that is, the factor matrix related to (R+r)th dimension is equal to Ar∗, and Nr=NR+rr=1,…,R. The mth frontal slice of Rs is a diagonal matrix whose main diagonal is given by the mth column of the covariance matrix Rs. For instance, considering R=2 for the sake of notation, the following expression satisfies the relationship previously cited:(14)Rs:,:,m,m=diagRs:,mm=1,…,M,where the matrix Rs:,:,m,m∈CM×M denotes the mth frontal slice of the covariance tensor Rs obtained by fixing its last two modes. The tensor Ro follows a dual-symmetric Tucker decomposition of order 2R with factor matrices Ar and Ar∗,r=1,…,R, and core tensor Rs.Considering the case in which the sources are uncorrelated and have unitary variance, we can rewrite (13) as(15)Ro=I2R,M×1A1×2A2⋯×RAR×R+1A1∗×R+2A2∗⋯×2RAR∗,where I2R,M is the identity tensor of order 2R in which each dimension has size M. In this case, the covariance tensor Ro follows a dual-symmetric PARAFAC decomposition of order 2R.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 covarianceRs 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. ALS-Tucker Algorithm

Our goal is to blindly estimate the spatial signature matricesAr and Ar∗,r=1,…,R which refer to the different dimensions of the sensor array from the covariance tensor Ro. For the sake of simplicity, from this point on, we consider A(R+r)=Ar∗. In matrix-based notation, the Tucker decomposition (13) allows the following factorization in terms of its factor matrices and core tensor in accordance with (5):(16)Ror=ArΔr,where(17)Δr=RsrA2R⊗⋯⊗Ar+1⊗Ar-1⊗⋯⊗A1T,while Ror and Rsr,  r=1,…,2R, denote the n-mode unfolding of the covariance tensor Ro and the core tensor Rs, respectively.From the matrix unfoldings ofRo, an ALS based algorithm is formulated to estimate the desired factor matrices. An estimate of the spatial signature matrix A^rr=1,…,2R, associated with the rth dimension of the covariance tensor, is obtained by solving the following least squares (LS) problem:(18)A^r=argminArRor-ArΔrF2, whose analytic solution is given by(19)A^r=RorΔr†.As discussed in Section4.1, the Tucker decomposition does not impose restrictions on the structure of the core tensor and its estimation becomes necessary. Let Rs be an unknown matrix of arbitrary structure. The following LS problem is formulated from the vectorization of the sample covariance matrix R^:(20)vecR^s=argminRsvecR^-A∗⊗AvecRsF2, from which an estimate of R^s can be obtained as(21)vecR^s=A∗⊗A†vecR^, where A=A1⋄A2⋯⋄AR∈CN×M.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 ALS-Tucker algorithm for R-D sensor arrays is summarized in Algorithm 1.Algorithm 1:Summary of the ALS-Tucker algorithm.
(1) Initialize A^(r)∈CNr×M for r=2,…,2R and the core tensorR^s randomly. (2) for  r=1,…,2R  do According to (19), obtain an estimate for the matrix A^(r): A^(r)=R(r)Δ(r)†. Note. The matrix Δr, described in (17), is updated by fixingA^r calculated previously. end (3) According to (21), obtain an estimate for R^s: vecR^s=A∗⊗A†vecR^; R^s=unvecM×MvecR^s. (4)  Using (14), reconstruct the core tensor R^s from R^s. (5)  Repeat Steps (2)–(4) until convergence.In this approach the factor matrices are treated as independent variables; that is, the dual-symmetry property of the covariance tensor is not exploited. In this case, a final estimate of the spatial signature matrix associated with therth dimension of the array is given by(22)A^finalr=A^r+A^R+r∗2.


## 4.3. ALS-ProKRaft Algorithm

In this section, a link is established between the ALS-ProKRaft algorithm proposed initially in [12] and blind spatial signature estimation in array signal processing. The main idea behind this algorithm is to exploit the dual-symmetry property of the PARAFAC decomposition described in (15). Next, the ALS-ProKRaft 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 as(23)Rmm=A1⋄⋯⋄ARA1⋄⋯⋄ARH,which assumes the following factorization:(24)Rmm=Rmm1/2·Rmm1/2H,where Rmm1/2∈CN×M can be obtained from the singular value decomposition of Rmm given by(25)Rmm=U·Σ·VH, obeying the following structure:(26)Rmm1/2=UM·ΣM·TH=A1⋄⋯⋄AR,where U[M]∈CN×M is formed by the first M columns of U and ΣM∈CM×M is a diagonal matrix which contains the M dominant singular values of Rmm. The matrix T represents a unitary rotation factor.Equation (26) describes an orthogonal Procrustes problem (OPP) [23], in which T is a transformation matrix that maps UM·ΣM to A1⋄⋯⋄AR such that UM·Σ[M]·TH=A1⋄⋯⋄AR. An efficient estimate for T is obtained minimizing the Frobenius norm of the residual error:(27)argminTUM·ΣM·TH-A1⋄⋯⋄ARF.This problem can be solved using a change of basis from the singular value decomposition of the matrix(28)A1⋄⋯⋄ARH·UM·ΣM=UP·ΣP·VPH,which leads to the following solution [23]:(29)T^=UP·VPH.From (26) and (29), an ALS-based iterative algorithm is formulated to estimate the spatial signature matrices from the PARAFAC decomposition (15). Firstly, individual estimates for each factor matrix A^r,  r=1,…,R, are obtained by applying the multidimensional LS Khatri-Rao factorization (LS-KRF) algorithm on the composite spatial signature matrix A^=A^1⋄A^2⋯⋄A^R. Then, the matrix T^ is obtained from (29). For more details and access to the pseudocode of this algorithm, we refer the interested reader to [12]. The ALS-ProKRaft algorithm for blind spatial signature estimation in R-D sensor arrays is summarized in Algorithm 2.Algorithm 2:Summary of the ALS-ProKRaft algorithm.
i = 0; Initialize T^i=0=IM. (1)  According to (25), obtain UM and ΣM from the SVD of the multimode unfolding matrixRmm. (2)  i=i+1; (3)  According to (26), obtain estimates for A^(i)(r) for r=1,…,R by applying the multidimensional LS-KRF algorithm on UM⋅ΣM⋅T^i-1H. (4)  According to (29), compute the SVD for the matrix A^i1⋄⋯⋄A^iRH⋅UM⋅ΣM=UP⋅ΣP⋅VPH and obtainT^(i)=UP⋅VPH. (5)  Repeat Steps (2)–(4) until convergence.Note that when compared with conventional ALS-based PARAFAC solutions [21] formulated from the unfolding matrices (8), the ALS-ProKRaft algorithm becomes preferred since only half of the factors matrices needs to be estimated by exploiting the dual-symmetry property of the covariance tensor. This generally leads to a fast convergence rate of the algorithm.


## 5. Spatial Signature Estimation inL-Shaped Sensor Arrays

In this section, the blind spatial signature estimation problem is formulated forL-shaped array configuration. Considering that the receiver array is divided into smaller subarrays, the Tucker decomposition of a fourth-order tensor is formulated from the sample cross-correlation matrix of the data received by the different subarrays. From this multidimensional structure the proposed generalized ALS-Tucker algorithm previously presented in Section 4.2 can be used to estimate the source’s spatial signatures.


### 5.1. Cross-Correlation Tensor forL-Shaped Sensor Arrays

In this approach, we consider anL-shaped sensor array equipped with N1+N2-1 sensors positioned in the x-z plane, as illustrated in Figure 1. Each linear array contains N1 and N2 sensors equally spaced at distances d1 and d2, respectively. We consider that each linear array is divided into P and W smaller subarrays, respectively. Each subarray contains N1sub=N1-P+1 and N2sub=N2-W+1 sensors. The signal received at the pth subarray, for p=1,…,P, is given by(30)Xp1=As1DpΦ1S+Vp1∈CN1sub×K,and the signal received at the wth subarray, w=1,…,W, is given by(31)Xw2=As2DwΦ2S+Vw2∈CN2sub×K,where (i)
Asr∈CNrsub×M is the spatial signature matrix of the first subarray (or reference subarray) for the rth array dimension, r=1,2;(ii)
DpΦ1 and DwΦ2 are diagonal matrices whose main diagonal is given by the pth and wth row of the matrices Φ(1)∈CP×M and Φ(2)∈CW×M, respectively. The rows of Φ1 and Φ2 capture the delays suffered by the signals impinging the pth and wth subarrays with respect to the reference subarray, which are defined based on the following spatial frequencies:(32)μm1=2π·d1·cos⁡αmλ,μm2=2π·d2·sin⁡αm·cos⁡βmλ,where αm and βm are the azimuth and elevation angles of the mth source, respectively.Figure 1
L-shaped array configuration with N1+N2-1 sensors. The distance between the sensors in the z-axis is d1 while the distance between the sensors in the x-axis is d2. The mth wavefront has elevation and azimuth angles equal to αm and βm, respectively.From (30) and (31), let us introduce the following extended data matrices:(33)X1=X11,…,XP1T∈CN1subP×K,X2=X12,…,XW2T∈CN2subW×K,or, more compactly,(34)Xr=Φr⋄AsrS+Vr,r=1,2.In contrast with (12), in this approach we shall work with the following sample cross-correlation matrix:(35)R^=Φ1⋄As1RsΦ2⋄As2H+σv2I.As mentioned in Section4.1, we can see that the noiseless term in (35) denotes a multimode unfolding of the following cross-correlation tensor of size N1sub×P×N2sub×W:(36)Ro=Rs×1As1×2Φ1×3As2∗×4Φ2∗,where Rs∈CM×M×M×M. In this modeling, the tensor Ro follows a fourth-order Tucker decomposition. By analogy with (4), we deduce the following correspondences: (37)I1,I2,I3,I4⟷N1sub,P,N2sub,WR1,R2,R3,R4⟷MG,A1,A2,A3,A4⟷Rs,As1,Φ1,As2∗,Φ2∗.The spatial signatures of the sources can be estimated from (36) by using the proposed generalized ALS-Tucker algorithm. Note that, in this case, the ALS-Tucker algorithm is simplified to a fourth-order 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 asei-ei-1≤10-6, where ei denotes the residual error of the ith iteration, defined as(38)ei=R-R^iF2,where R=Ro+V is a noisy version of R, V is an additive complex-valued white Gaussian noise tensor, and R^i is the covariance tensor reconstructed from the estimated factor matrices and core tensor. Since the ALS-ProKRaft algorithm exploits the dual-symmetry 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 matricesA^finalr, r=1,…,R, the final step is to estimate the spatial frequencies of the sources μ^mr, m=1,…,M. The final estimates can be computed from the average over the values obtained in each row of A^finalr as follows:(39)μ^mr=1Nr-1∑n=2NrargA^finalrn,mn-1.


## 5.1. Cross-Correlation Tensor forL-Shaped Sensor Arrays

In this approach, we consider anL-shaped sensor array equipped with N1+N2-1 sensors positioned in the x-z plane, as illustrated in Figure 1. Each linear array contains N1 and N2 sensors equally spaced at distances d1 and d2, respectively. We consider that each linear array is divided into P and W smaller subarrays, respectively. Each subarray contains N1sub=N1-P+1 and N2sub=N2-W+1 sensors. The signal received at the pth subarray, for p=1,…,P, is given by(30)Xp1=As1DpΦ1S+Vp1∈CN1sub×K,and the signal received at the wth subarray, w=1,…,W, is given by(31)Xw2=As2DwΦ2S+Vw2∈CN2sub×K,where (i)
Asr∈CNrsub×M is the spatial signature matrix of the first subarray (or reference subarray) for the rth array dimension, r=1,2;(ii)
DpΦ1 and DwΦ2 are diagonal matrices whose main diagonal is given by the pth and wth row of the matrices Φ(1)∈CP×M and Φ(2)∈CW×M, respectively. The rows of Φ1 and Φ2 capture the delays suffered by the signals impinging the pth and wth subarrays with respect to the reference subarray, which are defined based on the following spatial frequencies:(32)μm1=2π·d1·cos⁡αmλ,μm2=2π·d2·sin⁡αm·cos⁡βmλ,where αm and βm are the azimuth and elevation angles of the mth source, respectively.Figure 1
L-shaped array configuration with N1+N2-1 sensors. The distance between the sensors in the z-axis is d1 while the distance between the sensors in the x-axis is d2. The mth wavefront has elevation and azimuth angles equal to αm and βm, respectively.From (30) and (31), let us introduce the following extended data matrices:(33)X1=X11,…,XP1T∈CN1subP×K,X2=X12,…,XW2T∈CN2subW×K,or, more compactly,(34)Xr=Φr⋄AsrS+Vr,r=1,2.In contrast with (12), in this approach we shall work with the following sample cross-correlation matrix:(35)R^=Φ1⋄As1RsΦ2⋄As2H+σv2I.As mentioned in Section4.1, we can see that the noiseless term in (35) denotes a multimode unfolding of the following cross-correlation tensor of size N1sub×P×N2sub×W:(36)Ro=Rs×1As1×2Φ1×3As2∗×4Φ2∗,where Rs∈CM×M×M×M. In this modeling, the tensor Ro follows a fourth-order Tucker decomposition. By analogy with (4), we deduce the following correspondences: (37)I1,I2,I3,I4⟷N1sub,P,N2sub,WR1,R2,R3,R4⟷MG,A1,A2,A3,A4⟷Rs,As1,Φ1,As2∗,Φ2∗.The spatial signatures of the sources can be estimated from (36) by using the proposed generalized ALS-Tucker algorithm. Note that, in this case, the ALS-Tucker algorithm is simplified to a fourth-order 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 asei-ei-1≤10-6, where ei denotes the residual error of the ith iteration, defined as(38)ei=R-R^iF2,where R=Ro+V is a noisy version of R, V is an additive complex-valued white Gaussian noise tensor, and R^i is the covariance tensor reconstructed from the estimated factor matrices and core tensor. Since the ALS-ProKRaft algorithm exploits the dual-symmetry 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 matricesA^finalr, r=1,…,R, the final step is to estimate the spatial frequencies of the sources μ^mr, m=1,…,M. The final estimates can be computed from the average over the values obtained in each row of A^finalr as follows:(39)μ^mr=1Nr-1∑n=2NrargA^finalrn,mn-1.


## 6. Computational Complexity

In the following, we discuss the computational complexity of the iterative ALS-Tucker and ALS-ProKRaft 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 sizeI1×I2 the number of floating-point operations associated with the SVD computation is OI1·I2·min(I1·I2) [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 ALS-Tucker algorithm equals the complexity of 2R 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 ALS-ProKRaft algorithm equals the complexity of M(R-1) matrix SVDs associated with the application of the multidimensional LS-KRF algorithm in (26) plus the complexity of one additional matrix SVD associated with the update of the unknown unitary rotation factor matrix T 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 inR-D sensor arrays. As previously stated in Section 4.3, the ALS-ProKRaft algorithm works on the assumption that the sample covariance matrix of the sources Rs 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., K→∞), as well as when the source signals are perfectly uncorrelated. In practice, this assumption is not guaranteed. On the other hand, the ALS-Tucker algorithm previously formulated in Section 4.2 naturally captures any structure for the sources covariance into the core tensor Rs. Therefore, the assumption of uncorrelated source signals is not necessary for the ALS-Tucker 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 ALS-Tucker algorithm does not exploit the dual-symmetry property of the data covariance tensor and all factor matrices need to be estimated as independent variables. However, in the ALS-ProKRaft algorithm, only half of the factor matrices are estimated by exploiting the dual-symmetry property of the covariance tensor. Therefore, the ALS-ProKRaft algorithm is more computationally attractive than the ALS-Tucker algorithm. When compared with the state-of-the-art matrix-based algorithms such as MUSIC, ESPRIT, and Propagator Method, the proposed tensor-based 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 ALS-Tucker and ALS-ProKRaft algorithms forR-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 L-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 x-z plane. The mth wavefront has direction of arrival αm,βmm=1M, where αm and βm are elevation and azimuth angles, respectively.In Figures2 and 3, the performance is measured in terms of the root mean square error (RMSE) of the estimated spatial frequencies μ^mr in terms of Signal to Noise Ratio (SNR). The relations between directions of arrival and spatial frequencies are given by (32), where d(r) denotes the distance between sensors in the rth array dimension, which is assumed here equal to λ/2. In Figure 2, we consider Hadamard sequences for the sources, while in Figure 3 we consider BPSK modulated sources. In both cases, we have N=64 sensors (i.e., N1 and N2 equal to 8) and the sample covariance matrix of the received data (12) is calculated from a reduced number K=10 of samples. The total RMSE is defined as(40)RMSE=E∑r=1R∑m=1Mμmr-μ^mr2.Figure 2
Total RMSE versus SNR forN=64 sensors, K=10 samples, and DoAs: 30°,55° and 45°,60° for Hadamard sequences.Figure 3
Total RMSE versus SNR forN=64 sensors, K=10 samples, and DoAs: 30°,55° and 45°,60° for BPSK sequences.From Figure2, it can be seen that both ALS-Tucker and ALS-ProKRaft 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 ALS-ProKraft 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 ALS-Tucker 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 matrix-based standard ESPRIT (SE) [2], matrix-based MUSIC algorithm to planar array configuration [25], and tensor-based standard ESPRIT (STE) [10], the proposed algorithms have improved accuracy in all the considered scenarios.Figure4 shows the convergence performance of the iterative algorithms. In this experiment, the median values of the normalized estimation error ei/N2R are plotted in terms of the number of iterations for different SNR. It is noteworthy that the ALS-ProKRaft algorithm has a faster convergence compared to the ALS-Tucker algorithm. This behavior is expected since ALS-ProKRaft exploits the dual-symmetry property of the data tensor, which results in estimating half as many factor matrices compared to the ALS-Tucker approach.Figure 4
Convergence of the ALS-Tucker and ALS-ProKRaft algorithms.In the second part of this section, we consider aL-shaped configuration array. In Figure 5, we set N=13 sensors (i.e., N1 and N2 equal to 7) and K=500 samples. Each uniform linear array is divided into P=2 and W=2 subarrays, respectively. In this experiment, the performance of the proposed ALS-Tucker algorithm is compared to the state-of-the-art matrix-based methods, namely, Propagator Method (PM) [26], MUSIC [27], and ESPRIT [28], all of them originally formulated for L-shaped arrays. Note that the ALS-Tucker algorithm presents an improved performance over its competitors, with more evidenced gains in the low-to-medium SNR range. For high SNR values, the performance of the MUSIC method comes close to that of our proposal. However, the ALS-Tucker algorithm dispenses any estimation procedure via bidimensional peak search as occurs with MUSIC, being the former more computationally attractive.Figure 5
Total RMSE versus SNR forN=13 sensors, K=500 samples, and DoAs: 30°,45° and 50°,55°.Figure6 shows the performance of the ALS-Tucker by assuming different number of sensors. In this experiment, we consider K=500 samples. We can observe a better performance in terms of RMSE when the number sensors of the L-shaped array is increased. This is valid for all the simulated SNR values.Figure 6
Total RMSE versus SNR (performance of the ALS-Tucker algorithm for different number of sensors).In Figure7, we analyze the influence of the number K of samples on the performance of the ALS-Tucker 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 ALS-Tucker algorithm outperforms the state-of-the-art PM, ESPRIT, and MUSIC methods.Figure 7
Total RMSE versus number of samples.


## 9. Conclusion

In this paper, two tensor-based 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 covariance-based generalization of the Tucker decomposition forR-D sensor arrays. Then, a link between the ALS-ProKRaft algorithm and covariance-based array signal processing for blind spatial signatures estimation has been established. As another contribution, we have formulated a cross-correlation-based fourth-order Tucker decomposition which makes the proposed ALS-Tucker algorithm applicable in scenarios composed by L-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 ALS-Tucker 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 dual-symmetry property of the covariance tensor, which makes the ALS-ProKRaft algorithm preferable since it provides good estimation accuracy with a smaller number of ALS iterations. Finally, when compared with other state-of-the-art matrix-based and tensor-based techniques, the proposed tensor-based iterative algorithms have shown their effectiveness with remarkable gains in terms of estimation error.



---

*Source: 1615962-2017-02-15.xml*