Deep Transformer Model with Pre-Layer Normalization for COVID-19 Growth Prediction

Coronavirus disease or COVID-19 is an infectious disease caused by the SARS-CoV-2 virus. The first confirmed case caused by this virus was found at the end of December 2019 in Wuhan City, China. This case then spread throughout the world, including Indonesia. Therefore, the COVID-19 case was designated as a global pandemic by WHO. The growth of COVID-19 cases, especially in Indonesia, can be predicted using several approaches, such as the Deep Neural Network (DNN). One of the DNN models that can be used is Deep Transformer which can predict time series. The model is trained with several test scenarios to get the best model. The evaluation is finding the best hyperparameters. Then, further evaluation was carried out using the best hyperparameters setting of the number of prediction days, the optimizer, the number of features, and comparison with the former models of the Long Short-Term Memory (LSTM) and Recurrent Neural Network (RNN). All evaluations used metric of the Mean Absolute Percentage Error (MAPE). Based on the results of the evaluations, Deep Transformer produces the best results when using the Pre-Layer Normalization and predicting one day ahead with a MAPE value of 18.83. Furthermore, the model trained with the Adamax optimizer obtains the best performance among other tested optimizers. The performance of the Deep Transformer also exceeds other test models, which are LSTM and RNN.

gression and classification were carried out using several datasets, which were then compared with other models. The results of this study showed that the TST model outperformed the average evaluation value for time series regression and time series classification. Furthermore, Zhou et al. [25] proposed the Informer, a long-sequence time-series forecasting (LSTF) model. The Informer model gets test results that exceed the LogTrans, Reformer, LSTMa, Long-and Short-term Time-series network (LSTnet) models, DeepAR, ARIMA, and Prophet.
Xiong et al. [21] experimented on the Transformer model [14] by comparing the position of the normalization layer and the warm-up stage of the learning rate model. From the experimental results, the original Transformer model that puts the normalization layer outside the residual block gets the expected value of the gradient from the parameter near the output layer, which is large at initialization and causes unstable training when using a large learning rate. Furthermore, it is found that the Transformer model with the normalization layer position in the residual block can be trained without a warm-up stage and converges more quickly.
There are several studies conducted to predict positive cases of COVID-19. One of them is a study by Yudistira [22] which uses multivariate Long Short-Term Memory (LSTM) to predict the growth of positive cases of COVID-19 in the world. In this study, several test scenarios were carried out and the results showed that the LSTM model could exceed the performance of ordinary Recurrent Neural Network (RNN) and Vector Auto Regression (VAR). Then, Yudistira et al. [23] analyzed which variables contributed to the time series prediction in COVID-19 cases. These variables include demographics, economy, health, government, education, environment, and COVID-19. The analysis was carried out by making a saliency map using gradient-based visual attribution for each variable and its time dimension. The results of the test obtained a saliency map to be able to understand what factors affect the increase or decrease of cases.
Based on the research results described previously, no studies predict the growth in the number of COVID-19 cases in Indonesia using the Deep Transformer model. Experimental research that compares the position of the normalization layer has only been carried out on the original Transformer model, not yet on the Transformer for time series prediction. Experiments for hyperparameter optimization have also not been carried out by all the previously described studies. Therefore, this study uses the Deep Transformer model [20] to predict the growth in the number of COVID-19 cases in Indonesia, compare the position of the normalization layer to the predicted results, and optimize several hyperparameters. The results of this study are expected to be helpful in various scientific fields in making decisions during this pandemic.
In this research, several test scenarios will be carried out. First, we test to get the best hyperparameters from the Deep Transformer model. The hyperparameters are the size of the model, the number of encoder and decoder blocks, the dimensions of the Feed-Forward Network, Pre-Layer, and Post-Layer, and the amount of time lag. Then, after obtaining the best hyperparameters, several advanced tests were carried out: testing the predictions for the next few days, testing the number of input features, and finally comparing them with the Recurrent Neural Networks (RNN) and Long Short-Term Memory (LSTM). This paper will be divided into several sections. Section 1 explains the introduction and background of conducting this research. Section 2 describes previous research that is similar to this study. Finally, section 3 describes the methods used in the research, starting with the normalization of the data used, the Deep Transformer model, Mean Absolute Percentage Error (MAPE) to calculate errors and related to the dataset.

