# Tensor Block-Block Terms Decomposition for Matrix-Valued Imaging Applications

Saulo Cardoso Barreto, Julien Flamant, Sebastian Miron, David Brie

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

## To cite this version:

Saulo Cardoso Barreto, Julien Flamant, Sebastian Miron, David Brie. Tensor Block-Block Terms Decomposition for Matrix-Valued Imaging Applications. 2025. hal-05059862v2

HAL Id: hal-05059862

https://hal.science/hal-05059862v2

Preprint submitted on 2 Jul 2025

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Tensor Block-Block Terms Decomposition
for Matrix-Valued Imaging Applications*

Saulo Cardoso Barreto$^{a,,*}$, Julien Flamant$^a$, Sebastian Miron$^a$, David Brie$^a$

$^a$Université de Lorraine, CNRS, CRAN, F-54000 Nancy, France

**Abstract**

Matrix-valued images appear in many applications, ranging from polarimetric remote sensing to medical imaging. Such images can be represented as 4th-order tensors, where the first two dimensions correspond to spatial variables and the last two encode the matrix feature in each pixel. To efficiently analyze, decompose, and process these images, this paper introduces the block-block terms decomposition (BBTD), a novel low-rank tensor decomposition model that extends bilinear matrix factorization to 4th-order tensors by representing the latter as the sum of outer products of low-rank matrix blocks. Low-rank assumptions permit to greatly reduce the number of parameters to be estimated and enable the enforcement of key physical constraints on matrix sources. We establish both necessary and sufficient conditions for the uniqueness of the BBTD model. To enable the use of BBTD in covariance matrix-valued imaging, we develop an optimization framework that allows efficient handling of non-negativity and symmetry constraints together with low-rank assumptions on matrix blocks. This leads to a fast and computationally efficient algorithm. Numerical experiments on synthetic and real data from DTI illustrate the potential of the BBTD model in matrix-valued imaging.

Keywords: higher-order tensor decomposition, block-block terms decomposition, matrix-valued imaging

1. Introduction

Matrix-valued imaging is an emerging modality which consists in acquiring, for each pixel, a matrix encoding physical features that characterize the observed scene or medium. Important examples are optical polarimetric imaging [1, 2], where Mueller matrices permits a complete characterization of polarization properties of the medium in each pixel. Another application is polarimetric synthetic aperture radar (POL-SAR) imaging [3], where images are matrix-valued, each matrix corresponding to a covariance matrix built

*Research supported in part by French National Research Agency through projects RICOCHET-ANR-21-CE48-0013, ATEMPORAL-ANR-23-CE48-0007. For the purpose of Open Access, a CC-BY public copyright license has been applied by the authors to the present document and will be applied to all subsequent versions up to the Author Accepted Manuscript arising from this submission.

*Corresponding author

Email addresses: sbarreto.cardoso@gmail.com (Saulo Cardoso Barreto), julien.flamant@cnrs.fr (Julien Flamant), sebastian.miron@univ-lorraine.fr (Sebastian Miron), david.brie@univ-lorraine.fr (David Brie)
from scattering vector responses of the target. In medical imaging, a popular MRI technique is Diffusion Tensor Imaging (DTI), where at each pixel, a symmetric matrix characterizes the anisotropy of tissues [4] by quantifying the diffusion of water molecules in the brain. Last, the concept of Fluorescence Excitation-Emission-Matrix (FEEM) imaging has recently been proposed [5], which expands standard FEEM point-wise measurements to extended spatial characterization. In all of these applications, being able to analyze, decompose, or estimate latent factors from matrix-valued images with tailored data analysis tools is crucial to the understanding of underlying physical phenomena.

Matrix-valued images can be naturally represented as multidimensional arrays with four dimensions, or 4th-order tensors, where the first two dimensions encode *spatial dimensions* (e.g., the horizontal and vertical axes of an image) and where the last two dimensions represent that of the matrix feature, which is application dependent (e.g., 4 × 4 Mueller matrices in polarimetric imaging). Despite this natural algebraic representation, the topic has drawn little interest from the tensor community, and relevant tensor algebra tools remain to be developed. On the other hand, the specific case of covariance matrix-valued images (as in DTI or PolSAR) has attracted much interest owing to the interpretation of the space of covariance matrices as a Riemannian manifold. As such, it benefits from the development of Riemannian statistical learning tools, such as classification [6], clustering [7], geodesic PCA [8], Gaussian distributions [9], to cite only a few. In the context of DTI, several interpolation and smoothing methods have also been proposed [10, 11, 12]; see also [13] for a review of manifold-based image processing tools for covariance matrix data.

Tensor decompositions are ubiquitous for high-dimensional data analysis [14], as they allow to break down the complexity of analyses and methods by enforcing low-rank models while preserving the multidimensional structure of the data. Tensor methods are particularly suited for blind source separation. For instance, the standard Canonical Polyadic Decomposition (CPD) [15] leads to direct identification of chemical components in FEEM [16]. Another popular one is the Block Terms Decomposition (BTD) into rank-$(L_r, L_s, 1)$ terms [17] (also noted BTD-LL1), which can be tailored to tackle efficiently the linear unmixing problem in hyperspectral imaging [18, 19, 20]. However, in the context of matrix-valued imaging, these models suffer from several limitations: i) the CPD model enforces rank-one spatial maps – far from typical spatial structures encountered in imaging applications and ii) the BTD-LL1 model has been primarily developed for 3rd-order tensors.

This work introduces the *block-block terms decomposition* (BBTD), a novel tensor decomposition that leverages the representation of matrix-valued images as 4th-order tensors. It addresses the limitations of existing tensor decompositions by enabling the separate low-rank modelling of spatial maps (first block) and matrix features (second block). The resulting model can be viewed as a generalization of the usual linear mixing model (encountered in hyperspectral imaging, for instance), where low-rank assumptions permit: i) to reduce drastically the dimensionality of the problem and ii) to enforce key physical constraints on the matrix sources (e.g., in polarization imaging, a fully polarized source corresponds to a rank-one covariance
matrix [21]). Beyond matrix-valued imaging, which motivated its development, the BBTD model holds its own theoretical interest. Therefore, the first contribution of this paper lies in the theoretical study of the BBTD model, including trivial ambiguities, uniqueness conditions and equivalent representations in terms of existing tensor decompositions. The second contribution consists of tailoring this new model to the important case of covariance matrix imaging. This leads to the development of two algorithms (referred to as *exact* and *inexact*) that perform efficient constrained alternating optimization of each block, imposing low-rank and non-negativity or symmetry constraints simultaneously. Numerical experiments on synthetic and real DTI data demonstrate the potential of the proposed model for matrix-valued imaging applications.

The paper is organized as follows. Section 2 sets the notation and reviews relevant tensor decompositions. Section 3 introduces the BBTD model and discusses its first properties. In particular, we establish uniqueness conditions for this new tensor model, as well as an alternating least-squares algorithm that relies on the CONFAC/PARALIND interpretation of BBTD. Section 4 then describes the practical use of BBTD in covariance matrix imaging applications through a dedicated algorithmic framework. Section 5 provides numerical experiments on synthetic and real-data which illustrate the potential of the approach. Finally, Section 6 summarizes our contributions and outlines future research directions, while technical details are gathered in the Appendix.

## 2. Preliminaries

### 2.1. Notations

Let $\mathbb{R}$ be the set of real numbers and $\mathbb{C}$ the set of complex numbers, and let $\mathbb{F} = \mathbb{R}$ or $\mathbb{C}$ when the distinction is irrelevant. Vectors and matrices with entries in $\mathbb{F}$ are denoted in boldface using lowercase symbols for vectors (e.g., $\mathbf{v} \in \mathbb{F}^I$) and uppercase symbols for matrices (e.g., $\mathbf{M} \in \mathbb{F}^{I \times J}$). A column vector of ones of length $I$ is denoted by $\mathbf{1}_I$. The transpose, conjugate and conjugate-transpose (Hermitian) of a matrix $\mathbf{M}$ are represented as $\mathbf{M}^\top$, $\overline{\mathbf{M}}$, and $\mathbf{M}^\dagger$, respectively. Its inverse and pseudo-inverse are denoted as $\mathbf{M}^{-1}$ and $\mathbf{M}^\dagger$, respectively. Tensors of higher orders are represented with calligraphic letters, e.g., $\mathcal{T} \in \mathbb{F}^{I_1 \times \cdots \times I_N}$. The identity matrix of size $M$ is denoted $\mathbf{I}_M$ while the identity tensor of order $N$ and size $M$ in all of its modes is represented as $\mathcal{I}_{N,M}$. The symbol $\|\cdot\|_F$ denotes the usual Frobenius norm for matrices and tensors.

### 2.2. Basic matrix and tensor operations

The Kronecker product of two matrices $\mathbf{A} \in \mathbb{F}^{I \times J}$ and $\mathbf{B} \in \mathbb{F}^{K \times L}$ is denoted $\mathbf{A} \otimes \mathbf{B}$ such that

$$ \mathbf{A} \otimes \mathbf{B} = \begin{bmatrix} a_{11} \mathbf{B} & a_{12} \mathbf{B} & \cdots & a_{1J} \mathbf{B} \\ a_{21} \mathbf{B} & a_{22} \mathbf{B} & \cdots & a_{2J} \mathbf{B} \\ \vdots & \vdots & \ddots & \vdots \\ a_{I1} \mathbf{B} & a_{I2} \mathbf{B} & \cdots & a_{IJ} \mathbf{B} \end{bmatrix} \in \mathbb{F}^{IK \times JL}. \quad (1) $$
The Khatri-Rao product for matrices $A \in \mathbb{F}^{I \times K}$ and $B \in \mathbb{F}^{J \times K}$ is denoted by $A \odot B$ and defined as the column-wise Kronecker product

76
77

$$
A \odot B = [a_1 \otimes b_1 \ a_2 \otimes b_2 \ \dots \ a_K \otimes b_K] \in \mathbb{F}^{IJ \times K}. \qquad (2)
$$

For two tensors $\mathcal{T} \in \mathbb{F}^{I_1 \times \cdots \times I_N}$ and $\mathcal{S} \in \mathbb{F}^{J_1 \times \cdots \times J_M}$ their outer product is $\mathcal{Z} = \mathcal{T} \circ \mathcal{S} \in \mathbb{F}^{I_1 \times \cdots \times I_N \times J_1 \times \cdots \times J_M}$ with indices

78
79

$$
\mathcal{Z}_{i_1 \dots i_N j_1 \dots j_M} = \mathcal{T}_{i_1 \dots i_N} \mathcal{S}_{j_1 \dots j_M}. \qquad (3)
$$

The fiber of a tensor is a generalization of the concept of rows and columns in matrices. For a tensor of order $N$, a mode-$n$ fiber is a 1-dimensional subset of the tensor obtained by varying only the $n$-th index, with the remaining indices fixed. A tensor slice is obtained by fixing two dimensions in a tensor. The matricization of a tensor $\mathcal{T} \in \mathbb{F}^{I_1 \times \cdots \times I_N}$ along its mode-$n$ will be denoted as the $\mathcal{T}_{(n)} \in \mathbb{F}^{I_n \times (I_1 \cdots I_{n-1} I_{n+1} \cdots I_N)}$ where the ordering of tensor elements follows the convention of Kolda and Bader [22]. The mode-$n$ product of a tensor $\mathcal{T} \in \mathbb{R}^{I_1 \times \cdots \times I_N}$ with a matrix $M \in \mathbb{R}^{J \times I_n}$ yields a new tensor $\mathcal{S} = \mathcal{T} \bullet_n M \in \mathbb{F}^{I_1 \times \cdots \times I_{n-1} \times J \times I_{n+1} \times \cdots \times I_N}$ such that mode-$n$ unfoldings are transformed as $\mathcal{S}_{(n)} = M \mathcal{T}_{(n)}$.

