Quaternion Non-negative Matrix Factorization: definition, uniqueness and algorithm

TL;DR

This paper introduces quaternion non-negative matrix factorization (QNMF), a novel method for polarized signal analysis that leverages quaternion algebra and physical constraints, offering improved uniqueness and applicability in polarimetric imaging.

Contribution

The paper defines QNMF, establishes its uniqueness conditions, and proposes an efficient quaternion ALS algorithm, extending NMF to polarized signals with enhanced theoretical guarantees.

Findings

QNMF generalizes NMF for polarized signals.

The QALS algorithm efficiently solves QNMF.

Numerical experiments validate QNMF's effectiveness.

Abstract

This article introduces quaternion non-negative matrix factorization (QNMF), which generalizes the usual non-negative matrix factorization (NMF) to the case of polarized signals. Polarization information is represented by Stokes parameters, a set of 4 energetic parameters widely used in polarimetric imaging. QNMF relies on two key ingredients: (i) the algebraic representation of Stokes parameters thanks to quaternions and (ii) the exploitation of physical constraints on Stokes parameters. These constraints generalize non-negativity to the case of polarized signals, encoding positive semi-definiteness of the covariance matrix associated which each source. Uniqueness conditions for QNMF are presented. Remarkably, they encompass known sufficient uniqueness conditions from NMF. Meanwhile, QNMF further relaxes NMF uniqueness conditions requiring sources to exhibit a certain zero-pattern, by leveraging the complete polarization information. We introduce a simple yet efficient algorithm called quaternion alternating least squares (QALS) to solve the QNMF problem in practice. Closed-form quaternion updates are obtained using the recently introduced generalized HR calculus. Numerical experiments on synthetic data demonstrate the relevance of the approach. QNMF defines a promising generic low-rank approximation tool to handle polarization, notably for blind source separation problems arising in imaging applications.

Figures

Tables

Table 1. TABLE I: Matrix functions f : ℍ M × N → ℍ : 𝑓 → superscript ℍ 𝑀 𝑁 ℍ f:\mathbb{H}^{M\times N}\rightarrow\mathbb{H} and corresponding gradients. Matrices 𝐀 1 , 𝐀 2 subscript 𝐀 1 subscript 𝐀 2 \mathbf{A}_{1},\mathbf{A}_{2} are quaternion-valued with appropriate dimensions.

$f\left(𝐐\right)$	${\nabla }_{\overline{𝐐}}f$
$\mathrm{tr}{𝐀}_1{\mathrm{𝐐𝐀}}_2$	$-\frac{1}{2}{𝐀}_1^{\top }{𝐀}_2^{†}$
$\mathrm{tr}{𝐀}_1{𝐐}^{†}{𝐀}_2$	$\mathrm{Re}\left[{𝐀}_2\right]{𝐀}_1$
$\mathrm{tr}{𝐀}_1{𝐐}^{†}{\mathrm{𝐐𝐀}}_2$	$\mathrm{Re}\left[{\mathrm{𝐐𝐀}}_2\right]{𝐀}_1-\frac{1}{2}{\left({𝐀}_1{𝐐}^{†}\right)}^{\top }{𝐀}_2^{†}$

Equations

$$\Phi =\frac{\text{intensity of polarized part}}{\text{total intensity}}=\frac{\sqrt{S_1^2+S_2^2+S_3^2}}{S_0}.$$

$$w=S_0+\mathbf{i}{S}_3+\mathbf{j}{S}_1+\mathbf{k}{S}_2\in \mathbb{H}.$$

$$w=I+I\Phi \mathbf{\mu }$$

$$S_0\ge 0\text{ and }{S}_1^2+S_2^2+S_3^2\le S_0^2.$$

$$\mathbf{J}=\frac{1}{2}\left[\begin{array}{cc}{S}_0+S_1& S_2+\mathbf{i}{S}_3\\ S_2-\mathbf{i}{S}_3& S_0-S_1\end{array}\right]\in {\mathbb{C}}^{2\times 2}.$$

