Dictionary-based Tensor Canonical Polyadic Decomposition

TL;DR

This paper introduces a dictionary-based tensor canonical polyadic decomposition method that enhances interpretability and accuracy in tensor source extraction by enforcing known dictionary constraints, with demonstrated benefits in hyperspectral image unmixing.

Contribution

It proposes a novel formulation of sparse coding for high-dimensional tensors and explores the benefits of dictionary constraints in tensor decomposition.

Findings

Improved parameter identifiability and estimation accuracy.

Effective in hyperspectral image unmixing.

Validated on simulated data and real hyperspectral images.

Abstract

To ensure interpretability of extracted sources in tensor decomposition, we introduce in this paper a dictionary-based tensor canonical polyadic decomposition which enforces one factor to belong exactly to a known dictionary. A new formulation of sparse coding is proposed which enables high dimensional tensors dictionary-based canonical polyadic decomposition. The benefits of using a dictionary in tensor decomposition models are explored both in terms of parameter identifiability and estimation accuracy. Performances of the proposed algorithms are evaluated on the decomposition of simulated data and the unmixing of hyperspectral images.

Tables

Table 1. TABLE I: Basic definitions from linear algebra

	$ℰ\otimes ℱ$	$\text{: tensor product space, linear space mapped by }\otimes \text{ from }ℰ\times ℱ.$
	$𝒂\otimes 𝒃$	: tensor product of two vectors, i.e. an element of $ℰ\otimes ℱ$,
		can be understood as an outer product of vectors
		if the tensor space is an array space. [14]
	$𝑨\odot 𝑩$	: Khatri-Rao (columnwise Kronecker) product of matrices [19].
	$𝑨⊡𝑩$	: Hadamard product of matrices, i.e., element-wise product [19].

Table 2. TABLE II: Some definitions and properties of multilinear operators

	$𝑼\otimes 𝑽\otimes 𝑾$	: an operator acting on a third order tensor.
	$\left(𝑼\otimes 𝑽\otimes 𝑾\right)𝓣$	$\text{: application of }𝑼\text{ on the first mode, }𝑽\text{ on the}$
	$\text{second mode, }𝑾$	$\text{on the third mode, also noted }𝓣{\bullet }_1𝑼{\bullet }_2𝑽{\bullet }_3𝑾.$

Table 3. TABLE III: Average runtime of various algorithms over N = 100 𝑁 100 N=100 realization when R = R ​ e 𝑅 𝑅 𝑒 R=Re .

Alg.	ALS	MPALS	Flex-MPALS
Time (s)	0.21	0.56	0.72
Alg.	r-MPALS	SMPALS	FG
Time (s)	2.75	0.73	7.62

Table 4. TABLE IV: Reconstruction error ( % ) (\%) for the Urban and Terrain airport HSI. Best results are highlighted in bold. Computation time does not include 500 iterations of nonnegative least squares update for the abundances after each method (which takes about 10 seconds on average).

	Urban HSI		Terrain HSI
	time(s)	Rel.err.	time(s)	Rel.err.
SPA	0.2	9.58	0.3	5.89
d-SPA	20	4.57	48	3.50
ds-SPA	21	4.67	58	3.37
nmf-SPA	13	5.37	18	3.75
VCA	2.0	13.07	1.6	18.61
d-VCA	16	4.73	46	3.29
ds-VCA	20	4.66	35	3.37
nmf-VCA	12	6.63	17	4.27
SNPA	9.5	9.63	13	5.76
d-SNPA	17	5.02	47	3.78
ds-SNPA	21	4.91	59	3.88
nmf-SNPA	13	5.14	19	4.47
H2NMF	8	5.81	11.7	5.09
d-H2NMF	17	4.05	47	3.49
ds-H2NMF	21	4.05	64	3.49
nmf-H2NMF	11	5.03	17	4.51
FGNRS-100	2.2	5.58	1.6	3.34
d-FGNSR-100	17	4.47	45	3.01
ds-FGNSR-100	20	4.47	58	3.01
nmf-FGNSR-100	12	6.37	19	3.31

Equations

$$\mathcal{T}=\sum _{r=1}^R{\mathcal{D}}_r,$$

$$\mathcal{T}=\left(\mathbf{A}\otimes \mathbf{B}\otimes \mathbf{C}\right){\mathcal{I}}_R,$$

$$\begin{array}[]{l}\boldsymbol{T}_{1}=\boldsymbol{A}\left(\boldsymbol{B}\odot\boldsymbol{C}\right)^{T},\\ \boldsymbol{T}_{2}=\boldsymbol{B}\left(\boldsymbol{A}\odot\boldsymbol{C}\right)^{T},\\ \boldsymbol{T}_{3}=\boldsymbol{C}\left(\boldsymbol{A}\odot\boldsymbol{B}\right)^{T}.\end{array}$$

$$\mathbf{M}=\mathbf{A}{\mathbf{B}}^T=\mathbf{A}\mathbf{P}{\mathbf{P}}^{-1}{\mathbf{B}}^T$$

$$\mathbf{B}=\mathbf{D}\mathbf{S},\Vert {\mathbf{s}}_i{\Vert }_0=1\text{ for }i\in \{1,\dots ,R\},$$

