Empirically Accelerating Scaled Gradient Projection Using Deep Neural Network For Inverse Problems In Image Processing

TL;DR

This paper introduces a deep neural network-based iterative algorithm that accelerates traditional optimization methods for large-scale inverse problems in image processing, ensuring convergence and improved empirical performance.

Contribution

It presents a novel DNN approach to learn parameters for the scaled gradient projection method, enhancing convergence speed over conventional algorithms.

Findings

Significantly improved convergence rate in simulations.

Effective acceleration for various large-scale inverse problems.

Demonstrated convergence and empirical performance benefits.

Abstract

Recently, deep neural networks (DNNs) have shown advantages in accelerating optimization algorithms. One approach is to unfold finite number of iterations of conventional optimization algorithms and to learn parameters in the algorithms. However, these are forward methods and are indeed neither iterative nor convergent. Here, we present a novel DNN-based convergent iterative algorithm that accelerates conventional optimization algorithms. We train a DNN to yield parameters in scaled gradient projection method. So far, these parameters have been chosen heuristically, but have shown to be crucial for good empirical performance. In simulation results, the proposed method significantly improves the empirical convergence rate over conventional optimization methods for various large-scale inverse problems in image processing.

Figures

Tables

Table 1. Table 1: Averaged NMSE (dB) of all methods trained with 50%, tested on various cases for inpainting.

Method	Test-30%	Test-50%	Test-70%
FISTA@100	-5.20$\pm$1.01	-10.74$\pm$2.29	-17.69$\pm$2.85
U-Net-50%	-7.05$\pm$0.97	-19.40$\pm$2.91	-9.11$\pm$2.38
Step-learn-50%@20	-13.77$\pm$2.86	-17.73$\pm$3.01	-21.65$\pm$3.05
Diag-learn-50%@20	-13.85$\pm$2.71	-17.50$\pm$2.97	-21.35$\pm$3.05

Table 2. Table 2: Averaged NMSE (dB) of methods trained with 50%, tested on various partial Fourier cases.

Method	Test-30%	Test-50%	Test-70%
FISTA@100	-20.46$\pm$3.49	-24.84$\pm$3.99	-32.31$\pm$4.74
U-Net-50%	-18.22$\pm$3.43	-20.17$\pm$4.08	-19.62$\pm$4.57
Step-learn-50%@20	-20.80$\pm$3.67	-25.49$\pm$4.20	-33.04$\pm$4.79
Diag-learn-50%@20	-20.80$\pm$3.65	-25.52$\pm$4.17	-33.06$\pm$4.78

Table 3. Table 3: Averaged NMSE (dB) of methods trained with 50%, tested on various partial Fourier cases.

Method	Test-50%	Test-50% noisy
FISTA@100	-24.84$\pm$3.99	-22.79$\pm$2.85
U-Net-50%	-20.17$\pm$4.08	-19.85$\pm$3.79
Step-learn-50%@20	-25.49$\pm$4.20	-23.16$\pm$2.78
Diag-learn-50%@20	-25.52$\pm$4.17	-23.18$\pm$2.76

Equations

$$\underset{\mathbf{x}}{\min }F(\mathbf{x})=f(\mathbf{x})+g(\mathbf{x})$$

$${\mathbf{x}}^{(k+1)}={\text{prox}}_{t_k}(g)({\mathbf{x}}^{(k)}-t_k\nabla f({\mathbf{x}}^{(k)})),$$

$$\underset{\mathbf{x}\in S}{\min }f(\mathbf{x})$$

$${\mathbf{y}}^{(k)}={\mathbf{P}}_S({\mathbf{x}}^{(k)}-t_k\nabla f({\mathbf{x}}^{(k)}))$$

$$f({\mathbf{x}}^{(k)}+{\gamma }^{(k)}{\mathbf{z}}^{(k)})⩽f({\mathbf{x}}^{(k)})+\beta {\gamma }^{(k)}\nabla f({\mathbf{x}}^{(k)}{)}^T{\mathbf{z}}^{(k)}$$

$${\tilde{\mathbf{x}}}^{(k),n}={\mathcal{T}}_{\lambda t}\left({\mathbf{x}}^{(k),n}-t_{k,n}\nabla f({\mathbf{x}}^{(k),n})\right)$$

$$\underset{\mathbf{\Theta }}{\min }\frac{1}{2}{\Vert {\mathbf{x}}_{\ast }^n-{\mathcal{T}}_{\lambda t}\left({\mathbf{x}}^{(k),n}-t_{k,n}(\mathbf{\Theta })\nabla \mathbf{f}({\mathbf{x}}^{(\mathbf{k}),\mathbf{n}})\right)\Vert }^2$$

