Separable Quaternion Matrix Factorization for Polarization Images

Polarization is a unique characteristic of transverse wave and is represented by Stokes parameters. Analysis of polarization states can reveal valuable information about the sources. In this paper, we propose a separable low-rank quaternion linear mixing model to polarized signals: we assume each column of the source factor matrix equals a column of polarized data matrix and refer to the corresponding problem as separable quaternion matrix factorization (SQMF). We discuss some properties of the matrix that can be decomposed by SQMF. To determine the source factor matrix in quaternion space, we propose a heuristic algorithm called quaternion successive projection algorithm (QSPA) inspired by the successive projection algorithm. To guarantee the effectiveness of QSPA, a new normalization operator is proposed for the quaternion matrix. We use a block coordinate descent algorithm to compute nonnegative factor activation matrix in real number space. We test our method on the applications of polarization image representation and spectro-polarimetric imaging unmixing to verify its effectiveness.


Introduction
Polarization is one of the primary characteristics of transverse waves such as optics [7], seismology [1], radio [42], and microwaves [23]. It is a unique transverse wave phenomenon that describes the geometrical orientation of the oscillations. Different materials, the Earth's surface, and atmospheric targets in the process of reflection, scattering, and transmission and emission of electromagnetic radiation will produce their own nature of the polarization characteristics. Analysis of the polarization states imposed by these processes will potentially reveal useful information about the source. Hence it works well in many applications in remote sensing [27,40], agriculture [33], astronomy [21,24] and many other fields.
For example, hyperspectral imaging is a well-known technique collects and processes spectra intensity information from across the electromagnetic spectrum. But due to the limitations of optical resolution, ground structure, noise and other factors, different materials often exhibit "the same spectra" that are difficult to identify. To address this issue, the polarization modulation technique integrated into spectral imaging that can further lead to a more powerful technique, namely spectro-polarimetric imaging [6,31,32,43]. The spectro-polarimetric imaging data consists of spatial, frequency and polarization information. The emerging technique has been widely used in astronomy [2,38], climatology [22], medical image [8], among others.
The state of polarization can be entirely characterized by the four components of the Stokes vector [4,30]. The Stokes parameters were first defined by George Gabriel Stokes in 1852 and used widely to describe the intensity and polarization of light. They can be measured in experiments [18,39] and reformulated in quaternion algebra [25]. In the blind polarized source separation problems, low-rank approximation models based on quaternions play an important role and have been studied in many references, see for example [20,26,41]. Most of the research works rely on standard quaternion mixing model assumption, without considering physical properties of light. Until 2020, a quaternion linear mixing model named quaternion nonnegative matrix factorization(QNMF) has been proposed in reference [10], specifically for blind unmixing of spectro-polarimetric data. The QNMF model extends nonnegativity constraints based on physical properties of light, and a quaternion alternating least squares (QALS) algorithm is presented accordingly.
Since QNMF is extended from standard NMF, its solution could not be guaranteed to be unique in most cases. However, the uniqueness of a factorization is crucial in practice. For the NMF problem, additional constraints on the factor matrix are usually imposed to guarantee its uniqueness. Inspired by separable NMF, in this paper, we introduce separability into QNMF model, assuming that each column of source factor matrix equals a column of the polarized data matrix. We refer to it as separable quaternion matrix factorization (SQMF).
The separability assumption makes sense in practical situations. For example [28], in hyperspectral unmixing, each column of the data matrix is the spectral signature of a pixel. Separability, also known as pure-pixel assumption, means that for each material (i.e., endmember) within the hyperspectral image, there exists a pixel that only contains that material. For more separable NMF models, we refer the interested reader to [14,34,35].
Similarly, the SQMF model requires in the application of spectro-polarimetric unmixing, that each column of the data matrix consists of the spectral and polarization signature of a pixel. And for any pixels, data can be represented as a linearly weighted combination of several elementary sources. For each source corresponding to a material, there exists a pixel that only contains that source (material). Here the separability in spectro-polarimetric imaging can be seen as a natural extension of pure-pixel assumption in hyperspectral imaging.

