# An Inexact Penalty Decomposition Method for Sparse Optimization

Authors: Zhengshan Dong; Geng Lin; Niandong Chen
Journal: Computational Intelligence and Neuroscience
 (2021)

---

## Abstract

The penalty decomposition method is an effective and versatile method for sparse optimization and has been successfully applied to solve compressed sensing, sparse logistic regression, sparse inverse covariance selection, low rank minimization, image restoration, and so on. With increase in the penalty parameters, a sequence of penalty subproblems required being solved by the penalty decomposition method may be time consuming. In this paper, an acceleration of the penalty decomposition method is proposed for the sparse optimization problem. For each penalty parameter, this method just finds some inexact solutions to those subproblems. Computational experiments on a number of test instances demonstrate the effectiveness and efficiency of the proposed method in accurately generating sparse and redundant representations of one-dimensional random signals.


---

## Body

## 1. Introduction

In this paper, we consider solving the following sparse optimization problem by an inexact penalty decomposition (iPD) method:(1)minx∈Xlx+λx0,s.t.gx≤0,hx=0,where λ≥0 controls the sparsity of the solution, X⊂ℝn is a closed convex set in the n-dimensional Euclidean space ℝn, l:ℝn⟶ℝ,gx:ℝn⟶ℝp are continuously differentiable convex functions, h:ℝn⟶ℝm is an affine function, and x0 denotes the number of nonzero components of x.Sparse optimization is to solve some problems whose solutions are sparse or compressed. And it has attracted considerable attention in the past ten years since its broad applications, such as signal (image) processing [1–3], linear regression [4], inverse problem [5], model selection [6], and machine learning [6, 7]. In those applications, most information of interest has or can be coded by much low dimension though its own dimension is high.However, problem (1) is NP hard even though for some simple special cases [8]. Even so, many methods have been proposed for some special cases of problem (1). These methods can be classified into four categories: (1) greedy methods: matching pursuit [9, 10] and greedy coordinate descent [11]; (2) l1-norm relaxation methods: gradient projection [12, 13], iterative shrinkage-thresholding [5, 14], iterative reweighted method [15], alternating direction method [16], and homotopy method [17–20]; (3) lp-norm 0<p<1 relaxation methods [1, 2, 21]; and (4) l0-norm based methods, e.g., penalty decomposition method [22], block decomposition method [23], iterative hard thresholding method [22, 24–29], and so on. In this paper, we mainly discuss the PD method.The PD method was proposed for solving the generall0-norm minimization problem (1) by Lu and Zhang in [22]. And it had been successfully applied to solve compressed sensing [22], sparse logistic regression [22], sparse inverse covariance selection [22], low rank minimization [30], image restoration [3] problems, and so on. Moreover, the PD method is theoretically sound. Lu et al. stated that any accumulation point of the sequence generated by the PD method satisfies the first-order optimality conditions of problem (1) when the Robinson condition holds. Hence, the PD method is an effective and versatile method for sparse optimization. However, since the PD method found exact solutions of subproblems for each penalty parameter, it may be time consuming in practice.In this paper, an inexact penalty decomposition (iPD) method is proposed for the sparse optimization problem (1). The iPD method just finds some inexact solutions to those subproblems for each penalty parameter. In more detail, for the first convex subproblem, the iPD method just takes one gradient step and then goes to solve the second nonconvex subproblem. The second subproblem can be solved by the iterative hard thresholding method [26]. After the two steps, the penalty parameter is updated. Computational experiments on a number of random instances demonstrate the effectiveness of the proposed method in accurately generating sparse and redundant representations of one-dimensional random signals.The rest of this paper is organized as follows. Section2 is the preliminary, in which some notations and the basic method are described. Section 3 presents the iPD method. Computational experiments are presented in Section 4, and conclusions are drawn in Section 5.


## 2. Preliminaries

### 2.1. Notations

