Uncertainty Quantification of Hyperspectral Image Denoising Frameworks based on Sliding-Window Low-Rank Matrix Approximation

Fig. 1 demonstrates the routine process of the LMRA-based HSI restoration. The HSI is processed with a 3D sliding window to reduce the computational burden. For each 3D patch in the window, all the bands are vectorized and stacked into the 2D matrix. The permutated 2D matrix goes through the routine LMRA process, and the result is re-permutated to recover the 3D denoised patch. Finally, all the patches are summarized and averaged to create the restored HSI.

An HSI acquired by a sensor is typically modeled as a 3D image composed of Gaussian noise and random noise:

where $\mathbf{W}\in\mathbb{R}^{M\times N\times P}$ is the image acquired from the hyperspectral sensor with the size $M\times N$ and $P$ channels, $\mathbf{W}^{\star}\in\mathbb{R}^{M\times N\times P}$ is the noise-free image denoted as the ground truth, $\mathbf{G}\in\mathbb{R}^{M\times N\times P}$ is the element-wise i.i.d. Gaussian noise and $\mathbf{N}\in\mathbb{R}^{M\times N\times P}$ is the sparse and random noise (impulse signal). The random noise $\mathbf{N}$ is normally considered negligible with regard to $\mathbf{G}$, and the data quality quantification process is inferring the statistical distribution of $\mathbf{G}$. In practice, the ground truth HSI $\mathbf{W}^{\star}$ cannot be obtained. Thus, we denote $\hat{\mathbf{W}}$ as the estimated denoised HSI calculated from studies based on the low-rank assumption [21, 22].

The goal of HSI denoising is to recover the denoised $\hat{\mathbf{W}}$ with the 3D HSI noisy observation $\mathbf{W}$ and the exact rank prior. Assuming that the random noise $\mathbf{N}$ is negligible and the residual follows the i.i.d. zero-mean Gaussian distribution, the restored image $\hat{\mathbf{W}}$ is recovered with the following formulation:

The idea of Eq. (LABEL:Eq_image_recovery) is to convert the 3D patch into a typical 2D small-size matrix so that the conventional low-rank factorization algorithm can be executed with limited computational resources. The equation is often solved with a Lagrange multiplier or an augmented Lagrange multiplier (ALM); off-the-shelf low-rank solvers can be used. In this work, we choose the GoDec algorithm proposed by [34] to be consistent with the LRMA method [20]. However, it should be emphasized that the proposed uncertainty quantification method is independent of the solver since it only assumes that the estimation of the LRMA is close to the minimum. The better the performance of the LRMA solver, the more accurate the confidence interval of our method describes.

After all small overlapped patches have gone through the process of Eq. (LABEL:Eq_image_recovery), the optimized denoised patches are summed to restore the HSI:

Fig. 2 illustrates the general framework and the relationship with the routine LRMA approach. In the LRMA process, the temporary variable, i.e., the 2D sliding-window patch, is utilized to provide the uncertainty to the HSI patch based on our closed-form formulation. After the HSI patch and its associated closed-form uncertainty have been generated, the routine LRMA approach averages them to yield the restored HSI. We also provide a closed-form uncertainty propagation to quantify the final restored HSI. As shown in Fig. 2, in all processes, the proposed approach is independent of the choice of the low-rank solver. The method can also be easily integrated into the derivations of the LRMA process. In addition, since our method is a closed-form solution and only utilizes temporary variables, the processing time is relatively short compared to the routine LRMA process.

One major contribution of this work is that it provides an element-wise closed-form solution to quantify the uncertainty in a time-efficient manner. For the restored HSI $\hat{\mathbf{W}}$, one straightforward approach of the element-wise uncertainty quantification is conducting the Monte Carlo test with the given statistical distribution of the observation $\mathbf{W}$. However, the size of HSIs is notoriously large, resulting in a heavy computational burden to conduct numerous Monte Carlo tests. Roughly, a one-time process in [20] on a $100\times 100\times 100$ HSI requires 20 minutes with a typical commercial desktop computer. Since a robust uncertainty quantification requires a reasonable number of trials, the Monte Carlo test is inefficient for use in the element-wise HSI uncertainty quantification.

