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

HAL
open science

# Physically-constrained block-term tensor decomposition for polarimetric image recovery

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

► To cite this version:

Saulo Cardoso Barreto, Julien Flamant, Sebastian Miron, David Brie. Physically-constrained block-term tensor decomposition for polarimetric image recovery. International Conference on Acoustics, Speech and Signal Processing, ICASSP 2024, Apr 2024, Seoul, South Korea. hal-04344657v2

HAL Id: hal-04344657

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

Submitted on 12 Jan 2024

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
# PHYSICALLY-CONSTRAINED BLOCK-TERM TENSOR DECOMPOSITION FOR POLARIMETRIC IMAGE RECOVERY

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

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

## ABSTRACT

This paper introduces a complete approach for the recovery of polarimetric images from experimental intensity measurements. In many applications, such images collect, at each pixel, a Stokes vector encoding the polarization state of light. By representing a Stokes vector image as a third-order tensor, we propose a new physically-constrained block-term tensor decomposition called Stokes-BTD. The proposed model is flexible and comes with broad identifiability guarantees. Moreover, physical constraints ensure meaningful interpretation of low-rank terms as Stokes vectors. In practice, Stokes images must be recovered from indirect, intensity measurements. To this aim, we implement two recovery algorithms for Stokes-BTD based on constrained alternated optimization and highlight constraints related to Stokes vectors. Numerical experiments on synthetic and real data illustrate the potential of the approach.

**Index Terms**— block-term tensor decomposition, Stokes polarimetric imaging, alternated constrained optimization.

## 1. INTRODUCTION

Polarimetric imaging is crucial for many applications including biological tissue assessment [1], remote sensing [2], and material characterization [3]. It exploits light polarization, a fundamental vector property describing the trajectory of the electromagnetic (EM) field. The interaction of polarized light with polarization-sensitive media makes polarization properties change, therefore revealing key physical features (e.g., anisotropy) that are inaccessible to conventional imaging techniques. As a result, polarimetric imaging is often used to build highly contrasted images, which can be further processed depending on the task at hand.

Polarization information can be represented in many ways [4]. The choice of one representation over the others usually depends on the wavelength (e.g., radar or optical frequencies). This work focuses on the optical imaging context: in this setting, the polarization properties of the EM field are, at each pixel, encoded by a four-dimensional real vector called Stokes vector. However, such Stokes polarimetric images cannot be measured directly. Instead in practice, several polarimetric projections (i.e., intensity measurements through different polarizers) are performed, and the Stokes polarimetric image has to be recovered in a second step. Several methods for the recovery of Stokes images from intensity measurements have been proposed. These include Bayesian approaches [5, 6], non-local means [7], regularized [8] and min-max optimization [9]. In contrast, this paper exploits the natural representation of a Stokes image as a third-order tensor to leverage the powerful framework of multilinear algebra, in particular, low-rank tensor decompositions [10].

<figure><img src="image_3.png" /><figcaption>Fig. 1: Stokes vectors of linear, circular, and elliptical polarization states.</figcaption></figure>

This work proposes a new constrained block-term tensor decomposition (BTD) for the recovery of Stokes polarimetric images from noisy intensity measurements. The approach allows simultaneously: i) the reconstruction of a physically meaningful Stokes image and ii) the segmentation of the Stokes image into distinct areas with identical polarization properties. The proposed model, called Stokes-BTD, incorporates physical constraints related to the interpretation of Stokes parameters. It is flexible and benefits from broad identifiability guarantees. Two constrained optimization strategies for fitting Stokes-BTD are proposed and carefully evaluated on synthetic and real-world polarized microscopy data.

## 2. STOKES POLARIMETRIC IMAGING

### 2.1. Stokes parameters

In optical polarimetric imaging, the polarization state of light is described by a Stokes vector $s = [s_0, s_1, s_2, s_3]^\top \in \mathbb{R}^4$ where $s_0, s_1, s_2$ and $s_3$ are the Stokes parameters [4]. These four real values are interpreted as follows: the first parameter $s_0$ indicates the total intensity, i.e., the sum of intensities of the polarized and unpolared part of light. The three remaining parameters $s_1, s_2$, and $s_3$ describe the polarized part. The value of $s_1$ represents the amount of horizontal or vertical polarization, while $s_2$ gives the polarized component along the 45°/135° axis. Lastly, $s_3$ indicates the preponderance between left and right circular polarizations. Importantly, there is a one-to-one correspondence between a Stokes vector and an oriented ellipse, called *polarization ellipse* in polarimetric imaging. Fig. 1 shows this correspondence for some elementary polarization states. For $s_3 = 0$ this ellipse degenerates into a line segment (linear polarization) while for $s_1 = s_2 = 0$ it becomes a circle. For arbitrary $s_1, s_2, s_3$ polarization is in general elliptic. Stokes vectors must respect strong *physical constraints*. In particular, a Stokes vector $s = [s_0, s_1, s_2, s_3]^\top$ must belong to $\mathbb{R}_S^4$ where $S$ denotes the set of constraints given by

