Difference of two norms-regularizations for Q-Lasso

Abstract

The focus of this paper is in Q-Lasso introduced in Alghamdi et al. (2013) which extended the Lasso by Tibshirani (1996). The closed convex subset Q belonging in a Euclidean m-space, for m∈IN, is the set of errors when linear measurements are taken to recover a signal/image via the Lasso. Based on a recent work by Wang (2013), we are interested in two new penalty methods for Q-Lasso relying on two types of difference of convex functions (DC for short) programming where the DC objective functions are the difference of l1 and lσq norms and the difference of l1 and lr norms with r>1. By means of a generalized q-term shrinkage operator upon the special structure of lσq norm, we design a proximal gradient algorithm for handling the DC l1−lσq model. Then, based on the majorization scheme, we develop a majorized penalty algorithm for the DC l1−lr model. The convergence results of our new algorithms are presented as well. We would like to emphasize that extensive simulation results in the case Q={b} show that these two new algorithms offer improved signal recovery performance and require reduced computational effort relative to state-of-the-art l1 and lp (p∈(0,1)) models, see Wang (2013). We also devise two DC Algorithms on the spirit of a paper where exact DC representation of the cardinality constraint is investigated and which also used the largest-q norm of lσq and presented numerical results that show the efficiency of our DC Algorithm in comparison with other methods using other penalty terms in the context of quadratic programing, see Jun-ya et al. (2017).

Citation

Moudafi, A. (2021), "Difference of two norms-regularizations for Q-Lasso", Applied Computing and Informatics, Vol. 17 No. 1, pp. 79-89. https://doi.org/10.1016/j.aci.2018.07.002

1. Introduction and preliminaries

The process of compressive sensing (CS) [8], which consists of encoding and decoding, is rapidly consolidated year after year due to the blooming of large datasets which become increasingly important and available. The process of encoding involves taking a set of (linear) measurements, b=Ax, where A is a matrix of size m×n. If m<n, we can compress the signal x∈IRn, whereas the process of decoding is to recover x from b where x is assumed to be sparse. It can be formulated as an optimization problem, namely

(1.1)min∥x∥0 subject to Ax=b,

where ∥⋅∥0 is the l0 norm, which counts the number of nonzero entries of x; namely

(1.2)∥x∥0=|{xi;xi≠0}|

with |·| being here the cardinality, i.e., the number of elements of a set. Hence minimizing the l0 norm amounts to finding the sparsest solution. One of the difficulties in CS is solving the decoding problem above, since l0 optimization is NP-hard. An approach that has gained popularity is to replace l0 by the convex norm l1 since it often gives a satisfactory sparse solution and has been applied in many different fields such as geology and ultrasound imaging.

More recently, nonconvex metrics were used as alternative approaches to l1, especially the nonconvex metric lp for p∈(0,1) in [6] which can be interpreted as a continued approximation strategy of l0 as p→0. A great deal of research has been conducted into lp problems including all kinds of variants and related algorithms, as you can see in [4] and references therein. The convex l1 relaxation compared to the nonconvex problem (lp) is generally more difficult to handle. However, it was shown in [12] that the potential reduction method can solve this special nonconvex problem in polynomial time with arbitrarily given accuracy.

Most recently, the majority of such sparsity inducing functions are unified as the notion of DC programming in [9], including log-sum, smoothly clipped absolute deviation and capped-l1 penalty. Generally, DC programming problem can be solved through a primal–dual convex relaxations algorithm which is famous in the literature of DC Programming [11]. Other algorithms appeared as for solving application problems of DC programming in the area of finance and insurance, data analysis, machine learning as well as signal processing. However, as noted in [18], among the above mentioned DC programming approaches for sparse reconstruction, most of them are mainly preserving the separability properties of both l0 and l1 norms.

To begin with, let us recall that the lasso of Tibshirani [16] is given by the following minimization problem

(1.3)minx∈IRn12∥Ax−b∥22+γ∥x∥1,