In this subsection, some notations are presented to simplify presentation. The transpose of a vectorx∈ℝn is denoted by xT. If without special statement, all norms used are the Euclidean norm, denoted by ⋅2. PX⋅ denotes projection on a set X. Given a vector x∈ℝn, the nonnegative part of x is denoted by x+, i.e., x+=maxx,0. The index of nonzero components of a vector x is denoted by Sx=i:xi≠0 (called support set) and Sk:=Sxk. The size of Sx is denoted as s=Sx.Now, let us consider problem (1). It is easy to verify that problem (1) is equivalent to the following problem:(2)minx,y∈Xlx+λy0,s.t.gx≤0,hx=0,x=y.And the relative penalty function of problem (2) is defined as(3)pρx,y=lx+λy0+ρ2gx+22+hx22+x−y22,where ρ>0 is the penalty parameter.For simplicity, we also denote(4)Fρx=lx+λx0+ρ2gx+22+hx22,fx=12gx+22+hx22,qρx,y=lx+ρ2gx+22+hx22+x−y22.


### 2.2. The PD Method

In this subsection, we show the PD method proposed in [22]. First, the outline of the PD method is as presented in Algorithm 1. Then, we explain why the PD method is time consuming by a random example.Algorithm 1: The PD method [22].
Input:ρ0>0, x0, σ>1;Output:x^;(1)
initializationk⟵0, y0,0=x0;(2)
repeat(3)
i⟵0,yk,0=yk;(4)
repeat(5)
xk,i+1∈argminxpρkx,yk,i;(6)
yk,i+1∈argminypρkxk,i+1,y;(7)
i⟵i+1;(8)
untilPXxk,i−∇xqρkxk,i,yk,i2≤ϵk(9)
yk+1⟵yk,i;(10)
xk+1⟵xk,i;(11)
ρk+1⟵σρk;(12)
k⟵k+1;(13)
until some termination conditions reach(14)
x^⟵yk;Remark 1.
(i) The termination condition in Step 8 of Algorithm1 is used to establish the global convergence of the PD method. In practice, the termination criterion is based on the relative change of the sequence xk,i,yk,i such as the sequence satisfying(5)maxxk,i−xk,i−1∞maxxk,i∞,1,yk,i−yk,i−1∞maxyk,i∞,1≤ϵI,for some ϵI>0. In addition, the PD method terminates the outer iterations when(6)xk−yk∞≤ϵO,holds for some ϵO>0.
(ii) The second subproblem, i.e., in Step 6 of Algorithm1,(7)yk,i+1∈argminyλy0+ρk2y−xk,i+122,has a closed-form solution [26].(8)yjk,i+1=xjk,i+1,ifxjk,i+1>2λρk,0,otherwise,where ⋅j denotes the j-th entry of a vector, j∈1,2,…,n.In Step 5 of Algorithm1, minimizing the function pρkx,yk with respect to x is a convex problem. There exist many efficient methods for this purpose if X is simple. However, for each penalty parameter, the PD method solves the penalty subproblems a few times until some termination conditions are reached, which is time consuming.Consider a special case—compressed sensing [31]. One important task of compressed sensing is to find the sparsest solution to the underdetermined linear system, which is formulated as(9)minxx0,s.t.Ax=b,where A∈ℝm×n is the sensing matrix and b∈ℝm is the observation data. For this special problem, fx=1/2Ax−b22andFρx=x0+ρ/2Ax−b22. The value of fx is called data fidelity, and it can measure the feasibility of a solution x. Fρx is the penalty function of problem (9).Example 1.
We generate a sparse vectorx¯∗ with length n=5000 and s=100 nonzero components. These components independently follow the standard Gaussian distribution, and their locations are assigned randomly to x¯∗. Then, we create a Gaussian random matrix A with size 1000×5000, and let b=Ax¯∗. Then, we solve this instance by the PD method package, and the process data are presented as Figure 1.
Figure1 shows that the value of fx decreases slowly. It decreases steep just at the first few steps for each penalty parameter. There are many almost null steps during the process. And the value of the penalty function Fρx increases too much when updating the penalty parameter. Hence, we can just take one or a few iterations for each penalty parameter to save some time. In Section 3, we will improve the PD method by the above observations.Figure 1
Iteration process of penalty decomposition for solving compressed sensing with sizem=1000,n=5000,ands=100: (a) data fidelity at each iteration; (b) penalty function value at each iteration.
(a)(b)


## 2.1. Notations