$$s_0 \ge 0 \quad \text{and} \quad s_0^2 \ge s_1^2 + s_2^2 + s_3^2. \qquad (S)$$

The first constraint $s_0 \ge 0$ is classical and simply refers to the interpretation of $s_0$ as an intensity. The second constraint involves all

<footnote>Research supported by project ANR RICOCHET N° ANR-21-CE48-0013.</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 2: Direct measurement model of Stokes polarimetric data (2).</figcaption></figure>

four components of s and defines the set of physical polarization parameters $s_1, s_2$, and $s_3$. From a mathematical viewpoint, the subset $\mathbb{R}_S^4 \subset \mathbb{R}^4$ defines a second-order cone which can be interpreted as an extension of the classical non-negativity constraint to the vector case. To see this, one must exploit the bijective mapping between a Stokes vector s and the so-called coherence matrix $\mathbf{J}(s)$ given by [4]

$$ \mathbf{J}(s) = \frac{1}{2} \begin{bmatrix} s_0 + s_1 & s_2 + i s_3 \\ s_2 - i s_3 & s_0 - s_1 \end{bmatrix} \in \mathbb{C}^{2 \times 2}. \quad (1) $$

By physical construction [4], the matrix $\mathbf{J}$ must be a covariance matrix and hence positive semi-definite (or simply, non-negative definite). Moreover, straightforward calculations show that $s \in \mathbb{R}_S^4 \Leftrightarrow \mathbf{J}(s) \succeq 0$. This can be easily shown using the fact that a 2-by-2 complex matrix is positive semi-definite if and only if its trace and determinant are both non-negative. We will exploit this equivalence in the sequel for the practical implementation of (S).

## 2.2. Stokes polarimetric image measurement

A Stokes polarimetric image collects, for each pixel, the Stokes vector encoding the polarization state of light. Formally, for a grid of $I \times J$ pixels, a Stokes polarimetric image defines a three-dimensional array or third-order tensor $\mathcal{X} \in \mathbb{R}_S^{I \times J \times 4}$ such that the pixel $(i, j)$ reads $\mathcal{X}_{ij} = s(i, j) \in \mathbb{R}_S^4$. However, in optical imaging, only intensity measurements can be performed, making it impossible to retrieve the Stokes vector image $\mathcal{X}$ directly. Instead, the measurement of Stokes parameters relies on performing a sequence of $M$ polarimetric projections [11], leading to $M$ intensity images $\mathbf{Y}_1, \mathbf{Y}_2, \dots, \mathbf{Y}_M \in \mathbb{R}^{I \times J}$. This is easily implemented in practice by using optical components called polarizers. Formally, the $m$-th intensity image $\mathbf{Y}_m$ can be related to the Stokes tensor $\mathcal{X}$ by the pixel-wise operation $\mathbf{Y}_m(i, j) = \mathbf{p}_m^\top \mathcal{X}_{ij}$, where $\mathbf{p}_m \in \mathbb{R}_S^4$ is a reference Stokes vector supposed known. Physically, the image $\mathbf{Y}_m$ encodes the intensity of light polarized along $\mathbf{p}_m$ in $\mathcal{X}$. Finally, note that in experiments an overcomplete set of measurements is always considered ($M \ge 4$), such that the instrument matrix $\mathbf{P} \in \mathbb{R}^{M \times 4}$ gathering reference polarization states has full column rank.

The Stokes polarimetric image measurement model can be written efficiently using tensor algebra. Let $\mathcal{Y} \in \mathbb{R}^{I \times J \times M}$ be the tensor gathering the $M$ intensity measurements in its frontal slices $\mathbf{Y}_1, \mathbf{Y}_2, \dots, \mathbf{Y}_M \in \mathbb{R}^{I \times J}$. The measurement model then reads

$$ \mathcal{Y} = \mathcal{X} \bullet_3 \mathbf{P} + \mathcal{N}, \quad (2) $$

