BLIND SEPARATION OF EXPONENTIAL POLYNOMIALS AND THE DECOMPOSITION OF A TENSOR IN RANK-( L

BLIND SEPARATION OF EXPONENTIAL POLYNOMIALS AND THE DECOMPOSITION OF A TENSOR IN RANK-( L

, L

, 1) TERMS

LIEVEN DE LATHAUWER^†

Abstract. We present a new necessary and sufficient condition for essential uniqueness of the decomposition of a third-order tensor in rank-(Lr, Lr, 1) terms. We derive a new deterministic technique for blind signal separation that relies on this decomposition. The method assumes that the signals can be modeled as linear combinations of exponentials or, more generally, as exponential polynomials. The results are illustrated by means of numerical experiments.

Key words. multilinear algebra, higher-order tensor, singular value decomposition, canonical polyadic decomposition, block term decomposition, blind signal separation, exponential polynomial

AMS subject classifications. 15A18, 15A69 DOI. 10.1137/100805510

1. Introduction. In section 1.1 we explain our notation and give some basic definitions. In section 1.2 we briefly recall the multilinear singular value decomposition (MLSVD) and the canonical polyadic decomposition (CPD). In section 1.3 we set out the goals of this paper and explain how it is organized.

1.1. Preliminaries.

1.1.1. Notation. Scalars are denoted by lower-case letters (a, b, . . . ), vectors are written in boldface lower-case (a, b, . . . ), matrices correspond to boldface capitals (A, B, . . . ), and tensors are written as calligraphic letters ( A, B, . . . ). This notation is consistently used for lower-order parts of a given quantity. For instance, the entry with row index i and column index j in a matrix A, i.e., (A)

, is denoted by a

(also (a)

= a

and ( A)

= a

). The ith column vector of a matrix A is denoted by a

, i.e., A = [a

a

. . .]. Italic capitals indicate index upper bounds (e.g., i = 1, 2, . . . , I).

K stands for R or C. The symbol ⊗ denotes the Kronecker product,

A ⊗ B =

⎡

⎢ ⎣

a

B a

B · · · a

B a

B · · ·

.. . .. .

⎤

⎥ ⎦ .

The Khatri–Rao or columnwise Kronecker product is represented by :

A  B = [a

⊗ b

a

⊗ b

· · · ] .

∗Received by the editors August 16, 2010; accepted for publication (in revised form) September 13, 2011; published electronically December 8, 2011. This research was supported by the Research Coun- cil K.U. Leuven: GOA-MaNet, CoE EF/05/006 Optimization in Engineering (OPTEC), CIF1 and STRT1/08/023; the F.W.O.: (a) project G.0427.10N, (b) Research Communities ICCoS, ANMMM, and MLDM; the Belgian Federal Science Policy Office: IUAP P6/04 (DYSCO, “Dynamical systems, control and optimization,” 2007–2011); EU: ERNSI.

http://www.siam.org/journals/simax/32-4/80551.html

†Group Science, Engineering and Technology, Katholieke Universiteit Leuven Campus Kortrijk, E.

Sabbelaan 53, 8500 Kortrijk, Belgium (Lieven.DeLathauwer@kuleuven-kortrijk.be), and Department of Electrical Engineering (ESAT), Research Division SCD, Katholieke Universiteit Leuven, Kasteel- park Arenberg 10, B-3001 Leuven, Belgium (Lieven.DeLathauwer@esat.kuleuven.be, http://homes.

esat.kuleuven.be/∼delathau/home.html).

1451

The column space of a matrix is denoted by span(A). The rank of a matrix A is de- noted by rank(A). The superscripts ·

, ·

, ·

, and ·

denote the transpose, complex conjugate, complex conjugated transpose, and Moore–Penrose pseudoinverse, respec- tively. The operator diag( ·) stacks its scalar arguments in a square diagonal matrix.

For vectorization of a matrix A = [a

a

. . .] we adopt the following convention:

vec(A) = [a

a

. . .]

. The (N × N ) identity matrix is represented by I

. 1.1.2. Basic definitions.

Definition 1.1. Consider T ∈ K

and U

∈ K

, U

∈ K

, U

∈ K

. Then the products T ·

U

, T ·

U

, and T ·

U

are defined by

( T ·

U

)

=

I1



i1=1

t

u

∀j

, i

, i

,

( T ·

U

)

=

I2



i2=1

t

u

∀i

, j

, i

,

( T ·

U

)

=

I3



i3=1

t

u

∀i

, i

, j

,

respectively.

Definition 1.2. The outer product A ◦ B of a tensor A ∈ K

and a tensor B ∈ K

is the tensor defined by

( A ◦ B)

= a

b

for all values of the indices.

For instance, the outer product T of three vectors a, b, and c is defined by t

= a

b

c

for all values of the indices. The outer product T of a matrix E and a vector c is defined by t

= e

c

for all values of the indices.

Definition 1.3. In the matrix representations T

∈ K

, T

∈ K

, T

∈ K

of a tensor T ∈ K

, the entries are stacked as follows:

(T

)

= (T

)

= (T

)

= t

∀i

, i

, i

. (1.1)

Definition 1.4. A third-order tensor T ∈ K

has multilinear rank (R

, R

, R

) iff rank(T

) = R

, rank(T

) = R

, and rank(T

) = R

.

The values R

, R

, and R

are sometimes called the mode-1, mode-2, and mode- 3 rank, respectively. If the multilinear rank of a third-order tensor T is equal to (R

, R

, R

), then T is called rank-(R

, R

, R

). A rank-(1, 1, 1) tensor is briefly called rank-1. This implies that a third-order tensor has rank-1 iff it equals the outer product of three nonzero vectors.

Besides multilinear rank, a second generalization of matrix rank to higher-order tensors is the following.

Definition 1.5. The rank of a tensor T is the minimal number of rank-1 tensors that yield T in a linear combination.

An early discussion of the generalization of matrix rank to higher-order tensors

can be found in [12, 13].

1.2. Basic tensor decompositions.

Definition 1.6. An MLSVD of a rank-(R

, R

, R

) tensor T ∈ K

is a decomposition of T of the form

(1.2) T = S ·

U

·

U

·

U

, in which

• the matrices U

∈ K

, U

∈ K

, and U

∈ K

are orthogonal (unitary);

• the tensor S ∈ K

is – all-orthogonal (all-unitary):



i2=1 I3



i3=1

s

s

= σ

if j

= j

= i

= 0 if j

= j

,



i1=1 I3



i3=1

s

s

= σ

if j

= j

= i

= 0 if j

= j

,



i1=1 I2



i2=1

s

s

= σ

if j

= j

= i

= 0 if j

= j

. – ordered:

σ

 σ

 · · ·  σ

> σ

= · · · = σ

= 0, σ

 σ

 · · ·  σ

> σ

= · · · = σ

= 0, σ

 σ

 · · ·  σ

> σ

= · · · = σ

= 0.

The values σ

, σ

, and σ

are called mode-1, mode-2, and mode-3 singular values, respectively. The columns of U

, U

, and U

are called mode-1, mode-2, and mode-3 singular vectors, respectively.

The MLSVD has its roots in [30, 31] and was further studied in [5], where it was called higher-order singular value decomposition (HOSVD).

Definition 1.7. A canonical polyadic decomposition (CPD) of a rank-R tensor T ∈ K

is a decomposition of T in a linear combination of R rank-1 terms:

(1.3) T =



a

◦ b

◦ c

.

The fully symmetric variant, in which a

= b

= c

, r = 1, . . . , R, was studied

in the nineteenth century in the context of invariant theory [3]. The unsymmetric

decomposition was introduced in 1927 [12, 13]. Around 1970, the unsymmetric de-

composition was independently reintroduced in psychometrics [1] and phonetics [11],

where it was called canonical decomposition (CANDECOMP) and parallel factor de-

composition (PARAFAC), respectively. The terms CANDECOMP and PARAFAC

are sometimes used when the number of rank-1 terms is not minimal.

1.3. Goals and organization. In this paper we further study a tensor decom- position that we have recently introduced, namely the decomposition of a third-order tensor in rank-(L

, L

, 1) terms [6, 7, 8], and we develop a new technique for blind signal separation, based on this decomposition.

Section 2 contains the discussion of the tensor decomposition as such. In section 2.1 the definition is recalled. Section 2.2 concerns the uniqueness of the decomposition.

Theorem 2.4 is a new result. Section 2.3 is a note on computational methods.

Blind signal separation has been an active research area in the last twenty years;

see, for instance, the books [2, 4, 14]. A prominent technique is independent compo- nent analysis (ICA). In ICA, signals are separated on the basis of their hypothesized statistical independence. In this paper we develop a new deterministic technique for blind signal separation. Signals are separated under the assumption that they can be modeled as exponential polynomials. Such signals are ubiquitous; as a matter of fact, exponentials play a very fundamental role in signal processing and system theory [17, 24].

Section 3 derives the blind signal separation technique at a conceptual level.

The problem is formulated in terms of the decomposition in rank-(L

, L

, 1) terms of a third-order tensor with partial Hankel structure. Section 4 considers in detail signals that can be modeled as a linear combination of exponentials. The results are generalized to exponential polynomials in section 5. In section 6 the computational cost is reduced by means of tensor compression. Section 7 addresses the particular case of exponential polynomials that are mutually orthogonal. Section 8 illustrates the results by means of some numerical experiments. We conclude in section 9.

2. The decomposition of a third-order tensor in rank-( L

, L

, 1) terms.

2.1. Definition and representation. In [6, 7, 8] we introduced block term decompositions (BTD) of a higher-order tensor. If there is only one term, then the BTD of a given tensor reduces to its MLSVD. If all the terms have multilinear rank (1, 1, 1), then the BTD of a third-order tensor reduces to its CPD. A particular BTD is the decomposition in rank-(L

, L

, 1) terms.

Definition 2.1. A decomposition of a tensor T ∈ K

in a sum of rank- (L

, L

, 1) terms, 1  r  R, is a decomposition of T of the form

(2.1) T =



E

◦ c

,

in which matrix E

∈ K

is rank-L

and vector c

∈ K

is nonzero, 1  r  R.

We assume that R is minimal.

Note that the multilinear rank of the rth term in (2.1) is indeed equal to (L

, L

, 1).

It is clear that in (2.1) one can permute the rth and r

th terms when L

= L

. Also, one can scale E

, provided c

is counterscaled. We call the decomposition essentially unique when it is only subject to these trivial indeterminacies. If E

= U

· Σ

· V

denotes the SVD of E

, 1  r  R, then

T =



 U

· (c

 Σ

) · V

◦ (c

/c

) =



( c

 Σ

) ·

U

·

V

·

(c

/c

)

is a representation of (2.1) in which each term is in MLSVD form.

If we factorize E

as A

· B

, in which the matrix A

∈ K

and the matrix

B

∈ K

are rank-L

, r = 1, . . . , R, then we can write (2.1) as

(2.2) T = 

r=1

(A

· B

) ◦ c

.

The decomposition is visualized in Figure 2.1. Define A = [A

· · · A

], B = [B

· · · B

], C = [c

· · · c

]. In terms of the matrix representations of T defined in (1.1), (2.2) can be written as

T

= (C  B) · A

, (2.3)

T

= (C  A) · B

, (2.4)

T

= [vec(E

) · · · vec(E

)] · C

. (2.5)

T

I I

I

J J

J K K

K

= L

+ . . . + L

A

B

c

A

B

c

Fig. 2.1. Visual representation of the decomposition in rank-(Lr, Lr, 1) terms.

2.2. Uniqueness. In [7] several conditions were mentioned under which essential uniqueness of the decomposition is guaranteed. We recall [7, Theorem 4.1], which also appeared in a slightly different form in [26].

Theorem 2.2. Consider a decomposition of T ∈ K

in rank-(L

, L

, 1) terms as in (2.1)–(2.2), with I

, I



r=1

L

. If A = [A

A

· · · A

] and B = [B

B

· · · B

] have full column rank and C does not have proportional columns, then the decomposition is essentially unique.

Below we propose a new uniqueness theorem (Theorem 2.4). This theorem gener- alizes the CPD uniqueness result [15, Condition A]. We first reformulate [15, Condition A] as Theorem 2.3 and give a new proof. The reasoning will be generalized to prove Theorem 2.4.

Theorem 2.3. Consider the CPD (1.3) of T ∈ K

. Define E(w) =

r=1

w

a

· b

. Assume that the following conditions are satisfied:

(C1) For every w that has at least two nonzero entries, we have that rank(E(w)) >

1.

(C2) The columns of C = [c

c

· · · c

] are linearly independent.

Then the CPD of T is essentially unique. On the other hand, if condition (C1) is not satisfied, then the CPD of T is not essentially unique.

Proof. Assume that there exists an alternative decomposition

(2.6) T =



˜

a

◦ ˜b

◦ ˜c

.

Denote Z = C

, ˜ Z = ˜ C

, and Y = ˜ C

C.

We have

T ·

˜ z

= ˜ a

· ˜b

=



w

a

· b

,

where w

= y

. Because of condition (C1), we have that ˜ a

· ˜b

= d

a

·b

for some r ∈ {1, 2, . . . , R}, with d

= 0. Because of condition (C2), we have that ˜z

= d

z

. If r = 1, the first and rth terms in (1.3) can without loss of generality be switched. We conclude that ˜ a

· ˜b

= d

a

· b

and ˜ z

= d

z

, with d

= 0.

Induction step. Let us assume that ˜ a

· ˜b

= d

a

·b

and ˜ z

= d

z

with d

= 0, 1  r  ¯ R − 1. Then we prove that also ˜ a

· ˜b

= d

a

· b

and ˜ z

= d

z

with d

= 0. We have

T ·

˜ z

= ˜ a

· ˜b

=



w

a

· b

,

where w

= y

. Because of condition (C1), we have that ˜ a

· ˜b

= d

a

· b

for some r ∈ {1, 2, . . . , R}, with d

= 0. Because of condition (C2), we have that

˜

z

= d

z

. Because ˜ Z has full column rank, we can rule out the possibility that r < ¯ R. If r > ¯ R, the ¯ Rth and rth terms in (1.3) can without loss of generality be switched. The induction follows.

We now have that ˜ Z = ZD, in which D = diag(d

, d

, . . . , d

) ∈ K

is diagonal and nonsingular. This implies that ˜ C = CD

. On the other hand, we also have that

˜

a

· ˜b

= d

a

· b

, 1  r  R. We conclude that decompositions (1.3) and (2.6) are essentially equal.

Finally, we prove that condition (C1) is necessary. Let us assume that it is not satisfied. Then there exist nonzero vectors u and v such that u ·v

=

r=1

w

a

·b

, in which, say, w

= 0 = w

. Without loss of generality, we assume that w

= 1. Then we have that a

· b

= u · v

−

r=2

w

a

· b

. Hence, an alternative decomposition of T is

T = u ◦ v ◦ c

+



a

◦ b

◦ (c

− w

c

),

in which c

− w

c

is not proportional to c

.

We now generalize Theorem 2.3 to the decomposition in rank-(L

, L

, 1) terms.

Theorem 2.4. Consider a decomposition of T ∈ K

in rank-(L

, L

, 1) terms as in (2.1)–(2.2). Define E(w) =

r=1

w

E

. Assume that the following conditions are satisfied:

(C1) For every w that has at least two nonzero entries, we have that rank(E(w)) >

max

(L

).

(C2) The columns of C are linearly independent.

Then decomposition (2.1)–(2.2) is essentially unique. On the other hand, if condition (C1) is not satisfied, then decomposition (2.1)–(2.2) is not essentially unique.

Proof. Without loss of generality, we assume that the terms in (2.1)–(2.2) are ordered such that L

 L

 · · ·  L

. Condition (C1) is then equivalent to the following set of conditions:

(C1.1) If w

= 0 and ∃s ∈ {2, 3, . . . , R} such that w

= 0, then rank(E(w)) > L

.

(C1.2) If w

= 0 and ∃s ∈ {1, 3, . . . , R} such that w

= 0, then rank(E(w)) > L

. .. .

(C1.R) If w

= 0 and ∃s ∈ {1, 2, . . . , R−1} such that w

= 0, then rank(E(w)) > L

. Assume that there exists an alternative decomposition

(2.7) T =



E ˜

◦ ˜c

=



( ˜ A

· ˜ B

) ◦ ˜c

,

which is also ordered such that L

 L

 · · ·  L

. Denote Z = C

, ˜ Z = ˜ C

, and Y = ˜ C

C.

We have

T ·

˜ z

= ˜ E

=



w

E

,

where w

= y

. Since rank( ˜ E

) = L

, we have according to condition (C1.1) that (i) w

= w

= · · · = w

= 0 or that (ii) w

= 0 and ∃s ∈ {2, 3, . . . , R} such that w

= 0. In the special case that w

= w

= · · · = w

= w

= · · · = w

= 0 while w

= 0 and L

= L

, the first and ¯ sth terms in (2.1) can without loss of generality be switched. In all other cases, possibility (ii) can be ruled out because of (C1.2)–(C1.R) and the ordering constraint. We conclude that w

= w

= · · · = w

= 0, i.e., ˜ z

= d

z

with d

= 0.

Next, we have

T ·

˜ z

= ˜ E

=



w

E

,

where w

= y

. Since rank( ˜ E

) = L

, we have according to condition (C1.2) that (i) w

= w

= · · · = w

= 0 or that (ii) w

= 0 and ∃s ∈ {1, 3, . . . , R} such that w

= 0. Let us examine the latter possibility. First, it is impossible that w

= 0 while w

= w

= · · · = w

= 0. Indeed, this would mean that ˜ z

is in the orthogonal complement of span([c

c

· · · c

]). On the other hand, ˜ z

∈ span( ˜ Z) = span(Z) = span( ˜ C) = span(C) = span(T

2I1×I3

). This would imply that

˜

z

, just like ˜ z

, is proportional to z

, which is impossible since ˜ Z has full column rank.

Second, in the special case that w

= w

= · · · = w

= w

= · · · = w

= 0 while w

= 0 and L

= L

, the second and ¯ sth terms in (2.1) can without loss of generality be switched. Third, rank( ˜ E

) = L

rules out the special case w

= w

= · · · = w

= w

= · · · = w

= 0, w

= 0 and L

> L

. Fourth, because of (C1.3)–(C1.R) it is impossible that w

= 0, w

= 0, and ∃s ∈ {3, . . . , R} such that w

= 0. Fifth, because of (C1.3)–(C1.R) it is impossible that w

= w

= 0 and ∃s

, s

∈ {3, . . . , R}, s

= s

, such that w

= 0 = w

. We conclude that w

= w

= · · · = w

= 0, i.e., ˜ z

= d

z

with d

= 0.

Induction step. Let us assume that ˜ z

= d

z

with d

= 0, 1  r  ¯ R − 1. Then we prove that also ˜ z

= d

z

with d

= 0. We have

T ·

˜ z

= ˜ E

=



w

E

,

where w

= y

. Since rank( ˜ E

) = L

, we have according to condition (C1. ¯ R) that (i) w

= · · · = w

= w

= · · · = w

= 0 or that (ii) w

= 0 and ∃s ∈ {1, . . . , ¯ R − 1, ¯ R + 1, . . . , R} such that w

= 0. Let us examine the latter possibility.

First, it is impossible that w

= w

= · · · = w

= 0. Indeed, this would mean that ˜ z

is in the orthogonal complement of span([c

c

· · · c

]). On the other hand, ˜ z

∈ span(˜Z) = span(Z) = span( ˜ C) = span(C) = span(T

2I1×I3

). This would imply that ˜ z

∈ span([z

· · · z

]) = span([˜ z

· · · ˜z

]), which is impossible since Z has full column rank. Second, in the special case that w ˜

= · · · = w

(= · · · ) = w

= w

= · · · = w

= 0 while w

= 0 and L

= L

, the ¯ Rth and ¯ sth terms in (2.1) can without loss of generality be switched. Third, rank( ˜ E

) = L

rules out the special case w

= w

= · · · = w

= w

= · · · = w

= 0, w

= 0 and L

>

L

. Fourth, because of (C1. ¯ R + 1)–(C1.R) it is impossible that ∃s

∈ {1, . . . , ¯ R − 1}

such that w

= 0, w

= 0 and ∃s

∈ { ¯ R + 1, . . . , R} such that w

= 0. Fifth, because of (C1. ¯ R + 1)–(C1.R) it is impossible that w

= w

= · · · = w

= 0 and ∃s

, s

∈ { ¯ R + 1, . . . , R}, s

= s

, such that w

= 0 = w

. We conclude that w

= · · · = w

= w

= · · · = w

= 0, i.e., ˜ z

= d

z

with d

= 0.

After induction we have that ˜ Z = ZD, in which D ∈ K

is diagonal and nonsingular. This implies that ˜ C = CD

. From (2.5) we now obtain

T

· ˜Z = [vec(˜E

) · · · vec(˜E

)] = [vec(E

) · · · vec(E

)] · D.

We conclude that decompositions (2.1) and (2.7) are essentially equal.

Finally, we prove that conditions (C1.1)–(C1.R) are necessary. Let us assume that (C1. ¯ R) is not satisfied, i.e., for some w

= 0 and w

= 0 (s = ¯ R) we have that rank(E(w))  L

. Without loss of generality we assume that w

= 1. Then we have that E

= E(w) −

r= ¯R

w

E

. Hence, an alternative decomposition of T is T = E(w) ◦ c

+



E

◦ (c

− w

c

),

in which c

− w

c

is not proportional to c

.

We illustrate Theorem 2.4 by means of two examples. Note that the essential uniqueness in Example 2 does not follow from Theorem 2.2, nor from any other theorem in [7].

Example 1. Consider full-rank matrices X, Y ∈ K

and C ∈ K

. Define the rank-2 matrices E

= x

y

+ x

y

and E

= x

y

+ x

y

and define a tensor T ∈ K

by the decomposition in rank-(2,2,1) terms T = E

◦ c

+ E

◦ c

. Obviously, for this decomposition condition (C1) in Theorem 2.4 is not satisfied, since E([1 − 1]) = E

− E

= x

y

− x

y

is also rank-2. The proof of the theorem explains how alternative decompositions can be found. We have for instance T = E([1 − 1]) ◦ c

+ E

◦ (c

+ c

).

Example 2. Consider full-rank matrices X, Y ∈ K

, D

, D

∈ K

, and

C ∈ K

. Define the rank-3 matrices E

= [x

x

x

] · D

· [y

y

y

]

and

E

= [x

x

x

] ·D

·[y

y

y

]

, and define a tensor T ∈ K

by the decomposition

in rank-(3,3,1) terms T = E

◦ c

+ E

◦ c

. It is clear that w

E

+ w

E

has at

least rank 4 if w

= 0 = w

. Since the conditions in Theorem 2.4 are satisfied, the

decomposition of T is essentially unique.

2.3. Computation. An algorithm of the alternating least squares (ALS) type was proposed in [8]. This generalization of the most popular CPD algorithm [11, 19, 20, 27] is sometimes rather slow. Line search algorithms form an interesting alternative. In such schemes the multilinearity of the problem may be exploited to find a good or even the optimal step size [22, 25]. Enhanced line search was used to compute a BTD in [22]. In [23] a Levenberg–Marquardt algorithm was discussed.

A Levenberg–Marquardt iteration step is quite expensive, but the convergence is quadratic. We also mention that the proof of Theorem 2.2 is constructive [7, 26]. It shows that under the conditions of the theorem the decomposition can be computed via a matrix generalized eigenvalue decomposition (GEVD).

3. Blind separation and tensor decomposition. Consider the following data model:

(3.1) Y = M · S + N,

in which Y ∈ K

is the matrix of observed data, M ∈ K

is an unknown mixing matrix, S ∈ K

is a matrix that contains R unknown source signals, and N ∈ K

represents additive noise. The goal of blind signal separation is the estimation of M and/or S, given only Y. N is considered as a perturbation of the equation and will for convenience be ignored in the presentation.

Since the factorization M ·S is not unique, we need to make some assumptions on M and/or S. In principal component analysis (PCA) it is assumed that the columns of M are mutually orthogonal and that also the rows of S are mutually orthogonal [16].

The problem then reduces to the computation of the singular value decomposition (SVD) of Y. In independent component analysis (ICA) it is assumed that the rows of S are mutually statistically independent [4]. Algebraic methods for ICA rely in some way on a CPD of a higher-order tensor derived from X [9]. The latter can be an observed higher-order cumulant tensor or a third-order tensor in which a set of observed covariance matrices is stacked, to mention just two popular alternatives.

In this paper we propose a new technique for blind signal separation, which is based on a different assumption on S. Namely, in section 4 we assume that the sources can be modeled as linear combinations of exponentials. In section 5 we more generally consider exponential polynomials. This structure allows us to formulate the problem in terms of the decomposition of a third-order tensor, derived from Y, in a sum of rank-(L

, L

, 1) terms. If this decomposition is essentially unique, then its computation allows us to blindly separate the signals.

We work as follows. Each row of Y is mapped to an (I × J) Hankel matrix, with I + J − 1 = N . These matrices are stacked in a tensor Y ∈ K

. Formally, we have

(3.2) ( Y)

= (Y)

, 1  i  I, 1  j  J, 1  k  K.

Since the mapping is linear, the (I×J) slices of Y are linear combinations of the Hankel representations of the sources. The linear coefficients correspond to the entries of M.

We have

(3.3) Y =



H

◦ m

,

in which H

∈ K

is the Hankel matrix derived from the rth row of S, 1  r  R.

Under certain conditions on the sources, discussed in sections 4 and 5, the associated

Hankel matrices have low rank. That is, under certain conditions on the sources, decomposition (3.3) is of the type (2.1), and the uniqueness conditions mentioned in section 2.2 apply.

If decomposition (3.3) is essentially unique, then it directly yields the mixing matrix up to column permutation and scaling. Let ˆ M be the estimate of M. We ideally have ˆ M = M · D · P, in which D ∈ K

is a nonsingular diagonal matrix and P ∈ K

a permutation matrix.

If M has full column rank with K  R, then an estimate of the source matrix S can be obtained from (3.1) as ˆ S = ˆ M

· Y. If M is rank-deficient and/or K < R, then this approach is not possible. However, the source values can still be estimated by av- eraging the entries of the estimates ˆ H

of H

along the antidiagonals. Consistent with the indeterminacies of M, we ideally have ˆ S = P

· D

· S. We call the factorization M · S essentially unique when it is only subject to these trivial indeterminacies.

4. Linear combinations of exponentials. In section 4.1 we explain our source model. In section 4.2 we explain how the source structure leads to rank-(L

, L

, 1) terms in (3.3). In section 4.3 we discuss the uniqueness of the decomposition Y = M · S.

4.1. Source model. In this section we assume that the sources can be expressed as linear combinations of exponentials:

(4.1) s

(n)

= s

=

Lr



lr=1

c

z

, 0  n  N − 1, 1  r  R.

Because of Euler’s formula, this model subsumes linear combinations of possibly ex- ponentially damped sinusoids:

e

cos(ωn + φ) = cz

+ c

(z

)

,

where c =

e

and z = e

. Model (4.1) can also include hyperbolic sines and cosines:

sinh(n) = e

− e

2 , cosh(n) = e

+ e

2 .

4.2. Tensor decomposition. If s

(n) can be written as in (4.1), then its asso- ciated Hankel matrix H

∈ K

admits the Vandermonde decomposition

(4.2) H

= V

· diag(c

, c

, . . . , c

) · ˜ V

,

in which the Vandermonde matrices V

∈ K

and ˜ V

∈ K

are defined by (4.3)

V

=

⎡

⎢ ⎢

⎢ ⎣

1 1 · · · 1

z

z

· · · z

.. . .. . .. . z

z

· · · z

⎤

⎥ ⎥

⎥ ⎦ , V ˜

=

⎡

⎢ ⎢

⎢ ⎣

1 1 · · · 1

z

z

· · · z

.. . .. . .. . z

z

· · · z

⎤

⎥ ⎥

⎥ ⎦ ;

see [21]. Let us assume that I, J  max(L

, L

, . . . , L

). Since a Vandermonde matrix

generated by distinct poles has full rank, H

is rank-L

and the rth term in (3.3) is

rank-(L

, L

, 1), 1  r  R. Hence, (3.3) is a decomposition of Y in rank-(L

, L

, 1)

terms.

4.3. Uniqueness.

Theorem 4.1. Consider a matrix M ∈ K

that does not have proportional columns and a matrix S ∈ K

with structure (4.1). Assume that

 

r=1

L

. If all the poles z

are distinct, 1  l

 L

, 1  r  R, then the decomposition Y = M · S is essentially unique.

Proof. The constraint

 

r=1

L

allows us to map the rows of Y to (I × J) Hankel matrices with I, J 

r=1

L

. If all the poles are distinct, then the matrices V = [V

V

· · · V

] · diag(c

, c

, . . . , c

) ∈ K

and V = [ ˜ ˜ V

V ˜

· · · ˜ V

] ∈ K

, defined by (4.3), have full column rank. The result then follows from Theorem 2.2.

Remark 1. Note that under the conditions of Theorem 4.1 exact blind signal separation may be achieved by computing a matrix GEVD.

In the case of repeated poles, essential uniqueness may sometimes by proved by invoking Theorem 2.4. We discuss three examples.

Example 3. Consider a mixture of two sources. The mixing matrix M ∈ K

has full column rank, with K  2. The two rows of S ∈ K

are generated as in (4.1), with poles z

, z

, z

and z

, z

, z

, respectively. Apart from the common pole z

, all poles are distinct. The number of samples N  9. The two sources are associated with rank-3 Hankel matrices H

, H

∈ K

, where I, J  5. Since there is only one common pole, rank(w

H

+ w

H

)  4 when w

= 0 = w

. Theorem 2.4 now implies that decomposition (3.3) is essentially unique. Consequently, the factorization Y = M · S is essentially unique.

Example 4. Consider again a mixture of two sources. The mixing matrix M ∈ K

has full column rank, with K  2. The two rows of S ∈ K

are generated as in (4.1), with poles z

, z

, z

and z

, z

, z

, respectively. The poles z

, z

, z

, and z

are all different. The number of samples N  9. The two sources are associated with rank-3 Hankel matrices H

, H

∈ K

, where I, J  5. Theorem 2.4 indicates that decomposition (3.3) is not essentially unique since, for instance, rank(c

H

− c

H

)  3. As a matter of fact,

Y =

M ·

 c

−c

0 1



·

 c

−c

0 1



· S



is a valid decomposition too.

Example 5. Consider a mixture of three sources. The mixing matrix M ∈ K

has full column rank, with K  3. The three rows of S ∈ K

are generated as in (4.1), with poles z

, z

, z

; z

, z

, z

; z

, z

, z

, respectively. The poles z

, z

, . . . , z

are distinct. The number of samples N  11. The three sources are associated with rank-3 Hankel matrices H

, H

, H

∈ K

, where I, J  6. We first consider linear combinations w

H

+ w

H

+ w

H

in which exactly one of the coefficients is equal to zero. The rank of such linear combinations is at least equal to 4. Next, we consider linear combinations in which none of the coefficients is equal to zero. The rank of such linear combinations is always strictly greater than 3, if c

c

c

+ c

c

c

= 0.

Hence, under the latter condition Theorem 2.4 implies that decomposition (3.3) is essentially unique. The factorization Y = M · S is then essentially unique.

5. Exponential polynomials. In section 5.1 we present a source model that is

more general than the one in 4.1. In section 5.2 we explain how the source structure

leads to rank-(L

, L

, 1) terms in (3.3). In section 5.3 we discuss the uniqueness of

the decomposition Y = M · S.