80
81
82
83
84
85
86

2.3. Canonical polyadic decomposition

A tensor $\mathcal{T} \in \mathbb{F}^{I_1 \times \cdots \times I_N}$ of order $N \ge 3$ is said to follow a Canonical Polyadic Decomposition (CPD) if it
can be written as the sum of rank-one tensors

87
88
89

$$
\mathcal{T} = \sum_{r=1}^{R} z_{r}^{(1)} \circ z_{r}^{(2)} \circ \dots \circ z_{r}^{(N)}. \qquad (4)
$$

The rank $R$ denotes the minimal number of terms in (4) for which equality holds. By stacking individual vectors in factors matrices $\mathbf{Z}^{(n)} = [\mathbf{z}_{1}^{(n)} \cdots \mathbf{z}_{R}^{(n)}] \in \mathbb{F}^{I_n \times R}$ the CPD is written compactly as $\mathcal{T} = [[\mathbf{Z}^{(1)}, \dots, \mathbf{Z}^{(N)}]]_{\text{CPD}}$.

90
91
92
93
94
95

$$
\mathcal{T}_{(n)} = \mathbf{Z}^{(n)} (\mathbf{Z}^{(N)} \odot \dots \odot \mathbf{Z}^{(n+1)} \odot \mathbf{Z}^{(n-1)} \odot \dots \odot \mathbf{Z}^{(1)})^{\top}. \quad (5)
$$

A fundamental property of the CPD is its essential uniqueness, i.e., uniqueness up to trivial scaling and
permutation of columns of the factor matrices, provided certain rank conditions are met. One of the most
well-known conditions ensuring the essential uniqueness of the CPD (4) is Kruskal’s condition [23, 24]:

$$
\sum_{n=1}^{N} k_{\mathcal{Z}^{(n)}} \geq 2R + (N-1), \qquad (6)
$$

where the *k*-rank of a matrix *M*, denoted as *k*<sub>*M*</sub>, is the largest integer *r* such that every subset of *r* columns
of *M* is linearly independent.

96
97
## 2.4. Block Terms Decomposition for 3rd-order tensors

For 3rd-order tensors, a popular decomposition in signal and image processing is the Block Terms Decomposition (BTD) into rank-(L<sub>r</sub>, L<sub>r</sub>, 1) terms [25, 17]. Formally, a tensor T ∈ F<sup>I<sub>1</sub>×I<sub>2</sub>×I<sub>3</sub></sup> admits a BTD into rank-(L<sub>r</sub>, L<sub>r</sub>, 1) terms if it can be expressed as

$$ \mathcal{T} = \sum_{r=1}^{R} (\mathbf{Z}_{r}^{(1)} \mathbf{Z}_{r}^{(2)\top}) \circ \mathbf{z}_{r}^{(3)}, \qquad (7) $$

where R is the rank of the decomposition and where for each 1 ≤ r ≤ R, we have Z<sub>r</sub><sup>(1)</sup> ∈ F<sup>I<sub>1</sub>×L<sub>r</sub></sup>, Z<sub>r</sub><sup>(2)</sup> ∈ F<sup>I<sub>2</sub>×L<sub>r</sub></sup>, and z<sub>r</sub><sup>(3)</sup> ∈ F<sup>I<sub>3</sub></sup>. In the remainder of this paper, we refer to the decomposition (7) and its trivially permuted variants simply as BTD-LL1. Such decomposition is particularly useful in tensor imaging applications when two of the dimensions correspond, for example, to spatial coordinates. In such cases, the latent factors obtained from the classical CPD lack direct interpretability. In contrast, BTD-LL1 provides interpretable low-rank representations of images with uniqueness guarantees. One of the most successful application domains of BTD-LL1 is hyperspectral imaging, where this tensor model proves especially well-suited [19, 18, 20].

Crucial to the results derived in this paper is the interpretation of the BTD-LL1 model (7) as a CPD (4) with colinearities. Let L<sub>r</sub> = ∑<sub>r=1</sub><sup>R</sup> L<sub>r</sub> and build factor matrices Z<sup>(1)</sup> = [Z<sub>1</sub><sup>(1)</sup> ... Z<sub>R</sub><sup>(1)</sup>] ∈ F<sup>I<sub>1</sub>×L</sup> and Z<sup>(2)</sup> = [Z<sub>1</sub><sup>(2)</sup> ... Z<sub>R</sub><sup>(2)</sup>] ∈ F<sup>I<sub>2</sub>×L</sup> by horizontally stacking individual block factors in (7). Define Z<sup>(3)</sup> = [z<sub>1</sub><sup>(3)</sup> ... z<sub>R</sub><sup>(3)</sup>] ∈ F<sup>I<sub>3</sub>×R</sup>. Then classical calculations [26] show that

$$ \mathcal{T} = \sum_{r=1}^{R} (\mathbf{Z}_{r}^{(1)} \mathbf{Z}_{r}^{(2)\top}) \circ \mathbf{z}_{r}^{(3)} = \llbracket \mathbf{Z}^{(1)}, \mathbf{Z}^{(2)}, \mathbf{Z}^{(3)} \mathbf{\Theta} \rrbracket_{\text{CPD}}, \qquad (8) $$

where Θ = bdiag(1<sub>L<sub>1</sub></sub><sup>T</sup>, ..., 1<sub>L<sub>R</sub></sub><sup>T</sup>) is the constraint matrix associated to ranks L<sub>1</sub>, ..., L<sub>R</sub> in the BTD-LL1 model (7). Eq. (8) is also sometimes referred to as the PARALIND or CONFAC-(1,3) interpretation [26] of BTD-LL1.

Similarly to CPD, the BTD-LL1 model is unique up to certain trivial ambiguities: permutation of the R terms, linear transformations of the column spaces of Z<sub>r</sub><sup>(1)</sup> and Z<sub>r</sub><sup>(2)</sup>, as well as scalar scaling and counterscaling of the block {Z<sub>r</sub><sup>(1)</sup> Z<sub>r</sub><sup>(2)\top</sup>} and factor z<sub>r</sub><sup>(3)</sup> in each r-th term.

The uniqueness analysis of BTD-LL1 is based on an extension of the classical CPD uniqueness results to the case of partitioned factor matrices. Several sufficient uniqueness conditions are presented in [17]: for instance, the decomposition T = [Z<sup>(1)</sup>, Z<sup>(2)</sup>, Z<sup>(3)</sup>Θ]_{CPD} is essentially unique if:

$$ I_1 I_2 \geq \sum_{r=1}^{R} L_r^2 \quad \text{and} \quad k'_{Z^{(1)}} + k'_{Z^{(2)}} + k_{Z^{(3)}} \geq 2R + 2, \qquad (9) $$

in which the k'-rank (k'<sub>M</sub>), for a partitioned matrix M, is the largest r for which any set of r submatrices of M forms a linearly independent set of columns.
### 3. Block-block terms decomposition

This section introduces a novel decomposition for 4th-order tensors, named as block-block terms decomposition (BBTD). This model generalizes the classical bilinear factorization of matrices to 4th-order tensors. More precisely, while matrices can be decomposed as a sum of outer products of vectors, the BBTD model decomposes a 4th-order tensor as the sum of outer products of low-rank matrices – the so-called block-block terms. To the best of our knowledge, this decomposition is novel, and while its introduction is strongly motivated by applications to matrix-valued imaging, it also holds an independent theoretical interest.

This section is organized as follows: Section 3.1 introduces formally the BBTD decomposition for 4th-order tensors and studies first properties related to trivial ambiguities and dimension reduction. Section 3.2 then develops uniqueness conditions for BBTD, providing a necessary condition and a sufficient condition for identifiability. Section 3.3 discusses the interpretation of BBTD as a PARALIND / CONFAC model, which is useful for the design of an alternating least squares (ALS) algorithm for BBTD in Section 3.4.

#### 3.1. BBTD: definition, ambiguities and number of parameters

Let $\mathcal{T} \in \mathbb{F}^{I_1 \times I_2 \times I_3 \times I_4}$ be a 4th-order tensor. We say that $\mathcal{T}$ admits a block-block terms decomposition (BBTD) if it can be represented as a sum of outer products of matrix blocks $\{U_r\}$ and $\{V_r\}$ such that

$$ \mathcal{T} = \sum_{r=1}^{R} U_r \circ V_r, \text{ where rank } U_r = L_{1r}, \text{ rank } V_r = L_{2r}, r = 1, \dots, R. \quad (10) $$

We refer to the triplet $(R, \{L_{1r}\}, \{L_{2r}\})$ as the parameters of the BBTD model. Note that, by definition of the matrix rank, one has $L_{1r} \le \min(I_1, I_2)$ and $L_{2r} \le \min(I_3, I_4)$ for $r = 1, \dots, R$. An alternative equivalent definition for BBTD consists in parameterizing the r-th block-block term $U_r \circ V_r$ explicitly. Let $U_r = A_r B_r^\top$ with $A_r \in \mathbb{F}^{I_1 \times L_{1r}}$ and $B_r \in \mathbb{F}^{I_2 \times L_{2r}}$. Similarly, write $V_r = C_r D_r^\top$ where $C_r \in \mathbb{F}^{I_3 \times L_{2r}}$ and $D_r \in \mathbb{F}^{I_4 \times L_{2r}}$. Then the BBTD model (10) can be equivalently expressed as

$$ \mathcal{T} = \sum_{r=1}^{R} A_r B_r^\top \circ C_r D_r^\top = \llbracket A, B, C, D \rrbracket_{\text{BBTD}}, \quad (11) $$

where we have introduced the factors matrices $A, B, C$ and $D$ of the BBTD model such that

$$ A = [A_1 \cdots A_R] \in \mathbb{F}^{I_1 \times L_1}, \quad (12) $$

$$ B = [B_1 \cdots B_R] \in \mathbb{F}^{I_2 \times L_1}, \quad (13) $$

$$ C = [C_1 \cdots C_R] \in \mathbb{F}^{I_3 \times L_2}, \quad (14) $$

$$ D = [D_1 \cdots D_R] \in \mathbb{F}^{I_4 \times L_2}, \quad (15) $$

where $L_1 = \sum_{r=1}^{R} L_{1r}$ and $L_2 = \sum_{r=1}^{R} L_{2r}$. This parametrization will be useful later for establishing uniqueness conditions, studying equivalent models, and practical algorithm development.
**Trivial ambiguities** The BBTD model (11), like any tensor decomposition model, exhibits several trivial ambiguities. Specifically, the factor matrices of the BBTD model can only be reconstructed up to block-block term permutations, block-wise scaling, and up to subspaces associated with individual blocks {$U_r = A_r B_r^\top$} and {$V_r = C_r D_r^\top$}. More precisely, let $\mathcal{T} = [A, B, C, D]_{\text{BBTD}}$ a BBTD model and consider an alternative BBTD decomposition $\mathcal{T} = [\tilde{A}, \tilde{B}, \tilde{C}, \tilde{D}]_{\text{BBTD}}$. We say that the two decompositions are trivially related to one another if their factors can be related through a combination of the following transformations:

* Block-block terms permutation: Permutation $\sigma$ of the $R$ block-block terms such that $(\tilde{U}_r, \tilde{V}_r) = \sigma((U_r, V_r))$, with $1 \le r \le R$.

* Block scaling: Scaling/counterscaling of the blocks of each term by scaling factors $\lambda_1, \dots, \lambda_R \in \mathbb{F} \setminus \{0\}$ such that $\tilde{U}_r = \lambda_r U_r$ and $\tilde{V}_r = \lambda_r^{-1} V_r$, with $1 \le r \le R$.