Related Work
Research that has been done by Vaswani et al. [14] proposes a Transformer model that is used to perform language translation. The transformer is a model architecture that relies entirely on attention mechanisms to draw global dependencies between input and output. The Transformers enable significantly more parallelization and greatly improve translation quality after training for just twelve hours on eight P100 GPUs. The attention mechanism used is Scaled Dot-Product Attention and Multi-head Attention. Product Attention calculates the attention function on a set of queries simultaneously. Then it is entered into the matrix. Multi-head Attention allows the model to pay attention to information from different sub-spaces of representation at different positions.
Li et al. [8] conducted a study using the Transformer model for time series prediction. The research obtained satisfactory results, but there are two shortcomings: locality-agnostic and memory bottleneck. The first shortcoming can be overcome by the convolutional self-attention method, during the second shortcoming by the LogSparse method, so the proposed model name is LogSparse Transformer. This model is then compared with several popular models for time series prediction. These are Autoregressive Integrated Moving Average (ARIMA), Seasonal Error Trend (ETS), Temporally Regularized Matrix Factorization (TRMF), RNN-based autoregressive model (DeepAR), RNN-based state, and Space model (DeepState). It is found that the proposed model can exceed the performance of other models.
Furthermore, Wu et al. [20] also conducted research that developed a Transformer model so that it can be used to make time-series predictions. Furthermore, Zerveas et al. [24] conducted a study that developed the Transformer model into a Time Series Transformer (TST) model to predict multivariate time series. In this study, time series regression and time series classification were carried out using several datasets, which were then compared with other models such as Support Vector Regression (SVR), Random Forest, XGBoost, Rocket, LSTM, and several others. The TST model is divided into two tests: TST with supervised training and TST with unsupervised pre-training. The results of these tests showed that the average value of the evaluation for time series regression and time series classification was superior to the TST model with unsupervised pre-training and continued with TST with supervised training.
Furthermore, Zhou et al. [25] proposed a Transformer-based model that can predict a long sequence of time series or Long Sequence Time-series Forecasting (LSTF) called the Informer. This model applies the ProbSparse self-attention mechanism. The performance of the Informer was tested by comparing it with several models: LogTrans, Reformer, LSTMa, Long-and Short-term Time-series network (LSTnet), DeepAR, ARIMA, and Prophet. From the test results, it was found that the Informer model almost surpassed all the models tested.
Xiong et al. [21] experimented on the Transformer model by comparing the position of the normalization layer and the warm-up stage of the learning rate model. The experiment was carried out by combining the location of the normalization layer before or after the Transformer layer and using a warm-up stage, which is used or not on the translator machine. From the experimental results, the original Transformer model, which puts the normalization layer outside the residual block, gets the expected gradient value of the parameter near the output layer, which is large at initialization. Therefore, it causes unstable training when using a large learning rate. Furthermore, it is also found that the Transformer model with the normalization layer position in the residual block can be trained without a warm-up stage and converges more quickly.
Based on previous studies, it can be said that the Transformer model produces a time series prediction model with robust performance. This study conducts time-series predictions using the Deep Transformer model on Indonesia's daily growth dataset of COVID-19 data. Moreover, we assume that the position of the normalization layer influences the optimal parameters generated during learning.

Min-Max Normalization
Min-max normalization according to Patro and Sahu [11] is a scaling technique or mapping technique with the aim of creating a new range specifically with predetermined limits. Min-max normalization can be calculated using the formula in Equation 1. Variable A' is the min-max normalized data, A is the original data, D and C are the lower and upper ranges of the specified limits or can be written as [C,D].