Contributions and outline of this paper
This paper considers the quaternion matrix factorization under separability conditions, referred to as separable quaternion matrix factorization (SQMF). The contributions are summarized as follows.
1) To investigate the SQMF problem, we present several properties of SQ-matrices that can be decomposed by SQMF. We show the relationship between SQ-matirces and its corresponding intensity (real) component and polarization (imaginary) component. The uniqueness of SQMF has been proven from a linear algebra perspective.
2) In the SQMF problem, to find the source matrix in quaternion space, we propose an algorithm inspired by the successive projection algorithm (SPA) [3,17] which we refer to the quaternion successive projection algorithm (QSPA). It only requires O(mnr) operations. We demonstrate that QSPA can identify the quaternion-valued factor matrix correctly for SQ-matrices. To compute the activation factor matrix in real number space, we provide a simple algorithm extended from hierarchical nonnegative least squares.
3) We test the proposed method in the polarization image representation and spectro-polarized imaging unmixing. Compared to state-of-the-art techniques, SQMF can capture the main characteristic of polarized data and attain a high approximation to the data. Its computational time is very competitive.
The rest paper is organized as follows. In Section 2, we briefly review quaternion, polarization and Stokes parameters. Section 3 introduces the separable quaternion matrix factorization problem and presents some related properties. We also discuss the uniqueness of the SQMF problem. Section 4 proposes the quaternion successive projection algorithm to determine column subset to form the quaternion factor matrix, and quaternion hierarchical nonnegative least squares to compute the real factor matrix. Section 5 validates the effectiveness of the proposed method in numerical experiments on realistic polarization image and simulated spectro-polarimetric data set. The concluding remarks are given in Section 6.

Preliminaries
Let us start with notations used throughout this paper. The real number field, the complex number field and the quaternion algebra are defined by R, C and H respectively. Unless otherwise specified, lowercase letters represent real numbers, for example, a ∈ R. The bold lowercase letters represent real vectors, such as, a ∈ R n . Real matrices are denoted by bold capital letters, like A ∈ R m×n . The numbers, vectors, and matrices under the quaternion field are represented by the corresponding symbols wearing a hat, for exampleȃ ∈ H,ȃ ∈ H n andȂ ∈ H m×n .
We define the inner product of two quaternion vectors as follows. Givenȃ ∈ H m ,b ∈ H m , Their inner product is defined as, And the l 2 norm of any quaternion vectorq ∈ H m is given below, We use Mat(·) to represent all the components of a quaternion vectorq ∈ H m or matrixQ ∈ H m×n , that is It is easy to verify that In the following, we present an inequality based on the l 2 norm of quaternion vector which is useful in the paper. Lemma 1. For any quaternion matrixQ . = [q 1 , · · · ,q n ] ∈ H m×n , and any real vector

Stokes parameters
Stokes parameters S 0 , S 1 , S 2 , S 3 are used to describe the polarization state of electromagnetic radiation. The first parameter S 0 ≥ 0 represents the total intensity of light, equals sum of the intensity of polarized and un-polarized light. Parameters S 1 , S 2 , S 3 describe the polarization ellipse of light, specifically, S 1 = φS 0 cos 2ψ cos 2χ, φS 0 , 2ψ and 2χ are the spherical coordinates of the three-dimensional vector of cartesian coordinates (S 1 , S 2 , S 3 ). φ is the degree of polarization, the metric of relative contribution of polarized part to the total intensity, φ ∈ [0, 1] implies the polarization state. φ = 1 indicates fully polarized, and φ = 0 means unpolarized. For φ ∈ (0, 1), we say it is partially polarized. These four Stokes parameters (S 0 , S 1 , S 2 , S 3 ) can be expressed by using quaternion algebra as Due to the physical meaning indicated in (3), the Stokes vector (S 0 , S 1 , S 2 , S 3 ) T ∈ R 4 satisfies that Worth noting that, condition (4) is equivalent to the positive semi-definiteness of a 2 × 2 Hermitian matrix J, More precisely, (4) ⇔ trace (J) ≥ 0 and det(J) ≥ 0. The condition (4) implies, for any quaternion q ∈ H satisfies Req ≥ 0 and |Imq| 2 ≤ (Req) 2 . Define the set of constrained quaternions [10] H S ⊂ H such that H S . = {q ∈ H|Req ≥ 0 and |Imq| 2 ≤ (Req) 2 }.
For more details on polarization and Stokes parameters, we refer the interested reader to [12] and the references therein.

