IPGM: Inertial Proximal Gradient Method for Convolutional Dictionary Learning

Abstract

Inspired by the recent success of the proximal gradient method (PGM) and recent efforts to develop an inertial algorithm, we propose an inertial PGM (IPGM) for convolutional dictionary learning (CDL) by jointly optimizing both an ${\ell }_2$-norm data fidelity term and a sparsity term that enforces an ${\ell }_1$ penalty. Contrary to other CDL methods, in the proposed approach, the dictionary and needles are updated with an inertial force by the PGM. We obtain a novel derivative formula for the needles and dictionary with respect to the data fidelity term. At the same time, a gradient descent step is designed to add an inertial term. The proximal operation uses the thresholding operation for needles and projects the dictionary to a unit-norm sphere. We prove the convergence property of the proposed IPGM algorithm in a backtracking case. Simulation results show that the proposed IPGM achieves better performance than the PGM and slice-based methods that possess the same structure and are optimized using the alternating-direction method of multipliers (ADMM).

1. Introduction

Sparse representation is a popular a priori mathematical modeling approach in various signal and image processing applications. Inspired by deep learning-based convolutional operations, convolutional sparse representation is a hot topic at present. In the convolutional sparse representation model, the convolutional dictionary learning (CDL) process plays an important role, but the associated algorithm and convergence proof are difficult problems [1]. CDL involves training dictionaries and estimating codes from multiple signals. In the past decade, the convolutional sparse representation model has achieved impressive results in various applications. The resulting model can achieve excellent performance in image denoising and repair [2], image decomposition [3], image reconstruction [4], medical imaging [5,6], trajectory reconstruction [7], image segmentation [8], superresolution [9], audio processing [10], and other applications. CDL research is a typical interdisciplinary research subject combining mathematics and artificial intelligence. The research results in this field have important theoretical and practical value for a variety of signal and image processing applications. In such a case, a typical assumption is that a signal $y\in {ℝ}^N$, $y=D\Gamma$ is a linear combination of columns, also known as atoms. The coefficient $\Gamma \in {ℝ}^M$ is a sparse vector. A matrix $D\in {ℝ}^{N\times M}$ is called a dictionary. Given $y$ and $D$, this task of finding its sparsest representation is equivalent to solving the following problem:

$$\underset{\Gamma }{\min }{\Vert \Gamma \Vert }_0\mathrm{s}.\mathrm{t}.{\Vert y-D\Gamma \Vert }_2\le \mathsf{\epsilon }$$

(1)

where $\mathsf{\epsilon }$ represents the degree of model mismatch or the additive noise intensity. The solutions to such problems can be approximated by greedy algorithms (such as orthogonal matching pursuit (OMP) [11]) or convex formulas (such as basis pursuit (BP) [12]). The task of the developed learning model is to identify the dictionary $D$ that can best represent a set of training signals, called dictionary learning. Several processing methods have been proposed, including K-singular value decomposition (K-SVD) [13], the method of optimal directions (MOD) [14], online dictionary learning [15], and trainlet learning [16]. Unfortunately, many real-world signals, such as images and audio, are high-dimensional signals, making sparse coding computationally challenging. To avoid the curse of dimensionality, high-dimensional signals are decomposed into low-dimensional overlapping signal patches, and sparse coding is performed independently on each signal block [17]. Because of its simplicity and high performance, this method has been widely utilized in a successful manner. However, this method has some limitations that its patch-based process ignores the correlations between patches.

Another option that contrasts with this local paradigm is a global model [18] called the convolutional sparse representation model. The input signal is represented as the superposition of the convolutional results of some local atoms and sparse feature mappings. The problem is solved by applying a specific structure to the global dictionary involved. In particular, the dictionary in this model is constrained to a banded circulant matrix formed by atomic cascades, which is called convolutional sparse representation.

The convolutional sparse representation model is as follows:

$$y=D\Gamma ={\displaystyle \sum _m{d}_m\ast {\gamma }_m}$$

(2)

where a signal $y\in {ℝ}^N$, the banded global convolutional dictionary $D\in {ℝ}^{N\times N\mathrm{M}}$ is composed of a local convolutional dictionary filter, ${\left\{d_m\right\}}_{m=1}^M\in {ℝ}^n$, a cyclic shift composed of the atom $D_L\in {ℝ}^{n\times M}$, and the coding coefficient $\Gamma \in {ℝ}^{NM}$ is composed of ${\left\{{\gamma }_i\right\}}_{i=1}^M$, where ${\gamma }_i\in {ℝ}^N$.

This work proposes a theoretical analysis of a new global convolutional sparse representation model, which guarantees the local sparsity metric. The convolutional dictionary directly determines translation invariance to compensate for the problem that the correlations between adjacent signal patches are ignored. By using this convolutional dictionary structure, convolutional sparse coding (CSC) and CDL have arisen. For the convenience of calculation, the ${\ell }_1$ norm of the coefficient is used instead of the ${\ell }_0$ norm, and a spherical constraint is used for the dictionary instead of the sphere constraint. Several convex and relaxed CDL algorithms have been proposed [19,20,21,22]. They mainly utilize the alternating-direction method of multipliers (ADMM) solver based on the Fourier domain. However, the ADMM loses its connection to the patch-based processing paradigm that is widely used in many signal and image processing applications [23].

Papyan proposed a new convex relaxation algorithm [23] via slice-based dictionary learning (SBDL), where the sum of the slices forms the signal patches. The developed algorithm is an ADMM solver based on the signal domain. It adopts a local point of view and trains a convolutional dictionary only according to the local calculation of the signal domain. It runs on local slices and faithfully solves the problem of global convolutional sparse coding. This local–global method and its resulting decomposition follow a recent work [18] that compared it with the Fourier-based method. The SBDL algorithm is easy to implement, is intuitive, achieves the most advanced performance, and converges faster than other approaches. Moreover, it provides a better model that can naturally allow a different number of nonzero values to be contained in each spatial location according to the local signal complexity. The ADMM parameters strongly depend on the given problem. Although its corresponding pursuit algorithm can be proven to be convergent, when the convolutional dictionary is iteratively updated in the ADMM algorithm, it is difficult to prove the convergence of the tracking algorithm [24].

The convex relaxation and ADMM-based optimization algorithm produce a nonsparse coding solution, and it is challenging to prove the algorithmic convergence [25]. In addition, in nonconvex optimization, the greedy algorithm for CDL problems has a high computational cost, poor performance, and convergence proof difficulty [26]. Chun and Fessler [27] recently proposed an algorithm that could achieve full convergence. However, the method involves approximate convex constraints, and its overall performance is slightly better than that of the ADMM. The CDL problem is essentially a nonconvex and nonsmooth optimization problem, as shown in the following formula, making it difficult to propose an optimization algorithm and prove its convergence.

$$\underset{\left\{d_m\right\},\left\{{\gamma }_{l,m}\right\}}{\arg \min }\frac{1}{2}{{\displaystyle \sum _l\Vert {\displaystyle \sum _{m=1}^M{d}_m\ast {\gamma }_{l,m}-y_l}\Vert }}_2^2+{\lambda }_1{\displaystyle \sum _l{\displaystyle \sum _{m=1}^M{\Omega }_1\left({\gamma }_{l,m}\right)}}+{\lambda }_2{\displaystyle \sum _{m=1}^M{\Omega }_2\left(d_m\right)}$$

(3)

where ${\lambda }_i,i=1,2$ is the equilibrium coefficient and ${\Omega }_i,i=1,2$ is the nonconvex constraint of the coefficient and the convolutional dictionary. For example, in a case with typical nonconvex constraints, ${\Omega }_1\left(x\right)={\Vert x\Vert }_0$ is a zero norm, and ${\Omega }_2(d_m)={\Vert d_m\Vert }_2^2=1$ is a spherical constraint.

