In sparse linear regression, we seek to estimate a sparse vector β ∈ R p from
where A ∈ R n×p is the design matrix and w denotes the observation noise. In this paper, we study the sorted 1 penalization estimator (SLOPE) [1] (see also [2], [3]). Given a non-decreasing regularization sequence λ = [λ 1 , λ 2 , . . . , λ p ] with 0 ≤ λ 1 ≤ λ 2 ≤ • • • ≤ λ p , SLOPE estimates β by solving the following optimization problem
where
is a reordering of the absolute values |x 1 | , |x 2 | , . . . , |x p | in increasing order. In [1], the regularization term J λ (x) def = p i=1 λ i |x| (i) is referred to as the "sorted 1 norm" of x. The same regularizer was independently developed in a different line of work [2]- [4], where the motivation is to promote group selection in the presence of correlated covariates.
The classical LASSO estimator is a special case of SLOPE. It corresponds to using a constant regularization sequence, i.e., λ 1 = λ 2 = • • • = λ p = λ. However, with more general λ-sequences, SLOPE has the flexibility to penalize different coordinates of the estimate according to their magnitudes. This adaptivity endows SLOPE with some nice statistical properties that are not possessed by LASSO. For example, it is shown in [5], [6] that SLOPE achieves the minimax 2 estimation rate with high probability. In terms of testing, the authors of [1] show that SLOPE controls false discovery rate (FDR) for orthogonal design matrices, which is not the case in LASSO [7]. In addition, we note that the new regularizer J λ (x) is still a norm [1], [3]. Thus, the optimization problem associated with SLOPE remains convex, and it can be efficiently solved by using e.g., proximal gradient descent [1], [3].
In the aforementioned studies on analyzing SLOPE, the performance of the estimator is given in terms of non-asymptotic probabilistic bounds. Such bounds provide very limited information about how to optimally design the regularization sequence λ in different settings, an important open question in the literature. In this paper, we provide two main contributions:
1) We obtain a characterization of SLOPE in the asymptotic regime: n, p → ∞ and n/p → δ. Compared with the probabilistic bounds derived in previous work, our results are asymptotically exact. Similar asymptotic analysis has been done for LASSO [8] and many other regularized linear regression problems [9]- [11], but the main technical challenge in analyzing SLOPE is the nonseparability of J λ (x): it cannot be written as a sum of componentwise functions, i.e., J λ (x) = p i=1 J i (x i ). In our work, we overcome this challenge by showing that the proximal operator of J λ (x) is asymptotically separable.
2) Using our asymptotic characterization, we derive oracle optimal λ in two settings: (1) the optimal regularization sequence that minimizes the MSE E ββ 2 ; and (2) the optimal sequence that achieves the highest possible power in testing and variable selection under a given level of Type-I error. In both cases, we show that the optimal design can be recast as certain infinite-dimensional convex optimization problems, which have efficient and accurate finite-dimensional approximations.
A caveat of our optimal design is that it requires knowing the limiting empirical measure of β (e.g., the sparsity level and the distribution of its nonzero coefficients). For this reason, our results are oracle optimal. It provides the first step towards more practical optimal designs that are completely blind to β.
The rest of the paper is organized as follows. In Sec. II, we first prove the asymptotic separability of the proximal operator associated with J λ (x). This property allows us to derive our main asymptotic characterization of SLOPE, summarized as Theorem 1. Based on this analysis, we present the optimal design of the regularization sequence in Sec. III. Numerical simulations verify our asymptotic characterizations. They also demonstrate the superiority of our optimal design over LASSO and a previous sequence design in the literature [5]. Due to space constraints, we only state and illustrate the main results in this paper, and leave the technical proofs to [12].
There are four main objects in the description of our model and algorithm: (1) the unknown sparse vector β; (2) the design matrix A; (3) the noise vector w; and (4) the regularization sequence λ. Since we study the asymptotic limit (with p → ∞), we will consider a sequence of instances β (p) , A (p) , w (p) , λ (p) p∈N with increasing dimensions p, where β (p) , λ (p) ∈ R p , A (p) ∈ R n×p and w (p) ∈ R n . A sequence of vectors x (p) ∈ R p indexed by p is called a converging sequence [8] if its empirical measure µ p (x)
i ) converges weakly to a probability measure on R.
Our results are proved under the following assumptions: (A.1) The number of observations grows in proportion to p:
2) The number of nonzero elements in β (p) grows in proportion to p: k/p → ρ ∈ (0, 1].
3) The elements of A (p) are i.i.d. Gaussian distribution:
∼ N (0, 1 n ). (A.4) {β (p) } p , {w (p) } p and {λ (p) } p are converging sequences.
The distribution functions of the limiting measures are denoted by F β , F w and F λ , respectively. Moreover, we have
, where the probability P(•) and the expectations E[•] are all computed with respect to the limiting measures.
We start by studying the proximal operator associated with the sorted 1 norm J λ (x). Given y ∈ R p and a regularization sequence λ with 0 ≤ λ 1 ≤ λ 2 ≤ • • • ≤ λ p , the proximal operator is defined as the solution to the following optimization problem:
where
In the case of LASSO, which corresponds to choosing
the proximal operator is easy to characterize, as it is separable: [Prox λ (y)] i = sign(y i ) max(|y i |λ, 0). In other words, the ith element of Prox λ (y) is solely determined by y i . However, this separability property does not hold for a general regularization sequence. When p is finite, [Prox λ (y)] i depends not only on y i but also on other elements of y. This coupling makes it much harder to analyze the proximal operator. Fortunately, as we show below, when p → ∞, Prox λ (•) becomes asymptotically separable. Proposition 1: Let {y (p) } p and {λ (p) } p be two converging sequences. Denote by F y and F λ the distribution functions of their respective limiting measures. It holds that
where η(•; F y , F λ ) is a scalar function that is determined by F y and F λ , and η(y (p) ; F y , F λ ) denotes a coordinate-wise application of η(•; F y , F λ ) on y (p) . The asymptotic separability of Prox λ (•) greatly facilitates our asymptotic analysis and design of SLOPE, since it allows us to reduce the original high-dimensional problem to an equivalent one-dimensional problem. In what follows, we refer to η(•; F y , F λ ) as the limiting scalar function. In the appendix, we briefly present a procedure for constructing η(•; F y , F λ ) from the limiting distribution functions F y and F λ . More details can be found in [12].
Example 1: We compare the proximal operator Prox λ (y) and the limiting scalar function η(y; F y , F λ ), for two different λ-sequences shown in Fig. 1(a) and Fig. 1(c). The red curves represent the limiting scalar functions obtained in Proposition 1, whereas the blue circles are sample points of (y i , [Prox λ (y)] i ), with y i ∼ N (0, 1). For better visualization, we randomly sample 3% of all (y i , [Prox λ (y)] i ). It can be seen that under a moderate dimension p = 1024, the proximal operator can already be very accurately approximated by the limiting scalar function.
We are now ready to tackle the original optimization problem (2) associated with SLOPE. Our goal is to characterize the joint empirical measure of ( β, β):
Indeed, many quantities of interest, such as the MSE, type-I error, and FDR, are all functionals of this joint empirical measure. A function
where L is a positive constant. As in [8], we will depict the limit of µ p ( β, β) through its action on pseudo-Lipschiz functions.
Theorem 1: Assume (A.1) -(A.4) hold. For any pseudo-Lipschiz function ψ, we have
(5) Here, B, Z are two independent random variables with B ∼ F β and Z ∼ N (0, 1); η(• ; F y , F τ λ ) is the limiting scalar function defined in Proposition 1, with F y denoting the distribution function of B + σZ and F τ λ denoting that of τ λ for some τ ≥ 0. Moreover, the scalar pair (σ, τ ) is the unique solution of the following equations:
Remark 1: Readers familiar with the asymptotic analysis of LASSO will recognize that the forms of ( 6) and ( 7) look identical to the results of LASSO obtained in [8], [11]. Indeed, the first part of our proof directly applies the framework of analyzing LASSO asymptotics using convex Gaussian min-max theorem (CMGT) [10], [11]. Following [11], in the asymptotic regime, the limiting measure of SLOPE is determined by the following fixed point equations:
Note that ( 8) and ( 9) are similar to ( 6) and ( 7), except that they involve an R p → R p proximal mapping: Prox τ λ (β+σZ). This is where Proposition 1 becomes useful. Using the asymptotic separability stated in that proposition, we can simplify ( 8) and ( 9) to the scalar equations given in ( 6) and (7). Theorem 1 essentially says that the joint empirical measure of ( β, β) converges weakly to the law of (η(B + σZ; F y , F τ λ ), B). This means that although the original problem ( 2) is high-dimensional, its asymptotic performance can be succinctly captured by merely scalars random variables. In (5), if we let ψ(x, y) = (xy) 2 , we obtain the asymptotic MSE; by setting ψ(x, y) = 1 y=0, x =0 , we can recover the type-I error. (Technically, 1 y=0,x =0 is not pseudo-Lipschiz. However, with additional justifications [1], one can show that the conclusion is still correct.)
In this section, we will study the optimal design of the regularization sequence in SLOPE. Using the asymptotic characterization presented in Sec. II, we will derive the optimal limiting distribution F λ to achieve the best estimation or testing performance, given the oracle knowledge of F β .
We first turn to the problem of finding the optimal λsequence which minimizes the MSE of slope estimator. Since we work in the asymptotic regime, it boils down to finding the optimal distribution F * λ such that
= arg min
where B ∼ F β and the second equality follows from Theorem 1. From (6), this is further equivalent to finding F λ to minimize σ. However, directly searching over F λ appears unwieldy, since σ as a functional of F λ is defined indirectly through a nonlinear fixed point equation.
To simplify this problem, we first note that in ( 6) and ( 7), the influence of F λ on the solution (σ, τ ) is only through the limiting scalar function η. Therefore, instead of optimizing over F λ , we can find the optimal η * and then calculate the corresponding F * λ . The next question then becomes finding all possible η that can be realized by some F λ . In fact, we can compute all possible η(•) associated with any converging sequence y (p) p∈N . Let
) holds be the functional space that η belongs to. We have the following result: Proposition 2: For any converging sequence y (p) p∈N , we have
and M is a convex set. Moreover, for any η ∈ M, the corresponding distribution of the λ-sequence that yields η can be represented by: λ ∼ |Y |η(|Y |), where Y follows the limiting distribution of y (p) p∈N . Remark 2: Proposition 2 is the key ingredient in our optimal design. It shows that, with different choices of F λ , we can reach any nondecreasing and odd function that is Lipschitz continuous with constant 1. Clearly, the soft-thresholding functions associated with LASSO belongs to M, but the set M is much richer. This is the essence of how SLOPE generalizes LASSO: it allows for more degrees of freedom in the regularization.
Due to Proposition 2, the optimization problem can be simplified to that of finding the optimal η ∈ M such that σ as obtained from ( 6) and ( 7) is minimized. Specifically, we need to find 6) and ( 7)} . (10) Note that equations ( 6) and ( 7) involve two variables σ and τ . It is not easy to handle them simultaneously. A simplification we can make is to first set τ to 1 and find the minimum σ such that the first equation ( 6) and the inequality E B,Z [η (B + σZ; F y , F λ )] ≤ δ hold. Once we get σ min , the corresponding τ can then be obtained via (7):
) -1 and λ is in turn updated to be λ/τ . After this replacement, ( 6) and ( 7) will be both satisfied. It is not difficult to show that this procedure will lead to the same σ min as defined in (10). Therefore, the remaining task is to solve, for every candidate σ, the following problem:
Thanks to the convexity of M, we can show that ( 11) is an infinite-dimensional convex optimization problem. In practice, we can discretize over R to solve a finite-dimensional approximation. If the minimum value of ( 11) is smaller than δ(σ 2σ 2 w ), then it can be shown that σ min < σ and vice versa. Clearly, we only need to search for σ min over a compact set:
, since from (6), we know σ 2 ≥ σ 2 w and also
] to find the minimum σ. Once we find the optimal η(y) and (σ min , τ ), we know from Proposition 2 that the corresponding λ can be represented as:
In Fig. 2, we compare the MSEs achieved by our optimal design with those obtained by LASSO and the BHq sequences proposed [5], at different SNR and sparsity levels. For fair comparison, we optimize the parameters of the BHq and LASSO sequences. It can be seen from the figure that the empirical minimum MSEs match well with theoretical ones. We observe from Fig. 2a that, under low SNRs, the BHq sequence can lead to very similar performance as the oracle optimal design. However, at higher SNR levels, the optimal design outperforms the BHq sequence, while it gets close to LASSO. To unravel the underlying reason for this, we plot in Fig. 3 the distributions of the λ-sequences associated with the optimal design and the BHq design, respectively. It turns out that, in the low SNR case, the optimal design and BHq have similar distributions; at higher SNRs, the distribution of the optimal design is close to a delta-like distribution similar to LASSO.
Note that for small sparsity-level ρ, LASSO can outperform BHq and achieve performance close to that of the optimal sequence, but it is prone to higher bias when ρ grows. From Fig. 2b, we can find that LASSO's performance degrades much faster than the other two as ρ increases, This is because LASSO's penalization is not adaptive to the underlying sparsity levels [5].
Next we consider using SLOPE for variable selection. For a given level of type-I error, we want to find the optimal regularization sequence to achieve the highest possible power.
As we have shown in the last section, in the asymptotic region, this is equivalent to optimizing over η ∈ M. Let y thresh = sup y≥0 {y | η(y) = 0}. It follows from Theorem 1 that, in order to ensure P type-I = α, we need to have
, where Φ(•) is the CDF of the standard normal distribution. Similarly, we can compute the power of the test as
)). It can be shown that for any fixed β,
)) is a nonincreasing function of σ. Thus, under a given type-I error rate α, maximizing the power is equivalent to minimizing σ.
Similar to the procedure described in Sec. III-A, we can traverse through a bounded set of σ to find σ min . The difference here is that we need to enforce additional constraints that guarantee P type-I = α. We omit further details, which can be found in [12]. In Fig. 4, we compare the FDR curve of the optimal design with that of the BHq sequence. We verify that the optimal design indeed dominates the BHq sequence and that the empirical FDR curve matches well with theoretical prediction. ∼ (1ρ)δ(0) + ρN (µ 0 , σ 2 0 ) with ρ = 0.25, µ 0 = 2.125 and σ 0 = 0, w i i.i.d.
∼ N (0, σ 2 ) with σ = 0.25. The results are average over 100 realizations.
In this appendix, we briefly describe the procedure for constructing the limiting scalar function η(•; F y , F λ ) in Proposition 1. The procedure, as summarized in Algorithm 1, can be viewed as the asymptotic limit of a fast algorithm proposed in [1], [2] for solving (3). Here, λ(y) can be intuitively understood as a function that assigns each y a regularization strength λ and G(y) just thresholds y using the assigned λ. This is exactly in the same spirit of magnitude-dependent regularization applied in SLOPE. The WHILE LOOP in Step 2-14 is essentially an adjustment of G(y) obtained in Step 1: within certain intervals [y L , y R ], the original G(y) is replaced by E (G(y) | y ∈ [y L , y R ]). The WHILE LOOP ends until G(y) is nondecreasing.
For illustrations, we consider the simplest scenario when F Y is continuous and when G(y) as obtained in Step 1 is nondecreasing. In this case, WHILE LOOP in Step 2-14 will not be executed. It follows that η(y; F y , F λ ) = sign(y) max{0, |y| -F -1 λ [F |y| (|y|)]}. Clearly, this reduces to the soft-thresholding function associate with LASSO, with the regularization parameter given by F -1 λ [F |y| (|y|)] ≡ λ. Illustrations of several more complicated examples can be found in [12]. b Maximal nondecreasing interval (MNDI) is defined in the same way as MDI.