**1. Introduction**

Due to the use of hundreds of spectral bands, hyperspectral imaging (HSI) generally has enormous data volume and contains vast amount of information. This special characteristic results in several issues. First, the inter-band correlation of HSI is very high, and adjacent bands may contain redundant spectral information. This would lead to the well-known problem called the "curse of dimensionality" in data analysis. Second, the computational complexity of processing HSI data is very high. Third, storing HSI data usually requires large amount of disc space. Finally, transmitting HSI data requires higher bandwidths. Under such circumstances, removing partial data without significant loss of an image's spectral information is necessary. One commonly used approach is band selection (BS) [1]. BS takes advantage of such high-band correlation to remove the redundant bands, in order to achieve a wide range of applications, such as dimensionality reduction, data storage, data transmission, target detection, and classification.

Many different types of BS algorithms [2–15] have been proposed in the past two decades. Most of them make an assumption that the BS problem is an optimization problem, which maximizes or minimizes a pre-defined objective function that can measure the amount of information or inter-band

redundancy contained by the currently selected bands. For instance, Keshava et al. [2] adopts common distance measures, such as Euclidean distance and a spectral angle mapper, to measure the similarity of bands. Du et al. [3] uses linear regression and orthogonal subspace projection to sequentially select a new band based on the complementary vector space. Chang et al. [4] regards a spectral band as a signal, and designed a BS method based on a constrained energy minimization algorithm. Martínez-Usó et al. [5] introduces a hierarchical clustering-based method to select bands from the band groups that were pre-clustered by a specific method, with some informational measures. In addition, there has been a lot of research that adopts the perspective of information theory. The works in [6–8] adopt mutual information as the similarity index to find the maximum information spanned by the bands. In recent years, sparse regression models were used to perform BS. For instance, Sun et al. [9–11] uses a self-sparse regression (SSR) model to select bands. In an SSR problem, finding a new basis is equivalent to selecting the most representation bands in an HSI image. Lai et al. [12] proposed a SSR-based BS method, which adopts an orthogonal matching pursuit (OMP) algorithm to sequentially find the next band in an efficient manner. For the purpose of HSI data transmission, Chang et al. [13,14] proposed a progressive BS (PBS) method for spectral unmixing. In PBS, the bands are first prioritized by a certain criterion, and are transmitted by the order of prioritization. Du et al. [15] proposed a BS-based dimensional reduction method for change detection in multi-temporal hyperspectral images. Except for the above-mentioned work, there were still various kinds of BS research published in the literature [16–25].

Since many existing HSI algorithms were designed to deal with off-line problems, they usually require more computing time for pursuing extremely high accuracy. Undoubtedly, those sophisticated methods may not be appropriate for dealing with timeliness applications. Hence, developing real-time and efficient approaches is urgently needed. Based on that point, some research started to develop on-board computing approaches [26–29], where the data compression could be carried out on satellites to reduce the demand of downlink bandwidth, focusing additionally on real-time image classification with the aid of GPU acceleration [30,31].

On the other hand, processing HSI data in a real-time progressive fashion under the data transmission held from satellite to ground station has become another attractive topic. The term "progressive processing" for HSI data originated from [32]. According to [32], the term refers to a type of signal processing method. It decomposes entire data processing into a finite number of stages, and processes data progressively stage by stage, in the sense that the results obtained by previous stages can be used to update or improve the results to be processed in subsequent stages. Recently, the concept of progressive processing was utilized to process HSI data. For example, in order to achieve real-time spectral unmixing during band transmission, a revolutionary concept called progressive band processing (PBP) [33], was developed. PBP assumes that the HSI data is transmitted by band sequential (BSQ) format, and processes the image immediately after receiving a new band. Under the PBP framework, many common topics, such as spectral unmixing, target detection, and anomaly detection could be realized [33–35]. Unlike the traditional algorithms, which can only be implemented on the collected data, PBP makes use of recursive processing, to instantly process the current data segment. More specifically, PBP methods preserve the past information obtained at previous stage, and use it to accelerate the computation at the current stage. The most attractive feature of PBP is that the computing time will not significantly increase as the growth of number of received bands. Therefore, PBP algorithms have significant potential for the uses of applications in data transmission and communication.