Deep Transformer Model
Based on previous research by Wu et al. [20], the architecture of the Deep Transformer model follows the original architecture of the Transformer by Vaswani et al. [14]. Deep Transformer consists of two layers, those are encoder and decoder. Unlike the original architecture, in this model the input and output embedding layers are changed to input layers, which are fully connected layers (see Fig. 1). Encoder The encoder consists of an input layer, a positional encoding layer, and a stack of four identical encoder layers. At the input layer, input in the form of time series data is mapped into d dimensional vector model through a fullyconnected network. In this step, a multi-head attention mechanism is carried out which is an essential step for the model. The positional encoding stage is carried out using the sine and cosine functions which aim to encode sequential information in time series data by adding element-wised from the input vector to the positional encoding vector. The resulting vectors are then fed into the four encoder layers. Each encoder layer consists of two sub-layers, namely the self-attention sub-layer and the fully-connected feed-forward sub-layer. Each sublayer is followed by a normalization layer. The encoder generates d dimensional vector to proceed to the decoder [20].
Decoder The decoder consists of an input layer, one or several identical decoder layers, and an output layer. The input to the decoder starts from the last data point on the input of the encoder. At the input layer, the decoder input is mapped into a dimensional vector. In addition to the two sub-layers in each encoder layer, a third sub-layer is inserted in the decoder, namely the self-attention mechanism at the encoder output. Finally, there is an output layer which maps the output of the last decoder layer to the targeted time sequence. The decoder also uses look-ahead masking and a one-position offset between the decoder input and the target output in the decoder to ensure that the prediction of the time series data point only depends on the previous data point [20].
Positional Encoding Positional Encoding (PE) is entering some relative or absolute position information of each value in the sequence [14]. In this model, PE is added to the Pre-Layer and Post-Layer input embedding at the bottom of the encoder and decoder stack. PE has the same dimensions as the embedding, so they can be added together [20].
Multi Head Attention Attention is a function that maps a query and a set of key-value pairs to output. Query, key, value, and output are in vector form. The output value is calculated as the sum of the weights of the values. The weight for each value is calculated by a query function with the corresponding key [14,8,9]. Multi Head Attention is a module on the attention mechanism that calculates the Attention function several times in parallel. The output results from each head are combined by adding up the matrix values (concat) of each and linearly converted into the expected dimensions [14,8,9].
Based on the Deep Transformer model by Wu et al [20]. There are three types of Multi Head Attention, those are Multi-Head Self-Attention (MHSA), Multi-Head Causal-Attention (MHCA), and Encoder-Decoder Attention. MHCA is different from MHSA, where MHCA will do masking for the next prediction result, unlike MHSA where the value for each data is known. Then, in the Encoder-Decoder Attention key and value are obtained from the encoder results while the query is obtained from the MHCA results. Figure 1, it can be seen that between the sub-layers of the encoder and decoder layers there is an Add and Normalize layer [20]. The add layer is the residual connection layer [6], which is a skip-connection that learns the residual function by referring to the input layer, rather than studying the unreferenced function.

Residual Connection and Layer Normalization In
The normalization layer is designed to overcome the disadvantages of batch normalization [7] such as computationally expensive and batch size sensitive performance. Unlike batch normalization, layer normalization is a normalization method that directly estimates the normalization statistics from the inputs added to the neurons in the hidden layer so that normalization does not create new dependencies between training cases. The normalization layer can be calculated using Equation   where the normalization layer is placed between the residual connections, while in Pre-Layer Normalization (Pre-LN), the normalization layer is placed inside the residual connection (see Fig. 2). Xiong et al. [21] compared the performance of Pre-LN and Post-LN on language translation testing and found that the BLEU score and validation loss of the Pre-LN Transformer were better than post-LN. Furthermore, it was also found that Pre-LN converges faster than Post-LN and also found that changes in the position of the normalization layer affect changes in the optimizer.
Feed Fordward Networks Each layer on the encoder and decoder consists of fully connected feed-forward networks called Position-wise Feed-Forward Networks (FFN). Position-wise FFN is a network with a feed-forward layer consisting of two dense layers used in the last dimension, which means the same dense layer is used for each position item in an order called point-wise. Between the two dense layers, there is a ReLU activation function [14].

