## 2 Materials and Methods
### Estimation of Blur Kernel and the Intermediate Latent Image
#### 2.2.2 Solve \(\mathbf{k}\) with a given \(\mathbf{f}\)  
As the differential operations are linear, the blurs in the image gradient domains also satisfy the convolution model in Equation (1). Moreover, the blur kernel \(\mathbf{k}\) can be estimated with more accuracy in the gradient domain than directly solving the second subproblem in Equation (5) with pixel intensities [28]. In the following contexts, the analytical solution for the \(\mathbf{k}\)-subproblem is discussed in the matrix form, and the second part of Equation (5) can be rewritten as:  
\[\hat{\mathbf{K}}=\arg\min_{\mathbf{K}}\left(\mathbf{F}\mathbf{K}-\mathbf{G}\right)^{T}\left(\mathbf{F }\mathbf{K}-\mathbf{G}\right)+\alpha_{2}\mathbf{K}^{T}\mathbf{K} \tag{13}\]  
where \(\mathbf{F}=\left(\mathbf{F}_{1}+\mathbf{F}_{2}+\mathbf{F}_{3}+\mathbf{F}_{4}\right)\) and \(\mathbf{G}=\left(\mathbf{G}_{1}+\mathbf{G}_{2}+\mathbf{G}_{3}+\mathbf{G}_{4}\right)\); \(\mathbf{G}_{n}\) and \(\mathbf{G}_{n}\) with \(n\in N\) represent the block circulant matrices with circulant blocks generated from \(\
abla\mathbf{f}\) and \(\
abla\mathbf{g}\); and \(\mathbf{K}\) is rearranged from \(\mathbf{k}\) in a vector form so as to transform the convolution operation to the matrix multiplication [26,29].  
The optimization of the above equation is a least squares problem, and its closed-form solution can be achieved with the help of fast Fourier transform. However, the division and truncation in the frequency domain will lead to the amplification of both noise and estimation error. However, in this paper, we solve the optimization problem in the pixel domain directly with an effective conjugate gradient method combined with [PERSON]'s method.  
Let \(\mathcal{B}(\mathbf{K})\) stand for the \"\(\arg min\)\" function in the right-hand side of Equation (13), then its gradient \(\
abla\mathcal{B}(\mathbf{K})\) and Hessian matrix \(\mathbf{H}\) are calculated in Equation (14) as follows:\[\left\{\begin{array}{l}\
abla B(\@vec{K})=2\@vec{F}^{T}\@vec{F}K-2\@vec{F}^{T} \@vec{G}+2\alpha_{2}\@vec{K}\\ \\ \@vec{H}=2\@vec{F}^{T}\@vec{F}+2\alpha_{2}\@vec{I}\end{array}\right. \tag{14}\]  
where \(\@vec{I}\) is an identity matrix of the same size as \(\@vec{F}^{T}\@vec{F}\), and the variable \(\@vec{K}\) is eliminated in the Hessian matrix when calculating their second-order partial derivatives. The solution of \(\@vec{K}\) can be achieved by [PERSON]'s descent method. Furthermore, since \(\@vec{B}(K)\) is quadratic, the searching process converges within only one iteration [30]. Assume the starting point is \(\@vec{K}_{0}=\@vec{0}\), then the minimum is achieved for:  
\[\begin{array}{ll}\@vec{K}&=K_{0}-\@vec{H}^{-1}\
abla\@vec{B}(K_{0})\\ &=2\@vec{H}^{-1}\@vec{F}^{T}\@vec{G}\end{array} \tag{15}\]  
Considering that the computation of the inverse of the Hessian matrix \(\@vec{H}\) directly is very complicated and impracticable, and computing the inversion \(\@vec{H}^{-1}\) is not the final aim, the solution of Equation (15) can be converted to the following form:  
\[\@vec{H}\@vec{K}=2\@vec{F}^{T}\@vec{G}. \tag{16}\]  
It can be seen that Equation (16) is a linear system of equations; it can be solved by the conjugate descent method effectively. During the alternating optimization, it should be noted that the blur kernel must be kept positive and normalized so as to satisfy the property of a PSF. As a result, the negative elements of the estimated \(\@vec{K}\) are set to zero, and its elements, \(K_{i}\), are then normalized to hold \(\sum K_{i}=1\).