where $\bullet_3$ denotes the tensor mode product along the third dimension and $\mathcal{N}$ encodes measurement noise with i.i.d. Gaussian entries

$\mathcal{N}_{ijm} \sim \mathcal{N}(0, \sigma^2)$ with noise variance $\sigma^2$ – a reasonable assumption for many applications. See Fig. 2 for illustration.

It is important to note that a common practice in optics is to directly invert the model (2) by using the unconstrained least square solution $\hat{\mathcal{X}}_{\text{ULS}} = \mathcal{Y} \bullet_3 \mathbf{P}^\dagger$, where $\mathbf{P}^\dagger$ is the pseudo-inverse of $\mathbf{P}$. However, this approach does not guarantee at all that estimated Stokes vectors satisfy the physical constraints (S). Moreover, the analysis and decomposition of the Stokes polarimetric image into elementary regions has to be performed in a second step, e.g., using clustering approaches such as k-means or spatial mixture models [6]. This paper proposes to address these two limitations by taking advantage of the tensor nature of Stokes polarimetric images.

# 3. CONSTRAINED BLOCK-TERM DECOMPOSITION

## 3.1. Stokes-BTD model

In many practical settings, it is reasonable to suppose a linear mixture model for $\mathcal{X}$. Under this assumption, the Stokes vector $\mathcal{X}_{ij} = s(i, j)$ at a pixel $(i, j)$ can be written as the linear combination of $R$ elementary Stokes vectors $s_1, \dots, s_R \in \mathbb{R}_S^4$ such that $\mathcal{X}_{ij} = \sum_{r=1}^R (\mathbf{X}_r)_{ij} \cdot \mathbf{s}_r$, where $\mathbf{X}_1, \dots, \mathbf{X}_R \in \mathbb{R}_+^{I \times J}$ are the corresponding $R$ spatial mixing matrices. This model corresponds to the incoherent superposition of Stokes vectors associated with the different regions of the polarimetric image [4]. It can also be viewed as the polarimetric counterpart of the standard linear mixing model in hyperspectral imaging [12]. The model can be written in a compact way as

$$ \mathcal{X} \approx \sum_{r=1}^{R} \mathbf{X}_r \circ \mathbf{s}_r, \quad (3) $$

where the symbol $\circ$ denotes the tensor dot product. The matrices $\mathbf{X}_r$ are assumed to be non-negative, which ensures their interpretation as activation maps of the corresponding Stokes parameters vector $\mathbf{s}_r$. This also guarantees that (3) represents a combination of elements in $\mathbb{R}_S^4$ that lie within a cone, assuring that each pixel in the reconstructed polarimetric image is a Stokes vector, i.e. it satisfies the constraints specified in (S).

To reduce the number of parameters appearing in (3), a common approach is to model each spatial mixture matrix $\mathbf{X}_r$ by a (low-) rank $L$ approximation where $L \ll I, J$. Formally, one has $\mathbf{X}_r \approx \mathbf{A}_r \mathbf{B}_r^\top$, where $\mathbf{A}_r \in \mathbb{R}_+^{I \times L}$ and $\mathbf{B}_r \in \mathbb{R}_+^{J \times L}$ inherit non-negativity constraints from $\mathbf{X}_r$. Finally, the proposed Stokes-BTD model for a Stokes tensor $\mathcal{X}$ reads

$$ \mathcal{X} \approx \sum_{r=1}^{R} (\mathbf{A}_r \mathbf{B}_r^\top) \circ \mathbf{s}_r := \llbracket \mathbf{A}, \mathbf{B}, \mathbf{S} \rrbracket_{\text{SBTD}} \quad (4) $$

where the matrices $\mathbf{A} = [\mathbf{A}_1 \cdots \mathbf{A}_R] \in \mathbb{R}_+^{I \times LR}$ and $\mathbf{B} = [\mathbf{B}_1 \cdots \mathbf{B}_R] \in \mathbb{R}_+^{J \times LR}$ collect the $R$ blocks of size $I \times L$ and $J \times L$, respectively, and where $\mathbf{S} = [\mathbf{s}_1 \cdots \mathbf{s}_R] \in \mathbb{R}_S^{4 \times R}$ gathers the $R$ elementary Stokes vectors. The Stokes-BTD model defines a constrained block-term tensor decomposition (BTD) (also known as a rank-(L, L, 1) decomposition [13]). These constraints encode the physics of Stokes polarimetric images and ensure the meaningful interpretation of the low-rank terms of the decomposition (4) [14].

