Toward Designing Convergent Deep Operator Splitting Methods for Task-specific Nonconvex Optimization

1 Introduction

In this work, we consider the optimization problem

where x = {x1, · · · , xN} ∈ RD has N blocks (N ≥ 1), f(x) is continuously differentiable, but a series of gn(x) is not necessarily differentiable. Notice that convexity is not assumed for f, gn or Xn. By considering gn as an extended value function (i.e., take +∞ value), we can incorporate the set constraint xn ∈ Xn into gn(xn) since this is equivalent to minimize the indicator function of Xn. Therefore, we will not include this set constraint in our following analysis. Typically, f captures the loss of data fitting from the specific task modeling and g = ∑Nn=1 gn is the regularization that promotes desired structures on the variable x. Problems appearing in many learning and vision applications, such as sparse coding [Lin et al., 2011], tensor factorization [Xu and Yin, 2017], and image restoration [Liu et al., 2018] can all be (re)formulated in the form of Eq. (1).

∗Corresponding Author.

1.1 RelatedWorks

One of the most prevalent algorithms to solve Eq. (1) is the operator splitting approach. The main idea behind such kind of schemes is to reduce complex problems built from simple pieces into a series smaller subproblems which can be solved sequentially or in parallel. In the past several decades, a variety of splitting methods have been designed and analyzed. For example, [Passty, 1979] provided a prototype of the ForwardBackward Splitting (FBS) and proved its ergodic convergence. [Beck and Teboulle, 2009] presented the convergence rate of proximal gradient (PG) and accelerated proximal gradient (APG, also known as FISTA). Recently, [Davis and Yin, 2016] provided a unified way to analyze the convergence rate of Peaceman-Rachford Splitting (PRS) and Douglas-Rachford Splitting (DRS). It is also known that the widely used Alternating Direction Method of Multipliers (ADMM) can be reformulated within the operator splitting (e.g., DRS) framework in the dual space [Lin et al., 2015]. Though with mathematically proved convergence properties, the generally designed algorithms may still fail on some particular nonconvex optimization models in real scenarios. This is mainly because that due to their fixed updating schemes, it is hard to escape the unwanted saddle points during iterations.

To improve the performance in practical real-world applications, some researches tried to parameterize exiting iteration schemes and learned their parameters in the resulted propagation models. For example, [Liu et al., 2016] learned the parameters of a parameterized partial differential equation for various image and video processing tasks. Similarly, [Chen et al., 2015] introduced a higher-order diffusion system to perform data-dependent gradient descent for image denoising and super-resolution. The studies in [Uwe and Stefan, 2014] parametrized the half-quadratic splitting for practical applications. Very recently, inspired by the success of deep networks in different application fields, some works also tried to replace the standard iterations by existing network architectures. By considering convolutional neural networks (CNNs) as special image prior, [Zhang et al., 2017] proposed an iterative CNN scheme to address image restoration problems.

However, we have to point out that although with relatively better performance in some specific tasks, the nice convergence properties proved from theoretical side are completely missing in these methods. That is, neither the adaptive parameterization nor the replaced CNNs mentioned above can preserve the convergence results proved for the original iteration schemes. Moreover, it is even impossible to investigate and control the iterative behaviors (e.g., descent) of these methods, since their learned iterations actually no longer solve the original optimization model.

1.2 Contributions

In this work, we propose Learnable Bregman Splitting (LBS), a novel deep operator splitting algorithm for nonconvex optimization in real-world scenarios. Specifically, we first introduce a Bregman distance function to penalize the variables at each iteration. Then the basic LBS updating scheme is established based on a relaxed Krasonselskii-Mann iteration [Davis and Yin, 2016]. By introducing a novel triple operator splitting strategy, we can successfully combine the task-model-inspired and data-learning-driven operators within the LBS algorithmic framework. In summary, our contributions mainly include:

• LBS provides a novel learning strategy to extend prevalent mathematically designed operator splitting schemes for task-specific nonconvex optimization. Thanks to the learnable deep architectures, we can learn our iterations on collected training data to avoid unwanted solutions in particular applications.
• Different frommost existing learning-based optimization algorithms (e.g., iteration parameterization and CNN incorporation methods mentioned above), in which there is no theoretical guarantee, we provide rich investigations on the iterative behaviors, prove the global convergence and estimate the convergence rate of our LBS.
• We also demonstrate how to apply our algorithm for different computer vision applications and extensive results verify that LBS outperforms state-of-the-art methods on all the compared problems.

2 Learnable Bregman SplittingMethod

In this section, a learning-based operator splitting method, named Learnable Bregman Splitting (LBS), is developed for the nonconvex optimization model in Eq. (1).

2.1 Bregman Distance Penalization

As a fundamental proximity measure, the Bregman distance 1 plays important roles in various iteration algorithms. However, since it does not satisfy the triangle inequality nor symmetry, this function is not a real metric. Given a convex differential function h, the associated Bregman distance can be written as

1 The use of Bregman distance in optimization within various contexts is well spread. Many interesting properties of this function can be found in the comprehensive work [Bauschke et al., 1997].

Clearly, ∆h is strictly convex with respect to the first argument. Moreover, ∆h ≥ 0 for all (x, y) and is equal to zero if and only if x = y. So ∆h actually provides a natural (asymmetric) proximity measure between points in the domain of h.

In this work, we introduce ∆hn as a penalty term for each xn at t-th iteration. That is, we actually minimize the following energy to update xt+1 :

where we denote ∆h(x, xt ) = ∑Nn=1 ∆hn (xn, xtn ) and λ > 0 is the penalty parameter. It will be demonstrated that {∆hn (xn, xtn )} brings nice convergence properties for the proposed optimization model when h is µ-strong convex [Bauschke et al., 1997].

2.2 UniformCoordinate Updating Scheme

In this work, we consider the following general coordinate update scheme to minimize the energy function in Eq. (3):

where D(xt ) denotes the update direction (regarding to the problem) on xt , γ t > 0 is a step size and [·]n denotes the n-th block of the given variable. It should be pointed out that by formulating D := I-T (here I denotes the identity mapping), Eq. (4) can be further recognized as a relaxed KrasonselskiiMann iteration [Shi et al., 2016] with the operator T (i.e., xt+1 n = (1 - γ t )xtn + γ t [T (xt )]n) and then various existing first-order schemes can be reformulated in the form of Eq. (4).

Specifically, by defining JF = (I - F)-1 (resolvent), RF = 2JF - I (reflection) for F (operator about f), we can obtain a variety of prevalent splitting schemes, such as FBS, PRS, and DRS. As for the operator T in our work, if setting T = I-JG◦(I-F) and γ t = 1, we obtain FBS from Eq. (4), i.e., xt+1 = JG ◦ (I - F)(xt ), where ◦ denotes the operator composition. By considering JG = proxg 2 and F = ∇f , we further have the well-known proximal (or projected) gradient scheme from FBS. Setting T = RG ◦RF and γ t = 1, Eq. (4) reduces to xt+1 = (RG ◦ RF) (xt ), which is just the standard PRS iteration. Similarly, with the same T in PRS and γ = 1/2, we can also deduce DRS.

Additionally, it should be pointed out that the well-known ADMM [Lin et al., 2011] can also be deduced by applying DRS on its Lagrange dual space [Davis and Yin, 2016]. Therefore, although the original ADMM is designed for linearly constraint models, we can still reformulate it as a special case of Eq. (4) in the dual variable space. Thus Eq. (4) actually can also be utilized to address the constrained problems.

2 The proximal operation with respect to g (denoted as proxg ) is defined as proxg (x) := argminy g(y) + 1 2 ‖y - x‖2 .

In consideration of the nice properties and high flexibility of Eq. (4) mentioned above, we would like to utilize the general updating scheme in Eq. (4) with D = I - T , (named Uniform Coordinate Updating Scheme, UCUS for short) as our fundamental block-wise iteration rule.

2.3 Splitting with Learnable Architecture

