Parameter Estimation with the Ordered ℓ₂ Regularization via an Alternating Direction Method of Multipliers

Abstract

Regularization is a popular technique in machine learning for model estimation and for avoiding overfitting. Prior studies have found that modern ordered regularization can be more effective in handling highly correlated, high-dimensional data than traditional regularization. The reason stems from the fact that the ordered regularization can reject irrelevant variables and yield an accurate estimation of the parameters. How to scale up the ordered regularization problems when facing large-scale training data remains an unanswered question. This paper explores the problem of parameter estimation with the ordered ${\ell }_2$-regularization via Alternating Direction Method of Multipliers (ADMM), called ADMM-O${\ell }_2$. The advantages of ADMM-O${\ell }_2$ include (i) scaling up the ordered ${\ell }_2$ to a large-scale dataset, (ii) predicting parameters correctly by excluding irrelevant variables automatically, and (iii) having a fast convergence rate. Experimental results on both synthetic data and real data indicate that ADMM-O${\ell }_2$ can perform better than or comparable to several state-of-the-art baselines.

1. Introduction

In the machine learning literature, one of the most important challenges involves estimating parameters accurately and selecting relevant variables from highly correlated, high-dimensional data. Researchers have noticed many highly correlated features in high-dimensional data [1]. Models often overfit or underfit high-dimensional data because they have a large number of variables but only a few of them are actually relevant; most others are irrelevant or redundant. An underfitting model contributes to estimation bias (i.e., high bias and low variance) in the model fitting because it keeps out relevant variables, whereas an overfitting raises estimation error (i.e., low bias and high variance) since it includes irrelevant variables in the model.

To illustrate an application of our proposed method, consider a study of gene expression data. This is a high-dimensional dataset and contains highly correlated genes. The geneticist always likes to determine which variants/genes contribute to changes in biological phenomena (e.g., increases in blood cholesterol level, etc.) [2]. Therefore, the aim is to explicitly identify all relevant variants. The penalized regularization models such as ${\ell }_1$, ${\ell }_2$, and so forth have recently become a topic of great interest within machine learning, statistics [1], and optimization [3] communities as classic approaches to estimate parameters. The ${\ell }_1$-based method is not a preferred selection method for groups of variables among which pairwise correlations are significant because the lasso arbitrarily selects a single variable from the group without any consideration of which one to select [4]. Furthermore, if the selected value of a parameter is too small, the ${\ell }_1$-based method would select many irrelevant variables, thus degrading its performance. On the other hand, a large value of the parameter would yield a large bias [5]. Another point worth noting is that few ${\ell }_1$ regularization methods are either adaptive, computationally tractable, or distributed, but no such method contains all three properties together. Therefore, the aim of this study is to develop a model for parameter estimation and to determine relevant variables in highly correlated, high-dimensional data based on the ordered ${\ell }_2$. This model has the following three properties all together: adaptive (Our method is adaptive in the sense that it reduces the cost of including new relevant variables as more variables are added to the model due to rank-based penalization properties.), tractable (A computationally intractable method is a computer algorithm that takes a very long time to execute a mathematical solution. A computationally tractable method is exactly the opposite of intractable method.), and distributed.

Several adaptive and nonadaptive methods have been proposed for parameter estimation and variable selection in large-scale datasets. Different principles are adopted in these procedures to estimate parameters. For example, an adaptive solution (i.e., an ordered ${\ell }_1$) [5] is a norm and, therefore, convex. Regularization parameters are sorted in non-increasing order in the ordered ${\ell }_1$, in which the ordered regularization penalizes regression coefficients according to their order, with higher orders closer to the top and having larger penalties. Pan et al. [6] proposed a partial sorted ${\ell }_p$ norm, which is non-convex and non-smooth. In contrast, the ordered ${\ell }_2$ regularization is convex and smooth, just as the standard ${\ell }_2$ norm is convex and smooth in Reference [7]. Pan et al. [6] considered p-values between $0<p\le 1$ that do not cover ${\ell }_2$, ${\ell }_{\infty }$ norms, and so forth. Pan et al. [6] also did not provide details of other partially sorted norms when $p\ge 2$ and used random projection and the partial sorted ${\ell }_p$ norm to complete the parameter estimation, whereas we have used ADMM with the ordered ${\ell }_2$. A nonadaptive solution (i.e., an elastic net) [8] is a mixture of both ordinary ${\ell }_1$ and ${\ell }_2$. In particular, it is a useful model when the number of predictors (p) is much larger than the number of observations (n) or in any situation where the predictor variables are correlated.

Table 1. Properties of the different regularizers: Correlation is a method that shows how strong the relationship between explanatory variables is. In the correlation column, “Yes” means a strong correlation between variables and “No” means weaker (or no) correlation between variables. Stable means that estimates of the ${\ell }_2$-based methods are more stable when the explanatory variables are strongly correlated.

Regularizers	Promoting	Convex	Smooth	Adaptive	Tractable	Correlation	Stable
${\ell }_0$ [10,11]	Sparsity	No	No	No	No	No	No
${\ell }_1$ [12]	Sparsity	Yes	No	No	Yes	No	No
The Ordered ${\ell }_1$ [5]	Sparsity	Yes	No	Yes	Yes	No	No
${\ell }_2$ [13]	Grouping	Yes	Yes	No	Yes	Yes	Yes
The Ordered ${\ell }_2$	Grouping	Yes	Yes	Yes	Yes	Yes	Yes
Partial Sorted ${\ell }_p$ [6]	Sparsity	No	No	Yes	No	No	No

presents the important properties of the regularizers. As seen in Table 1, ${\ell }_2$ and the ordered ${\ell }_2$ regularizers are suitable methods for highly correlated, high-dimensional grouping data rather than ${\ell }_1$ and the ordered ${\ell }_1$ regularizers. The ordered ${\ell }_2$ encourages grouping, whereas most of ${\ell }_1$-based methods promote sparsity. Here, grouping signifies a group of strongly correlated variables in high-dimensional data. We used the ordered ${\ell }_2$ regularization in our method instead of ${\ell }_2$ regularization because the ordered ${\ell }_2$ regularization is adaptive. Finally, ADMM has a parallel behavior for solving large-scale convex optimization problems. Our model also employs ADMM and inherits distributed properties of native ADMM [9]. Hence, our model is also distributed. Bogdan et al. [5] did not provide any details about how they applied ADMM in the ordered ${\ell }_1$ regularization.