Peng [28] realized the joint and direct optimization of the CDL problem under a nonconvex and nonsmooth optimization scheme and proposed a forward–backward splitting algorithm based on the Fourier domain. The developed approach is better than the ADMM algorithm. To be more precise, the forward step involves estimating the smooth part of the objective function through a partial gradient. In contrast, the backward step counts the degree of nonsmoothness of the objective function through its proximal operator. Peng proved the convergence of the solution sequence of the proposed algorithm by using the semialgebraic property of reference [29] and the Kurdyka–Lojasiewicz (KL) property. Peng [28] used the gradient descent algorithm in the dictionary and coding process. Although this dramatically reduces the computational complexity of the overall algorithm, in theory, the gradient descent algorithm can only guarantee that it will reach the lowest local point, not the lowest global point. Many minimum points are contained in many complex functions. In many cases, we can only obtain the optimal local solution by using the gradient descent method and not the optimal global solution. In addition, when the sample size of the given dataset is large, the convergence speed of the gradient descent algorithm is slow.

Polyak [30] proposed the heavy-ball method that an inertia term is added to the standard gradient descent method. This method has a faster convergence rate than the standard gradient method based on an unchanging number of required calculations. Peter Ochs [31] applied this method to a convex optimization scheme. He proposed an iPiano algorithm combining an inertia term and a forward–backward splitting frame to address minimization problems consisting of differentiable (possibly nonconvex) and convex (possibly nondifferentiable) functions. These problems were strictly analyzed, and the global convergence of the function values and parameters was determined. Simultaneously, the efficiency of convergence was greatly improved. We will apply this inertia term algorithm to the CDL problem.

In this paper, an inertial forward–backward splitting algorithm is proposed for the CDL problem. The convergence of the algorithm is given and proven. The optimal convergence rate is derived. Finally, a large number of experiments show that the proposed algorithm has high efficiency and effectiveness.

The rest of this article is organized as follows. In Section 2, we introduce some related knowledge. In Section 3, the inertia forward and backward splitting algorithm is proposed to restate the CDL problem. In Section 4, the complexity of the proposed algorithm is analyzed. In Section 5, the convergence of the proposed algorithm is analyzed and proven. The convergence rate of the proposed algorithm is derived in Section 6. In Section 7, the performance of the proposed algorithm is evaluated through experiments and compared with other existing methods. Section 8 summarizes the full text.

2. Related Knowledge

2.1. Convex CDL in the Time Domain via Local Processing

The CSC model [19] assumes that a global signal $y\in {ℝ}^N$ can be decomposed as

$$y=D\Gamma ={\displaystyle {\sum }_{i=1}^N{P}_i^T{D}_L{\mathsf{\alpha }}_i}$$

(4)

The matrix $D\in {ℝ}^{N\times NM}$ is a banded convolutional dictionary. This matrix consists of all shifted versions of a local dictionary $D_L\in {ℝ}^{n\times m}$ whose columns are atoms. L represents the initial of the word “local.” The global sparse vector $\Gamma \in {ℝ}^{NM}$ can be decomposed into N non-overlapping, m-dimensional local sparse vectors ${\left\{{\mathsf{\alpha }}_i\right\}}_{i=1}^N$, where ${\mathsf{\alpha }}_i\in {ℝ}^M$ are called needles. The operator $P_i^T$ places $D_L{\mathsf{\alpha }}_i$ in the ith position of the signal. The above formula can be used to obtain sparse coefficients using the basic pursuit problem, as shown below:

$$\underset{{\left\{{\mathsf{\alpha }}_i\right\}}_{i=1}^N}{\min }\frac{1}{2}{\Vert y-{\displaystyle \sum _{i=1}^N{P}_i^T{D}_L{\mathsf{\alpha }}_i}\Vert }_2^2+\lambda {\displaystyle \sum _{i=1}^N{\Vert {\mathsf{\alpha }}_i\Vert }_1}$$

(5)

Papyan et al. [23] proposed slice-based local processing and defined $s_i=D_L{\mathsf{\alpha }}_i$ as the ith slice. Unlike other existing works in signal and image processing, which train a dictionary in the Fourier domain, they defined the learning problem based on the constructed slices, and CDL was carried out in the signal domain via local processing. Through local processing in the original signal domain, the global problem was solved completely. The global signal $y_l$ can be rewritten as $y={\displaystyle {\sum }_{i=1}^N{P}_i^T{s}_i}$.

The ADMM algorithm is used to solve the problem of minimizing the following augmented Lagrangian problem:

$$\underset{D_L,{\left\{{\mathsf{\alpha }}_i\right\}}_{i=1}^N,{\left\{s_i\right\}}_{i=1}^N,{\left\{u_i\right\}}_{i=1}^N}{\min }\frac{1}{2}{\Vert y-{\displaystyle \sum _{i=1}^N{P}_i^T{s}_i}\Vert }_2^2+{\displaystyle \sum _{i=1}^N\left(\lambda {\Vert {\mathsf{\alpha }}_i\Vert }_1+\frac{\mathsf{\rho }}{2}{\Vert s_i-D_L{\mathsf{\alpha }}_i+u_i\Vert }_2^2\right)}$$

(6)

where ${\left\{u_i\right\}}_{i=1}^N$ are the dual variables that satisfy the given constraint. $\mathsf{\rho }$ is the Lagrangian coefficient. References [23,24] call this method the SBDL algorithm.

In the SBDL algorithm, the CDL problem is transformed into a traditional dictionary learning problem solved by the K-SVD algorithm [13] or other dictionary learning algorithms. The sparse coding process is based on the ADMM algorithm for updating, and each coding update is regarded as a least absolute shrinkage and selection operator (LASSO) problem, which is solved by the least-angle regression and shrinkage (LARS) algorithm. The ADMM is used to solve the coding problem, which increases the numbers of auxiliary variables, calculations and redundant iterations. In addition, different solving methods are used to update dictionary and codes, which makes it difficult to prove the convergence of the algorithm.

2.2. Forward–Backward Splitting

Splitting algorithms for convex optimization problems usually originate from the proximal point algorithm [32]. The proximal point algorithm is very general, and the results regarding its convergence affect many other algorithms. However, in practice, an iteration of the computational algorithm can be as tricky as the original problem. The strategy to solve this problem involves splitting approaches such as the Douglas–Rachford method, several primal-dual algorithms, and forward–backward splitting.

Forward–backward splitting schemes have been used to solve a variety of problems in recent years. For example, forward–backward splitting algorithms are used to solve normal problems [33], to find generalized Nash equilibrium [34], to solve linear constraint problems [35], or to analyze related function problems in Banach space [36,37]. In particular, it is appealing to generalize forward–backward splitting schemes to nonconvex problems. This is due to their simplicity and simpler formulations in some exceptional cases, such as the gradient projection method, where the backward step is a projection onto a set. The backward step of the forward–backward algorithm studied in [28] is the solution of a proximal term of a nonconvex function.

The goal of a forward–backward splitting framework is to solve the following forms of an optimization problem:

$$\underset{x}{\arg \min }f\left(x\right)+g\left(x\right)$$

(7)

However, when a large amount of data is processed, even if misestimation is not considered during processing, the processing result of the forward–backward splitting algorithm becomes inaccurate. In this paper, the original algorithm is improved to increase the accuracy and reduce the induced error.

2.3. Inertia Item

Polyak studied a multistep scheme for the accelerated gradient method in [30] and proposed the heavy-ball method for the first time. Unlike the usual gradient method, this approach adds an inertia term, which is calculated by the difference between the values obtained during the previous two iterations. Compared with the standard gradient method, this method can accelerate the convergence speed of the algorithm while keeping the cost of each iteration unchanged. In addition, this method can obtain the optimal convergence rate without additional knowledge.