* Block subspace linear transformations: Transformations of the column subspaces of the matrices of each block by two non-singular matrices $P_r \in \mathbb{F}^{L_{1r} \times L_{1r}}$ and $Q_r \in \mathbb{F}^{L_{2r} \times L_{2r}}$ such that $\tilde{A}_r = A_r P_r$, $\tilde{B}_r = B_r (P_r^\top)^{-1}$, $\tilde{C}_r = C_r Q_r$ and $\tilde{D}_r = D_r (Q_r^\top)^{-1}$.

These ambiguities motivate the definition of essential uniqueness of the model BBTD presented in the next section.

**Number of parameters** A key benefit of the BBTD model (11) lies in its ability to greatly reduce the number of parameters to be estimated in $\mathcal{T}$. Indeed, for $\mathcal{T} \in \mathbb{F}^{I_1 \times I_2 \times I_3 \times I_4}$ one has $n_{\text{tensor}} = I_1 I_2 I_3 I_4$ variables. On the other hand, the BBTD model breaks down this exponential complexity by linearizing the dependence in terms of the tensor's dimensions. The number of free parameters in the BBTD model (11) is

$$n_{\text{BBTD}} = \sum_{r=1}^{R} [L_{1r}(I_1 + I_2 - L_{1r}) + L_{2r}(I_3 + I_4 - L_{2r})] - R. \quad (16)$$

Note that a naive computation (i.e., without taking into account the trivial ambiguities of the BBTD model) would yield $n'_{\text{BBTD}} = \sum_{r=1}^{R} [L_{1r}(I_1 + I_2) + L_{2r}(I_3 + I_4)]$. To obtain (16) one must subtract from $n'_{\text{BBTD}}$ the number of parameters that must be fixed to eliminate the block scaling and block subspace indeterminacies (permutations play no role in this context). Consider the $r$-th block-block term. To resolve the subspace linear ambiguities in the first block $A_r B_r^\top$, a total of $L_{1r}^2$ parameters must be fixed. This is equivalent to imposing the transformation matrix $P_r$ above. Consequently, the number of free parameters for this first block is given by $I_1 L_{1r} + I_2 L_{1r} - L_{1r}^2$. Similarly, for the second block $C_r D_r^\top$, the number of free parameters is $I_3 L_{2r} + I_4 L_{2r} - L_{2r}^2$. Additionally, to remove the scaling ambiguities between the two blocks, one extra parameter per block-block term must be fixed. This is equivalent to fixing the $\lambda_r$ scaling parameter above. This results in a total of $L_{1r}(I_1 + I_2 - L_{1r}) + L_{2r}(I_3 + I_4 - L_{2r}) - 1$ free parameters per block-block term. Summing over all $R$ terms one gets (16).
Eventually, note that the BBTD model can achieve remarkable compression. Indeed, for a tensor $\mathcal{T}$ with $I_1 = I_2 = I_3 = I_4 = 100$, BBTD with constant $(L_{1r}, L_{2r}) = (10, 10)$ leads to $n_{\text{BBTD}} = 3799R$, where $R$ is the number of block-block terms. For $R = 5$, the number of free parameters has thus been reduced by more than 5000 times using BBTD.

## 3.2. Uniqueness of BBTD

The trivial ambiguities of the BBTD model – see previous section – motivate the following definition for the essential uniqueness of the BBTD model.

**Definition 1** (Essential Uniqueness of BBTD). A BBTD decomposition $\mathcal{T} = [[A, B, C, D]]_{\text{BBTD}}$ is said to be essentially unique if any alternative BBTD decomposition $\mathcal{T} = [[\tilde{A}, \tilde{B}, \tilde{C}, \tilde{D}}]]_{\text{BBTD}}$ differs from the original only through a combination of trivial ambiguities: block-block terms permutation; block scaling; and block subspace linear transformations.

In the following, we will refer to the uniqueness of BBTD in the sense given by Definition 1. A necessary condition for uniqueness can be directly derived from the relationship between the number of data points and the number of free parameters in the BBTD model. The number of free parameters in the BBTD model must necessarily be less than or equal to the number of available data points. This ensures that the model is not over-parameterized, which is essential for achieving a unique decomposition. This result is summarized by the following proposition.

**Proposition 1** (Necessary uniqueness condition for BBTD). Suppose that $\mathcal{T} = [[A, B, C, D]]_{\text{BBTD}}$ is unique. Then $n_{\text{BBTD}} \le n_{\text{tensor}}$, where $n_{\text{BBTD}}$ is given by (16) and $n_{\text{tensor}} = I_1 I_2 I_3 I_4$.

This result is useful for obtaining upper bounds on the attainable rank for a given BBTD model. Consider the particular case, common in practice, where $L_{1r} = L_1$ and $L_{2r} = L_2$ for all $r = 1, \dots, R$. In this case, the necessary condition of Proposition 1 simplifies to:

$$R \le \left[ \frac{I_1 I_2 I_3 I_4}{L_1(I_1 + I_2 - L_1) + L_2(I_3 + I_4 - L_2) - 1} \right]. \quad (17)$$

Thus, for a tensor with equal dimensions $I_1 = I_2 = I_3 = I_4 = 10$ and subspace sizes $L_1 = L_2 = 5$, the maximum attainable rank according to this necessary condition is $R = 67$. We will see for this same example that this upper bound is quite optimistic, and far from what sufficient uniqueness condition for BBTD can guarantee.

The uniqueness properties of the BBTD model are closely related to those of the BTD-LL1 model. In fact, by vectorizing either the first two or the last two modes of (11), one obtains a BTD-LL1 decomposition. We leverage this relationship to establish sufficient uniqueness conditions for the proposed BBTD decomposition. However, it is important to highlight that, due to this additional flexibility, BBTD generally exhibits more favorable uniqueness properties compared to BTD-LL1. Thus, deriving sufficient uniqueness conditions for
BBTD relies on two key principles: (i) interpreting its 3rd-order unfoldings as BTD-LL1 decompositions and (ii) applying state-of-the-art BTD-LL1 uniqueness results to these unfoldings and selecting the most favorable one. We first establish the equivalent BTD-LL1 models for BBTD, followed by an analysis of their properties and identifiability results.

**Proposition 2.** Let $\mathcal{T} = [[A, B, C, D]]_{\text{BBTD}}$. Then, it is possible to partially unfold $\mathcal{T}$ as a 3rd-order tensor by combining two different modes. In particular, combining the first two modes, one obtains a 3rd-order tensor $\mathcal{T}^{1:2} \in \mathbb{F}^{I_1 I_2 \times I_3 \times I_4}$ such that

$$ \mathcal{T}^{1:2} = \sum_{r=1}^{R} \mathrm{vec}(A_r B_r^\top) \circ C_r D_r^\top, \qquad (18) $$

i.e., $\mathcal{T}^{1:2}$ follows a BTD into $R$ rank-(1, $L_{2r}$, $L_{2r}$) terms. Alternatively, combining the two last modes results in the 3rd-order tensor $\mathcal{T}^{3:4} \in \mathbb{F}^{I_1 \times I_2 \times I_3 I_4}$ such that

$$ \mathcal{T}^{3:4} = \sum_{r=1}^{R} A_r B_r^\top \circ \mathrm{vec}(C_r D_r^\top), \qquad (19) $$

i.e., $\mathcal{T}^{3:4}$ admits a BTD into $R$ rank-($L_{1r}$, $L_{1r}$, 1) terms.

*Proof.* The proof is direct by applying vectorizations onto pairs of modes (1, 2) and (3, 4), respectively. □

These 3rd-order interpretations are crucial to study the identifiability properties of the BBTD decomposition of a 4th-order tensor. Theorem 1 below provides a sufficient uniqueness condition in terms of 3rd-order tensorizations.

**Theorem 1** (Sufficient condition for BBTD uniqueness). Let $\mathcal{T} = [[A, B, C, D]]_{\text{BBTD}}$. Suppose that at least one of these conditions is satisfied:

1. the BTD decomposition (18) into rank-(1, $L_{2r}$, $L_{2r}$) terms of $\mathcal{T}^{1:2}$ is essentially unique;

2. the BTD decomposition (19) into rank-($L_{1r}$, $L_{1r}$, 1) terms of $\mathcal{T}^{3:4}$ is essentially unique.

Then the BBTD of $\mathcal{T}$ is unique.

*Proof.* Let $\mathcal{T}$ be decomposed in BBTD form, $\mathcal{T} = [[A, B, C, D]]_{\text{BBTD}}$. By Proposition 2, $\mathcal{T}^{1:2}$ admits a (3rd-order) BTD-LL1 model:

$$ \mathcal{T}^{1:2} = \sum_{r=1}^{R} \mathrm{vec}(A_r B_r^\top) \circ C_r D_r^\top. \qquad (20) $$

Assume that the decomposition (18) of $\mathcal{T}^{1:2}$ is unique up to trivial ambiguities of the BTD-LL1 model. Then, both $\mathrm{vec}(A_r B_r^\top)$ and $(C_r D_r^\top)$ are identifiable. As a result, all blocks $U_r = A_r B_r^\top$ and $V_r = C_r D_r^\top$ in the 4th-order BBTD are uniquely recoverable. Consequently the BBTD of $\mathcal{T}$ into $(A, B, C, D)$ is essentially unique according to Definition 1. Similarly, assuming the essential uniqueness of the decomposition of $\mathcal{T}^{3:4}$ leads to the same conclusion. Thus, the 4th-order BBTD is unique if either $\mathcal{T}^{1:2}$ or $\mathcal{T}^{3:4}$ is essentially unique. □
<figure><img src="image_3.png" /><figcaption>Figure 1: The heatmap provides the maximum attainable ranks *R* that guarantee uniqueness of the BBTD of a tensor *T* ∈ F¹⁰×¹⁰×¹⁰×¹⁰, for various combinations of *L*₁ and *L*₂. Color indicates which tensor unfolding yields the highest *R*: yellow for the vectorization of *U*<sub>r</sub> blocks (T<sup>1:2</sup>) and blue for the vectorization of the *V*<sub>r</sub> blocks (T<sup>3:4</sup>). Red indicates that both models give the same maximum *R*.</figcaption></figure>

The 3rd-order BTD-LL1 equivalent models of BBTD, each associated with distinct uniqueness conditions, provide two complementary perspectives, allowing us to select the model with the most favorable conditions. To illustrate the sufficient uniqueness conditions, we consider the specific case where $L_{1r} = L_1$ and $L_{2r} = L_2$ for all $r = 1, \dots, R$. Based on the BTD-LL1 uniqueness conditions [17] (see also (9) herein for a particular condition), Figure 1 depicts the relationships between the ranks $L_1$ and $L_2$ of the two BBTD blocks and the maximum attainable rank $R$ that ensures essential uniqueness of the BBTD of a tensor $T \in F^{10 \times 10 \times 10 \times 10}$.

The heatmap uses color to indicate the “dominant” unfolding that achieves the highest rank $R$: yellow corresponds to the $T^{1:2}$ unfolding, blue signifies the $T^{3:4}$ unfolding, and red highlights cases where both 3rd-order models yield the same maximum $R$. As expected, given that $I_1 = I_2 = I_3 = I_4$, the heatmap clearly illustrates the two complementary perspectives of the BBTD model. The observed symmetry reflects this duality, with $T^{1:2}$ prevailing for larger values of $L_1$ and $T^{3:4}$ dominating when $L_2$ is higher. Finally, in the top-right region, where both block ranks are large, the two models yield the same maximum $R$.

