Blind Multilinear Identification

I Introduction

There are two simple ideas for reducing the complexity or dimension of a problem that are widely applicable because of their simplicity and generality:

Sparsity: resolving a complicated entity, represented by a function $f$, into a sum of a small number of simple or elemental constituents:

Separability: decoupling a complicated entity, represented by a function $g$, that depends on multiple factors into a product of simpler constituents, each depending only on one factor:

The two ideas underlie some of the most useful techniques in engineering and science — Fourier, wavelets, and other orthogonal or sparse representations of signals and images, singular value and eigenvalue decompositions of matrices, separation-of-variables, Fast Fourier Transform, mean field approximation, etc. This article examines the model that combines these two simple ideas:

and we are primarily interested in its inverse problem, i.e., identification of the factors $\varphi_{kp}$ based on noisy measurements of $f$. We shall see that this is a surprisingly effective method for a wide range of identification problems.

Let $\mu_{k}=\max_{p\neq q}\lvert\langle\varphi_{kp},\varphi_{kq}\rangle\rvert$ and define the relative incoherence $\omega_{k}=(1-\mu_{k})/\mu_{k}$ for $k=1,\dots,d$. Note that $\mu_{k}\in$ and $\omega_{k}\in[0,\infty]$. We will show that if $d\geq 3$, and

then the decomposition in (1) is essentially unique and sparsest possible, i.e., $r$ is minimal. Hence we may in principle identify $\varphi_{kp}$ based only on measurements of the mixture $f$.

One of the keys in the identifiability requirement is that $d\geq 3$ or otherwise (when $d=1$ or $2$) the result would not hold. We will show that the condition $d\geq 3$ however leads to a difficulty (that does not happen when $d=1$ or $2$). Since it is rarely, if not never, the case that one has the exact values of $f$, the decomposition (1) is only useful in an idealized scenario. In reality, one has $\hat{f}=f+\varepsilon$, an estimate of $f$ corrupted by noise $\varepsilon$. Solving the inverse problem to (1) would require that we solve a best approximation problem. For example, with the appropriate noise models (see Section V), the best approximation problem often takes the form

with $\lVert\,\cdot\,\rVert$ an $L^{2}$-norm. Now, the trouble is that when $d\geq 3$, this best approximation problem may not have a solution — because the infimum of the loss function is unattainable in general, as we will discuss in Section VIII-A. In view of this, our next result is that when

the infimum in (3) is always attainable, thereby alleviating the aforementioned difficulty. A condition that meets both (2) and (4) follows from the arithmetic-geometric mean inequality

II Sparse separable decompositions

The notion of sparsity dates back to harmonic analysis and approximation theory , and has received a lot of recent attention in compressive sensing . The notion of separability is also classical — the basis behind the separation-of-variables technique in partial differential equations and special functions , fast Fourier transforms on arbitrary groups , mean field approximations in statistical physics , and the naïve Bayes model in machine learning . We describe a simple model that incorporates the two notions.

We will briefly examine the decompositions and approximations of our target function into a sum or integral of separable functions, adopting a tripartite notation for simplicity. There are three cases:

Here, we assume that $\nu$ is some given Borel measure and that $T$ is compact.

This may be viewed as a discretization of the continuous case in the $\mathbf{t}$ variable, i.e., $\theta_{p}(\mathbf{x})=\theta(\mathbf{x},\mathbf{t}_{p})$, $\varphi_{p}(\mathbf{y})=\varphi(\mathbf{y},\mathbf{t}_{p})$, $\psi_{p}(\mathbf{z})=\psi(\mathbf{z},\mathbf{t}_{p})$.

This may be viewed as a further discretization of the semidiscrete case, i.e.$~{}a_{ijk}=f(\mathbf{x}_{i},\mathbf{y}_{j},\mathbf{z}_{k})$, $u_{ip}=\theta_{p}(\mathbf{x}_{i})$, $v_{jp}=\varphi_{p}(\mathbf{y}_{j})$, $w_{kp}=\psi_{p}(\mathbf{z}_{k})$.

It is clear that when $i,j,k$ take finitely many values, the discrete decomposition (7) is always possible with a finite $r$ since the space is of finite dimension. If $i,j,k$ could take infinitely many values, then the finiteness of $r$ requires that equality be replaced by approximation to any arbitrary precision $\varepsilon>0$ in some suitable norm. This follows from the following observation about the semidiscrete decomposition: The space of functions with a semidiscrete representation as in (6), with $r$ finite, is dense in $C^{0}(\Omega)$, the space of continuous functions. This is just a consequence of the Stone-Weierstrass theorem . Discussion of the most general case (5) would require us to go into integral operators, which we will not do as in the present framework we are interested in applications that rely on the inverse problems corresponding to (6) and (7). Nonetheless (5) is expected to be useful and we state it here for completeness. Henceforth, we will simply refer to (6) or (7) as a multilinear decomposition, by which we mean a decomposition into a linear combination of separable functions. We note here that the finite-dimensional discrete version has been studied under several different names — see Section IX. Our emphasis in this paper is the semidiscrete version (6) that applies to multipartite functions on arbitrary domains and are not necessarily finite-dimensional. As such, we will frame most of our discussions in terms of the semidiscrete case, which of course includes the discrete version (7) as a special case (when $\mathbf{x},\mathbf{y},\mathbf{z}$ take only finite discrete values).

Multilinear decompositions arise in many contexts. In machine learning or nonparametric statistics, a fact of note is that Gaussians are separable

under a linear change of coordinates $\mathbf{z}=Q^{\mathsf{T}}\mathbf{x}$ where $A=Q\Lambda Q^{\mathsf{T}}$. Hence Gaussian mixture models of the form

where $A_{i}A_{j}=A_{j}A_{i}$ for all $i\neq j$ (and therefore $A_{1},\dots,A_{m}$ have a common eigenbasis) may likewise be transformed with a suitable linear change of coordinates into a multilinear decomposition as in (6).

We will later see several more examples from signal processing, telecommunications, and spectroscopy.

The multilinear decomposition — an additive decomposition into multiplicatively decomposable components — is extremely simple but models a wide range of phenomena in signal processing and spectroscopy. The main message of this article is that the corresponding inverse problem — recovering the factors $\theta_{p},\varphi_{p},\psi_{p}$ from noisy measurements of $f$ — can be solved under mild assumptions and yields a class of techniques for a range of applications (cf. Section IX) that we shall collectively call multilinear identification. We wish to highlight in particular that multilinear identification gives a deterministic approach for solving the problem of joint localization and estimation of radiating sources with short data lengths. This is superior to previous cumulants-based approaches , which require (i) longer data lengths; and (ii) statistically independent sources.

The experienced reader would probably guess that such a powerful technique must be fraught with difficulties and he would be right. The inverse problem to (6), like most other inverse problems, faces issues of existence, uniqueness, and computability. The approximation problem involved can be ill-posed in the worst possible way (cf. Section III). Fortunately, in part prompted by recent work in compressed sensing and matrix completion ), we show that mild assumptions on coherence allows one to overcome most of these difficulties (cf. Section VIII).

III Finite rank multipartite functions

In this section, we will discuss the notion of rank, which measures the sparsity of a multilinear decomposition, and the notion of Kruskal rank, which measures the uniqueness of a multilinear decomposition in a somewhat more restrictive sense. Why is uniqueness important? It can be answered in one word: Identifiability. More specifically, a unique decomposition means that we may in principle identify the factors. To be completely precise, we will first need to define the terms in the previous sentence, namely, ‘unique’, ‘decomposition’, and ‘factor’. Before we do that, we will introduce the tensor product notation. It is not necessary to know anything about tensor product of Hilbert spaces to follow what we present below. We shall assume that all our Hilbert spaces are separable and so there is no loss of generality in assuming at the outset that they are just $L^{2}(X)$ for some $\sigma$-finite $X$.

Let $X_{1},\dots,X_{d}$ be $\sigma$-finite measurable spaces. There is a natural Hilbert space isomorphism

with $\varphi_{kp}\in L^{2}(X_{k})$. The tensor product of functions $\varphi_{1}\in L^{2}(X_{1}),\dots,\varphi_{d}\in L^{2}(X_{d})$ is denoted by $\varphi_{1}\otimes\dots\otimes\varphi_{d}$ and is the function in $L^{2}(X_{1}\times\dots\times X_{d})$ defined by

With this notation, we may rewrite (9) as

“Multipartite functions are infinite-dimensional tensors.”

In this paper, functions having a finite decomposition will play a central role; for these we define

provided $f\neq 0$. The zero function is defined to have rank and we say $\operatorname{rank}(f)=\infty$ if such a decomposition is not possible.