$${\mathbf{x}}^{(k+1),n}={\mathcal{T}}_{\lambda /L}\left({\tilde{\mathbf{x}}}^{(k),n}-t_L\nabla f({\tilde{\mathbf{x}}}^{(k),n}\right)$$

$${\mathcal{L}}_k={\sum }_{i=0}^k{\tilde{\mathcal{L}}}_i$$

$${\hat{\mathbf{x}}}^n=\arg \underset{{\mathbf{x}}^n}{\min }\frac{1}{2}||{\mathbf{A}}_n{\mathbf{x}}^n-{\mathbf{y}}^n|{|}^2+\lambda ||{\mathbf{x}}^n|{|}_1$$

Full text

1. forward models and image priors implicitly for direct mappings from given measurements to solutions,
2. data-driven priors as proximal operators in conventional iterative algorithms, or
3. forward models, priors and/or static stepsizes in unfolded structures of optimization iterations. Here we investigate another way of utilizing convolutional neural network (CNN) for empirically accelerating conventional optimization for solving inverse problems in image processing. We propose a CNN to yield parameters in optimization algorithms that have been chosen heuristically, but have shown to be crucial for good empirical performance. Our CNN-incorporated scaled gradient projection methods, without compromising theoretical properties, significantly improve empirical convergence rate over conventional optimization based methods in large-scale inverse problems such as image inpainting, compressive image recovery with partial Fourier samples, deblurring and sparse view CT. During testing, our proposed methods dynamically select parameters every iterations to speed up convergence robustly for different degradation levels, noise, or regularization parameters as compared to direct mapping approach.

1 Introduction

Optimization based solvers for inverse problems have been widely investigated for image processing applications such as denoising [23], inpainting [28], deblurring [36], image recovery from incomplete Fourier samples [26], and image reconstruction from noisy Radon transformed measurements [35]. A typical pipeline for them is to construct an objective function with accurate forward modeling of image degradation processes and with reasonable image priors such as minimum total variation (TV) and/or sparsity in wavelet (or learned transformed) domain, and then to optimize the objective function to yield a solution, a recovered image, using theoretically well-grounded optimization algorithms such as iterative shrinkage-thresholding algorithm (ISTA) [12], fast ISTA (FISTA) [4], approximate message passing (AMP) based algorithm [11], or alternating directional method of multipliers (ADMM) [7]. Many related works have been proposed to improve theoretical convergence rates of algorithms and/or to design good image priors to regularize ill-posed inverse problems.

Deep neural networks (DNNs) have revolutionized the ways of solving inverse problems in image processing. These methods have significantly outperformed conventional optimization based methods by significantly improving the image quality of solutions and computation speed. There are largely three ways of using DNNs for inverse problems in image processing:

1. direct mapping DNNs from measurements (or analytic reconstructions) to solutions by implicitly learning forward models and image priors [32, 34, 21, 18, 17],
2. DNN based proximal operators for iterative algorithms by explicitly learning image priors [9, 24, 31, 16, 29], or
3. unfolded structure DNNs inspired by conventional optimization iterations by learning forward models, priors and/or stepsizes in optimizations [14, 30, 25, 13, 33, 10, 22]. Recently, it has been theoretically shown that unfolded LISTA has asymptotic linear convergence property for compressive sensing (CS) recovery problems [10, 22]. However, these methods have been applied only to small-scale CS recovery problems due to large number of parameters to train or to determine.

Most previous works determine static stepsizes with Lipschitz constant / training processes or dynamic stepsizes with backtracking. However, we argue that for empirically fast convergence rate, stepsizes in optimization iterations must be determined dynamically for different problems every iterations. Here we investigate an alternative way of utilizing CNNs for empirically accelerating optimization for solving large-scale inverse problems in image processing. We propose CNNs estimating a near optimal stepsize (or a diagonal matrix) per iteration for given problems to accelerate empirical convergence rates of conventional optimizations. These parameters in optimization algorithms have been selected heuristically, but have shown to be crucial for good empirical performance. Our CNN-incorporated scaled gradient projection (SGP) methods, without compromising theoretical properties, significantly improve empirical convergence rate over conventional optimization based methods such as ISTA [12] / FISTA [4] with backtracking in large-scale inverse problems such as image inpainting, CS image recovery with partial Fourier samples, image deblurring and sparse view CT reconstruction. Our ways of using DNN select parameters every iteration to speed up empirical convergence robustly for different degradation levels, noise, or regularization parameters as compared to typical direct mapping CNN ($e.g.$, U-Net [27]).

Here are our contributions: We

1. propose small CNNs to dynamically determine stepsizes in optimizations every iterations,
2. propose CNN-incorporated SGP methods without compromising convergence properties, and
3. demonstrate the performance and robustness of our methods for large-scale inverse problems in image processing.

2 Related Works

DNN-based direct mapping approaches have yielded state-of-the-art performance in image quality and computation speed for inverse problems in image processing such as image inpainting [32], image denoising [34], single image super resolution [21], sparse image recovery [18] and medical image reconstruction [17]. However, they also have limitations such as no mechanism to ensure that a solution is corresponding to a given measurement through forward models by correcting for intermediate errors in a current solution. Moreover, DNN-based direct mapping methods are based on black-box models with limited interpretability on the solutions for inverse problems. In contrast, conventional optimization based approaches often have theoretical guarantees for the exact recovery [8] or have interpretable converged solutions such as K-sparse images in wavelet domain with minimum distances to given measurements through forward models.

DNNs inspired by unfolded optimization iterations were also investigated with finite number of iterations. Learned ISTA (LISTA) [14] was the first work of this type to propose DNNs that implicitly learn forward models, image priors and stepsizes in optimizations from data. LISTA has been extended to unfolded LISTA with learned weights [10] and ALISTA with analytically determined weights [22]. The original LISTA compromised ISTA’s theoretical convergence properties, but there have been efforts to understand convergence properties in unfolded structure of DNNs [25, 13]. Recently, it has been theoretically shown that unfolded LISTA has asymptotic linear convergence property for CS recovery problems [10]. However, these methods have been applied to small-scale CS recovery problems due to large number of parameters to determine. If image size is 256$\times$256 and the compression ratio is about 20%, LISTA [14] and LISTA-CPSS [10] require about 4,465M, 2,727M parameters to train, while our proposed U-Net based CNN requires only 7M parameters for 256$\times$256 images. Recently, ALISTA [22] proposed methods for determining analytical weights using dictionary learning or convolutional sparse coding, but it does not seem to work for general forward models yet.

ADMM-net [30] and ISTA-Net [33] are also based on unfolded optimization iterations of ADMM and ISTA, respectively, but unlike LISTA approaches, they utilized forward models in their networks and trained convolutional neural networks (CNNs) for image priors such as transformations, parametrized non-linear functions and for optimization parameters such as stepsize. Using forward models explicitly allows these methods to deal with large-scale inverse problems. Similarly, there have been works using DNNs only for proximal operators in iterative algorithms [9, 24, 31, 16, 29]. Unlike ADMM-net and ISTA-Net with fixed number of iterations, these methods have flexibility of running any number of iterations for different cases. However, they focused on using DNNs as image priors within conventional optimization framework, rather than investigating acceleration of convergence. Most works in this category have static stepsizes that have been determined heuristically, selected conservatively to ensure convergence (e.g., 1/Lipschitz constant), or learned from data.

Lastly, there have been a few recent attempts to learn optimization algorithms as DNN based functions of gradients [1] and as policies of selecting algorithms using reinforcement learning [20]. These are similar to our proposed methods in dynamically determining algorithms every iterations. However, our methods are within the framework of SGP methods with theoretical convergence properties.

3 Background

3.1 Proximal gradient method

Consider an optimization problem of the form

[TABLE]

where $f$ is convex, differentiable and $g$ is convex, subdifferentiable for $\mathbf{x}$. Then, the following update equation at the $k$th iteration is called the proximal gradient method (PGM):

[TABLE]

where $\text{prox}_{t}(g)(\mathbf{x})=\arg\min_{\mathbf{z}}g(\mathbf{z})+\nicefrac{{1}}{{2t}}\|\mathbf{z}-\mathbf{x}\|^{2}$ and $\|\cdot\|$ denotes $l_{2}$-norm. PGM guarantees the convergence of $F(\mathbf{x}^{(k)})$ to $F^{*}$ at the solution $\mathbf{x}^{*}$ with the rate of $O(1/k)$.

One way to determine $t_{k}$ is based on a majorization-minimization technique for the cost function $F$ and its quadratic surrogate function for each iteration. ${\nabla}f$ is usually assumed to be Lipschitz continuous on a given domain and the reciprocal of a Lipschitz constant of ${\nabla}f$ is used for $t_{k}$. Another popular method to choose $t_{k}$ is a backtracking method. Note that both methods do not seek the largest possible stepsize since it is often more efficient to calculate the next iteration with conservative sub-optimal stepsize than to perform time-consuming stepsize optimization. Thus, if there is a way to quickly calculate near optimal stepsizes, it could help to accelerate empirical convergence rates.

3.2 Scaled gradient projection method

Recently, scaled gradient projection (SGP) methods with the convergence rate $O(1/k)$ have been proposed and have empirically demonstrated their general convergence speed improvements over FISTA with $O(1/k^{2})$ [6, 5]. The problem (1) can be seen as a constrained optimization problem

[TABLE]

where $S=\{\mathbf{x}\in\mathbb{R}^{n}:g(\mathbf{x})\leqslant\epsilon\}$ that is a convex set due to the convexity of $g$ and $\epsilon$ is determined by $g$. Then, an iterative algorithm can be formulated as the PGM given by

[TABLE]

where $\mathbf{P}_{S}(\mathbf{x})=\arg\min_{\mathbf{z}\in S}\|\mathbf{z}-\mathbf{x}\|^{2}.$ Whenever $\mathbf{y}^{(k)}\neq\mathbf{x}^{(k)}$, $\mathbf{y}^{(k)}-\mathbf{x}^{(k)}$ is a descent direction at $\mathbf{x}^{(k)}$ for the problem (3) and thus its inner product with ${\nabla}f(\mathbf{x}^{(k)})$ becomes negative. $\mathbf{y}^{(k)}=\mathbf{x}^{(k)}$ implies that $\mathbf{x}^{(k)}$ is a stationary point. Since $\mathbf{y}^{(k)}-\mathbf{x}^{(k)}$ is a descent direction at $\mathbf{x}^{(k)}$, Armijo line search can generate a convergent sequence of $\{\mathbf{x}^{(k)}\}_{k=0}^{\infty}$ that satisfy Armijo condition [3]:

[TABLE]

where $\mathbf{z}^{(k)}=\mathbf{y}^{(k)}-\mathbf{x}^{(k)}$ and $\beta,\gamma^{(k)}\in(0,1]$ for $\forall k$.

In the problem (3), SGP methods introduced an additional symmetric positive definite matrix $D_{k}$ in front of ${\nabla}f(\mathbf{x}^{(k)})$. Symmetry and positive definiteness are necessary conditions for a Hessian matrix of $f$ for Newton’s method and they are also important conditions for quasi-Newton methods. Newton-type methods usually converge with fewer iterations than first-order optimization methods, but they are computationally demanding especially for large-scale input data. SGP methods exploit symmetry and positive definiteness with the aim of less computational burdens while they can refine a direction vector $\mathbf{d}^{(k)}$ of the diagonal elements in $D_{k}$ to accelerate convergence rate. SGP methods are based on Armijo line search since $\mathbf{d}^{(k)}$ remains as a descent direction with the conditions on $D_{k}$ at $\mathbf{x}^{(k)}$. They are also applicable for proximal operators.

For the convergence of SGP methods, it requires additional condition for $D_{k}$. Define $\mathcal{D}_{\delta}$ for $\delta\geqslant 1$ as the set of all symmetric positive definite matrices whose eigenvalues are in the interval $[1/\delta,\delta]$. Then, for $\delta_{k}$ such that $D_{k}\in\mathcal{D}_{\delta_{k}}$, the condition $\Sigma_{k=0}^{\infty}(\delta_{k}^{2}-1)<\infty$, $\delta_{k}\geqslant 1$ should be satisfied. $\lim_{k\rightarrow\infty}D_{k}$ becomes an identity matrix and an iteration becomes similar to PGM. Appropriate $D_{k}$ [5] accelerated empirical convergence over fast PGMs such as FISTA.

For a given proximal operator $\text{prox}_{D}(g)(\mathbf{x})=\arg\min_{\mathbf{z}}g(\mathbf{z})+\nicefrac{{1}}{{2}}(\mathbf{z}-\mathbf{x})^{T}D^{-1}(\mathbf{z}-\mathbf{x}),$ the SGP method is summarized in Algorithm 1. Finding $\{D_{k}\}_{k=0}^{\infty}$ that can accelerate convergence still remains as an open problem. We propose to replace these heuristic decisions with DNNs.

4 Learning-based stepsize selection

We conjecture that DNN can be trained to generate a near-optimal stepsize per iteration if a current estimate and the gradient of a cost function at that estimate are given. Since no ground truth is available for the optimal sequence of stepsizes over all iterations, we propose to train a stepsize DNN to yield a greedy near-optimal stepsize per iteration by minimizing the distance between the estimated vector at the next iteration and the converged solution at each iteration.

4.1 Learning a stepsize for an iteration

To learn stepsizes by a DNN, a set of $N$ solution vectors $\{\mathbf{x}_{*}^{n}\}_{n=1}^{N}$ of optimization problems was generated and used as ground truth data. Solution vectors can be obtained by optimizing the original problems using any convex optimization algorithm ($e.g.$, FISTA with 1200 iterations). Suppose that the estimates at the $k$th iteration form a set of training data or $\{\mathbf{x}^{(k),n}\}_{n=1}^{N}$ that will be fed into the DNN. We denote the output of the DNN as a set of positive real numbers $\{t_{k,n}\}_{n=1}^{N}$ for stepsizes. Then, a set of vectors $\{\tilde{\mathbf{x}}^{(k),n}\}_{n=1}^{N}$ can be obtained at the next iteration:

[TABLE]

where $\mathcal{T}_{\lambda t}(\mathbf{x})_{i}=\max\{|x_{i}|-{\lambda t},0\}\text{sgn}(x_{i})$ for the $i$th element in the vector (soft threshold).

The desired stepsizes $\{t_{k,n}\}_{n=1}^{N}$ for the $n$th images at the $k$th iteration can be obtained by training the DNN $\Psi(\cdot;\bf\Theta)$ to minimize the following loss function with respect to $\bf\Theta$:

[TABLE]

where $t_{k,n}(\bf\Theta)=\Psi\left(\mathbf{x}^{(k),n},{\nabla}f(\mathbf{x}^{(k),n});\bf\Theta\right)$ and then by evaluating $t_{k,n}=\Psi\left(\mathbf{x}^{(k),n},{\nabla}f(\mathbf{x}^{(k),n});\bf\Theta\right).$ After $\{\tilde{\mathbf{x}}^{(k),n}\}_{n=1}^{N}$ are evaluated by the learned stepsizes $\{t_{k,n}\}_{n=1}^{N}$ using (6), we propose to generate another set of vectors for the next iteration $\{\mathbf{x}^{(k+1),n}\}_{n=1}^{N}$ by using a conventional stepsize based on Lipschitz constant $L$:

[TABLE]

where $t_{L}=1/L$. This additional step was necessary since $\{\tilde{\mathbf{x}}^{(k),n}\}_{n=1}^{N}$ were not often improved over $\{\mathbf{x}^{(k),n}\}_{n=1}^{N}$ when the DNN training was not done yet.

To sum, one iteration of our proposed method consists of two steps: 1) the first operation moves a current estimate towards its solution, 2) the second operation is applied for initial training of DNNs. In our simulations, our proposed training method worked well to reduce the loss quickly.