As discussed above, most existing splitting algorithms (e.g., FBS, PRS and DRS) specify the operator only based on the optimization model. However, due to the nonconvex nature of the model, it is hard for these schemes to escape undesired local minimum. Moreover, the complex data distribution in real applications will also slow down the redesigned iterations.

To partially address these issues, we provide a new splitting strategy, in which a learnable operator Td is introduced to extract information from the data. That is, we consider the following triple splitting scheme:

where Tgn and Tfλ are operators related to Ψth in Eq. (3). Here we just follow a FBS-like strategy to define Tfλ = I -

∇(f + 1λ ∆h) and Tgn = proxρgn . As for Td, we would like to build it as a learnable network architecture and train its parameters from collected training data set 3 . In this way, we can successfully incorporate data information to improve the iterative performance of the proposed algorithm.

Notice that it is challenging to analyze the convergence issues for the existing network-incorporated iterations (e.g., [Zhang et al., 2017]), since all their schemes are built in heuristic manners. In contrast, we will demonstrate in the following section that the convergence of our LBS can be strictly proved.

2.4 The Complete Algorithm

It can be seen that the learnable operator Td are not deduced from strict optimization rule, there may exists iteration errors when calculating T at each stage. Thus we introduce a new condition to control the inexactness of our updating scheme at each iteration. Specifically, we define the optimality errors of a given variable u at t-th iteration based on the first order subdifferential of Ψt+1 h , i.e.,

where dgn ∈ ∂gn(un) (here we denote ∂gn as the limiting Ferchet subdifferential of gn [Xu and Yin, 2017]) and n = 1, 2, · · · ,N. Then we consider the following so-called Relaxed Optimality Condition (ROC) for the given u ∈ RD .

Condition 1. (Relaxed Optimality Condition) Given any u ∈ RD , we define the relaxed optimality condition of Ψt+1 h for un ∈ u (n = 1, 2, · · · ,N) as ‖et+1 un ‖ ≤ c‖un - xtn ‖, where c is a fixed positive constant.

Based on the above condition, we are ready to propose our LBS algorithm for solving Eq. (1) in Alg. 1. Notice that the UCUS iteration, denoted as ucus(·), is independently stated in Alg. 2. It can be seen that if ROC is satisfied, the LBS iterations are fully based on the learnable network operator. While for some iterations, which do not satisfy ROC, we may still perform the model-based operators Tgn ◦ Tfλt to guarantee the final convergence. For convenience, hereafter the subvectors {x1, . . . , xn-1} and {xn+1, . . . , xN} are denoted as xn- and xn+ for short, respectively. We also denote ψt+1 n (xn) = f t+1 n (xn) + gn(xn), in which f t+1 n (xn) = f({xt+1 n- , xn, xtn+}).

3 See Sec. 4 for the details of this operator and its training strategy.

Algorithm 1 Learnable Bregman Splitting (LBS)

Require: x0 , Td, ∆h, c, µ, {γ t |0 < γ t ≤ 1}, {λt |λt > 0}.

1: while not converged do 2: zt+1 = Tfλt ◦ Td(xt ). 3: for n = 1, 2, . . . , N do 4: ut+1 n = Tgn (zt+1 n ). 5: if ut+1 n satisfies ROC then 6: vt+1 n = ut+1 n . 7: else 8: vt+1 n = Tgn ◦ Tfλt ({xt+1 n- , xtn , xtn+}). 9: end if 10: xt+1 n = ucus({xt+1 n- , vt+1 n , xtn+}, γ t ). 11: end for 12: end while

Algorithm 2 xt+1 n = ucus({xt+1 n- , vt+1 n , xtn+}, γ t )

1: wt+1 n = xtn -γ t (xtn - vt+1 n ). 2: if ψt+1 n (wt+1 n ) ≤ ψt+1 n (vt+1 n ) then 3: xt+1 n = wt+1 n . 4: else 5: xt+1 n = vt+1 n . 6: end if

3 Convergence Analysis