$$(\mathcal{S})\iff \{\begin{array}{ll}\mathrm{tr}\mathbf{J}& \ge 0\\ \det \mathbf{J}& \ge 0\end{array}.$$

$${\mathbb{H}}_{\mathcal{S}}\triangleq \left\{q\in \mathbb{H}|\mathrm{Re}q\ge 0\text{ and }|\mathrm{Im}q{|}^2\le (\mathrm{Re}q{)}^2\right\}.$$

$$\begin{array}{rl}x({\nu }_m,u_n)& =S_0({\nu }_m,u_n)+\mathbf{i}{S}_3({\nu }_m,u_n)\\ & +\mathbf{j}{S}_1({\nu }_m,u_n)+\mathbf{k}{S}_2({\nu }_m,u_n)\in {\mathbb{H}}_{\mathcal{S}}\end{array}$$

$$x({\nu }_m,u_n)=\sum _{p=1}^P{w}_p({\nu }_m)h_p(u_n).$$

$$\mathbf{X}=\mathbf{WH}$$

$$\mathbf{X}=\mathbf{WH}\iff \{\begin{array}{ll}\mathrm{Re}\mathbf{X}=\left[\mathrm{Re}\mathbf{W}\right]\mathbf{H}& \text{(NMF)}\\ \mathrm{Im}\mathbf{X}=\left[\mathrm{Im}\mathbf{W}\right]\mathbf{H}& \text{(polarization)}\end{array}.$$

$$\tilde{\mathbf{W}}$$

$$\tilde{\mathbf{H}}$$

$$\tilde{\mathbf{W}}\in {\mathbb{H}}_{\mathcal{S}}^{M\times P}\text{ and }\tilde{\mathbf{H}}\in {\mathbb{R}}_{+}^{P\times N}.$$

$$\mathbf{T}=\mathrm{diag}\left(t_1,t_2,\dots ,t_P\right),t_i>0,1\le i\le P$$

$$\mathbf{T}\text{ is a permutation matrix}$$

$$(\mathbf{W}{)}_{mp}=w_{mp}\in {\mathbb{H}}_{\mathcal{S}}\text{ and }(\mathbf{H}{)}_{pn}=h_{pn}\ge 0$$

$$w_{mp}=I_{mp}+I_{mp}{\Phi }_{mp}{\mathbf{\mu }}_{mp}.$$

$$\mathbf{T}(\alpha ,\beta )=\left[\begin{array}{cc}1-\alpha & \beta \\ \alpha & 1-\beta \end{array}\right].$$

$${\mathbf{T}}^{-1}(\alpha ,\beta )=\frac{1}{1-\alpha -\beta }\left[\begin{array}{cc}1-\beta & -\beta \\ -\alpha & 1-\alpha \end{array}\right].$$

$${\tilde{w}}_{m1}$$

$${\tilde{w}}_{m2}$$

$${\tilde{h}}_{1n}$$

$${\tilde{h}}_{2n}$$

$$\{\begin{array}{ll}\mathrm{Re}{\tilde{w}}_{m1}& \ge 0\\ \mathrm{Re}{\tilde{w}}_{m2}& \ge 0\end{array},$$

$$\{\begin{array}{l}|\mathrm{Im}{\tilde{w}}_{m1}{|}^2\le {\left(\mathrm{Re}{\tilde{w}}_{m1}\right)}^2\\ |\mathrm{Im}{\tilde{w}}_{m2}{|}^2\le {\left(\mathrm{Re}{\tilde{w}}_{m2}\right)}^2\end{array}.$$

$${\mathcal{I}}_{\alpha }^{\mathrm{QNMF}}$$

$${\mathcal{I}}_{\beta }^{\mathrm{QNMF}}$$

$${\mathcal{I}}_{\alpha }^{\mathrm{NMF}}=({\alpha }_{\text{min}}^{\text{NMF}},{\alpha }_{\text{max}}^{\text{NMF}}]\text{ and }{\mathcal{I}}_{\beta }^{\mathrm{NMF}}=({\beta }_{\text{min}}^{\text{NMF}},{\beta }_{\text{max}}^{\text{NMF}}].$$

$${\alpha }_{\text{min}}^{\text{NMF}}=-\underset{m\in {\mathbb{M}}_1}{\min }\frac{I_{m1}}{I_{m2}-I_{m1}},{\beta }_{\text{min}}^{\text{NMF}}=-\underset{m\in {\mathbb{M}}_2}{\min }\frac{I_{m2}}{I_{m1}-I_{m2}}$$

$${\alpha }_{\text{max}}^{\text{NMF}}=\underset{n}{\min }\frac{h_{2n}}{h_{1n}+h_{2n}},{\beta }_{\text{max}}^{\text{NMF}}=\underset{n}{\min }\frac{h_{1n}}{h_{1n}+h_{2n}}.$$

$$\mathrm{Re}{\tilde{w}}_{m1}$$

$$\mathrm{Im}{\tilde{w}}_{m1}$$

Full text

I Introduction

Polarization information is essential to many fields ranging from seismology [2], optics [3] to gravitational astronomy [4], among others. Polarization encodes the geometry of wave oscillations. It carries crucial morphological and physical insights about the medium which transmitted or reflected the polarized wave. For instance, polarization permits the detection of hidden targets in remote sensing [5, 6], it provides insights into solar magnetic fields [7] and reveals the anisotropy of biological tissues [8]. In summary, polarimetric imaging proves to be useful in a wide range of situations as it provides unique features of the observed scene, unaccessible to conventional intensity imaging.

Recent years have seen an increased interest in exploiting polarization diversity in hyperspectral imaging systems [9, 10, 11]. This modality, called spectro-polarimetric imaging, was first popularized in astronomy [12]. In full generality, it consists in acquiring, for a collection of wavelengths and locations, a 4-dimensional real vector gathering the 4 Stokes parameters. These energetic parameters are widely used in optics [13] to describe the complete properties of light (intensity and polarization). The joint acquisition of spatial, spectral and polarization diversities raises important and original challenges, in particular regarding blind polarized source separation problems.

Spectro-polarimetric data can be represented as a third order real-valued tensor, where the rows, columns and fibers correspond to frequency, space and polarization diversities, respectively. Note that, since polarization is described here by Stokes parameters, the fiber is of size 4. Solving the blind polarized source separation problem using usual low-rank approximation techniques raises two major issues. First, Stokes parameters being energetic quantities, they obey structural constraints [13] which need to be taken into account to guarantee the physical interpretation of the recovered rank-1 terms. Second, low-rank decompositions of spectro-polarimetric images such as the Canonical Polyadic Decomposition (CPD) [14] necessarily yield rank-1 terms with constant polarization. It thus restricts the physical validity of usual tensor decompositions to specific settings, e.g. narrow band sources with constant polarization. The general case of wideband sources with frequency-dependent polarization, as observed in spectro-polarimetric data, requires the development of new low-rank approximation tools to deal with the specificities of polarization.

To this aim, this paper introduces quaternion non-negative matrix factorization (QNMF). It extends the concept of non-negative matrix factorization (NMF) to polarized sources. It relies on two key ingredients: (i) the algebraic representation of polarization information using quaternions and (ii) the exploitation of physical constraints on Stokes parameters, which generalize non-negativity to the case of polarized signals. Non-negativity is to be understood here in a vectorial sense: it encodes the non-negativity, or positive semi-definiteness, of the 2-by-2 complex covariance matrix associated with each polarized source, as will be explained later.

This generalized non-negativity constraint on the source factor makes QNMF fundamentally different from NMF. In fact, the three Stokes parameters describing polarization are not required to be real non-negative: instead, Stokes parameters must lie inside a particular 4D second-order convex cone. This means that QNMF cannot be simply reformulated as an augmented real NMF problem by stacking Stokes parameters in a taller matrix. Instead, QNMF preserves the intrinsic structure of Stokes parameters by leveraging their quaternion representation: it allows the associated generalized non-negativity constraints to be handled straightforwardly.

QNMF encodes a quaternion blind source separation problem. However, unlike quaternion ICA [15, 16] which relies on a standard quaternion linear mixing model and assumes statistical independence of the sources, QNMF does not require any statistical model; it defines a quaternion linear mixing model with non-negativity constraints only. Such non-negativity constraints are essential to physical interpretation: they guarantee the general relevance of QNMF to solve blind unmixing of spectro-polarimetric data, just like NMF offers a standard approach towards hyperspectral image unmixing [17]. More generally, recent works [18, 19, 20, 21] have shown that Stokes parameters provide a natural reparameterization for the second-order properties of bivariate signals. As a result, QNMF defines a new low-rank approximation tool for arbitrary bivariate signals, such as current velocities in oceanography [22] or polarized gravitational waves [23].

Our contributions are threefold. First, we introduce the QNMF problem from a spectro-polarimetric linear mixing model perspective and demonstrate how it generalizes the usual NMF. Second, we study uniqueness conditions for the QNMF problem and reveal the key disambiguating role played by polarization. QNMF is shown to encompass known NMF sufficient uniqueness conditions, while relaxing NMF necessary uniqueness conditions. Finally, we provide a simple yet efficient algorithm in the quaternion domain that effectively solves the QNMF problem. We believe that quaternions, beyond offering an elegant way to handle polarization diversity, permit numerous theoretical (e.g. uniqueness conditions) and methodological (e.g. algorithmic design) developments that would have been cumbersome to obtain otherwise.

This paper is organized as follows. Section II collects necessary material regarding quaternions, polarization and associated constraints. Section III defines the concept of quaternion non-negative matrix factorization (QNMF). Section IV studies uniqueness conditions for QNMF, starting with the simpler 2 sources case, followed by a discussion of the general case. Section V introduces the quaternion alternating least squares (QALS) algorithm for solving the QNMF problem in practice. Section VI provides numerical evidence to the relevance of the approach, while Section VII presents concluding remarks. Appendices gather technical details and proofs.

II Preliminaries

This paper relies on the algebraic representation of polarization using quaternions together with the exploitation of structural non-negativity constraints. To this end, Section II-A briefly reviews quaternion algebra and Section II-B describes polarization using Stokes parameters. Section II-C then introduces the quaternion representation of Stokes parameters, whereas Section II-D details structural non-negative like constraints on Stokes parameters.

II-A Quaternions

Quaternions $\mathbb{H}$ forms a 4-dimensional normed division algebra over the real numbers $\mathbb{R}$ with canonical basis $\{1,\bm{i},\bm{j},\bm{k}\}$ where $\bm{i},\bm{j},\bm{k}$ are imaginary units $\bm{i}^{2}=\bm{j}^{2}=\bm{k}^{2}=-1$ such that $\bm{i}\bm{j}=\bm{k},\>\bm{i}\bm{j}=-\bm{j}\bm{i}$ and $\bm{i}\bm{j}\bm{k}=-1$. Importantly, quaternion multiplication is noncommutative, i.e. for $q,q^{\prime}\in\mathbb{H}$, $qq^{\prime}\neq q^{\prime}q$ in general. Any quaternion $q\in\mathbb{H}$ reads in Cartesian form $q=a+b\bm{i}+c\bm{j}+d\bm{k}$ where $a,b,c,d\in\mathbb{R}$. The real or scalar part of $q$ is $\mathrm{Re}\,q=a$ and its imaginary or vector part is $\mathrm{Im}\,q=b\bm{i}+c\bm{j}+d\bm{k}$. We note $\mathrm{Im}_{\bm{i}}\,q=b$, $\mathrm{Im}_{\bm{j}}\,q=c$ and $\mathrm{Im}_{\bm{k}}\,q=d$. When $\mathrm{Re}\,q=0$, $q$ is said to be pure imaginary, or simply pure. The quaternion conjugate of $q$ is denoted by $\overline{q}=\mathrm{Re}\,q-\mathrm{Im}\,q$. The modulus of $q$ is $|q|=\sqrt{q\overline{q}}=\sqrt{\overline{q}q}$.

A quaternion matrix $\mathbf{Q}\in\mathbb{H}^{M\times N}$ has elements $\mathbf{Q}_{mn}=q_{mn}\in\mathbb{H}$. Notation $\overline{\mathbf{Q}}$ denotes the entry-wise conjugate of $\mathbf{Q}$ whereas its conjugate-transpose is denoted by $\mathbf{Q}^{\dagger}=\overline{\mathbf{Q}}^{\top}$. Usual matrix operations require special care due to non-commutativity of quaternion multiplication: e.g. given two matrices $\mathbf{Q}_{1}\in\mathbb{H}^{M\times N}$ and $\mathbf{Q}_{2} \in \mathbb{H}^{N\times P}$, in general one has $\overline{\mathbf{Q}_{1}\mathbf{Q}_{2}}\neq\overline{\mathbf{Q}}_{1}\>\overline{\mathbf{Q}}_{2}$ and $(\mathbf{Q}_{1}\mathbf{Q}_{2})^{\top}\neq\mathbf{Q}_{2}^{\top}\mathbf{Q}_{1}^{\top}$. However, $(\mathbf{Q}_{1}\mathbf{Q}_{2})^{\dagger}=\mathbf{Q}_{2}^{\dagger}\mathbf{Q}_{1}^{\dagger}$ is always true. More details on quaternions and their properties can be found in [24]. See also [25, 26, 27, 28, 29] and references therein for more on quaternion matrix algebra.

II-B Polarization and Stokes parameters

For electromagnetic (EM) waves such as light traveling in free space, the instantaneous electric field vector $\mathbf{E}(t)$ lives in the 2D plane orthogonal to the direction of propagation. In such transverse waves, the end tip of $\mathbf{E}(t)\in\mathbb{R}^{2}$ describes a two-dimensional elliptical trajectory whose properties can vary with time, space or frequency. This phenomenon is called polarization: it describes the geometric nature of wave oscillations. This behaviour is not limited to electromagnetic radiation: transverse elastic waves and gravitational waves also exhibit polarization.

In optics [3, 30, 13], the polarization state of light is commonly described by four Stokes parameters $S_{0},S_{1},S_{2},S_{3}$. These real-valued parameters are energetic quantities that can be measured in experiments [31, 32]. The first Stokes parameter $S_{0}\geq 0$ is classical and measures the total intensity of light, i.e. the sum of intensities from the polarized and unpolarized parts of light. The three remaining Stokes parameters $S_{1},S_{2},S_{3}$ describe the polarized part of light, i.e. the properties of the polarization ellipse encoded by the trajectory of the electric field vector. The degree of polarization $\Phi\in[0,1]$ rules the relative contribution of the polarized part to the total intensity:

[TABLE]

When $\Phi=1$, light is said to be fully polarized, whereas for $\Phi=0$ it is said unpolarized. For $\Phi\in(0,1)$, light is said to be partially polarized.

II-C Quaternion representation of Stokes parameters

The Stokes vector $(S_{0},S_{1},S_{2},S_{3})^{\top}\in\mathbb{R}^{4}$ can be conveniently expressed using quaternion algebra as

[TABLE]

This algebraic representation of Stokes parameters has been long standing in the field of optics [33, 34, 35, 36, 37]. More recently [21, 18], it has been shown that (2) exhibits a general relevance for arbitrary bivariate signals. Our choice for the ordering of Stokes parameters in (2) adopts that of [21].

The quaternion representation of Stokes parameters enables straightforward physical and geometric interpretations while simplifying computations [18]. Eq. (2) can be rewritten for convenience as

[TABLE]

where $I:=S_{0}\geq 0$ is the total intensity and $\bm{\mu}$ is the polarization axis – a pure unit quaternion such that $\bm{\mu}^{2}=-1$. With this notation, the quantity $\Phi\bm{\mu}$ is a pure quaternion that can be identified with a vector of $\mathbb{R}^{3}$ inside the Poincaré sphere [3] which completely encodes the polarization state of light. The expression (3) also emphasizes one of the benefits of the quaternion formalism: it allows a straightforward separation between pure energetic information, conveyed by $\mathrm{Re}\,w=I$, and the geometric vector information encoded by $\mathrm{Im}\,w=I\Phi\bm{\mu}$.

II-D Polarization constraint and quaternion non-negativity

Stokes parameters are energetic quantities which are subject to admissibility conditions. A vector $(S_{0},S_{1},S_{2},S_{3})^{\top}\in\mathbb{R}^{4}$ is admissible as a Stokes vector if it satisfies the following constraints $(\mathcal{S})$

[TABLE]

Those two constraints bear a strong physical interpretation. The first one is classical and indicates that the total intensity $S_{0}$ is a non-negative real quantity. The second one means that the intensity of the polarized part cannot be larger than the total intensity. As a result, these constraints equivalently encode the range $0\leq\Phi\leq 1$ of the degree of polarization (1). Note that ($\mathcal{S}$) does not impose at all that $S_{1},S_{2}$ and $S_{3}$ are non-negative.

From a mathematical perspective, ($\mathcal{S}$) defines a closed second-order convex cone in $\mathbb{R}^{4}$, denoted by $\mathbb{R}^{4}_{\mathcal{S}}$. Conditions ($\mathcal{S}$) are referred to as non-negativity constraints for polarized signals by the following argument. Given $(S_{0},S_{1},S_{2},S_{3})^{\top}\in\mathbb{R}^{4}$, consider the 2-by-2 Hermitian matrix111In optics, the matrix $\mathbf{J}$ is known as the polarization or coherency matrix [13, Section 1.4] associated with the bivariate electric field $\mathbf{E}$. It can be statistically defined as a covariance matrix [13, Section 1.4.1], from which Hermitian positive semi-definiteness follows.:

[TABLE]

This mapping is bijective. Then, a quick check shows that imposing positive semi-definiteness of $\mathbf{J}$ – or simply, non-negativity of $\mathbf{J}$, is equivalent to ($\mathcal{S}$), i.e.

[TABLE]

This shows that $\mathbb{R}^{4}_{\mathcal{S}}$ and the set of 2-by-2 Hermitian non-negative matrices are isomorphic. This result will prove to be essential later on for imposing $\eqref{eq:polarizationConstraint}$ in practice, see Section V-A. Since the proposed approach relies on the quaternion representation (2) of Stokes parameters, we define by extension the set of non-negative quaternions $\mathbb{H}_{\mathcal{S}}\subset\mathbb{H}$ such that

[TABLE]

By definition, the sets $\mathbb{H}_{\mathcal{S}}$ and $\mathbb{R}^{4}_{\mathcal{S}}$ are isomorphic. They share the same properties: in particular, this means that $\mathbb{H}_{\mathcal{S}}$ is a closed second-order convex cone of $\mathbb{H}$.

III Quaternion Non-negative Matrix Factorization

This section introduces Quaternion Non-negative Matrix Factorization (QNMF), a new tool which generalizes the notion of Non-negative Matrix Factorization (NMF) to the case of polarized signals. It exploits two key features: (i) the polarization constraint ($\mathcal{S}$), which extends the notion of non-negativity to polarized signals and (ii) the algebraic representation (2) of Stokes parameters using quaternions. QNMF establishes a new and generic tool for low-rank approximations of polarized signals, with many potential applications in source separation or data completion.

III-A Definition

Without loss of generality, let us consider the setting of a typical spectro-polarimetric experiment. In such experiment, the four Stokes parameters $S_{i}(\nu_{m},u_{n})$, $i=0,1,2,3$ are acquired for a collection of sampled frequencies (or wavelengths) $(\nu_{m})_{1\leq m\leq M}$ and spatial locations $(u_{n})_{1\leq n\leq N}$. Using (2), the collected data can be written in quaternion form as

[TABLE]

for $1\leq m\leq M$, $1\leq n\leq N$. Note that spectro-polarimetric imaging includes standard hyperspectral imaging as a special case, in which only intensity measurements are performed i.e. $x(\nu_{m},u_{n})=S_{0}(\nu_{m},u_{n})$ .

A standard assumption in hyperspectral imaging is to consider that, at a given pixel $u_{n}$, data can be represented as a linearly weighted combination of $P$ elementary sources (or endmembers) [38, 17]. While being simple, the linear mixing model appears reasonable for many real-world settings [39, 40]. Generalizing this model to the case of polarized sources, we decompose spectro-polarimetric data (7) as the linear superposition of $P$ elementary polarized sources $w_{p}(\cdot)\in\mathbb{H}_{\mathcal{S}}$:

[TABLE]

The quantity $h_{p}(\cdot)\geq 0$ denotes the spatial activation coefficients attached to the source $p$. Considering that $h_{p}(\cdot)$ is a non-negative real quantity has two key implications. First, it means that $h_{p}(u_{n})$ can be interpreted as the abundance of the source $p$ at location $u_{n}$. Second, it ensures that, for any elementary source $w_{p}(\cdot)\in\mathbb{H}_{\mathcal{S}}$ the linear combination $\sum_{p=1}^{P}w_{p}(\nu_{m})h_{p}(u_{n})$ is in $\mathbb{H}_{\mathcal{S}}$ since $\mathbb{H}_{\mathcal{S}}$ is a convex cone.

Rewriting (8) into matrix form introduces Quaternion Non-Negative Matrix Factorization (QNMF) like

[TABLE]

where $\mathbf{X}\in\mathbb{H}_{\mathcal{S}}^{M\times N}$ is the data matrix with coefficients $(\mathbf{X})_{mn}=x_{mn}=x(\nu_{m},u_{n})$. The matrix $\mathbf{W}\in\mathbb{H}_{\mathcal{S}}^{M\times P}$ is the source matrix with coefficients $(\mathbf{W})_{mp}=w_{mp}=w_{p}(\nu_{m})$, whereas $\mathbf{H}\in\mathbb{R}_{+}^{P\times N}$ is called the activation matrix such that $(\mathbf{H})_{pn}=h_{pn}=h_{p}(u_{n})$.

III-B Relation with NMF

The QNMF problem extends the well-known NMF problem to the case of polarized signals. Compared to NMF, QNMF features a quaternion-valued sources factor $\mathbf{W}$ which exploits the constraint ($\mathcal{S}$) instead of the usual non-negativity constraint. Using quaternion algebra rules, the QNMF (9) can also be rewritten as

[TABLE]

Eq. (10) shows that QNMF can be seen as a co-factorization problem with common activation factor $\mathbf{H}$. The first factorization problem is an usual NMF on the real part of $\mathbf{X}$, i.e. on intensity data (Stokes parameter $S_{0}\geq 0$) only. The second one corresponds to a factorization problem on the imaginary part of $\mathbf{X}$ encoding polarization properties (Stokes parameters $S_{1},S_{2},S_{3}$). These two factorization problems are not independent for two reasons: (i) the activation factor $\mathbf{H}$ appears in both and (ii) for each coefficient $(m,p)$ of the source factor $\mathbf{W}$, the constraint ($\mathcal{S}$) connects the modulus of the imaginary part (polarization factorization problem) to its real part (NMF problem) by an inequality constraint.

The relationship (10) provides another illustration of how QNMF extends NMF to account for polarization diversity. It allows a precise quantification of the role played by the polarization information and its associated constraint ($\mathcal{S}$). This will be particularly useful for comparing uniqueness conditions of QNMF and NMF, see Section IV.

IV Uniqueness conditions for QNMF

This section deals with a fundamental question: upon which conditions on the source $\mathbf{W}$ and activation $\mathbf{H}$ factors is the QNMF $\mathbf{X}=\mathbf{W}\mathbf{H}$ unique? The uniqueness or identifiability question is central in any factorization or decomposition problem, as it is intimately related to the interpretability of the underlying model. Many authors have thus considered this question for the NMF problem, either by an analytical study of the properties of the NMF model [41, 42, 43, 44] or by the design of suitable identification criteria to guarantee identifiability [45, 46, 47, 48]. In this paper, we address the QNMF uniqueness question by a theoretical study of the properties of the QNMF model. To this aim, Section IV-A starts by introducing the intrinsic ambiguities of the QNMF model. We provide sufficient uniqueness conditions for QNMF, directly inherited from known NMF ones. To reveal the key role played by polarization information in model identifiability, Section IV-B studies in detail the $P=2$ sources case and provides a sufficient condition for QNMF uniqueness. Section IV-C considers the general $P\geq 2$ case and gives a necessary uniqueness condition that relaxes the requirement for sources $\mathbf{W}$ to have some zero entries.

IV-A Trivial ambiguities and NMF-based uniqueness

Let $\mathbf{X}\in\mathbb{H}_{\mathcal{S}}^{M\times N}$ and suppose that there exists two matrices $\mathbf{W}\in\mathbb{H}_{\mathcal{S}}^{M\times P}$ and $\mathbf{H}\in\mathbb{R}_{+}^{P\times N}$ such that $\mathbf{X}=\mathbf{W}\mathbf{H}$ holds. Given a nonsingular matrix $\mathbf{T}\in{\mathbb{H}}^{P\times P}$, any pair $(\widetilde{\mathbf{W}},\widetilde{\mathbf{H}})$ defined as

[TABLE]

leaves the data matrix unchanged i.e. $\mathbf{X}=\mathbf{W}\mathbf{H}=\widetilde{\mathbf{W}}\widetilde{\mathbf{H}}$. However, to be admissible, the linear transformation $\mathbf{T}$ should yield matrices $\widetilde{\mathbf{W}}$ and $\widetilde{\mathbf{H}}$ such that

[TABLE]

The real-valuedness constraint on $\widetilde{\mathbf{H}}$ directly imposes that $\mathbf{T}$ is a real matrix. Indeed, if the entries $\mathbf{T}$ are quaternion-valued, so are the entries of $\mathbf{T}^{-1}$, and thus $\widetilde{\mathbf{H}}=\mathbf{T}^{-1}\mathbf{H}$ is also quaternion-valued. As a result, only linear transformations $\mathbf{T}\in\mathbb{R}^{P\times P}$ may yield alternative factors $\widetilde{\mathbf{W}}$ and $\widetilde{\mathbf{H}}$ satisfying (13), just as in the standard NMF case. As a result, QNMF exhibits the same trivial ambiguities as NMF, summarized by Proposition 1 below.

Proposition 1**.**

The QNMF $\mathbf{X}=\mathbf{W}\mathbf{H}$ exhibits two intrinsic ambiguities, namely,

• •

scale indeterminacy:

[TABLE]

• •

order indeterminacy:

[TABLE]

for which $\widetilde{\mathbf{W}}=\mathbf{W}\mathbf{T}$ and $\widetilde{\mathbf{H}}=\mathbf{T}^{-1}\mathbf{H}$ define equally valid QNMF factors.

These inevitable indeterminacies motivate the following definition for the identifiability of the QNMF problem.

Definition 1**.**

The QNMF $\mathbf{X}=\mathbf{W}\mathbf{H}$ is said to be essentially unique if $\mathbf{X}=\widetilde{\mathbf{W}}\widetilde{\mathbf{H}}$ implies $\widetilde{\mathbf{W}}=\mathbf{W}\mathbf{D}\mathbf{P}$ and $\widetilde{\mathbf{H}}=\left(\mathbf{D}\mathbf{P}\right)^{-1}\mathbf{H}$ where $\mathbf{D}$ is a diagonal matrix with strictly positive entries and $\mathbf{P}$ is a permutation matrix.

QNMF features the same trivial ambiguities as NMF: those can be handled by standard techniques, e.g. by imposing normalization on the source coefficients to get rid of the scale indeterminacy.

As explained in Section III-B, QNMF generalizes NMF to the case of polarized signals. The co-factorization perspective given by (10) shows that QNMF inherits sufficient uniqueness conditions from NMF.

Proposition 2** (NMF-based sufficient condition for uniqueness).**

If the NMF $\mathrm{Re}\,\mathbf{X}=\left[\mathrm{Re}\,\mathbf{W}\right]\mathbf{H}$ is essentially unique, then the QNMF $\mathbf{X}=\mathbf{W}\mathbf{H}$ is essentially unique.

Proof.

Suppose that $\mathrm{Re}\,\mathbf{X}=\left[\mathrm{Re}\,\mathbf{W}\right]\mathbf{H}$ is essentially unique. Thus $\mathbf{H}$ can always be written as $\mathbf{H}=\mathbf{D}\mathbf{P}\widetilde{\mathbf{H}}$ where $\mathbf{D}$ is a diagonal matrix with strictly positive entries and $\mathbf{P}$ is a permutation matrix. Then $\mathbf{X}=\mathbf{W}\mathbf{H}=\widetilde{\mathbf{W}}\widetilde{\mathbf{H}}$, with $\widetilde{\mathbf{W}}=\mathbf{W}\mathbf{D}\mathbf{P}$ and $\widetilde{\mathbf{H}}=\left(\mathbf{D}\mathbf{P}\right)^{-1}\mathbf{H}$ which shows that the QNMF is essentially unique by Definition 1. ∎

Proposition 2 shows that known sufficient conditions for NMF directly extend to the QNMF case. For instance, if $\mathrm{Re}\,\mathbf{W}$ and $\mathbf{H}$ satisfy the separability related assumptions of [43, 42] or the more relaxed sufficiently scattered condition [44], then the QNMF $\mathbf{X}=\mathbf{W}\mathbf{H}$ is essentially unique. While this result is reassuring – QNMF encompasses known NMF uniqueness conditions – it does not take advantage of polarization information contained in $\mathrm{Im}\,\mathbf{W}$. The two remaining sections thus explore how this supplementary information together with the polarization constraint ($\mathcal{S}$) impact QNMF model identifiability.

IV-B Two sources case: range of admissible solutions and sufficient uniqueness condition

We first consider the simpler case of two sources ($P=2$), which is particularly instrumental in understanding the role played by the polarization constraint ($\mathcal{S}$). The approach follows closely the one presented in [41] for the NMF case.

For $1\leq m\leq M$ and $1\leq n\leq N$, the corresponding entries of each of the factors $\mathbf{W}$ and $\mathbf{H}$ are given by

[TABLE]

where $p=1,2$ denotes the source index. We use the convenient quaternion form (3) for $w_{mp}$, i.e.

[TABLE]

To model QNMF indeterminacies, we parameterize the matrix $\mathbf{T}$ in Eqs. (11) – (12) as $\mathbf{T}=\mathbf{T}(\alpha,\beta)$, $\alpha,\beta\in\mathbb{R}$ such that

[TABLE]

This form for $\mathbf{T}(\alpha,\beta)$, together with the assumption $\alpha+\beta<1$ ensures that $\mathbf{T}(\alpha,\beta)$ models only non-trivial ambiguities of QNMF, i.e. other than those given in Proposition 1. It follows that $\mathbf{T}$ is invertible, with inverse

[TABLE]

Now consider the linear transformed factors $\widetilde{\mathbf{W}}=\mathbf{W}\mathbf{T}(\alpha,\beta)$ and $\widetilde{\mathbf{H}}=\mathbf{T}^{-1}(\alpha,\beta)\mathbf{H}$. Using (18), the mixing of sources explicitly reads for $1\leq m\leq M$

[TABLE]

With (19) one gets the new activations for $1\leq n\leq N$

[TABLE]

Non-negativity of $\widetilde{\mathbf{H}}$ imposes that $\widetilde{h}_{1n}\geq 0$ and $\widetilde{h}_{2n}\geq 0$ for every $n$. The polarization constraints ($\mathcal{S}$) impose that $\widetilde{w}_{m1},\widetilde{w}_{m2}\in\mathbb{H}_{\mathcal{S}}$. In other words, for every $m$, one has

[TABLE]

Eq. (24) corresponds also to the usual NMF case [41]. However, condition (25) reveals the specifity of QNMF.

In order to find the range of admissible solutions for QNMF, we search for intervals $\mathcal{I}^{\mathrm{QNMF}}_{\alpha}$ and $\mathcal{I}^{\mathrm{QNMF}}_{\beta}$ such that, for every $(\alpha,\beta)\in\mathcal{I}^{\mathrm{QNMF}}_{\alpha}\times\mathcal{I}^{\mathrm{QNMF}}_{\beta}$ the linear transform $\mathbf{T}(\alpha,\beta)$ yields valid $\widetilde{\mathbf{W}}$ and $\widetilde{\mathbf{H}}$ factors. Remark that, since Eqs. (22) – (24) are identical to the NMF case, we can express the intervals $\mathcal{I}^{\mathrm{QNMF}}_{\alpha}$ and $\mathcal{I}^{\mathrm{QNMF}}_{\beta}$ as

[TABLE]

where $\mathcal{I}^{\mathrm{NMF}}_{\alpha}$ (resp. $\mathcal{I}^{\mathrm{NMF}}_{\beta}$) is the interval obtained for $\alpha$ (resp. $\beta$) using non-negativity constraints only, as explained below. For every $m$, $\mathcal{I}^{\mathrm{pol}}_{\alpha,m}$ and $\mathcal{I}^{\mathrm{pol}}_{\beta,m}$ indicate intervals for $\alpha$ and $\beta$ such that the polarization constraint (25) holds. Thanks to this rewriting, Eqs. (26)–(27) permit to separate the contributions of the non-negativity and polarization constraints, respectively. This allows to precisely quantify the role played by polarization constraint in QNMF to improve NMF identifiability.

Computation of $\mathcal{I}^{\mathrm{NMF}}_{\alpha}$ and $\mathcal{I}^{\mathrm{NMF}}_{\beta}$. These intervals encode 2 non-negativity conditions: the real part of the sources and the activation coefficients. Their expression follow directly from NMF results [41], that is:

[TABLE]

Lower bounds depend on the ratio of the real parts of the sources like

[TABLE]

where for convenience $\mathrm{Re}\,w_{mp}=I_{mp}$ as in (17) and $\mathbb{M}_{1}=\{m|I_{m2}>I_{m1}\}$, $\mathbb{M}_{2}=\{m|I_{m1}>I_{m2}\}$ are supposed to be non-empty sets. If $\mathbb{M}_{1}=\emptyset$ (resp. $\mathbb{M}_{2}=\emptyset$), one has $\alpha_{\text{min}}^{\text{NMF}}=-\infty$ (resp. $\beta_{\text{min}}^{\text{NMF}}=-\infty$). Activations coefficients control the upper bounds of $\mathcal{I}^{\mathrm{NMF}}_{\alpha}$ and $\mathcal{I}^{\mathrm{NMF}}_{\beta}$:

[TABLE]

Computation of $\mathcal{I}^{\mathrm{pol}}_{\alpha}$ and $\mathcal{I}^{\mathrm{pol}}_{\beta}$. Let us now exploit the polarization condition (25), on which relies QNMF. Fix $m$ $(1\leq m\leq M)$ and start with parameter $\alpha$. Plugging (17) into (20) ones gets:

[TABLE]

By developing all terms appearing in the polarization constraint (25) and reorganizing the expression as a second-order polynomial in $\alpha$, we get the following inequality

[TABLE]

where $\left\langle\bm{\mu}_{1},\bm{\mu}_{2}\right\rangle=-\mathrm{Re}\,\bm{\mu}_{1}\bm{\mu}_{2}$ can be identified with the usual inner product of $\mathbb{R}^{3}$. Starting from (21), the same approach is used for $\beta$, leading to

[TABLE]

For a given $m$, intervals $\mathcal{I}^{\mathrm{pol}}_{\alpha,m}$ and $\mathcal{I}^{\mathrm{pol}}_{\beta,m}$ are obtained by solving the corresponding polynomial inequalities (33) and (34). It involves finding the roots of the associated second-order polynomials in $\alpha$ and $\beta$, respectively. An easy but lengthy computation of polynomial discriminants shows that (33) and (34) always admit two real-valued solutions (possibly degenerate). Unfortunately, the amount of parameters involved (at least 4: $\Phi_{m1},\Phi_{m2},\left\langle\bm{\mu}_{m1},\bm{\mu}_{m2}\right\rangle$ and the ratio $I_{m1}/I_{m2}$) prevents from performing a general theoretical study of the roots behaviour. Such a study would require to account for numerous particular cases, notably due to the fact that the sign of the second-order term is not constant and can even cancel out for some values of parameters. Nonetheless, in practice when values of sources parameters are given, Eqs. (33) and (34) can be solved very efficiently numerically to yield desired intervals $\mathcal{I}^{\mathrm{pol}}_{\alpha,m}$ and $\mathcal{I}^{\mathrm{pol}}_{\beta,m}$. See Section VI for numerical illustrations on synthetic and real-world data.

Sufficient uniqueness condition. Despite the apparent complexity of the polarization conditions (33) and (34), a simple and interpretable sufficient condition for the essential uniqueness of QNMF with $P=2$ can be formulated.

Proposition 3**.**

If the following conditions are satisfied:

• •

$\exists\>m_{1},m_{2}\in\{1,2,\ldots M\}$* s.t.*

[TABLE]

• •

$\exists\>n_{1},n_{2}\in\{1,2,\ldots N\},n_{1}\neq n_{2}$* s.t.*

[TABLE]

then the QNMF $\mathbf{X}=\mathbf{W}\mathbf{H}$ is essentially unique.

Proof.

See Appendix A-A. ∎

On the one hand, condition (C2) is identical to the one found for the activation factor in the standard NMF case [41, 43, 44]. On the other hand, condition (C1) illustrates the key role played by polarization information. Compared to the usual NMF sufficient conditions for the 2 sources case [41, 43, 44], it does not require each source to vanish alternatively. Condition (C1) shows that it is sufficient that there exists two indices $m_{1},m_{2}$ such that the first source is fully polarized at $m_{1}$ (and the second one exhibits a different polarization state at $m_{1}$) and the second one is fully polarized at $m_{2}$ (and the first one exhibits a different polarization state at $m_{2}$), and that respective intensities of each source should be larger than the other one by a factor in $[0,1]$.

Note that (C1) does not require at all $m_{1}\neq m_{2}$. In fact, when $m_{1}=m_{2}=m$, it becomes condition (CS1’)

[TABLE]

In other words, if there exist $m$ such that both sources exhibit some arbitrary energy and are fully polarized with different polarization axes and if (C2) is satisfied, then the QNMF $\mathbf{X}=\mathbf{W}\mathbf{H}$ is essentially unique.

The sufficient condition given in Proposition 3 for the uniqueness of QNMF in the two source case appears remarkably broad compared to the NMF case [41]. By incorporating polarization information and its associated constraints, QNMF permits to achieve uniqueness even when the sources never vanish – a typical case where NMF is known to be not unique [44, 41, 48].

IV-C General case $(P\geq 2)$

The study of the uniqueness of QNMF $\mathbf{X}=\mathbf{W}\mathbf{H}$ for an arbitrary number of $P$ sources is much more cumbersome than the $P=2$ sources case, as the previous explicit modeling of indeterminacies is no longer amenable. However, the relationship between QNMF and NMF (see Section III-B) allows to formulate a necessary condition for uniqueness where standard NMF would fail.

Uniqueness of the NMF $\mathrm{Re}\,\mathbf{X}=\left[\mathrm{Re}\,\mathbf{W}\right]\mathbf{H}$ requires each column of $\mathrm{Re}\,\mathbf{W}$ to contain at least one entry equal to zero222Note that requiring some entries of $\mathrm{Re}\,\mathbf{W}$ to be zero is equivalent to require the corresponding entries of $\mathbf{W}$ to be zero due to the nature of ($\mathcal{S}$). [43, 41, 44]. Unfortunately, this condition is often violated, notably in hyperspectral imaging applications, and thus much effort has been devoted to the design of suitable optimization criteria [49, 45, 48, 46] to recover uniqueness of the (algorithmic) solution.

In contrast, QNMF relaxes the NMF zero-entry necessary condition on $\mathbf{W}$, stated by Proposition 4 below.

Proposition 4** (Necessary condition for uniqueness with non-vanishing sources).**

Suppose that the QNMF $\mathbf{X}=\mathbf{W}\mathbf{H}$ is essentially unique such that $\mathrm{Re}\,w_{mp}>0$ for every $m,p$. Then the following conditions are satisfied:

• •

$\forall(p,q),\>p\neq q,$**

[TABLE]

• •

$\forall(p,q),\>p\neq q,$**

[TABLE]

Proof.

See Appendix A-B. ∎

Proposition 4 shows that a necessary condition for uniqueness of the QNMF with non-vanishing sources is that, for any distinct pair of sources $(p,q)$, two criteria are satisfied: (A1) there exists an index $m$ such that the source $p$ is fully polarized and the source $q$ exhibits a different polarization state (it can be fully polarized but with a different polarization axis); (A2) there is an index $n$ such that the source $q$ is active while the source $p$ is not. Finally, note that condition (A2) on the activation factor is exactly the same as in the NMF case [41, 43, 44].

In our opinion, the result stated by Proposition 4 is particularly interesting: it suggests that, in the general case of $P\geq 2$ sources, there may exist sufficient uniqueness conditions for QNMF that do not require at all the source factor $\mathbf{W}$ to contains zero-entries, unlike NMF. This conjecture is motivated by Proposition 3, which provides such a sufficient condition for the $P=2$ sources case. It is also supported by our numerical experiments in Section VI-B, where QNMF uniqueness is observed on a synthetic, fully polarized spectro-polarimetric dataset with $P=3$ sources. QNMF might offer a convenient framework to deal with situations where NMF alone does not guarantee uniqueness, such as with non-vanishing sources in hyperspectral imaging. This motivates us to search for general sufficient uniqueness conditions in future work.

V An algorithm for QNMF

This section deals with the practical resolution of the QNMF problem. For a given value $P$, its resolution can be seen as an optimization procedure,

[TABLE]

where $D:\mathbb{H}^{M\times N}\times\mathbb{H}^{M\times N}\rightarrow\mathbb{R}_{+}$ can be an arbitary cost function for quaternion matrices. For simplicity, we choose here the Frobenius distance between quaternion matrices, referred to as the Euclidean cost in the sequel:

[TABLE]

The generic QNMF optimization problem (35) becomes

[TABLE]

The formulation of the resolution of the QNMF problem (35) appears much alike the standard NMF with Euclidean cost. However, two fundamental questions need to be answered: (i) is the constraint ($\mathcal{S}$) easy to implement? and (ii) can we optimize w.r.t. $\mathbf{W}$ directly in the quaternion domain?

Fortunately, the answers to these two questions are affirmative. For (i), the answer relies on the key link between ($\mathcal{S}$) and the set of non-negative Hermitian 2-by-2 matrices. A positive answer to (ii) is made possible by the recent advent of the theory of quaternion-domain derivatives [50, 51, 52, 53] – the so-called generalized $\mathbb{H}\mathbb{R}$ calculus.

Section V-A below explains how constraints on $\mathbf{W}$ and $\mathbf{H}$ factors are implemented. Section V-B describes the proposed alternating least squares strategy to solve (37). For completeness, Appendix B-A provides a quick introduction to quaternion optimization. Detailed computation of QNMF updates are found in Appendices B-C and B-B.

V-A Projections onto constraints

Projections onto the constraint sets of $\mathbf{W}$ and $\mathbf{H}$ are a cornerstone of any method attempting to solve numerically the QNMF factorization problem (35) or (37). Projection operators onto constraints $\mathbb{R}_{+}$ and ($\mathcal{S}$) are denoted by $\Pi_{\mathbb{R}_{+}}$ and $\Pi_{\mathbb{H}_{\mathcal{S}}}$, respectively. The projection $\Pi_{\mathbb{R}_{+}}$ is classical [54, 55]:

[TABLE]

where $\mathbf{M}\in\mathbb{R}^{P\times N}$, $1\leq p\leq P$, $1\leq n\leq N$.

To compute $\Pi_{\mathbb{H}_{\mathcal{S}}}$, we use the fact that the constraint ($\mathcal{S}$) can be derived from non-negativity (i.e. semi-definiteness) of 2-by-2 complex Hermitian matrices – as described in Section II-D. Consider an arbitrary matrix $\mathbf{M}\in\mathbb{H}^{M\times P}$ with $(m,p)$ coefficient $(\mathbf{M})_{mp}=M_{mp}=\mathrm{Re}\,M_{mp}+\bm{i}\mathrm{Im}_{\bm{i}}\,M_{mp}+\bm{j}\mathrm{Im}_{\bm{j}}\,M_{mp}+\bm{k}\mathrm{Im}_{\bm{k}}\,M_{mp}$. Then define the following bijective mapping $f:\mathbb{H}\rightarrow\mathbb{C}^{2\times 2}$ such that

[TABLE]

By construction, $f\left(M_{mp}\right)$ is an Hermitian matrix. Its projection onto the set of non-negative Hermitian matrices $\mathcal{K}_{+}^{2}$ is well-known [56, p. 399], i.e.

[TABLE]

where $\eta_{i}$ and $\mathbf{v}_{i}$ are the $i$-th eigenvalue and eigenvector of the matrix $f\left(M_{mp}\right)$, respectively. Since $\Pi_{\mathcal{K}_{+}^{2}}\left[f\left(M_{mp}\right)\right]$ is a 2-by-2 non-negative Hermitian matrix, it can be expressed as

[TABLE]

where $a\geq 0$, $b\geq 0$ and $c\in\mathbb{C}$ such that $ab-|c|^{2}\geq 0$ on account of positive semi-definiteness of $\Pi_{\mathcal{K}_{+}^{2}}\left[f\left(M_{mp}\right)\right]$. Remark that (41) is of the form (39): since the mapping $f$ is bijective, the projection $\left[\Pi_{\mathbb{H}_{\mathcal{S}}}(\mathbf{M})\right]_{mp}$ is uniquely given by

[TABLE]

where $a,b$ and $c$ as defined above.

V-B Proposed algorithm: quaternion alternating least squares

The proposed algorithm adopts a popular strategy for solving the QNMF problem, based on the alternating constrained minimization of (37) w.r.t. $\mathbf{H}$ and $\mathbf{W}$. This choice is motivated by the fact that, whereas the Euclidean cost (37) is not convex in both $\mathbf{W}$ and $\mathbf{H}$, it is convex in each variable separately. After initialization of the factors, the iteration $r>0$ reads

[TABLE]

Equations (43)-(44) describe a two-block coordinate descent (CD) scheme in the quaternion domain. Two-block real-domain CD schemes are known to converge to a stationary point under mild conditions [57, 58, 59]. These convergence results can be readily transposed to the quaternion case by identifying $\mathbb{H}$ with $\mathbb{R}^{4}$. Formally, since $\mathbb{R}^{+}$ and $\mathbb{H}_{\mathcal{S}}$ are closed convex sets, every limit point of the sequence $(\mathbf{W}_{r},\mathbf{H}_{r})$ generated by the two-block CD framework (43)-(44) is a stationary point of (37).

Convergence guarantees of the two-block CD scheme require to solve exactly both subproblems (43) and (44). However, the associated constraints make it difficult to obtain these updates directly in a closed form. The quaternion nature of the data $\mathbf{X}$ and source factor $\mathbf{W}$ adds a supplementary degree of difficulty due to non-commutativity of quaternion multiplication. Solving exactly each subproblem (43)-(44) requires to properly extend constrained optimization algorithms over the set of quaternions, which is beyond the scope of this paper.

Instead, we propose a simpler, closed-form algorithm that approximately solves each subproblem at each iteration, at the price of losing convergence guarantees. The main idea is the following: at a given iteration $r$, for each factor, one solves the unconstrained least-squares problem and projects the obtained solution onto the corresponding constraint:

[TABLE]

Due to its resemblance with the usual alternating least squares (ALS) algorithm [60] for NMF, we call this strategy quaternion alternating least squares (QALS). Derivation of explicit updates for unconstrained least-squares problems requires special care because of the quaternion nature of $\mathbf{X}$ and $\mathbf{W}$. It involves the notion of quaternion derivative as well as cautious handling of quaternion non-commutivity, see detailed computations provided in Appendices B-B and Appendix B-C. As a result, one gets the explicit updates:

[TABLE]

Projections onto constraints sets are carried out as described in Section V-A.

The proposed algorithm is remarkably simple and cheap. Performing optimization directly in the quaternion-domain yields updates expressions that are very much alike the standard ALS algorithm for NMF [61]. However, this apparent simplicity should not conceal the important underlying technical details presented in Appendix B. Such computations would not have been amenable without using the powerful theory of quaternion derivatives introduced recently [50, 51, 52, 53]. Despite the lack of theoretical guarantees on its convergence, the quaternion ALS algorithm provides a good baseline for the resolution of QNMF with reasonably good results in most situations. Thus, it paves the way to further developments of more sophisticated algorithms for the QNMF problem.

VI Numerical experiments

This section illustrates the relevance of the proposed approach by performing numerical experiments on synthetic data. Section VI-A provides numerical evidence to the key role played by polarization properties onto the range of admissible solutions for the 2 sources case. Section VI-B demonstrates the effectiveness of the QALS algorithm to solve the QNMF problem.

VI-A Polarization information and identifiability

This section illustrates the theoretical study provided in Section IV-B regarding the identifiability of QNMF in the $P=2$ sources case. To emphasize how polarization impacts the range of admissible solutions, we choose a setting where the QNMF $\mathbf{X}=\mathbf{W}\mathbf{H}$ is not unique – and thus, neither is the associated NMF $\mathrm{Re}\,\mathbf{X}=\mathrm{Re}\,\left[\mathbf{W}\right]\mathbf{H}$, see (10).

Both sources and associated activations are non-vanishing, i.e. for all $m,n$, one has $w_{m1},w_{m2},h_{1n},h_{2n}\neq 0$. Intensity responses $\mathrm{Re}\,w_{m1},\mathrm{Re}\,w_{m2}$ encode realistic $S_{0}$ parameters, obtained for instance in hyperspectral imaging of wood components [62]. To simplify, both sources are chosen to exhibit constant polarization properties:

[TABLE]

Each source being partially polarized, it precludes the QNMF to be unique, since the necessary conditions stated by Proposition 4 are not fulfilled.

Numerical computations yield the range of admissible solutions for the QNMF and its associated NMF problem. Following Section IV-B, the range of admissible parameters $\alpha$ and $\beta$ defining the transformation matrix $\mathbf{T}(\alpha,\beta)$ are found for the QNMF

[TABLE]

and for the associated NMF problem:

[TABLE]

Note that the QNMF and NMF intervals share the same upper bounds on $\alpha$ and $\beta$, which arise from the activation factor ratios (30).

Fig. 1 represents the set of admissible solutions corresponding to these NMF and QNMF intervals, respectively. Thick lines indicate the original sources and activations factors. Comparing respective ranges of solutions for $S_{0}$, it appears that QNMF significantly improves identifiability over standard NMF, by taking advantage of polarization information. Improvements are found on both $S_{0}$ and activations factors, and importantly, QNMF permits to reconstruct Stokes parameters $S_{1},S_{2},S_{3}$ with limited uncertainty. These results illustrate how QNMF takes advantage of the strong discriminative power of polarization in order to improve model identifiability.

VI-B Spectro-polarimetric blind unmixing using QALS

This section provides a first numerical validation of the QALS algorithm on realistic simulated spectro-polarimetric data. It consists in a linear mixture of $P=3$ sources imaged at $M=128$ wavelengths and $N=512\times 512$ spatial locations. Each source exhibits broad-band behavior, featuring non-vanishing intensity $S_{0}$ spectral profile together with spectrally varying polarization $S_{1}/S_{0},S_{2}/S_{0},S_{3}/S_{0}$ parameters. For simplicity, we assume that each source is fully polarized. Two-dimensional activations maps of each source are designed to represent a typical spectro-polarimetric imaging setting, e.g. solar spectro-polarimetry [63, 7], and such that the resulting activation matrix $\mathbf{H}$ satisfies the separability condition [44]. As a result, the designed $\mathbf{W}$ and $\mathbf{H}$ matrices satisfy the necessary conditions for QNMF identifiability stated in Proposition 4. Fig. 2A shows simulated spectro-polarimetric data $\mathbf{X}=\mathbf{W}\mathbf{H}$ for two distinct wavelengths indices $m=38$ and $m=76$.

We show that QALS permits to effectively recover the QNMF factors. We adopt a random initialization strategy. The initial activation $\mathbf{H}_{0}$ is chosen as a matrix of i.i.d. entries drawn from $\mathcal{U}([0,1])$, the uniform distribution on $[0,1]$. The initial source factor $\mathbf{W}_{0}$ has i.i.d. entries drawn from the quaternion circular unit Gaussian distribution [64] projected onto the constraint ($\mathcal{S}$). We then alternate Eqs. (47)–(48) until convergence. The latter was assessed by monitoring the relative error $\varepsilon_{r}=\|\mathbf{X}-\mathbf{W}_{r}\mathbf{H}_{r}\|^{2}_{F}/\|\mathbf{X}\|^{2}_{F}$ and stopping whenever the improvement was below a given threshold, i.e. $|\varepsilon_{r}-\varepsilon_{r-1}|\leq 10^{-5}$ in our case. Repeating this procedure for $K=100$ independent initializations of the QALS algorithm, it took on average 29 iterations for the algorithm to converge. Out of $100$ initializations, we observed that the QALS algorithm recovered each time identical sources and activation factors, strongly indicating that the QNMF of the data $\mathbf{X}$ may be essentially unique.

Fig. 2B depicts the 3 sources and corresponding activations factors recovered by the QALS algorithm for one arbitrary initialization. Excellent match between reconstructed and ground truth factors is indicated by very small relative errors, i.e. $\varepsilon_{\mathbf{H}}=\|\hat{\mathbf{H}}-\mathbf{H}\|^{2}_{F}/\|\mathbf{H}\|^{2}_{F}=4.91\cdot 10^{-5}$ and $\varepsilon_{\mathbf{W}}=\|\hat{\mathbf{W}}-\mathbf{W}\|^{2}_{F}/\|\mathbf{W}\|^{2}_{F}=3.92\cdot 10^{-5}$. Looking at sources intensities $S_{0}$, we observe that they never vanish, i.e. $S_{0}>0$, plus that the first and second sources exhibit very similar peaks. This is a very challenging situation for NMF-based methods, which would require additional regularizations such as minimum volume [46, 47, 45] to ensure uniqueness of the solution. It is thus remarkable that the QNMF model with the QALS algorithm is able to solve this blind source problem in a simple manner by leveraging polarization information. One of the advantages of polarization in source separation problems can be directly seen by inspecting $\mathbf{X}_{76:}$ in Fig. 2A: while it is not possible to distinguish 2 sources using $S_{0}$ only, the distinctive polarization pattern on $S_{2}$ makes it easy to identify source 1 and source 2.

These simulation results demonstrate the ability of the proposed QALS algorithm to effectively solve the QNMF problem. Despite no associated provable convergence guarantees, QALS offers a simple, computationally efficient algorithm and provides a baseline for the development of more advanced algorithms for QNMF.

VII Conclusion

This paper has introduced a new powerful tool called quaternion non-negative matrix factorization (QNMF), which generalizes the well-known concept of non-negative matrix factorization (NMF) to the case of polarized signals. The algebraic representation of Stokes parameters using quaternions, together with the generalization of the non-negativity constraint on Stokes parameters made it possible to formulate and prove uniqueness results that may have been cumbersome to obtain otherwise. We unveiled the relation between NMF and QNMF, showing that NMF sufficient conditions are also sufficient uniqueness conditions for QNMF. Remarkably, we have shown that QNMF improves model identifiability over NMF by leveraging the complete information available in the four Stokes parameters. This key role played by polarization information is stated by Proposition 3 for the two sources case and by Proposition 4 for the general case. Furthermore, taking advantage of recent results in quaternion optimization, we have proposed a simple yet efficient algorithm for solving the QNMF problem in practice. These first results appear very promising and pave the way to future work on both theoretical and methodological aspects of QNMF. In particular, it is an open question whether one can generalize to the QNMF case the notion of non-negative rank [65] and its related properties. This fundamental question is intimately related to the definition of a polarized source; as such, it would greatly improve our understanding of the role played by polarization diversity in blind source separation problems.

Appendix A Uniqueness conditions for QNMF

A-A Proof of Proposition 3

Suppose that $P=2$ and that the QNMF $\mathbf{X}=\mathbf{W}\mathbf{H}$ exists. We suppose that conditions (C1) and (C2) are satisfied. To show that the QNMF is identifiable, we prove that intervals $\mathcal{I}^{\mathrm{QNMF}}_{\alpha}$ and $\mathcal{I}^{\mathrm{QNMF}}_{\beta}$ defined in (26)-(27) are restricted to $\{0\}$.

Condition (C1) permits to simplify (33) and (34) for $m_{1},m_{2}$ like

[TABLE]

Using inequalities linking $I_{m1},I_{m2}$ for $m=m_{1},m_{2}$ in (C1) yields the following domains of solutions to (53) and (54):

[TABLE]

where $\alpha_{0},\beta_{0}\geq 1$. Condition (C2) imply that

[TABLE]

where $\alpha_{0}^{\prime},\beta_{0}^{\prime}\leq 0$. By intersection of intervals, one gets that $\mathcal{I}^{\mathrm{QNMF}}_{\alpha}$ and $\mathcal{I}^{\mathrm{QNMF}}_{\beta}$ defined in (26)-(27) are restricted to $\{0\}$, so that the QNMF is unique.

A-B Proof of Proposition 4

Suppose that the QNMF $\mathbf{X}=\mathbf{W}\mathbf{H}$ exists and that sources never vanish, i.e. $\mathrm{Re}\,w_{mp}>0$ for every $m,p$. We obtain conditions (A1) and (A2) by contradiction: one at a time we suppose that either (A1) or (A2) is false and show that there exists a non-trivial matrix $\mathbf{T}$ leading to different factors $\widetilde{\mathbf{W}}=\mathbf{W}\mathbf{T}$ and $\widetilde{\mathbf{H}}=\mathbf{T}^{-1}\mathbf{H}$.

Suppose that (A1) is not satisfied. Then, $\exists(p_{0},q_{0}),p_{0}\neq q_{0}$ such that

[TABLE]

or equivalently,

[TABLE]

Now consider the transformation $\mathbf{T}_{p_{0}q_{0}}^{\alpha}\in\mathbb{R}^{P\times P}$ defined by

[TABLE]

By construction $\mathbf{T}_{p_{0}q_{0}}^{\alpha}$ does not correspond to a trivial ambiguity of the QNMF when $\alpha\neq 0$. Note that $(\mathbf{T}_{p_{0}q_{0}}^{\alpha})^{-1} =\mathbf{T}_{p_{0}q_{0}}^{-\alpha}$. Consider the new factors $\widetilde{\mathbf{W}}=\mathbf{W}\mathbf{T}_{p_{0}q_{0}}^{\alpha}$ and $\widetilde{\mathbf{H}}=\mathbf{T}_{p_{0}q_{0}}^{-\alpha}\mathbf{H}$ such that $\mathbf{X}=\widetilde{\mathbf{W}}\widetilde{\mathbf{H}}$. From (59), the source indexed by $p_{0}$ is the only one affected by $\mathbf{T}_{p_{0}q_{0}}^{\alpha}$ like

[TABLE]

The corresponding activation coefficients read

[TABLE]

Supposing $\alpha>0$, then $\widetilde{h}_{p_{0}n}\geq 0$ for every $n$ by non-negativity of the matrix $\mathbf{H}$. It remains to find at least one $\alpha>0$ such that $\forall m,\>\widetilde{w}_{mp_{0}}=w_{mp_{0}}-\alpha w_{mq_{0}}\in\mathbb{H}_{\mathcal{S}}$. Imposing $\widetilde{w}_{mp_{0}}\in\mathbb{H}_{\mathcal{S}}$ for every $m$, yields

[TABLE]

where we used notation (17) for convenience. Since by assumption $I_{mp_{0}},I_{q_{0}}>0$, the first condition implies that $\alpha\leq\min_{m}(I_{mp_{0}}/I_{mq_{0}})$. The second condition related to polarization can be rewritten as, for every $m$

[TABLE]

According to (58) two cases may occur, not exclusively from each other. Fix $m$ and suppose that $\Phi_{mp_{0}}\in[0,1)$. Two cases are possible: either $\Phi_{mq_{0}}=1$ or $\Phi_{mq_{0}}\in[0,1)$. Assume that $\Phi_{mq_{0}}=1$, then (63) is true for $\alpha\leq I_{mp_{0}}(1-\Phi_{mp_{0}}^{2})I_{mq_{0}}^{-1}((1-\Phi_{mp_{0}}\left\langle\bm{\mu}_{mp_{0}},\bm{\mu}_{mq_{0}}\right\rangle)^{-1}$. For $\Phi_{mq_{0}}\in[0,1)$, since $I_{mq_{0}}^{2}(1-\Phi_{mq_{0}}^{2})>0$ in virtue of assumptions, the discriminant $\Delta$ reads

[TABLE]

Ranges of parameters imply that $\Delta\geq 0$. As a result, (63) is true for $\alpha\in(-\infty,\alpha_{-}^{m}]\cup[\alpha_{+}^{m},+\infty)$ where $\alpha_{\pm}^{m}$ are the roots of the polynomial. Sign of polynomial coefficients in (63) together with $0\leq\Phi_{mp_{0}}<1$ imply that $\alpha_{\pm}^{m}>0$. Suppose now that $\Phi_{mq_{0}}=1$ and $\bm{\mu}_{mq_{0}}=\bm{\mu}_{mp_{0}}$. Taking $\Phi_{mp_{0}}=1$ Eq. (63) becomes trivial. The case $\Phi_{mq_{0}}\in[0,1)$ is included in the previous discussion.

Summarizing all cases, for every $m$ such that $\Phi_{mp_{0}}\in[0,1)$ or $\Phi_{mq_{0}}\bm{\mu}_{mq_{0}}=\bm{\mu}_{mp_{0}}$, there always exists $\alpha>0$ such that conditions (62) are satisfied, meaning that the QNMF is not unique. This leads to the first necessary condition (A1). Repeating the same approach with $\alpha<0$ yields condition (A2).

Appendix B Derivation of quaternion ALS updates

B-A Quaternion derivatives using generalized $\mathbb{H}\mathbb{R}$-calculus

Given a function $f:\mathbb{H}\rightarrow\mathbb{H}$ of the variable $q\in\mathbb{H}$, one outstanding question is: does its derivative $\partial f/\partial q$ exists and if so, how do we compute it? Quaternion analytic functions [66, 67] are known to be differentiable, unfortunately they form a very restricted class of functions – of little interest to signal processing. In fact, cost functions $f:\mathbb{H}\rightarrow\mathbb{R}$ are not analytic [66], so that other strategies need to be deployed. First, a pseudo-derivative approach can be used by treating $f(q)$ as a function $f(q_{0},q_{1},q_{2},q_{3})$ of the four real components of $q$. As pointed out in [50], such approach requires lengthy and cumbersome computations, thus limiting its applicatibility. In order to compute derivatives directly in the quaternion domain, crucial steps have been made recently with the development of the $\mathbb{H}\mathbb{R}$-calculus [53], and subsequently the generalized $\mathbb{H}\mathbb{R}$-calculus [52]. The latter provides a complete framework including chain rule and product rules, extending naturally the $\mathbb{C}\mathbb{R}$-calculus [68] of complex-valued optimization to the case of quaternion functions.

For the quadratic-type of quaternion functions encountered in this paper, only important results related to quaternion optimization are necessary. These are given in Proposition 5 and Table I below. We refer to [52, 50] for a complete description of the generalized $\mathbb{H}\mathbb{R}$-calculus.

Proposition 5 characterizes stationary points of cost functions of quaternion matrices.

Proposition 5** (Theorem 4.1 [50]).**

Let $f:\mathbb{H}^{M\times N}\rightarrow\mathbb{R}$ and denote by $\nabla_{\mathbf{Q}}f$ its gradient w.r.t. $\mathbf{Q}$. Then,

[TABLE]

For practical computations, Table I reprints from [50] some useful functions of the quaternion matrix variable $\mathbf{Q}$ and their corresponding gradients.

B-B Updates for $\mathbf{H}$

Due to non-commutativity of the quaternion product, standard rules matrix derivatives cannot be applied directly. Instead, let us write explicitly the Euclidean cost function

[TABLE]

where we have used that $\overline{w_{mk}h_{kn}} =\overline{w_{mk}}h_{kn}$ since $h_{kn}$ is real-valued. Using the chain rule, the partial derivative of (67) w.r.t. $h_{ij}$ is

[TABLE]

As a result, one gets

[TABLE]

The Euclidean cost function is convex in $\mathbf{H}$; therefore, the minimizer ${\mathbf{H}}_{\star}$ of this cost is obtained by cancelling out the gradient (73):

[TABLE]

where again we used that $\overline{\mathbf{W}{\mathbf{H}}_{\star}}=\overline{\mathbf{W}}{\mathbf{H}}_{\star}$ since ${\mathbf{H}}_{\star}$ is a real matrix. It also implies that $\mathrm{Re}\,\left[\mathbf{W}^{\top}\overline{\mathbf{W}}{\mathbf{H}}_{\star}\right]=\mathrm{Re}\,\left[\mathbf{W}^{\top}\overline{\mathbf{W}}\right]{\mathbf{H}}_{\star}$, leading to the following expression for ${\mathbf{H}}_{\star}$:

[TABLE]

B-C Updates for $\mathbf{W}$

Let us rewrite the Euclidean distance (37) in terms of a trace operator

[TABLE]

where $\mathbf{Q}^{\dagger}=\overline{\mathbf{Q}}^{\top}$ is the (quaternion) conjugate-transpose of $\mathbf{Q}$. Developping333We use that, for any quaternion matrices $\mathbf{A}$ and $\mathbf{B}$, $(\mathbf{A}\mathbf{B})^{\dagger}=\mathbf{B}^{\dagger}\mathbf{A}^{\dagger}$. See e.g. [27] for more complete properties of quaternion matrices. (77) one gets

[TABLE]

By Proposition 5, stationary points are defined by cancelling values of the gradient of (78) w.r.t. to $\overline{\mathbf{W}}$. Thanks to Table I one gets

[TABLE]

The Euclidean cost function is convex in $\mathbf{W}$. Thus, cancelling out the gradient (81) yields the global minimizer ${\mathbf{W}}_{\star}$ of this cost such that

[TABLE]

so that

[TABLE]