## 3 Proposed Method
By setting the first-order derivative of \(\mathbf{Z}\) in Equation (30) as zero, we have  
\[\mathbf{0}=\mathcal{F}(\mathbf{Z})-\mathcal{F}(\widehat{\mathbf{Z}})+\eta\sum_{i=1}^{4}\{[ \mathcal{F}(\mathbf{K}_{i})]^{*}\circ[\mathcal{F}(\mathbf{K}_{i})\circ\mathcal{F}(\bm {Z})-\mathcal{F}(\mathbf{T}_{i}^{(n)})-\frac{\mathcal{F}(\mathbf{V}_{i}^{(n)})}{\eta}]\}, \tag{31}\]  
where \([\mathcal{F}(\mathbf{K}_{i})]^{*}\) is the conjugate map of \(\mathcal{F}(\mathbf{K}_{i})\). Then, \(\mathbf{Z}\) can be obtained as follows:  
\[\mathbf{Z}^{(n+1)}=\mathcal{F}^{-1}\frac{\mathcal{F}(\widehat{\mathbf{Z}})+\eta\sum_{i= 1}^{4}\{[\mathcal{F}(\mathbf{K}_{i})]^{*}\circ[\mathcal{F}(\mathbf{T}_{i}^{(n)})+ \frac{\mathcal{F}(\mathbf{V}_{i}^{(n)})}{\eta}]\}}{\mathbf{I}+\eta\sum_{i=1}^{4}[ \mathcal{F}(\mathbf{K}_{i})]^{*}\circ\mathcal{F}(\mathbf{K}_{i})}, \tag{32}\]  
where \(\mathcal{F}^{-1}\) represents the 2D inverse FFT operator; the division in the Equation (32) is component-wise as well. In addition, the symbol \(\mathbf{I}\) represents the matrix with all entries being \(1\). The sub-problem \(\mathbf{T}\) is  
\[\mathbf{T}_{i}^{(n+1)}=\arg\min_{\mathbf{T}_{i}}\lambda\|\mathbf{T}_{i}\|_{p}^{p}+\frac{\eta} {2}\|\mathbf{T}_{i}^{(n)}-\mathbf{K}_{i}*\mathbf{Z}+\frac{V_{i}^{(n)}}{\eta}\|_{2}^{2}. \tag{33}\]  
It can be solved using the soft threshold operator [35] as follows:  
\[\mathbf{T}_{i}^{(n+1)}=shrink_{p}(\mathbf{K}_{i}*\mathbf{Z}^{(n+1)}-\frac{\mathbf{V}_{i}^{(n)} }{\eta},\frac{\lambda}{\eta}), \tag{34}\]  
where  
\[shrink_{p}(\xi,\frac{1}{\beta})=max\{|\xi|-\beta^{p-2}|\xi|^{p-1},0\}\cdot\frac {\xi}{|\xi|}. \tag{35}\]  
Lastly, the Lagrange multiplier can be updated as  
\[\mathbf{V}_{i}^{(n+1)}=\mathbf{V}_{i}^{(n)}+\eta(\mathbf{T}_{i}^{(n+1)}-\mathbf{K}_{i}*\mathbf{Z}^ {(n+1)}). \tag{36}\]  
Similar to Algorithm 2, we introduce the acceleration step \(b_{i}^{(n+1)}\), variables \(\mathbf{\tilde{T}}_{i}^{(n+1)}\), and dual variables \(\mathbf{\tilde{V}}_{i}^{(n+1)}\):  
\[\mathbf{\Upsilon}_{i}^{(n+1)}=\mathbf{T}_{i}^{(n+1)}+\frac{b_{i}^{(n)}-1}{b_{i}^{(n+1 )}}(\mathbf{T}_{i}^{(n+1)}-\mathbf{T}_{i}^{(n)}), \tag{37}\]  
\[\mathbf{\tilde{V}}_{i}^{(n+1)}=\mathbf{V}_{i}^{(n+1)}+\frac{b_{i}^{(n)}-1}{b_{i}^{(n+1 )}}(\mathbf{V}_{i}^{(n+1)}-\mathbf{V}_{i}^{(n)}). \tag{38}\]  
In this manner, all the sub-problems of Equation (22) are solved independently. In all iterations, the sub-problem of \(\mathbf{Z}\) is the most time-consuming problem. Considering the special structure of the differential matrices, we regard the differential operators as convolution operators. By introducing the convolution theorem, the sub-problem is skillfully solved in the frequency domain. The entire algorithm to solve Equation (22) is summarized in Algorithm 3. ```
Initialize:\(b_{i}^{(0)},\ \eta,\ \varphi\in(0,1)\quad(i=1,2,3,4)\) While\(\|\mathrm{PSNR}(\mathbf{Z}^{(n+1)},\mathbf{\tilde{Z}})-\mathrm{PSNR}(\mathbf{Z}^{(n)},\mathbf{ \tilde{Z}})\|_{2}/\|\mathrm{PSNR}(\mathbf{Z}^{(n+1)},\mathbf{\tilde{Z}})\|_{2}>\text{tol}\)do
1. \(\mathbf{Z}^{(n+1)}=\mathcal{F}^{-1}\frac{\mathcal{F}(\tilde{Z})+\eta\sum_{i=1}^{4 }\{[\mathcal{F}(\mathbf{K}_{i})]^{*}\diamond[\mathcal{F}(\mathbf{T}_{i}^{(n)})+\frac{ \mathcal{F}(\tilde{V}_{i}^{(n)})}{\eta}]\}}{I+\sum_{i=1}^{4}[\mathcal{F}(\mathbf{K }_{i})]^{*}\diamond\mathcal{F}(\mathbf{K}_{i})}\)
2. \(\mathbf{T}_{i}^{(n+1)}=shrink_{p}(\mathbf{K}_{i}*\mathbf{Z}^{(n+1)}-\frac{\mathbf{\tilde{V}}_{ i}^{(n)}}{\eta},\frac{\lambda}{\eta})\)
3. \(\mathbf{V}_{i}^{(n+1)}=\mathbf{\tilde{V}}_{i}^{(n)}+\eta\left(\mathbf{T}_{i}^{(n+1)}-\mathbf{K} _{i}*\mathbf{Z}^{(n+1)}\right)\)
4. \(d_{i}^{(n+1)}=\eta^{-1}\|\mathbf{V}_{i}^{(n+1)}-\mathbf{\tilde{V}}_{i}^{(n+1)}\|^{2}+ \eta\|\mathbf{T}_{i}^{(n+1)}-\mathbf{\tilde{T}}_{i}^{(n+1)}\|^{2}\)
5. if\(d_{i}^{(i+1)}<\varphi d_{i}^{(n)},\quad i=1,2,3,4,\quad\text{then}\)
6. \(b_{i}^{(n+1)}=\frac{1+\sqrt{1+4(b_{i}^{(n)})^{2}}}{2}\)
7. \(\mathbf{\tilde{T}}_{i}^{(n+1)}=\mathbf{T}_{i}^{(n+1)}+\frac{b_{i}^{(n)}-1}{b^{(n+1)}}( \mathbf{T}_{i}^{(n+1)}-\mathbf{T}_{i}^{(n)})\)
8.