We will call a function $f$ with $\operatorname{rank}(f)\leq r$ a rank-$r$ function. Such a function may be written as a sum of $r$ separable functions but possibly fewer. A decomposition of the form

will be called a rank-$r$ multilinear decomposition. Note that the qualificative ‘rank-$r$’ will always mean ‘rank not more than $r$’. If we wish to refer to a function $f$ with rank exactly $r$, we will just specify that $\operatorname*{rank}(f)=r$. In this case, the rank-$r$ multilinear decomposition in (11) is of mininum length and we call it a rank-retaining multilinear decomposition of $f$.

A rank-$1$ function is both non-zero and decomposable, i.e., of the form $\varphi_{1}\otimes\dots\otimes\varphi_{d}$ where $\varphi_{k}\in L^{2}(X_{k})$. This agrees precisely with the notion of a separable function. Observe that the inner product (and therefore the norm) on $L^{2}(X_{1}\times\dots\times X_{d})$ of a rank-$1$ function splits into a product

where $\langle\cdot,\cdot\rangle_{p}$ denotes the inner product of $L^{2}(X_{p})$. This inner product extends linearly to finite-rank elements of $L^{2}(X_{1}\times\dots\times X_{d})$: for $f=\sum_{p=1}^{r}\varphi_{1p}\otimes\dots\otimes\varphi_{dp}$ and $g=\sum_{q=1}^{s}\psi_{1q}\otimes\dots\otimes\psi_{dq}$, we have

In fact this is how a tensor product of Hilbert spaces (the right hand side of (8)) is usually defined, namely, as the completion of the set of finite-rank elements of $L^{2}(X_{1}\times\dots\times X_{d})$ under this inner product.

When $X_{1},\dots,X_{d}$ are finite sets, then all functions in $L^{2}(X_{1}\times\dots\times X_{d})$ are of finite rank (and may in fact be viewed as hypermatrices or tensors as discussed in Section II). Otherwise there will be functions in $L^{2}(X_{1}\times\dots\times X_{d})$ of infinite rank. However, since we have assumed that $X_{1},\dots,X_{d}$ are $\sigma$-finite measurable spaces, the set of all finite-rank $f$ will always be dense in $L^{2}(X_{1}\times\dots\times X_{d})$ by the Stone-Weierstrass theorem.

The next statement is a straightforward observation about multilinear decompositions of finite-rank functions but since it is central to this article we state it as a theorem. It is also tempting to call the decomposition a ‘singular value decomposition’ given its similarities with the usual matrix singular value decomposition (cf. Example 4).

Let $f\in L^{2}(X_{1}\times\dots\times X_{d})$ be of finite rank. Then there exists a rank-$r$ multilinear decomposition

the functions $\varphi_{kp}\in L^{2}(X_{p})$ are of unit norm,

the coefficients $\sigma_{1},\dots,\sigma_{r}$ are real positive, and

This requires nothing more than rewriting the sum in (11) as a linear combination with the positive $\sigma_{p}$’s accounting for the norms of the summands and then re-indexing them in descending order of magnitudes. ∎

While the usual singular value decomposition of a matrix would also have properties (14), (15), and (16), the one crucial difference here is that our ‘singular vectors’ $\varphi_{k1},\dots,\varphi_{kr}$ in (13) will only be of unit norms but will not in general be orthonormal. Given this, we will not expect properties like the Eckhart-Young theorem, or that $\sigma_{1}^{2}+\dots+\sigma_{r}^{2}=\lVert f\rVert^{2}$, etc, to hold for (13) (cf. Section VI for more details).

One may think of the multilinear decomposition (13) as being similar in spirit, although not in substance, to Kolmogorov’s superposition principle ; the main message of which is that:

“There are no true multivariate functions.”

More precisely, the principle states that continuous functions in multiple variables can be expressed as a composition of a univariate function with other univariate functions. For readers not familiar with this remarkable result, we state here a version of it due to Kahane

It is in general not easy to determine $g$ and $\varphi_{1},\dots,\varphi_{2d+1}$ given such a function $f$. A multilinear decomposition of the form (13) alleviates this by allowing $g$ to be the simplest multivariate function, namely, the product function,

At this stage, it would be instructive to give a few examples for concreteness.

which clearly is the same as (9). The singular value decomposition (svd) of $A$ yields one such decomposition, where $\{\mathbf{u}_{1},\dots,\mathbf{u}_{r}\}$ and $\{\mathbf{v}_{1},\dots,\mathbf{v}_{r}\}$ are both orthonormal. But in general a rank-retaining decomposition of the form (13) will not have such a property.

where $\psi_{p}^{\ast}$ denotes the dual form of $\psi_{p}$.

Examples 4 and 5 are well-known but they are bipartite examples, i.e. $d=2$ in (13). This article is primarily concerned with the $d$-partite case where $d\geq 3$, which has received far less attention. As we have alluded to in the previous section, the identification techniques in this article will rely crucially on the fact that $d\geq 3$.

and if $A=0$, then its rank is set to be . This agrees of course with our use of the word rank in (10), the only difference is notational, since (20) may be written in the form

IV Uniqueness of multilinear decompositions

In Theorem 2, we chose the coefficients to be in descending order of magnitude and require the factors in each separable term to be of unit norm. This is largely to ensure as much uniqueness in the multilinear decomposition as generally possible. However there remain two obvious ways to obtain trivially different multilinear decompositions: (i) one may scale the factors $\varphi_{1p},\dots,\varphi_{dp}$ by arbitrary unimodulus complex numbers as long as their product is $1$; (ii) when two or more successive coefficients are equal, their orders in the sum may be arbitrarily permuted. We will call a multilinear decomposition of $f$ that meets the conditions in Theorem 2 essentially unique if the only other such decompositions of $f$ differ in one or both of these manners.

It is perhaps astonishing that when $d>2$, a sufficient condition for essential uniqueness can be derived with relatively mild conditions on the factors. This relies on the notion of Kruskal rank, which we will now define.

We now generalize Kruskal’s famous result to tensor products of arbitrary Hilbert spaces, possibly of infinite dimensions. But first let us be more specific about essential uniqueness.

We shall say that a multilinear decomposition of the form (13) (satisfying both (16) and (15)) is essentially unique if given another such decomposition,

we must have (i) the coefficients $\sigma_{p}=\lambda_{p}$ for all $p=1,\dots,r$; and (ii) the factors $\varphi_{1p},\dots,\varphi_{dp}$ and $\psi_{1p},\dots,\psi_{dp}$ differ at most via unimodulus scaling, i.e.

where $\theta_{1p}+\dots+\theta_{dp}\equiv 0\operatorname{mod}2\pi$, for all $p=1,\dots,r$. In the event when successive coefficients are equal, $\sigma_{p-1}>\sigma_{p}=\sigma_{p+1}=\dots=\sigma_{p+q}>\sigma_{p+q+1}$, the uniqueness of the factors in (ii) is only up to relabelling of indices, i.e.$\ p,\dots,p+q$.

Let $f\in L^{2}(X_{1}\times\dots\times X_{d})$ be of finite rank. Then a multilinear decomposition of the form

is both essentially unique and rank-retaining, i.e., $r=\operatorname{rank}f$, if the following condition is satisfied:

where $\Phi_{k}=\{\varphi_{k1},\dots,\varphi_{kr}\}$ for $k=1,\dots,d$.

It follows immediately why we usually need $d\geq 3$ for identifiability.

A necessary condition for Kruskal’s inequality (25) to hold is that $d\geq 3$.

If $d=2$, then $2r+d-1=2r+1>\operatorname{krank}\Phi_{1}+\operatorname{krank}\Phi_{2}$ since the Kruskal rank of of $r$ vectors cannot exceed $r$. Likewise for $d=1$. ∎

Lemma 9 shows that the condition in (25) is sufficient to ensure uniqueness and it is known that the condition is not necessary. In an appropriate sense, the condition is sharp . We should note that the version of Lemma 9 that we state here for general $d\geq 3$ is due to Sidiropoulos and Bro . Kruskal’s original version is only for $d=3$.

The main problem with Lemma 9 is that the condition (25) is difficult to check since the right-hand side cannot be readily computed. See Section VIII-F for a discussion.

the expected rank of $L^{2}(X_{1}\times\dots\times X_{d})$, since it is heuristically the minimum $r$ expected for a multilinear decomposition (13).

Let the notations be as above. If $f\in L^{2}(X_{1}\times\dots\times X_{d})$ has rank smaller than the expected rank, i.e.

then $f$ admits at most a finite number of distinct rank-retaining decompositions.