Based on the importance of BS, and the urgen<sup>t</sup> demand of real-time algorithms, this paper proposes a real-time progressive BS method for applications related to data transmission. To our best knowledge, there were no BS methods proposed for this purpose before. In our case, we assume that HSI data is transmitted in band-interleaved-by-pixel (BIP) or band-interleaved-by-sample (BIS) format, and is processed pixel by pixel (sample by sample). According to [32], such a process operated on BIP/BIS transmission is called progressive sample processing (PSP). To sum up, we developed a novel

concept called progressive sample processing of band selection (PSP-BS), which combines BS with PSP. Figure 1 illustrates the difference of traditional BS and PSP-BS.

To realize PSP-BS, we had to choose an appropriate BS algorithm as the core of PSP-BS. In this paper, a self-sparse regression BS algorithm, called orthogonal matching pursuit-based BS (OMPBS) [12], has been adopted because (1) the OMPBS can sequentially find the most information-complementary bands from the image, without any prior knowledge or any complicated optimization scheme; and (2) The algorithm is easily re-formulated to fulfill recursive processing through mathematical decomposition of the optimization equations of OMPBS. In the later sections, we first introduce OMPBS, and then derive its recursive version, called recursive OMPBS (Re-OMPBS), as the key algorithm to realize PSP-BS. In addition, we make an assumption that the pixel transmission sequence (PTS) could be defined before transmission. We propose three PTSs to make BS results converge more quickly during the PSP-BS. We call Re-OMPBS combined with PTS PSP-OMPBS.

To evaluate the effectiveness of PSP-OMPBS, two real hyperspectral datasets were used in the experiments. We conducted a comparison of computing time and the corresponding BS accuracy at every time stage in a progressive manner. Land cover classification was also adopted to evaluate the quality of the selected bands. According to the experimental results, the computational efficiency of PSP-OMPBS is stably high so that the process could be carried out in real time. The BS accuracy can be improved early on by using particular PTSs.

**Figure 1.** An illustration of the difference of traditional band selection (BS) and progressive sample processing of band selection (PSP-BS). (**a**) The traditional BS algorithm is implemented on the collected hyperspectral image (HSI) data. (**b**) The PSP-BS is implemented in data sample transmission. In PSP-BS, the real-time BS monitoring can be achieved.

#### **2. Orthogonal Matching Pursuit-Based BS (OMPBS)**