One of the key features of the Stokes-BTD model is that it is identifiable under broad conditions. Recall that a tensor decomposition is identifiable if the only ambiguities are the trivial ambiguities of the model. In particular, for (4) these correspond to permutation of the $R$ different terms or a scale factor $\lambda_r \neq 0$ such that
($\lambda_r A_r B_r^\top$) ◦ ($\lambda_r^{-1} s_r$) = $A_r B_r^\top$ ◦ $s_r$. Moreover, only the sub-space generated by $A_r$ and $B_r$ is identifiable since it is always possible to find an invertible matrix $D$ so that $A_r B_r^\top = A_r D D^{-1} B_r^\top$. The interpretation of (4) as a specific BTD model makes it possible to leverage existing sufficient uniqueness guarantees (see, e.g., [13, Sec. 4]) to study the identifiability of the Stokes-BTD model. The simplest condition states that if $A$ and $B$ are full column-rank and $S$ does not have collinear columns, then the model (4) is identifiable. Finally, note that in practical experiments where the Stokes polarimetric image $\mathcal{X}$ has to be retrieved from measurements $\mathcal{Y} \approx \mathcal{X}_{\bullet 3} P$, identifiability of the Stokes-BTD model is preserved by the forward model (2) since the instrument matrix $P \in \mathbb{R}^{M \times 4}$ will always be designed to have full column rank.

## 3.2. Algorithms for Stokes-BTD

We now consider algorithms for reconstructing a Stokes-BTD polarimetric image $\hat{\mathcal{X}} = [\hat{\mathbf{A}}, \hat{\mathbf{B}}, \hat{\mathbf{S}}]_{\text{SBTD}} \in \mathbb{R}^{I \times J \times 4}$ given a (noisy) tensor $\mathcal{Y} \in \mathbb{R}^{I \times J \times M}$ of $M \ge 4$ intensity images given by the measurement model (2). Since noise is assumed i.i.d. Gaussian, the inverse problem admits the constrained least squares formulation

$$ (\hat{\mathbf{A}}, \hat{\mathbf{B}}, \hat{\mathbf{S}}) = \underset{\mathbf{A}, \mathbf{B} \ge 0, \mathbf{S} \in \mathbb{R}_S^{4 \times R}}{\operatorname{argmin}} \|\mathcal{Y} - [\mathbf{A}, \mathbf{B}, \mathbf{S}]_{\text{SBTD}} \bullet_3 \mathbf{P}\|^2. \quad (5) $$

The problem (5) is convex in each matrix variable separately, since the cost is quadratic and associated constraints are convex (non-negativity and polarization constraint). However it is not jointly convex in $(\mathbf{A}, \mathbf{B}, \mathbf{S})$. Therefore, we adopt a classical constrained alternating minimization approach which consists in solving (5) successively for $\mathbf{A}, \mathbf{B}$, and $\mathbf{S}$. We rely on the fact that the cost function in (5) can be conveniently expressed using tensor unfoldings and factor matrices of the BTD model. Let us denote by $\mathbf{Y}_{(1)} \in \mathbb{R}^{J M \times I}$, $\mathbf{Y}_{(2)} \in \mathbb{R}^{I M \times J}$, $\mathbf{Y}_{(3)} \in \mathbb{R}^{I J \times M}$ the unfoldings of $\mathcal{Y}$ along its three modes. Then, at a given iteration $k$, the alternated optimization strategy reads explicitly

$$ \mathbf{A}^{k+1} \leftarrow \underset{\mathbf{A} \in \mathbb{R}_+^{I \times RL}}{\operatorname{argmin}} \|\mathbf{Y}_{(1)} - (\mathbf{B}^k \odot \mathbf{P}\tilde{\mathbf{S}}^k)\mathbf{A}^\top\|_{\mathrm{F}}^2, \quad (6) $$

$$ \mathbf{B}^{k+1} \leftarrow \underset{\mathbf{B} \in \mathbb{R}_+^{J \times RL}}{\operatorname{argmin}} \|\mathbf{Y}_{(2)} - (\mathbf{P}\tilde{\mathbf{S}}^k \odot \mathbf{A}^{k+1})\mathbf{B}^\top\|_{\mathrm{F}}^2, \quad (7) $$