It is worth noting that when dealing with tensors where the dimensions of the two blocks differ signifi-
cantly, unfolding the block with smaller dimensions generally provides the best uniqueness conditions. For
instance, consider a tensor where $I_1 = I_2 = 75$ and $I_3 = I_4 = 2$. In this scenario, the unfolding $\mathcal{T}^{3:4}$ typically
achieves the highest attainable $R$, whereas the unfolding $\mathcal{T}^{1:2}$ is only relevant when $L_2 = 1$ and $L_1$ is large.
### 3.3. Equivalent models for BBTD: PARALIND/CONFAC and rank-(L<sub>1r</sub>, L<sub>1r</sub>, L<sub>2r</sub>, L<sub>2r</sub>) BTD model

The BBTD model (11) can be viewed as a 4th-order CPD with structured colinearities between the columns of factor matrices, usually referred to as the PARALIND [27] or CONFAC [28] model. This interpretation will be useful to derive a practical algorithm for the BBTD model (11), see Section 3.4. The block-block structure in (11) can be enforced by two constraint matrices $\Phi$ and $\Psi$ such that

$$ \mathcal{T} = \llbracket A, B, C, D \rrbracket_{\text{BBTD}} = \llbracket A\Phi, B\Phi, C\Psi, D\Psi \rrbracket_{\text{CPD}}, \qquad (21) $$

where the constraint matrices are given by

$$ \Phi = \text{bdiag}(I_{L_{11}} \otimes 1_{L_{21}^{\top}}, \dots, I_{L_{1R}} \otimes 1_{L_{2R}^{\top}}), \qquad (22) $$

$$ \Psi = \text{bdiag}(1_{L_{11}^{\top}} \otimes I_{L_{21}}, \dots, 1_{L_{1R}}^{\top} \otimes I_{L_{2R}}). \qquad (23) $$

For the r-th BBTD term, the constraint matrix $\Phi$ reproduces $L_{2r}$ times the columns of $A_r$ and $B_r$, while the constraint matrix $\Psi$ alternately repeats $L_{1r}$ times the columns of $C_r$ and $D_r$ to match the two block terms together. When $L_{1r} = L_{2r} = 1$ for every $r = 1, \dots, R$, the constraint matrices $\Phi = \Psi = I_R$ and the BBTD model degenerates into the standard CPD (4). The PARALIND/CONFAC interpretation of BBTD (21) together with the usual CP decomposition unfoldings (5) enables the derivation of the matrix unfoldings of the BBTD model as:

$$ T_{(1)} = A\Phi(D\Psi \odot C\Psi \odot B\Phi)^{\top} \in \mathbb{F}^{I_1 \times I_2 I_3 I_4}, \qquad (24) $$

$$ T_{(2)} = B\Phi(D\Psi \odot C\Psi \odot A\Phi)^{\top} \in \mathbb{F}^{I_2 \times I_1 I_3 I_4}, \qquad (25) $$

$$ T_{(3)} = C\Psi(D\Psi \odot B\Phi \odot A\Phi)^{\top} \in \mathbb{F}^{I_3 \times I_1 I_2 I_4}, \qquad (26) $$

$$ T_{(4)} = D\Psi(C\Psi \odot B\Phi \odot A\Phi)^{\top} \in \mathbb{F}^{I_4 \times I_1 I_2 I_3}. \qquad (27) $$

Alternatively, the BBTD model (11) can be interpreted as a 4th-order BTD into rank-(L<sub>1r</sub>, L<sub>1r</sub>, L<sub>2r</sub>, L<sub>2r</sub>) terms by adapting 3rd-order definitions introduced in [17]. Namely, exploiting the structure of constraint matrices $\Phi$ and $\Psi$, one has

$$ \mathcal{T} = \llbracket A, B, C, D \rrbracket_{\text{BBTD}} = \sum_{r=1}^{R} \mathcal{D}_r \bullet_1 A_r \bullet_2 B_r \bullet_3 C_r \bullet_4 D_r, \qquad (28) $$

where the r-th core tensor $\mathcal{D}_r$ is explicitly given by $\mathcal{D}_r = \mathcal{I}_{4,L_{1r},L_{2r}} \bullet_1 \Phi_r \bullet_2 \Phi_r \bullet_3 \Psi_r \bullet_4 \Psi_r$, with $\mathcal{I}_{4,L_{1r},L_{2r}}$ being the 4th-order identity tensor of size $L_{1r}L_{2r}$ in all its modes and $\Phi_r = I_{L_{1r}} \otimes 1_{L_{2r}^{\top}}$ and $\Psi_r = 1_{L_{1r}^{\top}} \otimes I_{L_{2r}}$ are deduced from (22) and (23), respectively. It is easy to verify that the core tensor $\mathcal{D}_r$ has multilinear rank $(L_{1r}, L_{1r}, L_{2r}, L_{2r})$. Hence Eq. (28) shows that the BBTD model can be equivalently interpreted as a 4th-order BTD into R rank-(L<sub>1r</sub>, L<sub>1r</sub>, L<sub>2r</sub>, L<sub>2r</sub>) terms with a specific structured core tensor $\mathcal{D}_r$.

Note, however, that, up to our knowledge, these two equivalent interpretations cannot be leveraged to derive identifiability conditions for the BBTD model. Indeed, existing identifiability results [29, 30, 31, 32] for
CONFAC/PARALIND models rely on the absence of colinearities in at least one of the factor matrices. This is not the case for BBTD, where the structural constraints imposed across all modes introduce redundancy, preventing the direct application of these identifiability results. On the other hand, existing results for BTD uniqueness (including rank-(L<sub>r</sub>, L<sub>r</sub>, 1) and rank-(L<sub>r</sub>, M<sub>r</sub>, N<sub>r</sub>) decompositions) [17] are currently limited to 3rd-order tensors, and their extension to 4th-order tensors is not straightforward. This will be addressed in future work. Nonetheless, available results for 3rd-order BTD-LL1 already permit to derive sufficient uniqueness conditions for the BBTD model, as explained in the previous section (see Theorem 1).

### 3.4. Alternating least squares for BBTD

The PARALIND / CONFAC interpretation of BBTD – see Section 3.3 – enables a simple reconstruction strategy of factor matrices based on alternating least squares (ALS). Namely, for a given tensor $\mathcal{T} \in \mathbb{F}^{I_1 \times I_2 \times I_3 \times I_4}$ and ranks $R$, $\{L_{1r}\}$ and $\{L_{2r}\}$ we aim at solving the optimization problem

Thanks to unfoldings (24) – (27), one can extend the classical ALS procedure used in many unconstrained tensor decompositions [33, 34, 26] to solve (29). Namely, the so-called ALS-BBTD algorithm consists of cyclically updating individual factors by solving matrix LS problems obtained by plugging unfoldings expressions in (29). Such an approach enjoys reasonable convergence guarantees (such as convergence to a stationary point of the objective function), is parameter-free, and provides a baseline algorithm for any method attempting to solve the BBTD reconstruction problem. The $(k+1)$-th iteration of ALS-BBTD reads

Updates (30)–(33) are repeated until a suitable convergence criterion is met, such as stopping whenever the decrease of the objective function goes below a pre-defined threshold. In practice, the implementation of least squares solutions (30) – (33) can be efficiently addressed through fast solvers of normal equations, such as the ones implemented in Python’s `scipy.linalg Solve`.

## 4. BBTD for covariance imaging applications

As explained in the introduction, one of the main motivations for the study of the novel Block-Block Terms Decomposition (BBTD) in Section 3 lies in covariance matrix imaging applications. In this setting, images are matrix-valued, i.e., each pixel carries a covariance matrix of dimension $K \times K$, which can be

$$ \min_{A,B,C,D} \| \mathcal{T} - [A, B, C, D]_{\text{BBTD}} \|_{F}^{2}. \qquad (29) $$

$$ A^{k+1} \leftarrow (T_{(1)} M_{1}^{H}) (M_{1} M_{1}^{H})^{-1}, \quad \text{with } M_{1} = \Phi (D^{k} \Psi \odot C^{k} \Psi \odot B^{k} \Phi)^{\top}, \qquad (30) $$

$$ B^{k+1} \leftarrow (T_{(2)} M_{2}^{H}) (M_{2} M_{2}^{H})^{-1}, \quad \text{with } M_{2} = \Phi (D^{k} \Psi \odot C^{k} \Psi \odot A^{k+1} \Phi)^{\top}, \qquad (31) $$

$$ C^{k+1} \leftarrow (T_{(3)} M_{3}^{H}) (M_{3} M_{3}^{H})^{-1}, \quad \text{with } M_{3} = \Psi (D^{k} \Psi \odot B^{k+1} \Phi \odot A^{k+1} \Phi)^{\top}, \qquad (32) $$

$$ D^{k+1} \leftarrow (T_{(4)} M_{4}^{H}) (M_{4} M_{4}^{H})^{-1}, \quad \text{with } M_{4} = \Psi (C^{k+1} \Psi \odot B^{k+1} \Phi \odot A^{k+1} \Phi)^{\top}. \qquad (33) $$
real or complex-valued depending on the application. Important examples include polarization imaging and diffusion tensor imaging. A covariance matrix image can thus be represented as a tensor $\mathcal{T} \in \mathbb{F}^{I \times J \times K \times K}$ where $I \times J$ is the number of pixels. In such applications, a typical assumption is to model the observed tensor $\mathcal{T}$ as a linear superposition of $R$ spatially distributed elementary covariance sources such that

$$ \mathcal{T} \approx \sum_{r=1}^{R} S_r \circ \Sigma_r, \qquad (34) $$

where for every $r = 1, \dots, R$, the matrix $S_r \ge 0 \in \mathbb{R}_+^{I \times J}$ is a non-negative spatial activation factor and where $\Sigma_r \ge 0 \in \mathbb{F}^{K \times K}$ is an elementary covariance matrix (which is positive semidefinite by definition). Since the space of $K \times K$ covariance matrices is a convex cone in $\mathbb{F}^{K \times K}$, these constraints ensure that the resulting tensor (34) can be interpreted as a covariance matrix image. From a statistical viewpoint, this formulation can be seen as a structured generalization of $k$-means clustering, where instead of assigning each pixel to a single centroid, we assign weights to multiple elementary covariances $\Sigma_r$, modulated spatially by the maps $S_r$ for each pixel. The novelty of the model lies in its ability to simultaneously capture spatial variability and second-order statistical dependencies in a compact and interpretable manner.

The linear mixture model (34) also has a natural interpretation in many covariance-imaging applications. For instance, in the Stokes polarization imaging context, it can be interpreted as the superposition of $R$ incoherent sources, such that the $r$-th source has spatial distribution $S_r$ and constant polarization state described by $\Sigma_r$. In DTI, (34) can be interpreted as a decomposition of complex diffusion patterns into a mixture of canonical diffusion profiles, where each pixel of a slice is modeled as a sum of $R$ elementary covariance matrices $\Sigma_r$, each representing a principal diffusion direction associated with the $r$-th component and modulated by its corresponding spatial activation map $S_r$.