The same training method with a diagonal matrix can be applied by replacing the stepsize $t_{k,n}$ in (6) with a diagonal matrix $D_{k,n}$. The output dimension of the DNN $\Psi(\cdot;\bf\Theta)$ must be changed from $1$ to $G$ with its backpropagations.

4.2 Learning stepsizes for further iterations

Now we propose to further train DNNs to generate stepsizes for multiple iterations. Inspired the training strategy in [15], we define the following cumulative loss function:

[TABLE]

where $\tilde{\mathcal{L}}_{k}=\nicefrac{{1}}{{2}}{\sum_{n=1}^{N}}\|\mathbf{x}_{*}^{n}-\tilde{\mathbf{x}}^{(k),n}\|^{2},$ $\tilde{\mathbf{x}}^{(k),n}$ is defined in (6), and new input datasets as well as ground truth labels are defined as $\mathcal{I}_{k}=\cup_{i=1}^{k}\{\mathbf{x}^{(i),n},{\nabla}f(\mathbf{x}^{(i),n})\}_{n=1}^{N}$ and $\mathcal{O}_{k}=\cup_{i=1}^{k}\{\mathbf{x}_{*}^{n}\}_{n=1}^{N}$ where $\mathcal{O}_{k}$ contains duplicated sets.

Suppose that the DNN is to learn stepsizes of the first $K$ iterations. Initially, the DNN is trained with the input data set $\mathcal{I}_{0}$ and the ground truth label $\mathcal{O}_{0}$ at the [math]th iteration using the procedure in Section 4.1. Then, in the next iteration, the DNN is re-trained with the input data set $\mathcal{I}_{1}$ and the ground truth label $\mathcal{O}_{1}$ at the first iteration. This training process is repeated $K+1$ times so that the DNN can be trained cumulatively as summarized in Algorithm 2.

