# Optimal estimation of the canonical polyadic decomposition from low-rank tensor trains

Clémence Prévost, Pierre Chainais

► To cite this version:

Clémence Prévost, Pierre Chainais. Optimal estimation of the canonical polyadic decomposition from low-rank tensor trains. Signal Processing, 2025, 234, pp.110001. 10.1016/j sigpro.2025.110001. hal-05479957

HAL Id: hal-05479957

https://hal.science/hal-05479957v1

Submitted on 27 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Optimal estimation of the canonical polyadic decomposition from low-rank tensor trains☆

Clémence Prévost<sup>ID*,1</sup>, Pierre Chainais

Univ. Lille, CNRS, Centrale Lille, UMR 9189 CRISTAL, F-59000, Lille, France

<figure><img src="image_2.png" /></figure>

## ARTICLE INFO

**MSC:**

15A23

15A69

**Keywords:**

Canonical polyadic decomposition

Tensor train decomposition

Low-rank models

Optimization

Factor retrieval

MIMO systems

## ABSTRACT

Tensor factorization has been steadily used to represent high-dimensional data. In particular, the canonical polyadic decomposition (CPD) is very appreciated for its remarkable uniqueness properties. However, computing the high-order CPD is challenging: numerical issues and high needs for storage and processing can make algorithms diverge. Furthermore, the recovery of the CP factors is an ill-posed problem. One way to circumvent this limitation is to exploit the equivalence between the CPD and the Tensor Train Decomposition (TTD). This paper formulates the CPD as a dimension reduction using a TTD followed by a global marginally convex optimization problem. This global optimization scheme estimates the CP factors with minimal error. The resulting approach, Dimensionality Reduction, joint Estimation of the Ambiguity Matrices and the CP FACTORS (DREAMFAC), relies on a block-coordinate descent that reaches a first-order stationary point when estimating the CP factors. DREAMFAC is also shown to be an optimal estimator that reaches the corresponding constrained Cramér-Rao bound. It therefore appears as a state-of-the-art solution to estimate the best rank-K CPD of a tensor (when it exists). Its performance is illustrated on the problem of parameter estimation in a dual-polarized MIMO system. Numerical experiments show the excellent practical performance of DREAMFAC, even with very low SNR.

## 1. Introduction

Tensors provide a faithful representation of higher-order observations by preserving their multidimensional structure [1,2]. For instance, color images can be seen as data cubes with two spatial dimensions (the pixels) and a spectral dimension (the colors) that encodes the red, green and blue channels. Low-rank tensor decompositions have been steadily used in the past decade to identify latent factors. They have been successfully applied to various problems in, e.g., signal processing and machine learning [3], brain signal processing [4], video completion [5] or telecommunications [6]. The notion of low rank for tensors is not unified: various tensor decompositions are available, coming with different definitions of the rank as well as different properties. The most popular generalization of the concept of rank from matrices to tensors is the canonical polyadic decomposition (CPD) [7]. Any tensor admits a polyadic decomposition, that is a linear combination of rank-one tensors. When the number of rank-one tensors necessary in this decomposition is minimal, it is called the tensor rank and the decomposition is called the CPD. The best rank-K CPD approximation of a

tensor has gained a lot of interest for data processing and analysis, due to its remarkable uniqueness properties up to scaling and permutation ambiguities [8]. However, retrieving this best rank-K approximation is an ill-posed problem [9]. In practice, this task is usually carried out by suboptimal iterative algorithms. Furthermore, the storage and processing cost of tensors increase exponentially with their order N. This limitation is known as the curse of the dimensionality [10]. To this respect, the Tensor Train (TT) decomposition (TTD) has been recently introduced [11] for dimension reduction. This decomposition has two advantages. First, it exploits a stable numerical estimation procedure which avoids the iterative algorithms used to compute the CPD. Second, it breaks the curse of dimensionality by operating on a set of matrices and third-order tensors, called TT-cores. In particular, its storage cost is linear in N. Numerous works were devoted to the task of dimensionality reduction in the CP model. In [12], the authors provided an equivalence between the CPD and the TTD through a set of non-singular change-of-basis matrices shared among the TT-cores. An algorithm, called Joint dimensionality Reduction And Factor

<footnote>☆ This work was partly supported by the ANR, France project "Chaire IA Sherlock" ANR-20-CHIA-0031-01 hold by P. Chainais, as well as by the national support within the programme d'investissements d'avenir, France ANR-16-IDEX-0004 ULNE and Région HDF, France.</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: clemence.prevost@univ-lille.fr (C. Prévost).</footnote>

<footnote n="1">¹ EURASIP member.</footnote>
Estimation (JIRAFE), was proposed, using the TT-SVD algorithm to estimate the TT-cores. It estimates CP factors corresponding to the TT-cores of a TT-SVD approximation of the initial high-order tensor. A JIRAFE-like procedure was also recently extended to constrained CPD in [13]. These algorithms were based on a series of local and sequential optimization problems instead of a global strategy. They permit to estimate the CP factors at a low computational cost, while only proposing an approximate solution to the estimation problem. JIRAFE estimates the change-of-basis matrices appearing in the TTD only once before propagation through the following optimization problems, without enforcing full coherence between the TT-cores. For this reason, JIRAFE performs a suboptimal estimation of the estimated CP factors, in particular for those estimated in last position of the tensor train. Its performance is highly sensitive to the initialization of the algorithm and may call for some post-processing to correct its defaults.

This work formulates the joint dimensionality reduction and estimation of the CP factors as a global optimization problem. Retrieving the best rank-K CPD approximation is transformed into a well-posed problem by estimating the TT-cores, at the cost of a small approximation error arising from the use of TT-SVD to first produce a best rank-K tensor train decomposition. Then, the equivalence between the CPD and the TTD is leveraged to retrieve CP factors that correspond to the best rank-K CPD approximation, provided that it exists. Although the resulting optimization problem is not convex, it is marginally convex with respect to each of the CP factors to estimate. Theoretical convergence guarantees to a first-order stationary point are available for this kind of problem in the framework of block-coordinate descent (BCD), see [14,15] and references therein. This motivates the introduction of a new algorithm for joint dimensionality reduction and estimation of the CP factors that exploits this strength. This new algorithm is called Dimensionality Reduction, joint Estimation of the Ambiguity Matrices and the CP FACTORS (DREAMFAC). In contrast with JIRAFE, it ensures the full coherence between TT-cores by considering coupled updates for the change-of-basis matrices. Therefore, DREAMFAC reaches a first-order stationary point for the estimation of the CP factors while JIRAFE does not. The optimality of DREAMFAC is confirmed by estimation theory using the Cramér–Rao Bounds (CRB) [16,17]. Cramér–Rao bounds for tensor CP models have been extensively studied in the literature, including performance bounds for uncoupled CP models [18–20], Bayesian frameworks [21], and constrained models [22]. For the considered global optimization model, the constrained Cramér–Rao bound (CCRB) can be used. Its versatility was shown by numerous works [23–27]. The experimental results confirm that the mean-squared error yielded by DREAMFAC reaches the CCRB, showcasing the efficiency of the proposed algorithm. An application to the problem of harmonic retrieval in MIMO systems [28] illustrates the optimality of DREAMFAC on a realistic low-rank decomposition problem. For modeling and estimation of MIMO channels, it is important to accurately estimate channels parameters at the base station (angles of arrival, angles of departure, path gains, and polarization parameters) to perform beamforming and deal with multiuser interferences. In [6], a tensor-based approach for dual-polarized MIMO channel estimation was proposed. The MIMO channel was recast as a fourth-order tensor admitting a CPD. The authors proposed an Alternating Least Squares (ALS) to estimate CP factors. In [29], the authors adapted the JIRAFE procedure to treat the MIMO channel estimation problem. Therefore this problem can be used as a relevant benchmark. The experimental validation highlights the expected optimality of DREAMFAC, whose performance is better than that of state-of-the-art methods.

The paper is organized as follows. Section 2 introduces important background on tensors and their low-rank decompositions. Section 3 describes the coupled optimization problem at hand and the proposed algorithm DREAMFAC. Section 4 contains the derivation of the Cramér–Rao bounds for the proposed coupled model. It describes the detailed probabilistic framework, and closed-form expressions for the matrices

to invert. Section 5 gathers an extensive set of numerical experiments to illustrate the optimality of DREAMFAC and its performance with respect to state of the art. They include simulations highlighting the efficiency of DREAMFAC compared to the CCRB, and experimental validation on the selected application to the estimation of MIMO channels.

**Notation.** The following notations [1,2] are used: lower (a) or uppercase (A) plain font for scalars, boldface lowercase (a) for vectors, boldface uppercase (A) for matrices and calligraphic (A) for tensors. The elements of vectors, matrices, and tensors are denoted by $a_d$, $A_{d_1,d_2}$ and $A_{d_1,...,d_N}$, respectively. The transpose of a matrix A is denoted by $A^T$. The matrix $I_N$ is the $N \times N$ identity matrix and $0_{L \times K}$ is the $L \times K$ matrix of zeros. The symbols $\boxtimes$, $\circled$, $\boxdot$ and $\otimes$ denote the Kronecker, Khatri–Rao, Hadamard (element-wise) and outer products. The operator $\mathrm{vec}\{\cdot\}$ stands for the standard column-major vectorization of a matrix or a tensor. The operation $\mathrm{Diag}\{\mathbf{A}, \mathbf{B}\}$ produces a block-diagonal matrix whose blocks are $\mathbf{A}$ and $\mathbf{B}$.

## 2. Background on low-rank tensor models

### 2.1. Preliminaries on tensors

A tensor $\mathcal{X} \in \mathbb{R}^{D_1 \times \cdots \times D_N}$ is an $N$-dimensional array indexed by the elements $\mathcal{X}_{d_1, \dots, d_N}$, for $d_n \in \{1, \dots, D_n\}$ ($n \in \{1, \dots, N\}$). Each dimension of a tensor is called a mode. A mode-$p$ fiber of $\mathcal{X}$ is a vector obtained by fixing all but the $p$th dimension.

**Definition 2.1 (Tensor Unfoldings).** The mode-$p$ unfolding of a tensor $\mathcal{X}$, denoted by $X^{(p)}$, is the matrix whose rows are the mode-$p$ fibers of $\mathcal{X}$, ordered according to the vectorization order. For a tensor $\mathcal{X} \in \mathbb{R}^{D_1 \times \cdots \times D_N}$, e.g., $X^{(1)} \in \mathbb{R}^{D_N \times D_2 \times D_1}$.

