# A Probabilistic Hyperspectral Imagery Restoration Method
## 4 Probabilistic LRMA Model for HSI Restoration
### Clean HSI Estimation Method  
Considering matrix decomposition has the advantage to represent a low-rank matrix both in speed and scalability, it can be adopted to represent \(\mathbf{L}\). Specifically,  
\[\mathbf{L}=\mathbf{M}\mathbf{H}^{\text{T}}=\sum\
olimits_{r=1}^{R}\mathbf{m} _{.r}^{\text{T}}. \tag{10}\]  
The \(r\)-th column of \(\mathbf{M}\) and \(\mathbf{H}\) is denoted as \(\mathbf{m}_{.r}\) and \(\mathbf{h}_{.r}\), respectively. \(\theta=(\theta_{1},\theta_{2},\cdots,\theta_{R})\) is precision parameter. To guarantee the generated matrix \(\mathbf{L}\) a low-rank matrix, \(\mathbf{M}\) and \(\mathbf{H}\) can be probabilistically represented by Gaussian distributions [42] as  
\[ \begin{split} p(\mathbf{M}|\theta)&\sim\prod \limits_{r=1}^{\mathbf{R}}N(\mathbf{m}_{.r}|0,\,\theta_{r}^{-1}\mathbf{I}_{n_{ b}})\\ p(\mathbf{H}|\theta)&\sim\prod\limits_{r=1}^{ \mathbf{R}}N(\mathbf{h}_{.r}|0,\,\theta_{r}^{-1}\mathbf{I}_{n_{p}})\end{split}  \tag{11}\]  
where \(\mathbf{I}_{n_{p}}\) is an identity matrix in the size of \(n_{p}\times n_{p}\). \(\theta_{r}\) obeys Gamma distribution with hyperparameters \(a_{0}\) and \(b_{0}\), viz. \(p(\theta_{r})=\text{Gamma}(\theta_{r}|a_{0},b_{0})\). With such a matrix decomposition \(\mathbf{L}=\mathbf{M}\mathbf{H}^{\text{T}}\), \(p(\mathbf{X}|\mathbf{L},\mathbf{\omega},\mathbf{\mu},\mathbf{\tau})\) in Equation (9) turns to be  
\[p(\mathbf{X}|\mathbf{M},\mathbf{H},\mathbf{\omega},\mathbf{\mu},\mathbf{\tau})=\prod \limits_{(i,j)\in\Omega}\Bigl{(}\sum\
olimits_{k=1}^{K}\omega_{k}N(\mathbf{x} _{ij}|\mathbf{m}_{i}\mathbf{h}_{j.}^{\text{T}}+\mu_{k},\tau_{k}^{-1})\Bigr{)}, \tag{12}\]  
and \(p(\mathbf{L},\mathbf{\omega},\mathbf{\mu},\mathbf{\tau},\theta|\mathbf{X})\) turns to be \(p(\mathbf{M},\mathbf{H},\mathbf{\omega},\mathbf{\mu},\mathbf{\tau},\theta|\mathbf{X})\). It means instead of estimating \(\mathbf{L}\) directly, we need to estimate \(\mathbf{M}\) and \(\mathbf{H}\) first, and then obtain \(\mathbf{L}\) using \(\mathbf{L}=\mathbf{M}\mathbf{H}^{\text{T}}\). According to the Bayesian rule, \(p(\mathbf{M},\mathbf{H},\mathbf{\omega},\mathbf{\mu},\mathbf{\tau},\theta|\mathbf{X})\) is proportional to \(p(\mathbf{M},\mathbf{H},\mathbf{\omega},\mathbf{\mu},\mathbf{\tau},\theta,\mathbf{X})\), which can be further factorized into some simple distributions to guarantee the easier estimation from the perspective of generative model. The above analysis on \(p(\mathbf{M},\mathbf{H},\mathbf{\omega},\mathbf{\mu},\mathbf{\tau},\theta|\mathbf{X})\) can be represented by  
\[ \begin{split} p(\mathbf{M},\mathbf{H},\mathbf{\omega},\mathbf{\mu},\mathbf{ \tau},\theta|\mathbf{X})\propto& p(\mathbf{M},\mathbf{H},\mathbf{ \omega},\mathbf{\mu},\mathbf{\tau},\theta,\mathbf{X})\\ &=p(\mathbf{X}|\mathbf{M},\mathbf{H},\mathbf{\omega},\mathbf{\mu},\mathbf{ \tau})p(\mathbf{\omega})p(\mathbf{\mu}|\mathbf{\tau})p(\mathbf{\tau})p(\mathbf{M}| \theta)p(\mathbf{H}|\theta)p(\mathbf{\Theta})\end{split}  \tag{13}\]  
Though the factorized items such as \(p(\mathbf{M}|\theta)\) and \(p(\mathbf{H}|\theta)\), etc., can facilitate the estimation of unknown parameters \(\{\mathbf{M},\mathbf{H},\mathbf{\omega},\mathbf{\mu},\mathbf{\tau},\theta\}\), directly estimating \(\mathbf{M}\) and \(\mathbf{H}\) using \(p(\mathbf{M},\mathbf{H},\mathbf{\omega},\mathbf{\mu},\mathbf{\tau},\theta,\mathbf{X})\) is still intractable. To address this problem effectively, variational method or sampling-based method can be used to calculate \(\mathbf{M}\) and \(\mathbf{H}\). In this study, the variational method proposed in [42] is used for this purpose. The basic idea of this method is given as follows.