In this section, we provide strict analysis on the convergence behaviors of LBS. The following assumptions on the functions f , g, and Ψ are necessary for our analysis. Notice that all these assumptions are fairly loose in optimization area and satisfied in most vision and learning problems.

Assumption 1. 1) f is Lipschitz smooth and gn is proximable 4 . 2) Ψ is coercive.

We first summarize the roadmap of our analysis as follows 5 . Firstly, some important iterative properties of the sequences are proved in Propositions 1, 2, and 3. Based on these results, we prove the global convergence of LBS in Theorem 1. That is, the sequences generated by LBS are Cauchy sequences and converge to the critical points of Eq. (1). Finally, the convergence rate of the sequences is estimated in Corollary 1.

Proposition 1. (Sufficient descent). If c < µ 2λt and ρ < 1L , both the learnable operators Tgn ◦ Tfλt ◦ Td and Tgn ◦ Tfλt in Alg. 1 can get the objective inequality:

4 A function g is proximable if it is easy to obtain the minimizer of g(x) + 1 2β ‖x - y‖ for any given y and β > 0.

5 The detailed proofs of the theoretical results can be found in our technical report: https://arxiv.org/abs/1804.10798.

where M = min t { µ 2λt - c, 1 2ρ - L2 }, L is Lipschitz moduli of ∇f. Together with the direct comparison of function values in Alg. 2, there exists a non-increasing objective sequence, i.e.,

Remark 1. The inequalities in Proposition 1 builds the relationship of ψt+1 n (xtn ) and ψt+1 n (vt+1 n ), thus we can obtain a series of useful inequalities:

where n = N - 1, . . . , 1. It implies the non-increasing property of {Ψ(xt )}t∈N.

Proposition 2. (Square summable). If γ t ∈ (0, 1], {xt }t∈N, {vtn }t∈N are the sequences by Alg. 1, we have

Proposition 3. (Subsequence convergence). Let {xt }t∈N be the sequence generalized by Alg. 1. If x∗ is any accumulation point of {xt }t∈N. Then we have

when j → ∞.

Remark 2. Indeed, Propositions 2 and 3 are the key points to prove that the sequence {xt }t∈N has critical points. Proposition 3 can be derived by combining the Lipschitz smoothness of f, lower semi-continuous of gn, µ-strong convexity of h, and property of the learnable operators Tgn ◦ Tfλt ◦ Td and Tgn ◦ Tfλt in Steps 4 and 8 of Alg. 1.

Theorem 1. (Critical point and Cauchy sequence). Let Assumption 1 hold for Eq. (1), then the sequences {xt }t∈N generated by Alg. 1 has critical points x∗ of Ψ, i.e., if Ψ∗ is the limit of sequence {Ψ(xt )}t∈N, We have

If Ψ is a Kurdyka-Łojasiewicz function 6 , we can further prove that {xt }t∈N is a Cauchy sequence, thus globally converges to a critical point of Ψ.

Based on the above theorem, we can estimate convergence rate as follows.

Corollary 1. Let φ(s) = tθ (s)θ be a desingularizing function with a constant t > 0 and a parameter θ ∈ (0, 1]. Then {xt }t∈N generated by Alg. 1 converges after finite iterations if θ = 1. The linear and sub-linear rates can be obtained if choosing θ ∈ [1/2, 1) and θ ∈ (0, 1/2), respectively.

6 It should be pointed out that many functions arising in learning and vision areas, including ℓ0 norm and rational ℓp norms (i.e., p = p1/p2) are all Kurdyka-Łojasiewicz functions [Liu et al., 2018].

4 Numerical Results