$$ \mathbf{S}^{k+1} \leftarrow \underset{\mathbf{S} \in \mathbb{R}_S^{4 \times R}}{\operatorname{argmin}} \|\mathbf{Y}_{(3)} - (\mathbf{A}^{k+1} \odot \mathbf{B}^{k+1})(\mathbf{P}\tilde{\mathbf{S}})^\top\|_{\mathrm{F}}^2, \quad (8) $$

repeated until convergence. Note that to simplify expressions in (6)–(8), we used the notation $\tilde{\mathbf{S}} := \mathbf{S} \text{blkdiag}(1_L^{I \times 1_L}) \in \mathbb{R}^{4 \times LR}$ where each column of $\mathbf{S}$ is repeated $L$ times.

The outer loop iterations (6)–(8) define the backbone of the algorithms for Stokes-BTD studied in this paper. Two iterative methods – proximal gradient and alternating direction method of multipliers (ADMM) – are considered to solve the three convex subproblems. This results in two standard (constrained) tensor factorization algorithms known as block-proximal gradient (BPG) [15] and AO-ADMM [16], which have to be tailored to impose the physical constraints of Stokes-BTD. We implement BPG for Stokes-BTD following guidelines presented in [17], including FISTA-like acceleration and backtracking step-size procedures to ensure sufficient decay of the objective function. The implementation of AO-ADMM is standard [16, 18] with the penalty parameter $\rho$ being determined a priori and subproblem convergence being assessed by monitoring primal and dual residues. Both BPG and AO-ADMM offer guarantees of convergence to a stationary point of the cost function (5).

Explicit updates of BPG and AO-ADMM for the Stokes-BTD optimization problem (5) can be easily derived from the original papers [15, 16] and are omitted for space considerations. However, both algorithms require the computation of a projection operator for each constraint which we provide below for completeness. The projection of matrices $\mathbf{A}$ and $\mathbf{B}$ onto the set of matrices with non-negative entries is classic, i.e., $\Pi_{\mathbf{A}}(\mathbf{A}) = \max(0, \mathbf{A})$ and $\Pi_{\mathbf{B}}(\mathbf{A}) = \max(0, \mathbf{B})$ with the maximum being applied elementwise. In contrast, the projection of $\mathbf{S}$ onto $\mathbb{R}_S^{4 \times R}$ is less common. It consists in projecting each column $\mathbf{s}_1, \dots, \mathbf{s}_R$ of $\mathbf{S}$ onto $\mathbb{R}_S^4$. As detailed in Section 2.1 this is equivalent to projecting the 2-by-2 complex matrix $\mathbf{J}(\mathbf{s}_r)$ given by (1) onto the set of positive semidefinite matrices $\mathbb{C}_+^{2 \times 2}$. Let $\{\eta_i, \mathbf{v}_i\}$ represent the eigenvalues and eigenvectors of $\mathbf{J}(\mathbf{s}_r)$ such that $\mathbf{J}(\mathbf{s}_r) = \sum_{i=1}^2 \eta_i \mathbf{v}_i \mathbf{v}_i^\mathrm{H}$, where $\mathbf{v}_i^\mathrm{H}$ is the conjugate-transpose of $\mathbf{v}_i$. The projection of $\mathbf{J}(\mathbf{s}_r)$ onto $\mathbb{C}_+^{2 \times 2}$ reads

$$ \Pi_{\mathbb{C}_+^{2 \times 2}}(\mathbf{J}(\mathbf{s}_r)) = \sum_{i=1}^2 \max(0, \eta_i) \mathbf{v}_i \mathbf{v}_i^\mathrm{H} = \begin{bmatrix} \alpha & \gamma \\ \bar{\gamma} & \beta \end{bmatrix} \quad (9) $$

where $\alpha, \beta \in \mathbb{R}$ and $\gamma \in \mathbb{C}$ can be determined explicitly. Since $\mathbf{s}_r \mapsto \mathbf{J}(\mathbf{s}_r)$ is bijective, the projection of a single Stokes vector $\mathbf{s}_r$ onto $\mathbb{R}_S^4$ can be expressed from (9) as

$$ \Pi_{\mathbb{R}_S^4}(\mathbf{s}_r) = [(\alpha + \beta), (\alpha - \beta), 2 \operatorname{Re}(\gamma), 2 \operatorname{Im}(\gamma)]^\top \in \mathbb{R}_S^4, \quad (10) $$

and as a result, the projection of $\mathbf{S}$ onto $\mathbb{R}_S^{4 \times R}$ is obtained columnwise as $\Pi_{\mathbf{S}}(\mathbf{S}) = [\Pi_{\mathbb{R}_S^4}(\mathbf{s}_1) \cdots \Pi_{\mathbb{R}_S^4}(\mathbf{s}_R)]$.

