Multi-dimensional model order selection

Multi-dimensional model order selection (MOS) techniques achieve an improved accuracy, reliability, and robustness, since they consider all dimensions jointly during the estimation of parameters. Additionally, from fundamental identifiability results of multi-dimensional decompositions, it is known that the number of main components can be larger when compared to matrix-based decompositions. In this article, we show how to use tensor calculus to extend matrix-based MOS schemes and we also present our proposed multi-dimensional model order selection scheme based on the closed-form PARAFAC algorithm, which is only applicable to multi-dimensional data. In general, as shown by means of simulations, the Probability of correct Detection (PoD) of our proposed multi-dimensional MOS schemes is much better than the PoD of matrix-based schemes.


Introduction
In the literature, matrix array signal processing techniques are extensively used in a variety of applications including radar, mobile communications, sonar, and seismology.To estimate geometrical/physical parameters such as direction of arrival, direction of departure, time of direction of arrival, and Doppler frequency, the first step is to estimate the model order, i.e., the number of signal components.
By taking into account only one dimension, the problem is seen from just one perspective, i.e., one projection.Consequently, parameters cannot be estimated properly for certain scenarios.To handle that, multidimensional array signal processing, which considers several dimensions, is studied.These dimensions can correspond to time, frequency, or polarization, but also spatial dimensions such as one-or two-dimensional arrays at the transmitter and the receiver.With multidimensional array signal processing, it is possible to estimate parameters using all the dimensions jointly, even if they are not resolvable for each dimension separately.Moreover, by considering all dimensions jointly, the accuracy, reliability, and robustness can be improved.Another important advantage of using multi-dimensional data, also known as tensors, is the identifiability, since with tensors the typical rank can be much higher than using matrices.Here, we focus particularly on the development of techniques for the estimation of the model order.
The estimation of the model order, also known as the number of principal components, has been investigated in several science fields, and usually model order selection schemes are proposed only for specific scenarios in the literature.Therefore, as a first important contribution, we have proposed in [1,2] the one-dimensional model order selection scheme called Modified Exponential Fitting Test (M-EFT), which outperforms all the other schemes for scenarios involving white Gaussian noise.Additionally, we have proposed in [1,2] improved versions of the Akaike's Information Criterion (AIC) and Minimum Description Length (MDL).
As reviewed in this article, the multi-dimensional structure of the data can be taken into account to improve further the estimation of the model order.As an example of such improvement, we show our proposed R-dimensional Exponential Fitting Test (R-D EFT) for multi-dimensional applications, where the noise is additive white Gaussian.The R-D EFT successfully outperforms the M-EFT confirming that even the technique with the best performance can be improved by taking into account the multi-dimensional structure of the data [1,3,4].In addition, we also extend our modified versions of AIC and MDL to their respective multi-dimensional versions R-D AIC and R-D MDL.For scenarios with colored noise, we present our proposed multi-dimensional model order selection technique called closed-form PARAFAC-based model order selection (CFP-MOS) scheme [3,5].
The remainder of this article is organized as follows.After reviewing the notation in second section, the data model is presented in third section.Then the R-dimensional exponential fitting test (R-D EFT) and closedform PARAFAC-based model order selection (CFP-MOS) scheme are reviewed in fourth section.The simulation results in fifth section confirm the improved performance of R-D EFT and CFP-MOS.Conclusions are drawn finally.