This proposition has been proved in several cases, including symmetric tensors , but the proof still remains incomplete for tensors of most general form .

V Estimation of multilinear decompositions

In practice we would only have at our disposal $\hat{f}$, a measurement of $f$ corrupted by noise. Recall that our model for $f$ takes the form

Then we would often have to solve an approximation problem corresponding to (26) of the form

which we will call a best rank-$r$ approximation problem. A solution to (27), if exists, will be called a best rank-$r$ approximation of $\hat{f}$.

We will give some motivations as to why such an approximation is reasonable. Assuming that the norm in (27) is the $L^{2}$-norm and that the factors $\varphi_{kp}$, $p=1,\dots r$ and $k=1,\dots d$, have been determined in advance and we are just trying to estimate the parameters $\alpha_{1},\dots,\alpha_{r}$ from $\hat{f}^{(1)},\dots,$ $\hat{f}^{(N)}$ a finite sample of size $N$ of measurements of $f$ corrupted by noise, then the solution of the approximation problem in (27) is in fact (i) a maximum likelihood estimator (mle) if the noise is zero mean Gaussian, and (ii) a best linear unbiased estimator (blue) if the noise has zero mean and finite variance. Of course in our identification problems, the factors $\varphi_{kp}$’s are not known and have to be estimated too. A probabilistic model in this situation would take us too far afield. Note that even for the case $d=2$ and where the domain of $f$ $X_{1}\times X_{2}$ is a finite set, a case that essentially reduces to principal components analysis (pca), a probabilistic model along the lines of requires several strong assumptions and was only developed as late as 1999. The lack of a formal probabilistic model has not stopped pca, proposed in 1901 , to be an invaluable tool in the intervening century.

VI Existence of best multilinear approximation

As we mentioned in the previous section, in realistic situation where measurements are corrupted by additive noise, one has to extract the factors $\varphi_{kp}$’s and $\alpha_{p}$ through solving an approximation problem (27), that we now write in a slightly different (but equivalent) form,

Note that by Theorem 2, we may assume that the coefficients $\boldsymbol{\alpha}=(\alpha_{1},\dots,\alpha_{r})$ are real and nonnegative valued without any loss of generality. Such a form is also natural in applications given that $\alpha_{p}$ usually captures the magnitude of whatever quantity that is represented by the $p$ summand.

We will see this problem, whether in the form (27) or (28), has no solution in general. We will first observe a somewhat unusual phenomenon in multilinear decomposition of $d$-partite functions where $d\geq 3$, namely, a sequence of rank-$r$ functions (each with an rank-$r$ multilinear decomposition) can converge to a limit that is not rank-$r$ (has no rank-$r$ multilinear decomposition).

In Example 12, $\operatorname*{rank}(\hat{f})=3$ iff $\varphi_{i},\psi_{i}$ are linearly independent, $i=1,2,3$. Furthermore, it is clear that $\operatorname*{rank}(f_{n})\leq 2$ and

Note that our fundamental approximation problem may be regarded as the approximation problem

which always exists for an $f$ with $\operatorname*{rank}(f)\leq r$. The discussion above shows that there are target functions $\hat{f}$ for which

and thus (28) or (30) does not need to have a solution in general. This is such a crucial point that we are obliged to formally state it.

For $d\geq 3$, the best approximation of a $d$-partite function by a sum of $p$ products of $d$ separable functions does not exist in general.

Take the tripartite function $\hat{f}\in L^{2}(X_{1}\times X_{2}\times X_{3})$ in Example 12. Suppose we seek a best rank-$2$ approximation, in other words, we seek to solve the minimization problem

so as make $\|\hat{f}-\gamma g_{1}\otimes g_{2}\otimes g_{3}-\eta h_{1}\otimes h_{2}\otimes h_{3}\|$ as small as we desired by virtue of (29). However there is no rank-$2$ function $\gamma g_{1}\otimes g_{2}\otimes g_{3}+\eta h_{1}\otimes h_{2}\otimes h_{3}$ for which

In other words, the zero infimum can never be attained. ∎

Our construction above is based on an earlier construction in . The first such example was given in , which also contains the very first definition of border rank. We will define it here for $d$-partite functions. When $X_{1},\dots,X_{d}$ are finite sets, this reduces to the original definition in for hypermatrices.

Let $f\in L^{2}(X_{1}\times\dots\times X_{d})$. The border rank of $f$ is defined as

We say $\overline{\operatorname*{rank}}(f)=\infty$ if such a finite $r$ does not exist.

The discussions above show that strict inequality can occur. In fact, for the $\hat{f}$ in Example 12, $\overline{\operatorname*{rank}}(\hat{f})=2$ while $\operatorname*{rank}(\hat{f})=3$.

We would like to mention here that this problem applies to operators too. Optimal approximation of an operator by a sum of tensor/Kronecker products of lower-dimensional operators, which arises in numerical operator calculus , is in general an ill-posed problem whose existence cannot be guaranteed.

For linearly independent operators $\Phi_{i},\Psi_{i}:V_{i}\rightarrow W_{i}$, $i=1,2,3$, let $\widehat{T}:V_{1}\otimes V_{2}\otimes V_{3}\rightarrow W_{1}\otimes W_{2}\otimes W_{3}$ be

An example of an operator that has the form in (31) is the $3m$-dimensional Laplacian $\Delta_{3m}$, which can be expressed in terms of the $m$-dimensional Laplacian $\Delta_{m}$ as

There are several simple but artificial ways to alleviate the issue of non-existent best approximant. Observe from the proof of Theorem 14 that the coefficients in the approximant $\gamma,\eta$ becomes unbounded in the limit. Likewise we see this happening in Example 16. In fact this must always happen — in the event when a function or operator is approximated by a rank-$r$ function, i.e.

and if a best approximation does not exist, then the $r$ coefficients $\alpha_{1},\dots,\alpha_{r}$ must all diverge in magnitude to $+\infty$ as the approximant converges to the infimum of the norm loss function in (32). This result was first established in [21, Proposition 4.9].

So a simple but artificial way to prevent the nonexistence issue is to simply limit the sizes of the coefficients $\alpha_{1},\dots,\alpha_{r}$ in the approximant. One way to achieve this is regularization — adding a regularization term to our objective function in (28) to penalize large coefficients. A common choice is Tychonoff regularization, which uses a sum-of-squares regularization term:

Here, $\lambda$ is an arbitrarily chosen regularization parameter. It can be seen that this is equivalent to constraining the sizes $\alpha_{1},\dots,\alpha_{r}$ to $\sum_{p=1}^{r}\lvert\alpha_{p}\rvert^{2}=\rho$, with $\rho$ being determined a posteriori from $\lambda$. The main drawback of such constraints is that $\rho$ and $\lambda$ are arbitrary, and that they generally have no physical meaning.

More generally, one may alleviate the nonexistence issue by restricting the optimization problem (30) to a compact subset of its non-compact feasible set

Limiting the sizes of $\alpha_{1},\dots,\alpha_{r}$ is a special case but there are several other simple (but also artificial) strategies. In , the factors $\varphi_{k1},\dots,\varphi_{kp}$ are required to be orthogonal for all $k\in\{1,\dots,d\}$, i.e.

This remedy is acceptable only in very restrictive conditions. In fact a necessary condition for this to work is that

It is also trivial to see that imposing orthogonality between the separable factors removes this problem

This constraint is slightly less restrictive — by (12), it is equivalent to requiring (34) for some $k\in\{1,\dots,d\}$. Both (34) and (35) are nonetheless so restrictive as to exclude the most useful circumstances for the model (13), which usually involves factors that have no reason to be orthogonal, as we will see in Section IX. In fact, Kruskal’s uniqueness condition is such a potent tool precisely because it does not require orthogonality.

The conditions (34), (35), and (33) all limit the feasible sets for the original approximation (28) to a much smaller compact subset of the original feasible set. This is not the case for nonnegative constraints. In it was shown that the following best rank-$r$ approximation problem for a nonnegative-valued $\hat{f}$ and where the coefficients $\alpha_{p}$ and factors $\varphi_{kp}$ of the approximants are also nonnegative valued, i.e.

always has a solution. The feasible set in this case is non-compact and has nonempty interior within the feasible set of our original problem (28). The nonnegativity constraints are natural in some applications, such as the fluorescence spectroscopy one described in Section IX-F, where $\varphi_{kp}$ represent intensities and concentrations, and are therefore nonnegative valued.