In this paper, we propose “Parameter Estimation with the Ordered ${\ell }_2$ Regularization via ADMM” called ADMM-O${\ell }_2$ to find relevant parameters from a model. ${\ell }_2$ is a ridge regression; similarly, the ordered ${\ell }_2$ becomes an ordered ridge regression. The main contribution of this paper is not to present a superior method but rather to introduce a quasi-version of the ${\ell }_2$ regularization method and to concurrently raise awareness of the existing methods. As part of this research, we introduced a modern ordered ${\ell }_2$ regularization method and proved that the square root of the ordered ${\ell }_2$ is a norm and, thus, convex. Therefore, it is also tractable. In addition, the regularization method used an ordered elastic net method to combine the widely used ordered ${\ell }_1$ penalty with modern ordered ${\ell }_2$ penalty for ridge regression. The ordered elastic net is also proposed by the scholars in this paper. To the best of our knowledge, this is one of the first method to use the ordered ${\ell }_2$ regularization with ADMM for parameter estimation and variable selection. In Section 3 and Section 4, we explain the integration of ADMM with the ordered ${\ell }_2$ and further details about it.

The rest of the paper is arranged as follows. Related works are discussed in Section 2, along with a presentation of the ordered ${\ell }_2$ regularization in Section 3. Section 4 describes the application of ADMM to the ordered ${\ell }_2$. Section 5 presents the experiments conducted. Finally, Section 6 closes the paper with a conclusion.

2. Related Work

2.1. ${\ell }_1$ and ${\ell }_2$ Regularization

Deng et al. [14] presented efficient algorithms for group sparse optimization with mixed ${\ell }_{2,1}$ regularization for the estimation and reconstruction of signals. Their technique is rooted in a variable splitting strategy and ADMM. Zou and Hastie [8] suggested that an elastic net, a generalization of the lasso, is a linear combination of both ${\ell }_1$ and ${\ell }_2$ norm. It contributes to sparsity without permitting a coefficient to become too large. However, Candes and Tao [15] introduced a new estimator called the Dantzig selector for a linear model when the parameters are larger than the number of observations for which they established optimal ${\ell }_2$ rate properties under a sparsity. Chen et al. [16] enforced sparse embedding to ridge regression, obtaining solutions $\tilde{x}$ with ${\parallel \tilde{x}-x^{\ast }\parallel }_2\le ϵ\parallel x^{\ast }\parallel$ small, where $x^{\ast }$ is optimal, and also did this in $\mathcal{O}(nnz(A)+n^3/{ϵ}^2)$ time, where nnz(A) is the number of nonzero entries of A. Recently, Bogdan et al. [5] have proposed an ordered ${\ell }_1$-regularization technique inspired from a statistical viewpoint, in particular, by a focus on controlling the false discovery rate (FDR) for variable selection in linear regressions. Our proposed method is similar but focused on parameter estimation based on the ordered ${\ell }_2$ regularization and ADMM. Several methods have been proposed based on Reference [5] and similar ideas. For example, Bogdan et al. [2] introduced a new model-fitting strategy called Sorted L-One Penalized Estimation (SLOPE).which regularizes least-squares estimates with rank-dependent penalty coefficients. Zeng and Figueiredo [17] proposed DWSL1 as a generalization of octagonal shrinkage and clustering algorithm (OSCAR) that aims to promote feature grouping without previous knowledge of the group structure. Pan et al. [6] have introduced an image restoration method based on a random projection and a partial sorted ${\ell }_p$ norm. In this method, an input signal is decomposed into two components: a low-rank component and a sparse component. The low-rank component is approximated by random projection, and the sparse one is recovered by the partial sorted ${\ell }_p$ norm. Our method can potentially be used in various other domains such as cyber security [18] and recommendation [19].

2.2. ADMM

Researchers have paid a significant amount of attention to ADMM because of its capability of dealing with objective functions independently and simultaneously. Second, ADMM has proved to be a genuine fit in the field of large-scale data-distributed optimization. However, ADMM is not a new algorithm; it was first introduced by References [20,21] in the mid-1970s, with roots as far back as the mid-1950s. In addition, ADMM originated from an augmented method with a Lagrangian multiplier [22]. It became more popular when Boyd et al. [9] published papers about ADMM. The classic ADMM algorithm applies to the following “ADMM-ready” form of problems.