In general, the approximation problem (34) is ill-posed, as the decomposition is not unique unless additional constraints are imposed on $S_r$ and $\Sigma_r$. Consider the matricized form of the tensor $\mathcal{T}$, denoted $\mathcal{T} = S\Sigma^\top$, where $S \in \mathbb{R}_+^{I \times J \times R}$ contains the flattened spatial maps and $\Sigma \in \mathbb{C}^{K^2 \times R}$ the vectorized covariance components. If $S > 0$, it is always possible to insert an invertible matrix $M$ with non-negative entries such that $\mathcal{T} = SM^{-1}M\Sigma^\top$, with $SM^{-1} \ge 0$ and $M\Sigma^\top$ defining a new set of covariance terms (by convex conic combination). This reparameterization preserves the data representation but alters the decomposition components, highlighting the existence of many valid combinations ($\{S_r\}, \{\Sigma_r\}$) that explain the same observed $\mathcal{T}$. To cope with this issue, we consider low-rank assumptions of the spatial and covariance factors in (34). To motivate such a viewpoint, Section 4.1 first reviews the relevance of low-rank covariance matrices in imaging applications. Then Section 4.2 deploys an efficient constrained BBTD framework to compute the decomposition (34), while Section 4.3 describes the proposed alternating optimization strategy. Throughout the section, the proposed approach is detailed for the most general case $\mathbb{F} = \mathbb{C}$. The case $\mathbb{F} = \mathbb{R}$ follows directly.
### 4.1. Low-rank covariance matrices in imaging applications

The use of low-rank covariance matrices plays a crucial role in covariance imaging applications. In polarization imaging, for instance, Stokes parameters can be represented using rank-one covariance matrices, which indicate fully polarized states [35]. In Mueller matrix imaging, which is used to characterize the polarization properties of materials, a covariance matrix is derived from the 4×4 Mueller matrix to facilitate further analysis. A key property of the Mueller matrix is that the rank of its corresponding covariance matrix relates to depolarization, which accounts for multiple scattering mechanisms in a medium. A rank-1 covariance matrix characterizes a non-depolarizing medium, while ranks 2, 3, or 4 indicate increasing degrees of depolarization [36]. Additionally, any physical Mueller matrix can be expressed as a weighted sum of non-depolarizing Mueller matrices [37], reinforcing its fundamental structure.

Beyond polarization imaging, low-rank covariance matrices have significant implications for diffusion tensor imaging (DTI), where they capture the anisotropy of water molecule diffusion and their low-rank structure highlights principal directions of strain [38]. A rank-2 or rank-3 structure respectively corresponds to planar or fully anisotropic diffusion (i.e., the diffusion varies in all three directions), which helps, for example, to identify white matter fiber orientations in the brain. Additionally, low-rank approximations can denoise DTI data while preserving dominant diffusion directions or isolate pathological features (e.g., edema) that disrupt normal anisotropy [39].

### 4.2. Constrained BBTD approximation for covariance imaging

Consider the equality case in (34) where in addition spatial factors and covariances matrices are low rank, such that $S_r = A_r B_r^T$, with rank $S_r = L_{1r} \le \min(I, J)$ and $\Sigma_r = C_r C_r^H$, with rank $\Sigma_r = L_{2r} \le K$. Then one observes directly that the covariance matrix image $\mathcal{T}$ follows a particular BBTD model (11) such that $\mathcal{T} = [[A, B, C, \bar{C}]_{\text{BBTD}}$, with two additional interpretability constraints: 1) non-negativity of the spatial maps $\{S_r\}$ and 2) (conjugate) symmetry of the last two modes, which enforces the positive semidefiniteness of elementary covariance matrices $\{\Sigma_r\}$. In particular, this constrained BBTD model enjoys the uniqueness properties of the unconstrained BBTD introduced in Section 3.

In practice, however, for an arbitrary covariance matrix image $\mathcal{T}$, one can only expect to compute an approximation (34) taking the form of a constrained BBTD model. In this work, we propose to handle this approximation problem by solving the constrained optimization problem

$$ \begin{gathered} \min_{\{S_r\}_{r=1}^R, \{\Sigma_r\}_{r=1}^R} \left\| \mathcal{T} - \sum_{r=1}^R S_r \circ \Sigma_r \right\|_F^2 \quad (35) \\ \text{such that} \quad S_r \ge 0, \quad \text{rank } S_r \le L_{1r}, r = 1, \dots, R \\ \Sigma_r \ge 0, \quad \text{rank } \Sigma_r \le L_{2r}, r = 1, \dots, R \end{gathered} $$

Eq. (35) defines a constrained BBTD approximation problem, where constraints ensure the physical interpretability of the solution as well as its uniqueness, upon proper choice of ranks $\{L_{1r}\}$ and $\{L_{2r}\}$. Eventually,
note that compared to its unconstrained counterpart (29), the constrained BBTD approximation problem (35) is not formulated in terms of factor matrices ($A, B, C, D$). A crucial reason is that in imaging applications, the interpretable factors are the spatial maps {$S_r$} and covariances matrices {$\Sigma_r$} and not the individual low-rank factorizations. Nonetheless, low-rank factors will play a crucial role to solve the optimization problem (35), as explained below in Section 4.3.

### 4.3. Proposed approach: alternating optimization

To efficiently exploit the structure of the optimization problem (35), we devise an optimization strategy inspired by the AO-ADMM framework [40]. The main idea consists of alternatively optimizing variables {$S_r$} and {$\Sigma_r$}, where sub-problems are solved through carefully designed ADMM algorithms, as detailed below. While the proposed approach will be referred to as AO-ADMM for clarity, it is important to note that such approach cannot be directly deduced from the one presented in [40], which requires separable constraints on the tensor factors. In contrast, the proposed algorithm imposes jointly the low-rank and non-negativity constraints on the spatial maps {$S_r$}, while for covariance matrices {$\Sigma_r$} imposing low-rank positive semidefiniteness requires a dedicated strategy.

At the outer level, the $(k+1)$-th iteration of the proposed algorithm reads

$$ \{S_r^{k+1}\}_{r=1}^R = \arg\min_{\substack{S_r \ge 0 \\ \mathrm{rank} S_r \le L_{1r}}} \left\| \mathcal{T} - \sum_{r=1}^R S_r \circ \Sigma_r^k \right\|_F^2, \quad (36) $$

$$ \{\Sigma_r^{k+1}\}_{r=1}^R = \arg\min_{\substack{\Sigma_r \ge 0 \\ \mathrm{rank} \Sigma_r \le L_{2r}}} \left\| \mathcal{T} - \sum_{r=1}^R S_r^{k+1} \circ \Sigma_r \right\|_F^2. \quad (37) $$

The sequence (36)–(37) defines a two-block-coordinate descent (BCD) scheme, where exact minimization is performed for each block. While the cost function is quadratic, the low-rank constraints make the sub-problems non-convex and non-smooth optimization problems. Therefore, standard results on convergence properties of BCD methods [41, 42], which require stronger assumptions such as convexity or smoothness, do not apply in this context. In this paper, we provide empirical evidence towards convergence properties of the proposed algorithm, on synthetic and real-world datasets. The theoretical analysis to support these findings is left for future work, where we anticipate that recent advances [43, 44] in non-convex BCD methods could be useful.

Algorithm 1 summarizes the full procedure for solving (35), where explicit updates for (36) and (37) are now given.

#### 4.3.1. Updates for spatial maps

For $r = 1, \dots, R$, we consider replacing the low-rank constraint on spatial maps by the explicit factorization $S_r = A_r B_r^\top$, where $A_r \in \mathbb{R}^{I \times L_{1r}}$ and $B_r \in \mathbb{R}^{I \times L_{2r}}$. We further introduce an augmented variable $\tilde{S}_r$ such
that $\tilde{S}_r = S_r$. Using this splitting, the optimization problem (36) becomes

$$ \min_{\{\tilde{S}_r = A_r B_r^\top\}} \left\| \mathcal{T} - \sum_{r=1}^{R} A_r B_r^\top \circ \Sigma_r^k \right\|_F^2 + \sum_{r=1}^{R} \iota_+(\tilde{S}_r) \quad (38) $$

where we have lifted the non-negativity constraint in the objective, thanks to the indicator function $\iota_+$ on the non-negative orthant. This equivalent formulation (38) can be efficiently solved by ADMM. Letting $A = [A_1 \cdots A_R]$, $B = [B_1 \cdots B_R]$, $U = [U_1 \cdots U_R]$ the dual variable associated to the constraints $\{\tilde{S}_r = A_r B_r^\top\}$ and $\gamma$ the augmented Lagrangian penalty parameter, the $(\ell+1)$-th ADMM iteration for (38) reads

$$ (A^{\ell+1}, B^{\ell+1}) = \arg \min_{A,B} \left\| \mathcal{T} - \sum_{r=1}^{R} A_r B_r^\top \circ \Sigma_r^k \right\|_F^2 + \gamma \sum_{r=1}^{R} \left\| \tilde{S}_r^\ell - A_r B_r^\top + U_r^\ell \right\|_F^2 \quad (39) $$

$$ \tilde{S}_r^{\ell+1} = \max(0, A_r^{\ell+1} B_r^{\ell+1\top} - U_r^{\ell}), \quad r = 1, \dots, R \quad (40) $$

$$ U_r^{\ell+1} = U_r^{\ell} + \tilde{S}_r^{\ell+1} - A_r^{\ell+1} B_r^{\ell+1\top}, \quad r = 1, \dots, R \quad (41) $$

where the last two lines can be computed in parallel. The solution to (39) can be computed by ALS; see explicit computations in Appendix A.1. Once a suitable convergence criterion is reached, we set $S_r^{k+1} := \tilde{S}_r^{\ell+1}$, $r = 1, \dots, R$ as the $(k+1)$-th estimate of spatial maps.

### 4.3.2. Updates for covariance matrices

Given $\{S_r^{k+1}\}$, one is now concerned with solving (37). We take advantage of the explicit parameterization of the constraints (low-rank and positive semidefiniteness) such that $\Sigma_r = C_r C_r^H \ge 0$, where $C_r \in \mathbb{C}^{K \times L_{2r}}$, $r = 1, \dots, R$. Then problem (37) is equivalent to the unconstrained optimization problem

$$ \min_{\{C_r\}_{r=1}^R} \left\| \mathcal{T} - \sum_{r=1}^R S_r^{k+1} \circ C_r C_r^H \right\|_F^2, \quad (42) $$

which is quartic in $\{C_r\}$ due to the symmetry in the last two modes of the BBTD model. To circumvent this issue, we use a splitting strategy to equivalently rewrite (42) as

$$ \min_{C=\bar{D}} \left\| \mathcal{T} - \sum_{r=1}^{R} S_r^{k+1} \circ C_r D_r^\top \right\|_F^2, \quad (43) $$

where we have used the shorthand notation $C = [C_1 \cdots C_R]$ and $D = [D_1 \cdots D_R]$. Observe that the problem (43) is now biconvex, i.e., it is convex in $C$ ($D$ being fixed) and convex in $D$ ($C$ being fixed). Therefore, we can devise an ADMM algorithm to solve (43). Let $V = [V_1 \cdots V_R]$ denote the dual variable associated to the constraint $\{C_r = \bar{D}_r\}$ and let $\rho$ be the Lagrangian penalty parameter. The $(\ell+1)$-th ADMM iterations for (43) read

$$ C^{\ell+1} = \arg \min_{C} \left\| \mathcal{T} - \sum_{r=1}^{R} S_r^{k+1} \circ C_r D_r^{\ell\top} \right\|_F^2 + \rho \left\| C - \bar{D}^{\ell} + V^{\ell} \right\|_F^2 \quad (44) $$

$$ D^{\ell+1} = \arg \min_{D} \left\| \mathcal{T} - \sum_{r=1}^{R} S_r^{k+1} \circ C_r^{\ell+1} D_r^\top \right\|_F^2 + \rho \left\| C^{\ell+1} - \bar{D} + V^{\ell} \right\|_F^2 \quad (45) $$

$$ V^{\ell+1} = V^{\ell} + C^{\ell+1} - \bar{D}^{\ell+1} \quad (46) $$
Explicit updates for Eqs. (44) and (45) are found in Appendix A.2. At convergence of inner loop iterations, we simply return the $(k+1)$-th update for covariances matrices as $\Sigma_r^{k+1} = \tilde{C}_r \tilde{C}_r^H$, $r = 1, \dots, R$ where $\tilde{C}$ denotes the symmetrization of the last updates of $C^\ell$ and $D^\ell$ to avoid propagation of numerical errors due to early stopping. See Algorithm 1 for further details.