## 3.3. Initialization strategies

Since the optimization problem (5) is non-convex, the choice of initial points plays an important role in reconstruction performance of the algorithms. In this paper, we consider two initialization strategies: 1) random initialization of the factor matrices $\mathbf{A}, \mathbf{B}, \mathbf{S}$ and 2) a $k$-means based approach exploiting an initial guess of $\mathcal{X}$.

Random initialization draws the entries of matrices $\mathbf{A}$ and $\mathbf{B}$ from i.i.d. uniform distribution on $[0, 1]$. For $\mathbf{S}$, we draw at random each column $\mathbf{s}_r$ as follows. We first generate a random complex Gaussian vector $\mathbf{b} \in \mathbb{C}^2$ with unit variance and form the rank-one matrix $\mathbf{b} \mathbf{b}^\mathrm{H} \in \mathbb{C}_+^{2 \times 2}$. The vector $\mathbf{s}_r$ is then obtained using the inverse mapping to $\mathbb{R}_S^4$ given by (10). This ensures that the initial points satisfy the physical constraints of the Stokes-BTD model.

The second initialization strategy starts by computing the unconstrained least square solution $\hat{\mathcal{X}}_{\text{ULS}} = \mathcal{Y}_{\bullet 3} \mathbf{P}^\mathrm{T}$. A $k$-means algorithm is then applied on the fibers of $\hat{\mathcal{X}}_{\text{ULS}}$ with a number of clusters set to the number of terms $R$ of the decomposition. Centroids of each cluster define the initial values of each column of $\mathbf{S}$, where projection $\Pi_{\mathbf{S}}$ was applied to ensure that $(\mathbf{S})$ is met. Then, for each $r$, one draws a random array $\mathbf{U}_r \in \mathbb{R}^{I \times J}$ with i.i.d. uniform entries in $[0, 1]$ and the corresponding feature map identified by the $k$-means algorithm is used as a binary mask for non-zero values. Finally, a non-negative matrix factorization (NMF) with rank $L$ is applied to $\mathbf{U}_r$ to get initial guesses for matrices $\mathbf{A}_r$ and $\mathbf{B}_r$.

# 4. NUMERICAL EXPERIMENTS AND DISCUSSION

## 4.1. Comparison of algorithms for Stokes-BTD

We first investigate the behavior of BPG and AO-ADMM algorithms for Stokes-BTD. We fix $L = 6$ and $R = 5$ and generate a synthetic ground-truth tensor $\mathcal{X}_0 = [\mathbf{A}, \mathbf{B}, \mathbf{S}]_{\text{SBTD}} \in \mathbb{R}^{100 \times 100 \times 4}$ where factor matrices $\mathbf{A}, \mathbf{B}, \mathbf{S}$ are drawn at random and satisfy
<figure><img src="image_5.png" /><figcaption>**Fig. 3:** Performance analysis of BPG and AO-ADMM algorithms. Thick lines correspond to averages over 50 independent runs.</figcaption></figure>

the physical constraints of the Stokes-BTD model. The instrument matrix $\mathbf{P}$ is chosen to represent the standard deterministic Stokes measurement scheme [11] which uses $M = 6$ polarimetric projections. Finally, we add i.i.d. Gaussian noise following the forward model (2) such that the resulting signal-to-noise ratio (SNR) is 25 dB, where SNR = $\|\mathcal{X}_{\bullet 3}\mathbf{P}\|_F^2/(IJM\sigma^2)$ with $\sigma^2$ being the noise variance. The noisy measurement tensor $\mathcal{Y} \in \mathbb{R}^{100 \times 100 \times 6}$ is fixed for all experiments. Both algorithms were run for 1000 outer loop iterations, starting from the same points (for fair comparisons). Two figures of merits are used to evaluate the performance of algorithms. The first is the normalized mean squared error (NMSE) defined

as NMSE($\hat{\mathbf{A}}, \hat{\mathbf{B}}, \hat{\mathbf{S}}$) = $\frac{1}{R} \sum_{r=1}^{R} \frac{\|\mathbf{A}_r \mathbf{B}_r^\top \circ \mathbf{s}_r - \hat{\mathbf{A}}_r \hat{\mathbf{B}}_r^\top \circ \hat{\mathbf{s}}_r\|_F^2}{\|\mathbf{A}_r \mathbf{B}_r^\top \circ \mathbf{s}_r\|_F^2}$, which quantifies the reconstruction quality of the $R$ block terms. The second is the value of the cost function (5) which monitors convergence.