We expect that our trained DNN should yield near-optimal stepsizes for the first $K$ iterations, but may not be able to yield good stepsizes in later iterations that are larger than $K$. Thus, $K$ should be selected based on the trade-off between image quality and computation time.

5 DNN-incorporated convergent algorithms

5.1 SGP method as framework

The SGP method is described in Algorithm 1. If $D_{k}$ is an identity matrix and $\gamma^{(k)}$ is equal to 1 for $\forall k$, the SGP method is reduced to the PGM. Thus, the SGP method is a generalized version of the PGM by additionally multiplying a symmetric positive definite matrix with the gradient of a loss function that guarantees $\mathbf{z}^{(k)}$ to be a descent direction and by enforcing the Armijo condition for convergence. However, there is no known method to determine $\{D_{k}\}_{k=0}^{\infty}$ that can accelerate and guarantee convergence. We propose DNN to determine $\{D_{k}\}_{k=0}^{\infty}$ that can be trained using the learning procedure in Section 4. Since it is also possible for the DNN to yield $\{D_{k}\}_{k=0}^{\infty}$ that may not satisfy necessary conditions, we proposed to relax the SGP method to selectively use DNN based stepsize (or diagonal matrix) estimation or conservative Lipschitz constant based stepsize to guarantee convergence as summarized in Algorithm 3. We call this the direction relaxation scheme (DRS).