Tensor and matrix notation
In order to facilitate the distinction between scalars, matrices, and tensors, the following notation is used: Scalars are denoted as italic letters (a, b, ..., A, B, ..., a, b, ...), column vectors as lower-case bold-face letters (a, b, ...), matrices as bold-face capitals (A, B, ...), and tensors are written as bold-face calligraphic letters (A, B, . ..).Lower-order parts are consistently named: the (i, j)-element of the matrix A is denoted as a i,j and the (i, j, k)-element of a third order tensor X as x i,j,k .The n-mode vectors of a tensor are obtained by varying the nth index within its range (1, 2, ..., I n ) and keeping all the other indices fixed.We use the superscripts T, H, -1, +, and * for transposition, Hermitian transposition, matrix inversion, the Moore-Penrose pseudo inverse of matrices, and complex conjugation, respectively.Moreover the Khatri-Rao product (columnwise Kronecker product) is denoted by A ◊ B.
The tensor operations we use are consistent with [6]: The r-mode product of a tensor A ∈ C I 1 ×I 2 ×•••×I R and a matrix U ∈ C J r ×I ralong the rth mode is denoted as It is obtained by multiplying all r-mode vectors of A from the left-hand side by the matrix U. A certain r-mode vector of a tensor is obtained by fixing the rth index and by varying all the other indices.
The higher-order SVD (HOSVD) of a tensor where S ∈ C I 1 ×I 2 ×•••×I R is the core-tensor which satis- fies the all-orthogonality conditions [6] and U r ∈ C I r ×I r , r = 1, 2, ..., R are the unitary matrices of r-mode singular vectors.
Finally, the r-mode unfolding of a tensor A is symbo- lized by [A] (r) ∈ C I r ×(I 1 I 2 ...I r−1 I r+1 ...I R ) , i.e., it represents the matrix of r-mode vectors of the tensor A .The order of the columns is chosen in accordance with [6].

Data model
To validate the general applicability of our proposed schemes, we adopt the PARAFAC data model below where f (r) n (m r ) is the m r th element of the nth factor of the rth mode for m r = 1, ..., M r and r = 1, 2, ..., R, R +1.The M R+1 can be alternatively represented by N, which stands for the number of snapshots.
By defining the vectors T and using the outer product operator ∘, another possible representation of ( 2) is given by where is composed of the sum of d rank one tensors.Therefore, the tensor rank of X 0 coincides with the model order d.
For applications, where the multi-dimensional data obeys a PARAFAC decomposition, it is important to estimate the factors of the tensor X 0 , which are defined as ×d , and we assume that the rank of each F (r) is equal to min(M r , d).This definition of the factor matrices allows us to rewrite (3) according to the notation proposed in [7] where × r is the r-mode product defined in Section 2, and the tensor I R+1,d represents the R-dimensional iden- tity tensor of size d × d... × d, whose elements are equal to one when the indices i 1 = i 2 ... = i R+1 and zero otherwise.
In practice, the data is contaminated by noise, which we represent by the following data model where is the additive noise tensor, whose elements are i.i.d.zero-mean circularly symmetric complex Gaussian (ZMCSCG) random variables.Thereby, the tensor rank is different from d and usually it assumes extremely large values as shown in [8].Hence, the problem we are solving can therefore be stated in the following fashion: given a noisy measurement tensor X , we desire to estimate the model order d.Note that according to Comon [8], the typical rank of X is much bigger than any of the dimensions M r for r = 1, ..., R + 1.
The objective of the PARAFAC decomposition is to compute the estimated factors F(r) such that Since F(r) ∈ C M r ×d one requirement to apply the PAR- AFAC decomposition is to estimate d.
We evaluate the performance of the model order selection scheme in the presence of colored noise, which is given by replacing the white Gaussian white noise tensor N by the colored Gaussian noise tensor N (c) in (5).Note that the data model used in this article is simply a linear superposition of rank-one components superimposed by additive noise.
Particularly, for multi-dimensional data, the colored noise with a Kronecker structure is present in several applications.For example, in EEG applications [9], the noise is correlated in both space and time dimensions, and it has been shown that a model of the noise combining these two correlation matrices using the Kronecker product can fit noise measurements.Moreover, for MIMO systems the noise covariance matrix is often assumed to be the Kronecker product of the temporal and spatial correlation matrices [10].
The multi-dimensional colored noise, which is assumed to have a Kronecker correlation structure, can be written as where ⊗ represents the Kronecker product.We can also rewrite (7) using the n-mode products in the following fashion where is a tensor with uncor- related ZMCSCG elements with variance σ 2 n , and L i ∈ C M i ×M i is the correlation factor of the ith dimension of the colored noise tensor.The noise covariance matrix in the ith mode is defined as where a is a normalization constant, such that The equivalence between (7), (8), and ( 9) is shown in [11].
To simplify the notation, let us define M = R r=1 M r .For the r-mode unfolding we compute the sample covariance matrix as The eigenvalues of these r-mode sample covariance matrices play a major role in the model order estimation step.Let us denote the ith eigenvalue of the sample covariance matrix of the r-mode unfolding as λ (r) i .Notice that R (r)  xx possesses M r eigenvalues, which we order in such a way that λ M r .The eigenvalues may be computed from the HOSVD of the measurement tensor Note that the eigenvalues λ (r) i are related to the . The r-mode singular values σ (r) i can also be computed via the SVD of the r-mode unfolding X as follows where