$$\{\begin{array}{c}minimizef(x)+g(z)\hfill \\ s.t.Ax+Bz=c\hfill \end{array}$$

(1)

The wide range of applications have also inspired the study of the convergence properties of ADMM. Under mild assumptions, ADMM can converge for all choices of the step size. Ghadimi et al. [23] provided some advice on tuning over-relaxed ADMM in the quadratic problems. Deng and Yin [24] have also suggested linear convergence results under the consideration of only a single strongly convex term, given that linear operator’s A and B are full-rank matrices. These convergence results bound error as measured by an approximation to the primal–dual gap. Goldstein et al. [25] created an accelerated version of ADMM that converges more quickly than traditional ADMM under an assumption that both objective functions are strongly convex. Yan and Yin [26] explained, in detail, the different kinds of convergence properties of ADMM and their prerequisite rules for converging. For further studies on ADMM, see Reference [9].

3. Ridge Regression with the Ordered ${\ell }_2$ Regularization

3.1. The Ordered ${\ell }_2$ Regularization

The proposed parameter estimation and variable selection method in this paper is computationally manageable and adaptive. This procedure depends on the ordered ${\ell }_2$ regularization. Let $\lambda =({\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_p)$ be a decreasing sequence of positive scalars that satisfy the following condition:

$${\lambda }_1\ge {\lambda }_2\ge {\lambda }_3\ge \dots \ge {\lambda }_p\ge 0$$

(2)

The ordered ${\ell }_2$ regularization of a vector $x\in {\mathbb{R}}^p$ when ${\lambda }_1>0$ can be defined as follows:

$$J_{\lambda }(x)={\lambda }_1{x}_{(1)}^2+{\lambda }_2{x}_{(2)}^2+\cdots +{\lambda }_p{x}_{(p)}^2=\sum _{k=1}^p{\lambda }_{{BH}^{(k)}}{x}_{(k)}^2$$

(3)

where ${\lambda }_{{BH}^{(k)}}$ is called a BHq method [27], which generates an adaptive and a non-increasing value for $\lambda$ (Reference [2], Section 1.1). The details of ${\lambda }_{{BH}^{(k)}}$ are available in Section 4.2. For ease of presentation, we have written ${\lambda }_k$ in place of ${\lambda }_{{BH}^{(k)}}$ in the rest of the paper. $x_{(1)}^2\ge x_{(2)}^2\ge x_{(3)}^2\ge \cdots \ge x_{(p)}^2$ is the order statistic of the magnitudes of x [28]. The subscript k of x enclosed in parentheses indicates the kth-order statistic of a sample. Suppose that x is a sample size of 4. Hence, four numbers are observed in x if the sample values of x = (−2.1, −0.5, 3.2, 7.2). The order statistics of x would be $x_{(1)}^2={7.2}^2,x_{(2)}^2={3.2}^2,x_{(3)}^2={2.1}^2,x_{(4)}^2={0.5}^2$. The ordered ${\ell }_2$ regularization is expressed as the first largest value of $\lambda$ times the square of the first largest entry of x, plus the second largest value of $\lambda$ times the square of the second largest entry of x, and so on. $A\in {\mathbb{R}}^{n\times p}$ and $b\in {\mathbb{R}}^n$ are a matrix and a vector, respectively. The ordered ${\ell }_2$ regularized loss minimization can be expressed as follows:

$${\displaystyle \stackrel{min}{x\in {\mathbb{R}}^p}}\frac{1}{2}{\parallel Ax-b\parallel }_2^2+\frac{1}{2}\{{\lambda }_1{x}_{(1)}^2+\cdots +{\lambda }_p{x}_{(p)}^2\}$$

(4)

Theorem 1.

If the square root of$J_{\lambda }(x)$(Equation (3)) is a norm on${\mathbb{R}}^p$and a function$\parallel .\parallel :{\mathbb{R}}^p\to \mathbb{R}$satisfying the following three properties, then the following Corollaries 1 and 2 are true.

i 

(Positivity)$\parallel x\parallel \ge 0$for any$x\in {\mathbb{R}}^p$and$\parallel x\parallel =0$if and only if$x=0$.

ii 

(Homogeneity)$\parallel cx\parallel =\left|c\right|\parallel x\parallel$for any$x\in {\mathbb{R}}^p$and$c\in \mathbb{R}$.

iii 

(Triangle inequality)$\parallel x+y\parallel \le \parallel x\parallel +\parallel y\parallel$for any$x,y\in {\mathbb{R}}^p$.

Note:$\parallel x\parallel$and${\parallel x\parallel }_2$are used interchangeably.

Corollary 1.

When all the${\lambda }_k^{\prime }s$take on an equal positive value,$J_{\lambda }(x)$reduces to the square of the usual${\ell }_2$norm.

Corollary 2.

When${\lambda }_1>0$and${\lambda }_2$ = … = ${\lambda }_p=0$, the square root of $J_{\lambda }(x)$ reduces to the ${\ell }_{\infty }$ norm.

Proofs of the theorem and corollaries are provided in Appendix A.

Table 2. Notations and explanations.

Notations	Explanations	Notations	Explanations
Matrix	denoted by uppercase letter	f	loss convex function
Vector	denoted by lowercase letter	g	Regularizer part (${\ell }_1$ or ${\ell }_2$ etc.)
${\parallel .\parallel }_1$	${\ell }_1$ norm	$\partial f(x)$	subdifferential of convex function f at x
${\parallel .\parallel }_2$	${\ell }_2$ norm	$\partial g(z)$	subdifferential of convex function g at z
${\parallel x\parallel }_1={\displaystyle \sum _{k=1}^n}\left|x_k\right|$	${\ell }_1$ norm	L1	the ${\ell }_1$ norm or the lasso
${\parallel x\parallel }_2=\sqrt{{\displaystyle \sum _{k=1}^n}{\left|x_k\right|}^2}$	${\ell }_2$ norm	OL1	ordered ${\ell }_1$ norm or ordered lasso
${\parallel x\parallel }_2^2={\displaystyle \sum _{k=1}^n}{\left|x_k\right|}^2$	square of ${\ell }_2$ norm	OL2	ordered ${\ell }_2$ norm or ordered ridge regression
$J_{\lambda }(.)$	the ordered norm	Eq.	equation

Note: We often use the ordered ${\ell }_2$ norm/regularization, OL2, and ADMM-O${\ell }_2$ interchangeably.

shows the notations used in this paper and their meanings.

3.2. The Ordered Ridge Regression

We propose an ordered ridge regression in Equation (5) and call it the ordered ridge regression because we used an ordered ${\ell }_2$ regularization with an objective function instead of a standard ${\ell }_2$ regularization. The ordered ridge regression is commonly used for parameter estimation and variable selection, particularly when data are strongly correlated and highly dimensional. The ordered ridge regression can be defined as follows:

$$\begin{array}{c}\hfill {\displaystyle \stackrel{min}{x\in {\mathbb{R}}^p}}\frac{1}{2}{\parallel Ax-b\parallel }_2^2+\frac{1}{2}{J}_{\lambda }(x)\\ \hfill ={\displaystyle \stackrel{min}{x\in {\mathbb{R}}^p}}\frac{1}{2}{\parallel Ax-b\parallel }_2^2+\frac{1}{2}\sum _{k=1}^p{\lambda }_k{\left|x_{(k)}\right|}^2\end{array}$$

(5)

where $x\in {\mathbb{R}}^p$ denotes an unknown regression coefficient, $A\in {\mathbb{R}}^{n\times p}(p\gg n)$ is a known matrix, $b\in {\mathbb{R}}^n$ represents a response vector, and $J_{\lambda }(x)$ is the ordered ${\ell }_2$ regularization. The optimal parameter choice for the ordered ridge regression is much more stable than that for a regular lasso; also, it achieves adaptivity in the following senses.

i

For decreasing $({\lambda }_k)$, each parameter ${\lambda }_k$ marks the entry or removal of some variable from the current model (therefore, its coefficient becomes either zero or nonzero); thus, coefficients remain constant in the model. We achieved this by putting some threshold values for ${\lambda }_k$ (Reference [5], Section 1.4).

ii

We observed that the price of including new variables declines as more variables are added to the model when the ${\lambda }_k$ decreases.

4. Applying ADMM to the Ordered Ridge Regression

In order to apply ADMM to the problem in Equation (5), we first transform it into an equivalent form of the problem in Equation (1) by introducing an auxiliary variable z.

$${\displaystyle \stackrel{min}{x,z\in {\mathbb{R}}^p}}\frac{1}{2}{\parallel Ax-b\parallel }_2^2+\frac{1}{2}{J}_{\lambda }(z)s.t.x-z=0$$

(6)

We can see that Equation (6) has two blocks of variables (i.e., x and z). Its objective function is separable in the form of Equation (1) since $f(x)=\frac{1}{2}{\parallel Ax-b\parallel }_2^2$ and $g(z)=\frac{1}{2}{J}_{\lambda }(z)=\frac{1}{2}\lambda {\parallel z\parallel }_2^2=\frac{1}{2}{\displaystyle \sum _{k=1}^p}{\lambda }_k{\left|x_{(k)}\right|}^2$, where $A=I$ and $B=-I$. Therefore, ADMM is applicable to Equation (6). An augmented Lagrangian form of Equation (6) can be defined as follows:

$$\begin{array}{cc}\hfill {\mathcal{L}}_p(x,z,y)& =\frac{1}{2}{\parallel Ax-b\parallel }_2^2+\frac{1}{2}\lambda {\parallel z\parallel }_2^2+y^T(x-z)\hfill \\ & +\frac{\rho }{2}{\parallel x-z\parallel }_2^2\hfill \end{array}$$

(7)

where $y\in {\mathbb{R}}^p$ is a Lagrangian multiplier and $\rho >0$ denotes a penalty parameter. Next, we apply ADMM to the augmented Lagrangian equation of Equation (7) (Reference [9], Section 3.1), which renders ADMM iterations as follows:

$$\begin{array}{c}{x}^{k+1}:={\displaystyle \stackrel{argmin}{x\in {\mathbb{R}}^p}}{\mathcal{L}}_p(x,z^k,y^k)\hfill \\ z^{k+1}:={\displaystyle \stackrel{argmin}{z\in {\mathbb{R}}^p}}{\mathcal{L}}_p(x^{k+1},z,y^k)\hfill \\ y^{k+1}:=y^k+\rho (x^{k+1}-z^{k+1})\hfill \end{array}$$

(8)

Proximal gradient methods are well known for solving convex optimization problems for which the objective function is the sum of a smooth loss function and a non-smooth penalty function [9,29,30]. A well-studied example is ${\ell }_1$ regularized least squares [1,5]. It should be noted that an ordered ${\ell }_1$ norm is convex but not smooth. Therefore, these researchers used a proximal gradient method. In contrast, we have employed an ADMM method because ADMM can solve convex optimization problems for which the objective function is either the sum of a smooth loss function and a non-smooth penalty function or both loss and penalty function are smooth and ADMM also supports parallelism. In the ordered ridge regression, both loss and penalty function are smooth, whereas, in the ordered elastic net, loss function is smooth and penalty function is non-smooth.

4.1. Scaled Form

We can also define ADMM in scaled form by merging a linear and a quadratic term in augmented Lagrangian and then a scaled dual variable, which is shorter and more appropriate. The scaled dual form of ADMM iterations in Equation (8) can be expressed as follows:

$$\begin{array}{}(9\mathrm{a})& x^{k+1}:={\displaystyle \stackrel{argmin}{x\in {\mathbb{R}}^p}}(f(x)+(\rho /2){\parallel x-z^k+u^k\parallel }_2^2)\hfill (9\mathrm{b})& z^{k+1}:={\displaystyle \stackrel{argmin}{z\in {\mathbb{R}}^p}}(g(z)+(\rho /2){\parallel x^{k+1}-z+u^k\parallel }_2^2)\hfill (9\mathrm{c})& u^{k+1}:=u^k+x^k-z^k\hfill \end{array}$$

where $u=\frac{1}{\rho }y$ and u is the scaled dual variable. Next, we can minimize the augmented Lagrangian in Equation (7) with respect to x and z, successively. Minimizing Equation (7) with respect to x becomes the x subproblem of Equation (9a), and it can be expressed as follows:

$$\begin{array}{}& \hfill {\displaystyle \stackrel{min}{x\in {\mathbb{R}}^p}}\frac{1}{2}{\parallel Ax-b\parallel }_2^2+(\rho /2){\parallel x-z^k+u^k\parallel }_2^2(10\mathrm{a})& \hfill ={\displaystyle \stackrel{min}{x\in {\mathbb{R}}^p}}\frac{1}{2}\{x^T{A}^{T}Ax-2b^{T}Ax\}+\frac{\rho }{2}\{x^2-2{(z^k-u^k)}^{T}x\}\end{array}$$

After computing a derivative of Equation (10a) with respect to x, then the setting of the derivative of x becomes equal to zero. Notice that this is a convex problem; therefore, it minimizes to solve the following linear system of Equation (10b):

$$\begin{array}{}& \hfill \iff A^{T}Ax+\rho x-b^{T}A-\rho {(z^k-u^k)}^T=0(10\mathrm{b})& \hfill x^{k+1}={(A^{T}A+\rho \ast I)}^{-1}{(A^{T}b+\rho (z^k-u^k))}^T\end{array}$$

Minimizing problem Equation (7) w.r.t. z, we obtain Equation (9b), and it results in the following z subproblem:

$$\begin{array}{}& \hfill {\displaystyle \stackrel{min}{z\in {\mathbb{R}}^p}}\frac{1}{2}\lambda {\parallel z\parallel }_2^2+\frac{\rho }{2}{\parallel x^{k+1}+u^k-z\parallel }_2^2(11\mathrm{a})& \hfill ={\displaystyle \stackrel{min}{z\in {\mathbb{R}}^p}}\frac{1}{2}{\lambda }_k{z}^{T}z+\frac{\rho }{2}\{{(x^{k+1}+u^k-z)}^T(x^{k+1}+u^k-z)\}\end{array}$$

After computing a derivative of Equation (11a) with respect to z, then the setting of the derivative of z becomes equal to zero. Notice that this is a convex problem; therefore, it minimizes to solve the following linear system of Equation (11b):

$$\begin{array}{}& \hfill \iff \frac{1}{2}2{\lambda }_{k}z+\frac{\rho }{2}\{2z-2{(x^{k+1}+u^k)}^T\}=0(11\mathrm{b})& \hfill z^{k+1}={({\lambda }_k+\rho \ast I)}^{-1}\rho (x^{k+1}+u^k)\end{array}$$

Finally, the multiplier (i.e., the scaled dual variable u) is updated in the following way:

$$u^{k+1}=u^k+(x^k-z^k)$$

(12)

Optimality conditions: Primal and dual feasibility are essential and adequate optimality conditions for ADMM in Equation (6) [9]. Dual residual ($S^{k+1}$) and primal residual (${\gamma }^{k+1}$) can be defined as follows:

$$\begin{array}{cccc}\mathrm{Dual}\mathrm{residual}(S^{k+1})\mathrm{at}\mathrm{iteration}k+1\hfill & & & \hfill =\rho (z^k-z^{k+1})\\ \mathrm{Primal}\mathrm{residual}({\gamma }^{k+1})\mathrm{at}\mathrm{iteration}k+1\hfill & & & \hfill =x^{k+1}-z^{k+1}\end{array}$$

Stopping criteria: The stopping criterion for an ordered ridge regression is that primal and dual residuals must be small:

$$\begin{array}{ccc}{\parallel {\gamma }^k\parallel }_2\le {ϵ}^{pri}\hfill & & {\displaystyle where\text{−}{ϵ}^{pri}=\sqrt{p}{ϵ}^{abs}+{ϵ}^{rel}max\{{\parallel x^k\parallel }_2,{\parallel z^k\parallel }_2\}}\\ {\parallel S^k\parallel }_2\le {ϵ}^{dual}\hfill & & {\displaystyle where{ϵ}^{dual}=\sqrt{n}{ϵ}^{abs}+{ϵ}^{rel}{\parallel \rho \ast u^k\parallel }_2}\end{array}$$

We set ${ϵ}^{abs}={10}^{-4}$ and ${ϵ}^{rel}={10}^{-2}$. For further details about this choice, see Reference [9], Section 3).

4.2. Over-Relaxed ADMM Algorithm

By comparing Equations (1) and (6), we can write Equation (6) in the over-relaxation form as follows:

$$\alpha x^{k+1}+(1-\alpha )z^k{\displaystyle /*\mathrm{where}\mathrm{A}=\mathrm{I},\mathrm{B}=-\mathrm{I}\mathrm{and}\mathrm{c}=0\mathrm{in}\mathrm{our}\mathrm{case}\ast /}$$

(13)

Substituting $x^{k+1}$ with Equation (13) into z of Equation (11b) and u of Equation (12) updates the results in a relaxation form. Algorithm 1 presents an ADMM iteration for the ordered ridge regression of Equation (6).

Algorithm 1: Over-relaxed ADMM for the ordered ridge regression

1: Initialize $x^0\in {\mathbb{R}}^p$, $z^0\in {\mathbb{R}}^p$, $u^0\in {\mathbb{R}}^p$, $\rho >0$ 2: while (${\parallel {\gamma }^k\parallel }_2\le {ϵ}^{pri}\&\&{\parallel S^k\parallel }_2\le {ϵ}^{dual}$) do 3:  $x^{k+1}\leftarrow {(A^{T}A+\rho \ast I)}^{-1}{(A^{T}b+\rho (z^k-u^k))}^T$ 4:  ${\lambda }_k\leftarrow SortedLambda(\{{\lambda }_k\})$; $▹{\displaystyle \mathrm{refer}\mathrm{to}\mathrm{Algorithm}2}$ 5:  $z^{k+1}\leftarrow {({\lambda }_k+\rho \ast I)}^{-1}\rho (\alpha x^{k+1}+(1-\alpha )z^k+u^k)$ 6:  $u^{k+1}\leftarrow u^k+\alpha (x^{k+1}-z^{k+1})+(1-\alpha )(z^k-z^{k+1})$ 7: end while

We observed that ADMM Algorithm 1 computes an exact solution for each subproblem and that their convergence is guaranteed by existing ADMM theory [24,25,31]. The most important and computationally intensive operation here is matrix inversion in line 3 of Algorithm 1. Here, matrix A is high-dimensional ($p\gg n$) and $(A^{T}A+\rho \ast I)$ takes $\mathcal{O}(n{p}^2)$ and its inverse (i.e., ${(A^{T}A+\rho \ast I)}^{-1}$) takes $\mathcal{O}(p^3)$. We compute ${(A^{T}A+\rho \ast I)}^{-1}$ and $A^{T}b$ outside loop; then, we are left with (inverse * ${(A^{T}b+\rho (z^k-u^k))}^T$), which is $\mathcal{O}(p^2)$, while addition and subtraction take $\mathcal{O}(p)$. ${(A^{T}A+\rho \ast I)}^{-1}$ is also cacheable, so the complexity is just $\mathcal{O}(p^3)$ + k * $\mathcal{O}(n{p}^2+p)$ heuristically with k number of iteration.

Generating the ordered parameter$({\lambda }_k)$: As mentioned in the beginning, we set out to identify a computationally tractable and adaptive solution. The regularizing sequences play a vital role in achieving this goal. Therefore, we generated adaptive values of $({\lambda }_k)$ such that regressor coefficients are penalized according to their respective order. Our regularizing sequence procedure is motivated by the BHq procedure [27]. The BHq method generates $({\lambda }_k)$ sequences as follows:

$${\lambda }_{{BH}^{(k)}}={\Phi }^{-1}(1-q\ast \frac{k}{2p})$$

(14)

$${\lambda }_k={\lambda }_{{BH}^{(k)}}\sqrt{1+\frac{\sum _{j<k}{\lambda }_{{BH}^{(j)}}^2}{n-k}}$$

(15)

where $k>0$, ${\Phi }^{-1}(\alpha )$ is $\alpha$th quantile of a standard normal distribution, and q is a parameter, namely $q\in [0;1]$. We started with ${\lambda }_1={\lambda }_{{BH}^{(1)}}$ as an initial value of the ordered parameter $({\lambda }_k)$.

Algorithm 2 presents a method for generating sorted $({\lambda }_k)$. The difference between lines 5 and 6 in Algorithm 2 is that line 5 is for low-dimensional ($n\le p$) data and that line 6 is for high-dimensional data ($p\gg n$). Finally, we used the ordered $({\lambda }_k)$ from Algorithm 2 (i.e., the adaptive value of $({\lambda }_k)$) in the ordered ridge regression of Equations (6) and (7) instead of ordinary $\lambda$. This makes the ordered ${\ell }_2$ adaptive and different from standard ${\ell }_2$.

Algorithm 2: Sorted Lambda $(\{{\lambda }_k\})$

1: Initialize $q\in [0;1],k>0,p,n\in N$ 2: ${\lambda }_1\leftarrow {\lambda }_{{BH}^{(1)}}$; $▹{\displaystyle {\lambda }_{{BH}^{(1)}}isfromEquation(14)wherek=1}$ 3: for$k\in \{2,\cdots ,K\}$do 4:  ${\lambda }_{{BH}^{(k)}}\leftarrow {\mathsf{\Phi }}^{-1}(1-q\ast \frac{k}{2p})$ 5:  ${\lambda }_k\leftarrow {\lambda }_{{BH}^{(k)}}\ast \sqrt{1+\frac{{\displaystyle \sum _{j<k}}{\lambda }_{{BH}^{(j)}}^2}{p-k-1}}$; $▹{\displaystyle \mathrm{where}\mathrm{n}=\mathrm{p}}$ 6:  ${\lambda }_k\leftarrow {\lambda }_{{BH}^{(k)}}\ast \sqrt{1+\frac{{\displaystyle \sum _{j<k}}{\lambda }_{{BH}^{(j)}}^2}{2p-k-1}}$ ; $▹{\displaystyle \mathrm{where}\mathrm{n}=2\mathrm{p}}$ 7: end for

4.3. The Ordered Elastic Net

A standard ${\ell }_2$ (or an ordered ${\ell }_2$) regularization is a commonly used tool to estimate parameters for microarray datasets (strongly correlated grouping). However, a key drawback of the ${\ell }_2$ regularization is that it cannot automatically select relevant variables because the ${\ell }_2$ regularization shrinks coefficient estimates closer but not exactly equal to zero (Reference [12], Chapter 6.2). On the other hand, a standard ${\ell }_1$ (or an ordered ${\ell }_1$) regularization can automatically determine relevant variables due to its sparsity property. However, the ${\ell }_1$ regularization also has a limitation. Especially when different variables are highly correlated, the ${\ell }_1$ regularization tends to pick only a few of them and to remove the remaining ones—even important ones that might be better predictors. To overcome the limitations of both ${\ell }_1$ and ${\ell }_2$ regularization, we proposed another method called an ordered elastic net (the ordered ${\ell }_{1,2}$ regularization or ADMM-O${\ell }_{1,2}$ or -O${\ell }_{1,2}$), similar to a standard elastic net [8], by combining the ordered ${\ell }_2$ regularization with the ordered ${\ell }_1$ regularization and the elastic net. By doing so, the ordered ${\ell }_{1,2}$ regularization automatically selects relevant variables in a way similar to the ordered ${\ell }_1$ regularization. In addition, it can select groups of strongly correlated variables. The key difference between the ordered elastic net and the standard elastic net is a regularization term. We apply the ordered ${\ell }_1$ and ${\ell }_2$ regularization in the ordered elastic net instead of the standard ${\ell }_1$ and ${\ell }_2$ regularization. This approach means that the ordered elastic net inherits the sparsity, grouping, and adaptive properties of the ordered ${\ell }_1$ and ${\ell }_2$ regularization. We have also employed ADMM to solve the ordered ${\ell }_{1,2}$ regularized loss minimization as follows:

$$\iff {\displaystyle \stackrel{min}{x\in {\mathbb{R}}^p}}\frac{1}{2}{\parallel Ax-b\parallel }_2^2+\alpha {\lambda }_{BH}{\parallel x\parallel }_1+\frac{1}{2}(1-\alpha ){\lambda }_{BH}{\parallel x\parallel }_2^2$$

For simplicity, let ${\lambda }_1=\alpha {\lambda }_{BH}$ and ${\lambda }_2=(1-\alpha ){\lambda }_{BH}$. The ordered elastic net becomes

$$\iff {\displaystyle \stackrel{min}{x\in {\mathbb{R}}^p}}\frac{1}{2}{\parallel Ax-b\parallel }_2^2+{\lambda }_1{\parallel x\parallel }_1+\frac{1}{2}{\lambda }_2{\parallel x\parallel }_2^2$$

Now, we can transform the above ordered elastic net equation into an equivalent form of Equation (1) by introducing an auxiliary variable z.

$$\iff {\displaystyle \stackrel{min}{x,z\in {\mathbb{R}}^p}}\underset{f(x)}{\underbrace{\frac{1}{2}{\parallel Ax-b\parallel }_2^2}}+\underset{g(z)}{\underbrace{{\lambda }_1{\parallel z\parallel }_1+\frac{1}{2}{\lambda }_2{\parallel z\parallel }_2^2}}s.t.x-z=0$$

(16)

We can minimize Equation (16) w.r.t. x and z in the same way as we minimized the ordered ${\ell }_2$ regularization in Section 4 and Section 4.1 and Section 4.2. Therefore, we directly present the final results below without any details. The + sign means to select max (0,value).

$$\begin{array}{l}{x}^{k+1}={(A^{T}A+\rho \ast I)}^{-1}(A^{T}b+\rho (z^k-u^k))\\ z^{k+1}={(\frac{\rho (x^{k+1}+u^k)-{\lambda }_1}{{\lambda }_2+\rho })}_{+}-{(\frac{-\rho (x^{k+1}+u^k)-{\lambda }_1}{{\lambda }_2+\rho })}_{+}\\ u^{k+1}=u^k+(x^k-z^k)\end{array}$$

(17)

5. Experiments

A series of experiments were conducted on both simulated and real data to examine the performance of a proposed method. In this section, first, a concept to select a correct sequence of $({\lambda }_k)$s is discussed. Second, an experiment on synthetic data is presented that describes the convergence of the lasso, SortedL1, ADMM-O${\ell }_2$, and ADMM-O${\ell }_{1,2}$. Finally, the proposed method is applied to a real feature selection dataset. The performance of the ADMM-O${\ell }_{1,2}$ method is analyzed in comparison with state-of-the-art methods: the lasso and SortedL1. These two methods are chosen for comparison because they are very similar to the ADMM-O${\ell }_{1,2}$ method except they use the regular lasso and the ordered lasso, respectively, while ADMM-O${\ell }_{1,2}$ model employs the ordered ${\ell }_{1,2}$ regularization with ADMM.

Experimental setting: The algorithms were implemented on Scala Spark™ with Scala code in both distributed and non-distributed versions. A distributed version of experiments was carried out in a cluster of virtual machines with four nodes: one master and three slaves. Each node has 10 GB of memory, 8 cores, CentOS release 6.2, and amd64: core-4.0-noarch. Apache spark™ 1.5.1 was deployed on it. The scholars also used IntelliJ IDEA 15 ULTIMATE as a Scala editor, interactive build tool sbt version 0.13.8, and Scala version 2.10.4. The standalone machine is a Lenovo desktop running Windows 7 Ultimate with an Intel™ Core™ i₃ Duo 3.20 GHz CPU and 4 GB of memory. The scholars used MATLAB™ version 8.2.0.701 running on a single machine to draw all figures. The source codes of the lasso, SortedL1, and ADMM-O${\ell }_{1,2}$ are available at References [32,33,34], respectively.

5.1. Adjusting the Regularizing Sequence $({\lambda }_k)$ for the Ordered Ridge Regression

Figure 1 was drawn using Algorithm (2), where p = 5000. As seen in Figure 1, when the value of a parameter ($q=0.4$) becomes larger, the sequence $({\lambda }_k)$ decreases, while $({\lambda }_k)$ increases for a small value of $q=0.055$. However, the goal is to obtain a non-increasing order of sequence $({\lambda }_k)$ by adjusting the value of q, which stimulates convergence. Here, adjusting means tuning the value of the parameter q using the BHq procedure to yield a suitable sequence $({\lambda }_k)$ such that it improves performance.

5.2. Experimental Results of Synthetic Data

In this section, numerical examples show the convergences of ADMM-O${\ell }_{1,2}$, ADMM-O${\ell }_2$, and other methods. A tiny, dense example of an ordered ${\ell }_2$ regularization is examined, where the feature matrix A has n = 1500 examples and p = 5000 features. Synthetic data is generated as follows: create a matrix A and choose $A_{i,j}$ using $\mathcal{N}(0,1)$ and then normalize columns of the matrix A to have the unit ${\ell }_2$ norm. $x^0\in {\mathbb{R}}^p$ is generated such that each sampled from $x^0\sim$$\mathcal{N}(0,0.02)$ is a Gaussian distribution. Label b is calculated as b = A*$x^0$ + v, where v ∼$\mathcal{N}(0,{10}^{-3}\ast I)$, which is the Gaussian noise. A penalty parameter $\rho$ = 1.0, an over-relaxed parameter $\alpha$ = 1.0, and termination tolerances ${ϵ}^{abs}\le {10}^{-4}$ and ${ϵ}^{rel}\le {10}^{-2}$ are used. Variables $u^0\in {\mathbb{R}}^p$ and $z^0\in {\mathbb{R}}^p$ are initialized to be zero. $\lambda \in {\mathbb{R}}^p$ is a non-increasing ordered sequence according to Section 5.1 and Algorithm 2. Figure 2a,b indicates the convergence of ADMM-O${\ell }_2$ and ADMM-O${\ell }_{1,2}$, respectively. Figure 3a,b shows the convergence of the ordered ${\ell }_1$ regularization and the lasso, respectively. It can be seen from the Figure 2 and Figure 3 that the ordered ${\ell }_2$ regularization converges faster than all algorithms. The ordered ${\ell }_1$, lasso, ordered ${\ell }_{1,2}$, and ordered ${\ell }_2$ take less than 80, 30, 30, and 10 iterations, respectively to converge. Dual is not guaranteed to be feasible. Therefore, a level of infeasibility of dual is also needed to compute. A numerical experiment terminates whenever both the infeasibility ($\widehat{w}$) and relative primal-dual gap ($\delta (b)$) are less than equal to $\lambda (1)$(Tolinfeas (${ϵ}^{infeas}={10}^{-6}$) and TolRelGap (${ϵ}^{gap}={10}^{-6}$), respectively) the ordered ${\ell }_1$ regularization harness synthetic data provided by Reference [5]. The same data is generated for the lasso as for the ordered ${\ell }_2$ regularization except for an initial value of $\lambda$. For the lasso, set $\lambda =0.1\ast {\lambda }_{max}$, where ${\lambda }_{max}={\parallel A^T\ast b\parallel }_{\infty }$. The researchers also use 10-fold cross-validation (cv) with the lasso. For further details about this step, see Reference [9].

5.3. Experimental Results of Real Data

Variable selection difficulty arises when the number of features (p) are greater than the number of instances (n). The proposed method genuinely handles these types of issues. The practical application of the ADMM method is in many domains such as computer vision and graphics [35], analysis of biological data [36,37], and smart electric grids [38,39]. A biological leukaemia dataset [40] was used to demonstrate the performance of the proposed method. Leukemia is a type of cancer that impairs the body’s ability to build healthy blood cells. Leukemia begins in the bone marrow. There are many types of leukemia, such as acute lymphoblastic leukemia, acute myeloid leukemia, and chronic lymphocytic leukemia. The following two types of leukemia are used in this experiment: acute lymphoblastic leukemia and acute myeloid leukemia. The leukemia dataset consists of 7129 genes and 72 samples [41]. Randomly split the data into training and test sets. In the training set, there are 38 samples, among which 27 are type I ALL (acute lymphoblastic leukemia) and 11 are type II AML (acute myeloid leukemia). The remaining 34 samples allowed us to test the prediction accuracy. The test set contains 20 type I ALL and 14 type II AML. The data were labeled according to the type of leukemia (ALL or AML). Therefore, before applying an ordered elastic net, the type of leukemia (ALL = −1, or AML = 1) is converted as a (−1, 1) response y. Predicted response $\widehat{y}$ is set to 1 if $\widehat{y}>0$; otherwise, it is set to −1. $\lambda \in {\mathbb{R}}^p$ is a non-increasing, ordered sequence generated according to Section 5.1 and Algorithm 2. For the regular lasso, $\lambda$ is a single scalar value generated using Equation (14). $\alpha =0.1$ is used for the leukemia dataset. All other settings are the same as experiment with synthetic data.

Table 3. Summary of variable selection in leukaemia dataset.

q	Method	Test Error	#Genes	Time	#Iter
0.1	Lasso	2.352941	6	5.459234 s	10,000
	The ordered ${\ell }_1$	2.235294	56	16.208356 s	10,000
	The ordered ${\ell }_2$	2.117647	All	0.176351 s	216
	The ordered ${\ell }_{1,2}$	2.352941	109	1.391347 s	2104
0.2	Lasso	2.352941	6	5.419032 s	10,000
	The ordered ${\ell }_1$	2.352941	65	14.990179 s	10,000
	The ordered ${\ell }_2$	2.117647	All	0.167623 s	197
	The ordered ${\ell }_{1,2}$	2.352941	107	1.046763 s	1597
0.3	Lasso	2.352941	6	5.394828 s	10,000
	The ordered ${\ell }_1$	2.352941	85	15.477436 s	10,000
	The ordered ${\ell }_2$	2.117647	All	0.140347 s	185
	The ordered ${\ell }_{1,2}$	2.117647	108	0.820148 s	1276
0.4	Lasso	2.235294	7	5.470428 s	10,000
	The ordered ${\ell }_1$	2.352941	90	12.206446 s	10,000
	The ordered ${\ell }_2$	2.117647	All	0.135451 s	178
	The ordered ${\ell }_{1,2}$	2.0	107	0.685255 s	1055
0.5	Lasso	2.0	8	5.387423 s	10,000
	The ordered ${\ell }_1$	2.352941	126	13.226632 s	10,000
	The ordered ${\ell }_2$	2.117647	All	0.126034 s	172
	The ordered ${\ell }_{1,2}$	2.0	102	0.603964 s	871
Average	Lasso	2.2588234	6.6	5.426189 s	10,000
	The ordered ${\ell }_1$	2.3294116	84.4	14.4218098 s	10,000
	The ordered ${\ell }_2$	2.117647	All	0.1491612 s	189.6
	The ordered ${\ell }_{1,2}$	2.1647058	106.6	0.9094954 s	1380.6

illustrates the experiment results of the leukemia dataset for different types of regularization. The lowest average mean square error (MSE) of regularization is the ordered ${\ell }_2$, followed by the ordered ${\ell }_{1,2}$ and lasso, while the highest average MSE can be seen in the ordered ${\ell }_1$. Looking at Table 3 first, it is clear that the ordered ${\ell }_2$ converges the fastest among all the regularizations. The second fastest converging regularization is the ordered ${\ell }_{1,2}$, while the slowest converging regularization is the ordered ${\ell }_1$. The ordered ${\ell }_2$ takes an average iteration around 190 and an average time around 0.15 s to converge. On the other hand, the ordered ${\ell }_{1,2}$, the ordered ${\ell }_1$, and lasso take average iterations around 1381, 10,000, and 10,000, respectively, and average times around: 1.0, 14.0, and 5.0 s, respectively, to converge. It can also be seen from the data in Table 3 that the ordered ${\ell }_2$ selected all the variables but that the goal is to select only the relevant variables from strongly correlated, high-dimensional dataset. Therefore, the ordered elastic net was proposed, which only selects relevant variables and discards irrelevant variables. As can be seen from Table 3, average MSE, time, and iteration in the ordered ${\ell }_1$ regularization and lasso are significantly more than the ordered ${\ell }_{1,2}$ regularization, although an average gene selection in the ordered ${\ell }_{1,2}$ regularization is more than that of the ordered ${\ell }_1$ regularization and lasso. The ordered ${\ell }_1$ and lasso select averages around 84 and 7 variables, respectively, whereas the ordered ${\ell }_{1,2}$ selects an average around 107 variables. The lasso performs poorly on the leukemia dataset. The reason for this is that strongly correlated variables are present in the leukemia dataset. In general, the ordered elastic net performs better than the ordered ${\ell }_1$ and lasso. Figure 4 shows the ordered elastic net solution paths and the variable selection results.

6. Conclusions

In this paper, we showed a method for optimizing an ordered ${\ell }_2$ problem under an ADMM framework, called ADMM-O${\ell }_2$. As an implementation of ADMM-O${\ell }_2$, the ridge regression with the ordered ${\ell }_2$ regularization is shown. We also presented a method for variable selection, called ADMM-O${\ell }_{1,2}$ which employs the ordered ${\ell }_1$ and ${\ell }_2$. We see the ordered ${\ell }_{1,2}$ as a generalization of the ordered ${\ell }_1$, which is shown as an important tool for model fitting, feature selection, and parameter estimation.

Experimental results show that the ADMM-O${\ell }_{1,2}$ method correctly estimates parameter, selects relevant variables, and excludes irrelevant variables for microarray data. Our method is also computationally tractable, adaptive, and distributed. The ordered ${\ell }_2$ regularization is convex and can be optimized efficiently with faster convergence rate. Additionally, we have shown that our algorithm has complexity $\mathcal{O}(p^3)$ + k * $\mathcal{O}(n{p}^2+p)$ heuristically, where k is the number of iterations. In future work, we plan to apply our method to other regularization models with complex penalties.