**Definition 2.2 (Matrix Mode Product).** The matrix mode-$p$ product between a tensor $\mathcal{X}$ and a matrix $\mathbf{M}$ is denoted by $\mathcal{X} \circled_{p} \mathbf{M}$ and is constructed such that each mode-$p$ fiber of $\mathcal{X}$ is multiplied by $\mathbf{M}$, e.g., the elements of the mode-1 product between $\mathcal{X} \in \mathbb{R}^{D_1 \times \cdots \times D_N}$ and $\mathbf{M} \in \mathbb{R}^{L \times D_1}$ are accessed as

$$ (\mathcal{X} \circled_{1} \mathbf{M})_{\ell, d_2, \dots, d_N} = \sum_{i=1}^{D_1} \mathcal{X}_{i, d_2, \dots, d_N} M_{\ell, i}, \quad \ell \in \{1, \dots, L\}. \quad (1) $$

Moreover, it holds that $\mathcal{Y} = \mathcal{X} \circled_{p} \mathbf{M} \Leftrightarrow \mathcal{Y}^{(p)} = \mathcal{X}^{(p)} \mathbf{M}^{\top}$.

**Definition 2.3 (Tensor Contraction Product).** The contraction product on modes $p, q$ between two tensors $\mathcal{X} \in \mathbb{R}^{D_1 \times \cdots \times D_N}$ and $\mathcal{Y} \in \mathbb{R}^{J_1 \times \cdots \times J_M}$ with $D_p = J_q$ is denoted by $\mathcal{X} \circled_{p}^{q} \mathcal{Y}$. It produces a tensor of order $N + M - 2$ such that

$$ \begin{aligned} & (\mathcal{X} \circled_{p}^{q} \mathcal{Y})_{d_1 \dots d_{p-1}, d_{p+1} \dots d_N, J_1 \dots J_{q-1}, J_{q+1} \dots J_M} \\ & = \sum_{\ell=1}^{D_p=J_q} (\mathcal{X})_{d_1 \dots d_{p-1}, \ell, d_{p+1} \dots d_N} (\mathcal{Y})_{J_1 \dots J_{q-1}, \ell, J_{q+1} \dots J_M}. \end{aligned} \quad (2) $$

**Definition 2.4 (Outer Product).** The outer product between $N$ vectors $a_n \in \mathbb{R}^{D_n}$ ($n \in \{1, \dots, N\}$) is a rank-one tensor $\mathcal{X} = a_1 \otimes \cdots \otimes a_N \in \mathbb{R}^{D_1 \times \cdots \times D_N}$ whose elements are accessed as $\mathcal{X}_{d_1, \dots, d_N} = (a_1)_{d_1} \times \cdots \times (a_N)_{d_N}$.

### 2.2. The canonical polyadic decomposition

For all $N$-th order tensor $\mathcal{X} \in \mathbb{R}^{D_1 \times \cdots \times D_N}$, there exists an integer $K$ such that it admits a canonical polyadic (CP) decomposition (CPD) as
$$ \mathcal{X} = [[\mathbf{A}_1, \dots, \mathbf{A}_N]], \quad (3) $$
where $A_n \in \mathbb{R}^{D_n \times K}$ ($n \in \{1, \dots, N\}$) are called the CP factors. When minimal, the integer $K$ is the rank of the tensor $\mathcal{X}$. Then each entry of $\mathcal{X}$ can be expressed as

$$ \mathcal{X}_{d_1, \dots, d_N} = \sum_{k=1}^{K} (A_1)_{d_1, k} \times \dots \times (A_N)_{d_N, k}. \quad (4) $$

The CP factors are essentially unique up to scaling and permutation ambiguities, if the rank $K$ is not too large [1,2]. The permutation ambiguity means that the columns of the latent CP factors can be reordered arbitrarily by any permutation matrix $\Pi \in \mathbb{R}^{K \times K}$ as

$$ \mathcal{X} = [[A_1, \dots, A_N]] = [[A_1 \Pi, \dots, A_N \Pi]]. \quad (5) $$

The scaling ambiguity means that the individual factors can be scaled as

$$ \mathcal{X}_{d_1, \dots, d_N} = \sum_{k=1}^{K} \left( \lambda_{1,k} (A_1)_{d_1,k} \right) \times \dots \times \left( \lambda_{N,k} (A_N)_{d_N,k} \right). \quad (6) $$

where $\lambda_{1,k} \times \dots \times \lambda_{N,k} = 1$, for all $k \in \{1, \dots, K\}$.

### 2.3. The tensor train decomposition (TTD)

The tensor train decomposition (TTD) [11] factorizes a tensor $\mathcal{X} \in \mathbb{R}^{D_1 \times \dots \times D_N}$ as a series of contractions between third-order tensors and matrices as

$$ \mathcal{X} = G_1 \frac{1}{2} \mathcal{G}_2 \frac{1}{3} \dots \frac{1}{N-1} \mathcal{G}_{N-1} \frac{1}{N} \mathcal{G}_N, \quad (7) $$

where $G_1 \in \mathbb{R}^{D_1 \times K_1}$, $G_N \in \mathbb{R}^{K_N-1 \times D_N}$, and $G_n \in \mathbb{R}^{K_{n-1} \times D_n \times K_n}$ (for $n \in \{2, \dots, N_1\}$), are referred to as TT-cores. The integers $K_1, \dots, K_{N-1}$ are called the TT-ranks.

A fast and efficient way to estimate the TT-cores is to resort to the TT-SVD algorithm [11], a procedure that sequentially extracts dominant singular vectors from tensor unfoldings. As a result, the TTD in (7) is not unique. In fact, due to the use of the SVD, we can replace two successive TT-cores $G_n$ and $G_{n+1}$ by $G'_n$ and $G'_{n+1}$ such that

$$ G'_n = G_n \frac{1}{3} M_n^{-1}, \quad G'_{n+1} = M_n \frac{1}{2} G_{n+1}, \quad (8) $$

where $M_n \in \mathbb{R}^{K_n \times K_n}$ is a non-singular change-of-basis matrix. This means that the multiplicative ambiguities in the TTD correspond to post- and pre-multiplications by nonsingular matrices.

### 2.4. Equivalence between the CPD and the TTD

There exists a theoretical equivalence between the CP and TT decompositions. If the best rank-$K$ CPD exists, then TT-SVD recovers the best TTD, termed $\text{TTD}_{\text{best}}$. This best TTD is in one-to-one correspondence with the best rank-$K$ CPD. With TT-ranks $K_1 = \dots = K_{N-1} = K \le \min(D_1, \dots, D_N)$, the link between the TTD and the best rank-$K$ CPD is such that [12]