where $Y=\frac{1}{m}\sum_{i=1}^{m}\mathbf{y}_{i}\mathbf{y}_{i}^{\mathsf{T}}$. However the problem will not have a solution when the number of samples is smaller than the dimension, i.e., $m<n$, as the infimum of the loss function in (36) cannot be attained by any $X$ in the feasible set. This is an indication that we should seek more samples (so that we could get $m\geq n$, which will guarantee the attainment of the infimum) or use a different model (e.g., determine if $X^{-1}$ might have some a priori zero entries due to statistical independence of the variables). It is usually unwise to impose artificial constraints on the covariance matrix $X$ just so that the loss function in (36) would attain an infimum on a smaller feasible set — the thereby obtained ‘solution’ may bear no relation to the true solution that we want.

Our goal in Section VIII-A is to define a type of physically meaningful constraints via the notion of coherence. It ensures the existence of a unique minimum, but not via an artificial limitation of the optimization problem to a convenient subset of the feasible set. In the applications we discuss in Section IX, we will see that it is natural to expect existence of a solution when coherence is small enough, but not otherwise. So when our model is ill-posed or ill-conditioned, we are warned by the size of the coherence and could seek other remedies (collect more measurements, use a different model, etc) instead of forcing a ‘solution’ that bears no relation to reality. But before we get to that we will examine why, unlike in compressed sensing and matrix completion, the approximation of rank by a ratio of nuclear and spectral norms could not be expected to work here.

VII Nuclear and spectral norms

We introduce the notion of nuclear and spectral norms for multipartite functions. Our main purpose is to see if they could be used to alleviate the problem discussed in Section VI, namely, that a $d$-partite function may not have a best approximation by a sum of $r$ separable functions.

The definition of nuclear norm follows naturally from the definition of rank in Section III.

We define the nuclear norm (or Schatten $1$-norm) of $f\in L^{2}(X_{1}\times\dots\times X_{d})$ as

Note that for rank-$1$ functions, we always have that

A finite rank function always has finite nuclear norm but in general a function in $L^{2}(X_{1}\times\dots\times X_{d})$ need not have finite nuclear norm.

The definition of the spectral norm of a multipartite function is motivated by the fact the usual spectral norm of a matrix $A$ equals the maximal absolute value of its inner product $\operatorname{tr}(AX)$ with rank-$1$ unit-norm matrices $X=\mathbf{u}\mathbf{v}^{\ast}$, $\lVert\mathbf{u}\rVert_{2}=\lVert\mathbf{v}\rVert_{2}=1$.

We define the spectral norm of $f\in L^{2}(X_{1}\times\dots\times X_{d})$ as

Here we write $\lVert\,\cdot\,\rVert$ for the $L^{2}$-norms in $L^{2}(X_{k})$, $k=1,\dots,d$. Alternatively, we may also use $\operatorname{Re}\langle f,\varphi_{1}\otimes\dots\otimes\varphi_{d}\rangle$ in place of $\lvert\langle f,\varphi_{1}\otimes\dots\otimes\varphi_{d}\rangle\rvert$ in (39), which does not change its value. Note that a function in $L^{2}(X_{1}\times\dots\times X_{d})$ always has finite spectral norm.

We shall write $\lVert\,\cdot\,\rVert=\lVert\,\cdot\,\rVert_{2}$. The spectral norm of $T$ is

We have regarded $T$ as a trilinear functional defined by $T(\mathbf{x},\mathbf{y},\mathbf{z})=\sum_{i,j,k=1}^{n_{1},n_{2},n_{3}}t_{ijk}x_{i}y_{j}z_{k}$ and $\lVert T\rVert_{2,2,2}$ is its induced norm as defined in . Again, these clearly extend to any $d$-tensors. We will see in Lemma 21 that the nuclear and spectral norms for tensors are dual to each other.

Note that we have used the term tensors, as opposed to hypermatrices, in the above example. In fact, Definition 17 defines nuclear norms for the tensors, not just their coordinate representations as hypermatrices (see our discussion after Example 6), because of the following invariant properties.

has the property that if $T$ has a multilinear decomposition of the form

One may similarly show that the spectral norm is also unitarily invariant or deduce the fact from Lemma 21 below. ∎

If $\lVert\,\cdot\,\rVert$ is the norm induced by the inner product, then $\lVert\,\cdot\,\rVert^{\vee}=\lVert\,\cdot\,\rVert$; but in general they are different. Nevertheless one always have that $(\lVert\,\cdot\,\rVert^{\vee})^{\vee}=\lVert\,\cdot\,\rVert$ and

Let $X_{1},\dots,X_{d}$ be finite sets. Then nuclear and spectral norms on $L^{2}(X_{1}\times\dots\times X_{d})$ satisfy

We first need to establish (42) without invoking (41). Since $X_{1},\dots,X_{d}$ are finite, any $g\in L^{2}(X_{1}\times\dots\times X_{d})$ is of finite rank. Take any multilinear decomposition

by definition of spectral norm. Taking infimum over all finite-rank decompositions, we arrive at (42) by definition of nuclear norm. Hence

On the other hand, using (41) for $\lVert\,\cdot\,\rVert_{\ast}$ and $\lVert\,\cdot\,\rVert_{\ast}^{\vee}$, we get

where the last equality follows from (38). ∎

When $X_{1},\dots,X_{d}$ are only required to be $\sigma$-finite measurable spaces, we may use a limiting argument to show that (42) still holds for $f$ of finite spectral norm and $g$ of finite nuclear norm; a proper generalization of (43) is more subtle and we will leave this to future work since we do not require it in this article.

We had suspected that the following generalization might perhaps be true, namely, rank, nuclear, and spectral norms as defined in (10), (37), and (39) would also satisfy the same inequality:

If true, this would immediately imply the same for border rank

by a limiting argument. Furthermore, (44) would provide a simple way to remedy the nonexistence problem highlighted in Theorem 14: One may use the ratio $\lVert f\rVert_{\ast}/\lVert f\rVert_{\sigma}$ as a ‘proxy’ in place of $\operatorname*{rank}(f)$ and replace the condition $\operatorname*{rank}(f)\leq r$ by the weaker condition $\lVert f\rVert_{\ast}\leq r\lVert f\rVert_{\sigma}$. The discussion in Section VI shows that there are $\hat{f}$ for which

Unfortunately, (44) is not true when $d>2$. The following example shows that nuclear norm is not an underestimator of rank on the spectral norm unit ball, and can in fact be arbitrarily larger than rank on the spectral norm unit ball.

for some $c>0$ and for $n$ sufficiently large. On the other hand, Derksen has recently established the exact values for the nuclear and spectral norm of $T_{n}$:

Fortunately, we do not need to rely on (45) for the applications consider in this article. Instead another workaround that uses the notion of coherence, discussed in the next section, is more naturally applicable in our situtations.

VIII Coherence

We will show in this section that a simple measure of angular constraint called coherence, or rather, the closely related notion of relative incoherence, allows us to alleviate two problems simultaneously: the computational intractability of checking for uniqueness discussed in Section IV and the non-existence of a best approximant in Section VI.

where the supremum is taken over all distinct pairs $\varphi,\psi\in\Phi$. If $\Phi=\{\varphi_{1},\dots,\varphi_{r}\}$ is finite, we also write $\mu(\varphi_{1},\dots,\varphi_{r}):=\max_{p\neq q}\lvert\langle\varphi_{p},\varphi_{q}\rangle\rvert$.

We adopt the convention that whenever we write $\mu(\Phi)$ (resp.$\ \mu(\varphi_{1},\dots,\varphi_{r})$) as in Definition 23, it is implicitly implied that all elements of $\Phi$ (resp.$\ \varphi_{1},\dots,\varphi_{r}$) are of unit norm.

The notion of coherence has received different names in the literature: mutual incoherence of two dictionaries , mutual coherence of two dictionaries , the coherence of a subspace projection , etc. The version here follows that of . Usually, dictionaries are finite or countable, but we have here a continuum of atoms. Clearly, $0\leq\mu(\Phi)\leq 1$, and $\mu(\Phi)=0$ iff $\varphi_{1},\dots,\varphi_{r}$ are orthonormal. Also, $\mu(\Phi)=1$ iff $\Phi$ contains at least a pair of collinear elements, i.e., $\varphi_{p}=\lambda\varphi_{q}$ for some $p\neq q$, $\lambda\neq 0$.

We find it useful to introduce a closely related notion that we call relative incoherence. It allows us to formulate some of our results slightly more elegantly.

For a finite set of unit vectors $\Phi=\{\varphi_{1},\dots,\varphi_{r}\}$, we will also write $\omega(\varphi_{1},\dots,\varphi_{r})$ occasionally.