Multi-dimensional model order selection schemes
In this section, the multi-dimensional model order selection schemes are proposed based on the global eigenvalues, the R-D subspace, or tensor-based data model.First, we show the proposed definition of the global eigenvalues together with the presentation of the proposed R-D EFT.Then, we summarize our multidimensional extension of AIC and MDL.Besides the global eigenvalues-based schemes, we also propose a tensor data-based multi-dimensional model order selection scheme.Followed by the closed-form PARAFACbased model order selection scheme is proposed for white and also colored noise scenarios.For data sampled on a grid and an array with centro-symmetric symmetries, we show how to improve the performance of model order selection schemes for such data by incorporating forward-backward averaging (FBA).

R-D exponential fitting test (R-D EFT)
The global eigenvalues are based on the r-mode eigenvalues represented by λ (r) i for r = 1, ..., R and for i = 1, ..., M r .To obtain the r-mode eigenvalues, there are two ways.The first way shown in (10) is possible via the EVD of each r-mode sample covariance matrix, and the second way in ( 12) is given via an HOSVD.
According to Grouffaud et al. [12] and Quinlan et al. [13], the noise eigenvalues that exhibit a Wishart profile can have their profile approximated by an exponential curve.Therefore, by applying the exponential approximation for every r-mode, we obtain that where .., M r and r = 1, 2, ..., R + 1.The rate of the exponential profile q(a r , b r ) is defined as where a = min (M, N) and b = max (M, N).Note that ( 15) of the M-EFT is an extension of the EFT expression in [12,13].
In order to be even more precise in the computation of q, the following polynomial can be solved Although from ( 16) a + 1 solutions are possible, we select only the q that belongs to the interval (0, 1).For M ≤ N (15) is equal to the q of the EFT [12,13], which means that the PoD of the EFT and the PoD of the M-EFT are the same for M <N.Consequently, the M-EFT automatically inherits from the EFT the property that it outperforms the other matrix-based MOS techniques in the literature for M ≤ N in the presence of white Gaussian noise as shown in [2].
For the sake of simplicity, let us first assume that M 1 = M 2 = ... = M R .Then we can define global eigenvalues as being [1] Therefore, based on (14), it is straightforward that the noise global eigenvalues also follow an exponential profile, since where i = 1, ..., M R+1 .
In Figure 1, we show an example of the exponential profile property that is assumed for the noise eigenvalues.This exponential profile approximates the distribution of the noise eigenvalues and the distribution of the global noise eigenvalues.The exemplified data in Figure 1 have the model order equal to one, since the first eigenvalue does not fit the exponential profile.To estimate the model order, the noise eigenvalue profile gets predicted based on the exponential profile assumption starting from the smallest noise eigenvalue.When a significant gap is detected compared to this predicted exponential profile, the model order, i.e., the smallest signal eigenvalue, is found.
The product across modes increases the gap between the predicted and the actual eigenvalues as shown in Figure 1.We compare the gap between the actual eigenvalues and the predicted eigenvalues in the rth mode to the gap between the actual global eigenvalues and the predicted global eigenvalues.Here, we consider that X 0 is a rank one tensor, and noise is added according to (5) Then, in this case, d = 1.For the first gap, we have λ , while for the second one, we have λ . Therefore, the break in the profile is easier to detect via global eigenvalues than using only one mode eigenvalues Since all tensor dimensions may be not necessarily equal to each other, without loss of generality, let us consider the case in which M 1 ≥ M 2 ≥ ... ≥ M R+1 .In Figures 2, 3, and 4, we have sets of eigenvalues obtained from each r-mode of a tensor with sizes M 1 = 13, M 2 = 11, M 3 = 8 and M 4 = 3.The index i indicates the position of the eigenvalues in each rth eigenvalues set.
We start by estimating d with a certain eigenvaluebased model order selection method considering the first unfolding only, which in the example in Figure 2 has a size M 1 = 13.If d < M 2 , we could have taken advantage of the second mode as well.Therefore, we compute the global eigenvalues λ  the global eigenvalues stops using the three first modes.Clearly, the full potential of the proposed method can be achieved when all modes are used to compute the global eigenvalues.This happens when d < M R+1 , so that λ (G) i can be computed for 1 ≤ i ≤ M R+1 .Note that using the global eigenvalues, the assumptions of M-EFT, that the noise eigenvalues can be approximated by an exponential profile, and the assumptions of AIC and MDL, that the noise eigenvalues are constant, still hold.Moreover, the maximum model order is equal to max r M r, for r = 1, ..., R.
The R-D EFT is an extended version of the M-EFT operating on the λ (G) i .Therefore, 1) It exploits the fact that the noise global eigenvalues still exhibit an exponential profile; 2) The increase of the threshold between the actual signal global eigenvalue and the predicted noise global eigenvalue leads to a significant improvements in the performance; 3) It is applicable to arrays of arbitrary size and dimension through the sequential definition of the global eigenvalues as long as the data is arranged on a multi-dimensional grid.
To derive the proposed multi-dimensional extension of the M-EFT algorithm, namely the R-D EFT, we start by looking at an R-dimensional noise-only case.For the R-D EFT, it is our intention to predict the noise global eigenvalues defined in (18).Each r-mode eigenvalue can be estimated via λ(r) M−P = (P + 1) Equations ( 19) and ( 20) are the same expressions as in the case of the M-EFT in [2], however, in contrast to the M-EFT, here they are applied to each r-mode eigenvalue.
Let us apply the definition of the global eigenvalues according to (17) where in (18) the approximation by an exponential profile is assumed.Therefore, where a (G) is the minimum a r for all the r-modes considered in the sequential definition of the global eigenvalue.In (22), λ(G) i is a function of only the last global eigenvalue λ(G) α (G) , which is the smallest global eigenvalue and is assumed a noise eigenvalue, and of the rates q P + 1, M M r for all the r-modes considered in the sequential definition.Instead of using directly (22), we use λ(r) M−P according to (19) for all the r-modes considered in the sequential definition.Therefore, the previous eigenvalues that were already estimated as noise eigenvalues are taken into account in the prediction step.
Similarly to the M-EFT, using the predicted global eigenvalue expression (21) considering white Gaussian noise samples, we compute the global threshold coefficients η (G) P via the hypotheses for the tensor case Once all η (G) P are found for a certain higher order array of sizes M 1 , M 2 , ..., M R , and for a certain P fa , then   the model order can be estimated by applying the following cost function d = α (G) − min(P) where where a (G) is the total number of sequentially defined global eigenvalues.