$$ G_1 = A_1 M_1^{-1}, \quad G_n = [[M_{n-1}, A_n, M_n^{-T}], \quad G_N = M_{N-1} A_N^T, \quad (9) $$

where $M_n \in \mathbb{R}^{K_n \times K_n}$ are non-singular change-of-basis matrices. These matrices can be related to the permutation and scaling ambiguity matrices of the CPD, see the proof of [12, Theorem 6] for more details. If the best rank-$K$ CPD does not exist, TT-SVD will still yield an approximated TTD, whose distance to $\mathcal{X}$ in the Frobenius norm is less or equal to $\sqrt{N-1}\|\mathcal{X} - \text{TTD}_{\text{best}}\|_F$ [11]. The corresponding CPD inherits the same quality of approximation.

In practice, in many cases the best rank-$K$ CPD exists and therefore TT-SVD recovers the best TTD and the desired best rank-$K$ CPD. Instead of solving an ill-posed problem, a well-posed problem (estimating the TT-cores) is solved, at the cost of a small estimation error arising from the use of the TT-SVD algorithm. In the worst case where there is no best rank-$K$ CPD, the proposed pipeline is still expected to yield a CPD, but with the largest approximation error of $\sqrt{N-1}\|\mathcal{X} - \text{TTD}_{\text{best}}\|_F$.

## 3. Joint estimation and dimensionality reduction

This section introduces the problem of estimating the CP factors and its limitations. Assuming that a best rank-$K$ CPD exists, the equivalence between the CPD and the TTD is leveraged to estimate the CP-factors by solving a well-posed problem. The global optimization problem that corresponds to this task is introduced. Then, a method that performs optimal estimation of the CP factors using a block-coordinate descent (BCD) scheme is proposed. Finally, the proposed method is discussed with respect to the state-of-the-art, JIRAFE. We emphasize that this is a major contribution with respect to the state-of-the-art which is only suboptimal, contrary to the proposed approach that uses BCD.

### 3.1. Estimation of the CP factors: interest and limitations

Let $\mathcal{X}$ be a tensor of rank $K$ that admits a CPD of the form

$$ \mathcal{X} = [[A_1, \dots, A_N]], \quad \text{where } A_n \in \mathbb{R}^{D_n \times K} \quad \forall n \in \{1, \dots, N\}. \quad (10) $$

Let $\mathcal{Y} \in \mathbb{R}^{D_1 \times \dots \times D_N}$ be the noisy observation according to the following model

$$ \mathcal{Y} = \mathcal{X} + \mathcal{E}, \quad (11) $$

where the tensor $\mathcal{E}$ represents independent and identically distributed (i.i.d.) white Gaussian noise added to each entry of the tensor $\mathcal{X}$. The model (11) can be rewritten

$$ \mathcal{Y} = [[A_1, \dots, A_N]] + \mathcal{E}. \quad (12) $$

The problem of estimating the CP factors $A_n$ based on model (12) has found numerous interests. It has been applied, among others, to problems in signal processing and machine learning [3], brain signal processing [4], video completion [5] or telecommunications [6]. The popularity of the CPD lies in its high regularization power, mild uniqueness properties and versatility to model a variety of problems, from denoising to component analysis.

The most popular way to estimate the CP factors is an iterative Alternating Least-Squared (ALS) algorithm [30]. However, this procedure becomes exponentially expensive as the order $N$ increases. One iteration of CP-ALS requires $N \cdot \mathcal{O}(K^2 \prod_n D_n)$ flops. Note that this complexity can be reduced by resorting to stochastic-type optimization algorithms, see e.g., [31–33]. Nevertheless, estimators of the CP factors become less robust and are prone to convergence issues for high-order tensors. These limitations fall under the so-called curse of dimensionality [10]. Consequently, one may need to perform tensor dimensionality reduction to exploit the benefits of model (12).

### 3.2. Estimation of the CP factors by the TTD and coupled updates

For the sake of easier reading, the following notation will be adopted: $A = \{A_1, \dots, A_N\}$, $M = \{M_1, \dots, M_{N-1}\}$ and $G = \{G_1, G_2, \dots, G_{N-1}, G_N\}$. Estimating the CP factors by directly solving (12) using high-dimensional CP-ALS is computationally expensive. To circumvent this limitation, it is possible to perform low-rank denoising by first resorting to the TTD to perform dimensionality reduction.

A low rank-$K$ factorized approximation is obtained by solving

$$ \min_{G} \|\mathcal{Y} - \mathbf{G}_1 \frac{1}{2} \mathcal{G}_2 \dots \frac{1}{N-1} \mathcal{G}_{N-1} \frac{1}{N} \mathcal{G}_N\|_F^2. \quad (13) $$

This step can be interpreted as a low-rank denoising of $\mathcal{Y}$ using the TT-SVD with rank $K$. This operation provides the best rank-$K$ approximation of $\mathcal{Y}$ in the least-squared sense. Performing TT-SVD on $\mathcal{Y}$ requires $\mathcal{O}(K^2 \prod_n D_n)$ flops, which is $N$ times less costly than CP-ALS.
Then, the CP factors $\mathbf{A}$ can be recovered by exploiting their links with the TT-cores $\mathcal{G}$ (9) in the following optimization problem:

$$ \min_{\mathbf{A},\mathbf{M}} \sum_{n=2}^{N-1} (\|\mathcal{G}_n - [\mathbf{M}_{n-1}, \mathbf{A}_n, \mathbf{M}_n^\mathrm{T}]\|_F^2) + \|\mathbf{G}_1 - \mathbf{A}_1 \mathbf{M}_1^{-1}\|_F^2 + \|\mathbf{G}_N - \mathbf{M}_{N-1} \mathbf{A}_N^\mathrm{T}\|_F^2, \quad (14) $$

that accounts for the ambiguities of the TT-SVD through the $\mathbf{M}$ matrices. The problem (14) is a coupled optimization problem, because each nth ambiguity matrix $\mathbf{M}_n$ appears twice in (14): once in the term related to $\mathcal{G}_n$ and once in that related to $\mathcal{G}_{n+1}$. Although this problem is not convex, it is marginally convex with respect to each of the CP factors to estimate. For this kind of problem, theoretical convergence guarantees are available with the framework of block-coordinate descent (BCD), see [14,15] and references therein.

### 3.3. Proposed approach

This subsection proposes a new approach using BCD that exactly solves problem (14), in place of a sequential estimation of factors that only solves a series of local optimization sub-problems without reaching any optimum. The proposed approach fully takes into account the coupling induced by the $\mathbf{M}_n$ by carefully considering the global minimization of the cost function of (14). Each iteration of the algorithm will update every factor while keeping the others fixed. This is made possible thanks to a block-coordinate descent that alternates between all the parameters of interest. Within this framework, the proposed method is guaranteed to reach a stationary point of (14) [15]. The resulting algorithm, called Dimensionality Reduction, joint Estimation of the Ambiguity Matrices and the CP FACTors (DREAMFAC), is described below.

First, DREAMFAC estimates the $\mathbf{G}_1$, $\mathcal{G}_N$ and $\mathcal{G}_n$ $\forall n \in \{2, ..., N-1\}$ thanks to a TT-SVD applied to $\mathcal{Y}$. This step solves problem (13). Then, the coupling constraints on the ambiguity matrices lead to the following coupled model:

$$ \left\{ \begin{aligned} \mathbf{G}_1 \mathbf{M}_1 &= \mathbf{A}_1, \\ \mathbf{G}_2^{(1)\mathrm{T}} &= \mathbf{M}_1 (\mathbf{M}_2^{-\mathrm{T}} \odot \mathbf{A}_2)^{\mathrm{T}}, \\ \mathbf{G}_n^{(3)} \mathbf{M}_n &= (\mathbf{A}_n \odot \mathbf{M}_{n-1}), & \forall 2 \le n \le N-2, \\ \mathbf{G}_{n+1}^{(1)\mathrm{T}} &= \mathbf{M}_n (\mathbf{M}_{n+1}^{-\mathrm{T}} \odot \mathbf{A}_{n+1})^{\mathrm{T}}, & \forall 2 \le n \le N-2, \\ \mathbf{G}_{N-1}^{(3)} \mathbf{M}_{N-1} &= (\mathbf{A}_{N-1} \odot \mathbf{M}_{N-2}), \\ \mathbf{M}_{N-1} \mathbf{A}_N^{\mathrm{T}} &= \mathcal{G}_N. \end{aligned} \right. \quad (15) $$

These equations are obtained by unfolding the tensors $\mathcal{G}_n$ (resp. $\mathcal{G}_{n+1}$, $n \in \{1, ..., N-1\}$) in the first (resp. third dimension). Considering that the ambiguity matrices $\mathbf{M}_n$ are fixed, the $\mathbf{A}_n$ ($n \in \{2, ..., N-1\}$) can be updated from (15) by unfolding the $\mathcal{G}_n$ in the second dimension and solving a simple least-squares problem:

$$ \mathbf{A}_n = \arg \min_{\mathbf{A}_n} \|\mathbf{G}_n^{(2)\mathrm{T}} - \mathbf{A}_n (\mathbf{M}_n^{-\mathrm{T}} \odot \mathbf{M}_{n-1})^{\mathrm{T}}\|_F^2. \quad (16) $$

Considering that the CP factors $\mathbf{A}_n$ are fixed, the coherence between TT-cores leads to new coupled updates for the $\mathbf{M}_n$ matrices:

$$ \mathbf{M}_1 = \arg \min_{\mathbf{M}_1} \|\mathbf{G}_1 \mathbf{M}_1 - \mathbf{A}_1\|_F^2 + \|\mathbf{G}_2^{(1)\mathrm{T}} - \mathbf{M}_1 (\mathbf{M}_2^{-\mathrm{T}} \odot \mathbf{A}_2)^{\mathrm{T}}\|_F^2, \quad (17) \\ \forall 2 \le n \le N-2, $$

$$ \mathbf{M}_n = \arg \min_{\mathbf{M}_n} \|\mathbf{G}_n^{(3)} \mathbf{M}_n - (\mathbf{A}_n \odot \mathbf{M}_{n-1})\|_F^2 + \|\mathbf{G}_{n+1}^{(1)\mathrm{T}} - \mathbf{M}_n (\mathbf{M}_{n+1}^{-\mathrm{T}} \odot \mathbf{A}_{n+1})^{\mathrm{T}}\|_F^2, \quad (18) $$

$$ \mathbf{M}_{N-1} = \arg \min_{\mathbf{M}_{N-1}} \|\mathbf{G}_{N-1}^{(3)} \mathbf{M}_{N-1} - (\mathbf{A}_{N-1} \odot \mathbf{M}_{N-2})\|_F^2 + \|\mathbf{G}_N - \mathbf{M}_{N-1} \mathbf{A}_N^\mathrm{T}\|_F^2. \quad (19) $$

### Algorithm 1 DREAMFAC

**input:** Observation tensor $\mathcal{Y}$, CP-rank $K$ (assumed known)

**output:** CP factors $\mathbf{A}_1, ..., \mathbf{A}_N$

Estimate $\mathbf{G}_1$, $\mathcal{G}_N$ and $\mathcal{G}_n$ $\forall n \in \{2, ..., N-1\}$ using TT-SVD on $\mathcal{Y}$ repeat

Estimate the $\mathbf{A}_n$ using (16)

Estimate the $\mathbf{M}_n$ using (17)–(19) (Multi-ALS)

Estimate $\mathbf{A}_1$ and $\mathbf{A}_N$ using (9)

**until** convergence

Solutions to these least-squares problems are explicit. Appendix A gathers the resulting closed-form expressions for the ambiguity matrices. In the proposed approach, the $(\mathbf{M}_n)_{1 \le n \le N-1}$ and $(\mathbf{A})_{2 \le n \le N-1}$ matrices are all updated in the same loop, termed Multi-ALS. Algorithm 1 describes this new procedure, called Dimensionality Reduction, joint Estimation of the Ambiguity Matrices and the CP FACTors (DREAMFAC). The coherence between the TT-cores is ensured by the coupled estimation of the ambiguity matrices $\mathbf{M}_n$. Hence DREAMFAC fully exploits the information comprised in model (15), which constitutes an efficient solution to problem (14).

Turning to the algorithmic complexity, TT-SVD requires $\mathcal{O}(K^2 \Pi_n D_n)$ flops. Estimating the matrices $\mathbf{A}_n$ requires $\mathcal{O}(K^3 \sum_n D_n)$ flops. The coupled updates of the $\mathbf{M}_n$ can be handled by a fast solver for the Sylvester equations [34–36]. As a whole, these solvers require $\mathcal{O}(NK^3)$ flops so that the complexity is controlled by the (small) rank $K$. For comparison, JIRAFE uses TT-SVD followed by $\mathcal{O}(K^3 \sum_n D_n)$ flops to estimate the $\mathbf{A}_n$. Then, the estimation of the $\mathbf{M}_n$ requires $\mathcal{O}(K^3 \sum_n D_n)$ flops in place of $\mathcal{O}(NK^3)$ for DREAMFAC. Note that both methods are faster than CP-ALS, whose complexity is in $\mathcal{O}(NK^2 \prod_n D_n)$. This analysis of complexity is in favor of DREAMFAC for high-dimensional tensors.

### 3.4. Discussion with respect to state-of-the-art

In [12], the authors proposed the fast sequential optimization method JIRAFE that produced a sub-optimal solution to (14). After estimating the TT-cores by solving (13), JIRAFE performs CP-ALS on a single $\mathcal{G}_n$ ($n \in \{2, ..., N-1\}$), e.g. on $\mathcal{G}_2$, by minimizing the criterion

$$ \min_{\mathbf{M}_1, \mathbf{A}_2, \mathbf{M}_2} \|\mathcal{G}_2 - [\mathbf{M}_1, \mathbf{A}_2, \mathbf{M}_2^{-\mathrm{T}}]\|_F^2, \quad (20) $$

Since the three CP factors $\mathbf{M}_1$, $\mathbf{A}_2$, $\mathbf{M}_2$ underlying $\mathcal{G}_2$ are estimated thanks to an alternate least-squares algorithm, this step was referred to as Tri-ALS. In JIRAFE, the resulting estimate of $\mathbf{M}_2$ is then fixed and propagated to other $\mathcal{G}_n$, $n \in \{3, ..., N-1\}$. Then, for each $\mathcal{G}_n$, the cost function to minimize is

$$ \min_{\mathbf{A}_n, \mathbf{M}_n} \|\mathcal{G}_n - [\mathbf{M}_{n-1}, \mathbf{A}_n, \mathbf{M}_n^{-\mathrm{T}}]\|_F^2, \quad (21) $$

where $\mathbf{M}_{n-1}$ is supposed to be known. This operation is termed Bi-ALS. Finally, $\mathbf{A}_1$ and $\mathbf{A}_N$ are obtained using the links provided in (9).

By sequentially solving (20) then (21) instead of the global problem (14), JIRAFE only performs a suboptimal estimation of the ambiguity matrices $\mathbf{M}_n$. For instance, $\mathbf{M}_2$ is estimated only from $\mathcal{G}_2$ even though it is shared between the CPD of $\mathcal{G}_2$ and $\mathcal{G}_3$, as shown in the coupled model (15). More generally, each matrix $\mathbf{M}_n$ is shared between the CPD of $\mathcal{G}_n$ and $\mathcal{G}_{n+1}$. As such, some important information contained in (14) is ignored by JIRAFE so that the coherence between the TT-cores is lost. As a result, the performance of JIRAFE is highly dependent on the choice of the TT-core considered for the initial Tri-ALS step. By first estimating the CP factors underlying $\mathcal{G}_2$, the estimation error on the $\mathbf{M}_n$, $\mathbf{A}_n$ for $n > 2$ is expected to increase with $n$. As a result, $\mathbf{M}_{N-1}$ and $\mathbf{A}_N$ shall inherit from the estimation errors of previous CP factors. Conversely, if Tri-ALS was performed first on $\mathcal{G}_{N-1}$ before propagating estimations down to $n = 1$, $\mathbf{A}_1$ and $\mathbf{M}_1$ would not be estimated accurately. Fig. 1 shows a
<figure><img src="image_3.png" /><figcaption>Fig. 1. Graphical representation of JIRAFE (black) and DREAMFAC (blue).</figcaption></figure>

graphical representation of model (15) and summarizes the differences between JIRAFE and DREAMFAC. By resorting to a coupled strategy that actually solves (14), the main limitations of JIRAFE are corrected in the proposed algorithm DREAMFAC. This is done in practice without any additional cost, since the complexity of DREAMFAC is at most that of JIRAFE. Indeed, the TT-SVD step, and the estimation of the $A_n$ factor matrices have the same complexity for both DREAMFAC and JIRAFE. However, the number of flops necessary to estimate the $M_n$ factor matrices is reduced in DREAMFAC owing to fast solvers for Sylvester equations, see Section 3.3. Therefore, the complexity of DREAMFAC can be made even smaller in applications, e.g., harmonic retrieval in MIMO channels, where the number of elements $K$ is often very small. Furthermore, since DREAMFAC has convergence guarantees to a stationary point owing to the BCD framework, it is expected to recover the CP factors while JIRAFE is not.

## 4. Cramér-Rao bounds for the proposed approach

The goal of this Section is to highlight the capacity of the proposed approach to provide the best possible estimation of the CP factors by solving (14), after the TT-SVD step. Indeed, the proposed approach relies on problem (14) to recover the CP factors without having to solve (12), which is an ill-posed problem. In this Section, it is assumed that the observed tensor $\mathcal{Y}$ admits an exact CPD, therefore TT-SVD yields no approximation error. Should $\mathcal{Y}$ not admit such an exact CPD, a small residual approximation error arising from the TT-SVD approximation would appear in the results.

### 4.1. General probabilistic framework

To derive appropriate performance bounds, it is necessary to embed the problem in an appropriate probabilistic framework. This requires to properly define the probabilistic model and the parameters of interest. A general coupled model has the form

$$ \begin{cases} g_n \sim f_\omega(g_n) & \forall n \in \{1, \dots, N\}, \\ h_n(\omega) = 0, & \end{cases} \qquad (22) $$

where $g_n$ denote the observations. The functions $f_\omega(g_n)$ are the probability density functions (PDF) of the random real datasets $g_n$, parameterized by the unknown deterministic real parameter vector $\omega$. Several assumptions will be necessary: (i) the PDFs $f_\omega(g_n)$ are non-redundant functions differentiable w.r.t. $\omega$, and their support as functions of $g_n$ do not depend on $\omega$; (ii) the $g_n$ are statistically independent.<fnref n="2" /> The

constraints in the model are described by the functions $h_n$ that are non-redundant deterministic vector functions, everywhere differentiable with respect to $\omega$.

The model (15) can be rewritten equivalently under the form (22). First, the link between the TT-cores and the $A_n$, $M_n$ matrices is such that

$$ G_1 = A_1 M_{1,1}^{-1}, \quad G_n = [M_{n-1,n}, A_n, M_{n,n}^{-\mathrm{T}}], \quad G_N = M_{N-1,N} A_N^{\mathrm{T}}. \qquad (23) $$

This paper proposes an approach to solve the approximation problem (14) in the least-squares sense based on model (12). Therefore, it is assumed that the observations are Gaussian. To be more precise, let us note $g_n = \mathrm{vec}(G_n)$ for $n \in \{1, N\}$ and $g_n = \mathrm{vec}(G_n)$ for $n \in \{2, \dots, N-1\}$. The vectors $g_n$ are random real Gaussian distributed datasets parameterized by their mean, i.e.,

$$ g_n \sim \mathcal{N}(\mu_n(\omega), \sigma_n^2 I) \quad \text{where} \quad \mu_n = \mathrm{vec}\{[M_{n-1,n}, A_n, M_{n,n}^{-\mathrm{T}}]\}, \qquad (24) $$

and $\sigma_n^2$ denote the variances of the Gaussian noise on $g_n$.

The unknown real deterministic parameter $\omega$ is such that

$$ \omega = [\mathrm{vec}\{A_1\}^\mathrm{T} \ \mathrm{vec}\{M_{1,1}\}^\mathrm{T} \ \mathrm{vec}\{M_{1,2}\}^\mathrm{T} \ \mathrm{vec}\{A_2\}^\mathrm{T} \ \dots \ \mathrm{vec}\{A_N\}^\mathrm{T}]. \qquad (25) $$

Second, the coherence between the TT-cores is ensured by the set of constraints

$$ M_{n,n} = M_{n,n+1} \quad \forall n \in \{1, \dots, N-1\}. \qquad (26) $$

Reshaping (26) under the form of vector functions depending on $\omega$ yields

$$ h_n(\omega) = \mathrm{vec}\{M_{n,n}\} - \mathrm{vec}\{M_{n,n+1}\} \quad \forall n \in \{1, \dots, N-1\}. \qquad (27) $$

As a result, the model (24) under constraints (27) characterizes (15) as a general coupled model.

### 4.2. Calculation of the constrained Cramér-Rao bound

Evaluating the performance of the coupled model (24) under constraints (27) is necessary to illustrate the efficiency of DREAMFAC for the estimation of $\omega$. The standard tool for this task is the constrained Cramér-Rao bound (CCRB) [24,26], defined<fnref n="3" /> as

$$ \mathrm{CCRB}(\omega) = U [U^\mathrm{T} F U]^{-1} U^\mathrm{T}, \qquad (28) $$

where $F \triangleq F(\omega)$ is the Fisher Information Matrix (FIM) of the model, $J_h$ is the Jacobian matrix related to the constraints $h(\omega) = 0$ and $U$ is a basis of ker($J_h$).

<footnote n="2">² This assumption may not hold when using TT-SVD to estimate the $g_n$. Nevertheless, this assumption is necessary to derive the performance bounds.</footnote>

<footnote n="3">³ If $F$ is invertible, then (28) and the alternative expression for the CCRB provided in [37] are equivalent [24, Corollary 1].</footnote>
Additionally to the model parameter $\omega$, let us define $\mathbf{x} = \text{vec}\{\mathcal{X}\} \in \mathbb{R}^c$ ($c = \Pi_n D_n$), that represents the vectorized low-rank approximation of $\mathcal{Y}$. The parameter $\mathbf{x}$ can be linked to the model parameters through the relationship $\mathbf{x} - \mathbf{h}_{\mathcal{X}}(\omega) = 0$, where $\mathbf{h}_{\mathcal{X}}$ is a non-redundant deterministic vector function, everywhere differentiable with respect to $\omega$ such that

$$ \mathbf{h}_{\mathcal{X}}(\omega) = \text{vec}(\llbracket \mathbf{A}_1, \dots, \mathbf{A}_N \rrbracket). \quad (29) $$

The above parameterization allows us to compute the CCRB for $\mathbf{x}$ [38, p. 125] as

$$ \text{CCRB}(\mathbf{x}) = \left[ \frac{\partial \mathbf{h}_{\mathcal{X}}(\omega)}{\partial \omega^{\top}} \right]^{\top} \text{CCRB}(\omega) \left[ \frac{\partial \mathbf{h}_{\mathcal{X}}(\omega)}{\partial \omega^{\top}} \right]. \quad (30) $$

### 4.2.1. Expression of the FIM

The FIM for $\omega$ first needs to be computed in order to obtain the CCRB. The FIM for $\omega$ is a block-diagonal matrix of the form

$$ \mathbf{F}(\omega) = \begin{bmatrix} \mathbf{F}_1(\omega) & \mathbf{0} & \cdots & \mathbf{0} \\ \mathbf{0} & \ddots & \cdots & \mathbf{0} \\ \vdots & \ddots & \ddots & \vdots \\ \mathbf{0} & \mathbf{0} & \cdots & \mathbf{F}_N(\omega) \end{bmatrix}, \quad (31) $$

where the blocks $\mathbf{F}_n(\omega)$ encode the contributions of each $g_n$ to the estimation of $\omega$. They are obtained by using the Slepian-Bangs formula [39]:

$$ \mathbf{F}_n(\omega) = \frac{1}{\sigma_n^2} \left[ \frac{\partial \mu_n(\omega)}{\partial \omega^{\top}} \right]^{\top} \left[ \frac{\partial \mu_n(\omega)}{\partial \omega^{\top}} \right]. \quad (32) $$

The expression of $\left[ \frac{\partial \mu_n(\omega)}{\partial \omega^{\top}} \right]$ is obtained simply for $n=1$ and $n=N$:

$$ \frac{\partial \mu_1(\omega)}{\partial \text{vec}\{\mathbf{A}_1\}^{\top}} = \mathbf{M}_{1,1}^{\top} \boxtimes \mathbf{I}_{D_1}, \quad (33) $$

$$ \frac{\partial \mu_1(\omega)}{\partial \text{vec}\{\mathbf{M}_{1,1}\}^{\top}} = \mathbf{I}_K \boxtimes \mathbf{A}_1, \quad (34) $$

$$ \frac{\partial \mu_1(\omega)}{\partial \text{vec}\{\mathbf{M}_{N-1,N}\}^{\top}} = \mathbf{A}_N \boxtimes \mathbf{I}_K, \quad (35) $$

$$ \frac{\partial \mu_N(\omega)}{\partial \text{vec}\{\mathbf{A}_N\}^{\top}} = \Pi_N \left( \mathbf{M}_{N-1,N}^{-1} \boxtimes \mathbf{I}_{D_N} \right), \quad (36) $$

where $\Pi_N$ is a permutation matrix linking the entries of $\text{vec}\{\mathbf{G}_N^\top\}$ to those of $\text{vec}\{\mathbf{G}_N\}$. For $n \in \{2, \dots, N-1\}$, $\left[ \frac{\partial \mu_n(\omega)}{\partial \omega^{\top}} \right]$ can be computed using relationships between tensor unfoldings:

$$ g_n = \left[ (\mathbf{M}_{n,n}^{-\top} \odot \mathbf{A}_n) \boxtimes \mathbf{I}_K \right] \text{vec}\{\mathbf{M}_{n-1,n}\} \quad (37) $$

$$ = \Pi_n^{(2,1)} \left[ (\mathbf{M}_{n,n}^{-\top} \odot \mathbf{M}_{n-1,n}) \boxtimes \mathbf{I}_{D_n} \right] \text{vec}\{\mathbf{A}_n\} \quad (38) $$

$$ = \Pi_n^{(3,1)^{\top}} \left[ \mathbf{I}_K \boxtimes (\mathbf{A}_n \odot \mathbf{M}_{n-1,n}) \right] \text{vec}\{\mathbf{M}_{n,n}\}, \quad (39) $$

where $\Pi_n^{(2,1)}$ and $\Pi_n^{(3,1)}$ are permutation matrices that link the entries of $\text{vec}\{\mathbf{G}_n^{(2)}\}$ (resp. $\text{vec}\{\mathbf{G}_n^{(3)}\}$) to those of $g_n = \text{vec}\{\mathbf{G}_n^{(1)}\}$. To ease the notation, we define the matrices $\mathbf{S}_{1,n}$, $\mathbf{S}_{2,n}$, $\mathbf{S}_{3,n}$ ($n \in \{2, \dots, N-1\}$):

$$ \mathbf{S}_{1,n} = \left[ (\mathbf{M}_{n,n}^{-\top} \odot \mathbf{A}_n) \boxtimes \mathbf{I}_K \right], \quad (40) $$

$$ \mathbf{S}_{2,n} = \Pi_n^{(2,1)} \left[ (\mathbf{M}_{n,n}^{-\top} \odot \mathbf{M}_{n-1,n}) \boxtimes \mathbf{I}_{D_n} \right], \quad (41) $$

$$ \mathbf{S}_{3,n} = \Pi_n^{(3,1)^{\top}} \left[ \mathbf{I}_K \boxtimes (\mathbf{A}_n \odot \mathbf{M}_{n-1,n}) \right]. \quad (42) $$

As a result, we have

$$ \left[ \frac{\partial \mu_n(\omega)}{\partial \omega^{\top}} \right] = \begin{cases} \begin{bmatrix} \mathbf{M}_{1,1}^{\top} \boxtimes \mathbf{I}_{D_1} & \mathbf{I}_K \boxtimes \mathbf{A}_1 \\ \mathbf{S}_{1,n} & \mathbf{S}_{2,n} \\ \mathbf{A}_N \boxtimes \mathbf{I}_K & \Pi_N \left( \mathbf{M}_{N-1,N}^{-1} \boxtimes \mathbf{I}_{D_N} \right) \end{bmatrix} & \text{for } n=1, \\ \begin{bmatrix} \mathbf{M}_{n,n}^{-\top} \odot \mathbf{A}_n & \mathbf{S}_{n,n} \\ \mathbf{S}_{n,n} & \mathbf{S}_{n,n} \end{bmatrix} & \forall n \in \{2, \dots, N-1\}, \\ \begin{bmatrix} \mathbf{A}_n \boxtimes \mathbf{I}_K & \Pi_N \left( \mathbf{M}_{N-1,N}^{-1} \boxtimes \mathbf{I}_{D_N} \right) \end{bmatrix} & \text{for } n=N. \end{cases} \quad (43) $$

The closed-form expressions for the $\mathbf{F}_n(\omega)$ are provided in Appendix B.

### 4.2.2. Expression of the CCRB

The Jacobian matrix $\mathbf{J}_h$ is obtained by deriving the functions $\mathbf{h}_n(\omega)$ with respect to all the elements of $\omega$ into a block-matrix. Given the constraints on $\omega$ in (27), it holds that

$$ \frac{\partial \mathbf{h}_n(\omega)}{\partial \text{vec}\{\mathbf{A}_n\}^{\top}} = \mathbf{0}, \quad \forall n \in \{1, \dots, N-1\}, \quad (44) $$

$$ \frac{\partial \mathbf{h}_n(\omega)}{\partial \text{vec}\{\mathbf{M}_{m-1,m}\}^{\top}} = \begin{cases} \mathbf{I}_{K^2} & \text{if } m=n \\ \mathbf{0} & \text{if } m \neq n, \end{cases} \quad (45) $$

$$ \frac{\partial \mathbf{h}_n(\omega)}{\partial \text{vec}\{\mathbf{M}_{m,m}\}^{\top}} = \begin{cases} -\mathbf{I}_{K^2} & \text{if } m=n \\ \mathbf{0} & \text{if } m \neq n. \end{cases} \quad (46) $$

As a result, $\mathbf{J}_h$ is such that

$$ \mathbf{J}_h = \begin{bmatrix} \frac{\partial \mathbf{h}_1(\omega)}{\partial \omega^{\top}} \\ \vdots \\ \frac{\partial \mathbf{h}_{N-1}(\omega)}{\partial \omega^{\top}} \end{bmatrix} = \begin{bmatrix} \mathbf{0} & \mathbf{I}_{K^2} & -\mathbf{I}_{K^2} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} \\ \vdots & \mathbf{0} & \mathbf{0} & \vdots & \mathbf{I}_{K^2} & -\mathbf{I}_{K^2} & \mathbf{0} & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{I}_{K^2} & -\mathbf{I}_{K^2} \end{bmatrix}. \quad (47) $$

The matrix $\mathbf{U}$ can then be obtained by simply solving $\mathbf{J}_h\mathbf{U} = \mathbf{0}$. Therefore $\mathbf{U}$ is the identity matrix of size $K(2N - 2 + \sum_n D_n)$. The CCRB submatrices for the $\mathbf{A}_n$ and the $\mathbf{M}_{n-1,n}$ are obtained by developing expression (28):

$$ \text{CCRB}(\mathbf{A}_n) = \mathbf{F}(\mathbf{A}_n)^{-1}, \quad \text{CCRB}(\mathbf{M}_{n-1,n}) = \mathbf{F}(\mathbf{M}_{n-1,n})^{-1}. \quad (48) $$

## 4.3. Reparameterized CRB for reconstruction of the low-rank tensor

The expression for $\mathbf{h}_{\mathcal{X}}(\omega)$ is necessary to compute the reparametrized CRB in (30). The relationships between tensor unfoldings give

$$ \begin{aligned} \mathbf{x} &= [(\mathbf{A}_N \odot \dots \odot \mathbf{A}_2) \boxtimes \mathbf{I}_{D_1}] \text{vec}\{\mathbf{A}_1\} \\ &= \dots = \Pi_{\mathcal{X}}^{(n,1)} [(\mathbf{A}_{N-1} \odot \dots \odot \mathbf{A}_1) \boxtimes \mathbf{I}_{D_N}] \text{vec}\{\mathbf{A}_N\}, \end{aligned} \quad (49) $$

where the $\Pi_{\mathcal{X}}^{(n,1)}$ are permutation matrices that link the entries of the $\text{vec}\{\mathbf{X}^{(n)}\}$ vectors to those of $\text{vec}\{\mathbf{X}^{(1)}\} = \mathbf{x}$. The matrices $\mathbf{S}_{n,\mathcal{X}}$ are defined such that

$$ \mathbf{S}_{n,\mathcal{X}} = \Pi_{\mathcal{X}}^{(n,1)} \left[ (\mathbf{A}_N \odot \dots \odot \mathbf{A}_{n+1} \odot \mathbf{A}_{n-1} \odot \dots \odot \mathbf{A}_1) \boxtimes \mathbf{I}_{D_n} \right], \quad (50) $$

finally yielding

$$ \left[ \frac{\mathbf{h}_{\mathcal{X}}(\omega)}{\partial \omega^{\top}} \right] = \left[ \mathbf{S}_{1,\mathcal{X}} \quad \mathbf{0} \quad \mathbf{0} \quad \mathbf{S}_{2,\mathcal{X}} \quad \mathbf{0} \quad \mathbf{0} \quad \cdots \quad \mathbf{S}_{N,\mathcal{X}} \right]. \quad (51) $$

Therefore,

$$ \text{CCRB}(\mathbf{x}) = \text{Diag}\{\mathbf{S}_{1,\mathcal{X}}^{\top} \mathbf{F}(\mathbf{A}_1)^{-1} \mathbf{S}_{1,\mathcal{X}}, \dots, \mathbf{S}_{N,\mathcal{X}}^{\top} \mathbf{F}(\mathbf{A}_N)^{-1} \mathbf{S}_{N,\mathcal{X}}\}. \quad (52) $$

Given the expressions of the constrained Cramér-Rao bounds in Eq. (48) and Eq. (52), the performance of the proposed approach can now be evaluated numerically.

## 5. Simulations

This section illustrates the optimality of DREAMFAC in practical cases. Section 5.1 highlights the performance of DREAMFAC in synthetic test cases. Section 5.2 shows that, in all scenarios, DREAMFAC reaches the Cramér-Rao bounds derived in Section 4. Section 5.3 is an application to channel estimation in dual-polarized MIMO systems.
<figure><img src="image_4.png" /><figcaption>Fig. 2. Columns of the reference (black dashed lines) and estimated factors as returned by JIRAFE with forward propagation (red circles), JIRAFE with backward propagation (blue squares), and DREAMFAC (green diamonds).</figcaption></figure>

## 5.1. Recovery of the CP factors

A 7-order ($N = 7$) tensor $\mathcal{X}$ with $D_1 = \dots = D_N = 6$ and $K = 3$ is taken as a reference. This tensor admits a reference CPD as in (10) with i.i.d. entries generated from the normal distribution. The noisy tensor $\mathcal{Y}$ was generated using isotropic white Gaussian noise to yield a SNR of 20 dB.

The proposed approach is compared to two JIRAFE procedures. In the first one, Tri-ALS is initially performed on $\mathcal{G}_2$ and the factors is propagated towards the higher $n$. In the second one, Tri-ALS is performed on $\mathcal{G}_{N-1}$ with propagation towards the lower $n$. The TT-cores are estimated using TT-SVD with $K = 3$. The CP factors were initialized as $\mathbf{A}_n = \mathbf{U}_n (\mathbf{\Sigma}_n)^{\frac{1}{2}}$ using the SVD of the $n$th mode unfolding of $\mathcal{Y}$ with rank $K$, namely $\mathbf{Y}^{(n)} = \mathbf{U}_n \mathbf{\Sigma}_n \mathbf{V}_n$. The algorithms have a maximum of 1000 iterations, and the results are averaged over 100 noise realizations.

Fig. 2 shows columns of the reference and estimated factors. The columns of the $\mathbf{A}_n$ factors for $n \in \{2, \dots, N-1\}$ are correctly estimated by all the algorithms, except for a few outliers in $\mathbf{A}_6$. The JIRAFE-like procedures provided an incorrect estimation of the factors at the end of the tensor trains, respectively $\mathbf{A}_N$ and $\mathbf{A}_1$. Differently, Algorithm 1 correctly estimates $\mathbf{A}_1$ and $\mathbf{A}_N$ while usually providing a better estimate than the JIRAFE procedures.

To further illustrate the propagation of the error through the TT-cores by JIRAFE, the normalized mean square error between the estimated $\hat{\mathbf{A}}_n$ and the reference $\mathbf{A}_n$ is considered:

$$ \text{NMSE} = \frac{\|\hat{\mathbf{A}}_n - \mathbf{A}_n\|_F^2}{\|\mathbf{A}_n\|_F^2}. \qquad (53) $$

Fig. 3 shows in semi-log scale the NMSE for each factor, as provided by the three algorithms. As a comparison, the value of the uniform CCRB, normalized by $\|\mathbf{A}_n\|_F^2$, is also displayed. For the two JIRAFE procedures, the highest NMSE corresponds to the last estimated factor, resp. $\mathbf{A}_N$ for the forward procedure and $\mathbf{A}_1$ for the backward procedure. The NMSE provided by DREAMFAC is smaller than that provided by the JIRAFE algorithms. Since DREAMFAC converges to a stationary point of problem (14) (contrary to JIRAFE), it is the MLE associated with problem (14) under the assumption of i.i.d. white Gaussian noise on the observations. As such, it is expected to reach the CCRB, which is verified in the experiments. DREAMFAC indeed reaches the NMSE depicted by the normalized CRB.

## 5.2. Optimality of DREAMFAC using the CCRB

This subsection compares the performance of DREAMFAC to the constrained Cramér-Rao bound obtained in Section 4.

The reference tensor $\mathcal{X}$ is a 7th-order tensor with $D_1 = \dots = D_N = 6$ and $K = 3$. The entries of the true CP factors $\mathbf{A}_n$ were generated once as

<figure><img src="image_5.png" /><figcaption>Fig. 3. NMSE provided by JIRAFE with forward propagation (red dots), backward propagation (blue squares), and DREAMFAC (green diamonds), and normalized CRB.</figcaption></figure>

i.i.d. real standard Gaussian variables. White Gaussian noise is added to $\mathcal{X}$ to form the observation. We assume that the noise level was the same for all TT-cores: this assumption is reasonable since the $\mathcal{G}_n$ are all estimated from $\mathcal{Y}$ using TT-SVD. on the TT-cores. The SNR on the observed tensors in dB is defined as $\text{SNR}_i = 10 \log_{10} (\|\mathcal{G}_i\|_F^2 / \|\mathcal{E}_i\|_F^2)$, ($i = 1, \dots, N$).

The model parameters are retrieved using JIRAFE and DREAMFAC. The CP factors and ambiguity matrices are initialized randomly. The permutation ambiguities in the estimated factors are corrected using the Hungarian algorithm [40] to make comparisons aligned with the reference. The experiments show the uniform MSE and uniform CCRB obtained from the MSE and CCRB matrix traces, as widely considered in, e.g., [41–43]. The expressions for the bounds proposed in this paper permit the computation of the uniform CCRB by taking the trace of these matrices.

Fig. 4 shows in semi-log scale the uniform bounds and MSEs for all the entries of the $\mathbf{A}_n$ and $\mathbf{M}_n$, $\omega$ and $\mathcal{X}$ as a function of the SNR.

It is noticeable that the uniform MSE produced by DREAMFAC reaches the uniform CCRB. As expected, the proposed approach is optimal for estimation of the parameters and reconstruction of $\mathcal{X}$. JIRAFE yields a higher MSE. These observations were to be expected, JIRAFE is a suboptimal optimization method, as highlighted in Section 3.4. The two algorithms depict the same kind of behavior with respect to the SNR. This is particularly visible for estimation of the ambiguity matrices. DREAMFAC permits to gain 5 dB for the estimation of the $\mathbf{A}_n$, and 3 dB for the reconstruction of $\mathcal{X}$ compared to JIRAFE.
<figure><img src="image_6.png" /><figcaption>Fig. 4. Uniform CCRB and MSE provided by JIRAFE and DREAMFAC, for (a) estimation of all the $A_n$, (b) estimation of all the $M_n$, (c) estimation of $\omega$ and (d) reconstruction of $X$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5. Uniform CCRB and uniform MSE provided by JIRAFE and DREAMFAC for the estimation of the $A_n$.</figcaption></figure>

Fig. 5 shows in semi-log scale the uniform CCRB for $A_1, A_3, A_5$ and $A_7$, as well as the uniform MSEs, as a function of the SNR. The results for $A_2, A_4$ and $A_6$ are similar. The MSE obtained with JIRAFE was higher, and progressively increased for high values of $n$. The uniform MSE obtained with DREAMFAC also reaches the CCRB for each factor, and its scale barely varies trough the $A_n$. Therefore DREAMFAC is also optimal for the estimation of each CP factor individually.

### 5.3. Application to channel estimation in dual-polarized MIMO systems

The problem of channel estimation in dual-polarized massive MIMO aims at recovering the channel parameters at the base station (angles of

arrival and departure, path gains, polarization parameters). A tensor-based approach was proposed in [6]. The MIMO channel was recast as a fourth-order tensor admitting a CPD. In [29], the JIRAFE procedure was adapted to treat the channel estimation problem with receiver and transmitter rectangular arrays, hence the channel was viewed as a fifth-order tensor.

#### 5.3.1. Model description

The steering vectors for the $k$th path between a Uniform Rectangular Array (URA) transmitter of size $M_T^x \times M_T^y$ and a URA receiver of size $M_R^x \times M_R^y$ are such that

$$a_T(k) = a_T^x(k) \otimes a_T^y(k), \quad a_R(k) = a_R^x(k) \otimes a_R^y(k), \qquad (54)$$

with e.g., $a_T^x(k) = [1, \exp(j\omega_T^x(k)), \dots, \exp(j\omega_T^x(k)(M_T^x - 1))]^T$, and likewise for $a_T^y(k), a_R^x(k), a_R^y(k)$. As a result, for $K$ paths, the steering matrices in transmission and in reception are

$$A_T = A_T^x \odot A_T^y, \quad A_R = A_R^x \odot A_R^y, \qquad (55)$$

with e.g., $A_T^x = [a_T^x(1), \dots, a_T^x(K)]$.

The path-loss matrix $B \in \mathbb{C}^{4 \times K}$ contains the path-loss parameters. Let the $(p, q)$-th subchannel correspond to $p \in \{V_R, H_R\}$ for the vertical (V) polarized and horizontal (H) polarized receive antennas, and $q \in \{V_T, H_T\}$ for the V-polarized and H-polarized transmit antennas. For the $k$th path and for the $(p, q)$-th subchannel, $\beta_k^{(p,q)}$ is the path-loss parameter. Therefore, noting $\beta_k^{(p,q)} = [\beta_1^{(p,q)}, \dots, \beta_K^{(p,q)}]$,

$$B = [\beta^{(V_R, V_T)} \quad \beta^{(V_R, H_T)} \quad \beta^{(H_R, V_T)} \quad \beta^{(H_R, H_T)}], \qquad (56)$$

As a result, the channel tensor $H \in \mathbb{C}^{M_T^x \times M_T^y \times M_R^x \times M_R^y \times 4}$ can be written as

$$H = \llbracket A_R^x, A_R^y, (A_T^x)^*, (A_T^y)^*, B \rrbracket + \mathcal{E}, \qquad (57)$$

which corresponds to a fifth-order CPD of rank $K$. The noise term $\mathcal{E}$ encompasses the background noise and the estimation error due to the pre-estimation of the unstructured channel, and can be modeled as zero-mean circularly complex Gaussian random variables.
<figure><img src="image_8.png" /><figcaption>Fig. 6. Columns of the reference factors (black dashed lines) and estimated factors as returned by JIRAFE (red circles) and DREAMFAC (green diamonds).</figcaption></figure>

The main assumption of (57) is that $K \le \min(M_T^x, M_T^y, M_R^x, M_R^y)$, i.e., the steering matrices are full column-rank. Furthermore, it is assumed that there are few dominant paths, i.e., $K < 4$ and $B$ is a full column-rank matrix. The assumption of a small $K$ is usually made in massive MIMO scenarios.

### 5.3.2. Results

The matrices $A_T^x, A_T^y, A_R^x, A_R^y$ were generated once based on single random realizations of the angular frequencies $\omega_T^x(k), \omega_T^y(k), \omega_R^x(k), \omega_R^y(k)$ following a uniform distribution on $[0, \pi]$. The factor $B$ was drawn from a complex Gaussian distribution with zero mean and unit variance. The dimensions $M_T^x = M_R^x = 10$ and $M_T^y = M_R^y = 8$, and $K = 3$ were considered.

Fig. 6 shows in semi-log scale true factors and the estimated factors provided by DREAMFAC and JIRAFE averaged over 100 noise realizations with 30 dB SNR.

The factors estimated by JIRAFE seem coherent with the reference factors, but with the wrong scale and angle. DREAMFAC recovers the CP factors correctly, with their original scale and orientation, without the need for a post-processing step.

To further assess the performance of the DREAMFAC, the MSE between the true and estimated factors, and between the angular frequencies, were considered:

$$ \text{MSE}_{R}^{x} = \sum_{k=1}^{K} (A_{R}^{x}(\cdot, k) - \hat{A}_{R}^{x}(\cdot, k))^{2}, \quad \text{MSE}_{R}^{y} = \sum_{k=1}^{K} (A_{R}^{y}(\cdot, k) - \hat{A}_{R}^{y}(\cdot, k))^{2}, \qquad (58) $$

$$ \text{MSE}_{T}^{x} = \sum_{k=1}^{K} (A_{T}^{x}(\cdot, k) - \hat{A}_{T}^{x}(\cdot, k))^{2}, \quad \text{MSE}_{T}^{y} = \sum_{k=1}^{K} (A_{T}^{y}(\cdot, k) - \hat{A}_{T}^{y}(\cdot, k))^{2}, \qquad (59) $$

$$ \text{MSE}_{\omega} = \sum_{k=1}^{K} (\omega_{T}^{x}(k) - \hat{\omega}_{T}^{x}(k))^{2} + (\omega_{T}^{y}(k) - \hat{\omega}_{T}^{y}(k))^{2} + (\omega_{R}^{x}(k) - \hat{\omega}_{R}^{x}(k))^{2} + (\omega_{R}^{y}(k) - \hat{\omega}_{R}^{y}(k))^{2}. \qquad (60) $$

The MSE were evaluated over 10 values of the SNR in $[-20, 30]$ dB, and calculated by averaging the results over 500 independent noise realizations. Computing the complex CCRB is a delicate task that requires the calculation of Wirtinger derivatives [44]. For this reason, calculation of the CRB associated with model (57) is relegated to future works.

Four tensor-based methods were benchmarked. The first one was CP-ALS [30] followed by closed-form solutions to estimate the parameters from the factors. The second one was the so-called CP-VDM, for CPD with Vandermonde factor matrix, proposed in [45]. The third one was based on the generalized eigenvalue decomposition (GEVD) [46]. The fourth one was JIRAFE followed by a Vandermonde rectification

strategy [47] (termed JIRAFE-VDM) to enforce the structure of the steering matrices. This step is crucial to estimate the factors with the correct scale and angle, as exemplified in Fig. 6. The structure rectification in JIRAFE-VDM remains a patch to improve the approximate optimization performed by JIRAFE.

Fig. 7 shows in semi-log scale the averaged MSE for recovery of the steering matrices. DREAMFAC has better performance than all other approaches, especially for low SNR. For these SNR, its performance is better than that of JIRAFE-VDM, with a gain of approximately 3 dB on the quality of the estimation. JIRAFE-VDM could achieve similar performance to DREAMFAC for higher SNR, only thanks to the Vandermonde rectification. Without this correction, the MSE yielded by JIRAFE would have been significantly higher, due to the incorrect scale in the estimated factors. DREAMFAC yields the lowest MSE over the full range of the SNR values.

Fig. 8 shows in semi-log scale the averaged MSE$_\omega$ provided by the algorithms as a function of the SNR. As for the steering matrices, DREAMFAC performs better than the benchmarked algorithms for the estimation of the angular frequencies, especially for low SNR. For instance, DREAMFAC performs as well at an SNR of -15 dB as JIRAFE-VDM at -10 dB, with lower complexity and without the need for any post-processing. These results confirm the optimality of DREAMFAC in practice when other benchmarked methods suffer from a loss of performance. This general good performance of DREAMFAC stems from its theoretical guarantees to reach a stationary point of problem (14) as well as the CCRB.

## 6. Conclusion

This paper proposes a new approach called DREAMFAC for Dimensionality Reduction, joint Estimation of the Ambiguity Matrices and the CP FACTors. Assuming that the best rank-K CPD exists, dimensionality reduction is performed by exploiting the equivalence between the low-cost TTD and the expensive CPD. It relies on a global coupled optimization scheme, instead of a local and sequential strategy adopted in previous approaches. DREAMFAC yields the optimal solution to the problem (14) of estimating the CP factors from a TT decomposition. Indeed, the use of a block-coordinate descent algorithm guarantees an exact solution to this marginally convex problem. From an estimation theory perspective, constrained Cramér-Rao bounds were derived to evaluate the optimum performance under mild assumptions. Numerical experiments show that DREAMFAC reaches these Cramér-Rao bounds. The proposed approach, DREAMFAC, is an optimal solution to the problem (14), both from an optimization perspective and as a statistical estimator. The optimality of DREAMFAC is exemplified on the realistic problem of complex harmonics retrieval in dual-polarized MIMO
<figure><img src="image_9.png" /><figcaption>Fig. 7. (a) MSE<sub>R</sub><sup>x</sup>, (b) MSE<sub>R</sub><sup>y</sup>, (c) MSE<sub>T</sub><sup>x</sup>, (d) MSE<sub>T</sub><sup>y</sup> as a function of the SNR.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8. MSE on the estimation of the angular frequencies.</figcaption></figure>

channels. DREAMFAC provides better estimation of the base station parameters than state-of-the-art methods, in particular in the difficult case of low SNR.

As a conclusion, the proposed DREAMFAC yields the optimal solution to the identification of a rank-K CPD from a TT-SVD-based TTD. The only assumption needed is that the best rank-K CPD of the tensor exists. In the quest for a better estimation of the CP factors, new alternatives to the use of the TTD, or TT-SVD within TTD, may be necessary.

CRediT authorship contribution statement

**Clémence Prévost:** Writing – review & editing, Writing – original draft, Validation, Software, Methodology, Investigation, Conceptualization. **Pierre Chainais:** Writing – review & editing, Writing – original draft, Validation, Supervision, Software, Methodology.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Closed-form expression for the coupled M<sub>n</sub> in DREAMFAC

Solutions the least-squared problems (17)–(19) can be obtained by solving Sylvester equations:

$$
M_1 (M_2^{-1} M_2^{-\mathrm{T}} \Box A_2^{\mathrm{T}} A_2) + (G_1^{\mathrm{T}} G_1) M_1 = G_1^{\mathrm{T}} A_1 + G_2^{(1)^{\mathrm{T}}} (A_2 \odot M_2), \quad (A.1)
$$

$$
\begin{align}
M_n (M_{n+1}^{-1} M_{n+1}^{-\mathrm{T}} \Box A_{n+1}^{\mathrm{T}} A_{n+1}) + (G_n^{(3)^{\mathrm{T}}} G_n^{(3)}) M_n \nonumber \\
&= G_{n+1}^{(1)^{\mathrm{T}}} (M_{n+1}^{-\mathrm{T}} \odot A_{n+1}) + G_n^{(3)^{\mathrm{T}}} (A_n \odot M_{n-1}), \tag{A.2}
\end{align}
$$

$$
\begin{align}
& M_{N-1} (A_N^\top A_N) + (G_{N-1}^{(3)} {}^\top G_{N-1}) M_{N-1} \nonumber \\
& \quad = G_N A_N + G_{N-1}^{(3)\top} (A_{N-1} \odot M_{N-2}). \tag{A.3}
\end{align}
$$

Appendix B. Closed-form expressions for the uncoupled FIM

The FIM matrices $F_n(\omega)$ introduced in (32) are block-matrices such as

$$
F_1(\omega) = \frac{1}{\sigma_1^2} \begin{bmatrix} M_{1,1}^{-1} M_{1,1}^{-\mathrm{T}} \boxtimes I_{d_1} & M_{1,1}^{-1} \boxtimes A_1 \\ M_{1,1}^{-\mathrm{T}} \boxtimes A_1^{\mathrm{T}} & I_R \boxtimes A_1^{\mathrm{T}} A_1 \end{bmatrix}, \quad (B.1)
$$

$$
F_n(\omega) = \frac{1}{\sigma_n^2} \begin{bmatrix} S_{1,n}^\top S_{1,n} & S_{1,n}^\top S_{2,n} & S_{1,n}^\top S_{3,n} \\ S_{2,n}^\top S_{1,n} & S_{2,n}^\top S_{2,n} & S_{2,n}^\top S_{3,n} \\ S_{3,n}^\top S_{1,n} & S_{3,n}^\top S_{2,n} & S_{3,n}^\top S_{3,n} \end{bmatrix} \quad \text{for } n \in \{2, \dots, N-1\}, \quad (B.2)
$$

$$
F_N(\omega) = \frac{1}{\sigma_N^2} \times \begin{bmatrix} A_N^\top A_N \boxtimes I_R & (A_N^\top \boxtimes I_R) \Pi_N (M_{N-1,N}^{-1} \boxtimes I_{D_N}) \\ (M_{N-1,N}^{-1} \boxtimes I_{D_N}) \Pi_N (A_N \boxtimes I_R) & M_{N-1,N}^{-1} M_{N-1,N}^{-1} \boxtimes I_{D_N} \end{bmatrix} \quad (B.3)
$$
Let $F(A_n)$ and $F(M_{n-1,n})$ denote the blocks containing the contributions of the $A_n$ and the $M_{n-1,n}$ to the estimation. The matrices $F(A_n)$ and $F(M_{n-1,n})$ appear in (B.1)–(B.3) such that

$$F(A_n) = \begin{cases} \mathbf{I}_R \boxtimes \mathbf{A}_1^\top \mathbf{A}_1 & \text{for } n=1, \\ \mathbf{S}_{2,n}^\top \mathbf{S}_{2,n} & \text{for } n \in \{2, \dots, N-1\}, \\ \mathbf{A}_N^\top \mathbf{A}_N \boxtimes \mathbf{I}_R & \text{for } n=N, \end{cases} \quad (B.4)$$

$$F(M_{n-1,n}) = \begin{cases} \mathbf{S}_{1,n}^\top \mathbf{S}_{1,n} & \text{for } n \in \{2, \dots, N-1\}, \\ \mathbf{M}_{N-1,N}^{-\top} \mathbf{M}_{N-1,N}^{-1} \boxtimes \mathbf{I}_{D_N} & \text{for } n = N. \end{cases} \quad (B.5)$$

## Data availability

A software reproducing the results will be released upon acceptance of the manuscript.

## References

[1] P. Comon, Tensors: A brief introduction, IEEE Signal Process. Mag. 31 (3) (2014) 44–53.

[2] T.G. Kolda, B.W. Bader, Tensor Decompositions and Applications, SIAM Rev. 51 (3) (2009) 455–500.

[3] N.D. Sidiropoulos, L. De Lathauwer, X. Fu, K. Huang, E.E. Papalexakis, C. Faloutsos, Tensor decomposition for signal processing and machine learning, IEEE Trans. Signal Process. 65 (13) (2017) 3551–3582.

[4] E. Acar, C. Aykut-Bingol, H. Bingol, R. Bro, B. Yener, Multiway analysis of epilepsy tensors, Bioinformatics 23 (13) (2007) i10–i18.

[5] Q. Zhao, G. Zhou, L. Zhang, A. Cichocki, S.-I. Amari, Bayesian robust tensor factorization for incomplete multiway data, IEEE Trans. Neural Networks Learn. Syst. 27 (4) (2015) 736–748.

[6] C. Qian, X. Fu, N.D. Sidiropoulos, Y. Yang, Tensor-based channel estimation for dual-polarized massive MIMO systems, IEEE Trans. Signal Process. 66 (24) (2018) 6390–6403.

[7] L. Domanov, Canonical polyadic decomposition of third-order tensors: Reduction to generalized eigenvalue decomposition, SIAM J. Matrix Anal. Appl. 35 (2) (2014) 636–660.

[8] A. Stegeman, N.D. Sidiropoulos, On Kruskal’s uniqueness condition for the candecomp/parafac decomposition, Linear Algebra Appl. 420 (2–3) (2007) 540–552.

[9] N. Li, S. Kindermann, C. Navasca, Some convergence results on the regularized alternating least-squares method for tensor decomposition, Linear Algebra Appl. 438 (2) (2013) 796–812.

[10] I.V. Oseledets, E.E. Tyrtyshnikov, Breaking the curse of dimensionality, or how to use SVD in many dimensions, SIAM J. Sci. Comput. 31 (5) (2009) 3744–3759.

[11] I.V. Oseledets, Tensor-train decomposition, SIAM J. Sci. Comput. 33 (5) (2011) 2295–2317.

[12] Y. Zniyed, R. Boyer, A.F. De Almeida, G. Favier, High-order tensor estimation via trains of coupled third-order CP and tucker decompositions, Linear Algebra Appl. 588 (2020) 304–337.

[13] M. Giraud, V. Itier, R. Boyer, Y. Zniyed, A.L.F. de Almeida, Tucker decomposition based on a tensor train of coupled and constrained CP cores, IEEE Signal Process. Lett. (2023).

[14] J. Bolte, S. Sabach, M. Teboulle, Proximal alternating linearized minimization for nonconvex and nonsmooth problems, Math. Program. (1–2) (2014) 459–494, http://dx.doi.org/10.1007/s10107-013-0701-9.

[15] P. Tseng, Convergence of a block coordinate descent method for nondifferentiable minimization, J. Optim. Theory Appl. 475–494, http://dx.doi.org/10.1023/A: 1017501703105.

[16] H. Cramér, Mathematical Methods of Statistics, Univ. Press, Princeton, 1946.

[17] C.R. Rao, Information and accuracy attainable in the estimation of statistical parameters, Bull. Calcutta Math. Soc 37 (1945) 81–91.

[18] S. Sahnoun, P. Comon, Joint source estimation and localization, IEEE Trans. Signal Process. 63 (10) (2015) 2485–2595.

[19] X. Liu, N. Sidiropoulos, Cramér-Rao lower bounds for low-rank decomposition of multidimensional arrays, IEEE Trans. Signal Process. 49 (9) (2001) 2074–2086.

[20] M. Boizard, R. Boyer, G. Favier, J. Cohen, P. Comon, Performance estimation for tensor CP decomposition with structured factors, in: Proc. ICASSP, 2015.

[21] R. Cabral Farias, J. Cohen, P. Comon, Exploring multimodal data fusion through joint decompositions with flexible couplings, IEEE Trans. Signal Process. 64 (18) (2016) 4830–4844.

[22] C. Ren, R. Cabral Farias, P.-O. Amblard, P. Comon, Performance bounds for coupled models, in: Proc. 2016 IEEE SAM, 2016, event-place: Rio de Janeiro, Brazil.

[23] P. Stoica, T.L. Marzetta, Parameter estimation problems with singular information matrices, IEEE Trans. SP 49 (1) (2001) 87–90.

[24] P. Stoica, B.C. Ng, On the Cramér-Rao bound under parametric constraints, IEEE SP Lett. 5 (7) (1998) 177–179.

[25] T.J. Moore, B.M. Sadler, R.J. Kozick, Maximum-likelihood estimation, the Cramér-Rao bound, and the method of scoring with parameter constraints, IEEE Trans. SP 56 (3) (2008) 895–908.

[26] T. Menni, E. Chaumette, P. Larzabal, J.P. Barbot, New results on Deterministic Cramér-Rao bounds for real and complex parameters, IEEE Trans. SP 60 (3) (2012) 1032–1049.

[27] T. Menni, J. Galy, E. Chaumette, P. Larzabal, Versatility of constrained CRB for system analysis and design, IEEE Trans. AES 50 (3) (2014) 1841–1863.

[28] D.C. Araújo, T. Maksymyuk, A.L.F. de Almeida, T. Maciel, J.C.M. Mota, M. Jo, Massive MIMO: survey and future research topics, IET Commun. 10 (15) (2016) 1938–1946.

[29] Y. Zniyed, R. Boyer, A. De Almeida, G. Favier, Tensor train representation of MIMO channels using the JIRAFE method, Signal Process. 171 (2020) 107479.

[30] L. De Lathauwer, D. Nion, Decompositions of a higher-order tensor in block terms–Part III: Alternating least squares algorithms, SIAM J. Matrix Anal. Appl. 30 (3) (2008) 1067–1083.

[31] W. Pu, S. Ibrahim, X. Fu, M. Hong, Stochastic mirror descent for low-rank tensor decomposition under non-euclidean losses, IEEE Trans. Signal Process. 70 (2022) 1803–1818.

[32] X. Fu, S. Ibrahim, H.-T. Wai, C. Gao, K. Huang, Block-randomized stochastic proximal gradient for low-rank tensor factorization, IEEE Trans. Signal Process. 68 (2020) 2170–2185.

[33] T.G. Kolda, D. Hong, Stochastic gradients for large-scale tensor decomposition, SIAM J. Math. Data Sci. 2 (4) (2020) 1066–1095.

[34] R. Bartels, G. Stewart, Solution of the matrix equation AX+XB=C, Commun. ACM 15 (9) (1972) 820–826.

[35] G. Golub, S. Nash, C.V. Loan, A Hessenberg-Schur method for the problem AX+XB=C, IEEE Trans. Autom. Control. 24 (6) (1979) 909–913.

[36] V. Simoncini, Computational methods for linear matrix equations, SIAM Rev. 58 (3) (2016) 377–441.

[37] J. Gorman, A. Hero, Lower bounds for parametric estimation with constraints, IEEE Trans. Inform. Theory 36 (6) (1990) 1285–1301.

[38] E. Lehmann, G. Casella, Theory of Point Estimation, second ed., Springer, 1998.

[39] D. Slepian, Estimation of signal parameters in the presence of noise, Trans. IRE Prof. Group Inf. Theory 3 (3) (1954) 68–69.

[40] J. Munkres, Algorithms for the assignment and transportation problems, J. Soc. Ind. Appl. Math. 5 (1) (1957) 32–38.

[41] Y.C. Eldar, Minimum variance in biased estimation: Bounds and asymptotically optimal estimators, IEEE Trans. Signal Process. 52 (7) (2004) 1915–1930.

[42] A.O. Hero, A Cramér-Rao Type Lower Bound for Essentially Unbiased Parameter Estimation, Tech. Rep., Massachussets Inst. of tech., Lexington Lincoln Lab, 1992.

[43] A.O. Hero, J.A. Fessler, M. Usman, Exploring estimator bias-variance tradeoffs using the uniform CR bound, IEEE Trans. Signal Process. 44 (8) (1996) 2026–2041.

[44] E. Ollila, V. Koivunen, J. Eriksson, On the cramér-rao bound for the constrained and unconstrained complex parameters, in: 2008 5th IEEE Sensor Array and Multichannel Signal Processing Workshop, IEEE, 2008, pp. 414–418.

[45] M. Sørensen, L. De Lathauwer, Blind signal separation via tensor decomposition with vandermonde factor: Canonical polyadic decomposition, IEEE Trans. Signal Process. 61 (22) (2013) 5507–5519.

[46] E. Sanchez, B.R. Kowalski, Tensorial resolution: a direct trilinear decomposition, J. Chemom. 4 (1) (1990) 29–45.

[47] R. Boyer, P. Comon, Rectified ALS algorithm for multidimensional harmonic retrieval, in: 2016 IEEE SAM Signal Processing Workshop, IEEE, 2016, pp. 1–5.