## 4 Solution
### Solving the OSG-Lp Model
\(C(u)\) in the above equation is independent of \(v\), and \(\mathbf{D}(u)\in\mathbb{R}^{M^{2}\times M^{2}}\) is a diagonal matrix with its elements defined as  
\[\left[\mathbf{D}(u)\right]_{m,m}=\sqrt{\sum_{i=-N_{i}}^{N_{r}}\sum_{j=-N_{l}}^{N_ {r}}\left[\sum_{k_{1}=-N_{l}}^{N_{r}}\sum_{k_{2}=-N_{l}}^{N_{r}}\left|u_{m-i+k _{1},m-i+k_{2}}\right|^{2}\right]^{-\frac{1}{2}}}\qquad\left(m=1,2,\cdots,M^{ 2}\right), \tag{22}\]  
The entries of \(\mathbf{D}\) can be easily computed by using MATLAB built-in function \"conv2\". Putting Equations (17), (21), and (22) together, the optimization problem \(P(v)\) can be written as  
\[Q(v,u)=\frac{\alpha}{2}\left\|v-v_{0}\right\|_{2}^{2}+S(v,u)=\frac{\alpha}{2} \left\|v-v_{0}\right\|_{2}^{2}+\frac{1}{2}\left\|\mathbf{D}(u)v\right\|_{2}^{2}+C( u). \tag{23}\]  
When \(v=u\), \(Q(u,u)=P(u)\), to minimize \(P(v)\), the MM aims to iteratively solve  
\[v^{(n+1)}=\underset{v}{\text{argmin}}\frac{\alpha}{2}\left\|v-v_{0}\right\|_{ 2}^{2}+\frac{1}{2}\left\|\mathbf{D}\!\left(v^{(n)}\right)v\right\|_{2}^{2},n=1,2,3\ldots \tag{24}\]  
with the solution  
\[v^{(n+1)}=\left(\mathbf{I}+\frac{1}{\alpha}\mathbf{D}^{2}\!\left(v^{(n)}\right)\right) ^{-1}\!\!v_{0},n=1,2,3\ldots, \tag{25}\]  
where \(\mathbf{I}\in\mathbb{R}^{M^{2}\times M^{2}}\) is an identity matrix with the same size of \(\mathbf{D}^{2}\!\left(v^{(n)}\right)\), \(\bm(2) \(\mathbf{Z}_{3}^{(k+1)}\) could be solved by the famous soft thresholding shrinkage method [45] as shown below:  
\[ \begin{split}&\mathbf{Z}_{3}^{(k+1)}=\underset{\mathbf{Z}_{3}}{\text{ argmin}}\|\mathbf{Z}_{3}\|_{p}^{p}+\frac{\delta_{3}}{2}\|\mathbf{Z}_{3}-\mathbf{F}^{(k)}+\mathbf{G}- \widetilde{\mathbf{Z}}_{3}^{(k)}\|_{2}^{2}\\ &=\text{\emph{shrinkage}}(\mathbf{F}^{(k)}-\mathbf{G}+\widetilde{\mathbf{Z}}_ {3}^{(k)},\frac{1}{\beta_{3}},p)\\ &=\text{\emph{sign}}\Big{(}\mathbf{F}^{(k)}-\mathbf{G}+\widetilde{\mathbf{Z} }_{3}^{(k)}\Big{)}\\ &\cdot*\max\Big{(}\text{\emph{abs}}\Big{(}\mathbf{F}^{(k)}-\mathbf{G}+ \widetilde{\mathbf{Z}}_{3}^{(k)}\Big{)}-\Big{(}\frac{1}{\beta_{3}}\Big{)}\gamma(p -2).*\text{\emph{abs}}\Big{(}\mathbf{F}^{(k)}-\mathbf{G}+\widetilde{\mathbf{Z}}_{3}^{(k)} \Big{)}\gamma(p-1),0\Big{)}.\end{split}  \tag{28}\]  
(3) The \(\mathbf{Z}_{4}^{(k+1)}\) subproblem is computed as  
\[\mathbf{Z}_{4}^{(k+1)}=\underset{\mathbf{Z}_{4}}{\text{argmin}}\frac{\beta_{4}}{2}\| \mathbf{Z}_{4}-\mathbf{F}^{(k)}-\widetilde{\mathbf{Z}}_{4}^{(k)}\|_{2}^{2}=\min\Big{(}1, \max\Big{(}\mathbf{F}^{(k)}+\widetilde{\mathbf{Z}}_{4}^{(k)},0\Big{)}\Big{)},. \tag{29}\]  
(4) The \(\mathbf{F}^{(k+1)}\) subproblem is solved by substituting \(\mathbf{Z}_{1}^{(k+1)},\mathbf{Z}_{2}^{(k+1)},\mathbf{Z}_{3}^{(k+1)},\mathbf{Z}_{4}^{(k+1)}\) into Equation (12) and calculating for the variable \(\mathbf{F}^{(k+1)}\). Under the assumption of periodic boundary conditions, fast Fourier transform is applied to both sides of the equation to perform the computation in the frequency domain instead of a spatial domain, in order to reduce the computational complexity caused by matrix multiplication. Matrix multiplication is converted to dot product operation, in other words, solving the following normal equation:  
\[ \begin{split}&\Big{(}\beta_{1}\overline{\mathbf{K}}_{\mathbf{h}}^{*}* \overline{\mathbf{K}}_{\mathbf{h}}+\beta_{2}\overline{\mathbf{K}}_{\mathbf{v}}^{*}*\overline{ \mathbf{K}}_{\mathbf{v}}+\beta_{3}\mathbf{1}+\beta_{4}\mathbf{1}\Big{)}.*\overline{\mathbf{F}}^{(k +1)}\\ &=\beta_{1}\overline{\mathbf{K}}_{\mathbf{h}}^{*}*\left(\overline{\mathbf{Z}}_ {1}^{(k+1)}-\overline{\mathbf{Z}}_{1}^{(k)}\right)+\beta_{2}\overline{\mathbf{K}}_{ \mathbf{v}}^{*}*\left(\overline{\mathbf{Z}}_{2}^{(k+1)}-\overline{\mathbf{Z}}_{2}^{(k)} \right)\\ &+\beta_{3}\left(\overline{\mathbf{Z}}_{3}^{(k+1)}+\overline{\mathbf{G}}- \overline{\mathbf{Z}}_{3}^{(k)}\right)+\beta_{4}\left(\overline{\mathbf{Z}}_{4}^{(k+ 1)}-\overline{\mathbf{Z}}_{4}^{(k)}\right),\end{split}