The OMPBS [12] is a sequential BS method based on a self-sparse regression (SSR) model [36]. It aims to find a set of representative bands that can represent all bands based on the minimization of reconstructed errors in the SSR model, by using an orthogonal matching pursuit (OMP) algorithm [37]. Suppose *L* is the number of bands and *N* is number of pixels in a hyperspectral image. Let **b***i* ∈ *R N*×1 presents the *i*th band vector. In an SSR model, both the observation matrix as well as the dictionary matrix are the set of band vectors denoted by **B** = [**b**1, **b**2,..., **b***L*] and the coefficient matrix is denoted by **C**. Obviously, the goal is to find the optimal coefficient matrix **C** that minimizes the reconstructed error ||**B** − **BC**||<sup>2</sup> *F*. Since only *k* bands need to be selected, we impose the sparsity constraint ||**C**||0,2 *p* where *l*0/*l*2 norm ||**C**||0,2 counts the number of the non-zero rows in **C**. Thus, our sparse-BS optimization problem is

$$\hat{\mathbf{C}} = \arg\min\_{\mathbf{c}} ||\mathbf{B} - \mathbf{B}\mathbf{C}||\_F^2 \dots \\ \text{s.t.} \ ||\mathbf{C}||\_{0,2} \le p\_\prime \tag{1}$$

where ||**B** − **BC**||<sup>2</sup> *F* is the reconstructed error, ||**C**||0,2 *p* is the sparsity constraint, and *p* is the sparsity level. The sub-optimization problem (find **C**) is solved by least square estimator **Qˆ** = **P***T***P** −1 **<sup>P</sup>***T***B**,

where **P** is temporary basis matrix composed by the candidate bands. The original OMPBS algorithm is composed of a doubly-nested loop. So, it is required to solve *L*!/(*L* − *p*)! times of (1), so the overall computation is time-consuming. To improve its efficiency, we revised the original algorithm by residual perspective. The revised algorithm (Algorithm 1) is shown as follows.

## **Algorithm 1. OMPBS**

Objective: Find *p* representative bands for a HSI cube Input: **B** (*N* × *L* matrix), sparsity = *p* Output: Band indices set <sup>Ω</sup>*p* = .*<sup>b</sup>*1, *b*2,..., *bp*/ Steps: 1. Initialization: Set Ω0 = ∅, **R**0 = **B**, **P**0 = ∅, *j* = 1. 2. At *j*-th iteration, calculate the row norm vector **k** ∈ *L* × 1 of matrix **R***Tj* **<sup>R</sup>***j*. 3. Find the row index *bj* whose value is maximum at **k**. Then set **P***j* = \***<sup>P</sup>***j*−<sup>1</sup> **b***bj* +, Ω*j* = <sup>Ω</sup>*j*−<sup>1</sup> ∩ ,*bj*-. 4. Solve the sub optimization problem **Q ˆ** *j* = arg min**Q***j*------**B** − **<sup>P</sup>***j***Q**---|2*F* by **Q***j* = (**P***Tj* **<sup>P</sup>***j*)<sup>−</sup>1**P***Tj* **B**. 5. Update the residual matrix **<sup>R</sup>***j*+<sup>1</sup> = **B** − **<sup>P</sup>***j***Q***j*. 6. If *j* = *p*, break; otherwise, set *j* ← *j* + 1 and go to Step 2.

#### **3. Progressive Sample Processing-Based OMPBS (PSP-OMPBS)**

The proposed PSP-OMPBS method is fully introduced in this section. It consists of two parts, as shown in Figure 2. The first part controls the data transmission order, in which the pixel transmission sequence will be formed. The second part is the real-time BS computing core, recursive OMPBS (Re-OMPBS). Section 3.1 introduces the derivation of Re-OMPBS, where we convert some steps in OMPBS to their recursive forms. Section 3.2 describes the condition for implementing recursive equations introduced in Section 3.1. Section 3.3 summaries the details of the Re-OMPBS algorithm. Finally, the formation of the pixel transmission sequence is introduced in Section 3.4.

**Figure 2.** Conceptual flowchart of progressive sample processing-based OMPBS (PSP-OMPBS). The center part is the proposed PSP-OMPBS method which includes two blocks: formation of the pixel transmission sequence (PTS), and implementation of the recursive OMPBS (Re-OMPBS) algorithm. The transmission is carried out from the PTS block to the Re-OMPBS block.

#### *3.1. Derivation of Recursive OMPBS (Re-OMPBS)*

In this section, we aim to integrate OMPBS into PSP. According to the definition, PSP-BS can instantly acquire the current BS result during the data transmission by using two types of information: the processed information provided in previous stage (i.e., the (*n* − 1)th stage), and the new information provided in current stage (i.e., *n*th stage). More specifically, PSP-BS can immediately update the BS result when a new pixel (the *n*th pixel) is received. Assume that *n* presents the number of currently received pixels. We would like to utilize the past information obtained in the (*n* − 1)th stage to reduce the computational burden of the *n*th stage. In doing so, the optimization process of OMPBS must be broken down.

#### 3.1.1. Decomposition of the Least Square Estimator

In the following derivation, we assume that *j* is a fixed number and do not annotate it for each term. In the OMPBS algorithm, there are two calculations that can be decomposed. The first one is the least square estimator in Step 4:

$$\mathbf{Q}(n) = \left(\mathbf{P}(n)^T \times \mathbf{P}(n)\right)^{-1} \times \mathbf{P}(n)^T \times \mathbf{B}(n). \tag{2}$$