5.2 Proposed relaxation algorithms with DNN

For the DRS in Algorithm 3, note that $\text{prox}_{D}(g)=\mathcal{T}_{\mathbf{d}}$ when $g$ is $l_{1}$-norm and $D$ is a diagonal matrix whose diagonal elements form a vector $\mathbf{d}$. $\mathbf{z}^{(k)}_{1}$ represents a wild search direction generated by the trained DNN and $\mathbf{z}^{(k)}_{2}$ is a conservative search direction from conventional Lipschitz constant based stepsize. Then, depending on the relationship between $\mathbf{z}^{(k)}_{1}$ and $\mathbf{z}^{(k)}_{2}$, the final search direction will be either a linear combination of both of them or $\mathbf{z}^{(k)}_{2}$ alone. For the DNN to generate a single stepsize, $D_{k}$ will be an identity matrix multiplied by that stepsize.

We propose to incorporate the DNN based DRS method into the SGP algorithm as detailed in Algorithm 4. $\mathbf{z}^{(k)}$ is a search direction to yield the estimated vector for the next iteration. As in Algorithm 3, $\mathbf{z}^{(k)}$ is either the weighted average of $\mathbf{z}^{(k)}_{1}$ and $\mathbf{z}^{(k)}_{2}$ with $\gamma_{1}^{(k)}$ or $\mathbf{z}^{(k)}_{2}$ itself. The value of $\gamma_{1}^{(k)}$ was initially set to be $1$ and it remains the same or decreases by a factor $\eta_{1}$ over iterations depending on the Armijo condition at each iteration $k$. The ratio of the weight for $\mathbf{z}^{(k)}_{1}$ to the weight for $\mathbf{z}^{k}_{2}$ in $\mathbf{z}^{k}$ is evaluated at each iteration. Initially, $\mathbf{z}^{(k)}_{1}$ using the trained DNN is dominant in $\mathbf{z}^{(k)}$, but $\mathbf{z}^{(k)}=\mathbf{z}^{(k)}_{2}$ for later iterations and the DNN will not be used eventually. Thus, our proposed algorithm is initially the SGP with DNN search directions and becomes the PGM for later iterations.