R-D AIC and R-D MDL
In AIC and MDL, it is assumed that the noise eigenvalues are all equal.Therefore, once this assumption is valid for all r-mode eigenvalues, it is straightforward that it is also valid for our global eigenvalue definition.Moreover, since we have shown in [2] that 1-D AIC and 1-D MDL are more general and superior in terms of performance than AIC and MDL, respectively, we extend 1-D AIC and 1-D MDL to the multi-dimensional form using the global eigenvalues.Note that the PoD of 1-D AIC and 1-D MDL is only greater than the PoD of AIC and MDL for cases where M M r > M r, which cannot be fulfilled for one-dimensional data.
The corresponding R-dimensional versions of 1-D AIC and 1-D MDL are obtained by first replacing the eigenvalues Rxx by the global eigenvalues λ (G) i defined in (17).Additionally, to compute the number of free parameters for the 1-D AIC and 1-D MDL methods and their R-D extensions, we propose to set the parameter N = max r M r and a (G) is the total number of sequentially defined global eigenvalues similarly as we propose in [1].Therefore, the optimization problem for the R-D AIC and R-D MDL is given by d = arg min P J (G) (P) where where d represents an estimate of the model order d, and g (G) (P) and a (G) (P) are the geometric and arithmetic means of the P smallest global eigenvalues, respectively.The penalty functions p(P, N a (G) ) for R-D AIC and R-D MDL are given in Table 1.
Note that the R-dimensional extension described in this section can be applied to any model order selection scheme that is based on the profile of eigenvalues, i.e., also to the 1-D MDL and the 1-D AIC methods.