Recurrent Neural Networks
Recurrent Neural Networks (RNN) is a neural network architecture that aims to recognize patterns from sequence data. RNN can be used on time series data, genomes, handwriting, images and even to a higher level can be used in speech recognition. Mathematically, the hidden state and output of the RNN can be written in the form of an equation, those are Equation 5 and Equation 6.
The hidden state and the input at time t, respectively, are assumed to be H t and X t . Next, the weights will be used, those are W hh for the hidden state and W xh for input, while the bias will be b h . The output at time t is assumed to be O t with a weight of W ho and a bias of b o . Finally, the results of the calculations for both will be passed to the activation function φ such as logistic sigmoid or tanh.

Long Short-Term Memory
LSTM is a neural network that can handle vanishing gradient problems in RNN by using gated cells that can store information. In LSTM, there are three gates, namely input gate I t to read data into the cell, output gate O t to read cell entries, and forget gate F t to reset the contents of the cell. The three gates can be written into equations respectively in Equation 7, Equation 8, and Equation 9. The variables W xi , W xo , and W xf are the data weights of each gate. The variables W ho , W ho , and W hf are the hidden state weights of each gate. While the variables b i , b o , and b f are the bias. The calculation results from the three gates will be passed to the sigmoid activation function.
Next, there is a cellC t as a memory candidate which has the same calculation as before, but uses the activation function tanh. the calculation ofC t will use the weights W xc and W hc as well as the bias b c which can be seen in Equation 10. Finally, the value of the new memory will be calculated, namely C t and H t which can be seen in Equation 11 and Equation 12, where denotes elementwise multiplicationC

Mean Absolute Percentage Error
According to Adhikari and Agrawal [1], Mean Absolute Percentage Error (MAPE) is a performance measure to measure the performance of a time series prediction or forecasting model. It can be seen in Equation 13, MAPE represents the percentage of the average absolute error. The variable e is the error value of the prediction, y is the actual value, and n is the measure of the test data. The error value is obtained from the difference between the predicted value and the actual value.
M AP E = 1 n n t=1 e t y t * 100 (13)

Data
This study uses secondary data as the object of research. The data was taken from the public API on the official page of the COVID-19 Handling Task Force in Indonesia, that is covid.go.id. The data taken is in the form of a JSON (JavaScript Object Notation) file containing overall information on COVID-19 cases, daily updates, and accumulation from the beginning of the case to the present. The data will be normalized using min-max normalization with a value range of -1 to 1. The data is then divided into three, those are training data, evaluation data, and test data. Test data is the data used to test the model that has been trained. The test data were taken from the last 60 days of the total data. Training data is the data used to train the model and is determined as much as 70% of the data from the beginning except for test data. Evaluation data is used to evaluate the model in each epoch and is determined as much as 30% of the data from the back except for test data.

Experimental Results
The test was carried out on 750 initial COVID-19 data in Indonesia with the last 60 days used as test data. First, testing is carried out to obtain optimal hyperparameters. The hyperparameters are: embedding size; the number of encoder and decoder blocks; FFN, pre-layer, and post-layer dimensions; and the amount of time lag. Furthermore, further testing is carried out, those are testing the number of predictions for the next day, optimizer used in training, the number of features, and finally comparing the Deep Transformer model with the LSTM and RNN models.