It follows from our observation about coherence that $0\leq\omega(\Phi)\leq\infty$, $\omega(\Phi)=\infty$ iff $\varphi_{1},\dots,\varphi_{r}$ are orthonormal, and $\omega(\Phi)=0$ iff $\Phi$ contains at least a pair of collinear elements.

In the next few subsections, we will see respectively how coherence can inform us about the existence (Section VIII-A), uniqueness (Section VIII-B), as well as both existence and uniqueness (Section VIII-C) of a solution to the best rank-$r$ multilinear approximation problem (28). We will also see how it can be used for establishing exact recoverability (Section VIII-D) and approximation bounds (Section VIII-E) in greedy algorithms.

The goal is to prevent the phenomenon we observed in Example 12 to occur, by imposing natural and weak constraints; we do not want to reduce the search to a compact set. It is clear that the objective is not coercive, which explains why the minimum may not exist. But with an additional condition on the coherence, we shall be able to prove existence thanks to coercivity.

The following shows that a solution to the bounded coherence best rank-$r$ approximation problem always exists:

Let $f\in L^{2}(X_{1}\times\dots\times X_{d})$ be a $d$-partite function. If

where $\mu_{k}$ denotes the coherence as in Definition 23, then

The equivalence between (46) and (47) follows from Definition 24. We show that if either of these conditions are met, then the loss function is coercive. We have the following inequalities

which is true because $\sum_{p\neq q}(\lvert\lambda_{p}\rvert-\lvert\lambda_{q}\rvert)^{2}\geq 0$. This yields

Since by assumption $(r-1)\prod_{k=1}^{d}\mu_{k}<1$, it is clear that the left hand side of (49) tends to infinity as $\|\boldsymbol{\lambda}\|_{2}\rightarrow\infty$. And because $f$ is fixed, $\left\|f-\sum_{p=1}^{r}\lambda_{p}\varphi_{1p}\otimes\dots\otimes\varphi_{dp}\right\|$ also tends to infinity as $\|\boldsymbol{\lambda}\|_{2}\rightarrow\infty$. This proves coercivity of the loss function and hence the existential statement. ∎

The condition (46) or, equivalently, (47), in Theorem 25 is sharp in an appropriate sense. Theorem 25 shows that the condition (47) is sufficient in the sense that it guarantees a best rank-$r$ approximation when the condition is met. We show that it is also necessary in the sense that if (47) does not hold, then there are examples where a best rank-$r$ approximation fails to exist.

In fact, let $\hat{f}$ be as in Example 12. As demonstrated in the proof of Theorem 14, the infimum for the case $d=3$ and $r=2$,

for $k=1,2,3$, the corresponding coherence

as $n\rightarrow\infty$. For any values of $\mu_{1},\mu_{2},\mu_{3}\in$ such that (47) holds, i.e.$\ \mu_{1}\mu_{2}\mu_{3}<1/(r-1)=1$, we cannot possibly have $\mu(g_{k},h_{k})\leq\mu_{k}$ for all $k=1,2,3$ since

VIII-B Uniqueness and minimality via coherence

In order to relate uniqueness and minimality of multilinear decompositions to coherence, we need a simple observation about the notion of Kruskal rank introduced in Definition 7.

Let $\Phi\subseteq L^{2}(X_{1}\times\dots\times X_{d})$ be finite and $\operatorname{krank}\Phi<\dim\operatorname{span}\Phi$. Then

Let $s=\operatorname{krank}{\Phi}+1$. Then there exists a subset of $s$ distinct unit vectors in $\Phi$, $\{\varphi_{1},\dots,\varphi_{s}\}$ such that $\alpha_{1}\varphi_{1}+\dots+\alpha_{s}\varphi_{s}=0$ with $\lvert\alpha_{1}\rvert=\max\{\lvert\alpha_{1}\rvert,\dots,\lvert\alpha_{s}\rvert\}>0$. Taking inner product with $\varphi_{1}$ we get $\alpha_{1}=-\alpha_{2}\langle\varphi_{2},\varphi_{1}\rangle-\dots-\alpha_{s}\langle\varphi_{s},\varphi_{1}\rangle$ and so $\lvert\alpha_{1}\rvert\leq(\lvert\alpha_{2}\rvert+\dots+\lvert\alpha_{s}\rvert)\mu(\Phi)$. Dividing by $\lvert\alpha_{1}\rvert$ then yields $1\leq(s-1)\mu(\Phi)$. The condition $\operatorname{krank}{\Phi}<\dim\operatorname{span}\Phi$ prevents $\Phi$ from being orthonormal, so $\mu(\Phi)>0$ and we obtain (50). ∎

We now characterize the uniqueness of the rank-retaining decomposition in terms of coherence introduced in Definition 23.

Suppose $f\in L^{2}(X_{1}\times\dots\times X_{d})$ has a multilinear decomposition

where $\Phi_{k}:=\{\varphi_{k1},\dots,\varphi_{kr}\}$ are elements in $L^{2}(X_{k})$ of unit norm and $\operatorname{krank}\Phi_{k}<\dim\operatorname{span}\Phi_{k}$ for all $k=1,\dots,d$. Let $\omega_{k}=\omega(\Phi_{k})$. If

then $r=\operatorname{rank}(f)$ and the decomposition is essentially unique. In terms of coherence, (51) takes the form

Inequality (52) implies that $\sum_{k=1}^{d}\mu_{k}^{-1}\geq 2r+d-1$, where $\mu_{k}$ denotes $\mu(\Phi_{k})$. If it is satisfied, then so is Kruskal’s condition (25) thanks to Lemma 26. The result hence directly follows from Lemma 9 and Definition 24. ∎

Note that unlike the Kruskal ranks in (25), the coherences in (52) are trivial to compute. In addition to uniqueness, an easy but important consequence of Theorem 27 is that it provides a readily checkable sufficient condition for tensor rank, which is NP-hard over any field .

Since the purpose of Theorem 27 is to provide a computationally feasible alternative of Lemma 9, excluding the case $\operatorname{krank}\Phi_{k}=\dim\operatorname{span}\Phi_{k}$ is not an issue. Note that $\operatorname{krank}\Phi_{k}=\dim\operatorname{span}\Phi_{k}$ iff $\Phi_{k}$ comprises linearly independent elements, and the latter can be checked in polynomial time. So this is a case where Lemma 9 can be readily checked and one does not need Theorem 27.

VIII-C Existence and uniqueness via coherence

The following existence and uniqueness sufficient condition may now be deduced from Theorems 25 and 27.

If $d\geq 3$ and if coherences $\mu_{k}$ satisfy

then the bounded coherence best rank-$r$ approximation problem has a unique solution up to unimodulus scaling.

The existence in the case $r=1$ is assured, because the set of separable functions $\{\varphi_{1}\otimes\dots\otimes\varphi_{d}:\varphi_{k}\in L^{2}(X_{k})\}$ is closed. Consider thus the case $r\geq 2$. Since the function $f(x)=\frac{1}{x}-\left(\frac{d}{2x+d-1}\right)^{d}$ is strictly positive for $x\geq 2$ and $d\geq 3$, condition (53) implies that $\prod_{k=1}^{d}\mu_{k}$ is smaller than $1/r$, which permits to claim that the solution exists by calling for Theorem 25. Next in order to prove uniqueness, we use the inequality between harmonic and geometric means: if (53) is verified, then we also necessarily have $d\left(\sum_{k=1}^{d}\mu_{k}^{-1}\right)^{-1}\leq\frac{d}{2r+d-1}$. Hence $\sum_{k=1}^{d}\mu_{k}^{-1}\geq 2r+d-1$ and we can apply Theorem 27. ∎

In practice, simpler expressions than (53) can be more attractive for computational purposes. These can be derived from the inequalities between means:

Examples of stronger sufficient conditions that could be used in place of (53) include

VIII-D Exact recoverability via coherence

We now describe a result that follows from the remarkable work of Temlyakov. It allows us to in principle determine the multilinear decomposition meeting the type of coherence conditions in Section VIII-A.

for some $\mu\in[0,1)$ to be chosen later. Recall that the elements of $\Phi$ are implicitly assumed to be of unit norm (cf. remark after Definition 23).

$h_{m}\in L^{2}(X_{1}\times\dots\times X_{d})$ is a projection of $f$ onto $\operatorname{span}(g_{1},\dots,g_{m})$, i.e.

$f_{m}\in L^{2}(X_{1}\times\dots\times X_{d})$ is a deflation of $f$ by $h_{m}$, i.e.

Note that deflation alone, without the coherence requirement, generally does not work for computing multilinear decompositions . The following result, adapted here for our purpose, was proved for any arbitrary dictionary in .