Closed-form PARAFAC-based model order selection (CFP-MOS) scheme
In this section, we present the Closed-form PARAFACbased model order selection (CFP-MOS) technique proposed in [5].The major motivation of CFP-MOS is the fact that R-D AIC, R-D MDL, and R-D EFT are applicable only in the presence of white Gaussian noise.Therefore, it is very appealing to apply CFP-MOS, since it has a performance close to R-D EFT in the presence of white Gaussian noise, and at the same time it is also applicable in the presence of colored Gaussian noise.
According to Roemer and Haardt [14], the estimation of the factors F (r) via the PARAFAC decomposition is transformed into a set of simultaneous diagonalization problems based on the relation between the truncated HOSVD [6]-based low-rank approximation of X X ≈ S [s] and the PARAFAC decomposition of X where r ∈ C M r ×p r, p r = min (M r , d), and F(r) = U The closed-form PARAFAC (CFP) [14] decomposition constructs two simultaneous diagonalization problems for every tuple (k,ℓ), such that k, ∈ R, and k < ℓ.In order to reference each simultaneous matrix diagonalization (SMD) problem, we define the enumerator function e(k, ℓ, i) that assigns the triple (k, ℓ, i) to a sequence of consecutive integer numbers in the range 1, 2, ..., T.Here i = 1, 2 refers to the two simultaneous matrix diagonalizations (SMD) for our specific k and ℓ.Consequently, SMD (e (k, ℓ, 1), P) represents the first SMD for a given k and ℓ, which is associated to the simultaneous diagonalization of the matrices S rhs k, ,(n) by T k .Initially, we consider that the candidate value of the model order P = d, which is the model order.Similarly, SMD (e (k, ℓ, 2), P) corresponds to the second SMD for a given k and ℓ referring to the simultaneous diagonalizations of S lhs k, ,(n) by T ℓ .S rhs k, ,(n) and S lhs k, ,(n) are defined in [14].Note that each SMD(e(k, ℓ, i), P) yields an estimate of all factors F (r) [14,15], where r = 1, ..., R. Consequently, for each factor F (r) there are T estimates.
For instance, consider a 4-D tensor, where the third mode is degenerate, i.e., M 3 <d.Then, the set R + 1 is given by {1, 2, 4}, and the possible (k, ℓ)-tuples are (1,2), (1,4), and (2,4).Consequently, the six possible SMDs are enumerated via e(k, ℓ, i) as follows: e(1, 2, 1) = 1, e(1, 2, 2) = 2, e(1, 4, 1) = 3, e(1, 4, 2) = 4, e(2, 4, 1) = 5, and e (2, 4, 2) = 6.In general, the total number of SMD problems T is equal to There are different heuristics to select the best estimates of each factor F (r) as shown in [14].We define the function to compute the residuals (RESID) of the simultaneous matrix diagonalizations (SMD) as RESID (SMD(•)).For instance, we apply it to e(k, ℓ, 1) and for e(k, ℓ,2) where Since each residual is a positive real-valued number, we can order the SMDs by the magnitude of the corresponding residual.For the sake of simplicity, we represent the ordered sequence of SMDs to e(k, ℓ, i) by a single e (t) for t = 1, 2, ..., T, such that RESID(SMD (e (t) , P)) ≤ RESID(SMD(e (t+1) , P)).Since in practice d is not known, P denotes a candidate value for d, which is our estimate of the model order d.Our task is to select P from the interval dmin ≤ P ≤ dmax , where dmin is a lower bound and dmax is an upper bound for our candidate values.For instance, dmin equal to 1 is used, and dmax is chosen such that no dimension is degenerate [14], i.e. d ≤ M r for r = 1, ..., R. We define RESID(SMD (e (t) , P)) as being the tth lowest residual of the SMD considering the number of components per factor equal to P. Based on the definition of RESID(SMD(e (t) ,P)), one first direct way to estimate the model order d can be performed using the following properties 1) If there is no noise and P <d, then RESID(SMD(e (t) , P)) > RESID(SMD(e (t) , d)), since the matrices generated are composed of mixed components as shown in [16].
2) If noise is present and P >d, then RESID(SMD(e (t) , P)) > RESID(SMD(e (t) , d)), since the matrices generated with the noise components are not diagonalizable commuting matrices.Therefore, the simultaneous diagonalizations are not valid anymore.
Based on these properties, a first model order selection scheme can be proposed d = arg min P RESID(SMD(e (1) , P)). (29) However, the model order selection scheme in (29) yields a Probability of correct Detection (PoD) inferior to the some MOS techniques found in the literature.Therefore, to improve the PoD of (29), we propose to exploit the redundant information provided only by the closed-form PARAFAC (CFP) [14].
Let F(r) e (t) ,P denote the ordered sequence of estimates for F (r) assuming that the model order is P.In order to combine factors estimated in different diagonalizations processes, the permutation and scaling ambiguities should be solved.For this task, we apply the amplitude approach according to Weis et al. [15].For the correct model order and in the absence of noise, the subspaces of F (r) e (t) ,P should not depend on t.Consequently, a measure for the reliability of the estimate is given by comparing the angle between the vectors f (r) v,e (t) ,P for different t, where f (r) v,e (t) ,P corresponds to the estimate of the vth column of F v,e (t) ,P , f (r) v,e (1) ,P , where the operator ∢ gives the angle between two vectors and T lim represents the total number of simultaneous matrix diagonalizations taken into account.T lim , a design parameter of the CFP-MOS algorithm, can be chosen between 2 and T. Similar to the Threshold Core Consistency Analysis (T-CORCONDIA) in [4], the CFP-MOS requires weights Δ(P), otherwise the Probabilities of correct Dectection (PoD) for different values of d have a significant gap from each other.Therefore, to have a fair estimation for all candidates P, we introduce the weights Δ(P), which are calibrated in a scenario with white Gaussian noise, where the number of sources d varies.For the calibration of weights, we use the probability of correct detection (PoD) of the R-D EFT [1,4]  where E{PoD R -D EFT SNR (P)} returns the averaged probability of correct detection over a certain predefined SNR range using the R-D EFT for a given scenario assuming P as the model order, d max is defined as being the maximum candidate value of P, and Δ var is the vector with the threshold coefficients for each value of P. Note that the elements of the vector of weights Δ vary according to a certain defined range and interval and that the averaged PoD of the CFP-MOS is compared to the averaged PoD of the R-D EFT.When the cost function is minimized, then we have the desired Δ var .
Up to this point, the CFP-MOS is applicable to scenarios without any specific structure in the factor matrices.If the vectors f (r) v,e (t) ,P have a Vandermonde structure, we can propose another expression.Again let F(r) e (t) ,P be the estimate for the factor matrix obtained from SMD(e (t) , P).Using the Vandermonde structure of each factor we can estimate the scalars μ (r) v,e (t) ,P corresponding to the vth column of F(r) e (t) ,P As already proposed previously, for the correct model order and in the absence of noise, the estimated spatial frequencies should not depend on t.Consequently, a measure for the reliability of the estimate is given by comparing the estimates for different t.
Similar to the cost function in (30), to have a fair estimation for all candidates P, we introduce the weights Δ(P), which are calculated in a similar fashion as for T-CORCONDIA Var in [4] by considering data contaminated by white Gaussian noise.