As we know, computing Equation (2) is required for each iteration *j* in the inner loop of OMPBS. This computation would be time-consuming when *n* is large. This would be unfavorable to progressive processing. To resolve this issue, we must convert Equation (2) to its recursive version. To do that, we define two matrices:

$$\mathbf{B}(n) = \begin{bmatrix} b\_{11} & b\_{12} & \dots & b\_{1L} \\ b\_{21} & b\_{22} & \dots & b\_{2L} \\ \vdots & \vdots & \ddots & \vdots \\ b\_{n1} & b\_{n2} & \dots & b\_{nL} \end{bmatrix}, \quad \mathbf{P}(n) = \begin{bmatrix} p\_{11} & p\_{12} & \dots & p\_{1p} \\ p\_{21} & p\_{22} & \dots & p\_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ p\_{n1} & p\_{n2} & \dots & p\_{np} \end{bmatrix}. \tag{3}$$

**<sup>B</sup>**(*n*) represents the total pixel data after the *n*th pixel is received, and **<sup>P</sup>**(*n*) presents the data's temporary basis matrix, composed by selected band vectors. Let **r***n* = [*bn*1, *bn*2,..., *bnL*] be the new pixel and **p***n* = %*pn*1, *pn*2,..., *pnp*&. Next, we define auxiliary matrix **<sup>H</sup>**(*n*)

$$\mathbf{H}(n) = \mathbf{P}(n)^T \times \mathbf{P}(n) = \left[ \begin{array}{c} \mathbf{P}(n-1)^T & \widetilde{\mathbf{p}}\_n^T \end{array} \right] \left[ \begin{array}{c} \mathbf{P}(n-1) \\ \widetilde{\mathbf{p}}\_n \end{array} \right] = \mathbf{P}(n-1)^T \times \mathbf{P}(n-1) + \widetilde{\mathbf{p}}\_n^T \widetilde{\mathbf{p}}\_n. \tag{4}$$

Substituting (4) into (2), we ge<sup>t</sup>

$$\mathbf{Q}(n) = \left(\mathbf{P}(n-1)^T \times \mathbf{P}(n-1) + \breve{\mathbf{p}}\_n^T \breve{\mathbf{p}}\_n\right)^{-1} \times \mathbf{P}(n)^T \times \mathbf{B}(n). \tag{5}$$

Now use Woodberry's Identity [38]

$$\left(\mathbf{A} + \mathbf{c}\mathbf{c}^{T}\right)^{-1} = \mathbf{A}^{-1} - \frac{\left(\mathbf{A}^{-1}\mathbf{c}\right)\left(\mathbf{c}^{T}\mathbf{A}^{-1}\right)}{1 + \mathbf{c}^{T}\mathbf{A}^{-1}\mathbf{c}}\tag{6}$$

to decompose Equation (5). Let **A** = **<sup>P</sup>**(*n* − <sup>1</sup>)*<sup>T</sup>***P**(*n* − <sup>1</sup>), **c** = **p***Tn* , *v* = **A**−<sup>1</sup> and *ρ* = **<sup>c</sup>***T***A**−1**c**, and Equation (6) becomes

$$\mathbf{Q}(n) = \left( \mathbf{H}(n-1)^{-1} - \frac{\mathbf{v}\_{n\mid(n-1)} \mathbf{v}\_{n\mid(n-1)}^T}{1 + \rho\_{n\mid(n-1)}} \right) \left[ \begin{array}{c} \mathbf{P}(n-1) \\ \widetilde{\mathbf{p}}\_n \end{array} \right]^T \left[ \begin{array}{c} \mathbf{B}(n-1) \\ \mathbf{r}\_n \end{array} \right] \tag{7}$$