**Algorithm 1** Exact AO-ADMM for solving (35)

**Input:** covariance matrix image $\mathcal{T}$, initial $\{S_r^0\}$, $\{\Sigma_r^0\}$.

**Parameters:** ADMM parameters $\gamma, \rho$

1: **while** exit criterion not satisfied **do**

2: Update of spatial maps – subproblem (36)

3: init $A^0, B^0$ from $\{S_r^k\}, \tilde{S}_r^0 = S_r^k, U^0 = 0$

4: **while** not converged **do**

5: compute $(A^{\ell+1}, B^{\ell+1})$ by ALS for (39)

6: compute $\{\tilde{S}_r^{\ell+1}\}$ by (40)

7: compute $U^{\ell+1}$ by (41)

8: $\ell \leftarrow \ell + 1$

9: **end while**

10: set $S_r^{k+1} = \tilde{S}_r^\ell, r = 1, \dots, R$

11: Update of covariance matrices – subproblem (37)

12: init $C^0, D^0$ from $\{\Sigma_r^k\}, V^0 = 0$

13: **while** not converged **do**

14: compute $C^{\ell+1}$ by solving normal equations (A.9)

15: compute $D^{\ell+1}$ by solving normal equations (A.10)

16: compute $V^{\ell+1}$ with (46)

17: $\ell \leftarrow \ell + 1$

18: **end while**

19: set $\tilde{C} = (C^\ell + \bar{D}^\ell)/2$

20: set $\Sigma_r^{k+1} = \tilde{C}_r \tilde{C}_r^H, r = 1, \dots, R$

21: $k \leftarrow k + 1$

22: **end while**

**Output:** estimated $\{\hat{S}_r\}, \{\hat{\Sigma}_r\}$

### 4.4. Speed-up with inexact ADMM iterations

A crucial bottleneck in Algorithm 1 lies in the update of spatial maps (sub-problem (36)), which consists in two nested loops due to: 1) the ADMM formulation of the subproblem and 2) the joint estimation of
factors ($A^{\ell+1}, B^{\ell+1}$) in (39) with ALS. The computational burden associated with these two operations is inevitably high and prevents the direct use of Algorithm 1 for large-scale applications. However, it is possible to circumvent this issue by relaxing the spatial maps ADMM iterations. More precisely, if one replaces the ADMM loop (lines 3-10 in Algorithm 1) and the inner ALS loop (line 5) by a single update of each variable, one gets the simple sequence

$$A^{k+1} = \arg \min_{A} \left\| \mathcal{T} - \sum_{r=1}^{R} A_r B_r^{k\top} \circ \Sigma_r^k \right\|_F^2 + \gamma \sum_{r=1}^{R} \| S_r^k - A_r B_r^{k\top} \|_F^2 \quad (47)$$

$$B^{k+1} = \arg \min_{B} \left\| \mathcal{T} - \sum_{r=1}^{R} A_r^{k+1} B_r^{\top} \circ \Sigma_r^k \right\|_F^2 + \gamma \sum_{r=1}^{R} \| S_r^k - A_r^{k+1} B_r^{\top} \|_F^2 \quad (48)$$

$$S_r^{k+1} = \max(0, A_r^{k+1} B_r^{k+1\top}), \quad r = 1, \dots, R \quad (49)$$

where computations in Eqs. (47) and (48) can be carried out explicitly following derivations in Appendix A.1. Comparing exact ADMM iterations (39) – (41) with “inexact” ones (47)–(49), remark that the dual variable $U$ has been eliminated (being initialized to zero in ADMM [45]). The intermediate variable $\tilde{S}$ has also disappeared, the update of spatial maps being done directly in (49). Overall, this relaxation greatly simplifies the original algorithm; see Algorithm 2 for details. As demonstrated in Section 5.1, it offers tremendous speed-up, while preserving empirical convergence performance.

### 4.5. Initialization strategies

Non-convex problems are known for their sensitivity to initial points; the problem (35) is no exception. Three different initialization strategies are considered hereafter:

* random initialization: For each $r = 1, \dots, R$, the entries of $S_r^0$ are drawn i.i.d. from the uniform distribution on $[0, 1]$, while $\Sigma_r^0 = C_r^0 C_r^{0\top}$, where $C_r^0$ has i.i.d. entries drawn from the standard normal distribution on $\mathbb{C}$.

* k-means based initialization: given data tensor $\mathcal{T} \in \mathbb{C}^{I \times J \times K \times K}$, one first vectorizes each pair of blocks to get a matrix $M \in \mathbb{C}^{I \times J \times K^2}$. Then k-means is run to determine $R$ vectorized covariance space clusters centers $\sigma_1, \dots, \sigma_R \in \mathbb{C}^{K^2}$, in the column space of $M$. Once determined, we obtain the $r$-th matrix $\Sigma_r^0$ by, in order: matricizing $\sigma_r$, projecting onto the SPD cone, and performing truncated SVD with rank $L_{2r}$. Corresponding spatial maps are initialized using k-means labels, reshaped as an $I \times J$ matrix, and binarized. Then, an NMF with rank $L_{1r}$ is performed to ensure that $S_r^0$ satisfies both non-negativity and low-rank constraints.

* projected ALS-BBTD: given data tensor $\mathcal{T}$, one first runs the ALS-BBTD algorithm (Eqs. (30)–(33)) to get $\mathcal{T} \approx [\hat{A}, \hat{B}, \hat{C}, \hat{D}]$. Then one forms, for each $r = 1, \dots, R$, the spatial maps $S_r^0 = \max(0, \hat{A}_r \hat{B}_r^{\top})$ and covariance matrices $\Sigma_r^0 = \tilde{C}_r \tilde{C}_r^{\top}$ where $\tilde{C}_r = (\hat{C}_r + \overline{\hat{D}_r})/2$.

Comparison of such strategies is presented in Section 5.2 below.
**Algorithm 2** Inexact AO-ADMM for solving (35)

**Input:** covariance matrix image $\mathcal{T}$, initial $\{S_r^0\}$, $\{\Sigma_r^0\}$.

**Parameters:** ADMM parameters $\gamma, \rho$

1: init $A^0, B^0$ from $\{S_r^0\}$

2: **while** exit criterion not satisfied **do**

3: (inexact) update of spatial maps – subproblem (36)

4: update $A^{k+1}$ by solving (47)

5: update $B^{k+1}$ by solving (48)

6: set $S_r^{k+1} = \max(0, A_r^{k+1}B_r^{k+1\top}), r = 1, \dots, R$

7: Update of covariance matrices – subproblem (37)

8: init $C^0, D^0$ from $\{\Sigma_r^k\}, V^0 = 0$

9: **while** not converged **do**

10: compute $C^{\ell+1}$ by solving normal equations (A.9)

11: compute $D^{\ell+1}$ by solving normal equations (A.10)

12: compute $V^{\ell+1}$ with (46)

13: $\ell \leftarrow \ell + 1$

14: **end while**

15: set $\tilde{C} = (C^\ell + \bar{D}^\ell)/2$

16: set $\Sigma_r^{k+1} = \tilde{C}_r \tilde{C}_r^H, r = 1, \dots, R$

17: $k \leftarrow k + 1$

18: **end while**

**Output:** estimated $\{\hat{S}_r\}$, $\{\hat{\Sigma}_r\}$
<figure><img src="image_4.png" /><figcaption>Figure 2: Difference between the cost function in Equation (35) and its optimal value for a synthetic tensor $T \in \mathbb{C}^{110 \times 120 \times 4 \times 4}$, with $R = 3$, $L_1 = 5$, and $L_2 = 2$ in a 45 dB noise scenario, comparing the performance of Algorithm 1 (represented in orange) and Algorithm 2 (represented in green). A dashed line indicates where the difference between the cost function and the optimal value reaches 90% of its final value. The markers on each curve correspond to the first and then every 5th outer loop iteration.</figcaption></figure>

# 5. Numerical experiments

## 5.1. Exact vs. Inexact algorithms

We start by comparing the exact AO-ADMM (Algorithm 1) with its non-exact counterpart (Algorithm 2). For this comparison, we construct a ground truth, complex-valued tensor $T_0$ of dimensions $110 \times 120 \times 4 \times 4$ following a BBTD model with $R = 3$, $L_{1r} = 5$ and $L_{2r} = 2$, for $r = 1, \dots, R$. Then noisy observations are constructed as $T = T_0 + \mathcal{E}$, where $\mathcal{E}$ has i.i.d. Gaussian entries where variance $\sigma^2$ is chosen such that the signal-to-noise ratio SNR = $\frac{\|T_0\|_F^2}{\lambda K^2 \sigma^2}$ is 45 dB. For fair comparison, both algorithms employ identical hyperparameters and the same exiting criteria: iterations were terminated when the relative change in the cost function fell below $10^{-6}$. ALS-BBTD initialization was used in both cases.

Figure 2 illustrates the difference between the cost function and its optimal value, in which we can see the impact of the two nested loops in the exact ADMM algorithm. The exact algorithm required over 170 seconds to converge (with an average of 1.7s/iteration), while the non-exact version took only 6.5 seconds (0.4s/iteration on average) on a M1 Max MacBook Pro using Python 3.12.4. Notably, the dashed line in the figure represents the point where the error reaches 90% of its final value. The inexact version reaches this threshold in approximately 2 seconds, whereas the exact version takes more than 54 seconds. Similar results (not shown) are obtained for other tensor sizes, noise level, initializations, and choices of BBTD model parameters. These results underscore the prohibitive computational cost of the exact Algorithm 1 in large-scale applications, whereas inexact AO-ADMM iterations preserve empirical convergence and exhibit great scalability. Thus, subsequent experiments will focus solely on Algorithm 2.
<figure><img src="image_5.png" /><figcaption>Figure 3: Performance analysis for initialization strategies in a scenario with 45 dB for Algorithm 2 evaluating the cost function represented in Equation (35). Blue lines correspond to 50 independent runs with random initialization. Orange dashed line with circles corresponds to *k*-means initialization, and the pink dashdotted line represents the performance using as initial guess the projected solution from ALS-BBTD as described in Section 4.5.</figcaption></figure>

## 5.2. Initialization strategies