In this subsection, some notations are presented to simplify presentation. The transpose of a vectorx∈ℝn is denoted by xT. If without special statement, all norms used are the Euclidean norm, denoted by ⋅2. PX⋅ denotes projection on a set X. Given a vector x∈ℝn, the nonnegative part of x is denoted by x+, i.e., x+=maxx,0. The index of nonzero components of a vector x is denoted by Sx=i:xi≠0 (called support set) and Sk:=Sxk. The size of Sx is denoted as s=Sx.Now, let us consider problem (1). It is easy to verify that problem (1) is equivalent to the following problem:(2)minx,y∈Xlx+λy0,s.t.gx≤0,hx=0,x=y.And the relative penalty function of problem (2) is defined as(3)pρx,y=lx+λy0+ρ2gx+22+hx22+x−y22,where ρ>0 is the penalty parameter.For simplicity, we also denote(4)Fρx=lx+λx0+ρ2gx+22+hx22,fx=12gx+22+hx22,qρx,y=lx+ρ2gx+22+hx22+x−y22.


## 2.2. The PD Method

In this subsection, we show the PD method proposed in [22]. First, the outline of the PD method is as presented in Algorithm 1. Then, we explain why the PD method is time consuming by a random example.Algorithm 1: The PD method [22].
Input:ρ0>0, x0, σ>1;Output:x^;(1)
initializationk⟵0, y0,0=x0;(2)
repeat(3)
i⟵0,yk,0=yk;(4)
repeat(5)
xk,i+1∈argminxpρkx,yk,i;(6)
yk,i+1∈argminypρkxk,i+1,y;(7)
i⟵i+1;(8)
untilPXxk,i−∇xqρkxk,i,yk,i2≤ϵk(9)
yk+1⟵yk,i;(10)
xk+1⟵xk,i;(11)
ρk+1⟵σρk;(12)
k⟵k+1;(13)
until some termination conditions reach(14)
x^⟵yk;Remark 1.
(i) The termination condition in Step 8 of Algorithm1 is used to establish the global convergence of the PD method. In practice, the termination criterion is based on the relative change of the sequence xk,i,yk,i such as the sequence satisfying(5)maxxk,i−xk,i−1∞maxxk,i∞,1,yk,i−yk,i−1∞maxyk,i∞,1≤ϵI,for some ϵI>0. In addition, the PD method terminates the outer iterations when(6)xk−yk∞≤ϵO,holds for some ϵO>0.
(ii) The second subproblem, i.e., in Step 6 of Algorithm1,(7)yk,i+1∈argminyλy0+ρk2y−xk,i+122,has a closed-form solution [26].(8)yjk,i+1=xjk,i+1,ifxjk,i+1>2λρk,0,otherwise,where ⋅j denotes the j-th entry of a vector, j∈1,2,…,n.In Step 5 of Algorithm1, minimizing the function pρkx,yk with respect to x is a convex problem. There exist many efficient methods for this purpose if X is simple. However, for each penalty parameter, the PD method solves the penalty subproblems a few times until some termination conditions are reached, which is time consuming.Consider a special case—compressed sensing [31]. One important task of compressed sensing is to find the sparsest solution to the underdetermined linear system, which is formulated as(9)minxx0,s.t.Ax=b,where A∈ℝm×n is the sensing matrix and b∈ℝm is the observation data. For this special problem, fx=1/2Ax−b22andFρx=x0+ρ/2Ax−b22. The value of fx is called data fidelity, and it can measure the feasibility of a solution x. Fρx is the penalty function of problem (9).Example 1.
We generate a sparse vectorx¯∗ with length n=5000 and s=100 nonzero components. These components independently follow the standard Gaussian distribution, and their locations are assigned randomly to x¯∗. Then, we create a Gaussian random matrix A with size 1000×5000, and let b=Ax¯∗. Then, we solve this instance by the PD method package, and the process data are presented as Figure 1.
Figure1 shows that the value of fx decreases slowly. It decreases steep just at the first few steps for each penalty parameter. There are many almost null steps during the process. And the value of the penalty function Fρx increases too much when updating the penalty parameter. Hence, we can just take one or a few iterations for each penalty parameter to save some time. In Section 3, we will improve the PD method by the above observations.Figure 1
Iteration process of penalty decomposition for solving compressed sensing with sizem=1000,n=5000,ands=100: (a) data fidelity at each iteration; (b) penalty function value at each iteration.
(a)(b)