Note that the proposed DRS method with relaxation only determines a search direction for the next estimate in a descent direction for inverse problems. $\gamma^{(k)}_{2}$ is the final stepsize parameter, starting from $1$ and decreases its value by a factor $\eta_{2}$ until it satisfies the Armijo condition. Therefore, our proposed method in Algorithms 3 and 4 using the trained DNN is converging theoretically.

6 Simulation results

6.1 Inverse problem settings

We performed various image processing simulations such as image inpainting, CS image recovery with partial Fourier samples, image deblurring and large-scale sparse-view CT to evaluate our proposed methods. An optimization problem for inverse problems in image processing has the following form:

[TABLE]

where $n$ is an index for image, $\mathbf{A}_{n}$ is a matrix that describes an image degradation forward process, $\mathbf{y}^{n}$ is a measurement vector and $\lambda$ is a regularization parameter to balance between data fidelity and image prior. An image is modeled to be in wavelet domain (three level symlet-4) so that $\mathbf{x}^{n}$ is a wavelet coefficient vector for an image. So, the linear operator $\mathbf{A}_{n}$ is a measurement matrix (different for image by image) with an inverse sparsifying transform $\mathbf{W}^{T}$. For normalized measurement matrices, their Lipschitz constants are less than $1$ and normalized gradients helped to yield better results. BSDS500 dataset [2] with $256\times 256$ images was used for all simulations where 450 / 50 images were used for training / testing, respectively. $\lambda$ was set to be 0.1.