Another benefit of this research is that it solves the uncertainty propagation of the sliding-window process. The uniqueness of HSIs is that the permutated 2D matrix after dimensionality reduction is still relatively large. The solution cannot be achieved with a one-time LRMA process. Instead, researchers have adopted the sliding-window strategy to segment the HSI into many small patches and process them parallel. The different statistical properties of each overlapped patch make it challenging to quantify the uncertainty of the final restored HSI. We solve this problem and present a closed-form solution in the proposed sliding-window element-wise closed-form HSI uncertainty quantification.

This section focuses on solving the uncertainty quantification of the patch-wise low-rank factorization process (Eq. (LABEL:Eq_image_recovery)). We follow the example of [31] and draw a similar conclusion for the uncertainty quantification of the exact low-rank patch factorization. A concise and clear proof is also provided. HSI patch-based averaging (Eq. (5)) will be discussed in the next subsection. A recent study [31] proposed an uncertainty quantification method for the exact low-rank 2D matrix. This element-wise algorithm can quantify the confidence interval for a common low-rank based 2D matrix denoising algorithm. Synthetic and real-world tests validate the efficiency of this method. We introduce this approach and modify it to enable the quantification of the uncertainty of each estimated 3D HSI patch $\hat{\mathbf{W}}^{(ijk)}$ using Eq. (LABEL:Eq_image_recovery). We also provide concise proof different from [31].

In the patch optimization process shown in Eq. (LABEL:Eq_image_recovery), we first aim to quantify the uncertainty of the temporary 2D square matrix $\hat{\mathbf{W}}^{\sharp{(ijk)}}$ ($\hat{\mathbf{W}}^{\sharp{(ijk)}}=g(\hat{\mathbf{W}}^{{(ijk)}})$ and it can be any patch in $\Omega$). We denote the ground truth of $\hat{\mathbf{W}}^{\sharp{(ijk)}}$ and $\hat{\mathbf{W}}^{{(ijk)}}$ as ${\mathbf{W}}^{\sharp\star{(ijk)}}$ and ${\mathbf{W}}^{\star{(ijk)}}$. For the sake of conciseness, the rank-r SVD decomposition is ${\mathbf{W}}^{\sharp\star{(ijk)}}=\mathbf{U}^{\star(ijk)}\Sigma^{\star(ijk)}\mathbf{V}^{\star(ijk)\top}$. For conciseness, we simplify $(\mathbf{U}^{\star(ijk)},\Sigma^{\star(ijk)},\mathbf{V}^{\star(ijk)})$ by ignoring the patch index $(ijk)$ and rewrite it as $(\mathbf{U}^{\star},\Sigma^{\star},\mathbf{V}^{\star})$. Two matrices are defined as the product of the unitary matrix and the squared singular matrix is: $\mathbf{X}^{\star}\triangleq\mathbf{U}^{\star}\Sigma^{\star 1/2}\in\mathbb{R}^{K\times r}$, $\mathbf{Y}^{\star}\triangleq\mathbf{V}^{\star}\Sigma^{\star 1/2}\in\mathbb{R}^{L\times r}$. The following rules apply:

Following the above definition, the ground truth 2D temporary matrix is $\mathbf{W}^{\sharp\star{(ijk)}}=\mathbf{X}^{\star}\mathbf{Y}^{\star\top}$. Similarly, we define the estimated version $\hat{\mathbf{W}}^{\sharp{(ijk)}}\approx\hat{\mathbf{X}}\hat{\mathbf{Y}}^{\top}$ where $\hat{\mathbf{X}}$ and $\hat{\mathbf{Y}}$ are the estimated matrices regarding $\mathbf{X}^{\star}$ and $\mathbf{Y}^{\star}$. We also obtain the estimated version of $\Sigma^{\star}$ as $\hat{\Sigma}$:

It should be noticed that the factorization of $\hat{\mathbf{W}}$ is not unique. For each pair of estimated auxiliary matrices ($\hat{\mathbf{X}},\hat{\mathbf{Y}}$), there are infinite unitary matrices $\mathbf{T}$ that satisfy:

Among the unitary matrices, there is a rectification matrix $\hat{\mathbf{R}}$ that minimizes the distance between the auxilliary matrices ($\hat{\mathbf{X}},\hat{\mathbf{Y}}$) and the ground truth (${\mathbf{X}}^{\star},{\mathbf{Y}}^{\star}$):

We provide the proof of theorem 1 in Section A. This proof is specifically designed for the full HSI patch denoising scenario and is clearer and more concise than the one in [31].

In theorem 1, we make two assumptions, i.e., the Gaussian noise distribution and the global minimum obtained from the low-rank solver. The former ensures the correctness of the error propagation equations. The latter ensures that the residual noises $\Psi_{X}$ and $\Psi_{Y}$ are negligible compared to $\mathbf{Z}_{\mathbf{X}}$ and $\mathbf{Z}_{\mathbf{Y}}$, and Eq. (10) is approximated to:

Assuming that the first-order expansion is reasonably tight and omitting the high order expansion, the element-wise error between the estimation $\hat{\mathbf{W}}_{uv}^{\sharp{(ijk)}}$ and the ground truth ${\mathbf{W}}_{uv}^{\sharp\star{(ijk)}}$ at the position $uv$ is:

In practice, the rows of the ground truth $\mathbf{U}_{u,.}^{\star}$ and $\mathbf{V}_{v,.}^{\star}$ are impossible to obtain. Assuming that the solver achieves a solution close to the global minimum, we approximate them with the estimated version $\hat{\mathbf{U}}_{u,.}$ and $\hat{\mathbf{V}}_{v,.}$, which are the $u$/$v$th row of the estimations $\hat{\mathbf{U}}$ and $\hat{\mathbf{V}}$. Thus the element-wise variance can be approximated as:

We provide a toy example as the weighted average of events $X$ and $Y$ as $F=\mathrm{a}X+\mathrm{b}Y$. Assuming the variance of $X$ and $Y$ are $\sigma_{X}^{2}$ and $\sigma_{Y}^{2}$, the propagated uncertainty of $F$ is $\sigma_{F}^{2}$:

For convenience, we rewrite Eq. (5) to the element-wise form:

We provide a practical solution to estimate the optimal ratio $\eta^{(ijk)(i^{\prime}j^{\prime}k^{\prime})}$. Similar to theorem 1, we assume that the adopted low-rank solver (GoDec algorithm in this paper, but it could be any solver) achieves the global minimum, the information on the denoised element $\xi$ in the patch $\hat{\mathbf{W}}^{(ijk)}$ is determined by the sum of the information embedded in all the elements in the patch $\hat{\mathbf{W}}^{(ijk)}$. The mutual information encoded in two overlapping patches is linearly related to the number of elements shared by the overlapping patches [36]. Thus the mutual information of the element from two patches is linearly related to the number of shared elements in the two patches. Therefore, the optimal ratio between the two patches $\eta^{(ijk)(i^{\prime}j^{\prime}k^{\prime})}$ is close to:

Fig. 3 presents an example of Eq. (21). Each neighboring block pair has $50\%$ overlap. Based on the statistics of the number of overlapping elements, the $\eta$ of the pairs $(1,2)$, $(2,3)$, $(3,4)$ and $(4,1)$ is 0.5, whereas the $\eta$ of pairs $(1,3)$ and $(2,4)$ is 0.25.

