Citation/Reference M. Sorensen and L. De Lathauwer (2017),
Multidimensional Harmonic Retrieval via Coupled Canonical Polyadic Decomposition --- Part II: Algorithm and Multirate Sampling IEEE Transactions on Signal Processing, vol. 65, no. 2, pp. 528-539, Jan.15, 15 2017.
Archived version Author manuscript: the content is identical to the content of the published paper, but without the final typesetting by the publisher
Published version http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=7581082&isnumbe r=7748608
Journal homepage http://ieeexplore.ieee.org/
Author contact Mikael.Sorensen@kuleuven.be
Abstract
IR https://lirias.kuleuven.be/handle/123456789/550210
(article begins on next page)
Multidimensional Harmonic Retrieval via Coupled Canonical Polyadic Decomposition
— Part II: Algorithm and Multirate Sampling
Mikael Sørensen and Lieven De Lathauwer, Fellow, IEEE
Abstract—In Part I, we have presented a link between Multidimensional Harmonic Retrieval (MHR) and the re- cently proposed coupled Canonical Polyadic Decomposition (CPD), which implies new uniqueness conditions for MHR that are more relaxed than the existing results based on a Vandermonde constrained CPD. In Part II, we explain that the coupled CPD also provides a computational framework for MHR. In particular, we present an algebraic method for MHR based on simultaneous matrix diagonalization that is guaranteed to find the exact solution in the noiseless case, under conditions discussed in Part I. Since the simultaneous matrix diagonalization method reduces the MHR problem into an eigenvalue problem, the proposed algorithm can be interpreted as a MHR generalization of the classical ESPRIT method for one-dimensional harmonic retrieval. We also demonstrate that the presented coupled CPD framework for MHR can algebraically support multirate sampling. We develop an efficient implementation which has about the same computational complexity for single-rate and multirate sampling. Numerical experiments demonstrate that by simul- taneously exploiting the harmonic structure in all dimensions and making use of multirate sampling, the coupled CPD framework for MHR can lead to an improved performance compared to the conventional Vandermonde constrained CPD approaches.
Index Terms—coupled canonical polyadic decomposition, simultaneous matrix diagonalization, multidimensional har- monic retrieval, multirate sampling, direction of arrival estimation.
I. Introduction
In the signal processing literature, the uniqueness properties of Multidimensional Harmonic Retrieval (MHR) have mainly been studied using arguments that
M. Sørensen and L. De Lathauwer are with KU Leuven - E.E. Dept.
(ESAT) - STADIUS Center for Dynamical Systems, Signal Processing and Data Analytics, Kasteelpark Arenberg 10, B-3001 Leuven-Heverlee, Belgium, the Group Science, Engineering and Technology, KU Leuven Kulak, E. Sabbelaan 53, 8500 Kortrijk, Belgium, and iMinds Medi- cal IT, Kasteelpark Arenberg 10, B-3001 Leuven-Heverlee, Belgium, {Mikael.Sorensen, Lieven.DeLathauwer}@kuleuven.be.
Research supported by: (1) Research Council KU Leuven:
CoE EF/05/006 Optimization in Engineering (OPTEC), C1 project C16/15/059-nD, (2) F.W.O.: project G.0830.14N, G.0881.14N, (3) the Belgian Federal Science Policy Office: IUAP P7 (DYSCO II, Dynamical systems, control and optimization, 2012–2017), (4) EU: The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013) / ERC Advanced Grant: BIOTENSORS (no.
339804). This paper reflects only the authors’ views and the Union is not liable for any use that may be made of the contained information.
implicitly or explicitly relate to properties of a Vander- monde constrained Canonical Polyadic Decomposition (VDM-CPD), see [21] and references therein. Based on this connection several algorithms have been developed (e.g. [4], [7], [9], [10], [5]) that solve the MHR problem via a Generalized EigenValue Decomposition (GEVD).
In the companion paper [23] we presented a link be- tween MHR and the recently introduced coupled Canon- ical Polyadic Decomposition (CPD) modeling framework [22], [25], which led to improved identifiability results. In this paper we explain that the coupled CPD model also provides an algebraic framework from which algorithms for MHR can be derived. We develop an algorithm based on coupled CPD that, similar to ESPRIT [16] for one-dimensional (1D) Harmonic Retrieval (HR), reduces the MHR problem into a set of decoupled single-tone MHR problems via a GEVD. A closed-form solution to the single-tone MHR problem will be provided. In the inexact case, the algorithm can be used to initialize an optimization-based method for MHR (e.g. [2], [20]).
We argue that an additional reason to consider cou- pled CPD as opposed to the conventional VDM-CPD is that the former approach o↵ers a more flexible frame- work. In particular, we show that coupled CPD can support multirate sampling. This makes it interesting for large-scale MHR problems where sub-Nyquist sampling is often a desired feature in order to reduce hardware, data acquisition or computational costs.
The paper is organized as follows. The rest of the introduction will present the notation and the tensor decompositions used throughout the paper. Section II briefly reviews the connections between MHR, VDM- CPD and coupled CPD. From the coupled CPD inter- pretation of MHR, an algebraic ESPRIT-type algorithm is derived in Section III. In Section IV we extend the algorithm to multirate sampling. Numerical experiments are reported in section V. Section VI concludes the paper.
Notation: Vectors, matrices and tensors are denoted by lower case boldface, upper case boldface and upper case calligraphic letters, respectively. The rth column vector of A is denoted by ar. The symbols b and d denote the Kronecker and Khatri-Rao product:
AbB :“
»
—–
a11B a12B . . . a21B a22B . . . ... ... . ..
fi
ffifl , AdB :“ ra1bb1a2bb2. . .s,
in which pAqmn “ amn. The outer product of N vectors apnqP CInis denoted by ap1qbbb¨ ¨ ¨bbbapNqP CI1ˆI2ˆ¨¨¨ˆIN, such that pap1qbbb¨ ¨ ¨bbbapNqqi1,i2,...,iN “ ap1qi1 ap2qi2 ¨ ¨ ¨ apNqiN .
The conjugate, transpose, conjugate-transpose, Frobe- nius norm, Moore-Penrose pseudo-inverse, range and kernel of the matrix A are denoted by A˚, AT, AH, }A}F, A:, rangepAq and kerpAq respectively. The real and imaginary parts of a number a are denoted by Re tau and Imtau, respectively. From the context it should be clear when i denotes the imaginary unit number, i.e., i “ ?
´1. Kronecker’s delta function, denoted by ij, is equal to one when i “ j and zero otherwise. Matlab index notation will occasionally be used for submatri- ces of a given matrix. We write Ap1:k,:) to denote the submatrix of A consisting of the rows from 1 to k of A.
Let A P CIˆJ, then A “ A p1 : I ´ 1, :q P CpI´1qˆJ and A “ A p2 : I, :q P CpI´1qˆJ, i.e., A and A are obtained by deleting the bottom and top row of A, respectively.
DkpAq P CJˆJ denotes the diagonal matrix holding row k of A P CIˆJ on its diagonal.
The k-th compound matrix of A P CIˆR is denoted by CkpAq P CCkIˆCkR, where Ckm“ k!pm´kq!m! . It is the matrix containing the determinants of all kˆk submatrices of A, arranged with the submatrix index sets in lexicographic order. For more elementary notation, we refer to [23].
Tensor decompositions: Let us briefly present the tensor decompositions that are used in this paper. See the references in [23] for further details.
1) Polyadic Decomposition (PD) and CPD: A PD is a decomposition of a tensor X P CIˆJˆK into rank-1 terms:
X “ ÿR r“1
arbbbbrbbbcr, (1) with factor matrices A “ ra1, . . . ,aRs P CIˆR, B “ rb1, . . . ,bRs P CJˆR and C “ rc1, . . . ,cRs P CKˆR. The rank of a tensor X is equal to the minimal number of rank-1 tensors arbbbbrbbbcr that yield X in a linear combination.
Assume that the rank of X is R, then (1) is the CPD of X.2) Coupled PD and coupled CPD: More generally, we say that a collection of tensors Xpnq P CInˆJnˆK, n P t1, . . . , Nu, admits an R-term coupled PD if each tensor Xpnq can be written as:
Xpnq“ ÿR r“1
apnqr bbbbpnqr bbbcr, n P t1, . . . , Nu, (2) with factor matrices Apnq“ rapnq1 , . . . ,apnqR s P CInˆR, Bpnq“ rbpnq1 , . . . ,bpnqR s P CJnˆR and C “ rc1, . . . ,cRs P CKˆR. We define the coupled rank of tXpnqu as the minimal number of coupled rank-1 tensors apnqr bbbbpnqr bbbcrthat yield tXpnqu in a linear combination. Assume that the coupled rank of tXpnqu is R, then (2) is the coupled CPD of tXpnqu.
II. Review of links between multidimensional harmonic retrieval and(coupled) CPD
The connections between MHR, CPD and coupled CPD will be reviewed in this section. More details can
be found in the companion paper [23]. In order to understand why coupled CPD is a natural framework for MHR, we will first review the link between 1D HR and CPD.
A. 1D HR and CPD
1) CPD interpretation of 1D HR: Consider the 1D HR factorization
X “ ACTP CIˆM, (3) where A P CIˆR is a Vandermonde matrix with columns ar “ r1, zr,z2r, . . . ,zI´1r sT and where C P CMˆR is an unconstrained matrix that has full column rank. Let us stack X and X as slices of Y P C2ˆJˆM. Because of the shift-invariance A “ AD2pAq, we have
Y “„ X X
⇢
“
„ BCT
BD2pAp2qqCT
⇢
“ pAp2qd BqCTP C2JˆM, (4) where Ap2q “ “1 ¨¨¨ 1
z1¨¨¨ zR
‰ P C2ˆR and B “ A P CJˆR with J “ I ´ 1. Equation (4) can be seen as the matrix representation of the Vandermonde constrained CPD of Y:
Y “ ÿR r“1
ap2qr bbbbrbbbcrP C2ˆJˆM. (5) It was explained in [23] that the Vandermonde structure in (5) can be relaxed without a↵ecting the uniqueness of the decomposition. In other words, the 1D HR factoriza- tion of X is unique if and only if the unconstrained CPD of Y is unique. The latter decomposition is known to be unique if B and C have full column rank and if zr, zs,
@r , s. It is known that under this condition, the CPD of Y can be determined via the GEVD of the matrix pencil pX, Xq (e.g., [17], [16], [6], [8]). This is in fact precisely what happens in ESPRIT, see the next subsection.
2) ESPRIT: A particular method for computing the 1D HR factorization of X via the matrix pencil pX, Xq in (4) is ESPRIT [16]. It was originally formulated in terms of the correlation matrix XXH “ ASTS˚AH, but it can also be formulated in terms of an SVD. Briefly, let the columns of U P CIˆR constitute a basis for rangepXq, numerically obtained from the SVD of X. Since B and C have full column rank, there exists a nonsingular matrix F P CRˆR such that U “A FT. The shift-invariance property A “ AD2pAp2qq now results in the EVD
U:U “ F´1D2pAp2qqF. (6) As the generators tzruRr“1 form the second row of Ap2q, they can be found as the eigenvalues of U:U, i.e., the generalized eigenvalues of the pencil pU, Uq. We have now found the Vandermonde column vectors of A that form a basis for rangepXq. Finally, the linear coefficients of the columns of X in this Vandermonde basis can be determined as C “ pA:XqT.
B. MHR and Vandermonde constrained CPD
1) VDM-CPD interpretation of MHR: Consider now an N-dimensional HR problem that involves R complex exponentials. The generator of the r-th exponential in the n-th mode is denoted by zr,n, 1 § r § R, 1 § n § N.
The data tensor X P CI1ˆ¨¨¨ˆINˆM admits a constrained PD
X “ ÿR r“1
ap1qr bbb¨ ¨ ¨bbbapNqr bbbcr, (7) with Vandermonde factor matrices Apnq “ rapnq1 , . . . ,apnqR s P CInˆR with apnqr “ r1, zr,n,z2r,n, . . . ,zIr,nn´1sT and unconstrained C “ rc1, . . . ,cRs P CMˆR, where M is the number of snapshots. The PD of X in (7) has the following matrix representation
X “ pAp1qd ¨ ¨ ¨ d ApNqqCTP Cp±Nn“1InqˆM. (8) Note that the Vandermonde constraints in (7) can in- crease the number of rank-1 terms that are required to fit X, compared to unconstrained CPD. The mini- mal number of Vandermonde constrained rank-1 terms ap1qr bbb¨ ¨ ¨bbbapNqr bbbcrthat yield X will be denoted by rMHRpXq.
2) Exploiting harmonic structure in one mode: If the Van- dermonde structure of only one of the factor matrices, say Apnq, is taken into account, then (7) corresponds to a VDM-CPD problem that can be reduced to a CPD as in Subsection II-A. In detail, using the row-selection matrices
SpIpnq1,...,INq“ I±n´1
p“1Ipb IInb I±Nq“n`1Iq , (9) SpIpnq1,...,INq“ I±n´1p“1Ipb IInb I±Nq“n`1Iq , (10) we build the tensor YpnqP C2ˆpˆn´1p“1IpqˆpIn´1qˆpˆNq“n`1IqqˆM:
Ypnq“ ÿR r“1
ap2,nqr bbbbpnqr bbbcr, (11) with matrix representation
Ypnq“
« SpIpnq1,...,INqX SpIpnq1,...,INqX
ff
“ pAp2,nqd BpnqqCTP C2JnˆM, (12) where Jn“ p±N
p“1,p,nIpqpIn´ 1q and Ap2,nq““ 1 ¨¨¨ 1
z1,n¨¨¨ zR,n
‰P C2ˆR, (13)
Bpnq“n´1ä
p“1
Appqd Apnqd äN
q“n`1
Apqq. (14) Ignoring the structure of Bpnq, (11) can be interpreted as an unconstrained two-slice CPD. Hence, if Bpnq and C have full column rank and if zr,n , zs,n, @r , s, then the generators tzr,nuRr“1 can be found via a (G)EVD. In a similar way as in ESPRIT for 1D HR, the remaining unknown factor matrices tApnqu and C in (8) can now be deduced. Details and variants of this procedure can be found in [4], [7], [9], [10], [27], [21], [5] and references therein.
C. MHR and coupled CPD
1) Coupled CPD interpretation of MHR: In contrast to the existing algebraic VDM-CPD approaches to MHR, coupled CPD simultaneously considers all tensors in the set tYp1q, . . . ,YpNqu. Overall, we obtain
»
—– Yp1q
... YpNq
fi ffifl “
»
—–
Ap2,1qd Bp1q ... Ap2,Nqd BpNq
fi
ffifl CT, (15)
where Ypnqis given by (12). The structure of (15) is such that X must have the N-dimensional harmonic structure (7), even when the structure of tBp1q, . . . ,BpNqu is ignored [23]. In fact, we make the assumption that Ap2,nq, n P t1, . . . , Nu do not have zero entries, i.e., we assume that Ap2,nqhave a (trivial) Vandermonde structure. Hence, in principle any numerical method for coupled CPD can be used to compute the MHR decomposition (7) via (15).
2) Exploiting several harmonic structures: Both neces- sary and sufficient MHR uniqueness conditions have been obtained in [23]. In this paper we develop a lin- ear algebra based method. Similar to ESPRIT, it finds the pure harmonics from rangepXq, and it does so via (G)EVD. Like ESPRIT, it is guaranteed to find the so- lution in the noise-free case under some algebraic as- sumptions. It relies on the sufficient and “almost nec- essary” MHR uniqueness condition stated in Theorem II.1, which jointly takes the Vandermonde structures of Ap1q, . . . ,ApNq into account. Consequently, it works under more relaxed conditions than the methods that essentially rely on the harmonic structure in one mode.
Theorem II.1 makes use of the matrix
GpNq“
»
——
—– C2
´Ap2,1q¯ d C2´
Bp1q¯ ...
C2
´Ap2,Nq¯ d C2´
BpNq¯ fi ffiffi
ffifl. (16)
Theorem II.1. Consider the PD of X P CI1ˆ¨¨¨ˆINˆM in (7) where the factor matrices tApnqu are Vandermonde. If
#Cin (8) has full column rank,
GpNq in (16) has full column rank, (17) then rMHRpXq “ R and the VDM-CPD of X is unique.
III. Coupled CPD algorithms for MHR
In this section we use coupled CPD to obtain ESPRIT- type algorithms for MHR that work under the assump- tions in Theorem II.1.
A. Find a basis for rangepAp1qd ¨ ¨ ¨ d ApNqq
Similar to ESPRIT [14], [16] we will find the individual multidimensional harmonics from rangepAp1q d ¨ ¨ ¨ d ApNqq. We first explain that under condition (17), the latter is readily available as rangepXq. First we note that GpNqhaving full column rank implies that Ap1qd¨ ¨ ¨dApNq
has full column rank. Indeed, if Ap1qd ¨ ¨ ¨ d ApNq does not have full column rank, then there exists a vector y P kerpAp1qd ¨ ¨ ¨ d ApNqq and a vector x P CM such that we obtain an alternative decomposition X “`Ap1qd
¨ ¨ ¨ d ApNq˘ `
C ` xyT˘T, contradicting Theorem II.1. Since C is assumed to have full column rank, we have that rangepXq “ rangepAp1qd ¨ ¨ ¨ d ApNqq. Numerically, we obtain rangepXq from the compact SVD X “ U⌃⌃⌃VH, where U P Cp±Nn“1InqˆR and V P CMˆR are columnwise orthonormal and ⌃⌃⌃ P CRˆR is diagonal. We know that rangepXq “ rangepU⌃⌃⌃q and that there exists a nonsingu- lar matrix F P CRˆR such that
U⌃⌃⌃“ pAp1qd ¨ ¨ ¨ d ApNqqFT. (18) Compared to (8), (18) involves the smaller pRˆRq matrix F. Note that this dimensionality reduction step is similar to the one in classical ESPRIT for 1D HR, see Subsection II-A. Note also that if R is unknown, then it can be estimated as the number of significant singular values of X.
B. From MHR to coupled CPD
From (18) we build the tensors Upnqred P C2ˆpˆn´1p“1IpqˆpIn´1qˆpˆNq“n`1IqqˆR, n P t1, . . . , Nu, each admitting a factorization
Upnqred:“
„ Up1,nq Up2,nq
⇢
“ pAp2,nqd BpnqqFT, (19) where ’red’ stands for reduced, and
Up1,nq“ SpIpnq1,...,INqU⌃⌃⌃P CMnpIn´1qˆR, (20) Up2,nq“ SpIpnq1,...,INqU⌃⌃⌃P CMnpIn´1qˆR, (21) in which Mn “ p±N
m“1,m,nImq such that Jn “ MnpIn ´ 1q. Note that (19), in combination with the way Upnqred is constructed, implements the harmonic structure in the nth mode. Overall, we obtain the compressed variant of (15):
»
—— –
Up1qred ... UpNqred
fi ffiffi fl “
»
—–
Ap2,1qd Bp1q ... Ap2,Nqd BpNq
fi
ffifl FT, (22)
in which each submatrix implements the harmonic struc- ture in the corresponding mode. In the next subsection we adapt the algebraic SD method in [1] to our MHR problem. In particular, a more efficient implementation will be developed that takes the harmonic structure of the MHR problem into account.
C. From coupled CPD to SD
From (22) we build R symmetric matrices tMprqu, admitting the factorizations
Mprq“ G⇤⇤⇤prqGT, r P t1, . . . , Ru, (23)
where G “ F´1 and ⇤⇤⇤prq P CRˆR are diagonal matrices.
The construction will be explained below. Note that (23) can be seen as the slice-wise representation of the CPD of the partially symmetric tensor M “∞R
r“1grbbbgrbbb⇣⇣⇣r
with factor matrices G and “ r⇣⇣⇣1, . . . , ⇣⇣⇣Rs in which
⇣⇣⇣r “ r prq11, . . . , prqRRsT. It can verified that both G and have full column rank [1]. This implies that in the exact case, it suffices to exploit the structure of only two of the matrix slices, denoted by Mpr1q and Mpr2q. More precisely, as in the ESPRIT method, the problem can be solved by means of a GEVD of the matrix pencil pMpr1q,Mpr2qq in which the columns of G´Tare the gen- eralized eigenvectors. In the inexact case, the robustness can be increased by simultaneously exploiting all R equations in (23), using for instance optimization-based methods (e.g., [28], [20]). The CPD of M can equivalently be seen as the simultaneous diagonalization of its matrix slices, which is why (23) is referred to as SD.
Following the steps in [25] we will now explain how to construct the matrices tMprqu in (23). Thereafter, we explain that the complexity of the construction can be reduced by exploiting the MHR structure.
Construction of matrices tMprqu: Consider the bilinear mappings pnq: CJnˆ2ˆ CJnˆ2Ñ C4J2n defined by [1]:
pnqpX, Yqpi´1q4Jn`4pj´1q`2pk´1q`l
“ xikyjl` yikxjl´ xilyjk´ yilxjk, n P t1, . . . , Nu.
Denote pnqr,s “ pnqprup1,nqr ,up2,nqr s, rup1,nqs ,up2,nqs sq, in which up1,nqt and up2,nqt correspond to the t-th column of Up1,nq and Up2,nq, respectively. Define pnq “ r pnq1 , 2¨ pnq2 s P C4J2nˆC2R`1 with pnq1 “ r pnq1,1, pnq2,2, . . . , pnqR,Rs P C4Jn2ˆR and
pnq
2 “ r pnq1,2, pnq1,3, pnq2,3, . . . , pnqR´1,Rs P C4J2nˆC2R. It can be verified that if condition (17) in Theorem II.1 is satisfied, then
“
»
—–
p1q
...
pNq
fi
ffifl P C4p∞Nn“1J2nqˆC2R`1 (24)
has an R-dimensional kernel [25]. Furthermore, the ker- nel vectors can be associated with the matrices tMprqu.
Namely, Mprqsatisfies
mprq“ 0 , @r P t1, . . . , Ru,
where mprq “ rmprq1,1,mprq2,2, . . . ,mprqR,R,mprq1,2,mprq1,3, . . . ,mprqR´1,RsT. Since kerp q “ kerp H q, a basis for the kernel of (24) can be found from the smaller (C2R`1ˆ C2R`1) positive semi-definite matrix
H “
ÿN n“1
pnqH pnq. (25)
Numerically, we can take tmprqu equal to the eigen- vectors associated with the R smallest eigenvalues of
H . In the computation of the smaller matrix (25), the structure of t pnqH pnqu can be exploited. Concretely, define the set SR “ tps1,s2q : 1 § s1 † s2 § Ru in which the C2R elements are ordered lexicographically.
Consider the mapping fR : N2 Ñ t1, 2, . . . , C2Ru that returns the position of its argument in the set SR, e.g., fRpp1, 3qq “ 3. Define also Yp1,nq P CRˆR, Yp2,nq P CRˆC2R and Yp3,nqP CC2RˆC2R as follows
pYp1,nqqr,s“ 8pwpnqrs ´ bpnqHrs bpnqsr q , (26) pYp2,nqqr, fRps1,s2q“ 16pwpnqrs1wpnqrs2 ´ bpnqHrs1 bpnqrs2q , (27) pYp3,nqqfRpr1,r2q, fRps1,s2q“
4pwpnqr1s1wpnqr2s2` wpnqr2s1wpnqr1s2´ bpnqHr1s1 bpnqr2s2´ bpnqHr2s1 bpnqr1s2q , (28) where
wpnqrs :“ wp1,nqrs ` wp2,nqrs ,
wp1,nqrs “ up1,nqHr up1,nqs , (29) wp2,nqrs “ up2,nqHr up2,nqs , (30)
bpnqrs :“ rwp1,nqrs ,zpnqrs ,zpnq˚sr ,wp2,nqrs sT,
zpnqrs “ up2,nqHr up1,nqs . (31) It can be verified that pnqH pnq““ Yp1,nq 2¨Yp2,nq
2¨Yp2,nqH4¨Yp3,nq
‰, so that (25) can be expressed in terms of (26)–(28):
H “
ÿN n“1
pnqH pnq“ ÿN n“1
„ Yp1,nq 2 ¨ Yp2,nq 2 ¨ Yp2,nqH 4 ¨ Yp3,nq
⇢ . (32) Note that the matrices Yp1,nq and Yp3,nq are Hermitian, implying that it suffices to compute their upper trian- gular part. Note also that the vector products bpnqHr1s1 bpnqr2s2 that appear in (26)–(28) can be computed prior to the construction of Yp1,nq, Yp2,nq and Yp3,nq. In addition, due to the MHR structure, the involved scalars twp1,nqrs ,wp2,nqrs u can be computed efficiently, as explained next.
Exploiting MHR structure: By construction the matrix slices tUp1,nqu of the tensors tUredpnqu are largely overlap- ping. (Similarly for tUp2,nqu.) More concretely, observe that for every n P t1, . . . , Nu, the vector up1,nqr “ SpIpnq1,...,INqur
consists of all except Mnentries of ur. To put it di↵erently, the vectors in the set tup1,1qr , . . . ,up1,Nqr u are almost identi- cal. (Similarly for the vectors in the set tup2,1qr , . . . ,up2,Nqr u, which all are related via up2,nqr “ SpIpnq1,...,INqur.) Moreover, Up1,nqand Up2,nqin (19) are taken from the matrix U⌃⌃⌃, the columns of which are mutually orthogonal, see (20)–(21).
These properties enable a faster computation of wp1,nqrs
and wp2,nqrs in (29)–(30). We work as follows. Consider the following row-selection matrices complementary to SpIpnq1,...,INq and SpIpnq1,...,INq:
TpIpnq1,...,INq“ I±n´1
p“1Ipb IInpIn, :q b I±Nq“n`1Iq P CMnˆp±Nn“1Inq , TpIpnq1,...,INq“ I±n´1p“1Ipb IInp1, :q b I±Nq“n`1IqP CMnˆp±Nn“1Inq . The selection is complementary in the sense that SpIpnq1,...,INq` TpIpnq1,...,INq “ SpIpnq1,...,INq ` TpIpnq1,...,INq “ I±Nn“1In, i.e., TpIpnq1,...,INq and TpIpnq1,...,INq select the Mn rows not selected by SpIpnq1,...,INqand SpIpnq1,...,INq, respectively. This property together
with relations (20)–(21) allows us to express (29)–(30) as follows:
wp1,nqrs “ up1,nqHr up1,nqs “ r suHrSpIpnq1,...,INqSpIpnq1,...,INqus
“ r sp rs´ uHrTpIpnq1,...,INqTpIpnq1,...,INqusq, (33) wp2,nqrs “ up2,nqHr up2,nqs “ r suHrSpIpnq1,...,INqSpIpnq1,...,INqus
“ r sp rs´ uHrTpIpnq1,...,INqTpIpnq1,...,INqusq, (34) where the orthogonality uHrus“ rshas been taken into account and where r denotes the r-th singular value of X.
Denoting the Hermitian matrices
ÄWp1,nq“ pTpIpnq1,...,INqU⌃⌃⌃qH¨ pTpIpnq1,...,INqU⌃⌃⌃q P CRˆR, (35) ÄWp2,nq“ pTpIpnq1,...,INqU⌃⌃⌃qH¨ pTpIpnq1,...,INqU⌃⌃⌃q P CRˆR, (36) the products in (29)–(31) can be computed as the entries of
Wp1,nq“ Up1,nqHUp1,nq“ ⌃⌃⌃2´ ÄWp1,nqP CRˆR, (37) Wp2,nq“ Up2,nqHUp2,nq“ ⌃⌃⌃2´ ÄWp2,nqP CRˆR, (38)
Zpnq“ Up2,nqHUp1,nqP CRˆR, (39)
where n P t1, . . . , Nu. Note that by exploiting the struc- ture of up1,nqHr up1,nqs and up2,nqHr up2,nqs , the computational cost of wp1,nqrs and wp2,nqrs have both been reduced from 2MnpIn´ 1q ´ 1 flops to 2Mn flops.
The di↵erent steps in the transformation of the original MHR decomposition (8) into the SD factorization (23) are listed in Table I, together with an estimate of their com- plexity (in flops). t is assumed that the R-SVD and QZ methods are used for the computation of the SVD and GEVD, respectively [3]. The most expensive steps are the initial dimensionality reduction step, which has a com- plexity of the size of 6p±N
n“1InqM2, and the computation of the matrices tZpnqu used in the construction of H , with a complexity of the order of 2Np±N
n“1InqR2. The to- tal complexity is of the order of p6M`2NRM2qp±N
n“1InqM.
If there are many snapshots (M °° R), then the cost is dominated by the conventional SVD reduction in Step 1. If M w R, then the cost is about p6 ` 2NqR flops per tensor entry.
D. From SD to single-tone MHR problems
After G “ F´1 has been found from (23) we obtain U⌃⌃⌃GT “ Ap1qd ¨ ¨ ¨ d ApNq from relation (18). The gen- erators tzr,nuNn“1associated with the rth source can now be found from the rth column of U⌃⌃⌃GT. In other words, coupled CPD has turned the MHR problem (7) into a set of decoupled single-tone MHR problems. Denote
pnq“ SpIpnq1,...,INqU⌃⌃⌃GT , (40)
pnq“ SpIpnq1,...,INqU⌃⌃⌃GT . (41)
Algorithm operations Eq. Complexity order (flops)
SVD of X (8) 6p±N
n“1InqM2` 11M3 Compute tÄWp1,nq, ÄWp2,nquNn“1 (35)–(36) 2p2p∞N
n“1Mnq ´ NqC2R`1 Compute tWp1,nq,Wp2,nquNn“1 (37)–(38) 2NR
Compute tZpnquNn“1 (39) p2p∞N
n“1MnpIn´ 1qq ´ NqR2 Compute
tbpnqHr1s1 bpnqr2s2u1§n§N1§r1,r2,s1,s2§R
(26)–(28) 7NR4 Compute
tYp1,nq,Yp2,nq,Yp3,nquNn“1 (26)–(28) Np2C2R`1` 3RC2R` 6C2R`1q Compute H (32) NpC2R`1` 2RC2R` 2C2R`1q
EVD of H (32) 43pC2R`1q3
GEVD of pair pMpr1q,Mpr2qq (23) 30R2 TABLE I
An estimate of the complexity of the different steps in the transformation from the originalMHR decomposition (8) into the
SD factorization (23). We denote Mn“ p±
m“1,m,nImq.
Due to the shift-invariance property Apnq “ ApnqD1przr,1,, . . . ,zR,nsq, we have pnqr “ pnqr zr,n
so that the generator can be obtained as zr,n“ p pnqHr pnq
r q{p pnqHr pnq r q.
In several applications, such as DOA estimation, zr,n
is located on the unit circle, i.e., zr,n“ ei¨ r,n“ cosp r,nq ` i sinp r,nq for some r,nP r0, 2⇡q. In that case, we can find the generator by minimizing
f p r,nq “››› pnqr zr,n´ pnqr
››
›2F “››› pnqr ei¨ r,n´ pnqr
››
›2F
“ ↵r,n´ 2Re t r,nu cosp r,nq ´ 2Imt nu sinp r,nq
“ ↵r,n´ 2Re t r,nu`cos2p r,n{2q ´ sin2p r,n{2q˘
´ 4Imt r,nu sinp r,n{2q cosp r,n{2q
“ ↵r,n´ 2
„ cosp r,n{2q sinp r,n{2q
⇢T
Gpr,nq
„ cosp r,n{2q sinp r,n{2q
⇢ , where ↵r,n “ pnqHr pnqr ` pnqHr pnq
r , r,n “ pnqHr pnqr and Gpr,nq “ “´Retr,nu Imtr,nu
Imtr,nu Retr,nu
‰. Let xpr,nq P R2 denote the normalized eigenvector associated with the dominant eigenvalue of Gpr,nq, then the unit-norm generator is obtained as zr,n“ pxpr,nq1 ` i ¨ xpr,nq2 q2.
E. Summary and discussion of overall method
The overall procedure is summarized as Algorithm 1.
If condition (17) in Theorem II.1 is satisfied, then Algo- rithm 1 is guaranteed to find the MHR decomposition in the exact case. On the other hand, if condition (17) is not satisfied, then Algorithm 1 does not apply. In the latter case, what goes wrong is that we cannot find the R symmetric matrices tMprqu in (23) from H , as the dimension of its kernel is larger than R [25]. It is also clear that Algorithm 1 works under more relaxed conditions than the approaches in Subsection II-B since it takes all harmonic structures into account at once. See [23] for a thorough discussion of the conditions.
In the inexact case, Algorithm 1 does not necessarily provide the Least Squares (LS) estimate of the factors of the original MHR decomposition (7) as they would be
obtained from direct fitting. However, some robustness is built into the algorithm. More precisely, in Steps 1 and 5, LS estimates of rangepAp1qd¨ ¨ ¨dApNqq and the kernel of H can be obtained via SVD and EVD, respectively.
Likewise, in Step 8 a LS estimate of the generators tzr,nu can be obtained via the closed-form solutions discussed in Subsection III-D. As mentioned in Subsection III-C, an algebraic solution to the SD problem (23) in Step 6 can be obtained via a GEVD. The result can serve as a good initialization for an SD method that in turn takes all the matrix slices in (23) into account in LS sense. Obviously, Algorithm 1 as a whole can be used to initialize a LS fitting method for MHR (e.g., [20]), i.e., the optimization algorithm is only used to refine the estimates.
Algorithm 1 SD method for MHR (Theorem II.1).
Input: X “ pAp1qd ¨ ¨ ¨ d ApNqqCT of the form (8).
1.Compute SVD X “ U⌃⌃⌃VH(and find R from ⌃⌃⌃).
2. Build Upnqred in (19) for n P t1, . . . , Nu.
3. Construct tWp1,nq,Wp2,nq,ZpnquNn“1given by (37)–(39).
4. Construct H in (32) viatWp1,nq,Wp2,nq,ZpnquNn“1. 5. Obtain matrices tMprqu from the kernel of H . 6. Solve SD problem Mprq“ G⇤⇤⇤prqGT, r P t1, . . . , Ru.
7. Compute N “ U⌃⌃⌃GT.
8. Obtain tzr,nu from N via single-tone MHR.
Output: tzr,nu
F. Variant 1: Compression
Note that the submatrices in (22) have the size (JnˆR), with Jn w ±N
n“1In. We may consider compressing the tensors Upnqred. Ideally, Up1,nq and Up2,nq in (20)–(21) will be replaced by matrices of size pR ˆ Rq, as there are only R harmonics. It makes sense to project on the R- dimensional dominant subspace of rUp1,nq,Up2,nqs, which may be computed by means of an SVD. We may equiva- lently compute the R-dimensional dominant eigenspace of
rUp1,nq,Up2,nqsH¨ rUp1,nq,Up2,nqs (42)
“
„ Wp1,nq ZpnqH Zpnq Wp2,nq
⇢
P C2Rˆ2R, (43)
where Wp1,nq, Wp2,nq and Zpnq are the matrices that al- ready appeared in (37)–(39). While the diagonal blocks can be computed cheaply thanks to the orthogonality of U, the computation of Zpnq is expensive. Recall from Subsection III-C that the computation of Zpnq was also the most expensive part of the construction of H . Consequently, although replacing Up1,nq and Up2,nq by (RˆR) matrices makes the computational cost of Steps 3- 4 in Algorithm 1 negligible, the savings are compensated by the cost of the compression itself, i.e., the complexity is about the same as that of Algorithm 1. However, the “compressed variant” of Algorithm 1 has another interesting feature. Recall from Subsection III-B that we use Uredpnq to impose the harmonic structure in the nth