where **<sup>v</sup>***n*|(*<sup>n</sup>*−<sup>1</sup>) = **P***Tn*−1**P***<sup>n</sup>*−<sup>1</sup>−1**p***Tn* = **<sup>H</sup>**(*n* − <sup>1</sup>)−1**p***Tn* , *ρn*|(*<sup>n</sup>*−<sup>1</sup>) = **<sup>p</sup>***n***P***Tn*−1**P***<sup>n</sup>*−<sup>1</sup>−1**p***Tn* = **<sup>p</sup>***n***v***n*|(*<sup>n</sup>*−<sup>1</sup>). −<sup>1</sup>),

We expand Equation (7), and further simplify it by **Q**(*n* − 1) = **<sup>H</sup>**(*n* − <sup>1</sup>)−1**P**(*n* − <sup>1</sup>)*<sup>T</sup>***B**(*n* to ge<sup>t</sup>

$$\begin{split} \mathbf{Q}(n) &= \mathbf{H}(n-1)^{-1} \left( \mathbf{P}(n-1)^{T} \mathbf{B}(n-1) \right) + \widetilde{\mathbf{p}}\_{n}^{T} \mathbf{r}\_{n} \\ & \quad - \frac{\mathbf{v}\_{n|(n-1)} \mathbf{v}\_{n|(n-1)}^{T}}{1 + \rho\_{n|(n-1)}} \left( \mathbf{P}(n-1)^{T} \mathbf{B}(n-1) + \widetilde{\mathbf{p}}\_{n}^{T} \mathbf{r}\_{n} \right) \\ &= \mathbf{Q}(n-1) + \mathbf{H}(n-1)^{-1} \widetilde{\mathbf{p}}\_{n}^{T} \mathbf{r}\_{n} - \frac{\mathbf{v}\_{n|(n-1)} \mathbf{v}\_{n|(n-1)}^{T}}{1 + \rho\_{n|(n-1)}} \left( \mathbf{P}(n-1)^{T} \mathbf{B}(n-1) + \widetilde{\mathbf{p}}\_{n}^{T} \mathbf{r}\_{n} \right) . \end{split} \tag{8}$$

*Remote Sens.* **2018**, *10*, 367

Multiply the 3rd term of Equation (8) by **<sup>H</sup>**(*n* − <sup>1</sup>)**H**(*n* − 1)−1, and then re-organize it. Finally we obtain

$$\begin{split} \mathbf{Q}(n) &= \mathbf{Q}(n-1) - \frac{\mathbf{v}\_{n[(n-1)} \mathbf{v}\_{n[(n-1)}^{\mathrm{T}}]}{1 + \rho\_{n[(n-1)}} \mathbf{H}(n-1) \mathbf{Q}(n-1) \\ &\quad + \left( \mathbf{H}(n-1)^{-1} - \frac{\mathbf{v}\_{n[(n-1)} \mathbf{v}\_{n[(n-1)}^{\mathrm{T}}]}{1 + \rho\_{n[(n-1)}} \right)} \widetilde{\mathbf{p}}\_{n}^{\mathrm{T}} \mathbf{r}\_{n} \\ &= \left( \mathbf{I}\_{L \times L} - \frac{\mathbf{v}\_{n[(n-1)} \mathbf{v}\_{n[(n-1)}^{\mathrm{T}}]}{1 + \rho\_{n[(n-1)}} \mathbf{H}(n-1) \right) \mathbf{Q}(n-1) \\ &\quad + \left( \mathbf{H}(n-1)^{-1} - \frac{\mathbf{v}\_{n[(n-1)} \mathbf{v}\_{n[(n-1)}^{\mathrm{T}}]}{1 + \rho\_{n[(n-1)}} \right)} \widetilde{\mathbf{p}}\_{n}^{\mathrm{T}} \mathbf{r}\_{n} \end{split} \tag{9}$$

Equation (9) is the recursive version of Equation (2). In Equations (7)–(9), **Q**(*n* − 1) and **<sup>H</sup>**(*n* − 1) are the past information obtained at the (*n* − 1)th stage; **<sup>v</sup>***n*|(*<sup>n</sup>*−<sup>1</sup>) and *ρn*|(*<sup>n</sup>*−<sup>1</sup>) are the so-called innovation information [32], whose calculation is involved with both previous information and current information; and **p***n* and **r***n* are provided by the new incoming pixel. Once **Q**(*n*) is obtained, we calculate **<sup>H</sup>**(*n*) = **<sup>H</sup>**(*n* − 1) + **<sup>p</sup>***Tn***p***n* for the use in the next stage. It should be noted that Equation (9) can only be used when the updating condition is met. This detail will be presented in Section 3.2.

#### 3.1.2. Decomposition of the Residual Multiplication Term

On the other hand, in OMPBS, the other calculation that can be simplified is the residual multiplication term **<sup>R</sup>**(*n*)*Tj*+1**<sup>R</sup>**(*n*)*j*+<sup>1</sup> in Step 5 of OMPBS. Calculating this term would be time-consuming when *n* becomes large. According to OMPBS, the residual term is formed by

$$\mathbf{R}\_{j+1}(n) = \mathbf{B}(n) - \mathbf{P}\_j(n)\mathbf{Q}\_j(n). \tag{10}$$

We expand **<sup>R</sup>**(*n*)*Tj*+1**<sup>R</sup>**(*n*)*j*+<sup>1</sup> in terms of Equation (10):

$$\begin{split} \mathbf{R}\_{j+1}(n)^T \times \mathbf{R}\_{j+1}(n) &= \left(\mathbf{B}(n) - \mathbf{P}\_{j}(n)\mathbf{Q}\_{j}(n)\right)^T \left(\mathbf{B}(n) - \mathbf{P}\_{j}(n)\mathbf{Q}\_{j}(n)\right) \\ &= \mathbf{B}(n)^T \mathbf{B}(n) - \mathbf{Q}\_{j}(n)^T \mathbf{P}\_{j}(n)^T \mathbf{B}(n) - \mathbf{B}(n)^T \mathbf{P}\_{j}(n)\mathbf{Q}\_{j}(n) + \mathbf{Q}\_{j}(n)^T \mathbf{P}\_{j}(n)^T \mathbf{P}\_{j}(n) \mathbf{Q}\_{j}(n) . \end{split} \tag{11}$$

There are four terms required to be computed in Equation (11). The first term is **<sup>B</sup>**(*n*)*<sup>T</sup>***B**(*n*), which can be obtained by the following recursive equation

$$\mathbf{B}(n)^{T}\mathbf{B}(n) = \left[ \begin{array}{c} \mathbf{B}(n-1)^{T} \\ \end{array} \mathbf{r}\_{n}^{T} \right] \left[ \begin{array}{c} \mathbf{B}(n-1) \\ \end{array} \right] = \mathbf{B}(n-1)^{T}\mathbf{B}(n-1) + \mathbf{r}\_{n}^{T}\mathbf{r}\_{n}.\tag{12}$$

Similarly, the second term can be represented by

$$\mathbf{Q}\_{\rangle}(n)^{\mathsf{T}}\mathbf{P}\_{\rangle}(n)^{\mathsf{T}}\mathbf{B}(n) = \mathbf{Q}\_{\rangle}(n)^{\mathsf{T}} \left[ \begin{array}{c} \mathbf{P}\_{\rangle}(n-1)^{\mathsf{T}} & \overline{\mathbf{p}}\_{n}^{\mathsf{T}} \end{array} \right] \left[ \begin{array}{c} \mathbf{B}(n-1) \\ \mathbf{r}\_{n}^{\mathsf{T}} \end{array} \right] \\ = \mathbf{Q}\_{\rangle}(n)^{\mathsf{T}} \left( \mathbf{P}\_{\rangle}(n-1)^{\mathsf{T}}\mathbf{B}(n-1) + \overline{\mathbf{p}}\_{n}^{\mathsf{T}}\mathbf{r}\_{\mathsf{R}} \right), \tag{13}$$

where **<sup>Q</sup>***j*(*n*)*<sup>T</sup>* is supposed to be known by Equation (9). The third term **<sup>B</sup>**(*n*)*<sup>T</sup>***P***j*(*n*)**Q***j*(*n*) is the transposition of the second term, so it is unnecessary to re-compute it. Finally, the fourth term can be calculated by

$$\mathbf{Q}\_{\dot{j}}(n)^T \mathbf{P}\_{\dot{j}}(n)^T \mathbf{P}\_{\dot{j}}(n) \mathbf{Q}\_{\dot{j}}(n) = \mathbf{Q}\_{\dot{j}}(n)^T \mathbf{H}\_{\dot{j}}(n) \mathbf{Q}\_{\dot{j}}(n) = \mathbf{Q}\_{\dot{j}}(n)^T \left(\mathbf{H}\_{\dot{j}}(n-1) + \widetilde{\mathbf{p}}\_n^T \widetilde{\mathbf{p}}\_n\right) \mathbf{Q}\_{\dot{j}}(n) \tag{14}$$

In conclusion, term **<sup>R</sup>**(*n*)*Tj*+1**<sup>R</sup>**(*n*)*j*+<sup>1</sup> can be quickly calculated by using Equations (11)–(14).

#### *3.2. Condition for Applying Recursive Equations*

Unlike related work about PBP [33–35], where the recursive calculation can be applied to each iteration of PBP, Equation (9) of Re-OMPBS can only be implemented under particular conditions. To understand it, we assume that <sup>Ω</sup>*j*(*n*) presents the BS result at the *j*th iteration in the *n*th stage, and the result was just obtained. Now we move to calculate **Q***j*(*n*) for finding the residual. If the condition <sup>Ω</sup>*j*(*n*) = <sup>Ω</sup>*j*(*<sup>n</sup>* − 1) holds, we can use Equation (9) to calculate **Q***j*(*n*). More specifically, if the selected bands' indices of the *j*th iteration of the previous stage are the same as the band indices of the *j*th iteration in the current stage, we can use the recursive equation to save computing time. This is simply because the sparse coefficient matrix **Q***j*(*n*) can only be "updated" when the currently-selected bases (i.e., band indices) are exactly the same as the previous ones. More specifically, the information in the horizontal dimension of temporary matrices **<sup>P</sup>***j*(*<sup>n</sup>* − 1) and **<sup>P</sup>***j*(*n*) are aligned so that the information of previous stage can be shared with the current stage. In the following, we call it the "updating condition".

Satisfying the updating condition not only reduces computing time for **Q***j*(*n*), but also slightly accelerates Equations (12)–(14). One key of implementing Re-OMPBS is that we have to store the computing results obtained at the (*n* − 1)th stage, such as terms . **Q***j*−<sup>1</sup>(*n*) /*p j*=1 , . **<sup>H</sup>***j*−<sup>1</sup>(*n*) /*p j*=1 , **<sup>B</sup>**(*n* − 1) *<sup>T</sup>***B**(*n* − <sup>1</sup>), and , **<sup>P</sup>**(*n* − 1) *<sup>T</sup>***B**(*n* − 1) -*p j*=1, for applying in Equation (9) and Equations (11)–(14), when the updating condition is held at the *n*th stage. One should note that meeting this condition is not necessary for Re-OMPBS. For those cases that do not meet the conditions, we directly use the original least square formula **Q***j*(*n*) = **<sup>P</sup>***j*(*n*) *<sup>T</sup>***P***j*(*n*) −1 **<sup>P</sup>***j*(*n*) *<sup>T</sup>***B**(*n*) to calculate **Q***j*(*n*).

One might think that the probability of reaching these conditions would be low. However, this is not the case. In fact, when the amount of received pixels *n* reaches a certain amount, a single new incoming pixel will not obviously influence the BS result. In other words, the BS result is expected to be consistent. This phenomenon will be demonstrated in our experiment section. Moreover, Re-OMPBS is regarded as the "approximate" version of OMPBS. Theoretically, there would be some tiny numerical errors accumulated in **Q***j*(*n*) if the updating condition has been reached continuously for a long period. In this case, an occasional lack of meeting the updating condition exactly gives the chance to reset the numerical error to zero.