## 3. The Proposed Method

In this section, we describe the process of the iPD method. From the outline of Algorithm1, we find that, for each penalty parameter ρk, the block coordinate descent method needs to alternately solve two minimization subproblems many times, and an example in Section 2 shows that there are many almost null step for each penalty parameter. Hence, the original PD method may be time consuming if convergence speed of the block coordinate descent is slow.Motivated by the analysis inSection 2 and the above demonstration, we accelerate the PD method by alteratively solving the two penalty subproblems once a time after updating the penalty parameter. For solving the first penalty subproblem, a gradient step is taken, and its step-length is searched by the backtracking line search method.Now, we present the outline of the accelerated penalty decomposition method as follows.Remark 2.
A practical termination criterion in Step 11 of Algorithm2 can be(10)xk−yk∞maxxk∞,1≤tol,for some tol>0.Algorithm 2: The inexact PD method.
Input:ρ0>0, x0, σ0>1, L0>0, γinc>1, γdec>1;Output:x^;(1)
initializationk⟵0, yk=xk;(2)
repeat(3)
whilepρkxk+1,yk+η/2xk+1−yk22>pρkyk,ykdo(4)
Lk=minγincLk,Lmax(5)
xk+1=yk−∇xpρkyk,yk/Lk;(6)
end while(7)
yk+1∈argminypρkxk+1,y;(8)
ρk+1=σkρk(9)
Lk+1=minLk/γdec,Lmin(10)
k⟵k+1;(11)
until some termination conditions reach(12)
x^⟵yk;Theorem 1.
If the gradient of the functionpρkx,yk with respect to x is Lipschitz continuous (its Lipschitz constant is denoted as Lp), then the line search between Steps 3–6 can be terminated in a finite number of iterations.Proof.
Sincepρkxk+1,yk satisfies(11)pρkxk+1,yk≤pρkyk,yk+∇xpρkyk,ykTxk+1−yk+Lp2xk+1−yk22,it together with xk+1=yk−∇xpρkyk,yk/Lk implies that(12)pρkxk+1,yk≤pρkx,yk−Lk−Lp2xk+1−xk22.
Then, ifLk≥Lp+η/2,(13)pρkxk+1,yk+η2xk+1−yk22≤pρkyk,yk,holds, which means that the while loop in Algorithm 2 terminates if Lk≥Lp+η/2. Let L¯k be the final value of Lk after the while loop. Then, L¯k/γinc≤Lp+η/2 holds, i.e., L¯k≤γincLp+η/2. Let n^k be the number of iterations in the while loop at the k-th iteration. Then, one can get that(14)Lminγincn^k−1≤Lk0γincn^k−1≤L¯k≤γincLp+η2,where Lk0 is the initial value of Lk in the line search. Therefore,(15)n^k≤N≔logLp+η−log2Lminlogγ+2.


## 4. Experiments

In this section, we implement the proposed accelerated PD method to solve the compressed sensing problem. To verify the efficiency of PD empirically, a large number of computational experiments are performed on one-dimensional random signals. We mainly compare the performance of our improved PD method with that of the original PD method [22]. All experiments were performed on a personal computer with an Intel(R) Core(TM)i7-7700HQ CPU (2.80 GHz) and 8 GB memory, using a MATLAB toolbox (version R2018b).We compare the performance of the compared methods by the CPU time (in seconds) required, the size of the support set of the reconstructed datax^, and the mean squared error (MSE) with respect to x¯∗, which is defined as(16)MSE=1nx^−x¯∗,and the data fidelity of Ax^−y is defined as(17)DF=12Ax^−y2,and NS as the number of successfully recovered instances. We say a signal x^ is successfully recovered if the positions of the nonzero components of x^ are the same as x¯∗ and the corresponding MSE value is less than 10−4.


### 4.1. Data Generation and Parameter Setting