Applying forward-backward averaging (FBA)
In many applications, the complex-valued data obeys additional symmetry relations that can be exploited to enhance resolution and accuracy.For instance, when sampling data uniformly or on centro-symmetric grids, the corresponding r-mode subspaces are invariant under flipping and conjugation.Such scenarios are known as having centro-symmetric symmetries.Also in such scenarios, we can incorporate FBA [17] to all model order selection schemes even with a multidimensional data model.First, let us present modifications in the data model, which should be considered to apply the FBA.Comparing the data model of ( 4) to the data model to be introduced in this section, we summarize two main differences.The first one is the size of X 0 , which has R + 1 dimensions instead of the R dimensions as in (4).Therefore, the noiseless data tensor is given by This additional (R + 1)th dimension is due to the fact that the (R + 1)th factor represents the source symbols matrix F (R+1) = S T .The second difference is the restriction of the factor matrices F (r) = for r = 1, ..., R of the tensor X 0 in (33) to a matrix, where each vector is a function of a certain scalar μ (r) i related to the rth dimension and the ith source.In many applications, these vectors have a Vandermonde structure.For the sake of notation, the factor matrices for r = 1, ..., R are represented by A (r) , and it can be written as a function of μ (r) i as follows In [18,19] it was demonstrated that in the tensor case, forward-backward averaging can be expressed in the following form where [A n B] represents the concatenation of two tensors A and B along the nth mode.Note that all the other modes of A and B should have exactly the same sizes.The matrix Π n is defined as In multi-dimensional model order selection schemes, forward-backward averaging is incorporated by replacing the data tensor X in (11) by Z .Moreover, we have to replace N by 2 • N in the subsequent formulas since the number of snapshots is virtually doubled.
In schemes like AIC, MDL, 1-D AIC, and 1-D MDL, which requires the information about the number of sensors and the number of snapshots for the computation of the free parameters, once FBA is applied, the number of snapshots in the free parameters should be updated from N to 2 • N.
To reduce the computational complexity, the forwardbackward averaged data matrix Z can be replaced by a real-valued data matrix {Z} ℝ M × 2N which has the same singular values as Z [20].This transformation can be extended to the tensor case where the forward-backward averaged data tensor Z is replaced by a real-valued data tensor ϕ{Z} ∈ R M 1 ×•••×M R ×2N possessing the same r-mode singular values for all r = 1, 2, ..., R + 1 (see [19] for details).
where Z is given in (35), and if p is odd, then Q p is given as and p = 2 • n + 1.On the other hand, if p is even, then Q p is given as and p = 2 • n.