Suppose $f\in L^{2}(X_{1}\times\dots\times X_{d})$ has a multilinear decomposition

with $\varphi_{1p}\otimes\dots\otimes\varphi_{dp}\in\Phi$ and the condition that

for some $t\in(0,1]$. Then the woga algorithm recovers the factors exactly, or more precisely, $f_{r}=0$ and thus $f=h_{r}$.

So $h_{r}$, by its definition in (57) and our choice of $\Phi$, is given in the form of a linear combination of rank-$1$ functions, i.e., an rank-$r$ multilinear decomposition.

VIII-E Greedy approximation bounds via coherence

This discussion in Section VIII-D pertains to exact recovery of a rank-$r$ multilinear decomposition although our main problem really takes the form of a best rank-$r$ approximation more often than not. We will describe some greedy approximation bounds for the approximation problem in this section.

By our definition of rank and border rank,

It would be wonderful if greedy algorithms along the lines of what we discussed in Section VIII-D could yield an approximant within some provable bounds that is a factor of $\sigma_{r}(\hat{f})$. However this is too much to hope for mainly because a dictionary comprising all separable functions, i.e., $\{f:\operatorname*{rank}(f)=1\}$ is far too large to be amenable to such analysis. This does not prevent us from considering somewhat more restrictive dictionaries like what we did in the previous section. So again, let $\Phi\subseteq\{f\in L^{2}(X_{1}\times\dots\times X_{d}):\operatorname*{rank}(f)=1\}$ be such that

for some given $\mu\in[0,1)$ to be chosen later. Let us instead define

since the infimum is taken over a smaller dictionary.

The special case where $t=1$ in the woga described in Section VIII-D is also called the orthogonal greedy algorithm (oga). The result we state next comes from the work of a number of people done over the last decade: (59) is due to Gilbert, Muthukrisnan, and Strauss in 2003 ; (60) is due to Tropp in 2004 ; (61) is due to Dohono, Elad, and Temlyakov in 2006 ; and (62) is due to Livshitz in 2012 . We merely apply these results to our approximation problem here.

Let $\hat{f}\in L^{2}(X_{1}\times\dots\times X_{d})$ and $f_{r}$ be the $r$th iterate as defined in woga with $t=1$ and input $\hat{f}$.

It would be marvelous if one could instead establish bounds in (59), (60), (61), and (62) with $\sigma_{r}(\hat{f})$ in place of $s_{r}(\hat{f})$ and $\{f:\operatorname*{rank}(f)=1\}$ in place of $\Phi$, dropping the coherence $\mu$ altogether. In which case one may estimate how well the $r$th oga iterates $f_{r}$ approximates the best rank-$r$ approximation. This appears to be beyond present capabilites.

We would to note that although the approximation theoretic technqiues (coherence, greedy approximation, redundant dictionaries, etc) used in this article owe their newfound popularity to compressive sensing, they owe their roots to works of the Russian school of approximation theorists (e.g., Boris Kashin, Vladimir Temlyakov, et al.) dating back to the 1980s. We refer readers to the bibliography of for more information.

VIII-F Checking coherence-based conditions

Since the conditions in Theorems 25, 27, 29, and Corollary 28 all involve coherence, we will say a brief word about its computation.

It has recently been established that computing the spark of a finite set of vectors $\Phi$, i.e., the size of the smallest linearly dependent subset of $\Phi$, is strongly NP-hard . Since $\operatorname{spark}\Phi=\operatorname{krank}\Phi+1$, it immediately follows that the same is true for Kruskal rank.

is strongly NP-hard. $\binom{\Phi}{k}=$ set of all $k$-element subsets of $\Phi$.

Given the NP-hardness of Kruskal rank, one expects that Lemma 9, as well as its finite-dimensional counterparts , would be computationally intractable to apply in reality. The reciprocal of coherence is therefore a useful surrogate for Kruskal rank by virtue of Lemma 26 and the fact that computing $\mu(\Phi)$ requires only $r(r-1)/2$ inner products $\langle\varphi_{i},\varphi_{j}\rangle$, $i\neq j$, where $r=\lvert\Phi\rvert$.

IX Applications

The goal of this section is two-fold. First we provide a selection of applications where the rank-$r$ multilinear decomposition (13) arises naturally via considerations of first principles (in electrodynamics, quantum mechanics, wave propagation, etc). Secondly, we demonstrate that the coherence conditions discussed extensively in Section VIII invariably have reasonable interpretations in terms of physical quantities.

The use of a rank-$r$ multilinear decomposition model in signal processing via higher-order statistics has a long history . Our signal processing applications here are of a different nature, they are based on geometrical properties of sensor arrays instead of considerations of higher-order statistics. This line of argument first appeared in the work of Sidiropoulos and Bro , which is innovative and well-motivated by first principles. However, like all other applications considered thus far, whether in data analysis, signal processing, psychometrics, or chemometrics, it does not address the serious nonexistence problem that we discussed at length in Section VIII-A. Without any guarantee that a solution to (28) exists, one can never be sure when the model would yield a solution. Another issue of concern is that the Kruskal uniqueness condition in Lemma 9 has often been invoked to provide evidence of a unique solution but we now know that this condition is practically impossible to check because of Corollary 31. The applications considered below would use the coherence conditions developed in Section VIII to avoid these difficulties. More precisely, Theorem 25, Theorem 27, and Corollary 28 are invoked to guarantee the existence of a solution to the approximation problem and provide readily checkable conditions for uniqueness of the solution, all via the notion of coherence. Note that unlike Kruskal’s condition, which applies only to an exact decomposition, Corollary 28 gives uniqueness of an approximation in noisy circumstances.

In this section, applications are presented in finite dimension. In order to avoid any confusion, $X^{\ast}$, $X^{\mathsf{H}}$ and $X^{\mathsf{T}}$ will denote complex conjugate, hermitian transpose, and transpose, of the matrix $X$ respectively.

Consider a narrow band transmission problem in the far field. We assume here that we are in the context of wireless telecommunications, but the same principle could also apply in other areas. Let $r$ signals impinge on an array, so that their mixture is recorded. We wish to recover the original signals and to estimate their directions of arrival and respective powers at the receiver. If the channel is specular, some of these signals can correspond to different propagation paths of the same radiating source, and are therefore correlated. In other words, $r$ does not denote the number of sources, but the total number of distinct paths viewed from the receiver.

In the present framework, we assume that channels can be time-varying, but that they can be regarded to be constant over a sufficiently short observation length. The goal is to be able to work with extremely short samples.

Under these assumptions, the signal received at discrete time $t_{k}$, $k=1,\dots,n_{3}$, on the $i$th sensor of the reference subarray can be written as

with $\psi_{i,j,p}=\jmath\frac{\omega}{C}(\mathbf{b}_{i}^{\mathsf{T}}\mathbf{d}_{p}+\boldsymbol{\Delta}_{j}^{\mathsf{T}}\mathbf{d}_{p})$. If we let $\boldsymbol{\Delta}_{1}=\mathbf{0}$, then (63) also applies to the reference subarray. The crucial feature of this structure is that variables $i$ and $j$ decouple in the function $\exp(\psi_{i,j,p})$, yielding a relation resembling the rank-retaining multilinear decomposition:

where $u_{ip}=\exp\left(\jmath\frac{\omega}{C}\mathbf{b}_{i}^{\top}\mathbf{d}_{p}\right)$, $v_{jp}=\exp\left(\jmath\frac{\omega}{C}\boldsymbol{\Delta}_{j}^{\mathsf{T}}\mathbf{d}_{p}\right)$ and $w_{kp}=\sigma_{p}(t_{k})/\|\boldsymbol{\sigma}_{p}\|$, $\lambda_{p}=\|\boldsymbol{\sigma}_{p}\|$.

However, the observation model (63) is not realistic, and an additional error term should be added in order to account for modeling inaccuracies and background noise. It is customary (and realistic thanks to the central limit theorem) to assume that this additive error has a continuous probability distribution, and that therefore the hypermatrix $S$ has the generic rank. Since the generic rank is at least as large as $\lceil n_{1}n_{2}n_{3}/(n_{1}+n_{2}+n_{3}-2)\rceil$, which is always larger than Kruskal’s bound , we are led to the problem of approximating the hypermatrix $S$ by another of rank $r$. We have seen that the angular constraint imposed in Section VIII permits us to deal with a well-posed problem. In order to see the physical meaning of this constraint, we need to first define the tensor product between sensor subarrays.

IX-B Tensor product between sensor subarrays

then we may view all measurements as the superimposition of decomposable hypermatrices $\lambda_{p}\mathbf{u}_{p}\otimes\mathbf{v}_{p}\otimes\mathbf{w}_{p}$.