Each instance is generated randomly with sizem,n,s, where m×n is the dimension of matrix A and s is the sparsity level, such as m=1000, n=5000, and s=100. The elements of matrix A follow the Gaussian distribution. The vector x¯∗ is generated with the same distribution at s randomly chosen coordinates. Finally, the vector b is generated by b=Ax¯∗.Unless otherwise stated, all parameters in the PD method are set as default, and parameters in the IPD package are set as in Table1.Table 1
Parameter settings in the acceleration of the PD method.
ParameterValuex00ρρ0=10,ρk+1=min1.1ρk,1015tol10−6η1LL0=0.1maxAj22γinc=2,γdec=3


### 4.2. Compare with the Original PD Method

Firstly, we compare the iteration process of the iPD method with that of the PD method on a random instance. All parameters are set as before, and the problem size ism=1000,n=5000,ands=100. Figure 2 describes the data fidelity and the penalty function value over the iteration process. From Figure 2(a), we find that the iPD method does not have many null steps, and the values of data fidelity generated by the iPD method decrease much fast than those of the original PD method. Furthermore, the iPD method just requires about 150 steps while the original PD method requires about 400 steps. And the running time of the iPD method is about 7 seconds, which is less than half of the time required by the original PD method. Moreover, the penalty function value generated by the iPD method is much stable than that by the original PD method.Figure 2
Iteration process of the compared methods for solving compressed sensing with sizem=1000,n=5000,ands=100: (a) data fidelity at each iteration; (b) penalty function value at each iteration.
(a)(b)In the second experiment, we compare the accelerated PD method with its original PD method at different sampling numbers. We fix the dimensionm=5000 and the sparsity level s=100. For each sampling number m, 100 instances are generated, and the averaged performance of the two methods is presented in Figure 3.Figure 3
Averaged results of the penalty decomposition methods for the compressed sensing problem with different sampling numbers on 100 instances: (a) CPU time over sampling number; (b) recovered rate over sampling number; (c) MSE over sampling number; (d) data fidelity over sampling number; (e) number of nonzero components over sampling number.
(a)(b)(c)(d)(e)From Figure3(a), we see that the accelerated PD method requires not more than 10 seconds while the original PD method requires much more time. And the time required by the accelerated PD method is stable at different sampling numbers. Figure 3(b) shows that the recovered rate by the accelerated PD method is higher than that by the original PD method when m is bigger than 600. When the sampling number is bigger than 700, the accelerated PD method can recover all signals successfully. We find that the MSE value and the DF value generated by the accelerated PD method are lower than those generated by the original PD method. The averaged number of nonzero components also shows that the accelerated PD method performs better.In the next experiment, we compare the accelerated PD method with its original version for solving the compressed sensing problem with different sparsity levelss. All parameters are set as the same value as those stated before. The averaged computational results on 100 instances are presented in Table 2.Table 2
Averaged results on 100 instances with sizem=1000andn=5000 for each sparsity level s
AlgorithmsTime (s)NNZMSEDFNSPD5018.2505.62×10−95.95×10−710010021.6999.951.06×10−82.61×10−69515026.36149.882.50×10−83.50×10−58920076.38297.656.72×10−58.19×10−273250102.541225.211.22×10−33.53×10−10270105.181283.921.58×10−31.93×10−10iPD506.72503.22×10−91.77×10−71001008.5699.995.20×10−93.62×10−7991509.781507.31×10−95.85×10−710020010.322009.32×10−97.66×10−710025011.75656.875.50×10−51.03×10−63027012.751235.93.67×10−49.72×10−71From Table2, we find that the PD method not works well when the sparsity level is greater than 150, especially when it is greater than 200. However, the sparsity level recovered by the iPD method can reach 200. When the two methods can recover sparse signals, the iPD method just needs about one third of the time required by the PD method. Moreover, the recovered rate of the iPD method is higher than that of the original PD method. From MSE and DF value, we see that the signals recovered by the iPD method are more exact than those recovered by the PD method. When s=100, there is one instance not recovered exactly by the iPD method since there exist several very small components and one of them is not recovered.


## 4.1. Data Generation and Parameter Setting

Each instance is generated randomly with sizem,n,s, where m×n is the dimension of matrix A and s is the sparsity level, such as m=1000, n=5000, and s=100. The elements of matrix A follow the Gaussian distribution. The vector x¯∗ is generated with the same distribution at s randomly chosen coordinates. Finally, the vector b is generated by b=Ax¯∗.Unless otherwise stated, all parameters in the PD method are set as default, and parameters in the IPD package are set as in Table1.Table 1
Parameter settings in the acceleration of the PD method.
ParameterValuex00ρρ0=10,ρk+1=min1.1ρk,1015tol10−6η1LL0=0.1maxAj22γinc=2,γdec=3