Simulation results
In this section, we evaluate the performance, in terms of the probability of correct detection (PoD), of all multidimensional model order selection techniques presented previously Monte Carlo simulations considering different scenarios.
Comparing the two versions of the CORCONDIA [4,21] and the HOSVD-based approaches, we can notice that the computational complexity is much lower in the R-D methods.Moreover, the HOSVD-based approaches outperform the iterative approaches, since none of them are close to the 100% Probability of correct Detection (PoD).The techniques based on global eigenvalues, R-D EFT, R-D AIC, and R-D MDL maintain a good performance even for lower SNR scenarios, and the R-D EFT shows the best performance if we compare all the techniques.
In Figures 5 and 6, we observe the performance of the classical methods and the R-D EFT, R-D AIC, and R-D MDL for a scenario with the following dimensions M 1 = 7, M 2 = 7, M 3 = 7, and M 4 = 7.The methods described as M-EFT, AIC, and MDL correspond to the simplified one-dimensional cases of the R-D methods, in which we consider only one unfolding for r = 4.
In Figures 7 and 8, we compare our proposed approach to all mentioned techniques for the case that white noise is present.To compare the performance of CFP-MOS for various values of the design parameter T lim , we select T lim = 2 for the legend CFP 2f and T lim = 4 for CFP 4f.In Figure 7, the model order d is equal to 2, while in Figure 8, d = 3.In these two scenarios, the proposed CFP-MOS has a performance very close to R-D EFT, which has the best performance.
In Figures 9 and 10, we assume the noise correlation structure of Equation ( 9), where W i of the ith factor for M i = 3 is given by where r i is the correlation coefficient.Note that also other types of correlation models different from (40) can be used.
In Figures 9 and 10, the noise is colored with a very high correlation, and the factors L i are computed based on ( 9) and (40) as a function of r i .As expected for this scenario, the R-D EFT, R-D AIC, and R-D MDL completely fail.In case of colored noise with high correlation, the noise power is much more concentrated in the signal components.Therefore, the smaller are the values of d, the worse is the PoD.The behavior of the CFP-MOS, AIC, MDL, and EFT are consistent with this effect.The PoD of AIC, MDL, and EFT increases from 0.85, 0.7, and 0.7 in Figure 9 to 0.9, 0.85, and 0.85 in Figure 10.CFP-MOS 4f has a PoD = 0.98 for SNR = 20 dB in Figure 9, while a PoD = 0.98 for SNR = 15 dB in Figure 10.
In contrast to CFP-MOS, AIC, MDL, and EFT, the PoD of RADOI [22] degrades from Figures 9 and 10.In Figure 9, RADOI has a better performance than the CFP-MOS version, while in Figure 10, CFP-MOS outperforms RADOI.Note that the PoD for RADOI becomes constant for SNR ≤ 3 dB, which corresponds to a biased estimation.Therefore, for severely colored noise scenarios, the model order selection using CFP-MOS is more stable than the other approaches.
In Figure 11, no FBA is applied in all model order selection techniques, while in Figure 12 FBA is applied in all of them according section 4. In general, an improvement of approximately 3 dB is obtained when FBA is applied.
In Figure 12, d = 3.Therefore, using the sequential definition of the global eigenvalues from "R-D Exponential Fitting Test (R-D EFT)", we can estimate the model order considering four modes.By increasing the number of sources to 5 in Figure 13, the sequential definition of the global eigenvalues is computed considering the second, third, and fourth modes, which are related to M 2 , M 3 , and N.
By increasing the number of sources even more such that only one mode can be applied, the curves of the R-D EFT, R-D AIC and R-D MDL are the same as the curves of M-EFT, 1-D AIC, and 1-D MDL, as shown in Figure 14.