We analyze now the behavior of Algorithm 2 for the different initializations proposed in Section 4.5. We use a similar setting as before, starting from a ground truth BBTD tensor $T_0 = [[A, B, C, \bar{C}]_{BBTD} \in \mathbb{C}^{110 \times 120 \times 4 \times 4}$, in which the factors matrices are drawn at random for $R = 3$, $L_{1r} = 5$ and $L_{2r} = 2$, $r = 1, \dots, R$. Considering a noisy measurement scenario with SNR = 45 dB, we evaluated the influence of three initialization strategies: random initialization, projected ALS-BBTD, and a *k*-means based approach.

Figure 3 illustrates the algorithm’s performance based on the cost function defined in Eq. (35) across the proposed initialization strategies. For random initialization, 50 independent Monte Carlo runs were performed, whereas the other two strategies produced consistent results, eliminating the need for multiple runs. By the end of the outer loop iterations, all initialization methods yielded comparable results. But some Monte Carlo runs presented a slower convergence, as shown in the curves in blue. Although *k*-means and ALS-BBTD yield similar convergence, *k*-means is significantly faster (0.1 s vs. 5 s) since it avoids the extra computation of the unconstrained solution. Therefore, we use *k*-means initialization in the following experiments.

## 5.3. Decomposition of real DTI data with BBTD

To demonstrate the potential of the BBTD model, we use the human brain DTI dataset [46], where each voxel contains a 3 × 3 diffusion matrix. A coronal slice is selected, forming the tensor $T \in \mathbb{R}^{127 \times 160 \times 3 \times 3}$ which is used for all the experiments of this section.

In Figure 4, we set all $L_{1r}$ to be equal to $L_1$ and then evaluate different $L_1$ values with fixed $R = 4$ and $L_2 = 1$. As expected, increasing $L_1$ introduces a tradeoff: higher computational complexity extends
<figure><img src="image_6.png" /><figcaption>Figure 4: Tradeoff between relative reconstruction error and execution time for different values of $L_1$ with fixed $R = 4$, and $L_2 = 1$ for a tensor $T \in \mathbb{R}^{127 \times 160 \times 3 \times 3}$.</figcaption></figure>

execution time, while simultaneously reducing the relative reconstruction error. However, beyond $L_1 = 50$, the improvement in reconstruction error becomes negligible, indicating that the spatial structure is inherently of that rank.

We analyze three BBTD decompositions of the same sample $T \in \mathbb{R}^{127 \times 160 \times 3 \times 3}$ with different decomposition parameters but fixed $R = 4$ (see Fig. 5). The left panel illustrates the original data. In (I), activation maps correspond to each covariance matrix element, e.g., $c_{xx}$ is $T(\cdot, \cdot, 0, 0)$. The data was acquired with a minimal six diffusion-weighted images, and approximately 1% of pixels contain non-PSD $3 \times 3$ matrices due to noise. In (II), an RGB representation based on eigenvectors highlights anatomical structures despite known ambiguities [38, 47, 48] in the mapping [38, Eq. 3.7]. In the standard RGB representation for DTI, the color is determined by the orientation of the principal eigenvector relative to the RGB axes, while an additional isotropic white component modulates the intensity according to the anisotropy level.

The right panel in Fig. 5 shows decompositions with increasing rank parameters. In (III) ($L_1 = 10$, $L_2 = 1$), spatial maps do not clearly delineate structures, but main diffusion directions remain well-defined in each vector $c_r$. In (IV) ($L_1 = 30$, $L_2 = 1$), increasing $L_1$ enhances anatomical clarity. In the first term ($S_1$) cerebellar commissures appear with strong activation, and the fornix is distinguishable between the lateral ventricles. White matter regions with similar diffusion properties show lower activation. The second and fourth terms represent the posterior limb of the internal capsule, respectively, while the third term captures the corpus callosum. Finally, (V) ($L_1 = 30$, $L_2 = 3$) provides the most refined segmentation, with all anatomical structures more distinctly visible and separated, and we can see that with rank-3 covariance matrices, the abundance maps are less prone to spatial overlapping. Each full-rank covariance matrix $\Sigma_r$ is
represented as an ellipsoid, where the eigenvectors define the orientation and the eigenvalues determine the scaling. When $L_2 = 1$, each pixel can be viewed as a weighted sum of $R = 4$ rank-1 which, when combined, produce a rank-3 matrix due to their dimensionality and structure. While each rank-1 term captures a single dominant direction, the resulting covariance reflects the contribution of multiple elementary components. It is also worth noting that, although the original noisy data included non-PSD covariance matrices, each estimated matrix in the BBTD decomposition is guaranteed to be PSD. Relative reconstruction errors are (III) 0.34, (IV) 0.22, and (V) 0.21, confirming that low-rank covariance matrices in BBTD effectively segment regions with similar diffusion properties while maintaining a low reconstruction error when compared to full-rank covariance matrices.

When comparing the BBTD-derived representation with the conventional RGB encoding shown in Fig. 5(II), we know that, while each color channel is nominally tied to one anatomical axis, the resulting hue at each pixel is modulated by both the principal diffusion direction and the anisotropy weighting. As a result, a pixel that appears green is not necessarily purely aligned along the corresponding axis; it may also contain contributions from other directions. Furthermore, the desaturation step (weighted by linear anisotropy) introduces an intensity-dependent bias: isotropic regions appear gray, and visual interpretation becomes less reliable. This is further compounded by the fact that human perception of color is nonlinear: small changes in one RGB component can dominate how the color is perceived. In contrast, our approach avoids these ambiguities by using activation maps directly tied to $rank-L_2$ diffusion profiles, ensuring a clearer and more interpretable representation of directional information without relying on arbitrary color encoding.

# 6. Conclusion and perspectives

This paper introduces the Block-Block Terms Decomposition (BBTD), a novel tool for decomposing 4th-order tensors in terms constructed from outer products of low-rank matrices. While the original motivation lies in matrix-valued imaging applications, this new model presents a standalone theoretical interest. To this aim, we studied its theoretical properties, including ambiguities, the number of free parameters, and an equivalent CONFAC model. Using third-order tensorizations, Theorem 1 provides sufficient uniqueness conditions for BBTD. These demonstrate the relevance of the 4th-order modeling to provide further relaxed uniqueness guarantees. To enable the use of BBTD as a meaningful low-complexity model for covariance matrix imaging applications, we propose a dedicated optimization framework that allows to cope with non-negativity and symmetry constraints – which are key to physical interpretation – as well as with low-rank constraints – which are key to guaranteeing the identifiability of the model. Eventually, numerical experiments on synthetic and real-data from DTI experiments illustrate the potential of the BBTD model.

This work triggers several future research directions. The first one will focus on evaluating the potential of
<figure><img src="image_7.png" /><figcaption>Figure 5: Three BBTD decompositions of a coronal slice from a human brain DTI dataset [46], where each pixel contains a 3 × 3 covariance matrix representing diffusion properties along three axes. (I) Activation maps for each covariance matrix element. (II) RGB representation derived from principal diffusion eigenvectors, with key anatomical structures labeled. (III) Decomposition with $L_1 = 10$, $L_2 = 1$, and $R = 4$, where each spatial map $S_r$ is associated with a direction vector $c_r$. (IV) Increasing $L_1$ to 30 while keeping $L_2 = 1$ illustrates the effect of higher spatial rank. (V) Setting $L_2 = 3$ associates each component $S_r$ with a full-rank covariance matrix, visualized as ellipsoids defined by eigenvalues (scaling) and eigenvectors (orientation).</figcaption></figure>

the BBTD model in other relevant matrix-valued imaging scenarios, such as the recently proposed excitation-emission matrix fluorescence imaging [5]. This will require different constraints on the BBTD factors and, therefore, the development of dedicated algorithms. The second perspective will target the extension of the BBTD framework to arbitrary order-N tensors: indeed, one can define a general (p,q)-BBTD model such that $p+q=N$ and where block-block terms are formed by outer products of low-rank tensors of order $p$ and $q$, respectively. Such decomposition requires further analysis and characterization, yet it opens exciting research avenues to break the curse of dimensionality of higher-order tensors.

## Acknowledgements

The authors would like to thank Dr. Ricardo Borsoi for stimulating discussions and feedback during the writing of this manuscript.
## Appendix A. Computational details for Algorithms 1 and 2

### Appendix A.1. Explicit ALS updates for (39)

The main challenge in deriving ALS-type updates for (39) lies in handling properly the real-valuedness of matrices $A$ and $B$ when the remaining factors are complex-valued. We first consider the estimation of $A$ given fixed $B$. We drop iteration indices $k, \ell$ for convenience. Remark that the BBTD fitting error can be expressed using mode-1 unfolding of the 3-rd order tensorization $T^{3:4}$ of $T$ such that

$$ \left\| \mathcal{T} - \sum_{r=1}^{R} A_r B_r^\top \circ \Sigma_r \right\|_F^2 = \left\| T_{(1)}^{3:4} - A (\Sigma \Theta_1 \odot B)^\top \right\|_F^2 \quad (A.1) $$

with $\Sigma = [\text{vec } \Sigma_1 \dots \text{vec } \Sigma_R]$ and $\Theta_1 = \text{diag}(\mathbf{1}_{L_{11}}, \dots, \mathbf{1}_{L_{1R}})$ the constraint matrix associated to the decomposition into rank-(1, $L_{1r}$, $L_{1r}$) terms. Let $M_1 = (\Sigma\Theta_1 \odot B)^\top$. Then, expand (A.1) as

$$ \begin{aligned} \|T_{(1)}^{3:4} - AM_1\|_F^2 &= \operatorname{Tr}(T_{(1)}^{3:4} - AM_1)^H (T_{(1)}^{3:4} - AM_1) \\ &= \operatorname{Tr} T_{(1)}^{3:4H} T_{(1)}^{3:4} - \operatorname{Tr} M_1^H A^\top T_{(1)}^{3:4} \\ &\quad - \operatorname{Tr} T_{(1)}^{3:4H} AM_1 + \operatorname{Tr} M_1^H A^\top AM_1 \end{aligned} \quad (A.2) $$

where we have used the fact that $A$ is real. Similarly, the second term in (39) is expanded as

$$ \sum_{r}^{R} \| \tilde{S}_r - A_r B_r^\top + U_r \|_F^2 = \sum_{r=1}^{R} \operatorname{Tr} \tilde{S}_r^\top \tilde{S}_r - 2 \operatorname{Tr} A_r B_r^\top \tilde{S}_r^\top + \operatorname{Tr} B_r A_r^\top A_r B_r^\top. \quad (A.3) $$

where we used shorthand notation $\tilde{S}_r = \tilde{S}_r + U_r$. Derivating (A.2) with respect to $A$ one gets

$$ \begin{aligned} \nabla_A \|T_{(1)}^{3:4} - AM_1\|_F^2 &= -T_{(1)}^{3:4} M_1^H - \overline{T_{(1)}^{3:4}} M_1^\top + A [M_1 M_1^H + (M_1 M_1^H)^\top] \\ &= -2 \operatorname{Re}\{T_{(1)}^{3:4} M_1^H\} + 2 A \operatorname{Re}\{M_1 M_1^H\}. \end{aligned} \quad (A.4) $$

To compute the gradient of (A.3), we proceed blockwise and obtain

$$ \frac{1}{2} \nabla_A \sum_{r}^{R} \| A_r B_r^\top - \tilde{S}_r \|_F^2 = A_r \text{ bdiag}(B_1^\top B_1, \dots, B_R^\top B_R) - [\tilde{S}_1 B_1 \dots \tilde{S}_R B_R]. \quad (A.5) $$

Eventually, from Eqs. (A.4) and (A.5) one gets the normal equations for the $A$-update

$$ A (\operatorname{Re}\{M_1 M_1^H\} + \gamma \operatorname{diag}(B_1^\top B_1, \dots, B_R^\top B_R)) = \operatorname{Re}\{T_{(1)}^{3:4} M_1^H\} + \gamma [\tilde{S}_1 B_1 \dots \tilde{S}_R B_R]. \quad (A.6) $$

The normal equations for $B$ can be obtained in a similar way:

$$ B (\operatorname{Re}\{M_2 M_2^H\} + \gamma \operatorname{diag}(A_1^\top A_1, \dots, A_R^\top A_R)) = \operatorname{Re}\{T_{(2)}^{3:4} M_2^H\} + \gamma [\tilde{S}_1^\top A_1 \dots \tilde{S}_R^\top A_R] \quad (A.7) $$

where $M_2 = (\Sigma\Theta_1 \odot A)^\top$. Solving alternatively the normal equations (A.6) – (A.7) yields the ALS-update for problem (39).
Appendix A.2. Explicit updates for Eqs. (44) and (45)

Start with the computation of $C^{\ell+1}$ in (44). We drop the outer iteration index $k$ for notational convenience. First remark that the BBTD fitting error can be expressed using mode-2 unfolding of the 3-rd order tensorization $\mathcal{T}^{1:2}$ of $\mathcal{T}$ such that

$$ \left\| \mathcal{T} - \sum_{r=1}^{R} S_r \circ C_r D_r^{\ell \top} \right\|_F^2 = \left\| T_{(2)}^{1:2} - C(D^{\ell} \odot S\Theta_2)^{\top} \right\|_F^2 \quad (A.8) $$

with $S = [\text{vec } S_1 \dots \text{vec } S_R]$ and where $\Theta_2 = \text{bdiag}(\mathbf{1}_{L_{21}^\top}, \dots, \mathbf{1}_{L_{2R}^\top})$ is the constraint matrix associated to the decomposition into rank-(1, $L_{2r}$, $L_{2r}$) terms. Since $C$ is complex-valued, computing the gradient of the cost function in (44) requires the computation of Wirtinger derivatives [49] with respect to $\bar{C}$. Setting the gradient to zero, one obtains the normal equations

$$ C(M_3 M_3^H + \rho I) = T_{(2)}^{1:2} M_3^H + \rho (\bar{D}^\ell - V^\ell), \quad (A.9) $$

where $M_3 = (D^\ell \odot S\Theta_2)^\top$. Similar computations for (45) yield the normal equations for the update of $D$

$$ D(M_4 M_4^H + \rho I) = T_{(3)}^{1:2} M_4^H + \rho (\bar{V}^\ell + \bar{C}^{\ell+1}), \quad (A.10) $$

where $M_4 = (C^{\ell+1} \odot S\Theta_2)^\top$.

References

[1] R. M. Azzam, Stokes-vector and mueller-matrix polarimetry, Journal of the Optical Society of America A 33 (7) (2016) 1396–1408.

[2] R. Chipman, W. S. T. Lam, G. Young, Polarized light and optical systems, CRC press, 2018.

[3] J.-S. Lee, E. Pottier, Polarimetric radar imaging: from basics to applications, CRC press, 2017.

[4] J.-D. Tournier, S. Mori, A. Leemans, Diffusion tensor imaging and beyond, Magnetic resonance in medicine 65 (6) (2011) 1532.

[5] O. Katz, T. Ferguson, E. Abbey, S.-J. Klose, C. Prüfert, H.-P. Loock, Fluorescence excitation-emission-matrix imaging, Analytical Chemistry 95 (34) (2023) 12631–12639.

[6] A. Collas, F. Bouchard, A. Breloy, G. Ginolhac, C. Ren, J.-P. Ovarlez, Probabilistic pca from heteroscedastic signals: geometric framework and application to clustering, IEEE Transactions on Signal Processing 69 (2021) 6546–6560.

[7] A. Goh, R. Vidal, Clustering and dimensionality reduction on riemannian manifolds, in: 2008 IEEE Conference on computer vision and pattern recognition, IEEE, 2008, pp. 1–7.

[8] P. T. Fletcher, S. Joshi, Principal geodesic analysis on symmetric spaces: Statistics of diffusion tensors, in: International Workshop on Mathematical Methods in Medical and Biomedical Image Analysis, Springer, 2004, pp. 87–98.

[9] S. Said, L. Bombrun, Y. Berthoumieu, J. H. Manton, Riemannian gaussian distributions on the space of symmetric positive definite matrices, IEEE Transactions on Information Theory 63 (4) (2017) 2153–2170.

[10] V. Arsigny, P. Fillard, X. Pennec, N. Ayache, Log-euclidean metrics for fast and simple calculus on diffusion tensors, Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine 56 (2) (2006) 411–421.

[11] P. Fillard, X. Pennec, V. Arsigny, N. Ayache, Clinical dt-mri estimation, smoothing, and fiber tracking with log-euclidean metrics, IEEE transactions on medical imaging 26 (11) (2007) 1472–1482.
[12] C. Chefd’Hotel, D. Tschumperlé, R. Deriche, O. Faugeras, Regularizing flows for constrained matrix-valued images, Journal of Mathematical Imaging and Vision 20 (2004) 147–162.

[13] X. Pennec, Manifold-valued image processing with spd matrices, in: Riemannian geometric statistics in medical image analysis, Elsevier, 2020, pp. 75–134.

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

[15] F. L. Hitchcock, The expression of a tensor or a polyadic as a sum of products, Journal of Mathematics and Physics 6 (1-4) (1927) 164–189.

[16] C. M. Andersen, R. Bro, Practical aspects of parafac modeling of fluorescence excitation-emission data, Journal of Chemometrics: A Journal of the Chemometrics Society 17 (4) (2003) 200–215.

[17] L. De Lathauwer, Decompositions of a higher-order tensor in block terms—Part II: Definitions and uniqueness, SIAM Journal on Matrix Analysis and Applications 30 (3) (2008) 1033–1066.

[18] C. Prévost, R. A. Borsoi, K. Usevich, D. Brie, J. C. Bermudez, C. Richard, Hyperspectral super-resolution accounting for spectral variability: Coupled tensor LL1-based recovery and blind unmixing of the unknown super-resolution image, SIAM Journal on Imaging Sciences 15 (1) (2022) 110–138.

[19] M. Ding, X. Fu, T.-Z. Huang, J. Wang, X.-L. Zhao, Hyperspectral super-resolution via interpretable block-term tensor modeling, IEEE Journal of Selected Topics in Signal Processing 15 (3) (2020) 641–656.

[20] M. Ding, X. Fu, X.-L. Zhao, Fast and structured block-term tensor decomposition for hyperspectral unmixing, IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 16 (2023) 1691–1709.

[21] J. J. Gil, Polarimetric characterization of light and media: Physical quantities involved in polarimetric phenomena, The European Physical Journal Applied Physics 40 (1) (2007) 1–47. doi:10.1051/epjap:2007153.

[22] T. G. Kolda, B. W. Bader, Tensor Decompositions and Applications, SIAM Review 51 (3) (2009) 455–500. doi:10.1137/07070111X.

[23] J. B. Kruskal, Three-way arrays: rank and uniqueness of trilinear decompositions, with application to arithmetic complexity and statistics, Linear algebra and its applications 18 (2) (1977) 95–138.

[24] N. D. Sidiropoulos, R. Bro, On the uniqueness of multilinear decomposition of N-way arrays, Journal of Chemometrics: A Journal of the Chemometrics Society 14 (3) (2000) 229–239.

[25] L. De Lathauwer, Decompositions of a higher-order tensor in block terms—Part I: Lemmas for partitioned matrices, SIAM Journal on Matrix Analysis and Applications 30 (3) (2008) 1022–1032.

[26] G. Favier, A. L. de Almeida, Overview of constrained PARAFAC models, EURASIP Journal on Advances in Signal Processing 2014 (2014) 1–25.

[27] R. Bro, R. A. Harshman, N. D. Sidiropoulos, M. E. Lundy, Modeling multi-way data with linearly dependent loadings, Journal of Chemometrics: A Journal of the Chemometrics Society 23 (7-8) (2009) 324–340.

[28] A. L. de Almeida, G. Favier, J. C. M. Mota, A constrained factor decomposition with application to MIMO antenna systems, IEEE Transactions on Signal Processing 56 (6) (2008) 2429–2442.

[29] A. Stegeman, A. L. De Almeida, Uniqueness conditions for constrained three-way factor decompositions with linearly dependent loadings, SIAM Journal on Matrix Analysis and Applications 31 (3) (2010) 1469–1490.

[30] D. Brie, S. Miron, F. Caland, C. Mustin, An uniqueness condition for the 4-way candecomp/parafac model with collinear loadings in three modes, in: 2011 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), IEEE, 2011, pp. 4108–4111.

