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.

1. Introduction

In this paper, we consider solving the following sparse optimization problem by an inexact penalty decomposition (iPD) method:where controls the sparsity of the solution, is a closed convex set in the -dimensional Euclidean space , are continuously differentiable convex functions, is an affine function, and denotes the number of nonzero components of .

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 [13], 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) -norm relaxation methods: gradient projection [12, 13], iterative shrinkage-thresholding [5, 14], iterative reweighted method [15], alternating direction method [16], and homotopy method [1720]; (3) -norm relaxation methods [1, 2, 21]; and (4) -norm based methods, e.g., penalty decomposition method [22], block decomposition method [23], iterative hard thresholding method [22, 2429], and so on. In this paper, we mainly discuss the PD method.

The PD method was proposed for solving the general -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. Section 2 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 vector is denoted by . If without special statement, all norms used are the Euclidean norm, denoted by . denotes projection on a set . Given a vector , the nonnegative part of is denoted by , i.e., . The index of nonzero components of a vector is denoted by (called support set) and . The size of is denoted as .

Now, let us consider problem (1). It is easy to verify that problem (1) is equivalent to the following problem:

And the relative penalty function of problem (2) is defined aswhere is the penalty parameter.

For simplicity, we also denote

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.

Input:, , ;
Output:;
(1)initialization , ;
(2)repeat
(3);
(4)repeat
(5)  ;
(6)  ;
(7)  ;
(8)until
(9);
(10);
(11);
(12);
(13)until some termination conditions reach
(14);

Remark 1. (i) The termination condition in Step 8 of Algorithm 1 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 such as the sequence satisfyingfor some . In addition, the PD method terminates the outer iterations whenholds for some .
(ii) The second subproblem, i.e., in Step 6 of Algorithm 1,has a closed-form solution [26].where denotes the -th entry of a vector, .

In Step 5 of Algorithm 1, minimizing the function with respect to is a convex problem. There exist many efficient methods for this purpose if 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 aswhere is the sensing matrix and is the observation data. For this special problem, . The value of is called data fidelity, and it can measure the feasibility of a solution . is the penalty function of problem (9).

Example 1. We generate a sparse vector with length and nonzero components. These components independently follow the standard Gaussian distribution, and their locations are assigned randomly to . Then, we create a Gaussian random matrix with size , and let . Then, we solve this instance by the PD method package, and the process data are presented as Figure 1.
Figure 1 shows that the value of 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 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.

3. The Proposed Method

In this section, we describe the process of the iPD method. From the outline of Algorithm 1, we find that, for each penalty parameter , 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 in Section 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 Algorithm 2 can befor some .

Input:, , , , , ;
Output:;
(1)initialization , ;
(2)repeat
(3)whiledo
(4)  
(5)  ;
(6)end while
(7);
(8)
(9)
(10);
(11)until some termination conditions reach
(12);

Theorem 1. If the gradient of the function with respect to is Lipschitz continuous (its Lipschitz constant is denoted as ), then the line search between Steps 3–6 can be terminated in a finite number of iterations.

Proof. Since satisfiesit together with implies thatThen, if ,holds, which means that the while loop in Algorithm 2 terminates if . Let be the final value of after the while loop. Then, holds, i.e., . Let be the number of iterations in the while loop at the -th iteration. Then, one can get thatwhere is the initial value of in the line search. Therefore,

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 data , and the mean squared error (MSE) with respect to , which is defined asand the data fidelity of is defined asand NS as the number of successfully recovered instances. We say a signal is successfully recovered if the positions of the nonzero components of are the same as and the corresponding MSE value is less than .

4.1. Data Generation and Parameter Setting

Each instance is generated randomly with size , where is the dimension of matrix and is the sparsity level, such as , , and . The elements of matrix follow the Gaussian distribution. The vector is generated with the same distribution at randomly chosen coordinates. Finally, the vector is generated by .

Unless otherwise stated, all parameters in the PD method are set as default, and parameters in the IPD package are set as in Table 1.

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 is . 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.

In the second experiment, we compare the accelerated PD method with its original PD method at different sampling numbers. We fix the dimension and the sparsity level . For each sampling number , 100 instances are generated, and the averaged performance of the two methods is presented in Figure 3.

From Figure 3(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 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 levels . All parameters are set as the same value as those stated before. The averaged computational results on 100 instances are presented in Table 2.

From Table 2, 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 , 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.

Data Availability

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This research was funded partly by the Natural Science Foundation of Fujian Province of China, under grant 2020J01843, and the Science and Technology Project of the Education Bureau of Fujian, China, under grant JAT200403.