To verify the convergence and performance of LBS for nonconvex optimization, we apply it on two widely researched vision problems, i.e., image completion and deblurring. In our algorithm, we adopt residual network as the learnable network architecture for Td, which can well describe the sparse priors. Specially, there are 19 layers in our network which includes 7 convolution layers, 6 ReLU layers, 5 batch normalization layers and one loss layer. Every convolution layer has 64 kernels of size 3 × 3, and possesses the dilation attribute. In training stage, we randomly select 800 natural images from ImageNet database [Deng et al., 2009]. The chosen pictures are cropped into small patches of size 35 × 35 and Gaussian noise is imposed to these patches. In the following, we always consider Mahalanobis distance ‖ · ‖A (with weight matrix A) [Bauschke et al., 1997] as ∆h in our applications. All experiments are performed on a PC with Intel Core i7 CPU @ 3.4 GHz, 32 RAM and NVIDIA GeForce GTX 1050 Ti GPU.

4.1 ℓp-Sparse Coding for Image Completion

We first consider to solve a ℓp-sparse coding model to address the problem of image completion (also known as image inpainting). The purpose of this task is to restore a visually plausible image in which data are missing due to damage or occlusions. This problem can be formulated as:

where y is the observed image, M denotes a mask, B is the dictionary, α is its corresponding sparse coefficients and ρ > 0 is a parameter. Following [Beck and Teboulle, 2009], we consider B as a inverse wavelet basis (i.e., multiplying by B corresponds to performing inverse wavelet transform) and thus Bα is just the latent image (denoted as x). To enforce the sparsity of α, we utilize nonconvex ℓ0.8 norm in the above coding model. As for the Bergman distance, we set A = µI, where I denotes the unit matrix and µ = 0.01.

It is easy to check that Eq. (6) is just a specific case of Eq. (1) with single variable. In this following, we first verify the theoretical results proposed in this work, and then test the performance of LBS on challenging benchmark datasets.

Iteration Behaviors Analysis: We first choose example images from CBSD68 dataset [Zhang et al., 2017] to demonstrate the iterative behaviors of LBS together with some other widely used splitting schemes (e,g, FBS, FISTA, and ADMM). For fair comparisons, the stopping criterion of all the compared methods are set in the same manner. That is, we denote x = Bα and consider ‖xt+1 - xt ‖/‖xt ‖ ≤ 10-4 as the stopping criterion in all these methods.

Fig. 1 showed the convergence curves from different aspects, including iteration error ("Iter. Error", defined as log(‖xt+1 - xt ‖/‖xt ‖)) and reconstruction error ("Rec. Error", defined as log(‖xt -xgt‖/‖xgt‖)). Our LBS have superiority against traditional FBS, FISTA, and ADMM on both convergence rate and final reconstruction. LBS can achieve the convergence precision only almost a dozen steps while FBS and FISTA need few hundreds steps and ADMM needs four dozens of steps. Since introducing the network as Td, our strategies have less reconstruction error than others obviously. The PSNR and SSIM of the final results also verify that our LBS has better performance. Concretely, our PSNR is approximately higher 3dB than the compared methods.

We also compared the curves of objective function value errors ("Func. Error", based on Φ(xt )) for different settings of LBS, including naive LBS (nLBS, do not check the ROC and monotone conditions) and the complete LBS in Alg. 1. From the left subfigure of Fig. 2, it is easy to observe that the proposed criterion can lead to very fast convergence, while there are severe oscillations on the curves of nLBS. Furthermore, we plotted the bars of ROC (i.e., the error etut and the threshold c‖ut - xt ‖2 ) on the right part of Fig. 2. It can be seen that the ROC condition is always satisfied except at the last two iterations. Thus deep networks are performed at most of our iterations. Only at the last stages, LBS tended to perform model-inspired iterations (i.e., Step 8 in Alg. 1) to obtain accurate solution for the given optimization model.

Comparisons on Benchmarks: To further express the superiority of LBS, we generated random masks of different levels (including 20%, 40%, 60% and 80% missing pixels) on CBSD68 dataset [Zhang et al., 2017] for comparison, which contains 68 images with the size of 481×321. Then we compared LBS with four state-of-the-art methods, namely, FoE [Roth and Black, 2009], VNL [Arias et al., 2011], ISDSB [He and Wang, 2014], and JSM [Zhang et al., 2014]. Tab. 1 reports the averaged quantitative results, including PSNR, SSIM, and time (in second). It can be seen that regardless the proportion of masks, LBS can achieve better performance against the state-of-the-art approaches. This is mainly due to our superior strategy which uses learnable network operator.

