A unified consensus-based parallel ADMM algorithm for high-dimensional regression with combined regularizations

The advancement of modern science and technology has made data collection increasingly easy, resulting in the generation of massive amounts of data. However, due to the sheer volume of data and other factors such as privacy concerns, it has become necessary to store this data in a distributed manner. Therefore, it is essential to design parallel algorithms that can adapt to this massive and distributed storage dataset. Interested readers can refer to Boyd et al. (2010), Parikh and Boyd (2014), Lin et al. (2022), and their references for further information.

In this paper, our focus is on high-dimensional regression problems with massive data. When not considering distributed storage, the data required for a regression model typically includes the following,

Here, $y_{i}$ denotes the target value for the $i$-th observation $\bm{x}_{i}$, and $\bm{x}_{i}$ represents a $p$-dimensional vector that signifies the $i$-th observation value of $p$ features. These features could be the original observations and/or selected functions constructed from them. The objective of regression is to establish a relationship between the features in $\bm{X}$ and the target values in $\bm{y}$ by estimating coefficients that minimize the difference between the predicted and actual values of $\bm{y}$. To avoid overfitting and improve interpretability, some regularization terms are considered to be added to the objective optimization function, that is,

where $\mathcal{L}(\bm{y}-\bm{X}\bm{\beta})=\frac{1}{n}\sum_{i=1}^{n}{L}(y_{i}-\bm{x}_{i}^{\top}\bm{\beta})$. Here, the loss function ${L}$ is a generic function that includes least squares, Huber loss in Huber (1964) and quantile loss in Koenker and Basset (1970). We also consider situations where $\mathcal{L}(\bm{y}-\bm{X}\bm{\beta})$ cannot be separated, such as square root in Belloni et al. (2011). The regularization term $P_{\lambda}(|\bm{\beta}|)$ with tuning parameter $\lambda>0$ can be a single regularization term, such as $\ell_{2}$ or ridge in Hoerl and Kennard (1970), total variation in Rudin et al. (1992), lasso in (Tibshirani (1996)), group lasso in Yuan and Lin (2008); or a combination of them, such as sparse fused lasso (lasso plus total variation) in Tibshirani et al. (2005), elastic-net (lasso plus $\ell_{2}$) in Zou (2006), sparse group lasso (lasso plus group lasso) in Wu and Lange (2008); or their nonconvex variants. Due to the lack of unified algorithms for these combined regularization regressions, this paper focuses primarily on addressing this issue.

The alternating direction method of multipliers (ADMM) is an iterative optimization algorithm commonly used to solve convex minimization problems with linear constraint. It decomposes these problems into subproblems and can be applied to various statistical learning models. Boyd et al. (2010) provided a good overview of its development and detailed applications in various fields. Because some of the functions $L$ are nonsmooth and nondifferentiable, it is necessary to introduce the linear constraint $\bm{r}=\bm{y}-\bm{X}\bm{\beta}$ in order to use ADMM to uniformly solve (1.2). The constrained optimization problem is then formulated as follows,

Thus, the augmented Lagrangian form of (1) is

where $\bm{d}$ is dual variable corresponding to the linear constraint, and $\mu>0$ is a given penalty parameter. The iterative scheme of ADMM for (1) is

In general, when designing algorithms for distributed parallel processing, it requires a central machine and several local machines. Assume that

are stored in a distributed manner across $M$ local machines. The central machine receives information transmitted by the local machines, consolidates the information, and then forwards it to the respective local machines. The decomposition in (1.6) enables the algorithm to be naturally parallelizable across multiple local machines. Each local machine can independently solve its designated subproblem, and the solutions are then communicated and coordinated between the central machine and the local machines through variable updates. This distributed storage allows for parallel processing and efficient handling of large-scale datasets. The traditional ADMM algorithm in (1.5) can only handle problems on a single machine and is not sufficient for data stored in a distributed manner. Recently, many parallel ADMM algorithms have been proposed to solve statistical learning models when data is stored in a distributed manner, such as linear regression models (see Boyd et al. (2010), Yu et al. (2017) and Fan et al. (2020)). The existing parallel ADMM algorithms employ strategies similar to the regularized consensus problem described by Boyd et al. (2010). The consensus problem in distributed computing refers to the need for multiple nodes in a distributed system to agree on a common decision.

These existing parallel algorithms are specifically designed for individual regularization terms which can achieve sparse feature selection, such as lasso, group lasso, scad (Fan and Li (2001)) and mcp (Zhang (2010)). However, practical applications often involve other types of structural information beyond sparse feature structures. These include high correlation, feature grouping, or smooth trends among features. Combined regularization terms can simultaneously utilize different regularization techniques to integrate multiple structural information in the model, enhancing its performance and generalization ability. Thus, in such cases, many studies have suggested replacing a single regularization term with a combined regularization term. For example, elastic-net in Zou and Hastie (2005), fused lasso in Tibshirani et al. (2005), and sparse group lasso in Wu and Lange (2008). Moreover, many nonconvex variants of these combined regularization terms have been developed. For example, elastic-net variants such as snet (scad plus $\ell_{2}$) in Wang et al. (2010) and mnet (mcp plus $\ell_{2}$) in Huang et al. (2015); sparse fused lasso variants like sctv (scad plus total variation) and mctv (mcp plus total variation) in Xiu et al. (2019); sparse group lasso variants such as scgl (scad plus group lasso) and mcgl (mcp plus group lasso) in Huang et al. (2012) and Tugnait (2021b).