A being an m×n real matrix, b∈IRm and γ>0 is a tuning parameter. The latter is nothing else than the basic pursuit (BP) of Chen et al. [7], namely

(1.4)minx∈IRn∥x∥1 such that Ax=b.

However, the constraint Ax=b being inexact due to errors of measurements, the problem (1.4) can be reformulated as

(1.5)minx∈IRn∥x∥1 subject to∥Ax−b∥p≤ ε,

where ε>0 is the tolerance level of errors and p is often 1,2 or ∞. It is noticed in [1] that (1.5) can be rewritten as

(1.6)minx∈IRn∥x∥1 subject to Ax∈Q,

in the case when Q:=Bε(b), the closed ball in IRn with center b and radius ε.

Now, when Q is a nonempty closed convex set of IRm and PQ the orthogonal projection from IRm onto the set Q and by observing that the constraint is equivalent to the condition Ax−PQ(Ax)=0, this leads to the following Lagrangian formulation

(1.7)minx∈IRn12∥(I−PQ)Ax∥22+γ∥x∥1,

γ>0 being a Lagrangian multiplier.

A link is also made in [1] with split feasibility problems [5] which consist in finding x satisfying

(1.8)x∈C, Ax∈Q,

with C and Q two nonempty closed convex subsets of IRn and IRm, respectively. An equivalent formulation of (1.8) as a minimization problem is given by

(1.9)minx∈C12∥(I−PQ)Ax∥22,

and its l1-regularization is

(1.10)minx∈C12∥(I−PQ)Ax∥22+γ∥x∥1,

with γ>0 a regularization parameter.

This convex relaxation approach was frequently employed, see for example [1,20] and references there in. As the level curves of l1-l2 are closer to l0 than those of l1, this motivated us in [14] to propose a regularization of split feasibility problems by means of the nonconvex l1-l2, namely

(1.11)minx∈C12∥(I−PQ)Ax∥22+γ(∥x∥1−∥x∥2),

and present three algorithms with their convergence properties [14]. Unlike the separable sparsity inducing functions involved in the aforementioned DC programming for problem (l0), we are interested in the two first sections of this work to two specific types of DC programming with un-separable objective functions, which are in the form of difference functions between two norms, namely the new notion lσq denoting the sum of q largest elements of a vector in magnitude (i.e., the l1 norm of q-term best approximation of a vector) introduced in [18] and the classical lr norm with r>1. Obviously lσq and lr (r>1) are regular convex norms. The corresponding DC programs are as follows:

(1.12)minx∈IRn(∥x∥1−ε∥x∥σq:Ax=b),

and

(1.13)minx∈IRn(∥x∥1−ε∥x∥r:Ax=b),

where ε∈(0,1],∥x∥σq is defined as the sum of the q largest elements of x in magnitude, q∈{1,2,···n} and r>1. We would like to emphasize that the following least-squares variant of (1.12) and (1.13), were studied in the recent work by Wang [18]:

(1.14)minx(f(x):=12∥Ax−b∥2+μ(∥x∥1−ε∥x∥σq)),

where μ>0 and ε∈(0,1), and

(1.15)minx(f˜(x):=12∥Ax−b∥2+μ(∥x∥1−ε∥x∥r)),

where r>0 and ε∈(0, 1).

This paper proposes generalizations to Q-Lasso, namely

minx(f(x):=12∥(I−PQ)Ax∥2+μ(∥x∥1−ε∥x∥σq)),

where μ>0 and ε∈(0, 1), as well as

minx(f˜(x):=12∥(I−PQ)Ax∥2+μ(∥x∥1−ε∥x∥r)),

where r>0 and ε∈(0,1), and our attention will be focused on the algorithmic aspect.

The rest of the paper is organized as follows. In Sections 2 and 3, two DC-penalty methods instead of conventional methods such as l1 or l1−l2 minimization are proposed. Their convergence to a stationary point are also analyzed. The first iterative minimization method is based on the gradient proximal algorithm and the second one is designed by means of the majored penalty strategy. Furthermore, relying on DCA (difference of convex algorithm) two other algorithms are proposed and their convergence results are established in Section 3 and 4.