$$\left\{\begin{array}[]{l}\boldsymbol{\mathcal{T}}=\left(\boldsymbol{A}\mathop{\otimes}\boldsymbol{D}\boldsymbol{S}\mathop{\otimes}\boldsymbol{C}\right)\boldsymbol{\mathcal{I}}_{R}+\boldsymbol{\mathcal{E}},\\ \|\boldsymbol{s}_{i}\|_{0}=1\text{ for }i\in\{1,\dots,R\},~{}\boldsymbol{S}\in\{0,1\}^{d\times R}.\end{array}\right.$$

$$\begin{array}[]{l}\underset{\boldsymbol{A},\boldsymbol{S},\boldsymbol{C}}{\text{argmin}}\,\|\boldsymbol{\mathcal{T}}-\left(\boldsymbol{A}\mathop{\otimes}\boldsymbol{D}\boldsymbol{S}\mathop{\otimes}\boldsymbol{C}\right)\boldsymbol{\mathcal{I}}_{R}\|_{F}^{2}\\ \text{such that }\|\boldsymbol{s}_{i}\|_{0}=1\text{ for }i\in\{1,\dots,R\}\text{ and }\boldsymbol{S}\in\{0,1\}^{d\times R}.\end{array}$$

$$\mathbf{B}=f(\mathbf{D}\mathbf{S},\theta ),$$

$$\begin{array}[]{l}\underset{\boldsymbol{A},\boldsymbol{B},\boldsymbol{S},\boldsymbol{C}}{\text{argmin}}\,\|\boldsymbol{\mathcal{T}}-\left(\boldsymbol{A}\mathop{\otimes}\boldsymbol{B}\mathop{\otimes}\boldsymbol{C}\right)\boldsymbol{\mathcal{I}}_{R}\|_{F}^{2}+\frac{1}{\sigma_{c}^{2}}\|\boldsymbol{B}-\boldsymbol{D}\boldsymbol{S}\|_{F}^{2}\\ \text{such that }\|\boldsymbol{s}_{i}\|_{0}=1\text{ for }i\in\{1,\dots,R\}\text{ and }\boldsymbol{S}\in\{0,1\}^{d\times R}.\end{array}$$

$$\left\{\begin{array}[]{l}\boldsymbol{M}_{m}=\boldsymbol{A}\text{Diag}(\boldsymbol{c}_{m})\boldsymbol{B}_{m}^{T},\\ \boldsymbol{B}_{m}=\boldsymbol{P}_{m}\boldsymbol{Q},\\ \boldsymbol{P}_{m}^{T}\boldsymbol{P}_{m}=\boldsymbol{I}_{R},\\ \boldsymbol{Q}=\boldsymbol{D}\boldsymbol{S},\\ \|\boldsymbol{s}_{i}\|_{0}=1\text{ for }i\in\{1,\dots,R\},~{}\boldsymbol{S}\in\{0,1\}^{d\times R},\end{array}\right.$$

$$\mathbf{M}=\mathbf{A}{\left({\mathbf{M}}^T\mathbf{S}\right)}^T,$$

$$\mathbf{M}=\mathbf{D}\mathbf{X},$$

$$||X|{|}_{{\ell }_p/{\ell }_q}=||v|{|}_p,\text{ where }{v}_i=||{\mathbf{x}}_i|{|}_q\forall i.$$

$$\begin{array}[]{l}\min_{\boldsymbol{X}\in\mathbb{R}_{+}^{d\times n}}||\boldsymbol{M}-\boldsymbol{X}\boldsymbol{M}||_{F}^{2}\\ \text{such that }\;\text{ $\boldsymbol{X}$ has $r$ non-zero columns, }\end{array}$$

$$\underset{\mathbf{A},\mathbf{C}}{\inf }\Vert \mathcal{T}-\left(\mathbf{A}\otimes \mathbf{B}\otimes \mathbf{C}\right){\mathcal{I}}_R{\Vert }_F^2$$

$$\underset{\mathbf{A},||{\mathbf{a}}_r|{|}_2=1\forall r,\mathbf{C}}{\inf }\Vert {\mathbf{T}}_3-\mathbf{C}{\left(\mathbf{A}\odot \mathbf{B}\right)}^T{\Vert }_F^2$$

$$\Vert {\mathbf{T}}_3-\mathbf{C}{\left(\mathbf{A}\odot \mathbf{B}\right)}^T{\Vert }_F\le \Vert {\mathbf{T}}_3{\Vert }_F$$

$${\sigma }_R(\mathbf{A}\odot \mathbf{B})\ge {\sigma }_R(\mathbf{B})>0,$$

$${\sigma }_R^2(\mathbf{A}\odot \mathbf{B})$$

$$=\underset{||x|{|}_2=1}{\min }\sum _j||\mathbf{B}(x⊡{\mathbf{a}}^j)|{|}_2^2$$

$$\ge {\sigma }_R^2(\mathbf{B})\underset{||x|{|}_2=1}{\min }\sum _j||x⊡{\mathbf{a}}^j|{|}_2^2$$

$$={\sigma }_R^2(\mathbf{B})\underset{||x|{|}_2=1}{\min }\sum _{r=1}^R|x_r|||{\mathbf{a}}_r|{|}_2^2={\sigma }_R^2(\mathbf{B}).$$

$$\hat{\mathbf{B}}={\mathbf{T}}_2{\left(\mathbf{A}\odot \mathbf{C}\right)}^{†}={\mathbf{T}}_2\left(\mathbf{A}\odot \mathbf{C}\right){\left({\mathbf{A}}^T\mathbf{A}⊡{\mathbf{C}}^T\mathbf{C}\right)}^{-1},$$

$$\hat{\mathbf{B}}=\left({\mathbf{T}}_2\left(\mathbf{A}\odot \mathbf{C}\right)+\lambda \mathbf{D}\mathbf{S}\right){\left({\mathbf{A}}^T\mathbf{A}⊡{\mathbf{C}}^T\mathbf{C}+\lambda {\mathbf{I}}_{R\times R}\right)}^{-1}$$

$$\mathbf{B}=\left({\mathcal{T}}_{(2)}\left(\mathbf{A}\odot \mathbf{C}\right)+\lambda \mathbf{D}\mathbf{S}\right){\left({\mathbf{A}}^{\mathsf{T}}\mathbf{A}⊡{\mathbf{C}}^{\mathsf{T}}\mathbf{C}+\lambda {\mathbf{I}}_{R\times R}\right)}^{-1}$$

$$\mathbf{B}=\mathbf{D}\mathbf{S}+O\left(\frac{1}{\lambda }\right).$$

$$\underset{x\in \mathcal{B}(1,{\ell }_2)}{\mathrm{argmin}}\Vert x{\Vert }_1$$

$$\left\{\begin{array}[]{l}\min_{\boldsymbol{S}}\gamma(\boldsymbol{S})=\frac{1}{2}\|\boldsymbol{\mathcal{T}}-\left(\boldsymbol{A}\mathop{\otimes}\boldsymbol{D}\boldsymbol{S}\mathop{\otimes}\boldsymbol{C}\right)\boldsymbol{\mathcal{I}}_{R}\|_{F}^{2}+\delta\|S\|_{1}\\ \text{ such that }\|\boldsymbol{s}_{i}\|_{2}^{2}=1~{}~{}\forall i\leq R,~{}~{}\boldsymbol{S}\geq 0,\end{array}\right.$$

$$\frac{\partial \gamma }{\partial \mathbf{S}}={\mathbf{D}}^T\mathbf{D}\mathbf{S}\left({\mathbf{A}}^T\mathbf{A}⊡{\mathbf{C}}^T\mathbf{C}\right)-{\mathbf{D}}^T{\mathbf{T}}_2\left(\mathbf{A}\odot \mathbf{C}\right)+\delta 1_{d\times R}$$

$$\text{step}=\min \left(\frac{1}{{ϵ}_S},\underset{j}{\min }\underset{i|g_S(i,j)>0}{\max }\frac{S(i,j)}{g_S(i,j)}-10^{-12}\right)$$

$$\begin{array}[]{ll}\boldsymbol{d}^{i}_{k}=&|a_{k}\boldsymbol{u}+b_{k}+\nu_{k}^{i}\text{sinc}(\frac{\pi}{6}\boldsymbol{u}-e_{k}^{i})\\ &+\mu_{k}^{i}\left(\text{tri}(\boldsymbol{u}-f_{k}^{i}+2)-\text{tri}(\boldsymbol{u}-f_{k}^{i}-2)\right)|\end{array}$$

$$\mathbf{C}\leftarrow {\mathbf{C}}^{(0)}\left(\rho {\mathbf{I}}_R+\frac{(1-\rho )}{R}{1}_{R\times R}\right),$$

Full text

I Introduction

Given a mixture of several components, a classic problem in signal processing is to separate the contribution of each of these components using solely the information contained in the data. This problem is known as blind source separation and has been a particularly widely studied topic for the last two decades [1]. When the available data are contained in a multiway array, that is, a table of three entries or more, blind source separation has been successfully achieved by means of tensor decomposition techniques in a large variety of applications, ranging from telecommunications [2] to chemometrics [3], spectral unmixing [4], neuroimaging [5], social sciences [6] and machine learning [7].

The key concept behind tensor decomposition methods is the linearity of blocks of parameters of interest such as spectra, concentrations or time signatures with respect to experimental parameters such as wavelength, pixel index, subject index or time. For the canonical polyadic decomposition (CPD) model studied in this paper, each block of parameters of interest depends respectively on only one experimental parameter.

As an illustration, tensor decomposition techniques can be used to perform source separation with hyperspectral images, a task often referred to as spectral unmixing. Hyperspectral images are 2D images collected for a large number of wavelengths, and possibly along time. In this scenario, the blocks of parameters of interest are spectral signatures (characteristic responses of materials to light stimulation), abundances (relative concentrations) and time evolution.

For most applications, tensor decomposition models are however not exactly following a physical models because of their strong multilinearity assumptions. As a consequence, modeling error has to be accounted for when processing the results of the source separation using tensors. This matters when the goal of the source separation is to identify the components in the mixture. The estimated parameters of interest have to be compared with benchmarks, and modeling error induces error on the estimated parameters which may in turn deteriorate identification performances.

In this paper, we want to avoid splitting the identification procedure and the source separation procedure. It is shown that using a formalism inspired from sparse coding [8], merging source separation and identification is not only possible but also offers advantages in term of uniqueness properties of the tensor decomposition, and may reduce estimation error on identified factors provided the a priori information is accurate.

Outline and contributions

In this work, after discussing notations and vocabulary, we make the following contributions:

• •

Formalize high order tensor sparse coding by modifying the usual sparse coding formulation for matrices [8, 9]. We call the obtained tensor model dictionary canonical polyadic decomposition (DCPD).

• •

Provide tools for introducing flexibility in the DCPD formulation. This makes the DCPD model more suitable for practical problems.

• •

Study the identifiability of parameters of the DCPD for matrices and higher order tensors, and study the existence of a best low rank DCPD approximation.

• •

Develop greedy and continuous algorithms to compute the DCPD and its flexible variants.

• •

Check the identification performances of the DCPD model on synthetic data with respect to the CPD model. We also use DCPD in the matrix case to perform spectral unmixing under pure-pixel assumption of the Urban and Terrain data sets111available at http://www.agc.army.mil/ and compare our results with state-of-the-art methods [10, 11, 12, 13].

Notation and Vocabulary

Among various notation habits in the multiway array processing community, we choose to follow notation from [14, 15], as presented in Tables I and II. Although our results will be applicable to tensors of any order, we focus in this paper on third order tensor in order to simplify the presentation. We call a third-order real $K\times L\times M$ tensor $\boldsymbol{\mathcal{T}}$ a vector from a tensor space $(\mathds{R}^{K}\mathop{\otimes}\mathds{R}^{L}\mathop{\otimes}\mathds{R}^{M},\mathop{\otimes})$ with $\mathop{\otimes}$ being a tensor product. A three way array is an element from $(\mathds{R}^{K\times L\times M},\mathop{\otimes})$ where $\mathop{\otimes}$ is the outer product. In other words, the set of arrays is a tensor space when the outer product is used as the tensor product. A higher order tensor is a vector from a tensor space featuring at least three linear subspaces.

In Table II, some useful properties of multilinear operators acting on tensors are specified. Multilinear operators generalize linear operators acting on vectors by defining linear operations on each vector space $\mathds{R}^{N}$ composing the tensor space. These multilinear operators also form a tensor space $(\mathds{R}^{R_{1}\times K}\mathop{\otimes}_{op}\mathds{R}^{R_{2}\times K}\mathop{\otimes}_{op}\mathds{R}^{R_{3}\times K},\mathop{\otimes}_{op})$ where $\mathop{\otimes}_{op}$ will be abusively denoted $\mathop{\otimes}$ although it is not the outer product. Indeed, $\mathop{\otimes}_{op}$ is defined by $(\boldsymbol{U}\mathop{\otimes}_{op}\boldsymbol{V})(\boldsymbol{a}\mathop{\otimes}\boldsymbol{b})=\boldsymbol{U}\boldsymbol{a}\mathop{\otimes}\boldsymbol{V}\boldsymbol{b}$.

Given a $K\times L\times M$ tensor $\boldsymbol{\mathcal{T}}$, its canonical polyadic decomposition (CPD) of rank $R$ can be written as follows:

[TABLE]

where $\boldsymbol{\mathcal{D}}_{r}$ are decomposable tensors of the form $\boldsymbol{\mathcal{D}}_{r}=\boldsymbol{a}_{r}\mathop{\otimes}\boldsymbol{b}_{r}\mathop{\otimes}\boldsymbol{c}_{r}$. The rank of $\boldsymbol{\mathcal{T}}$ is the minimal value of $R$ such that (1) holds exactly, while we define the rank of a CPD model as the number of component $R$ in that model. A CPD model or a tensor is said to be a low-rank model or a low-rank tensor if $R$ is small with respect to all dimensions of the data or tensor.

Finding the CPD of a third-order tensor means finding rank-one tensors $\boldsymbol{\mathcal{D}}_{r}$. Yet, each tensor $\boldsymbol{\mathcal{D}}_{r}$ may be defined by three vectors $\boldsymbol{a}_{r}$, $\boldsymbol{b}_{r}$ and $\boldsymbol{c}_{r}$, only up to two scaling ambiguities; in fact, $\boldsymbol{a}_{r}\mathop{\otimes}\boldsymbol{b}_{r}\mathop{\otimes}\boldsymbol{c}_{r}=\alpha\boldsymbol{a}_{r}\mathop{\otimes}\beta\boldsymbol{b}_{r}\mathop{\otimes}\boldsymbol{c}_{r}/\alpha\beta$, $\forall\alpha,\beta\neq 0$.

Next, it is often convenient to store these vectors in matrices called factors as $\boldsymbol{A}=[\boldsymbol{a}_{1},\dots,\boldsymbol{a}_{R}]$, $\boldsymbol{B}=[\boldsymbol{b}_{1},\dots,\boldsymbol{b}_{R}]$ and $\boldsymbol{C}=[\boldsymbol{c}_{1},\dots,\boldsymbol{c}_{R}]$. This leads to a convenient writing:

[TABLE]

where $\boldsymbol{A}\in\mathds{R}^{K\times R}$, $\boldsymbol{B}\in\mathds{R}^{L\times R}$ and $\boldsymbol{C}\in\mathds{R}^{M\times R}$ are called factor matrices, and $\boldsymbol{\mathcal{I}}_{R}=\sum\limits_{r=1}^{R}{\boldsymbol{e}_{i}\mathop{\otimes}\boldsymbol{e}_{i}\mathop{\otimes}\boldsymbol{e}_{i}}$ belongs to $\mathds{R}^{R\times R\times R}$ with $\boldsymbol{e}_{i}$ a canonical basis vector of $\mathds{R}^{R}$, that is, $\boldsymbol{\mathcal{I}}_{R}$ is a diagonal core tensor with only ones on the diagonal. A convenient interpretation of (2) is to see it as a change of basis. Indeed, (2) means that the vector $\boldsymbol{\mathcal{T}}$ is expressed by coefficients $\boldsymbol{\mathcal{I}}_{R}$ in the image of multilinear operator $\boldsymbol{A}\mathop{\otimes}\boldsymbol{B}\mathop{\otimes}\boldsymbol{C}$, with linear operators $\boldsymbol{A}$, $\boldsymbol{B}$ and $\boldsymbol{C}$ spanning respectively the first, second and third mode of $\boldsymbol{\mathcal{T}}$. Note that model (2) now contains $2R$ scaling indeterminacies whereas definition (1) did not contain any.

Conditions on the dimensions of the tensor and rank of the decomposition are given in the literature [16, 17, 18] to ensure uniqueness of the factors in an unconstrained CP model, but only when noise is absent. When these conditions are satisfied, the factors in the CPD model are unique up to permutation ambiguity and the scaling described above. The model is then said to be identifiable. For higher order tensors, that is, arrays of three ways and more, these conditions are mild.

Matrices are a particular case of tensors with only two modes. One of the main differences of the CPD for matrices compared with the CPD of higher order tensors is that it is never identifiable without additional constraints as soon as $R>1$. Note that, since the two modes of a matrix are the column space and the row space, the CPD of a matrix can be written as $\boldsymbol{M}=\boldsymbol{A}\boldsymbol{B}^{T}$.

A difficult problem not adressed in this manuscript is finding the rank of a tensor, in both an exact decomposition and an approximate decomposition scenario. A naive approach is to look at the singular value profiles of the unfoldings, but an interested reader may refer to [20, 21] and references therein.

A tool that needs to be introduced here is the unfolding of three-way arrays. By choosing three particular ways to cut the cube of data into slices and stacking the obtained matrices, it is possible to rewrite the CPD in a matrix format. We chose the unfolding defined in [15], yielding the following matrix format CPD :

[TABLE]

Unfolding $\boldsymbol{T}_{i}$ contains all the vectors along the mode $i$. For instance, $\boldsymbol{T}_{1}$ contains all the $L\times M$ columns of $\boldsymbol{\mathcal{T}}$, thus its column space is the subspace spanned on the first mode by $\boldsymbol{\mathcal{T}}$, that is, the span of columns of $\boldsymbol{A}$.

II Dictionnary CPD : Model Definition

An important property of low-rank models is that the computed factors bear physical meaning. Factors are interpretable when the multilinear relationship between the experimental parameters and the block of parameters of interest stems from a meaningful modeling, and when the model is identifiable.

The first condition is met in many applications of low-rank models, for instance in fluorescence spectroscopy where spectra and concentrations depend linearly, as a first-order approximation, on the emission wavelength, the excitation wavelength and the mixture index [22]. When mining a collection of text documents, low-rank factorization techniques can identify characteristic words and documents [23].

The identifiability condition may however not be verified in practice. For matrices, it is indeed well known that identifiability is never achieved without imposing additional constraints for $R$ greater than one, since

[TABLE]

for any invertible matrix $\boldsymbol{P}$. In the higher-order case, conditions on the rank of the model and the dimensions of the data have been reported in the literature with some variations, but are mild enough to ensure identifiability of the CPD in most applications [16, 17, 18].

If the factors are meaningful, it may be of crucial importance to identify them. In the example of fluorescence spectroscopy, once emission and excitation spectra have been extracted from the data using a CPD model, the final step is to recognize which chemical compound is present in the mixtures by matching the extracted spectra with known spectra of known chemicals. The goal of the models presented below is to merge this identification step with the source separation procedure using the a priori information available on the factors.

II-A Dictionary-based CPD

Formally, assume the following relationship between a matrix $\boldsymbol{D}\in\mathds{R}^{L\times d}$ called the (over-complete) dictionary and factor $\boldsymbol{B}$:

[TABLE]

where $\boldsymbol{S}\in\{0,1\}^{d\times R}$ is a binary matrix, $\|\boldsymbol{s}_{i}\|_{0}$ counts the number of non-zero values in the $i^{th}$ column of $\boldsymbol{S}$, and $d$ is much larger than $R$. Here $\boldsymbol{S}$ has exactly one $1$ in each column and is a selection matrix that identifies, among the $d$ atoms of the dictionary, $R$ atoms present in factor $\boldsymbol{B}$. The chosen columns of $\boldsymbol{D}$ that constitute $\boldsymbol{B}$ are given by the row index of all the ones in $\boldsymbol{S}$. A priori, no restrictions on the atoms are assumed although some technical conditions to ensure identifiability in particular settings are discussed in Section IV.

Viewing (5) as a parameterization of $\boldsymbol{B}$, a dictionary CPD model (DCPD) can be written as follows:

[TABLE]

In this model, the parameters are $\boldsymbol{A}$, $\boldsymbol{S}$ and $\boldsymbol{C}$. Tensor $\boldsymbol{\mathcal{E}}$ represent the noise, and typically its entries follow i.i.d. known distributions such as a zero-mean Gaussian distribution. In that case, it is possible to derive the maximum likelihood estimator of parameters $\boldsymbol{A}$, $\boldsymbol{S}$ and $\boldsymbol{C}$:

[TABLE]

Of course, the DCPD does not apply in any situation where the CPD works. It is far from obvious that for a given application, a dictionary containing exactly all the factors on one mode is available. Even if such a dictionary is available, it is likely that some variability has to be accounted for between the atoms and the factors. Below, we try to tackle respectively cases where the relation (5) should not be exact, and cases where no dictionary is available.

II-B Accounting for Variability

In spectral unmixing, two different approaches have been studied in the literature: the fully blind unmixing and the semi-blind approach, using already known spectra to compute regression [24]. However in this particular application, it is well known that the same material emits a slightly different spectrum depending on additional parameters not accounted for in a low-rank model. These variations are called spectral variability and are a major issue to using both blind unmixing and semi-blind unmixing [25]. Learning from this example, we emphasize that a naive DCPD as introduced above in (6) may be viewed as unrealistic.

Therefore, it may be necessary to introduce flexibility in the DCPD model. A first way is to generalize the relationship between factor $\boldsymbol{B}$ and the dictionary $\boldsymbol{D}$:

[TABLE]

for some function $f(x,\theta)$ mapping to $\mathds{R}^{L\times R}$ and where $\theta$ is a random variable following some known probability $p(\theta)$. A simple instance of (8) is obtained by setting $f(x,\theta)=x+\theta$ and $p(\theta)$ is a Gaussian distribution of zero mean and known white covariance $\sigma_{c}\boldsymbol{I}_{K}\mathop{\otimes}\boldsymbol{I}_{R}$. Then (8) is a noisy version of (5) and interpreting (8) as a priori information on $\boldsymbol{B}$, a maximum a posteriori estimator yields (MAP derivation similar to [26]):

[TABLE]

In the case where atoms can be grouped with sufficient group population, the user may want to identify the columns of factor $\boldsymbol{B}$ with a group of atoms rather than with a single element of $\boldsymbol{D}$. This is the method used in [27] for hyperspectral imaging, where the dictionary contains multiple spectra for each class of mineral. In this scenario, a solution is to cluster the available atoms and use the centroids as columns of $\boldsymbol{D}$ along with inter-class covariance. Then comparing $\boldsymbol{B}$ to multiple classes of atoms with known averages and covariances amounts to penalizing using the Mahalanobis distance [27] in (9).

Dictionary-based PARAFAC2

There are other ways to account for discrepancy between the dictionary and the true phenomena underlying the data. In particular, it is not unreasonable that this discrepancy depends on the other experimental parameters (that is, the first or third mode in our case). For instance, in hyperspectral imaging, a modification of the geometry of the ground due to weather conditions can modify the spectra. If the data is collected along time, then the spectrum of a particular material can relate to one atom in the dictionary but with a time-dependent discrepancy. Since modeling variability is a topic-dependent task, we only provide one instance of a modified DCPD that handles this task, to serve as a guide for further works.

A well-known modified CPD model that accounts for variability is the PARAFAC2 model [28, 29]. Here we show that the PARAFAC2 model can be adapted incorporating dictionary information. The model is the following: for all positive integer $m$ smaller than $M$,

[TABLE]

where $\boldsymbol{M}_{m}\in\mathds{R}^{K\times L}$ is the $m^{th}$ slice of data tensor $\boldsymbol{\mathcal{T}}$ along the third mode, $\boldsymbol{c}_{m}$ is the $m^{th}$ row of a matrix of parameters $\boldsymbol{C}$, matrices $\boldsymbol{P}_{m}$ with orthogonal columns are unknown and $\boldsymbol{Q}$ is an unknown latent factor related to all the factors $\boldsymbol{B}_{m}$. Dictionary $\boldsymbol{D}$ has sizes $D\times d$ where $D$ can be different from $M$.

Here the three-way data is seen as a collection of matrices with shared first mode factor $\boldsymbol{A}$ and similar factors $\boldsymbol{B}_{m}=\boldsymbol{P}_{m}\boldsymbol{Q}$. The matrix $\boldsymbol{Q}$ stands for a latent shared factor among matrices $\boldsymbol{B}_{m}$. To include additional knowledge provided by the dictionary, is it assumed that this $\boldsymbol{Q}$ follows equation (5). Because the columns of $\boldsymbol{P}$ are orthogonal, the underlying hypothesis here is that all $\boldsymbol{B}_{m}$ have the same covariance matrix, which is a relaxed assumption with respect to the underlying hypothesis in (6), that is, $\boldsymbol{B}_{m}=\boldsymbol{Q}$ for all $m$.

An hidden advantage of the PARAFAC2 version of DCPD is that row dimension of the dictionary can be different from the number of samples on the second mode, since $\boldsymbol{P}_{m}$ does not need to be square. In practice this means that when using the PARAFAC2 dictionary model, the measured spectral bands can be different than the spectral bands at which atoms are sampled.

II-C Self-Dictionnary

If no dictionary is available on any of the factor, it is still possible in some specific cases to obtain the dictionary from the data. The topic of learning a dictionary from auxiliary data will be dealt with in future works. Rather, we focus here on the separability assumption [30, 31], that is when the data itself can be used as a dictionary. In the matrix case, the separability has been well studied in the context of nonnegative matrix factorization (NMF), see next section for details. With the suggested formalism, the self-dictionary CPD for matrices is written as follows:

[TABLE]

where   $\|\boldsymbol{s}_{i}\|_{0}=1\text{ for }i\in\{1,\dots,R\},~{}\boldsymbol{S}\in\{0,1\}^{d\times R}$. A working hypothesis of this model for matrices is that the second mode factor is a subset of the rows of $\boldsymbol{M}$.

A generalization of the separability assumption to higher order tensors is however not straightforward. Recently, an attempt was made to define “pure slices” [32], but another possible generalization is obtained by supposing the columns of factor $\boldsymbol{B}$ are contained in all the $K\times M$ columns of the unfolding matrix $\boldsymbol{T}_{2}$. The drawback of this model is the possibly very large number of correlated atoms in the obtained dictionary $\boldsymbol{T}_{2}$. There are various approaches to extend the separability assumption to high-order tensors, but this paper only deals with the case where $\boldsymbol{D}=\boldsymbol{T}_{2}$. Considering other variants is a direction for further research.

III Related works

III-A Sparse coding

Dictionary CPD presented in this paper can be seen as an extension of sparse coding for higher order tensors, see for instance [8, 9]. In sparse coding, a dictionary $\boldsymbol{D}$ is available, and the data is to be expressed as linear combinations of a small number of $R$ atoms:

[TABLE]

where $\boldsymbol{X}$ is imposed to have at most $R$ non-zero rows. This constraint can be formulated using the $\ell_{0}/\ell_{q}$ pseudo norm of $\boldsymbol{X}^{T}$ for any $q\geq 0$, where the $\ell_{p}/\ell_{q}$ norm of a matrix is defined as

[TABLE]

In fact, $\boldsymbol{X}$ has at most $R$ non-zero rows if and only if $||\boldsymbol{X}^{T}||_{\ell_{0}/\ell_{q}}\leq R$.

This model is equivalent to dictionary CPD for matrices by setting $\boldsymbol{S}\boldsymbol{B}^{T}=\boldsymbol{X}$. Nevertheless there is a crucial difference between the two formulations:

• •

By splitting the variable $\boldsymbol{X}$ into two variables $\boldsymbol{S}\boldsymbol{B}^{T}$, the linear constraint $\boldsymbol{M}=\boldsymbol{D}\boldsymbol{X}$ becomes non-convex with respect to $\boldsymbol{S}$ or $\boldsymbol{B}$.

• •

The sparse coding formulation involves a large number of variables, $d\times m$ if $\boldsymbol{M}$ is $n$ by $m$, whereas the DCPD formulation for matrices involves $(d+m)\times R$ parameters.

In summary, sparse coding involves much more parameters but yields a convex parametrization of the data.

In practice, because of the presence of noice, the constraint $\boldsymbol{M}=\boldsymbol{D}\boldsymbol{X}$ is replaced by the minimization of $||\boldsymbol{M}-\boldsymbol{D}\boldsymbol{X}||$ (which is convex for any norm $||.||$). Since the mixed pseudo norm $\ell_{0}/\ell_{q}$ is non-convex, most methods for computing sparse coding rely on convex relaxations. Row-sparsity of $\boldsymbol{X}$ can be achieved in different ways; in particular using convexifications based on the $\ell_{1}$ norm, e.g., $\ell_{1}/\ell_{2}$ [9, 33], $\ell_{1}/\ell_{\infty}$ [34], or using more sophisticated models [35, 36]. As far as we know, it is the first time sparse coding is tackled using the reformulation $\boldsymbol{X}=\boldsymbol{S}\boldsymbol{B}^{T}$. Although the resulting problem cannot be relaxed easily, it involves significantly fewer variables hence will be applicable to large-scale problems.

Recently, Salhoun et al. suggested a direct extension of sparse coding for higher order tensors [37] in the particular context of harmonic retrievals, but an obvious difficulty is that the row-sparsity has to be imposed on Khatri-Rao products of factors.

III-B NMF with self dictionary in spectral unmixing

Similarly, self dictionaries have been studied for matrices using the sparse coding formalism, by setting $\boldsymbol{D}=\boldsymbol{M}$. This is often called the pure-pixel assumption in spectral unmixing [24], since the self-dictionary model assumes some columns of the matrix $\boldsymbol{M}$ are not mixtures of more than one column of $\boldsymbol{B}$. As far as we know, there are mainly two types of approaches to tackle (11):

• •

Geometric approaches that selects the atoms in the dictionary based on some geometric criteria, typically based on the volume of the convex hull of $\boldsymbol{M}^{T}\boldsymbol{S}$. These approaches include for example vertex component analysis (VCA) [11] and the successive projection algorithm (SPA) [38, 39, 40, 13]. They are usually fast, running in $\mathcal{O}(mnR)$ operations. However, they do not always select atoms leading to a small data fitting term $||\boldsymbol{M}-\boldsymbol{A}(\boldsymbol{M}^{T}\boldsymbol{S})^{T}||_{F}$ since, most of them do not take it into account directly, as they usually put an emphasis on some geometric properties of $\boldsymbol{M}^{T}\boldsymbol{S}$ (such as having a large volume). In particular, these methods are in general sensitive to outliers.

• •

Sparse regression approaches that are based on the sparse coding reformulation of (11)

[TABLE]

and achieve column sparsity constraints on the scores $\boldsymbol{X}$ in different ways [9, 34, 35, 36, 33]. These methods have the advantage to better model (11) than geometric approaches as they take into account the data fitting term explicitly. They usually provide good solutions but are rather costly as an optimization problem in $dn$ variables must be solved, where $n$ stands for the column dimension of $\boldsymbol{M}$. In particular, since here $\boldsymbol{D}=\boldsymbol{M}^{T}$, we have $d=n$ hence $n^{2}$ variables. In hyperspectral unmixing, $n$ is the order of millions and these approaches are impractical. Hence pixels have to be selected in a preprocessing step [41] (e.g., using a geometric approach). Moreover, the problem solved is an approximation of the original problem, which results may not be as close as desired to the solutions of the non-convex problem.

In section V, we describe several algorithms to tackle the proposed formulation (11). They combine the advantages of the two types of approaches described above: they are fast, running in $\mathcal{O}(mnR)$ operations, but taking explicitly the data fitting term $||\boldsymbol{M}-\boldsymbol{X}\boldsymbol{M}||_{F}^{2}$ into account.

III-C Constrained tensor decompositions

The DCPD can be understood as a constrained CPD model. It is similar in spirit to computing CPD when a basis of representation is known for one of the factors, which typically happens when using Tucker Decomposition [42] or a basis of splines [43] for compression. Other linearly constrained tensor decomposition model are obtained when the components are linearly dependent [44], or when the factors are stuctured, e.g. Hankel or Toeplitz matrices [45].

From a constrained tensor decomposition perspective, the novelty of the present work is that the dictionary is overcomplete and is therefore not a basis of the factor space. Thus sparsity constraints are imposed on the coefficients and our approach is rather combinatorial. We have already published some preliminary results focusing on matrices and applications to spectral unmixing[46].

III-D Parameterized factors in tensor decomposition

Using a dictionary to help with recovering factor $\boldsymbol{B}$ of the CPD is also closely related to parameterizing the columns of that factor. Parametrization is a viable option when an analytical formulation of the atoms of the dictionary is possible. The dictionary is then a continuous dictionary, see for instance [47, section IV] and references therein. This continuous parameterization may yet not always be achievable in real-life applications. The motivation behind the two methods is however the same: reduce the number of degrees of freedom in the tensor decomposition model by providing a set of admissible solutions to improve estimation accuracy, and to restore identifiability in some pathological cases, some of which are discussed in the next section.

III-E Abundances estimation in hyperspectral unmixing

In the context of spectral unmixing of hyperspectral images, using a known library of spectra to estimate the second factor $\boldsymbol{B}$ is a widely studied topic. Matrix $\boldsymbol{A}$ contains the abundances and refers to relative concentrations of materials on each pixel. Except for the methods described above, the most widely used techniques to compute $\boldsymbol{A}$ when $\boldsymbol{D}$ is known are, to the best of our knowledge, MESMA [48], MELSUM [49], BSMA [50] and AutoMCU [51], which all rely on a more or less exhaustive search of all combinations of atoms in $\boldsymbol{D}$, except AutoMCU which only draws randomly a subset of possible atoms. For all possible combinations, abundances are computed, and the best abundances are those that minimize reconstruction error in addition to satisfying some interpretability criteria.

The proposed approach differs from these techniques since it merges blind source separation techniques for spectral unmixing, namely NMF, but featuring atom selection. To the best of our knowledge, our approach has not been described yet in the spectral unmixing literature, and is bound to be computationally less expensive than an exhaustive search.

IV Identifiability

The following section contains partial results on the identifiability of the DCPD model, cast in particular cases of interest. In propositions 1, 3 and 4 below, the assumption that no atoms are picked twice is made. This assumption is necessary for proving the identifiability of the DCPD parameters, but should not be necessarily imposed in the higher-order case. For instance, time series of hyperspectral images may require multiple abundance and time components for a single material when decomposed with the CPD model.

IV-A Matrix case

It is well known that for matrices, the low rank CPD model is not identifiable because of the rotation ambiguity, see equation (4). However, when a dictionary is available for one of the modes, this rotation ambiguity may be fixed given some conditions on $\boldsymbol{D}$:

Proposition 1**.**

*Let $\boldsymbol{M}$ be a real $n\times m$ matrix of rank $R$, and let $\boldsymbol{D}$ be a real $n\times d$ matrix with $\text{spark}(\boldsymbol{D})>R$ where $\text{spark}(\boldsymbol{D})$ is the minimum integer $k$ such that at least one subset of $k$ columns of $\boldsymbol{D}$ is rank-defficient. If there exist a full column-rank $\boldsymbol{S}\in\{0,1\}^{d\times R}$ with column sparsity set to 1, and $\boldsymbol{A}\in\mathds{R}^{n\times R}$ with nonzero columns such that $\boldsymbol{M}=\boldsymbol{A}(\boldsymbol{D}\boldsymbol{S})^{T}$, then $\boldsymbol{S}$ and $\boldsymbol{A}$ are unique up to permutation ambiguity. *

Proof.

Since $\text{spark}(\boldsymbol{D})>R$ and $\boldsymbol{S}$ is full column rank, $\boldsymbol{B}=\boldsymbol{D}\boldsymbol{S}$ has full column rank hence $\boldsymbol{A}$ is unique up to permutations in the decomposition $\boldsymbol{M}=\boldsymbol{A}\boldsymbol{B}^{T}$ if $\boldsymbol{B}$ is unique up to permutations. Moreover, given $\boldsymbol{M}=\boldsymbol{A}\boldsymbol{B}^{T}$, because $\boldsymbol{A}$ has no zero columns, the column space of $\boldsymbol{B}$ and the row space of $\boldsymbol{M}$ are equal. Such a $\boldsymbol{B}$ is built by selecting $R$ atoms in $\boldsymbol{D}$ that both span and belong to the row space of $\boldsymbol{M}$. Because the span of $\boldsymbol{D}$ is strictly larger than $R$, there is only one such set of $R$ atoms. Thus $\boldsymbol{B}$ is unique up to permutation. ∎

IV-B CPD and DCPD uniqueness

For all applications where the identifiability of the CPD is usually verified, a natural question to ask is whether the DCPD will automatically be identifiable or if some additional conditions need to be checked before trying to use the DCPD model. It turns out that given the uniqueness of the CPD up to scaling and permutations of the factors, the only requirement to obtain uniqueness of the DPCD is the uniqueness of the factorization $\boldsymbol{B}=\boldsymbol{D}\boldsymbol{S}$, which itself is quite simple to check.

Proposition 2**.**

Let $\boldsymbol{\mathcal{T}}\in\mathds{R}^{K\times L\times M}$ be such that $\boldsymbol{\mathcal{T}}$ admits a unique rank $R$ CPD up to scaling and permutations. Then if $\boldsymbol{D}$ does not contain collinear atoms and there exist $\boldsymbol{S}$ verifying (5), the DCPD is also unique up to scaling and permutations.

Proof.

Because the CPD of $\boldsymbol{\mathcal{T}}=\left(\boldsymbol{A}\mathop{\otimes}\boldsymbol{B}\mathop{\otimes}\boldsymbol{C}\right)\boldsymbol{\mathcal{I}}$ is unique, the DCPD is unique if and only if $\Omega=\{\boldsymbol{S}|\boldsymbol{B}=\boldsymbol{D}\boldsymbol{S}\}$ is a singleton. The set $\Omega$ is not empty by assumption, and if $\boldsymbol{S}_{1},\boldsymbol{S}_{2}$ belong to $\Omega$, then the columns they select in $\boldsymbol{D}$ have to be collinear (or equal if there was no scaling ambiguity on $\boldsymbol{B}$). By hypothesis this implies $\boldsymbol{S}_{1}=\boldsymbol{S}_{2}\Pi$, where $\Pi$ is a column permutation matrix. ∎

The existence of $\boldsymbol{S}$ is what is really difficult to asses in practice. On the other hand, the dictionary may contain very correlated atoms, but not exactly collinear atoms, so that the DCPD as a model is identifiable in practice whenever the CPD is identifiable.

Moreover, using proposition 1, it is possible to derive a mild sufficient condition for uniqueness of DCPD.

Proposition 3**.**

Let $\boldsymbol{\mathcal{T}}$ be a real $K\times L\times M$ tensor of rank $R$ and $\boldsymbol{D}$ a real $L\times d$ matrix. If there exist a full column-rank $\boldsymbol{S}\in\{0,1\}^{d\times R}$ with column sparsity set to 1, and $\boldsymbol{A}\in\mathds{R}^{K\times R}$ and $\boldsymbol{C}\in\mathds{R}^{M\times R}$ such that $\boldsymbol{\mathcal{T}}=\left(\boldsymbol{A}\mathop{\otimes}\boldsymbol{D}\boldsymbol{S}\mathop{\otimes}\boldsymbol{C}\right)\boldsymbol{I}_{R}$, if $\text{spark}(\boldsymbol{D})>R$ and if $\boldsymbol{A}\odot\boldsymbol{C}$ is full column-rank and $\boldsymbol{A}$ and $\boldsymbol{C}$ do not have zero columns, then $\boldsymbol{S}$, $\boldsymbol{A}$ and $\boldsymbol{C}$ are unique up to permutation and scaling ambiguity.

Proof.

By applying Proposition 1 to the rank $R$ matrix $\boldsymbol{T}_{2}=\boldsymbol{D}\boldsymbol{S}\left(\boldsymbol{A}\odot\boldsymbol{C}\right)^{T}$, $\boldsymbol{S}$ and the Khatri-Rao product $\boldsymbol{P}=\boldsymbol{A}\odot\boldsymbol{C}$ are unique up to permutation ambiguity. It remains to prove that the decomposition of $\boldsymbol{P}$ into $\boldsymbol{A}\odot\boldsymbol{C}$ is unique provided that the dimensions of the problem are fixed. Column-wise, we need to check that the factorization $\text{matr}(\boldsymbol{p}_{i})=\boldsymbol{a}_{i}\boldsymbol{c}_{i}^{T}$ is unique. This rank-one approximation is moreover up to scaling, and because this fixed matricization operator is an isomorphism, the decomposition of $\boldsymbol{P}$ is therefore unique. ∎

An important remark is that Proposition 3 is true for any order, since rank-one decompositions are always unique up to scaling. Similar discussion on identifiability when one factor is unique can be found in [52] and references therein. Moreover, if $\boldsymbol{A}$ has many collinear columns as may be the case with spectra in spectral unmixing, if $\boldsymbol{A}$ has no zero column and $\boldsymbol{C}$ is full column rank, then $\boldsymbol{A}\odot\boldsymbol{C}$ is full column-rank and the DCPD model does not suffer from the rotation ambiguity inherent to the CPD with collinear columns in factors. Colinear columns in factors appear when multiple columns of a factor in one mode are necessary to express the evolution of only one physically meaningful component.

IV-C Existence of the best low-rank DCPD

For third-order tensors, another important advantage of DCPD is that it makes the optimization problem well-posed. In fact, for CPD, the optimal solution may not exist as the feasible set is open; see for example [53].

Proposition 4**.**

Let $\text{spark}(\boldsymbol{D})>R$ and impose that $\boldsymbol{S}$ is full column rank, then the optimal solution of (7) is attained.

Below we provide two different proofs, each shedding a different light on the existence of the best low rank DCPD approximation.

First proof.

Let us show that the set $\mathcal{E}_{\boldsymbol{B}}=\{\boldsymbol{B}(\boldsymbol{A}\odot\boldsymbol{C})~{}|~{}\boldsymbol{A}\in\mathds{R}^{K\times R},\boldsymbol{C}\in\mathds{R}^{M\times R}\}$ is closed. First, notice that $\mathcal{E}_{\boldsymbol{B}}$ is the image of $\mathcal{E}=\{\boldsymbol{A}\odot\boldsymbol{C}~{}|~{}\boldsymbol{A}\in\mathds{R}^{K\times R},\boldsymbol{C}\in\mathds{R}^{M\times R}\}$ by a full column rank linear operator. Indeed, $\text{spark}(\boldsymbol{D})>R$. Then it is sufficient to show that $\mathcal{E}$ is closed, since $\mathcal{E}_{\boldsymbol{B}}=(\boldsymbol{B}^{\dagger})^{-1}\left(\mathcal{E}\right)\bigcap\text{col}(\boldsymbol{B})$, where † is the left pseudo-inverse and $\text{col}(\boldsymbol{B})$ is the column space of $\boldsymbol{B}$. Finally, since the set of rank-one matrices is closed, and the matrices in $\mathcal{E}$ are columnwise vectorized rank-one matrices, $\mathcal{E}$ is closed. ∎

Second proof.

Since $\text{spark}(\boldsymbol{D})>R$, for any full column rank $\boldsymbol{S}$, $\boldsymbol{B}=\boldsymbol{D}\boldsymbol{S}$ has rank $R$. Note that there are a finite number of such $\boldsymbol{S}$. It remains to show that for any $\boldsymbol{B}$ of rank $R$, the infimum of

[TABLE]

is attained. Without loss of generality, we can assume that $||\boldsymbol{a}_{r}||_{2}=1$ for $1\leq r\leq R$ by the scaling degree of freedom of each rank-one tensor $\boldsymbol{\mathcal{D}}_{r}=\boldsymbol{a}_{r}\mathop{\otimes}\boldsymbol{b}_{r}\mathop{\otimes}\boldsymbol{c}_{r}$. Using unfolding, the infimum of the above problem is attained if and only if the infimum of

[TABLE]

is attained. Moreover, since $\boldsymbol{C}=0$ is a feasible solution, we can add the constraint

[TABLE]

implying $\|\boldsymbol{C}\left(\boldsymbol{A}\odot\boldsymbol{B}\right)^{T}\|_{F}\leq 2\|\boldsymbol{T}_{3}\|_{F}$. For $\text{rank}(\boldsymbol{B})=R$ and $||\boldsymbol{a}_{r}||_{2}=1\forall r$, we can show that

[TABLE]

where $\sigma_{R}(\boldsymbol{B})$ is the $R$th singular value of $B$. In fact, denoting $\boldsymbol{a}^{j}$ the $j$th row of $\boldsymbol{A}$,

[TABLE]

Finally, $\|\boldsymbol{C}\|_{F}\leq 2\frac{\|\boldsymbol{T}_{3}\|_{F}}{\sigma_{\min}(\boldsymbol{B})}$ hence the feasible set can be reduced to a compact set hence the infimum is attained since the objective function is continuous and bounded below. ∎

V Decomposition Algorithms : Greedy and non-greedy approaches

In the literature of sparse approximation, two families of algorithms have been studied extensively: greedy approaches based on matching pursuit, and continuous approaches based on convex relaxations of the $\ell_{0}$ pseudo norm [54, 55]. In the same spirit, we develop in the next two sections the two same kind of algorithms to attack (6). We will compare these approaches in Section VI.

This paper does not explicitly discuss algorithms for computing the CPD itself. In DCPD, factor matrices $\boldsymbol{A}$ and $\boldsymbol{C}$ can be estimated using any off-the-shelf CPD algorithm assuming $\boldsymbol{B}$ is fixed. This can be done for instance using the alternating least squares procedure in the unconstrained case, or exact non-negative least squares for non-negative CPD [56, 57]. Therefore, below, only the estimation of $\boldsymbol{S}$ and $\boldsymbol{B}$ for fixed $\boldsymbol{A}$ and $\boldsymbol{C}$ is discussed222Algorithms introduced in this section are available at https://jeremy-e-cohen.jimdo.com/downloads/.

V-A Greedy algorithms : matching pursuit

Let us first provide an algorithm to compute the minimum (7) with respect to $\boldsymbol{S}$. Since the set of solution for $\boldsymbol{S}$ is discrete, the underlying optimization problem is combinatorial. On the other hand, computing the unconstrained CPD can be done efficiently using alternating least squares (ALS). We want to take advantage of both unconstrained CPD and greedy algorithms for computation efficiency when estimating respectively factors $\boldsymbol{A},\boldsymbol{B},\boldsymbol{C}$ and the selection matrix $\boldsymbol{S}$.

This is what the matching pursuit-ALS (MPALS) described below does. The variable $\boldsymbol{B}$ is injected in the problem to be a proxy of $\boldsymbol{D}\boldsymbol{S}$. The matrix $\boldsymbol{B}$ is estimated using the least square update333The inverse is not actually computed, rather we solve the least square problem using any efficient solver. The bottleneck here is the large product $\boldsymbol{T}_{2}\left(\boldsymbol{A}\odot\boldsymbol{C}\right)$ since $R\ll\min(K,L,M)$.

[TABLE]

and $\widehat{\boldsymbol{S}}$ is then evaluated by choosing the closest atom in $\boldsymbol{D}$ up to a scaling factor. Finally $\boldsymbol{B}$ is reevaluated as $\widehat{\boldsymbol{B}}=\boldsymbol{D}\widehat{\boldsymbol{S}}$. In other words, $\boldsymbol{B}$ is estimated through a projected least squares update, where the projection space is the set spanned by $\boldsymbol{D}\boldsymbol{S}$ for all $\boldsymbol{S}$ with coefficients in $\{0,1\}$ and column sparsity set to 1.

Like most projected ALS algorithms such as the ALS algorithm for NMF444In Matlab, this is the ‘als’ algorithm of ‘nnmf’. [58], convergence of the global MPALS algorithm, that is including the alternating least squares estimation of $\boldsymbol{A}$ and $\boldsymbol{C}$, cannot be ensured since the cost function is not guaranteed to decrease at each step. Indeed, the least squares update of $\boldsymbol{B}$ decreases the cost function, but the projection step increases it. Therefore MPALS is bound to have few provable results in term of convergence, contrary to continuous algorithms presented below. It is however a very simple algorithm to implement with no parameter to tune and it provides good results in both simulated and real data experiments reported in Section VI. In particular, we observed convergence in practice, and the cost function decreases for most iterations in simulations. The complexity of $\mathcal{O}(RKLM+RLd)$ operations per iteration of MPALS inside an ALS algorithm is about the same as plain ALS if the dictionary is not excessively large.

Note that the term greedy is a bit abusive since an atom chosen to belong to $\boldsymbol{B}$ at some iteration of the global procedure may be discarded at a further iteration. What is greedy in MP-ALS is the procedure to choose atoms in the dictionary at each inner iteration. Also, MP-ALS can be easily adapted if the constraint on the number of elements in $\boldsymbol{S}$ is modified to allow for more than one atom to be used to approximate the columns of $\boldsymbol{B}$, thus the borrowed name matching pursuit.

MPALS can also easily be adapted to tackle other similar optimization problems. If the factors are constrained to be nonnegative, then the estimates of factors $\boldsymbol{A}$ and $\boldsymbol{C}$ can be obtained by nonnegative least squares as mentioned earlier, while factor $\boldsymbol{B}$ should be nonnegative since the dictionary should be nonnegative in this case. Also, if no column of $\boldsymbol{D}$ may be selected twice, then after computing the scores $\text{max}_{j}\frac{\langle\boldsymbol{B}_{i}|\boldsymbol{D}_{j}\rangle}{\|\boldsymbol{D}_{j}\|}$ for each column of $\boldsymbol{B}$, the atoms are assigned to each such column by solving an assignment problem in order to maximize the sum of the scores [59].

A smooth version of MPALS

MPALS features a projection of $\boldsymbol{B}$ on the atoms of the dictionary at each outer iteration, which is a very rough way to impose the sparsity constraint on the scores $\boldsymbol{S}$. To obtain a smoother optimization algorithm, we suggest to enforce sparsity constraints in a continuous manner using the flexible formulation of the DCPD for factor $\boldsymbol{B}$, while using the projected factor $\boldsymbol{D}\boldsymbol{S}$ when estimating the other factors. We call this algorithm smooth MPALS (SMPALS), it is summarized in Algorithm 2. This algorithm also finds the minimum of (9) solving the flexible dictionary problem with Gaussian noise. MPALS is modified in a straightforward manner, by making the least squares update of $\boldsymbol{B}$ depend on $\boldsymbol{S}$, and not reevaluate $\boldsymbol{B}$ after evaluating $\boldsymbol{S}$. This means removing $\boldsymbol{B}=\boldsymbol{D}\boldsymbol{S}$ in Algorithm 1, and replacing $\boldsymbol{B}$ first estimate with

[TABLE]

where $\lambda$ is a given parameter, set by the user at the beginning of the algorithm, meant to approach $\frac{1}{\sigma_{c}^{2}}$ in (9).

In the alternating outer loop, $\boldsymbol{B}$ can be set either to the exact $\boldsymbol{D}\boldsymbol{S}$, or to the approximate version computed by (15). The first choice provides the algorithm we called SMPALS, while the second is fully flexible, therefore called Flex-MPALS. For Flex-MPALS, $\lambda$ is also kept constant. In SMPALS, until $\|\boldsymbol{B}-\boldsymbol{D}\boldsymbol{S}\|_{F}^{2}\geq 0.01\|\boldsymbol{B}\|_{F}^{2}$ is reached, the coupling strength $\lambda$ increases by a multiplicative constant $p$. We chose a relatively aggressive choice for the increase $p=1.1$ for third-order tensors. For matrices, even a higher value of $p=1.5$ gave good results and allowed faster convergence.

Convergence of proposed algorithms

Although it is difficult to asses the convergence of MPALS, results for SMPALS and Flex-MPALS can be derived. The convergence of SMPALS is guaranteed if $\lambda$ is allowed to grow to infinity. Here is a sketch of the proof. The update of $B$ can be written as follows; see (15):

[TABLE]

This means that $\boldsymbol{B}$ gets arbitrarily close to $\boldsymbol{D}\boldsymbol{S}$ as $\lambda$ goes to infinity. This implies that, for $\lambda$ sufficiently large, $\boldsymbol{S}$ and $\boldsymbol{B}$ will no longer be modified: if $\boldsymbol{B}$ does not change sufficiently, $\boldsymbol{S}$ does not change because of the discrete nature of the problem. In the mean time, the updates of $\boldsymbol{A}$ and $\boldsymbol{C}$, that use $\boldsymbol{D}\boldsymbol{S}$, decrease the cost function and converge to a stationary point of the corresponding objective function (for $\boldsymbol{S}$ fixed) [60].

The convergence of Flex-MPALS, i.e. with a fixed $\lambda$ and no replacement $\boldsymbol{B}$ by $\boldsymbol{D}\boldsymbol{S}$, is also guarantied for a normalized dictionary. Indeed, Flex-MPALS is a block-coordinate descent algorithm where the blocks $\boldsymbol{A},\boldsymbol{B},\boldsymbol{C}$ are updated using an optimal least squares estimate, while the estimated $\boldsymbol{S}$ in Flex-MPALS minimizes $\|\boldsymbol{B}-\boldsymbol{D}\boldsymbol{S}\|_{F}^{2}$ if the atoms of the dictionary have unit norm. Therefore the cost function decreases at each step of Flex-MPALS hence convergence to some value (since it is bounded below by zero).

V-B Continuous approaches

As explained above, the original optimization problem underlying the DCPD model is combinatorial because of the $\ell_{0}$ pseudo-norm and the fact that $\boldsymbol{S}$ is binary. To develop continuous algorithms such as gradient descent, a first step is to derive a continuous relaxation for the DCPD formulation. The first relaxations that comes to mind when working on $\ell_{0}/\ell_{p}$ pseudo-norm is a mixed norm $\ell_{1}/\ell_{2}$ which encourages sparsity column-wise. However in the DCPD model, $\boldsymbol{S}$ can only have one non-zero coefficient per column, which makes a $\ell_{1}/\ell_{2}$ based optimization difficult to tune. Rather than mixed norm, we choose to use the $\ell_{1}$ norm as a sparsity enhancing penalization on all the entries of $\boldsymbol{S}$, but under the constraint that the $\ell_{2}$ norm of the columns of $\boldsymbol{S}$ are set to one (which is not convex). Indeed, if $\mathcal{B}(1,\ell_{2})$ is the unit ball of the $\ell_{2}$ norm, the solutions to

[TABLE]

are exactly located on the coordinate axes, which is also the constraint imposed on the columns of $\boldsymbol{S}$. In light of this remark, solely imposing a $\ell_{1}$ penalization on coefficients of $\boldsymbol{S}$, that could be easily done with for instance ADMM, may not be sufficient to ensure the high level of column sparsity that is sought.

We chose to use the fast gradient from Nesterov [61] which can be understood as a proximal gradient descent with averaged steps. Although our problem is not convex, this method can still be applied and has been shown to work well in the non-convex case [62]. Fast gradient is used to solve the following subproblem in $\boldsymbol{S}$

[TABLE]

for a given parameter $\delta$ fixed by the user, and where $\boldsymbol{s}_{i}$ is the $i^{th}$ column of $\boldsymbol{S}$. The non-negativity constraint on $\boldsymbol{S}$ makes $\gamma$ differentiable with respect to $\boldsymbol{S}$:

[TABLE]

and the gradients of (18) with respect to $\boldsymbol{A}$ and $\boldsymbol{C}$ can be found for instance in [63] (we set $\boldsymbol{B}=\boldsymbol{D}\boldsymbol{S}$) if an all-at-once optimization is sought. The normalization constraint on the columns of $\boldsymbol{S}$ is imposed by normalization of $\boldsymbol{S}$ at each iteration. To avoid all values of a column of $\boldsymbol{S}$ to be non-positive, which would make the normalization meaningless, the gradient step is constrained as follows to ensure all columns of $\boldsymbol{S}$ after the gradient update have at least one positive entry:

[TABLE]

where $\epsilon_{S}$ is the Lipschitz constant and $g_{S}$ is the gradient of the cost function with respect to $\boldsymbol{S}$. $M(i,j)$ refers to the entry of matrix $\boldsymbol{M}$ indexed by $(i,j)$.

The global optimization procedure we suggest is therefore a mixture of ALS and fast gradient, denoted ALS-FG and summarized in Algorithm 3.

An important remark is that this algorithm is to be used inside an outer loop consisting of least squares problems. To ensure a smooth transition from the unconstrained problem to enforcing the dictionary, the penalty coefficient $\delta$ is linearly increased up from zero to a maximum value at the last outer iteration specified by the user. We observed that using the fast gradient improves convergence speed with respect to a simple coordinate descent, but other methods can also be implemented to solve the subproblem. Again, a flexible version of ALS-FG tackling optimization problem (9) can be derived from Algorithm 3 by adapting the gradient of $\gamma$ with respect to $\boldsymbol{S}$ and using $\boldsymbol{B}$ as a variable.

Stopping criterion

Choosing the stopping criterion is application-dependent, and we let this parameter be tuned by interested users. Nevertheless, a baseline that we used in simulations below is to compute the residual error $E_{i}$ at iteration $i$ of the global optimization procedure, every few iterations, and to check whether the error is still decreasing enough by computing $\frac{|E_{i}-E_{i-1}|}{E_{i}}$. We set the number of iterations of the fast gradient inside ALS-FG to 10.

Normalization

An inherent ambiguity of the CPD model which is also present in the DCPD model is the scaling ambiguity, that is, the norm of columns of the factors is not determined solely by the model. It is possible to fix the norm of factors in both MPALS and ALS-FG by normalizing the columns of factor $\boldsymbol{A}$ after each update of this factor, but this is not mandatory.

Initialization

Because the global optimization problems underlying the DCPD model and its flexible counterparts are non-convex and the proposed optimization algorithms work locally, it is a priori crucial to use a good initial point. If any information is available on the factors like non-negativity, it should be used in the initialization procedure. However, in the general case where no such information is available, one possible strategy is to compute an unconstrained CPD and use the obtained factors as an initialization for MPALS and ALS-FG. In the experiments below, we observed that indeed MPALS is highly sensitive to initialization. In particular, using a single random initialization leads in many cases to poor results for third-order tensors. Note however that is worked well for hyperspectral images.

VI Experiments on simulated and hyperspectral data

VI-A Simulated data experiments

VI-A1 Methodology

A critical question yet unanswered is whether using a dictionary within the tensor factorization model improves on identification error, which is the percentage of columns of $\boldsymbol{B}$ correctly matched to atoms of $\boldsymbol{D}$, with respect to only projecting the result of an unconstrained factorization on the set of atoms. Because information is added about factor $\boldsymbol{B}$ in the model, the identification error is expected to be smaller by using DCPD, in particular if the simulations are made challenging in terms of signal to noise ratio (SNR), rank under-estimation and over-estimation, or conditioning of the factors.

If the SNR is high, if the tensor is well-conditioned and if the rank is known, unconstrained CPD may provide relatively low identification error because of the generic uniqueness of the solution. For this reason, we choose difficult scenarios to compare DCPD to unconstrained CPD with outputs identified using a projection on the set of atoms. The dictionary used is constituted of atoms as follows:

[TABLE]

where $\boldsymbol{u}=\{1,\dots,L\}$, $a_{k},b_{k},\nu_{k}^{i},\mu_{k}^{i},e_{k}^{i},f_{k}^{i}$ are realizations of uniformly distributed random variables on respectively $[-1/L,1/L],[-1,1],[-1/4,1/4]$, $[-1/4,1/4]$, $\{1,\dots,L\}$, $\{1,\dots,L\}$ and tri is a triangular pulse function, of support $[-2,2]$. Atoms are then normalized with the $\ell_{2}$ norm. Index $k\in\{1,\dots,c\}$ is a class index, $c$ being the number of classes, and $i\in\{1,\dots,\frac{d}{c}\}$ indexes individual atoms in each class. This means that atoms of the dictionary are a sum of a linear baseline common to groups of atoms and of two individual features; see Figure 1 for an example. This design is meant to bluntly echo spectral signatures of materials in hyperspectral imaging, where only a small number of features discriminate materials from the same family. For the experiments below, a single realization of the dictionary was used, with no collinear columns (as to satisfy Proposition 2), but it could not be checked whether its spark satisfied the condition stated in Proposition 1. However, $\text{spark}(\boldsymbol{D})$ is equal to the dimension of the atoms plus one with probability one, which ensures that the parameters of the DCPD models used in these simulations are identifiable almost surely.

The data are generated as follows: each entry of the factors $\boldsymbol{A}$ and $\boldsymbol{C}$ is drawn according to a unitary centered Gaussian distribution, and the columns of $\boldsymbol{A}$ and $\boldsymbol{C}$ are then normalized with the $\ell_{2}$ norm, which makes these factors well-conditioned matrices. Matrix $\boldsymbol{S}$ is fixed in all experiments so that only one atom of each class is contained in the columns of $\boldsymbol{B}$. White Gaussian noise of fixed variance $\sigma$ is added to the data tensor.

Given the structure of the dictionary, the correlation among atoms of the same class can be quite high. The dictionary is constituted of a large number of atoms ($d=1000$) subdivided in $c=50$ groups. The dimensions of the tensors are $20\times 50\times 7$, yielding an inter-group correlation of atoms that reaches $0.999$ at most. We choose a noise variance of $\sigma=0.01$ which leads to an average SNR of about $11.5$dB, and set the rank to $R=10$.

As explained above, the simulations focus on the impact on identification performance when the rank is wrongly estimated, and when one of the factor matrix, here $\boldsymbol{C}$, is ill-conditioned. Therefore, we chose to grid over an estimated rank $Re$ ranging from $7$ to $13$ with a good conditioning of $\boldsymbol{C}$. The impact of the conditionning of $\boldsymbol{C}$ is studied through a grid on a parameter $\rho$ such that given a randomly-drawn well conditioned $\boldsymbol{C}^{(0)}$, it is modified as follows

[TABLE]

and is then normalized. When $\rho$ is equal to 0, $\boldsymbol{C}$ has column rank equal to 1, while its entries follow i.i.d. Gaussian distributions when $\rho$ is 1. The estimated rank is set equal to the true rank when studying the ill-conditioned case. The number of realization in each setting is set to $N=100$.

The CPD algorithm used as a baseline for comparison and for initialization is the NWAY toolbox [64]. We compare projected results of the NWAY toolbox with MPALS, SMPALS, Flex-MPALS and ALS-FG. Coupling strength in Flex-MPALS is fixed to $\lambda=0.04$ in this experiment. The value 0.04 provided spectra close to the dictionary, while allowing some flexibility. In practice, $\lambda$ could be tuned for example in order to achieve a given relative error between the columns of $B$ of the corresponding matched atoms of the dictionary (for example 1% relative error). The number of iterations is fixed to a maximum of 1000, except for the rand-MPALS (MPALS initialized randomly) where the maximum iteration number is 5000. Algorithm also stops if the stopping criterion mentioned in the previous section goes below $10^{-4}$.

Finally, in both cases we use the NWAY toolbox for initialization, except for the randomly initialized MPALS which is initialized with factors following the same distribution as the true factors.

VI-A2 Results

Figures 2 and 3 report the mean identification rate that measures the percentage of well-matched columns between $\boldsymbol{B}$ and $\widehat{\boldsymbol{B}}$. In the case where the rank and the estimated rank are miss-matched, the unmatched factors of $\boldsymbol{B}$ and $\widehat{\boldsymbol{B}}$ are considered miss-matched. This leads to a best possible identification rate given by $\frac{|Re-R|}{\max(Re,R))}$ and refered to as “oracle”. Figure 4 reports the relative mean square error on factor $\boldsymbol{B}$ given by $\mathds{E}\left[\frac{\|\boldsymbol{B}-\widehat{\boldsymbol{B}}\Pi\|_{F}^{2}}{\|\boldsymbol{B}\|_{F}^{2}}\right]$ where $\Pi$ is a permutation matrix computed to match the estimated factors to the true factors. Table III reports the mean run time of each algorithm in the first experiment when $R$ equals $Re$.

As expected the identification rate increases with respect to projected NWAY when using MPALS and its variants, with a gain ranging from a few percent to over ten percent. This is a direct but non-obvious consequence of the increased estimation performances. It can be seen that MPALS performs the best in almost all cases, especially when the factor $\boldsymbol{C}$ is ill-conditioned. However, the wrong estimation of the rank does not impact the gain in performance between projected CPD and DCPD. Notably, SMPALS performs slightly worse than the naive MPALS.

An important result to observe is the steady performance of the randomly initialized MPALS over various values of $\rho$, both in terms of mean square error and identification rate. Because DCPD is supposed to be identifiable even when $\rho$ tends to [math], the performance of MPALS should not depend too much on it, as observed with the random initialization. However, when the unconstrained CPD model is used as an initialization method, the results of MPALS do depend heavily on $\rho$. This shows that MPALS and all other proposed algorithm are very sensitive to initialization, a fact also supported by the experiment in the next section. As a consequence, a good initialization method is crucial for the performance of MPALS, this is a topic for further research.

Even though the identification rates evolve similarly for all methods in the ill-conditioning experiment, the mean square error on $\boldsymbol{B}$ tells a different story. Since many atoms in the dictionary are very correlated, it is reasonable to assume that the mean square error on $\boldsymbol{B}$ could be small even though all atoms are wrongly identified. It can be observed on Figure 4 that on average, all proposed algorithms for computing the DCPD significantly outperform the projected ALS method when $\boldsymbol{C}$ is ill-conditioned. As a conclusion, the atoms picked by MPALS and its variants are on average much closer to the true atoms than the ones picked using a projected ALS method. This is especially true for the randomly initialized MPALS.

Finally, from this experiment, we observe that computing the DCPD with a combinatorial greedy approach performs better than a continuous approach, since the performance of ALS-FG is relatively poor. We believe the reason for this is that the greedy approaches can escape the basin of attraction of a local minimum during the projection step, which a (standard) continuous approach will not be able to do. In fact, we observed that, in most cases, once an entry is set to zero by ALS-FG, it remains zero in the course of the iterations, which is not the case for MPALS and its variants.

VI-B Spectral Unmixing with the self-dictionary model

For this communication on dictionary based tensor factorization models, we choose to try out the different models and algorithms on a well-known sparse coding problem, namely spectral unmixing under the pure-pixel assumption. As explained in Section II, the data itself can be used as the dictionary, and models described in this paper for matrices can be straightforwardly compared with state-of-the-art sparse coding approaches for spectral unmixing, namely the succesive projection algorithm (SPA) [38], the successive nonnegative projection algorithm (SNPA) [10], the Hierarchical Clustering algorithm (H2NMF) [65] and FGNSR [41]. Below, the spectra contained in the data are used as atoms, so that the dictionary has row dimension equals to the number of spectral bands ($\approx 150$), and number of atoms equals to the number of pixels in the hyperspectral images (HSI) ($\approx 10^{5}$).

The two data sets that will be used to compare DCPD with these methods are Urban and Terrain. These two HSI satisfy approximately the pure pixel assumption. The dimensions of the HSI after vectorization of the pixel dimensions are respectively $307^{2}\times 162$ and $(500\times 307)\times 166$. We chose to decompose the two data sets with respective ranks $R_{urban}=6$ and $R_{terrain}=5$ according to what has been done previously in the literature; see [41] and the references therein. The maximal number of iterations of the DCPD algorithms is set to 50. Also because it does not perform well, results for the continuous fast gradient algorithm are not presented below.

We initialize the MPALS algorithm and its variants with each state-of-the-art algorithm with 10 additional A-HALS steps. For instance MPALS and SMPALS initialized with SPA are respectively denoted as d-SPA and ds-SPA. Results are presented in Table IV. Additionally, considering the performance of projected ALS in the multiway experiment above, the performance of projected NMF is also shown, for instance nmf-SPA refers to the A-HALS algorithm initialized with SPA, which endmembers output are projected onto the data points and which abundances are re-estimated using non-negative least squares.

The relative reconstruction error is used as a performance metric, since it is not possible to assess the identification performance without a ground-truth. Therefore, this experiment only studies the efficiency of the MPALS algorithm for minimizing the objective function.

In all cases, using either SMPALS or MPALS improves on the initial values for identified spectra. As shown already in the previous experiment, initialization plays an important role in the final reconstruction error, but at least using the DCPD model always refines the solutions, even when it is initially low like for FGNSR-100. Moreover, the DCPD algorithms are not excessively costly with respect to other state-of-the-art methods. There is no significant difference between MPALS and its smooth counterpart. Figure 5 shows the estimated abundances and spectra with the H2NMF-MPALS algorithm for both Urban and Terrain, and materials can be identified for each component by the user (no ground truth is available). Finally, Figure 6 shows the reconstruction error map on Urban and Terrain HSIs for the H2NMF-MPALS. Clearly the remaining error is not distributed as an i.i.d. Gaussian noise, which means the Frobenius norm used as the distance metric is not adapted. Also, most of the remaining error comes from rooftops and roads in the Urban HSI, a zone probably corrupted by large spectral variability.

Because Flex-MPALS does not require exactly that the pure pixel assumption is verified, we do not include it in this simulation. Indeed, using a reconstruction error criteria, Flex-MPALS would outperform the other methods, but that would not mean that obtained abundances and endmembers are better. On the other hand, the flexible model better tackles spectral variability since it can modify to some extent the spectra extracted from the pure pixels.555Interested readers will find the code for matrix Flex-MPALS online.

VII Conclusions

To jointly separate and identify sources using tensor canonical polyadic decomposition and a known dictionary, we introduced in this paper the DCPD model along with some flexible variants. Identifiability of the DCPD model parameters was discussed in both the matrix and the higher-order tensor cases. We proposed both greedy and continuous algorithms for DCPD and compared them on synthetic data sets and hyperspectral images. We observed that (i) the greedy algorithms provide in most cases better results, and that (ii) our DCPD formulation improves results with respect to standard CPD on the synthetic data sets, and with respect to spectral unmixing approaches based either on convex relaxation or geometric methods. A particularly promising direction for further research would be the design of efficient initialization schemes for the proposed greedy algorithms which are particularly sensitive to initialization.

An interesting continuation of this work would also be to compute the Cramér Rao bounds of the parameters in the DCPD to formally support what was presented in the simulation section. Moreover, in most application cases the dictionary is either unknown, or describes only a subset of the components. Thus a partial DCPD model should be investigated, as well as a dictionary learning scheme for tensors following the DCPD model.

Acknowledgment

The authors would like to thank Pierre Comon, Rodrigo Cabral Farias and Miguel A. Veganzones for early fruitful discussions on this work.