We then compared the visual performance of LBS with all these methods. Fig. 3 presented the comparisons on an image from ImageNet database [Deng et al., 2009] with 60% missing pixels. It can be seen that LBS outperformed all the compared methods on both visualization and metrics (PSNR and SSIM). The edge of motorcycle wheels can be restored smoother and clearer by LBS, while other approaches exist some noises and masks to affect the visual effects.

4.2 Nonconvex TV for Image Deblurring

We further evaluate LBS on image deblurring, which is a challenging problem in computer vision area. Here we consider the following widely used total variation (TV) based formulation:

where k, x, y denote the blur kernel, latent image, and blurry observation, respectively. TVp(u) = ‖Dhu‖pp + ‖Dvu‖pp is the nonconvex TV regularization with gradient matrices Dh and Dv (here we also set p = 0.8 for the ℓp norm). χΩu (u) is the indicator function of the set Ωu := {u ∈ Rn : 0 ≤ ui ≤ 1}. Following the half-quadratic splitting technique, Eq. (7) (with auxiliary variables vh and vv) can be reformulated as

Obviously, Eq. (8) is a special case of Eq. (1) with three blocks, thus can be efficiently addressed by LBS. Here we define A = diag (µuIu, µvh Ivh , µvv Ivv ) with µu = 0.01, and µvh = µvv = 0.001 in Bergman distance ∆h.

Fig. 4 demonstrated the convergence behaviors of LBS on {u, vh, vv}. It can be seen from the left subfigure that "Iter. Error" of all blocks quickly decreased to log(10-3 ), notice that "Iter. Error" of u is even less than log(10-4 ). On the right subfigure, the "Rec. Error" of vh and vv also have dramatic decline trend, which are shown along with the right vertical ordinate. Due to the different range of values, we have to plot the curves of u w.r.t. the left vertical ordinate. We can see that it still obtained the least "Rec. Error".

We then reported results on the challenging image deblurring benchmark dataset collected by Sun et al. [Sun et al., 2013] (which includes 640 blurry images with 1% Gaussian noises) for quantitative evaluation. We compared LBS with plenty of competitive approaches, including TV [Wang et al., 2008], HL [Krishnan and Fergus, 2009], CSF [Uwe and Stefan, 2014], IDDBM3D [Danielyan et al., 2012], EPLL [Zoran and Weiss, 2011], RTF [Schmidt et al., 2016], MLP [Schuler et al., 2013], IRCNN [Zhang et al., 2017], FDN [Kruse et al., 2017], and PADNet [Liu et al., 2018]).

It is known that learning-based methods (e.g., CSF, RTF, MLP, IRCNN, FDN, and PADNet) can achieve better performance than other conventional approaches in terms of quantitative metrics (e.g., PSNR and SSIM). However, due to the weak theoretical guarantee, they are worse than LBS (see Tab. 2). Fig. 5 expressed the qualitative results of LBS against other methods (top 4 in Tab. 2) on an example blurry image, which is generated with a large scale blur kernel (75×75 pixels) on an image from ImageNet [Deng et al., 2009]. It can be seen that LBS can restore the text and windows more distinctly than others. Although IRCNN has relatively higher PSNR than others (but lower than LBS), its visual quality and SSIM are not satisfied.

5 Conclusions

This paper proposed Learnable Bregman Splitting (LBS), a novel deep architecture based operator splitting algorithm for task-specific nonconvex optimization. It is demonstrated that both the model-based operators and the data-dependent networks can be used in our iteration. We also provided solid theoretical analysis to guarantee the convergence of LBS. The experimental results verified that LBS can obtain better performance against most other state-of-the-art approaches.

Acknowledgments

This work is partially supported by the National Natural Science Foundation of China (Nos. 61672125, 61733002, 61572096, 61432003 and 61632019), and the Fundamental Research Funds for the Central Universities.