The result of Eq. (16) is based on the assumption that the noises $\Psi_{X}$ and $\Psi_{Y}$ are significantly smaller than the variables $\mathbf{Z}_{\mathbf{X}}$ and $\mathbf{Z}_{\mathbf{Y}}$ with a normal distribution. If $\Psi_{X}$ and $\Psi_{Y}$ are not ignored, the conclusion does not hold. Therefore, it remains challenging to quantify the impact of random noise versus Gaussian noise. However, as [33] pointed out, the RPCA-like formulation (Eq. (LABEL:Eq_image_recovery).) is robust against random noise since the low-rank solver will adjust the use of the norm-1 or norm-2 to regularize the state matrix. Correspondingly, elements not corrupted by random noise are unaffected, whereas the impulse noise can be identified. Therefore, regarding the impact of random noise, the robustness of our approach is closely related to the choice of the LRMA algorithm. It should be noted that most existing HSI denoising studies have focused on modeling Gaussian noise [21, 25]. A few LRMA-based methods have addressed non-Gaussian noise [24]. These approaches utilize a total-variance related term to filter out non-Gaussian noise in the HSI to maintain good boundaries after filtering. Similarly, the proposed uncertainty estimation method can also be integrated with the method in [24] since the total-variance term is only a constraint.

The accuracy of the proposed closed-form uncertainty quantification method is extensively validated with numerical Monte Carlo tests. We apply different levels of Gaussian noises to the HSI. The numerical experiments provide sufficient samples to test if the statistical distributions of the samples are described accurately with the output confidence interval. Specifically, $\mathrm{T}=100$ Monte Carlo tests are performed on the input of the same noise $\sigma_{0}^{2}$. We should emphasize that the 100 trials are chosen due to limited computational resources. Due to the large volume of the HSIs, 100 Monte Carlo trials are sufficient to obtain a statistic assessment of the efficiency of the proposed method. The code is provided for readers to conduct more tests. For each element of the estimated HSI in trial $l$ at position $\xi$ as $\hat{\mathbf{W}}_{\xi}^{\{l\}}$, we define its deviation from the corresponding $\mathrm{T}$ times mean value as $\overline{\hat{\mathbf{W}}}_{\xi}$. The deviation from the mean is defined as:

If the closed-form solution describes the statistical distribution of the estimation well, $95\%$ of $\hat{\mathbf{W}}_{\xi}^{\{l\}}$ should fall within the range of $\left[\overline{\hat{\mathbf{W}}}_{\xi}-1.96\hat{\sigma}_{\xi},\ \overline{\hat{\mathbf{W}}}_{\xi}+1.96\hat{\sigma}_{\xi}\right]$. Following Eq. (22), we define the ratio of $\mathrm{e}^{\{l\}}_{\xi}$ falling within the $1.96\hat{\sigma}_{\xi}$ bound as the element-wise coverage rate. The average coverage rate of all elements is the average coverage rate. Table I presents the average and standard deviation of the coverage rate. The result indicates that the uncertainty approach robustly and accurately describes the element-wise denoised HSI when the level of noise is low ($\sigma_{0}$ less than $0.100$) ⁴⁴4https://JingweiSong@bitbucket.org/JingweiSong/uncertainty-estimation-for-hyper-spectral-image-denoising.git.. In the case $\sigma_{0}\geq 0.100$, we notice an under-confident of the prediction. We hypothesize that the accuracy is related to the level of noise and the choice of the optimal rank. When the noise level increases, the fixed rank LRMA does not accurately filter the noise while retaining useful information. A lower rank is needed to ensure that the residual follows the assumption of the normal distribution of the noise. As depicted in Table I, we strictly follow [20] by setting $\mathrm{r}=7$. To determine the reason, we conduct an experiment by choosing different $\mathrm{r}$. The results are shown in Fig. (6), indicating that the performance (uncertainty quantification) of our approach improves by decreasing the arbitrary rank. This implies that the Monte Carlo experiment provides a method to validate the correct choice of the rank. Further studies are needed to determine the rank at high noise level.