Conclusions
In this article, we have compared different model order selection techniques for multi-dimensional high-resolution parameter estimation schemes.We have achieved the following results considering a multi-dimensional data model.
1) In case of white Gaussian noise scenarios, our R-D EFT outperforms the other techniques presented in the literature.
2) In the presence of colored noise, the CFP-MOS is the best technique, since it has a performance close to the R-D EFT in case of no correlation, and a performance more stable than RADOI, in case of severely correlated noise.
3) For researchers, which prefer to use information theoretic criteria (ITC) techniques, we have also proposed multi-dimensional extensions of AIC and MDL, called R-D AIC and R-D MDL, respectively.
In Table 2, we summarize the scenarios to apply the different techniques shown in this article.Also in Table 2, wht stands for white noise and clr stands for colored noise.Note that the PoD of the CFP-MOS is close to the one of the R-D EFT for white noise, which means that it has a multi-dimensional gain.Moreover, since the CFP-MOS is suitable for white and colored noise applications, we consider it the best general-purpose scheme.
as in(17) for 1 ≤ i ≤ M 2 , thus discarding the M 1 -M 2 last eigenvalues of the first mode.We can obtain a new estimate d.As illustrated in Figure3, we utilize only the first M 2 highest eigenvalues of the first and of the second modes to estimate the model order.If d < M 3 we could continue in the same fashion, by computing the global eigenvalues considering the first three modes.In the example in Figure4, since the model order is equal to 6, which is greater than M 4 , the sequential definition algorithm of

Figure 2
Figure 2 Sequential definition of the global eigenvalues-1st eigenvalue set.

Figure 3
Figure 3 Sequential definition of the global eigenvalues-1st and 2nd eigenvalue sets.

r
• T r for a nonsingular transformation matrix T r ℂ d × d for all modes r ∈ R where R = {r|M r ≥ d, r = 1, . . .R + 1} denotes the set of non-degenerate modes.As shown in (25) and in (26), the operator R+1 × r=1 r denotes a compact representation of R r-mode products between a tensor and R + 1 matrices.

e
(t) ,P .Hence, this gives rise to an expression to estimate the model order using CFP-MOS d = arg min P as a reference, since the R-D EFT achieves the best PoD in the literature even in the low SNR regime.Consequently, we propose the following expression to obtain the calibrated weights Δ var var = arg min J var ( ) where J var ( ) = dmax P=dmin E PoD CFP -MOS SNR ( (P)) − E{PoD R -D EFT SNR (P)} (31) da Costa et al.EURASIP Journal on Advances in Signal Processing 2011, 2011:26 http://asp.eurasipjournals.com/content/2011/1/26 Hence, this gives rise to the new cost function d = arg min P v,e (t) ,P − μ(r) v,e (1) ,P .

Figure 11
Figure 11 Probability of correct Detection (PoD) versus SNR for an array of size M 1 = 5, M 2 = 7, and M 3 = 9.The number of snapshots N is set to 10 and the number of sources d = 3.No FBA is applied.

Figure 12 Figure 13 Figure 14
Figure 12 Probability of correct Detection (PoD) versus SNR for an array of size M 1 = 5, M 2 = 7, and M 3 = 9.The number of snapshots N is set to 10 and the number of sources d = 3. FBA is applied.

Table 1
Penalty functions for R-D information theoretic criteria ApproachPenalty function p(P, N, a