We implemented our stepsize (or diagonal matrix) DNN based on U-Net architecture [27] and modified FBPConvNet [17] using MatConvNet on MATLAB. Note that the input and output for the DNN are in a sparsifying transform domain and they have improved the overall performance of inverse problems. $3\times 3$ convolution filters are used for all convolutional layers and batch normalization and rectified linear unit (ReLU) were used after each convolution layer. $2\times 2$ max pooling was applied in the first half of the DNN and deconvolution layers / skip connections were used for the second half of the DNN. We reduced the number of layers in the original FBPConvNet to lower computation time (7M params). For stepsize learning, one fully connected layer was added at the end of the DNN to generate a single number. All simulations were run on an NVIDIA Titan X.

We compared our proposed SGP methods with stepsize DNN (called Step-learned) and diagonal matrix DNN (called Diag-learned) to conventional algorithms such as ISTA and FISTA with backtracking. We also compared our proposed methods with the U-Net that was trained to yield ground truth converged images for (10) from input measurements. We chose U-Net to compare since it has been shown to yield good results in various image processing problems with large-scale inputs and with different forward models for images including inpainting and compressive image recovery [19]. Unfortunately, we were not able to compare ours with ReconNet [18] or LISTA-CPSS [10] that were limited to small-scale CS image recovery problems with the image patch sizes of $33\times 33$ and $16\times 16$, respectively, and empirically one fixed forward model $\mathbf{A}_{n}=\mathbf{A},\forall n$.

NMSE was used for evaluation criteria [10] using $10\log_{10}(\|\mathbf{W}^{T}(\hat{\mathbf{x}}-\mathbf{x}^{\mathrm{gt}})\|^{2}/\|\mathbf{W}^{T}\mathbf{x}^{\mathrm{gt}}\|^{2}).$ Note that all DNNs were trained with converged solutions $\mathbf{x}^{*}$ while the evaluations were done with the oracle solutions $\mathbf{x}^{\mathrm{gt}}$ (the original BSDS500 dataset) for better evaluating the robustness of all methods under significantly different forward models and additional measurement noise. All DNNs for our Step-learned/Diag-learned methods as well as U-Net to yield solutions directly were trained for the case of recovering from noiseless 50% samples with different forward models for different training images. Then, all methods were tested on noiseless 30%, 50%, and 70% samples for test dataset to evaluate performance as well as robustness. For image recovery with partial Fourier samples, noisy 50% samples (Gaussian with standard deviation 5) were also used for further evaluation. We also performed image deblurring for Gaussian blur kernel with $\sigma=2$ and sparse-view CT reconstruction with 144 views using $512\times 512$ CT images to show the feasibility on other applications.

6.2 Image inpainting

Our proposed methods were applied to inpainting problems with different sampling rates. $\mathbf{A}_{n}=\mathbf{M}_{n}\mathbf{W}^{T}$ was used where $\mathbf{M}_{n}$ is a sampling matrix that is different for image. All DNNs were trained with 50% samples (called U-Net-50%, Step-learned-50%, Diag-learned-50%). All results tested on different sampling rates are reported in Figures 1, 2 and Table 1.

When all methods were tested on the same sampling rate (50%), non-iterative U-Net instantly yielded the best image quality among all methods including our proposed methods and even the ground truth converged images (FISTA at 1200 iteration yielded -18.15dB). Both of our proposed methods at 20 iteration yielded image qualities comparable to the converged images, while FISTA at 100 iteration performed much worse than our proposed methods at 20 iteration.

However, U-Net did not show robust performance and yielded substantial artifacts for the tests with different sampling rates such as 30% and 70%. However, our proposed methods without re-training yielded robust accelerations for different test cases as illustrated in Figure 2.

6.3 Image recovery with partial Fourier samples

Similar simulations were performed for image recovery with partial Fourier samples. Note that the input image of the DNN has four channels such that the first two channels are the real and imaginary of the estimated image. Initial images in Figure 3 were obtained using inverse Fourier transform with zero padding. All results are reported in Figures 3, 4 and Tables 2, 3.