Hyperparameter Testing
Prior to testing to obtain optimal hyperparameters, the initial hyperparameters are set as benchmarks or defaults, that is: input step or time lag of 7; predict step or the number of predictions ahead of 1; d size of 64; epoch maximum of 300; the number of encoders and decoders is 2 each; the number of heads from each encoder and decoder is 1; there are 50 hidden layers in feed forward, prelayer, and post-layer each; drop out on pre-layer, post-layer, and model is 0.2; the position of the normalized layer is Pre-LN; At the beginning of the training the learning rate is determined to be 1 and will decrease as the training is carried out. The optimizer used is the Adam optimizer with the LambdaLR scheduler with a function created by Wu et al. [20].
Testing on Embedding Size The first test was carried out to obtain the best d size among several sizes, those are 32, 64, 128, and 258. Table 1  Based on these test results, it was found that the best MAPE values for both types of normalized layers were equally obtained from the d size of 64. Therefore, the d size of 64 was used for advanced testing. Testing on the Number of Encoder and Decoder Blocks Furthermore, testing is carried out to obtain the best combination of the number of encoder and decoder blocks. The combination of the number of encoder and decoder blocks are 1-1, 2-2, 4-4, 2-4, and 4-2. Table 2 shows the result of the test. Based on the results of these tests, it was found that the combination of the number of encoder and decoder blocks which both amounted to 2 got the best MAPE value for the second type of normalization layer. Therefore, the number of encoder and decoder blocks used for advanced testing is 2.
Testing the Dimensions of the Feed Forward Network, Pre-Layer, and Post-Layer The next test combines the dimensions of the Feed Fordward Network layer, Pre-Layer, and Post-Layer. The size dimensions are either 30, 50, or 100. Table 3 shows the result of the test.
Based on the test results, it was found that the Pre-Layer with dimensions of 100 and the rest with dimensions of 50 got the best MAPE value for the Pre-LN model. Meanwhile, the best post-LN model yields the best performance using dimensions of FFN of 100 and the rest has dimensions of 50. Then the average is calculated for models with both normalization layers to determine which combination will be used for forthcoming testing. Finally, the best average is obtained for the FFN with the dimension of 100 and the remaining of 50.
Testing the Amount of Time Lag Next, the number of time lags was tested, which were 4, 7, 14, and 30 days. Table 4 is the result of testing the amount of time lag. Based on these results, it was found that the 7-day time lag got the best MAPE value for both types of normalized layers. Therefore, a time lag of 7 days was used for further testing.

Testing and Analysis of the Number of Predicted Days in the Future
After obtaining the best hyperparameters from the Deep Transformer model, further testing is carried out, testing the number of predictions for the next day. The test is performed using the d size of 64, the number of encoder and decoder blocks of 2 each, FFN dimension of 100, Pre-Layer and Post-Layer dimensions of 50, and the amount of time lag of 7 days. Then, judging from the previous four tests, the Pre-LN normalization layer got a better error value than the post-LN. Therefore, Pre-LN will be used for further testing. First, testing is carried out for predictions one day ahead of the last 60 days. The test was carried out ten times in a row and the average with its standard deviation of the error was calculated. The average MAPE value obtained is 22.298 and the standard deviation is 1.363. Figure 3 shows a graph comparing the average prediction results to the original value along with the range of prediction results from 10 trials. Second, testing is carried out for the next two days with the same data and method as the previous test. These tests obtained an average of 31.15 and a standard deviation of 2.901. The comparison of the average prediction results to the original value can be seen in Figure 4. Next, testing is carried out for the next four days with the same data and method as the previous test. These tests obtained an average of 37.428 and a standard deviation of 1.635. The comparison of the average prediction results to the original value can be seen in Figure 5. Furthermore, testing is carried out for the next 7 days with the same data and method as the previous test. From these tests, obtained an average of 42.35 and a standard deviation of 2.893. The comparison of the average prediction results to the original value can be seen in Figure 6. Based on the results of the four tests, it can be seen that the prediction error value for the next 1 day is the smallest and predictions for the next 2, 4, and 7 days will increase along with the increase in the number of predictions for the next day. Therefore, it can be said that the more days that will be predicted in the future, the resulting error value will increase. The comparison of the average error value and standard deviation can be seen in Table 5.