The inertia term of the heavy-ball method was first applied to the minimization of differentiable convex functions, and then it was extended to subdifferential convex functions. Later, the heavy-ball method was used for forward–backward splitting [37,38]. Recently, it has been frequently applied to the minimization of nonconvex functions. In [39,40], the inertia term was introduced into the nonconvex optimization problem to improve the convergence speed. Refs. [31,41] used an inertia term to develop a nonconvex optimization CDL scheme. The authors of [31] aimed to minimize problems composed of differentiable (possibly nonconvex) and convex (possibly nondifferentiable) functions. The iPiano algorithm was proposed by combining forward–backward splitting with an inertia term. The global convergence of the function values and parameters was determined. In [41], an inertial version of the proximal alternating linearization minimization (PALM) algorithm was proposed, and its global convergence to the critical point of the input objective function was proven.

3. Proposed IPGM Algorithm

The CDL problem to be solved via local processing is given as follows:

$$\underset{{\left\{{\mathsf{\alpha }}_i\right\}}_{i=1}^N,D_L}{\arg \min }\frac{1}{2}{{\displaystyle \sum _l\Vert {\displaystyle \sum _{i=1}^N{P}_i^T{D}_L{\mathsf{\alpha }}_{l,i}-y_l}\Vert }}^2+{\lambda }_1{\displaystyle \sum _l\left({\displaystyle \sum _{i=1}^N{\Omega }_1\left({\mathsf{\alpha }}_{l,i}\right)}\right)}+{\lambda }_2{\Omega }_2\left(D_L\right)$$

(8)

$D_L$ is the local convolutional dictionary, which has n rows and m columns. ${\mathsf{\alpha }}_{l,i}$, which has m rows, is the sparse coding of each component i of each sample l. $P_i^T$, which has N rows and n columns, is the operator that puts $D_L{\mathsf{\alpha }}_{l,i}$ in the ith position and pads the rest of the entries with zeros. $y_l$ is the observed signal. $\Vert \cdot \Vert$ has different meanings; when the interior is a vector, $\Vert \cdot \Vert$ represents the ${\ell }_2$-norm, and when the interior is a matrix, $\Vert \cdot \Vert$ represents the Frobenius norm. ${\lambda }_1$ and ${\lambda }_2$ are hyperparameters. ${\Omega }_1$ is the sparse constraint imposed on the column vectors with the ${\ell }_0$ norm or ${\ell }_1$ norm, which is defined as follows: ${\Omega }_1\left(x\right)={\Vert x\Vert }_0$ or ${\Omega }_1\left(x\right)={\Vert x\Vert }_1$. ${\Omega }_2$ is the indicator function of unit-norm sphere.

To solve the above convex or nonconvex CDL optimization problem (8) via local processing, we use the inertial forward-backward splitting framework. An inertial forward-backward splitting algorithm via local processing called the inertial proximal gradient method (IPGM) algorithm is proposed.

To form the objective function of the proposed CDL optimization problem via local processing with the developed framework, based on the proposed Equation (8), supposing that $x=\left(D_{{}_L}^t,\left\{{\mathsf{\alpha }}_{l,i}^t\right\}\right)$, $\left\{x^t\right\}$ can be generated by iteratively implementing the following equation:

$$x^{t+1}\leftarrow prox\left((x^t-{\eta }_t\nabla f(x^t)+{\xi }_t(x^t-x^{t-1})\right)$$

(9)

where prox refers to the proximal mapping operation.

$f$ and $g$ are, respectively defined as follows:

$$f\left(D_L,\left\{{\mathsf{\alpha }}_{l,i}\right\}\right)=\frac{1}{2}{{\displaystyle \sum _l\Vert {\displaystyle \sum _{i=1}^N{P}_i^T{D}_L{\mathsf{\alpha }}_{l,i}-y_l}\Vert }}^2$$

(10)

$$g\left(D_L,\left\{{\mathsf{\alpha }}_{l,i}\right\}\right)={\displaystyle \sum _l\left({\displaystyle \sum _{i=1}^N{\lambda }_1{\Omega }_1\left({\mathsf{\alpha }}_{l,\mathrm{i}}\right)}\right)}+{\lambda }_2{\Omega }_2\left(D_L\right)$$

(11)

To generate the sequence $\left\{x^t\right\}=\left\{D_{{}_L}^t,\left\{{\mathsf{\alpha }}_{l,i}^t\right\}\right\}$ using iterative Equation (9) based on inertial forward-backward splitting, we first need to derive the gradient of $f$ and the proximal mapping of $g$.

Since $f$ is a function of the composite variable $x$ or $\left(D_{{}_L}^{},\left\{{\mathsf{\alpha }}_{l,i}^{}\right\}\right)$, we define the gradient of $f$ as follows:

$$\nabla f:=\left\{{\nabla }_{D_L}f,\left\{{\nabla }_{{\mathsf{\alpha }}_{l,i}}f\right\}\right\}$$

(12)

${\nabla }_{D_L}f$ and ${\nabla }_{{\mathsf{\alpha }}_{l,i}}f$ can be computed as follows:

$${\nabla }_{D_L}f={\displaystyle \sum _l{\displaystyle \sum _{i=1}^N\left(P_i\left({\displaystyle \sum _{i=1}^N{P}_i^T{D}_L{\mathsf{\alpha }}_{l,i}-y_l}\right){\mathsf{\alpha }}_{{}_{l,i}}^T\right)}}$$

(13)

$${\nabla }_{{\mathsf{\alpha }}_{a,b}}f=D_{{}_L}^T{P}_b\left({\displaystyle \sum _{i=1}^N{P}_i^T{D}_L{\mathsf{\alpha }}_{a,i}-y_a}\right)$$

(14)

To show the result of the descent process, we use an intermediate variable $\left(D_{{}_L}^{t+1/2},\left\{{\mathsf{\alpha }}_{l,i}^{t+1/2}\right\}\right)$:

$$D_L^{t+1/2}:=D_L^t-{\eta }_t{\nabla }_{D_L}f+{\xi }_t\left(D_L^t-D_L^{t-1}\right)$$

(15)

$${\mathsf{\alpha }}_{l,i}^{t+1/2}:={\mathsf{\alpha }}_{l,i}^t-{\eta }_t{\nabla }_{{\mathsf{\alpha }}_{l,i}}f+{\xi }_t\left({\mathsf{\alpha }}_{l,i}^t-{\mathsf{\alpha }}_{l,i}^{t-1}\right)$$

(16)

where ${\eta }_t$ is a step size or descent parameter and ${\xi }_t$ is an inertial parameter. Therefore, we compute the proximal mapping of $g$ at $\left(\left\{D_{{}_L}^{t+1/2}\right\},\left\{{\mathsf{\alpha }}_{l,i}^{t+1/2}\right\}\right)$ as follows:

$$\left(D_{{}_L}^{t+1},\left\{{\mathsf{\alpha }}_{l,i}^{t+1}\right\}\right)=pro{x}_{{\eta }_{t}g}\left(D_{{}_L}^{t+1/2},\left\{{\mathsf{\alpha }}_{l,i}^{t+1/2}\right\}\right)=(\left\{pro{x}_{{\eta }_t{\lambda }_2{\Omega }_2}\left(D_{{}_L}^{t+1/2}\right)\right\},\left\{pro{x}_{{\eta }_t{\lambda }_1{\Vert \cdot \Vert }_0}\left({\mathsf{\alpha }}_{l,\mathrm{i}}^{t+1/2}\right)\right\})$$

(17)

The function prox refers to the proximal mapping operation, which is defined with the following form:

$$pro{x}_{\eta p}(u)=\underset{v}{\arg \min }\frac{1}{2\eta }{\Vert u-v\Vert }^2+p(v)$$

(18)

Furthermore, in the inertial forward-backward splitting framework, to accelerate the convergence rate, we propose using ${\eta }_t$ to fit the local curvature of $f$ in each iteration and capturing the local curvature $f$ by estimating the local Lipschitz constant of the CDL problem via local processing. This constant is defined as follows:

$$L_t=\frac{\sqrt{{\Vert \nabla f\left(x^t\right)-\nabla f\left(x^{t-1}\right)\Vert }^2}}{\sqrt{{\Vert x^t-x^{t-1}\Vert }^2}}$$

(19)

However, the direct derivation of ${\eta }_t$ from $L_t$ does not satisfy the convergence requirement. It is suggested to insert a backtracking scheme in the inertial forward-backward splitting framework to restore its convergence.

We introduce a sequence $\left\{{\tau }_t\right\}$, where ${\tau }_t>1$ holds for all $t$, to represent the adaptive parameter such that the sequence $\left\{x^t\right\}$ satisfies the following equation:

$$f\left(x^{t+1}\right)\le f\left(x^t\right)+\langle \nabla f\left(x^t\right),x^{t+1}-x^t\rangle +\frac{{\tau }_t{L}_t}{2}{\Vert x^{t+1}-x^t\Vert }^2$$

(20)

We assume that each step size ${\eta }_t$ is maintained by ${\tau }_t$ using the inequality $0<{\eta }_t<\left(1/\left({\tau }_t{L}_t\right)\right)$ and that (20) can be reformulated as follows:

$$\begin{array}{l}f\left(x^{t+1}\right)=f\left(D_L^{t+1},\left\{{\mathsf{\alpha }}_{l,i}^{t+1}\right\}\right)\le f\left(x^t\right)+\langle \nabla f\left(x^t\right),x^{t+1}-x^t\rangle +\frac{1}{2{\eta }_t}{\Vert x^{t+1}-x^t\Vert }^2\\ =f\left(D_L^t,\left\{{\mathsf{\alpha }}_{l,i}^t\right\}\right)+{\displaystyle \sum _{p,q}\left[{\mathsf{\alpha }}_{{}_{p,q}}^{t+1}\odot \left({\nabla }_{{\mathsf{\alpha }}_{p,q}}f\right)\left(D_L^t,\left\{{\mathsf{\alpha }}_{l,i}^t\right\}\right)+\frac{1}{2{\eta }_t}{\Gamma }_r\left({\left({\mathsf{\alpha }}_{p,q}^{t+1}-{\mathsf{\alpha }}_{p,q}^t\right)}^T\left({\mathsf{\alpha }}_{p,q}^{t+1}-{\mathsf{\alpha }}_{p,q}^t\right)\right)\right]}\\ +D_{{}_L}^{t+1}\odot \left({\nabla }_{D_{{}_L}}f\right)\left(D_L^t,\left\{{\mathsf{\alpha }}_{l,i}^t\right\}\right)+\frac{1}{2{\eta }_t}{\Gamma }_r\left({\left(D_{{}_L}^{t+1}-D_{{}_L}^t\right)}^T\left(D_{{}_L}^{t+1}-D_{{}_L}^t\right)\right)\end{array}$$

(21)

where $\odot$ stands for the Hadamard product.

We present the IPGM algorithm that solves the proposed CDL problem via local processing using inertial forward–backward splitting in Algorithm 1.

Algorithm 1 The IPGM algorithm for solving the proposed CDL problem via local processing

Input: initial dictionary $D_{{}_L}^0$; initial dictionary $D_L^{-1}=D_L^0$; initial coefficients map $\left\{{\mathsf{\alpha }}_{l,i}^0\right\}$; initial coefficient maps $\left\{{\mathsf{\alpha }}_{l,i}^{-1}\right\}=\left\{{\mathsf{\alpha }}_{l,i}^0\right\}$; adaptive parameter $\tau >1$; choose $c_1>0$ close to 0; inertial parameter ${\xi }_0\in \left[0,\frac{1}{2}\right)$; and descent parameter $1<{\eta }_0<\left(1-2{\xi }_0\right)/\tau L$, where $L$ is the Lipschitz constant of $\nabla f$. Output: learned dictionary $D_{{}_L}^t$; learned coefficients $\left\{{\mathsf{\alpha }}_{l,i}^t\right\}$; 1: Initialization: $t\leftarrow 0$ 2: repeat 3:               if $t<2$, then 4:                               ${\eta }_t\leftarrow {\eta }_0$, ${\xi }_t\leftarrow {\xi }_0$ 5:                else 6:                               compute $L_t$ using Equation (19) 7:                               ${\eta }_t\leftarrow \frac{1}{L_t}$ 8:               end if 9:                               $k\leftarrow 0$ 10:               repeat 11:                               compute ${\nabla }_{D_L}f$ and ${\nabla }_{{\mathsf{\alpha }}_{l,i}}f$ using Equations (13) and (14) 12:                               compute $\left(D_{{}_L}^{t+1/2},\left\{{\mathsf{\alpha }}_{l,i}^{t+1/2}\right\}\right)$ using Equations (15) and (16) 13:                               update $\left(D_{{}_L}^{t+1},\left\{{\mathsf{\alpha }}_{l,i}^{t+1}\right\}\right)$ using Equation (17) 14:                               ${\tau }_t\leftarrow {\tau }^k$ 15:                               ${\eta }_t\leftarrow \frac{{\eta }_t}{\tau }$ 16:                               $k\leftarrow k+1$ 17:               until Equation (21) is satisfied 18:                               ${\eta }_t\ge c_1$, ${\xi }_t\ge 0$ 19:                               ${\delta }_t:=\frac{1}{2{\eta }_t}-\frac{{\tau }_t{L}_t}{2}-\frac{{\xi }_t}{2{\eta }_t}$ 20:                                 ${\gamma }_t:=\frac{1}{2{\eta }_t}-\frac{{\tau }_t{L}_t}{2}-\frac{{\xi }_t}{{\eta }_t}$ 21:       until ${\delta }_t\ge {\gamma }_t\ge c_1$ and ${\left({\delta }_t\right)}_{t=0}^{\infty }$ is monotonically decreasing 22:       $t\leftarrow t+1$ 23: until convergence 24: return $\left\{\left\{D_{{}_L}^t\right\},\left\{{\mathsf{\alpha }}_{l,i}^t\right\}\right\}$

4. Computational Complexity Analysis

The computational complexity of our algorithms is discussed as follows. We assume the following. $N$ is the signal dimension. $I$ is the number of signals. $n$ is the patch size. $m$ is the number of filters. $k$ is the maximal number of nonzeros in a needle ${\mathsf{\alpha }}_{l,i}$, which is very sparse $k\ll m$. $C$ is the number of backtracking loops, which is usually very small. The convolution is performed by the local processing operation. The operation $P_i^T$ is only an operator of a column to an image.

The dominant computational complexity of the signal reconstruction performed by the local processing operation is approximately $O\left(INnk\right)$. The computation of the signal residual requires $O\left(IN\right)$ operations. Neglecting some minor factors, the computational complexity of the gradient of needles and dictionary performed by the local processing operation is effectively $O\left(INnm\right)$ and $O\left(INnk\right)$, respectively. In addition, the inertia term includes an addition and a multiplication for each iteration which is negligible.

5. Proof of the Convergence of Algorithm 1

Before describing the convergence theorem, let us analyze the relevant lemmas and hypotheses.

Lemma 1.

(The abstract convergence result for KL functions) This convergence result is based on three abstract conditions for a sequence${\left(z^t\right)}_{t\in \mathbb{N}}:={\left(x^t,x^{t-1}\right)}_{t\in \mathbb{N}}\in {ℝ}^{2N}$, $x^t\in {ℝ}^N,x^{t-1}\in {ℝ}^N$. We fix two positive constants a > 0 and b > 0 and consider a proper lower semicontinuous function$F:{ℝ}^{2N}\to ℝ\cup \{\infty \}$. Then, the conditions for${\left(z^t\right)}_{t\in \mathbb{N}}$are as follows:

(H1) For each $k\in \mathbb{N}$, it holds that

$$F\left(z^{k+1}\right)+a\Vert x^k-x^{k-1}\Vert \le F\left(z^k\right)$$

(H2) For each $k\in \mathbb{N}$, there exists a $w^{k+1}\in \partial F\left(w^{k+1}\right)$ such that

$$\Vert w^{k+1}\Vert \le \frac{b}{2}\left(\Vert x^k-x^{k-1}\Vert +\Vert x^{k+1}-x^k\Vert \right)$$

(H3) There exists a subsequence ${\left(z^{k_j}\right)}_{j\in \mathbb{N}}$ such that

$$z^{k_j}\to \tilde{z} and F\left(z^{k_j}\right)\to F\left(\tilde{z}\right) as j\to \infty$$

When $F$ is a KL function and H1, H2, and H3 are satisfied, $F$ satisfies the convergence result.

Lemma 2.

For all$n\ge 0$, ${\delta }_t\ge {\gamma }_t$, ${\xi }_t\in \left[0,\frac{1}{2}\right)$, and${\eta }_t\le \left(1-2{\xi }_t\right)/{\tau }_t{L}_t$. Furthermore, given$L_t>0$, there exists a pair of parameters${\eta }_t$and${\xi }_t$such that${\left({\delta }_t\right)}_{t=0}^{\infty }$is monotonically decreasing.

Proof of Lemma 2.

By the algorithmic requirements,

$${\delta }_t=\frac{1}{2{\eta }_t}-\frac{{\tau }_t{L}_t}{2}-\frac{{\xi }_t}{2{\eta }_t}\ge \frac{1}{2{\eta }_t}-\frac{{\tau }_t{L}_t}{2}-\frac{{\xi }_t}{{\eta }_t}={\gamma }_t\ge 0$$

(22)

The upper bounds for ${\xi }_t$ and ${\eta }_t$ are obtained by rearranging ${\gamma }_t\ge c_1$ to ${\xi }_t\le \left(1-{\eta }_t{\tau }_t{L}_t-2{\eta }_t{c}_1\right)/2$ and ${\eta }_t\le \left(1-2{\xi }_t\right)/\left(2c_1+{\tau }_t{L}_t\right)$, respectively.

The last statement follows by incorporating the descent property of ${\delta }_t$. Let ${\delta }_{-1}\ge c_1$ be chosen initially. Then, the decent property of ${\left({\delta }_t\right)}_{t=0}^{\infty }$ requires one of the equivalent statements below

$${\delta }_{t-1}\ge {\delta }_t\iff {\delta }_{t-1}\ge \frac{1}{2{\eta }_t}-\frac{{\tau }_t{L}_t}{2}-\frac{{\xi }_t}{2{\eta }_t}\iff {\eta }_t\ge \left(1-{\xi }_t\right)/\left(2{\delta }_{t-1}+{\tau }_t{L}_t\right)$$

(23)

to be true. An upper bound on ${\eta }_t$ is obtained by

$${\gamma }_t\ge c_1\iff {\eta }_t\le \frac{1-2{\xi }_t}{2c_1+{\tau }_t{L}_t}$$

(24)

Consider the condition for a nonnegative gap between the upper and lower bounds of ${\eta }_t$:

$$\frac{1-2{\xi }_t}{2c_1+{\tau }_t{L}_t}-\frac{1-{\xi }_t}{2{\delta }_{n-1}+{\tau }_t{L}_t}\ge 0\Rightarrow \frac{2{\delta }_{n-1}+{\tau }_t{L}_t}{2c_1+{\tau }_t{L}_t}\ge \frac{1-{\xi }_t}{1-2{\xi }_t}$$

(25)

Defining $b:=\left(2\delta +{\tau }_t{L}_t\right)/\left(2c_1+{\tau }_t{L}_t\right)\ge 1$, it is easily verified that there exists ${\xi }_t\in \left[0,\frac{1}{2}\right)$ satisfying the equivalent condition

$$\frac{b-1}{2b-1}\ge {\xi }_t$$

(26)

As a consequence, the existence of a feasible ${\eta }_t$ follows, and the decent property for ${\delta }_t$ holds. □

Proposition 1.

(a) The sequence${\left(H_{{\delta }_t}\left(x^t,x^{t-1}\right)\right)}_{n=0}^{\infty }$is monotonically decreasing and thus convergent. In particular, it holds that

$$H_{{\delta }_{t+1}}\left(x^{t+1},x^t\right)\le H_{{\delta }_t}\left(x^t,x^{t-1}\right)-{\gamma }_t{\Delta }_t^2$$

(27)

(b) It holds that${\displaystyle {\sum }_{n=0}^{\infty }{\Delta }_t^2<\infty }$, and thus,${\lim }_{t\to \infty }{\Delta }_t=0$.

Proof of Proposition 1.

(a) For a more convenient notation, we abbreviate $h=f+g$ in the IPGM algorithm for the nonconvex nonsmooth CDL problem via local processing; this does not mean that the value of the function drops, so we construct a function $H_{\delta }\left(x,y\right)=h\left(x\right)+\delta {\Vert x-y\Vert }_{}^2,\delta \in ℝ$. Note that for $x=y$, $H_{\delta }\left(x,y\right)=h\left(x\right)$.

We show below that $H_{\delta }\left(x,y\right)$ satisfies the convergence requirement. In the inertial forward–backward splitting framework, the sequence $\left\{x^t\right\}$ is generated by iteratively implementing Equation (9).

Notably, in the IPGM algorithm, the function $f$ satisfies the adaptive descent formula in (20).

The function $g$ is a nonconvex, nonsmooth, and proper closed function. Combining the iterative mapping of the sequence $\left\{x^t\right\}$ in the inertial forward–backward splitting framework and the definition of the proximal operator, we can obtain

$$\begin{array}{l}\frac{1}{2{\eta }_t}{\Vert x^{t+1}-\left(x^t-{\eta }_t\nabla f\left(x^t\right)+{\xi }_t\left(x^t-x^{t-1}\right)\right)\Vert }^2+g\left(x^{t+1}\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \le \frac{1}{2{\eta }_t}{\Vert x^t-\left(x^t-{\eta }_t\nabla f\left(x^t\right)+{\xi }_t\left(x^t-x^{t-1}\right)\right)\Vert }^2+g\left(x^t\right)\end{array}$$

(28)

Thus, this simplifies to

$$g\left(x^{t+1}\right)+\frac{1}{2{\eta }_t}{\Vert x^{t+1}-x^t\Vert }^2+\frac{1}{{\eta }_t}\langle x^{t+1}-x^t,{\eta }_t\nabla f\left(x^t\right)-{\xi }_t\left(x^t-x^{t-1}\right)\rangle \le g\left(x^t\right)$$

(29)

Now, using (20) and (29) by summing both inequalities, it follows that

$$\begin{array}{l}h(x^{t+1})\le h(x^t)-\left(\frac{1}{2{\eta }_t}-\frac{{\tau }_t{L}_t}{2}\right){\Vert x^{t+1}-x^t\Vert }^2+\frac{{\xi }_t}{{\eta }_t}\langle x^{t+1}-x^t,x^t-x^{t-1}\rangle \\ \hspace{1em}\hspace{1em}\hspace{1em}\le h(x^t)-\left(\frac{1}{2{\eta }_t}-\frac{{\tau }_t{L}_t}{2}-\frac{{\xi }_t}{2{\eta }_t}\right){\Vert x^{t+1}-x^t\Vert }^2+\frac{{\xi }_t}{2{\eta }_t}{\Vert x^t-x^{t-1}\Vert }^2\end{array}$$

(30)

where the second line follows from $2\langle A,B\rangle \le {\Vert A\Vert }_2^2+{\Vert B\Vert }_2^2$ for vectors $A,B\in {ℝ}^N$. Then, a simple rearrangement of the terms yields:

$$h\left(x^{t+1}\right)+{\delta }_t{\Vert x^{t+1}-x^t\Vert }^2\le h\left(x^t\right)+{\delta }_t{\Vert x^t-x^{t-1}\Vert }^2-{\gamma }_t{\Vert x^t-x^{t-1}\Vert }^2$$

(31)

which establishes (27) as ${\delta }_t$ is monotonically decreasing. The sequence ${\left(H_{{\delta }_t}\left(x^t,x^{t-1}\right)\right)}_{n=0}^{\infty }$ monotonically decreases if and only if ${\gamma }_t\ge 0$, which is confirmed by the algorithmic requirements. By assumption, $h$ is bounded from below by some constant $\underset{\_}{h}>-\infty$; hence, ${\left(H_{{\delta }_t}\left(x^t,x^{t-1}\right)\right)}_{n=0}^{\infty }$ converges.

(b) Summing (27) from $t=0,\cdots ,T$ yields (note that $H_{\delta }{}_t\left(x^0,x^{-1}\right)=h(x^0)$)

$${\displaystyle \sum _{t=0}^T{r}_t{\Delta }_t^2}\le {\displaystyle \sum _{t=0}^T{H}_{{\delta }_t}\left(x^t,x^{t-1}\right)-}{H}_{{\delta }_{t+1}}\left(x^{t+1},x^t\right)=h\left(x^0\right)-H_{{\delta }_{T+1}}\left(x^{T+1},x^T\right)\le h\left(x^0\right)-\underset{\_}{h}<\infty$$

(32)

Letting $\mathrm{T}$ tend to $\infty$, it can be seen from (32) that $\underset{t\to T}{\lim }{\gamma }_{\mathrm{t}}{\Delta }_t^2=0$ since ${\gamma }_T\ge c_1>0$ implies the above statement. □

Proposition 2.

(a) The sequence${(h(x^t))}_{t=0}^{\infty }$converges.

(b) There exists a convergent subsequence${(x^{t_k})}_{k=0}^{\infty }.$

(c) Any limit point$x^{\ast }:=li{m}_{k\to \infty }{x}^{t_k}$is a critical point of$h(x)$, and$h(x^{t_k})\to h(x^{\ast })$as$k\to \infty$.

Proof of Proposition 2.

(a) This follows from the squeeze theorem as for all $t\ge 0$, the following holds:

$$H_{-{\delta }_t}\left(x^t,x^{t-1}\right)\le h\left(x^t\right)\le H_{{\delta }_t}\left(x^t,x^{t-1}\right)$$

(33)

and due to Propositions 1(a) and (b),

$$\underset{t\to \infty }{\lim }{H}_{-{\delta }_t}\left(x^t,x^{t-1}\right)=\underset{t\to \infty }{\lim }{H}_{{\delta }_t}\left(x^t,x^{t-1}\right)-2{\delta }_t{\Delta }_n^2=\underset{t\to \infty }{\lim }{H}_{{\delta }_t}\left(x^t,x^{t-1}\right)$$

(34)

(b) By Proposition 1(a) and the fact that $H_{{\delta }_0}\left(x^0,x^{-1}\right)=h\left(x^0\right)$, it is clear that the whole sequence ${(x^t)}_{t=0}^{\infty }$ is contained in the level set $\left\{x\in {ℝ}^N:h\le h\left(x\right)\le h\left(x^0\right)\right\}$, which is bounded due to the coercivity of $h$ and because $h={\inf }_{x\in {ℝ}^N}h\left(x\right)>-\infty$. Using the Bolzano-Weierstrass theorem, we deduce the existence of a convergent subsequence ${(x^{t_k})}_{k=0}^{\infty }.$

(c) To show that each limit point $x^{\ast }=li{m}_{j\to \infty }{x}^{t_j}$ is a critical point of $h\left(x\right)$, we recall that the subdifferential is closed. We define

$${\varsigma }_j:=\frac{x^{t_j}-x^{t_j+1}}{{\eta }_{t_j}}-\nabla f\left(x^{t_j}\right)+\frac{{\xi }_{t_j}}{{\eta }_{t_j}}\left(x^{t_j}-x^{t_j-1}\right)+\nabla f\left(x^{t_j+1}\right).$$

(35)

Then, the sequence $\left(x^{t_j},{\varsigma }_j\right)\in Graph\left(\partial h\right):=\{\left(x,\varsigma \right)\in {ℝ}^N\times {ℝ}^N|\varsigma \in \partial h\left(x\right)\}.$ Furthermore, it holds that $x^{\ast }=li{m}_{j\to \infty }{x}^{t_j}$, and due to Proposition 1(b), $\nabla f$ is Lipschitz-continuous, and

$$\Vert {\varsigma }^j-0\Vert \le \frac{1}{{\eta }_{t_j}}{\Delta }_{n_j+1}+\frac{{\xi }_{t_j}}{{\eta }_{t_j}}{\Delta }_{n_j}+\Vert \nabla f\left(x^{t_j+1}\right)-\nabla f\left(x^{t_j}\right)\Vert$$

(36)

It holds that $li{m}_{j\to \infty }{\varsigma }^j=0$. It remains to be shown that $li{m}_{j\to \infty }h\left(x^{t_j}\right)=h(x^{\ast }).$ By the closure property of the subdifferential $\partial h$, $(x^{\ast },0)\in Graph\left(\partial h\right)$, which means that $x^{\ast }$ is a critical point of $h$.

According to the iterative mapping (Equation (9)) of the sequence $\left\{x^t\right\}$ in the inertial forward-backward splitting framework, we can obtain:

$$\begin{array}{l}\frac{1}{2{\eta }_t}{\Vert x^{t+1}-\left(x^t-{\eta }_t\nabla f\left(x^t\right)+{\xi }_t\left(x^t-x^{t-1}\right)\right)\Vert }^2+g\left(x^{t+1}\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\le \frac{1}{2{\eta }_t}{\Vert x-\left(x^t-{\eta }_t\nabla f\left(x^t\right)+{\xi }_t\left(x^t-x^{t-1}\right)\right)\Vert }^2+g\left(x\right)\end{array}$$

(37)

which implies that

$$\begin{array}{l}g\left(x^{t_j+1}\right)+\frac{1}{{\eta }_{t_j}}\langle {\eta }_{t_j}\nabla f\left(x^{t_j}\right)-{\xi }_t\left(x^{t_j}-x^{t_j-1}\right),x^{t_j+1}-x\rangle \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\frac{1}{2{\eta }_{t_j}}\left({\Vert x^{t_j+1}-x^{t_j}\Vert }^2-{\Vert x-x^{t_j}\Vert }^2\right)\le g\left(x\right)\end{array}$$

(38)

Proposition 1(b) and the boundedness of ${\eta }_{t_j}\nabla f\left(x^{t_j}\right)-{\xi }_t\left(x^{t_j}-x^{t_j-1}\right)$ yield that $limsu{p}_{j\to \infty }g\left(x^{t_j}\right)\le g\left(x\right)$. Invoking the lower semicontinuity of $g$ yields $li{m}_{j\to \infty }g\left(x^{t_j}\right)=g(x^{\ast })$.

Moreover, $f$ is differentiable and continuous, thus, $li{m}_{j\to \infty }f\left(x^{t_j}\right)=f(x^{\ast })$. We imply that $li{m}_{j\to \infty }h\left(x^{t_j}\right)=h(x^{\ast }).$

Now, using Lemma 1, we can verify the convergence of the sequence ${\left(x^t\right)}_{t\in \mathbb{N}}$ generated by Algorithm 1.  □

Theorem 1

(Convergence of the IPGM algorithm to a critical point). Let${\left(x^t\right)}_{t\in \mathbb{N}}$be generated by Algorithm 1, and let${\delta }_t=\delta$for all$t\in \mathbb{N}$. Then, the sequence${\left(x^{t+1},x^t\right)}_{t\in \mathbb{N}}$satisfies H1, H2, and H3 for the function$H_{\delta }:{ℝ}^{2N}\to ℝ\cup \{\infty \},\left(x,y\right)\mapsto h\left(x\right)+\delta {\Vert x-y\Vert }_2^2.$Moreover, if the sequence possesses the KL property at a cluster point, then the sequence has a finite length and is a critical point$H_{\delta }$; hence,$x^{\ast }$is a critical point of h.

Proof of Proposition Theorem 1.

First, we prove that the function has the KL property.$h$ is a semialgebraic function. If we consider that ${\Vert x-y\Vert }^2={\displaystyle \sum _{i=1}^N{\left[x\left(i\right)-y\left(i\right)\right]}^2}$, then ${\Vert x-y\Vert }^2$ is a polynomial function. Therefore, it is semialgebraic, so $\delta {\Vert x-y\Vert }_{}^2$ is a semialgebraic function. $H_{\delta }\left(x,y\right)$ is semialgebraic and has the KL property.

Next, we verify that Assumptions H1, H2, and H3 are satisfied.

Condition H1 is proven in Proposition 1(a) with $a=c_2\le {\gamma }_t.$ To prove Condition H2, consider $w^{t+1}:={(w_x^{t+1},w_y^{t+1})}^T\in \partial H_{\delta }\left(x^{t+1},x^t\right)$ with $w_x^{t+1}\in \partial g\left(x^{t+1}\right)+\nabla f\left(x^{t+1}\right)+2\delta (x^{t+1}-x^t)$ and $w_y^{t+1}=-2\delta (x^{t+1}-x^t)$. The Lipschitz continuity of $\nabla f$ and the use of (9) to specify an element from $\partial g\left(x^{t+1}\right)$ imply that

$$\begin{array}{l}\Vert w^{t+1}\Vert \le \Vert w_x^{t+1}\Vert +\Vert w_y^{t+1}\Vert \\ \hspace{1em}\hspace{1em} \le \Vert \nabla f\left(x^{t+1}\right)-\nabla f\left(x^t\right)\Vert +\left(\frac{1}{{\eta }_t}+4\delta \right)\Vert x^{t+1}-x^t\Vert +\frac{{\xi }_t}{{\eta }_t}\Vert x^t-x^{t-1}\Vert \\ \hspace{1em}\hspace{1em} \le \frac{1}{{\eta }_t}\left({\eta }_t{\tau }_t{L}_t+1+4{\eta }_t\delta \right)\Vert x^{t+1}-x^t\Vert +\frac{1}{{\eta }_t}{\xi }_t\Vert x^t-x^{t-1}\Vert \end{array}$$

(39)

As ${\delta }_t\ge {\gamma }_t\ge c_1>0$, we can obtain ${\eta }_t{\tau }_t{L}_t\le 1-2{\xi }_t\le 1$, and $\delta {\eta }_t=\frac{1}{2}-\frac{{\eta }_t{\tau }_t{L}_t}{2}-\frac{{\xi }_t}{2}\le \frac{1}{2}$. Setting $b=\frac{4}{c_1}$ verifies condition H2, i.e., $\Vert w^{t+1}\Vert \le b\left(\Vert x^t-x^{t-1}\Vert +\Vert x^{t+1}-x^t\Vert \right)$. Condition H2 is proven.

In Proposition 2(c), it is proven that there exists a subsequence ${\left(x^{t_j+1}\right)}_{j\in N}$ of ${\left(x^t\right)}_{t\in N}$ such that $\underset{j\to \infty }{lim}h\left(x^{t_j+1}\right)=h\left(x^{\ast }\right)$. The following corollary uses the fact that semialgebraic functions have the KL property. Proposition 1(b) shows that $\Vert x^{t+1}-x^t\Vert \to 0$ as $t\to \infty$; hence, $li{m}_{j\to \infty }{x}_{j+1}^t=x^{\ast }$. As the term $\delta {\Vert x-y\Vert }^2$ is continuous in x and y, we deduce that

$$\underset{j\to \infty }{lim}H(x^{t_{j+1}},x^{t_j})=\underset{j\to \infty }{lim}h(x^{t_j+1})+\delta \Vert x^{t_j+1}-x^{t_j}\Vert =H(x^{\ast },x^{\ast })=h(x^{\ast })$$

(40)

Therefore, Condition H3 is proven.

Now, Theorem 1 concludes the proof. □

Remark 1.

IPGM algorithm is convergent under nonconvex optimization. It is easy to prove the convergence under convex optimization constraints, which is equivalent to a special case of the iPiano algorithm for the CDL problem.

6. Convergence Rate

We prove a global $O\left(1/k\right)$ convergence rate for ${\Vert x^{t+1}-x^t\Vert }^2$. We first define the error ${\mu }_N$ to be the smallest squared ${\ell }_2$ norm value of successive iterations:

$${\mu }_N:=\underset{0\le n\le N}{min}{\Vert x^t-x^{t-1}\Vert }^2$$

(41)

Theorem 2.

Algorithm 1 guarantees that for all N ≥ 0,

$${\mu }_N\le c_1^{-1}\frac{h\left(x^0\right)-h}{N+1}$$

(42)

Proof of Proposition Theorem 2.

In view of Proposition 1(a) and the definition of ${\gamma }_N$ in (21), summing both sides of (26) for $n=0,...,N$ and using the fact that ${\delta }_N>0$ from (21), we obtain

$$h\le h\left(x^0\right)-{\displaystyle \sum _{n=0}^N{\gamma }_n{\Vert x^t-x^{t-1}\Vert }^2}\le h\left(x^0\right)-\left(N+1\right)\underset{0\le n\le N}{min}{\gamma }_n{\mu }_N.$$

(43)

As ${\gamma }_n>c_1,$ a simple rearrangement concludes the proof. □

7. Experimental Results and Analysis

In this section, we compare the performance of the proposed IPGM with that of various existing methods with respect to solving CDL problems.

7.1. Experimental Setup

7.1.1. List of Compared Methods

• SBDL denotes the method proposed in [23] based on convex optimization with an ${\ell }_1$ norm-based sparsity-inducing function that uses slices via local processing.
• Local block coordinate descent (LoBCoD) denotes the method proposed in [24] based on convex optimization with an ${\ell }_1$ norm-based sparsity-inducing function that utilizes needle-based representation via local processing.
• The PGM denotes the method proposed in [1] based on convex optimization with an ${\ell }_1$ norm-based sparsity-inducing function that uses the Fast Fourier transform (FFT) operation.
• The IPGM denotes the method proposed in Algorithm 1 in this paper uses needle-based representation via local processing.

7.1.2. Parameter Settings

The parameter settings used for the comparison methods are described as follows:

• In SBDL [23], the model parameter is $\lambda$, and its initial value is 1. The maximum nonzero value of the LARS algorithm is 5. The filter size of the dictionary is 11*11*100. In addition, special techniques are used in the first few iterations of the dictionary learning process.
• In LoBCoD [24], the model parameter $\lambda$ is initialized to 1. The maximum nonzero value of the LARS algorithm is 5. The filter format of the dictionary is 11*11*100. In addition, the dictionary learning process first carries out 40 iterations and then uses the proximal gradient descent algorithm.
• In the PGM [1] model, the Lipschitz constant L is set to 1000 for fruit, and $\lambda =1$.
• In the IPGM, the model parameter $L_t$ is set at 1000 for fruit, and $\lambda =1$.

7.1.3. Implementation Details

The computations of the SBDL, LoBCoD, and PGM algorithms, as well as the IPGM algorithm, are performed using a PC with an Intel i5 CPU and 12 GB of memory.

7.2. Efficiency of Dictionary Learning Algorithms

7.2.1. Motivation and Evaluation Metrics

The efficiency and stability of objective function are the most important criteria for evaluating a numerical optimization algorithm. The sequence generated by an efficient dictionary learning algorithm converges quickly to the corresponding clustering point with a lower function value. The relevant evaluation indicators are listed below.

• Final Value of $f+g$ and the Algorithmic Convergence Properties: In the $f+g$ minimization-based analysis method, we use the function value of the generated sequence to evaluate its optimization efficiency and convergence.
• Computing Time: We compare the computing times of the SBDL, LoBCoD, PGM, and IPGM algorithms and compare the computing efficiency.

7.2.2. Training Data and Experimental Settings

A set of generic grayscale natural images is used as the training data in this experiment. The set includes 10 fruit images (100 × 100 pixels). It is decomposed into high- and low-frequency components. The means are subtracted. These images are normalized to a range of 0 to 1. The dictionary includes 100 elements, the sizes of which are 11 × 11. One thousand iterations are carried out for each method. The dataset comes from SBDL [23] and LOBCOD [24], and the experimental process refers to the experimental process in references SBDL and LOBCOD. Similar to references [23,24], there is a similar trend in the experimental results under different atomic numbers or different sizes. We give the final objection function value representing the convergence properties and computing time representing convergence efficiency using only the conditions of 100 atoms with a size of 11 × 11 from the fruit grayscale natural images. Each of the average results represents the four trials using different initial dictionaries, from which the initial coefficients were derived.

7.2.3. Results

The experimental results are described as follows. The set of training data described above is used to derive the results. Figure 1 shows the functional values of all compared methods in each iteration of the dictionary learning procedure. In addition,

Table 1. The experimental results obtained by all compared methods after 1000 iterations with the grayscale natural fruit images.

Fruit	SBDL	LoBCoD	PGM	IPGM
Obj	9.545 × 103	9.857 × 103	1.048 × 104	9.449 × 103
Data consistent	3.067 × 103	3.136 × 103	3.528 × 103	3.035 × 103
Regularization	6.478 × 103	6.722 × 103	6.957 × 103	6.415 × 103
Sparsity	0.146%	0.141%	1.696%	0.145%
Time (s)	1110.134	1879.400	1955.365	1222.722
PSNR (dB)	29.348	29.252	28.799	29.438

shows the experimental data yielded by all compared methods over 1000 iterations. Figure 2 and Figure 3 correspond to the trained dictionaries and reconstructed images of all compared methods after 1000 iterations, respectively.

7.2.4. Discussion

A discussion of the experimental results is as follows. Table 1 shows the experimental results produced by the four algorithms for the fruit dataset, which shows that the IPGM method yields better results in terms of performances. In terms of the objective function values of the four algorithms, the IPGM algorithm is lowest with 9.449 × 10³. In terms of sparsity, except for the poor sparsity of the PGM, there is little difference among the other algorithms. The time consumption of our algorithm is 122.722 s, which is better than that of LoBCoD and the PGM but worse than that of SBDL. The SBDL algorithm calls for the C++ function in the MATLAB program, so it requires less time than the IPGM method. The peak signal-to-noise ratio (PSNR) of the IPGM algorithm is best with 29.438 dB, which is the highest value among the four algorithms. After analyzing Figure 3 above, we see that compared with other algorithms, the IPGM algorithm provides the reconstructed image with the clearest texture details. Generally, the IPGM algorithm has the best performance and the best effect when compared with other algorithms.

7.3. Ablation Experiment

7.3.1. Motivation and Evaluation Metrics

In the proposed IPGM algorithm, the inertia parameters are assigned directly, and they are equivalent to fixed values. Nevertheless, the magnitudes of the inertia parameters have specific impacts on the performance of the algorithm. Therefore, this section uses ablation experiments to analyze the performance of IPGM under different inertia parameters and the performance of IPGM under different optimization conditions.

7.3.2. Training Data and Experimental Settings

The experimental setting of the ablation experiment is consistent with that of the experiment in the previous section. The fruit dataset is also used as the ablation experiment dataset. The dataset comes from SBDL [23] and LOBCOD [24], and the experimental process refers to the experimental process in references SBDL and LOBCOD. We experiment on IPGM algorithms with different inertia parameters under convex optimization constraints and nonconvex optimization constraints.

7.3.3. Results

A discussion of the experimental results is as follows.

Table 2. Objective function values and PSNR values of IPGM models with different inertia parameters under the constraint of nonconvex optimization.

	0.1	0.2	0.3	0.4
Obj	1.072 × 104	1.059 × 104	1.097 × 104	1.167 × 104
PSNR	31.960	32.139	32.335	32.734

shows the objective function and PSNR values obtained with different inertia parameters under the constraint of nonconvex optimization.

Table 3. Objective function values and PSNR values of IPGM models with different inertia parameters under the constraint of convex optimization.

	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
Obj	9.737 × 103	9.642 × 103	9.604 × 103	9.569 × 103	9.565 × 103	9.543 × 103	9.516 × 103	9.484 × 103	9.449 × 103
PSNR	29.350	29.394	29.401	29.427	29.414	29.422	29.433	29.431	29.438

shows the objective function and PSNR values obtained with different inertia parameters under the constraint of convex optimization.

7.3.4. Discussion

• The above table summarizes the performance of the IPGM algorithm under different inertia parameter settings. Under nonconvex optimization, with the increase in the inertia parameter, the objective function value decreases, and the PSNR value increases. The effect is the best when the inertia parameter is 0.4. The objective function value is 1.167 × 10⁴, and the PSNR value is 32.734. Under convex optimization, with the increase in the inertia parameter, the value of the objective function decreases, and the value of PSNR increases. The effect is the best when the inertia parameter is 0.9. At this time, the value of the objective function is 9.449 × 10³, and the value of PSNR is 29.438. The results show that the larger the inertia parameter setting is, the better the performance of IPGM is under both convex and nonconvex constraints.
• According to the ranges of inertia parameters allowed under different constraints, the results in the above table are obtained. By comparing Table 2 and Table 3, it can be found that with the increase in inertia parameters, the PSNR of the IPGM algorithm under convex optimization constraints is higher, the objective function value is low, and the performance improves. Under the constraint of nonconvex optimization, although the PSNR of the IPGM algorithm is high, the value of the objective function is also high, and the performance is not ideal. This result shows that the simple utilized inertia term cannot achieve the perfect effect under the nonconvex constraint. In the future, we will continue to study the IPGM algorithm with dynamic inertia parameters under nonconvex optimization.

8. Conclusions

For the CDL problem, an IPGM algorithm based on a forward–backward splitting algorithm with an inertial term is proposed. The complexity of the algorithm is analyzed, and the convergence of the algorithm is proven. Finally, the IPGM algorithm is compared with other algorithms in terms of performance and effect through experiments. The results show that the IPGM algorithm produces a lower objective function value, lower sparse value, and higher reconstruction PSNR. In summary, the IPGM algorithm has many good theoretical properties, and it is efficient and straightforward.

In addition, according to the allowable range of inertia parameters under different constraints, the ablation experiment of the IPGM algorithm is carried out. According to the experimental results, it can be found that the IPGM algorithm achieves better performance under convex optimization. However, under the nonconvex optimization constraint, the performance of the IPGM algorithm does not reach the ideal value, indicating that the perfect effect cannot be achieved by simply using the inertia term under the nonconvex constraint. In the future, we will further study the IPGM algorithm with dynamic inertia parameters under nonconvex optimization.