For all test cases, non-iterative U-Net yielded the worse results among all methods including FISTA at 100 iteration. Forward models for partial Fourier sampling is much more complicated than forward models for inpainting, and thus more complicated DNN with much more dataset seems desirable for better performance. Our proposed methods without re-training yielded robust and excellent performance at early iteration for all test cases including the same sampling rate (50%), different sampling rates (30, 70%), and additional measurement noise case (noisy 50%). Thus, our proposed DNNs does seem robust to different models and noise in compressive sensing recovery with partial Fourier samples that was inherited from conventional optimization based algorithms.

6.4 Robustness to regularization parameters

We investigate the robustness of our proposed methods by running them that were trained with the original regularization parameter on the test set with different regularization parameters that are half of the original value and twice (2x) of the original value.

Figures 2 (d), 4 (d) illustrate that our proposed methods were robust to small changes in regularization parameters such as half or twice. However, large changes such as 10 times smaller or larger than the original parameter seem to break fast empirical convergence properties of our proposed methods. These phenomena were expected since changing regularization parameters leads to changing ground truth images, thus our DNNs whose inputs are dependent on current estimate and its corresponding gradient should behave in a different way ($e.g.$, if the current estimate is the same as the converged solution with original regularization parameter, then, zero stepsize should be obtained for the problem with the same regularization parameter, but non-zero stepsize should be obtained for different problem with different regularization parameter). Thus, large changes in regularization parameter may require re-training the DNN.

6.5 Image deblurring

Proposed methods were applied to image deblurring problems. Images were blurred using Gaussian kernel with $\sigma=2$. Then, image deblurring was performed with the regularization parameter 0.00001. Note that the initial data fidelity term for deblurring problem is usually much larger than other inverse problems such as inpainting problems. Unlike other inverse problems in image processing, learned diagonal matrix based relaxed SGP yielded the best image quality among all compared methods as shown in Figure 5 qualitatively and quantitatively. It seems that large discrepancy in the data fidelity term was quickly compensated when using the learned diagonal matrix in SGP.

6.6 Sparse-view medical image reconstruction

Lastly, our proposed method was investigated for sparse-view CT image reconstruction. Initial image in Figure 6 (a) was obtained by filtered back-projection from 144 views of projections and had streaking artifacts. With the regularization parameter 0.0005, we ran FISTA-b and proposed diag-learned SGP. At the 10th iteration, our proposed method yielded visually better image than FISTA as illustrated in Figure 6 (b) and (c). Figure 6 (d) shows that our proposed method achieved faster convergence rate than FISTA-b.

6.7 Limited robustness to other measurements

Robustness of our trained DNN for determining stepsize or diagonal matrix when different forward models are used. Our proposed methods have shown robustness for the problems such as image inpainting. It shows that the trained DNN still achieved much faster convergence than FISTA. Similar tendency was observed for image inpainting with 70% sampling. The trained DNN also yielded robust performance for lower or higher (70%) sampling in partial Fourier image recovery. For image deblurring problems with different blur levels, the trained DNN yielded sub-optimal performance compared to FISTA. However, note that step-learned SGP yielded relatively robust performance to diag-learned SGP and it yielded better performance than FISTA for early iterations. For sparse-view image reconstruction, the trained DNN with 144 views did not yield good performance for the test with 45 views. Thus, the robustness of the trained DNN to other forward models seems application-dependent. However, many DNN based algorithms are not robust to other measurement models [17].

7 Conclusion

We proposed a new way of using CNNs for empirically accelerating convergence for inverse problems in image processing with dynamic parameter selections over iterations with different forward models and without breaking theoretical properties such as convergence and robustness. Our trained DNN enabled SGP to empirically outperform FISTA that is theoretically faster than SGP and yield robust performance compared to direct mapping DNNs.

Acknowledgments

This work was supported partly by Basic Science Research Program through the National Research Foundation of Korea(NRF) funded by the Ministry of Education(NRF-2017R1D1A1B05035810), the Technology Innovation Program or Industrial Strategic Technology Development Program (10077533, Development of robotic manipulation algorithm for grasping/assembling with the machine learning using visual and tactile sensing information) funded by the Ministry of Trade, Industry & Energy (MOTIE, Korea), and a grant of the Korea Health Technology R&D Project through the Korea Health Industry Development Institute (KHIDI), funded by the Ministry of Health & Welfare, Republic of Korea (grant number: HI18C0316).