Testing and Analysis of the Optimizer Used
Furthermore, research was conducted on the optimizer used in training. The optimizers tested are Adam, AdamW, Adamax, Adagrad, Adadelta, SGD, and RMSprop. The hyperparameters used and the method of testing is the same as testing the predictions for the next day. For each optimizer, a LambdaLR scheduler is used with a function created by Wu et al. [20]. The results of comparing the MAPE values of the tested optimizer can be seen in Figure 7. In contrast, the graph of the comparison of the average prediction results with the ground-truth along with the range of the prediction results can be seen in Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13, and Figure 14.   Based on the results of the seven tests, it can be seen that the prediction error value using the Adamax optimizer is the smallest followed by Adam, AdamW, SGD, Adagrad, RMSprop, and the worst is the Adadelta optimizer. Therefore, it can be said that the Adamax optimizer can train the model better with a MAPE value of 22.343 compared to other optimizers. The comparison of the average error value and standard deviation can be seen in Table 6.

Testing and Analysis of the Number of Features
Furthermore, the number of features is tested, those are 1 feature (univariate) and multivariate with 3 features. One feature is using the number of positive cases, while three features are the number of positive cases, the number of deaths,     and the number of recoveries. Both tests have one output, that is the number of positive cases. The hyperparameters used are the same as the previous advanced tests. Testing is done by repeating the training and testing process several times for the last 60 days of data until the best results are obtained.  After several iterations of testing, the best result is that the model with one feature gets a MAPE value of 18.83, while the three features get a MAPE value of 25.24. The comparison graph between the prediction results and the original value of the prediction of one feature and 3 features, respectively, can be seen in Figure 15 and Figure 16.

Testing and Analysis of Other Models
Finally, a test was conducted to compare the Deep Transformer model with other models, those are Long Short-Term Memory (LSTM) and Recurrent Neural Network (RNN). The Deep Transformer model uses hyperparameters as before with a Pre-LN normalization layer, while the LSTM and RNN models use 16 hidden layers, one layer, 2000 epochs, and use the Adam optimizer with a learning rate of 0.01. Testing is carried out by conducting the training and testing process several times for the data for the last 60 days until the best results are obtained.
After several training and testing processes were carried out, the best MAPE results were obtained for each model. The Deep Transformer model got the best MAPE result of 18.83, the LSTM model got the MAPE result of 24.33, while the RNN model got the MAPE result of 23.15. A more detailed comparison of the MAPE values of each model can be seen in Figure 17, while a comparison of the prediction results for each model can be seen in Figure 18.  Based on the test results of the three models, it can be said that the Deep Transformer model outperforms the other two models by a considerable difference. The performance of the LSTM and RNN models is almost the same, but the RNN model is superior to LSTM with the same number of hidden layers.

Conclusion
The design of the Deep Transformer model in predicting the growth of COVID-19 cases in Indonesia starts by taking datasets from the public API on the official page of the COVID-19 Handling Task Force in Indonesia. Then, from the dataset, the required column is taken. Those are the date and the positive number of 750 days since the beginning of the case. Furthermore, the Deep Transformer model was developed in such a way that it could be used to predict the number of positive cases from the COVID-19 dataset in Indonesia. After the model is completed, training is carried out to obtain a model that can optimally predict positive cases by testing several existing hyperparameters. After obtaining the best hyperparameters, the testing process is carried out to predict the number of positive cases from the last 60 days of data. Finally, the test results will calculate the error value using Mean Absolute Percentage Error (MAPE).
Based on the results of several tests, the best hyperparameters were obtained. These are: embedding size of 64, the number of encoder and decoder blocks of two each, dimensions of the FFN of 100, Pre-Layer and Post-Layer dimensions are 50, the amount of time lag of 7 days, drop out on Pre-Layer of 0.2, the normalized layer of Pre-LN, and using the Adam optimizer of Adam with a learning rate initialization of 1. The loss will decrease as the training progresses with a total of 300 epochs. Using these hyperparameters, the Deep Transformer model gets a MAPE value of 18.83.
Suggestions for further experiments are comparing the results of model training with several optimizers other than Adam. Furthermore, developing a proper Deep Transformer model with more than one feature or multivariate input that can predict more than one day is promising.


extfieldsofstudy:
Computer Science
Mathematics

text:
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.