In addition to the average coverage rate, we also adopt the Monte Carlo experiments to prove that the restored HSI $\hat{\mathbf{W}}^{\{l\}}$ has an innate Gaussian distribution. Quantile-Quantile (Q-Q) plots are created, and Shapiro-Wilk tests are conducted to validate the proposed approach. Fig. 5 shows the Q-Q plots of all the datasets. Note that we randomly plot one element in each dataset as an example, but all other elements are similar since the observation is isotropic. Fig. 5 indicates that the Monte Carlo estimation of all the trials $\hat{\mathbf{W}}^{\{l\}}$ is close to a line; thus, the sample follows the Gaussian distribution.

In the Shapiro-Wilk tests, it is not possible to present all the elements. Thus, we randomly choose several elements from the three datasets with different levels of Gaussian noise. The results are presented in Table II. The $p$-values indicate that all samples follow the normal distribution ($>0.05$). Since all pixels are acquired under the same conditions, and their statistics should be consistent (either obey or reject the normal distribution assumption), we can conclude that the restored HSI still follows the Gaussian distribution.

In addition to the accuracy tests, we determine the time consumption. We find that the closed-form uncertainty estimation method has similar time consumption as the LRMA and far less time consumption than the Monte Carlo numerical tests. We use the LRMA code used in [20] and add the uncertainty estimation module. The proposed algorithm, the LRMA, and the Monte Carlo tests are all implemented in Matlab 2019a. In the Monte Carlo tests, we utilize the Parallel Computing Toolbox of Matlab to speed up the steps. The hardware set is $3.1GHz$ Intel Xeon E5-2687wv3 and 64 GB DDR4-Ram. In the Monte Carlo experiments, all 20 CPU threads are utilized. Note that in the Monte Carlo tests, each element is covered by several sliding windows. Since the memory is not sufficient to store all sliding windows, we calculate the result of the patch, store it on the hard disk, and reload it to calculate each element in the final step.

Table III shows the time consumption of the Monte Carlo uncertainty estimation, the LRMA method [20], and the proposed method integrated into [20]. The 100 Monte Carlo trials are significantly speeded up because of the running in parallel using 20 CPU threads. The results indicate that the proposed closed-form estimation does not add a significant amount of time to the LRMA method, whereas the Monte Carlo test requires substantially more time. As indicated by Eq. (14), the proposed approach requires no more than one band-wise SVD decomposition of the given HSI in addition to the original LRMA solver.

In addition to the coverage rate tests, we assess if the proposed approach is robust to the random noise from the sensor. As Eq. (2) shows, our work ignores the random noise $\mathbf{N}$ and only assumes that the residuals follow the Gaussian noise $\mathbf{G}$. As described in Section III-D, we expect that the proposed method is robust to random noise such as impulse noise, dead lines, and stripes [37, 38]. We conduct the tests by deliberately corrupting the original HSI with random noise and test the performance of the proposed uncertainty quantification method. We simulate different proportions of impulse noise on all datasets and determine the coverage rate of the Monte Carlo experiments.

Table IV shows the results of the Monte Carlo experiments for different coverage rates. The results indicate a relationship between the robustness of the proposed closed-form method and the coverage rate. When the noise $\sigma_{0}$ is greater than $0.1$, the uncertainty quantification approach is less sensitive to pulse noises. When the noise $\sigma_{0}$ is less than $0.075$, the proposed method is susceptible to random noise. Theoretically, the proposed approach is only effective for Gaussian noise, which is the low-rank prior to Eq. (2). Therefore, less Gaussian noise means less reliable of the assumption that $\mathbf{N}$ is less than $\mathbf{G}$ in Eq. (1). Thus, the results from the closed-form approach are less reliable. When the Gaussian noise level of the observation is low, the proposed approach tends to be over-confident of the result, and the result is lower than the coverage rate with a higher $\sigma_{0}$. A higher $\sigma_{0}$ means that the distribution is closer to the impulse noise. Furthermore, since the solvers (either GoDec or ALM solver) are all based on SVD decomposition, our approach has the same computational complexity as the original LRMA methods.