Fig. 3 depicts the performances of the proposed algorithms with the initialization strategies proposed on Section 3.3, averaged over 50 independent Monte Carlo runs. For both initialization strategies, AO-ADMM provides faster convergence and smaller error, both in terms of NMSE and cost function value. Moreover, we observe that $k$-means initialization provides faster decrease of the NMSE for both BPG and AO-ADMM algorithms, while cost functions values show no significant differences. Runtimes for 1000 iterations are comparable, with an average 9.3s for BPG and 13.7s for AO-ADMM on a 2021 M1 Max MacbookPro using Python 3.11.

## 4.2. Validation on real data and discussion

We now demonstrate the ability of the Stokes-BTD approach to provide a meaningful decomposition of real-world polarimetric images. We consider data [19] acquired in a polarized microscopy experiment of a mouse brain tissue. This setting uses the minimal number $M = 4$ of polarimetric projections such that the instrument matrix $\mathbf{P} \in \mathbb{R}^{4 \times 4}$ is known and invertible. After downsampling, we obtain an intensity measurement tensor $\mathcal{Y}$ of size $250 \times 250 \times 4$. Following our analysis of Section 4.1, we only consider the AO-ADMM algorithm and use $k$-means initialization. The algorithm was run for a maximal number of 1000 outer loop iterations, stopping whenever the relative improvement of the solution did not exceed $10^{-4}$.

We first analyze the choice of $R$ and $L$ of the Stokes-BTD model. The value of $R$ depends on a priori information of the data i.e., the expected number of different constant polarization states. For this data, $R = 4$ elementary polarimetric regions can be easily identified, e.g., by $k$-means analysis, so that we focus on the choice of $L$ instead. Intuitively, a larger rank for $\mathbf{A}_r$ and $\mathbf{B}_r$ enables

<figure><img src="image_6.png" /><figcaption>**Fig. 4:** Analysis of the tradeoff reconstruction error / runtime for the different values of the rank $L$ in Stokes-BTD.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 5:** Stokes-BTD for mouse brain tissue data [19] with $R = 4$ and two values of $L = 10$ (top) and $L = 30$ (bottom). Each row shows the corresponding low-rank activation map $\mathbf{A}_r \mathbf{B}_r^\top$ associated with a single polarization state $\mathbf{s}_r$ (represented as a polarization ellipse).</figcaption></figure>

Stokes-BTD to represent more complex spatial information, thus improving the reconstruction error $\|\mathcal{Y} - \mathcal{X}_{\bullet 3}\mathbf{P}\|_F^2$. However, the dimensions of the factor matrices grow linearly with $L$, leading to an increased computational burden. Fig. 4 illustrates this trade-off for values of $L$ ranging from $L = 3$ to $L = 100$ (maximum value guaranteeing identifiability of Stokes-BTD). From Fig. 4a) we observe that the reconstruction error decreases rapidly with $L$, showing no significant improvements after $L = 50$. In contrast, Fig. 4b) shows that the average AO-ADMM runtime (over 50 realizations) increases with $L$ and might become prohibitive for large $L$ values.

Fig. 5 shows the decomposition of the real data with Stokes-BTD for $L = 10$ and $L = 30$. As expected, a greater value of $L$ allows the model to represent more complex spatial information. We also note that the polarization states are almost identical for both $L$ values, which shows the robustness of Stokes-BTD with respect to polarization information. From a physical standpoint, we observe that Stokes-BTD is able to recover the strongly circularly polarized background (first column), together with three other regions (last three columns) with different levels of linear birefringence (i.e. specific polarimetric properties that encode distinct biological features).

Real data experiments in polarimetric microscopy highlight the Stokes-BTD approach's potential for interpretable reconstructions and segmentation of Stokes polarimetric data. This physically-informed low-rank tensor model, leveraging tensor algebra, exhibits scalable resolution while enjoying broad identifiability guarantees. The AO-ADMM algorithm coupled with a $k$-means initialization approach provides an efficient reconstruction strategy. This new tool offers promising research avenues for practical applications relying on polarimetric imaging, e.g., cancerous tissue characterization [1].
5. REFERENCES

[1] Chao He, Honghui He, Jintao Chang, Binguo Chen, Hui Ma, and Martin J Booth, “Polarisation optics for biomedical and clinical applications: a review,” *Light: Science & Applications*, vol. 10, no. 1, pp. 194, 2021.