[31] X. Guo, S. Miron, D. Brie, A. Stegeman, Uni-mode and partial uniqueness conditions for candecomp/parafac of three-way arrays with linearly dependent loadings, SIAM Journal on Matrix Analysis and Applications 33 (1) (2012) 111–129.

[32] A. L. de Almeida, X. Luciani, A. Stegeman, P. Comon, CONFAC decomposition approach to blind identification of
underdetermined mixtures based on generating function derivatives, IEEE Transactions on Signal Processing 60 (11) (2012) 5698-5713.

[33] P. Comon, X. Luciani, A. L. De Almeida, Tensor decompositions, alternating least squares and other tales, Journal of Chemometrics: A Journal of the Chemometrics Society 23 (7-8) (2009) 393-405.

[34] L. De Lathauwer, D. Nion, Decompositions of a higher-order tensor in block terms—part iii: Alternating least squares algorithms, SIAM journal on Matrix Analysis and Applications 30 (3) (2008) 1067-1083.

[35] D. H. Goldstein, Polarized light, CRC press, 2017.

[36] R. Ossikovski, O. Arteaga, Complete Mueller matrix from a partial polarimetry experiment: the nine-element case, Journal of the Optical Society of America A 36 (3) (2019) 403-415.

[37] J. J. Gil, Characteristic properties of Mueller matrices, Journal of the Optical Society of America A 17 (2) (2000) 328-334.

[38] G. Kindlmann, Visualization and analysis of diffusion tensor fields, Ph.D. thesis, The University of Utah, Salt Lake City, Utah (Sep. 2004).

[39] M. Lara, A. Moll, A. Mas, M. J. Picado, C. Gassent, J. Pomar, J. A. Llompart-Pou, M. Brell, J. Ibáñez, J. Pérez-Bárcena, Use of diffusion tensor imaging to assess the vasogenic edema in traumatic pericontusional tissue, Neurocirugia 32 (4) (2021) 161-169.

[40] K. Huang, N. D. Sidiropoulos, A. P. Liavas, A flexible and efficient algorithmic framework for constrained matrix and tensor factorization, IEEE Transactions on Signal Processing 64 (19) (2016) 5052-5065.

[41] L. Grippo, M. Sciandrone, Globally convergent block-coordinate techniques for unconstrained optimization, Optimization Methods and Software 10 (4) (1999) 587-637. doi:10.1080/10556789908805730.

[42] D. P. Bertsekas, Nonlinear programming, 2nd Edition, Athena Scientific, 1999.

[43] Y. Xu, W. Yin, A Globally Convergent Algorithm for Nonconvex Optimization Based on Block Coordinate Update, Journal of Scientific Computing 72 (2) (2017) 700-734. doi:10.1007/s10915-017-0376-0.

[44] L. Peng, R. Vidal, Block Coordinate Descent on Smooth Manifolds: Convergence Theory and Twenty-One Examples, arXiv:2305.14744 (Oct. 2023). doi:10.48550/arXiv.2305.14744.

[45] S. Boyd, N. Parikh, E. Chu, B. Peleato, J. Eckstein, et al., Distributed optimization and statistical learning via the alternating direction method of multipliers, Foundations and Trends® in Machine learning 3 (1) (2011) 1-122.

[46] G. Kindlmann, A. Alexander, Brain dataset, Courtesy of Gordon Kindlmann at the Scientific Computing and Imaging Institute, University of Utah, and Andrew Alexander, W. M. Keck Laboratory for Functional Brain Imaging and Behavior, University of Wisconsin-Madison (February 2000).
URL https://www.sci.utah.edu/~gk/DTI-data/

[47] S. Pajevic, C. Pierpaoli, Color schemes to represent the orientation of anisotropic tissues from diffusion tensor data: application to white matter fiber tract mapping in the human brain, Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine 42 (3) (1999) 526-540.

[48] A. L. Alexander, J. E. Lee, M. Lazar, A. S. Field, Diffusion tensor imaging of the brain, Neurotherapeutics 4 (3) (2007) 316-329.

[49] A. Hjorungnes, D. Gesbert, Complex-valued matrix differentiation: Techniques and key results, IEEE Transactions on Signal Processing 55 (6) (2007) 2740-2746.