Geometrical information of the sensor array is contained in $\mathbf{u}_{p}\otimes\mathbf{v}_{p}$ while energy and time information on each path $p$ is contained in $\lambda_{p}$ and $\mathbf{w}_{p}$ respectively. Note that the reference subarray and the set of translations play symmetric roles, in the sense that $\mathbf{u}_{p}$ and $\mathbf{v}_{p}$ could be interchanged without changing the whole array. This will become clear with a few examples.

When we are given a structured sensor array, there can be several ways of splitting it into a tensor product of two (or more) subarrays, as shown in the following simple examples.

This subarray is depicted in Figure 1(b). By translating it via the translation in Figure 1(c) one obtains another subarray. The union of the two subarrays yields the array of Figure 1(a). The same array is obtained by interchanging roles of the two subarrays, i.e., three subarrays of two sensors deduced from each other by two translations.

This array, depicted in Figure 2(a), can either be obtained from the union of subarray of Figure 2(b) and its translation defined by Figure 2(c), or from the array of Figure 2(c) translated three times according to Figure 2(b). We express this relationship as

In fact, \raisebox{-5.16663pt}{\includegraphics[scale={0.4}]{subarray4}}=\raisebox{-5.16663pt}{\includegraphics[scale={0.4}]{subarray2v}}\otimes\raisebox{-5.16663pt}{\includegraphics[scale={0.4}]{subarray2h}} and \raisebox{0.0pt}{\includegraphics[scale={0.4}]{subarray3}}=\raisebox{0.0pt}{\includegraphics[scale={0.4}]{subarray2h}}\otimes\raisebox{0.0pt}{\includegraphics[scale={0.4}]{subarray2h}}. However, it is important to stress that the various decompositions of the whole array into tensor products of subarrays are not equivalent from the point of view of performance. In particular, the Kruskal bound can be different, as we will see next.

Similar observations can be made for grid arrays in general.

Take an array of $9$ sensors located at $(x,y)\in\{1,2,3\}\times\{1,2,3\}$. We have the relations

Let us now have a look at the maximal number of sources $r_{\max}$ that can be extracted from a $n_{1}\times n_{2}\times n_{3}$ hypermatrix in the absence of noise. A sufficient condition is that the total number of paths, $r$, is smaller than Kruskal’s bound (25). We shall simplify the bound by making two assumptions: (a) the loading matrices are generic, i.e., they are of full rank, and (b) the number of paths is larger than the sizes $n_{1}$ and $n_{2}$ of the two subarrays entering the array tensor product, and smaller than the number of time samples, $n_{3}$. Under these simplifying assumptions, Kruskal’s bound becomes $2r_{\max}\leq n_{1}+n_{2}+r_{\max}-2$, or:

The table below illustrates the fact that the choice of subarrays has an impact on this bound.

IX-C Significance of the angular constraint

We are now in a position to interpret the meanings of the various coherences in light of this application. According to the notations given in (64), the first coherence

corresponds to the angular separation viewed from the reference subarray. The vectors $\mathbf{b}_{i}$ and $\mathbf{d}_{p}$ have unit norms, as do the vectors $\mathbf{u}_{p}$. The quantity $\lvert\mathbf{u}_{p}^{\mathsf{H}}\mathbf{u}_{q}\rvert$ may thus be viewed as a measure of angular separation between $\mathbf{d}_{p}$ and $\mathbf{d}_{q}$, as we shall demonstrate in Proposition 36.

where $\lambda=2\pi C/\omega$ denotes the wavelength.

Let $\mathbf{b}_{i}$, $\mathbf{d}_{p}$ and $\mathbf{u}_{q}$ be defined as in (64), $i=1,\dots,n_{1}$, $p,q=1,\dots,n_{2}$. Then we have the following.

If $\{\mathbf{b}_{1},\dots,\mathbf{b}_{n}\}$ is resolvent with respect to three linearly independent directions, then

Note that the condition in Definition 35 is not very restrictive, since sensor arrays usually contain sensors separated by half a wavelength or less. Thanks to Proposition 36, we now know that uniqueness of the matrix factor $U=[\mathbf{u}_{1},\dots,\mathbf{u}_{r}]$ and the identifiability of the directions of arrival $\mathbf{d}_{p}$ are equivalent. By the results of Section VIII, the uniqueness can be ensured by a constraint on coherence such as (53).

As in Section IX-B, the second coherence may be interpreted as a measure of the minimal angular separation between paths, viewed from the subarray defining translations.

The third coherence is the maximal correlation coefficient between signals received from various paths on the array

In conclusion, the best rank-$r$ approximation exists and is unique if either signals propagating through various paths are not too correlated, or if their direction of arrival are not too close, where “not too” is taken to mean that the product of coherences satisfies inequality (53) of Corollary 28. In other words, one can separate paths with arbitrarily high correlation provided they are sufficiently well separated in space.

Hence, the decomposition of a sensor array into a tensor product of two (or more) sensor subarrays depends not only on Kruskal’s bound, as elaborated in Section IX-B, but also on the ability of the latter subarrays to separate two distinct directions of arrival (cf. Proposition 36).

IX-D CDMA communications

The application to antenna array processing we described in Section IX-A also applies to all source separation problems , provided an additional diversity is available. An example is the case of Code Division Multiple Access (CDMA) communications. In fact, as pointed out in , it is possible to distinguish between symbol and chip diversities. We will elaborate on the latter example.

where $B_{kp}=\sum_{t}H_{p}(k-t)C_{p}(t)$ denotes the output of the $p$th channel excited by the $p$th coding sequence, upon removal of the guard chips (which may be affected by two different consecutive symbols) .

In order to avoid multiple access interferences, spreading sequences are usually chosen to be uncorrelated for all delays, which implies that they are orthogonal. However, the results obtained in Section VIII show that spreading sequences do not need to be orthogonal, and symbol sequences need not be uncorrelated, as long as the directions of arrival are not collinear. In particular, shorter spreading sequences may be used for the same number of users, which increases throughput. Alternatively, for a given spreading gain, one may increase the number of users. These are possible because the coherence conditions in Section VIII allow one to relax the constraint of having almost orthogonal spreading sequences. On the other hand, some directions of arrival may be collinear if the corresponding spreading sequences are sufficiently well separated angularly. These conclusions are essentially valid when users are synchronized, i.e., for downlink communications.

IX-E Polarization

The use of polarization as an additional diversity has its roots in . Several attempts to use this diversity in the framework of tensor-based source localization and estimation can be found in the literature .

Coherences $\mu_{1}$ and $\mu_{3}$ are the same as in Section IX-C, and represent respectively the angular separation between directions of arrival, and correlation between arriving sources. It is slightly more difficult to see the significance of $\mu_{2}$, the coherence associated with polarization.

For this, we need to go into more details . Let $\alpha_{p}\in(-\pi/2,\pi/2]$ and $\beta_{p}\in(-\pi/4,0)\cup(0,\pi/4)$ denote respectively the orientation and ellipticity angles of the polarization of the $p$th wave. Let $\theta_{p}\in[0,2\pi)$ and $\phi_{p}\in(-\pi/2,\pi/2]$ denote respectively the azimuth and elevation of the direction of arrival of the $p$th path. We have

The unit vector defining the $p$th direction of arrival is

So the triplet $(\mathbf{d}_{p},\mathbf{e}_{p},\mathbf{f}_{p})$ forms a right orthonormal triad.

First note that $Q(\alpha_{p})^{\mathsf{H}}Q(\alpha_{q})=Q(\alpha_{q}-\alpha_{p})$. Hence $\mathbf{g}_{p}^{\mathsf{H}}\mathbf{g}_{q}$ can be of unit modulus only if $\mathbf{h}_{p}$ and $Q(\alpha_{q}-\alpha_{p})\mathbf{h}_{q}$ are collinear. But the first entry of $\mathbf{h}_{p}$ is real and the second is purely imaginary. So the corresponding imaginary and real parts of $Q(\alpha_{q}-\alpha_{p})\mathbf{h}_{q}$ must be zero, which implies that $\sin(\alpha_{q}-\alpha_{p})=0$. Consequently $Q(\alpha_{q}-\alpha_{p})=\pm I$, which yields $\mathbf{h}_{p}=\pm\mathbf{h}_{q}$. But because the angle $\beta$ lies in the interval $(-\pi/4,\pi/4)$, only the positive sign is acceptable. ∎

We have $\lvert\mathbf{v}_{p}^{\mathsf{H}}\mathbf{v}_{q}\rvert=\lvert\mathbf{g}_{p}^{\mathsf{H}}B_{p}^{\mathsf{T}}B_{q}\mathbf{g}_{q}\rvert$. Notice that the matrix $B_{p}^{\mathsf{T}}B_{q}$ is of the form