## 4.2. Compare with the Original PD Method

Firstly, we compare the iteration process of the iPD method with that of the PD method on a random instance. All parameters are set as before, and the problem size ism=1000,n=5000,ands=100. Figure 2 describes the data fidelity and the penalty function value over the iteration process. From Figure 2(a), we find that the iPD method does not have many null steps, and the values of data fidelity generated by the iPD method decrease much fast than those of the original PD method. Furthermore, the iPD method just requires about 150 steps while the original PD method requires about 400 steps. And the running time of the iPD method is about 7 seconds, which is less than half of the time required by the original PD method. Moreover, the penalty function value generated by the iPD method is much stable than that by the original PD method.Figure 2
Iteration process of the compared methods for solving compressed sensing with sizem=1000,n=5000,ands=100: (a) data fidelity at each iteration; (b) penalty function value at each iteration.
(a)(b)In the second experiment, we compare the accelerated PD method with its original PD method at different sampling numbers. We fix the dimensionm=5000 and the sparsity level s=100. For each sampling number m, 100 instances are generated, and the averaged performance of the two methods is presented in Figure 3.Figure 3
Averaged results of the penalty decomposition methods for the compressed sensing problem with different sampling numbers on 100 instances: (a) CPU time over sampling number; (b) recovered rate over sampling number; (c) MSE over sampling number; (d) data fidelity over sampling number; (e) number of nonzero components over sampling number.
(a)(b)(c)(d)(e)From Figure3(a), we see that the accelerated PD method requires not more than 10 seconds while the original PD method requires much more time. And the time required by the accelerated PD method is stable at different sampling numbers. Figure 3(b) shows that the recovered rate by the accelerated PD method is higher than that by the original PD method when m is bigger than 600. When the sampling number is bigger than 700, the accelerated PD method can recover all signals successfully. We find that the MSE value and the DF value generated by the accelerated PD method are lower than those generated by the original PD method. The averaged number of nonzero components also shows that the accelerated PD method performs better.In the next experiment, we compare the accelerated PD method with its original version for solving the compressed sensing problem with different sparsity levelss. All parameters are set as the same value as those stated before. The averaged computational results on 100 instances are presented in Table 2.Table 2
Averaged results on 100 instances with sizem=1000andn=5000 for each sparsity level s
AlgorithmsTime (s)NNZMSEDFNSPD5018.2505.62×10−95.95×10−710010021.6999.951.06×10−82.61×10−69515026.36149.882.50×10−83.50×10−58920076.38297.656.72×10−58.19×10−273250102.541225.211.22×10−33.53×10−10270105.181283.921.58×10−31.93×10−10iPD506.72503.22×10−91.77×10−71001008.5699.995.20×10−93.62×10−7991509.781507.31×10−95.85×10−710020010.322009.32×10−97.66×10−710025011.75656.875.50×10−51.03×10−63027012.751235.93.67×10−49.72×10−71From Table2, we find that the PD method not works well when the sparsity level is greater than 150, especially when it is greater than 200. However, the sparsity level recovered by the iPD method can reach 200. When the two methods can recover sparse signals, the iPD method just needs about one third of the time required by the PD method. Moreover, the recovered rate of the iPD method is higher than that of the original PD method. From MSE and DF value, we see that the signals recovered by the iPD method are more exact than those recovered by the PD method. When s=100, there is one instance not recovered exactly by the iPD method since there exist several very small components and one of them is not recovered.


## 5. Conclusions

In this paper, we have proposed an acceleration of the penalty decomposition for the sparse approximation problem. The proposed method does not solve the penalty subproblems exactly and alternately solve penalty subproblems once a time after updating penalty parameters. We show that this method enhances the performance of the penalty decomposition method by computational experiments on a number of random instances for solving the compressed sensing problem. The experiments demonstrate that the proposed method indeed improves the original PD method since it recovers better solutions with less running time.



---

*Source: 9943519-2021-07-15.xml*