However, when the individual regularization is replaced by the combined regularization, the parallel algorithms mentioned above cannot be directly applied to these combined regularization regressions due to the complexity of the combination penalty terms, as well as the non-separability of certain composite regularization terms (such as sparse fused lasso and its nonconvex variants), which do not have closed-form solutions for their proximal operators. Currently, there are no parallel algorithms proposed for solving these high-dimensional problems with combined regularizations. The main work of this paper is to fill this gap and develop a unified consensus-based parallel ADMM (CPADMM) algorithm for high-dimensional regression problems with combined regularizations. The proposed algorithm has three main advantages: 1. It can not only handle massive amounts of data, but also adapt to the structure of distributed data storage. 2. Our algorithm is applicable to a wide range of convex loss functions, as long as we can derive the closed-form solutions for the proximal operators corresponding to the losses. 3. It has a guarantee of global convergence and a linear convergence rate.

The rest of the paper is organized as follows. Section 2 introduces some widely used combined regularization terms, and introduces some proximal operators that can play a crucial role in the algorithm. In Section 3, we first review the existing parallel ADMM algorithms for high-dimensional linear regression, and then propose a unified regression parallel ADMM algorithm for solving combined regularization. Section 4 discusses the specific implementation of the algorithm for solving various combined regularization regression. The convergence analysis is provided in Section 5. Section 6 presents numerical results demonstrating the scalability, efficiency and accuracy of the proposed algorithm. Section 7 summarizes the findings and concludes the paper, with a discussion on future research directions. Technical proofs, algorithm extensions and additional experiments are provided in the Appendix. The R package CPADMM for implementing the parallel ADMM algorithms can be obtained at https://github.com/xfwu1016/CPADMM.

Notations: $\bm{0}_{n}$ and $\bm{1}_{n}$ represent $n$-dimensional vectors with all elements being 0 and 1, respectively. $\bm{F}$ is a $(p-1)\times p$ matrix with all elements being 0, except for 1 on the diagonal and -1 on the superdiagonal. $\bm{I}_{n}$ represents the $n$-dimensional identity matrix. The Hadamard product is denoted by $\odot$. The sign$(\cdot)$ function is defined component-wise such that sign$(t)=1$ if $t>0$, sign$(t)=0$ if $t=0$, and sign$(t)=-1$ if $t<0$. $(\cdot)_{+}$ signifies the element-wise operation of extracting the positive component, while $|\cdot|$ denotes the element-wise absolute value function. For any vector $\bm{u}$, $\|\bm{u}\|_{1}$, $\|\bm{u}\|_{2}$, and $\|\bm{u}\|_{2,1}$ denote the $\ell_{1}$ norm, the $\ell_{2}$ norm, and the $\ell_{2,1}$ norm of $\bm{u}$, respectively. $\|\bm{u}\|_{\bm{H}}:=\sqrt{\bm{u}^{\top}\bm{H}\bm{u}}$ is used to denote the norm of $\bm{u}$ under the matrix $\bm{H}$, where $\bm{H}$ is a matrix.

In this section, we first review some commonly used combined regularization techniques and loss functions. Then we derive their corresponding proximal operators, as they play a crucial role in our algorithm.

Assuming that there are a total of $M$ local machines available, the response vector $\bm{y}\in\mathbb{R}^{n}$ and the sample observation matrix $\bm{X}\in\mathbb{R}^{n\times p}$ can be divided into $M$ blocks as follows,

where $\bm{y}_{m}\in\mathbb{R}^{n_{m}}$, $\bm{X}_{m}\in\mathbb{R}^{n_{m}\times p}$ and $\sum_{m=1}^{M}n_{m}=n$. In this section, we briefly review the existing consensus parallel ADMM algorithms for solving high-dimensional regression models. Based on existing algorithms, we propose a unified parallel ADMM algorithm for high-dimensional regression models with combined regularization.

From (3.11) and Algorithm 1, it can be observed that the operations on each local machine do not need to vary with changes in the regularization term. In other words, we only need to modify the operations on the central machine to solve high-dimensional regression models with different combined regularizers. Here, we first discuss the updates of local machines (in parallel), and then discuss the updates of the central machine separately according to three types of combined regularization terms. Note that we only need to focus on updating the primal variables, while the updates of the dual variables, $\bm{d}_{m}$, $\bm{e}_{m}$ and $\bm{f}$ involve linear algebraic operations and are easy to implement.

For the update of the subproblem $\bm{r}^{k+1}_{m}$ in (3.11), it is evidently an $\mathcal{L}$ proximal operator ($\text{prox}_{\frac{1}{n}\mathcal{L}}(\bm{y}_{m}+\bm{d}_{m}^{k}/\mu-\bm{X}_{m}\bm{\beta}_{m}^{k})$), and its closed-form solution can be obtained from Table 2.

For the update of the subproblem $\bm{\beta}^{k+1}_{m}$ in (3.11), the minimization problem is quadratic and differentiable, allowing us to solve the subproblem by solving the following linear equations:

When $p$ is large than $n_{m}$, Yu et al. (2017) suggested using the Woodbury matrix identity $(\bm{X}_{m}^{T}\bm{X}_{m}+\bm{I}_{p})^{-1}=\bm{I}_{p}-\bm{X}_{m}^{\top}(\bm{I}_{n_{m}}+\bm{X}_{m}\bm{X}_{m}^{\top})^{-1}\bm{X}_{m}$. This method is actually very practical when the size of $n_{m}$ is small because throughout the iteration process of the entire algorithm, the inverse only needs to be computed once. However, when both $n_{m}$ and $p$ are large, computing the inverse can be time-consuming, and in such cases, the conjugate gradient method performs better in solving the equation in (4.2). The conjugate gradient method is an efficient iterative algorithm used to solve symmetric positive definite linear systems of equations, and it has advantages in handling large-scale sparse problems.

Next, we will provide a detailed description of the update steps for the central machine for different combined regularization terms.