[2] J Scott Tyo, Dennis L Goldstein, David B Chenault, and Joseph A Shaw, “Review of passive imaging polarimetry for remote sensing applications,” *Applied optics*, vol. 45, no. 22, pp. 5453–5469, 2006.

[3] Maria Losurdo, Michael Bergmair, Giovanni Bruno, Denis Cattelan, Christoph Cobet, Antonello De Martino, Karsten Fleischer, Zorana Dohcevic-Mitrovic, Norbert Esser, Melanie Galliet, et al., “Spectroscopic ellipsometry and polarimetry for materials and systems analysis at the nanometer scale: state-of-the-art, potential, and perspectives,” *Journal of Nanoparticle Research*, vol. 11, pp. 1521–1554, 2009.

[4] José J Gil and Razvigor Ossikovski, *Polarized light and the Mueller matrix approach*, CRC press, 2022.

[5] Jihad Zallat and Christian Heinrich, “Polarimetric data reduction: a Bayesian approach,” *Optics express*, vol. 15, no. 1, pp. 83–96, 2007.

[6] Giorgos Sfikas, Christian Heinrich, Jihad Zallat, Christophoros Nikou, and Nikos Galatsanos, “Recovery of polarimetric Stokes images by spatial mixture models,” *Journal of the Optical Society of America A*, vol. 28, no. 3, pp. 465, Mar. 2011.

[7] Sylvain Faisan, Christian Heinrich, François Rousseau, Alex Lallement, and Jihad Zallat, “Joint filtering estimation of Stokes vector images based on a nonlocal means approach,” *Journal of the Optical Society of America A*, vol. 29, no. 9, pp. 2028, Sept. 2012.

[8] John R Valenzuela and Jeffrey A Fessler, “Joint reconstruction of Stokes images from polarimetric measurements,” *Journal of the Optical Society of America A*, vol. 26, no. 4, pp. 962–968, 2009.

[9] Carole Le Guyader, Samia Ainouz, and Stéphane Canu, “A physically admissible Stokes vector reconstruction in linear polarimetric imaging,” *Journal of Mathematical Imaging and Vision*, pp. 1–26, 2023.

[10] Tamara G Kolda and Brett W Bader, “Tensor decompositions and applications,” *SIAM review*, vol. 51, no. 3, pp. 455–500, 2009.

[11] Beth Schaefer, Edward Collett, Robert Smyth, Daniel Barrett, and Beth Fraher, “Measuring the Stokes polarization parameters,” *American Journal of Physics*, vol. 75, no. 2, pp. 163–168, 2007.

[12] Nirmal Keshava and John F Mustard, “Spectral unmixing,” *IEEE Signal Processing Magazine*, vol. 19, no. 1, pp. 44–57, 2002.

[13] Lieven De Lathauwer, “Decompositions of a higher-order tensor in block terms—Part II: Definitions and uniqueness,” *SIAM Journal on Matrix Analysis and Applications*, vol. 30, no. 3, pp. 1033–1066, 2008.

[14] Saulo Cardoso Barreto, Julien Flamant, Sebastian Miron, and David Brie, “Décomposition tensorielle de rang faible pour l’imagerie de Stokes polarimétrique,” in *GRETSI 2023-XXIXème Colloque Francophone de Traitement du Signal et des Images*, 2023, (in French).

[15] Yangyang Xu and Wotao Yin, “A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion,” *SIAM Journal on Imaging Sciences*, vol. 6, no. 3, pp. 1758–1789, Jan. 2013.

[16] Kejun Huang, Nicholas D. Sidiropoulos, and Athanasios P. Liavas, “A flexible and efficient algorithmic framework for constrained matrix and tensor factorization,” *IEEE Transactions on Signal Processing*, vol. 64, no. 19, pp. 5052–5065, Oct. 2016.

[17] Amir Beck, *First-order methods in optimization*, MOS-SIAM Series on Optimization. SIAM, 2017.

[18] Stephen Boyd, Neal Parikh, Eric Chu, Borja Peleato, Jonathan Eckstein, et al., “Distributed optimization and statistical learning via the alternating direction method of multipliers,” *Foundations and Trends® in Machine learning*, vol. 3, no. 1, pp. 1–122, 2011.

[19] Cameron Foltz and Shalin Mehta, “Dataset for demonstration of quantitative label-free imaging with phase and polarization,” June 2021, available online at https://dx.doi.org/10.5281/zenodo.5178487.