Definition and Properties
In this section, we will introduce a linear mixing model for polarized matrixM ∈ H m×n S . It is known that each pixel of a polarization image contains light intensity and its polarization information, To illustrate the motivation of the separable quaternion model, we will use spectro-polarimetric imaging as an example.
As we know, in the application of hyperspectral unmixing, the linear mixing model for hyperspectral imaging (HSI) data assumes that the spectrum of a pixel is a linear weighted combination of the pure spectra (endmembers) of the components present in that pixel. One commonly used assumption in endmember extraction is the presence of pure pixels in the given images, i.e., for each distinct material, at least one pixel (i.e., spectral signature) exists in the image only contains that material. We note that if the imaginary components of all pixels are equal to zero, all lights are unpolarized. The spectro-polarimetric image will degenerate to (HSI) data. From the physical interpretation of light intensity and its polarization, spectro-polarimetric imaging is supposed to share similar material and structural properties as HSI. We consider pure pixel assumption in the linear mixing model for spectro-polarimetric data represented as quaternion form.
Accordingly, we will generalize the standard separable nonnegative model to the quaternion matrix factorization model. Formally, and each column ofW is a column ofM. Precisely,W =M(:, K), where K is column set of cardinality r ofM.
For simplicity, we call a quaternion matrixM the separable quaternion matrix (SQ-matrix) if it has the form of (6). Note that the aim of the separable quaternion factorization is to find the minimal number of source columnsW from the givenM, r here refers to the minimal number of columns inW.
As we know that every entry of the quaternion matrixM is in the constrained quaternion set H S , it is necessary to investigate the plausibility of the assumption in SQMF model (6). The result is shown in the following proposition.
Easy to know that S 0 (m ij ) ≥ 0. To provem ij ∈ H S , in the following, we will show that The results follow.
Proposition 1 implies that any nonnegative weighted linear combination of quaternions in H S also lies in H S , i.e., H S is convex cone. In the following proposition, we will show the connection between a r-separable quaternion matrixM and its contained four components.
and at least one column set K l contains the others, i.e., for some p ∈ {0, 1, 2, 3}, and there exists matrix Q ∈ R r×(n−r) + , such that where Π is appropriate permutation matrix,K . = {1, · · · , n} − K, and I r l is identity matrix of r l . The factor matrices (W, H) : Proof. (=⇒) IfM is quaternion separable matrix, the results are obviously from Definition 1.
(⇐=) Without loss of generality, assume Similarly, We then constructW ∈ H m×r S , and H ∈ R r×n + , such thatM =WH. The results follow.
From Property 1, we remark that quaternion matrixM can not be guaranteed to be r-separable, even if its all four components are r-separable.
The following property shows that quaternion separable matrix factorization admits a unique solution up to permutation and scaling.
ThenW is unique up to permutation and scaling. And if Mat(M) = r, H is also unique up to permutation and scaling .
Proof.M is r quaternion separable matrix, then let From Proposition 1, X = YH and Y = X(:, K). If Y is not unique up to permutation and scaling, there exist (Y 1 , , and the columns of V 2 are from H 2 . Hence we have, Let G = V 1 V 2 , and assume that I r − G = 0. Without loss of generality, let the j-th column of G is non zero column, then Note that g ij = 0 for some i, and Y 1 (:, i) = 0 for all i, Notice that Y 1 (1 : m, :) = S 0 (W) ≥ 0, and G ≥ 0, easy to get 1 − g jj > 0. It implies that Y 1 (:, j) can be represented as linear combination of other columns of Y 1 , with nonnegative coefficients. This is contradict to minimal r-separability ofM. Hence, g jj = 1.
From (11), we get As V 1 and V 2 are nonnegative matrices, V 1 and V 2 can only be permutation matrix or scaling matrix. Therefore, Y is unique up to permutation and scaling.
Specifically, if rank(Mat(M)) = r, that is rank(X) = r, one can easily deduce that Y has full column rank, or else it will contradict to assumption of rank(X) = r. Since X = YH, we get that H is unique up to permutation and scaling from fundamental results of matrix equation.
Remark 1. Reference [10] indicates that if S 0 (M) is r-separable,W is unique up to permutation and scaling. Property 2 extends to the case that S 0 (M) is r 1 -separable when r 1 < r.
We remark that for a quaternion r separable matrixM ∈ H m×n S , in general, H is not unique up to permutations and scaling when r > rank(Mat(M)). In the following, we will give an example to show this fact.

Algorithms
This section proposes quaternion successive projection algorithm extended from SPA, to identify the column set K from a given matrixM. The extracted columns will form the source matrix W. For solving activation matrix H, a quaternion hierarchical nonnegative least squares method is proposed.

Quaternion SPA
Inspired from a well-known separable NMF algorithm named successive projection algorithm (SPA), we design a fast algorithm specifically for separable quaternion matrix factorization. SPA was introduced in [3] regarding spectral unmixing and has been well discussed in [13,29]. Moreover, it is robust in the presence of noise [17]. SPA assumes that the input matrix M ∈ R m×n has similar form as (1), that is M = M(:, K)[I r , H ]Π, where Π is a permutation matrix and H ≥ 0 and ||H (:, t)|| 1 ≤ 1 for all t. The main idea behind SPA which sequentially identifies the columns in K is given as follows: at each step, it first determines the column of M that has the largest 2 norm, and then projects all columns of M onto the orthogonal complement of the extracted column.
To start successive projection algorithm on quaternion field, given a quaternion separable matrix M ∈ H m×n S , there exists an index set K of cardinality r, and some H ∈ R r×n + , such thatM =WH, whereW =M(: Note that the assumption of ||H (:, t)|| 1 ≤ 1 for all t is necessary to guarantee SPA identifies the correct column set K. We remark that it also plays a critical role in the proposed quaternion SPA. From Lemma 1, one can easily verify that the column of largest 2 norm of quaternion separable matrixM will be identified from K if ||u t || 1 ≤ 1 for all t. However in general, ||u t || 1 ≤ 1 is not satisfied. To address this issue, a normalize operator can be applied on each column of SQ-matrix M. How to define a suitable normalize operator for a quaternion matrix will be the key to propose quaternion SPA.
the normalized vector ofq. Apply the normalize operator on SQ-matrixM, for any t, We have, Hence, we define the normalization function f (q) . = S 0 (q) 1 that can satisfy the condition (12). The quaternion successive projection algorithm is presented as follows. k = arg max j z j 2 .

5:
Compute inner product vector p p .
• No parameters need to be tuned.
The correctness of quaternion SPA in the noiseless case is proven in the following theorem in details. Proof. Without loss of generality, we let the column set K = {1, 2, · · · r}, M = w 1w2 · · ·w r |Wu 1 · · ·Wu t , whereW = w 1w2 · · ·w r . After the normalization, , thenM = w 1w 2 · · ·w r |W ū 1 · · ·W ū t . In the following, we will show that the first column will be chosen fromW . Note S 0 (W) ≥ 0 and U ≥ 0, for j ∈ {1, 2, · · · t}, Then we have Hence the first column will be selected fromW . Without loss of generality, assume thatw 1 is chosen, then in the projection procedure, Hence, Z = 0 |w . Similar result as (13) can be derived. The second column will be selected from K.
For simplicity, assume L columns have been selected from K, let P ⊥ W (L) be the projection onto the orthogonal complement of the columns ofW (L) . Easy get P ⊥ It implies that the L + 1-th column will be selected from K.
Since rank(Mat(M)) = r, we can deduce that rank(Mat(W)) = r. After select the r-th column from K, the residue matrix will degenerate into zero matrix. The results hence follow.

Algorithms for Solving H
After column set K has been determined from Quaternion-SPA, the factor matrixW =M(:, K). The factor matrix H ∈ R r×n + can be computed by solving the following quaternion nonnegative least square problem.
Let (14) can be solved by the following optimization model.
The following quaternion nonnegative least squares method is proposed in [10] for model (15).
Note that QNLS method is simple but rough. The projection onto R r×n + may lead to many zeros in H. It may cause algorithms fail, for example the quaternion nonnegative matrix factorization (QNMF) proposed in [10]. Hence, to address this disadvantage and improve the solution, an iterative method is given in the following section to compute H.

Quaternion Hierarchical Nonnegative Least Squares
As we know that objective function (14) can be further rewritten as,

Because of the nonnegative constrains on H, h pj is obtained by
, j = 1, 2, · · · , n.
Note that the row of H may equal to zero thanks to the projection Π R+ , which will lead to failure of algorithms. To address this problem, one strategy is to use a very small positive number ξ to replace lower bound 0 in the projection Π R+ in practice. Precisely, we define projection Π R +ξ as follows, where ξ ≪ 1 is a very small positive number. The row of H can be updated successively by where H (k+1) p,: .
Here H p,: refers to the p-th row of H. The iteration is presented in Algorithm 3.
are strictly convex. Here H p,: and H :,p refer to the p-th row and p-th column of H respectively. We can derive that the limit points of the iterates (17) are stationary points, based on the results in [5,36]. The convergence results are presented as follows.
Theorem 2. The points of algorithm 3 will converge to stationary points for every ξ > 0.
We remark that QHNLS is inspired by hierarchical nonnegative least squares (HNLS) [9,16] that is an exact block-coordinate descent method for solving nonnegative matrix factorization. The convergence results of Algorithm 3 are similar to HNLS (see Theorem 3 in reference [15]). In the proposed SQMF model,W is quaternion matrix formed by the columns of quaternion matrixM, we only apply QHNLS to compute H.

Numerical Experiments
In this section, we test the algorithms on realistic polarization image for application of image representation, and simulated spectro-polarimetric data set for blind unmixing. All experiments were run on Intel(R) Core(TM) i5-5200 CPU @2.20GHZ with 8GB of RAM using Matlab. We compare our method to the following state-of-the-art methods.
It selects the column with the largest l 2 norm then projects all columns of M ∈ R m×n on the orthogonal complement of the extracted column at each step. Note that SPA can only work on real matrix M ∈ R m×n , we consider the following variant for an input quaternion matrixM ∈ H m×n S .
• SPA * : We determine the column set by applying SPA on the real part of quaternion matrixM ∈ H m×n S , that is S 0 (M).
2. QNMF (quaternion NMF, [10]) is a state-of-the-art quaternion NMF model generalized from the standard NMF. Given a quaternion matrixM ∈ H m×n S , QNMF aims to findW ∈ H m×r S and H ∈ R r×n + such that M −WH is minimize. Quaternion alternating least squares algorithm (QALS) extended from standard ALS is proposed to solve QNMF model.

ImQNMF: We apply quaternion hierarchical ALS to replace quaternion ALS to compute H
for the quaternion NMF model, and refer it to as the improved quaternion NMF (ImQNMF).
The stopping criterion of QNMF and ImQNMF: We will use δ = 10 −4 , and the maximum number of iterations is 1000.

Polarization Image Representation
The polarization image "cover" is from polarization image dataset containing 40 images with several scenes, constructed by Qiu and etc. in [37]. The "cover" picture has 1024 × 1024 pixels, and each pixel contains information of intensity and polarization. That is, Figure 1: Top: image "cover"; Bottom: total intensity S 0 (T), polarization S 1 (T) and S 2 (T) where S 0 (T) ∈ R 1024×1024 + represents the total intensity, and (S 1 (T), S 2 (T), S 3 (T)) stands for its polarization. We show the image "cover" captured by a colour camera and polarization image sensor in Fig.1. The top row of Fig.1 only shows the intensity of "cover" captured by the colour camera with polarization filter at 0 • . The total intensity S 0 (T) and polarization (S 1 (T), S 2 (T)) measured by four intensities with linear polarizers oriented at 0 • , 45 • , 90 • , 135 • , are shown at the the bottom row of Fig.1. Since S 3 (T) = 0 here, we do not show it in Fig.1. We have the following observations: • The total intensities of background and the cover are almost the same, except that of the edge of the cover.
• The differences between polarization S 1 (T), S 2 (T) of background and cover are significant. Especially, the center of cover has higher polarization values than the other parts.
From the observation, the edge and center of the picture "cover" show the main characteristic. As we know that QSPA and SPA * are able to select data points that can be regarded as feature points, it will be interesting to see their identification ability and the reconstruction effect for the "cover" image. Hence, we first divide image "cover"T into 256 small image blocks, and each block has 64 × 64 pixels. The data matrixM ∈ H 4096×256 S is then generated by these 256 blocks; that is, each column ofM is generated by vectorizing every image block. We test QSPA and SPA * onM, and determine r columns, which correspond to the r key image blocks.
In Fig. 2, we show r = 10 images blocks determined by QSPA and SPA * . It is interesting to observe that the areas determined by QSPA are in the background and the center of the cover, while SPA * selects the edge of the cover. The results are surprising coincident with the observation of the characteristics of total intensity S 0 and polarization (S 1 , S 2 ). We remark that QNMF and ImQNMF are not able to identify critical regions that contain latent features from a polarization image.
In Table 1, we report the approximations of the reconstructed image of all the methods, which are defined as follows.
1. The total relative approximation in percentage: We remark that for QNMF, ten initials are generated randomly, and the results with the highest "Appro" value are reported. From Table 1, we observe that: • In terms of the relative approximation of all components, as r, the number of features increases, so make the approximations of all the methods. It is not surprising that ImQNMF gets the highest values, since ImQNMF does not consider extra constraints. QSPA receives the second best results, while SPA * has the worst ones. We note that QNMF fails when r = 30, 50 for all initial guesses.
• In terms of running time, SPA * is the fastest method, which needs less than 1 second. The proposed QSPA is the second fast, and ImQNMF is the slowest method. Especially when r = 50, ImQNMF is almost 10 5 times slower than QSPA, and 10 6 times slower than SPA * .
It is worth mentioning that, QSPA and SPA only need to store the indexes of important columns and the weight matrix H, to represent the polarization image. Precisely, for a givenM ∈ H m×n S and number of features r, to represent the factor matricesW ∈ H m×r S and H ∈ R r×n , QSPA and SPA only store r + rn parameters, while QNMF and ImQNMF need 4rm + rn.
Visualization of reconstructed images by r = 10 features is shown in Fig. 3. It is interesting to find the results of QSPA capture the main characteristics of the polarization image "cover". The results of QNMF and ImQNMF are more smooth, compared to that of SPA * and QSPA. We also notice that SPA * is not able to recover S 2 (T).
In Appendix, we show more visualization results on r = 30, 50.
We conclude the advantages of QSPA in the application of image representation as follows: QSPA can capture important characteristics from a polarization image. It can reconstruct original image well when the number of features is big enough. At the mean time, it runs very fast.

Spectro-Polarimetric Image Unmixing
We simulate spectro-polarimetric data sets from the Urban HSI [19]. Assume (W 0 , H 0 ) be the unmixing ground truth of HSI, W 0 ∈ R m×r0 + stands for the endmember matrix and H 0 ∈ R r0×n + represents the abundance matrix of all endmembers. The spectro-polarimetric data sets are simulated byM =M * + N whereM * =W * H * ∈ H m×n S , and N ∈ H m×n is noise matrix.
• H * ∈ R r×n + is generated from the abundance matrix H 0 of HSI data sets.
• The noise matrix N ∈ H m×n is generated at random normal distribution, that is are generated by the randn function of MATLAB. N is normalized such that N F = M * F . is the noise level defined as To evaluate the quality of solution (W, H) computed by a method, the total relative approximation "Appro" of (19) and relative approximation on the components "app-s l " (l = 0, 1, 2, 3) in (20) will be reported. We also present the relative approximation to ground truth (W * , H * ) in percentage defined as follows.
where π w , and π h are permutations. We test all the methods on three noise levels = [0, 0.05, 0.1] on the simulated data sets. For each noise level, ten such matrices are generated. For each matrix, ten initials are generated randomly for QNMF method, and we report the results of that has the best "Appro" value.

The Simulated Spectro-Polarimetric Urban Dataset
The Urban HSI [19] is taken from Hyper-spectral Digital Imagery Collection Experiment (HY-DICE) air-borne sensors and contains 162 clean spectral bands where each image has dimension 307 × 307. Therefore the "spectral × pixels" data matrix M 0 has dimensions 162 by 94249. The Urban data is mainly composed of 6 types of materials: Road; Grass; Tree; Roof; Metal; Dirt (for more details, see [19]). Its ground truth (W 0 , H 0 ) are , j = 1, 2, · · · , 6.
We simulate six sources spectro-polarimetric dataset by (21) based on Urban HSI with The average quality measures in percent is reported in Table 2. We have the following observations: − Regarding the relative approximation on all components, both QSPA and SPA * can achieve 100% in the noiseless case. As the noise level increases, the relative approximation of all methods decreases. On all noise levels, QSPA gets better results than SPA * . ImQNMF overperforms QNMF; in a certain sense it validates the effectiveness of QHNLS. It is reasonable that ImQNMF performs the best as the noise level ≥ 5%, since it aims to find factor matrices without separability constraints.
− In terms of recovering the ground truth (W * , H * ), QSPA performs the best at the noise level ≤ 5%. ImQNMF is the best at the noise level = 10%. The reason is that the matrix at the noise level = 10% is already far from a separable quaternion matrix. Worth noting that SPA * can recover ground truth in noiseless case; the reason is the real component S 0 (M * ) is 6-separable matrix, i.e., the spectrum already contains 6 endmembers.
− We also report average running time for all the methods. SPA * only needs an average of 0.11 seconds, and QSPA spends around 148 seconds. Both QNMF and ImQNMF are very slow, take more than 2000 seconds in the most cases.
In Appendix, we show some visual results at noise level = 5% for one arbitrary group of simulated polarimetric parameters {α, β} .
As we know that the "different objects with same intensity" problem is challenging in HSI blind unmixing. In the following simulated data set, we will show that the problem can be well solved with the help of polarization. We simulate a polarized Urban dataset containing ten objects, i.e., r = 10, with some different objects having the same intensity.  Table 3: Numerical results (in percent) for 10 sources spectro-polarimetric Urban dataset (r = 10).
• The intensity of source matrix S 0 (W * ) : where h * 7 takes the last 500 positions in h 0 1 where the values are 1, and the last 1000 positions in h 0 1 whose values are between 0 to 1. We report the average quality measures in percent in Table 3 and observe: − In terms of the relative approximation on all components, only QSPA can achieve 100% in a noiseless case. The relative approximation of all methods decreases as the noise level increases. When the noise level ≥ 5%, ImQNMF has the best results. It is also nice to see that QSPA obtains the second best at the noise level = 5%, 10%. We notice that SPA * does not perform well here.
− Regarding recovering the ground truth (W * , H * ), QSPA performs the best. We remark that in the noiseless case, QSPA can recover the factor matrices, while SPA * could not work because the rank of the real component ofM is 6, less than r = 10, the number of the simulated sources. However QNMF takes more than 5600 seconds and ImQNMF needs more than 13000 seconds, much slower than the other two methods.
Since QNMF fails for some initial guesses in this ten sources spectro-polarimetric dataset, we report the number of its failures among ten initials for each generated matrixM at all the noise levels in Fig. 4. We find that QNMF fails more frequently at a negligible noise level, thanks to the singularity ofW caused by too many zeros of H by QNLS. Especially in the noiseless case, it fails at 97 initials in 100 in total. On the right of Fig. 4, we report the identification accuracy of SPA * and QSPA defined as follows.
where K * is the true column indices of pure sources inM * . It shows the proportion of column indices that are correctly identified. It is not surprising to see the accuracies of QSPA are better than SPA * at all levels of noise. The visual results under one arbitrary group of simulated polarimetric parameters {α, β} are presented in Figs. 5-10. The simulated spectro-polarimetric image dataM =WH for three distinct wavelengths indices t = 30, 90, 150 are shown in Fig.5. The sources and the corresponding activations factors from all the methods at noise level 5% are shown in Figs.6-10. Visually from Fig.  6, it is exciting that QSPA can almost recover the activation matrix H at the noise level 5%. From Figs. 7-10, since QSPA and SPA * determineW from the noisy data matrixM, the components of the source matrixW from QSPA and SPA * have some noise compared to the ground truthW * . Compared to the other three methods, visually, the components of the source matrixW from QSPA are much closer to ground truthW * . The visual results indicate that QSPA overperforms the others to get factor matrices (W, H), which verify the results shown in Table 3.
We now make a conclusion on the advantages of QSPA in the application of spectro-polarimetric blind unmixing. QSPA can well identify the underlying objects and sources from spectro-polarimetric data, even in the situation of the same intensity but foreign objects. Its computational time is competitive.

Conclusion
In this paper, we introduced separability into the quaternion matrix factorization and referred the problem to as separable quaternion matrix factorization (SQMF). We studied some properties of SQ-matrices that can be decomposed by SQMF. We showed SQMF is unique up to scaling and permutations. To solve quaternion-valued sources of SQMF, we proposed a fast and efficient method called quaternion successive projection algorithm (QSPA) extended from the successive projection algorithm (SPA). We have demonstrated that QSPA can correctly identity the quaternion-valued sources factor in noiseless cases. To compute the activation factor matrix, we provided a simple algorithm named quaternion hierarchical nonnegative least squares (QHNLS). We tested the algorithms on polarization image for image representation, and simulated spectro-polarimetric data set for blind unmixing. The numerical results showed that the proposed method works promisingly.
Further work include a robust analysis of QSPA, and to design more efficient algorithms for the SQMF problem in the presence of noise.