2. Proximal gradient algorithm

First, we recall that the subdifferential of a convex function φ is given by

(2.1)∂φ(x):={u∈IRn;φ(y)≥φ(x)+〈u,y−x〉∀y∈IRn}.

Each element of ∂φ(x) is called subgradient. If φ(x)=12∥(I−PQ)Ax∥2, it is well-known that

(2.2)∂φ(x)=∇φ(x)=AT(I−PQ)Ax,

and when φ(x)=∥x∥1, we have

(2.3)(∂φ(x))i={sgn(xi)if   xi≠0;[−1,1]if   xi=0.

The indicator function of a set C ⊆ IRn is defined by

(2.4)iC(x)={0if x∈C;+∞otherwise.

Moreover, the normal cone of a set C at x∈C, denoted by NC(x) is defined as

(2.5)NC(x):={d∈IRn|〈d,y−x〉≤0,∀y∈C}.

Connection between the above definitions is given by the key relation ∂iC=NC.

In this section our interest is in solving the DC programming

(2.6)minx(f(x):=12∥(I−PQ)Ax∥2+μ(∥x∥1−ε∥x∥σq)),

where μ>0 and ε∈(0,1).

Similar to l1 norm, l2 norm, etc., we adopt the notation ∥x∥σq to denote the norm of lσq which is defined a line below (1.13) and we design an iterative algorithm based both on a generalized q-term shrinkage operator and on the proximal gradient algorithm framework.

At this stage, observe that the restriction on ε guarantees that f(x)≥0 for all x. To solve (2.6), we consider the following standard proximal gradient algorithm:

• 1. Initialization: Let x0 be given and set L>λmax(ATA) with λmax(ATA) the maximal eigenvalue.

• 2. For k=0,1,··· find

(2.7)xk+1∈Argminx∈IRn(〈AT(I−PQ)Axk,x−xk〉+L2∥x−xk∥2+μ(∥x∥1−ε∥x∥σq)).

Observe that subproblem (2.7) can equivalently formulated as

(2.8)minx(L2∥x−(xk−1LAT(I−PQ)Axk)∥2+μ(∥x∥1−ε∥x∥σq)).

Thus, it suffices to consider the solutions to the following minimization problem

(2.9)minx(12∥x−y∥2+λ1∥x∥1−λ2∥x∥σq),

with a given vector y and positive numbers λ1>λ2>0. An explicit solution of this problem is given by the following result, see [18].

Proposition 2.1 Let {i1,⋯,in} be the indices such that

|yi1|≥|yi2|≥⋯≥|yin|.

Then x∗:=proxλ1∥x∥1−λ2∥x∥σq(y) with

(2.10)xi∗={sign(yi)max{|yi|−(λ1−λ2),0}if  i=i1,i2,···,iq;sign(yi)max{|yi|−λ1,0}otherwise

is a solution of (2.9).

The proximal operator above (called the generalized q-term shrinkage operator in [18]) amounts to write the algorithm as follows:

Proximal Gradient Algorithm:

1. 1.

   Start: Let x0 be given and set L>λmax(ATA) with λmax(ATA) the maximal eigenvalue.

2. 2.

   For k=0,1,⋯ find

yk+1=(xk−1LAT(I−PQ)Axk),

Sort yk+1 as |yi1|≥|yi2|≥⋯≥|yin|,

(2.11)(xk+1)i={sign(yi)max{|yi|−μ(1−ε),0}if i=il;sign(yi)max{|yi|−μ,0}otherwise

where l=1,⋯,q.

End.

Now, we are in a position to show the following convergence result of the scheme (2.7):

Proposition 2.2 The sequence (xk) generated by the Proximal Gradient Algorithm above converges to a stationary point of problem (2.6).

Proof. Remember that h(x)=12∥(I−PQ)Ax∥2 is differentiable and its gradient ∇h(x)=AT(I−PQ)Ax is Lipschitz continuous with constant L˜:=λmax(ATA). By [3]-Proposition A.24, we have

f(xk+1)≤12∥(I−PQ)Axk∥2+〈AT(I−PQ)Axk,xk+1−xk〉+L˜2∥xk+1−xk∥2+μ(∥xk+1∥1−ε∥xk+1∥σq).

Combining this with definition of xk+1, we obtain

(2.12)f(xk+1)≤f(xk)−L−L˜2∥xk+1−xk∥2.

Since L>L˜, we see immediately that f(xk+1)≤f(xk) and thus the sequence (f(xk)) is convergent since f is a non-negative function. Furthermore, we obtain that ∑k∥xk+1−xk∥2<+∞ which follows by summing (2.12) from k=0 to ∞. As a further consequence, we note that

μ(1−ε)∥xk∥1≤μ(∥xk∥1−ε∥xk∥σq)≤f(xk)≤f(x0).

Since μ(1−ε)>0, we have that (xk) is bounded. Moreover, the objective function f is square term plus a piecewise linear function which ensures that f is semi-algebric and hence satisfies Kurdyka-Lojasiewicz inequality. [2]-Theorem 5.1 is then applicable and obtain that (xk) is convergent to a stationary point of (2.6). □

3. Majorized penalty algorithm

Consider the following minimization problem

(3.1)minx(f˜(x):=12∥(I−PQ)Ax∥2+μ(∥x∥1−ε∥x∥r),

where A∈IRm×n,Q a nonempty closed convex set of IRm,r>0 and ε∈(0,1).

First, observe again that conditions on ε guarantees that f˜(x)≥0 for all x. We will now describe an algorithm for solving (3.1), based on the majorized penalty approach see, for example, [18] and references therein. Following the same lines as in [18], we start by constructing a majorization of f˜. To that end let L>λmax(ATA), then for any x,y∈IRn, we have

12∥(I−PQ)Ax∥22≤12∥(I−PQ)Ay∥22+〈AT(I−PQ)Ay,x−y〉+L2∥x−y∥22.

Moreover, by invoking the convexity of the norm ∥x∥r and definition of its subdifferential, we also have

∥x∥r≥∥y∥r+〈g(y),x−y〉with g(y)∈∂∥y∥r,

where

(3.2)[g(y)]i={sign(yi)|yi|r−1∥y∥rr−1if y≠0;0otherwise.

Hence, if we define

F(x,y)=12∥(I−PQ)Ay∥22+〈AT(I−PQ)Ay,x−y〉+L2∥x−y∥2+μ(∥x∥1−ε∥y∥r−ε〈g(y),x−y〉),

hence, for every x,y∈IRn, we get

F(x,y)≥f˜(x)and F(y,y)=f˜(y).

Starting with an initial iterate x0, the majorized penalty approach above updates xk by solving

(3.3)xk+1=argminxF(x,xk).

This leads to the following explicit formulation of xk+1 by means of the proximity (shrinkage) operator of ∥x∥1:

xk+1=argminx(〈AT(I−PQ)Axk,x−xk〉+L2∥x−xk∥2+μ(∥x∥1−ε〈g(xk),x−xk〉))=argminx(L2∥x−xk+1L(AT(I−PQ)Axk−μεg(xk))∥2+μ∥x∥1)=proxμL∥x∥1(xk−1L(AT(I−pQ)Axk−μεg(xk)))=sgn(vk)∘max{|vk|−μL,0},

where

vk=xk−1L(AT(I−PQ)Axk−μεg(xk))with g(xk)∈∂∥xk∥r.

We summarize the algorithm as follows:

Majorized Penalty Algorithm:

• 1. Initialization: Let x0 be given and set L>λmax(ATA).

• 2. For k=0,1,⋯ find

(3.4)xk+1=sgn(xk−1L(AT(I−PQ)Axk−μεg(xk)))∘max{|vk|−μL,0}

End.

The following proposition contains the convergence result of this Penalty Algorithm.

Proposition 3.1 Let (xk) be the sequence generated by the Majorized Penalty Algorithm above. Then

(3.5)L2∥xk−xk+1∥2≤f˜(xk)−f˜(xk+1).

Furthermore, the sequence (xk) is bounded and any cluster point is a stationary point of problem (3.1).

Proof. Since xk+1 minimizes F(x,xk), thanks to the first-order optimality condition we can write

(3.6)0∈AT(I−PQ)Axk+L(xk+1−xk)+μ∥xk+1∥1−μεg(xk),

g(xk) being a subgradient of ∥x∥r at xk+1. This combined with the definition of the subdifferential of ∥x∥1 at xk+1 gives

μ∥xk∥1−μ∥xk+1∥1≥〈−AT(I−PQ)Axk−L(xk+1−xk)+μεg(xk),xk−xk+1〉=〈−AT(I−PQ)Axk+μεg(xk),xk−xk+1〉+L∥xk+1−xk∥22=〈AT(I−PQ)Axk+μεg(xk),xk+1−xk〉+L∥xk+1−xk∥22.

Hence

μ∥xk+1∥1−μ∥xk∥1+〈AT(I−PQ)Axk−μεg(xk),xk−xk+1〉≤−L∥xk+1−xk∥22.

This together with the definition of F, for any k≥1, leads to

f˜(xk+1)−f˜(xk)≤F(xk+1,xk)−f˜(xk)=〈AT(I−PQ)Axk,xk+1−xk〉+L2∥xk+1−xk∥22+ μ(∥xk+1∥1−∥xk∥1−εg(xk),xk−xk+1〉)=L2∥xk+1−xk∥22+μ∥xk+1∥1−μ∥xk∥1+〈AT(I−PQ)Axk−μεg(xk),xk+1−xk〉≤L2∥xk+1−xk∥22−L∥xk+1−xk∥22.

Consequently,

(3.7)f˜(xk+1)−f˜(xk)≤−L2∥xk+1−xk∥22.

Hence f˜(xk+1)≤f˜(xk) and thus the sequence (f˜(xk)) is convergent since f˜ is a non-negative function. Furthermore, the sequence (xk) is such that

∑k=0∞∥xk+1−xk∥2<+∞.

Indeed, by summing (3.7) from k=0 to ∞, we obtain that

L2∑k=0∞∥xk+1−xk∥22≤f˜(x0)−limk→+∞f˜(xk)≤f˜(x0)<+∞.

Consequently, the sequence (xk) is asymptotically regular, i.e., limk→+∞∥xk−xk+1∥=0. On the other hand, observe that the definition of f˜ for any k≥1, leads to

μ(∥xk∥1−ε∥xk∥r)≤12∥(I−PQ)Axk∥22+μ(∥xk∥1−ε∥xk∥r)=f˜(xk)≤ f˜(x0).

Since ∥xk∥1≥∥xk∥r, we obtain that μ(1−ε)∥xk∥r≤f˜(x0). This implies that (xk) is bounded since 0<ε<1. To conclude, we prove that every cluster point of (xk) is a stationary point of (3.1). Let x∗ be a cluster point of (xk), then x∗=limVxkv,(xkv) being subsequence of (xk). By passing to the limit in (3.6) along the subsequence (xkv) and in the light of the upper semicontinuity of (Clarke) subdifferentials, we obtain the desired result, namely

0∈AT(I−PQ)Ax∗+μ∂∥x∗∥1−μεg(x∗),

which is nothing else than the first-order optimality condition of (3.1). □

4. DCA algorithm

Now we turn our attention to a DC Algorithm (DCA), where the dual step at each iteration can be efficiently carried out due to the accessible subgradients of the largest- q-norm ∥·∥σq and ∥·∥r norm. Remember that to find critical points of f:=φ−ψ, the DCA consists in designing of sequences (xk) and (yk) by the following rules

(4.1){yk∈∂ψ(xk);xk+1=argminx∈IRn(φ(x)−(ψ(xk)+〈yk,x−xk〉)).

Note that by the definition of subdifferential, we can write

ψ(xk+1)≥ψ(xk)+〈yk,xk+1−xk〉.

Since xk+1 minimizes φ(x)−(ψ(xk)+〈yk,x−xk〉), we also have

φ(xk+1)−(ψ(xk)+〈yk,xk+1−xk〉)≤φ(xk)−ψ(xk).

Combining the last inequalities, we obtain

f(xk)=φ(xk)−ψ(xk)≥φ(xk+1)−(ψ(xk)+〈yk,xk+1−xk〉)≥f(xk+1).

Therefore, the DCA leads to a monotonically decreasing sequence (f(xk)) that converges as long as the objective function f is bounded below.

Now, we can decompose the objective function in (2.6) as follows

(4.2)minx(f(x):=(12∥(I−PQ)Ax∥2+μ∥x∥1)−(με∥x∥σq)),

where μ>0, ε∈(0,1), here φ(x)=12∥(I−PQ)Ax∥2+μ∥x∥1 and ψ(x)=με∥x∥σq.

At each iteration, DCA solves The convex subproblem defined by linearizing the concave term −ε∥x∥σq is solved by DCA at each iteration until a convergence condition is satisfied. More precisely, we have

(4.3){yk∈με ∂∥xk∥σq;xk+1=argminx∈IRn(12∥(I−PQ)Ax∥2+μ∥x∥1−(με∥xk∥σq+〈yk,x−xk〉)).

Especially, if either the function φ or ψ is polyhedral, the DCA is said to be polyhedral and terminates in finite iterations [15]. Note that the our proposed DCA is polyhedral since the largest-q norm term −ε∥x∥σq can be expressed as a pointwise maximum of 2qCnq linear functions, see [10]. On the other hand, the subdifferential of ∥x∥σq at a point xk is given in, see for example [19],

(4.4)∂∥xk∥σqargmaxy(∑i=1n|[xk]i|yi:∑i=1nyi=q,0≤yi≤1,i=1,⋯,n),

that is

∂∥xk∥σq={(y1,⋯,yn):yi1=⋯=yiq=1,yiq+1=0=⋯=yin=0},

where yij denotes the element of y corresponding to xij in the linear program (4.4). Observe that a subgradient y∈∂∥xk∥σq can be computed efficiently by first sorting the elements |xi| in decreasing order, namely |xi1|≥|xi2|≥⋯≥|xiq|. Then, assign 1 to yi which corresponds to xi1⋯,xiq.

To conclude, let us consider the following DC formulation of (3.1):

(4.5)minx(f(x):=(12∥(I−PQ)Ax∥2+μ∥x∥1)−(με∥x∥r)),

where r>0, ε∈(0,1), here φ(x)=12∥(I−PQ)Ax∥2+μ∥x∥1 and ψ(x)=με∥x∥r.

The subgradient y∈∂∥xk∥r is also available via the formula (3.2) and the DCA in this context take the following form

(4.6){yk∈με∂∥xk∥r;xk+1=argminx∈IRn(12∥(I−PQ)Ax∥2+μ∥x∥1−(με∥xk∥r+〈yk,x−xk〉)).

where

[yk]i={sign([xk]i)|[xk]i|r−1∥xk∥rr−1if  xk≠0;0otherwise.

For the details of DCA convergence properties, see [15].

5. Concluding remarks

The focus of this paper is on Q-Lasso relying on two new DC-penalty methods instead of conventional methods such as l1 or l1−l2 minimization developed in [13,17] and [21]. Two iterative minimization methods based on the gradient proximal algorithm as well as the majored penalty algorithm are designed and their convergence to a stationary point is proved. Furthermore, by means of DC (difference of convex) Algorithm, two other algorithms are devised and their convergence results are also stated.