where $\gamma$ and $\eta$ are real, $\gamma=\frac{1}{2}(\mathbf{e}_{p}^{\mathsf{T}}\mathbf{e}_{q}+\mathbf{f}_{p}^{\mathsf{T}}\mathbf{f}_{q})$ and $\eta=\frac{1}{2}(\mathbf{e}_{p}^{\mathsf{T}}\mathbf{f}_{q}-\mathbf{f}_{p}^{\mathsf{T}}\mathbf{e}_{q})$. Since $\mathbf{g}_{p}$ and $\mathbf{g}_{q}$ are of unit norms, $\lvert\mathbf{v}_{p}^{\mathsf{H}}\mathbf{v}_{q}\rvert$ can be of unit modulus only if $B_{p}^{\mathsf{T}}B_{q}$ has an eigenvalue of unit modulus, which requires that $\gamma^{2}+\eta^{2}=1$. We now prove that $\gamma^{2}+\eta^{2}\leq 1$ with equality if and only if the four sets of equalities hold.

With this goal in mind, define the $6$-dimensional vectors

Then $\gamma=\mathbf{z}^{\mathsf{T}}\mathbf{w}$ and $\gamma=\mathbf{z}^{\mathsf{T}}\mathbf{w}^{\prime}$. Decompose $\mathbf{z}$ into two orthogonal parts: $\mathbf{z}=\mathbf{z}_{0}+\mathbf{z}_{1}$, with $\mathbf{z}_{0}\in\operatorname{span}\{\mathbf{w},\mathbf{w}^{\prime}\}$ and $\mathbf{z}_{0}\bot\mathbf{z}_{1}$. Clearly, $\gamma^{2}+\eta^{2}=\|\mathbf{z}_{0}\|^{2}$. Moreover, $\|\mathbf{z}_{0}\|^{2}\leq\|\mathbf{z}\|^{2}=1$, with equality if and only if $\mathbf{z}\in\operatorname{span}\{\mathbf{w},\mathbf{w}^{\prime}\}$. By inspection of the definitions of $\mathbf{e}_{p}$ and $\mathbf{e}_{q}$, we see that the third entry of $\mathbf{z}$ and $\mathbf{w}$ is $\mathbf{0}$. Hence $\mathbf{z}\in\operatorname{span}\{\mathbf{w},\mathbf{w}^{\prime}\}$ is possible only if either $\mathbf{z}$ is collinear to $\mathbf{w}$ or if the third entry of $\mathbf{w}^{\prime}$ is $\mathbf{0}$. In the latter case, it means that $\phi_{q}=\pi/2$, and so $\phi_{p}=\pi/2$ and $\theta_{p}=\theta_{q}$. In the former case, it can be seen that $\sin\theta_{p}=\sin\theta_{q}$, and finally that $\phi_{p}=\phi_{q}$.

The last step is to rewrite $\gamma$ and $\eta$ as a function of angle $\theta_{p}-\theta_{q}$, using trigonometric relations: $\gamma=\cos(\theta_{p}-\theta_{q})(1+\sin\phi_{p}\sin\phi_{q})+\cos\phi_{p}\cos\phi_{q}$ and $\eta=\sin(\theta_{p}-\theta_{q})(\sin\phi_{p}+\sin\phi_{q})$. This eventually shows that $\gamma=1$ and $\eta=0$. As a consequence, $\lvert\mathbf{v}_{p}^{\mathsf{H}}\mathbf{v}_{q}\rvert=1$ only if $B_{p}^{\mathsf{T}}B_{q}=I$, and the result follows from Lemma 37. ∎

Proposition 38 shows that a constraint on the coherence $\mu_{2}$ compels source paths to have either different directions of arrival or different polarizations, giving $\mu_{2}$ physical meaning.

IX-F Fluorescence spectral analysis

This reveals the true chemical factors responsible for the data: $r=\operatorname{rank}(A)$ gives the number of pure substances in the mixtures, $\mathbf{x}_{p}=(x_{1p},\dots,x_{lp})$ gives the relative concentrations of $p$th substance in specimens $1,\dots,l$; $\mathbf{y}_{p}=(y_{1p},\dots,y_{mp})$ gives the excitation spectrum of $p$th substance; and $\mathbf{z}_{p}=(z_{1p},\dots,z_{np})$ gives the emission spectrum of $p$th substance. The emission and excitation spectra would then allow one to identify the pure substances.

Of course, this is only valid in an idealized situation when measurements are performed perfectly without error and noise. Under realistic noisy circumstances, one would then need to a find best rank-$r$ approximation, which is where the coherence results of Section VIII play a role. In this case, $\mu(\mathbf{x}_{1},\dots,\mathbf{x}_{r})$ measures the relative abundance of the pure substances in the samples while $\mu(\mathbf{y}_{1},\dots,\mathbf{y}_{r})$ and $\mu(\mathbf{z}_{1},\dots,\mathbf{z}_{r})$ measure the spectroscopic likeness of these pure substances in the sense of absorbance and fluorescence respectively.

IX-G Statistical independence induces diversity

We will now discuss a somewhat different way to achieve diversity. Assume the linear model below

where only the signal $\mathbf{x}(t)$ is observed, $U=[\mathbf{u}_{1},\dots,\mathbf{u}_{r}]$ is an unknown $n\times r$ mixing matrix, and $\mathbf{s}(t)=(s_{1}(t),\dots,s_{r}(t))$ has mutually statistically independent components. One may construct $K_{d}(\mathbf{x})$, the $d$th order cumulant hypermatrix of $\mathbf{x}(t)$, and it will satisfy the multilinear model

where $\lambda_{p}(\mathbf{s})$ denotes the $p$th diagonal entry of the $d$th cumulant hypermatrix of $\mathbf{s}$. Because of the statistical independence of $\mathbf{s}(t)$, the off-diagonal entries of the $d$th cumulant hypermatrix of $\mathbf{s}$ are zero . If $d\geq 3$, then the matrix $U$ and the entries $\lambda_{p}(\mathbf{s})$ can be identified . One may apply the results of Section VIII to deduce uniqueness of the solution.

Such problems generalize to convolutive mixtures and have applications in telecommunications, radar, sonar, speech processing, and biomedical engineering .

IX-H Nonstationarity induces diversity

If a signal $x(t)$ is nonstationary, its time-frequency transform, defined by

for some given kernel $\kappa$, bears information. If variables $t$ and $f$ are discretized, then the values of $X(t,f)$ can be stored in a matrix $X$; and the more nonstationary the signal $x(t)$, the larger the rank of $X$. A similar statement can be made on a signal $y(\mathbf{z})$ depending on a spatial variable $\mathbf{z}$. The discrete values of the space-wavevector transform $Y(\mathbf{z},\mathbf{w})$ of a field $y(\mathbf{z})$ can be stored in a matrix $Y$; and the less homogeneous the field $y(\mathbf{z})$, the larger the rank of $Y$. This is probably the reason why algorithms proposed in permit one to localize and extract dipole contributions in the brains using a multilinear model, provided that one has distinct time-frequency or space-wavevector patterns. Nevertheless, such localization is guaranteed to be successful only under restrictive assumptions.

X Further work

A separate article discussing practical algorithms for the bounded coherence best rank-$r$ multilinear approximation is under preparation with additional coauthors. These algorithms follow the general strategy of the greedy approximations woga and oga discussed in Sections VIII-D and VIII-E but contain other elements exploiting the special separable structure of our problem. Extensive numerical experiments will be provided in the forthcoming article.

Acknowledgement

We thank Ignat Domanov for pointing out that in Theorem 25, the factor on the right hand side may be improved from $r$ to $r-1$, and that the improvement is sharp. We owe special thanks to Harm Derksen for pointing out an error in an earlier version of Section VII and for very helpful discussions regarding nuclear norm of tensors. Sections VIII-D and VIII-E came from an enlightening series of lectures Vladimir Temlyakov gave at the IMA in Minneapolis and the useful pointers he graciously provided afterwards. We gratefully acknowledge Tom Luo, Nikos Sidiropoulos, Yuan Yao, and two anonymous reviewers for their helpful comments.

The work of LHL is partially supported by AFOSR Young Investigator Award FA9550-13-1-0133, NSF Collaborative Research Grant DMS 1209136, and NSF CAREER Award DMS 1057064. The work of PC is funded by the European Research Council under the European Community’s Seventh Framework Programme FP7/2007–2013 Grant Agreement no. 320594.

References