Parallel ADMM Algorithm with Gaussian Back Substitution for High-Dimensional Quantile Regression and Classification

The alternating direction method of multipliers (ADMM) is an iterative optimization method designed to solve complex convex minimization problems with linear constraints. It works by breaking down the original problem into smaller, more manageable subproblems that are easier to solve. ADMM alternates between optimizing these subproblems and updating dual variables to enforce the constraints, making it particularly useful for large-scale problems and various statistical learning applications (see Boyd et al. (2010) for more details). Since the quantile loss function $\rho_{\tau}(\cdot)$ is nonsmooth and nondifferentiable, it is necessary to introduce the linear constraint $\bm{r}=(r_{1},r_{2},\dots,r_{n})^{\top}=\tilde{\bm{y}}-\tilde{\bm{X}}\bm{\beta}$ to apply ADMM for solving problem (1.2). To better address the subproblem involving $\bm{\beta}$, one can introduce the equality constraint $\bm{z}=\bm{\beta}$. Consequently, the constrained optimization problem can be formulated as follows,

	$\displaystyle\min_{\bm{\beta},\bm{r},\bm{z}}$	$\displaystyle\quad\rho_{\tau}(\bm{r})+{P}_{\lambda}(|\bm{\beta}|),$
	s.t.	$\displaystyle\tilde{\bm{y}}-\tilde{\bm{X}}\bm{z}=\bm{r},\ \bm{z}=\bm{\beta},$		(2.1)

where $\rho_{\tau}(\bm{r})=\sum_{i=1}^{n}\rho(r_{i})$. The augmented Lagrangian form of (2.1) is

	$\displaystyle L_{\mu}(\bm{\beta},\bm{r},\bm{z},\bm{d}_{1},\bm{e}_{1})$	$\displaystyle=\rho_{\tau}(\bm{r})+{P}_{\lambda}(|\bm{\beta}|)-\bm{d}_{1}^{\top}(\tilde{\bm{y}}-\tilde{\bm{X}}\bm{z}-\bm{r})+\frac{\mu}{2}\|\tilde{\bm{y}}-\tilde{\bm{X}}\bm{z}-\bm{r}\|_{2}^{2}$
		$\displaystyle-\bm{e}_{1}^{\top}(\bm{z}-\bm{\beta})+\frac{\mu}{2}\|\bm{z}-\bm{\beta}\|_{2}^{2},$		(2.2)

where $\bm{d}_{1}$ and $\bm{e}_{1}$ are dual variables corresponding to the linear constraints, and $\mu>0$ is a given augmented parameter. Given $(\bm{r}^{0},\bm{z}^{0},\bm{d}_{1}^{0},\bm{e}_{1}^{0})$, the iterative scheme of ADMM for problem (2.1) is as follows,

$$\left\{\begin{array}[]{l}\bm{\beta}^{k+1}\ \leftarrow\mathop{\arg\min}\limits_{\bm{\beta}}\left\{L_{\mu}(\bm{\beta},\bm{r}^{k},\bm{z}^{k},\bm{d}_{1}^{k},\bm{e}_{1}^{k})\right\};\\ \bm{r}^{k+1}\ \leftarrow\mathop{\arg\min}\limits_{\bm{r}}\left\{L_{\mu}(\bm{\beta}^{k+1},\bm{r},\bm{z}^{k},\bm{d}_{1}^{k},\bm{e}_{1}^{k})\right\};\\ \bm{z}^{k+1}\ \leftarrow\mathop{\arg\min}\limits_{\bm{z}}\left\{L_{\mu}(\bm{\beta}^{k+1},\bm{r}^{k+1},\bm{z},\bm{d}_{1}^{k},\bm{e}_{1}^{k})\right\};\\ \bm{d}_{1}^{k+1}\ \leftarrow\bm{d}_{1}^{k}-\mu(\tilde{\bm{y}}-\tilde{\bm{X}}\bm{z}^{k+1}-\bm{r}^{k+1});\\ \bm{e}_{1}^{k+1}\ \leftarrow\bm{e}_{1}^{k}-\mu(\bm{z}^{k+1}-\bm{\beta}^{k+1}).\end{array}\right.$$

(2.3)

The entire iterative process of (2.3) is the non parallel QPADM algorithm proposed by Yu et al. (2017) and Wu et al. (2024). On the other hand, Gu et al. (2018) and Wu et al. (2025b) introduced the linearized ADMM algorithms for solving penalized quantile regression and SVM. These methods do not introduce an auxiliary variable $\bm{z}$, and instead linearize the quadratic function in the $\bm{\beta}$-subproblem to facilitate finding the solution for $\bm{\beta}$. However, their did not extend the algorithms to handle distributed storage data, meaning there is no parallel version of their ADMM algorithms. Therefore, we will focus on discussing the parallel QPADM algorithm from Yu et al. (2017) and Wu et al. (2025a).

When designing algorithms for distributed parallel processing, the setup generally involves a central machine and several local machines. Assume the data matrix $\tilde{\bm{X}}$ and the response vector $\tilde{\bm{y}}$ are distributed across $M$ local machines as follows,

$\displaystyle\tilde{\bm{X}}=(\tilde{\bm{X}}_{1}^{\top},\tilde{\bm{X}}_{2}^{\top},\dots,\tilde{\bm{X}}_{M}^{\top})^{\top}\ \text{and}\ \ \tilde{\bm{y}}=(\tilde{\bm{y}}_{1}^{\top},\tilde{\bm{y}}_{2}^{\top},\dots,\tilde{\bm{y}}_{M}^{\top})^{\top},$

(2.4)

where $\tilde{\bm{X}}_{m}\in\mathbb{R}^{n_{m}\times p}$, $\tilde{\bm{y}}_{m}\in\mathbb{R}^{n_{m}}$  and $\sum_{m=1}^{M}n_{m}=n$. To accommodate the structure of distributed storage data, Boyd et al. (2010) introduced $\bm{\beta}_{m}=\bm{\beta},m=1,2,\dots,M$, to enable parallel processing of the data. Then, the constrained optimization problem can be formulated as

	$\displaystyle\min_{\bm{\beta},\bm{r}_{m},\bm{\beta}_{m}}$	$\displaystyle\quad\sum_{m=1}^{M}\rho_{\tau}(\bm{r}_{m})+{P}_{\lambda}(|\bm{\beta}|),$
	$\displaystyle\text{s.t.}\ \bm{\beta}_{m}=\bm{\beta},\ \tilde{\bm{y}}_{m}-$	$\displaystyle\tilde{\bm{X}}_{m}\bm{\beta}_{m}=\bm{r}_{m},\ m=1,2,\dots,M,$		(2.5)

where $\rho(\bm{r}_{m})=\sum_{i=1}^{n_{m}}\rho(r_{i})$. The augmented Lagrangian form of (2.2) is

	$\displaystyle\small L_{\mu}(\bm{\beta},\bm{r}_{m},\bm{\beta}_{m},\bm{d}_{m},\bm{e}_{m})$	$\displaystyle=\ \sum_{m=1}^{M}\left[\rho(\bm{r}_{m})+\bm{d}_{m}^{\top}(\tilde{\bm{y}}_{m}-\tilde{\bm{X}}_{m}\bm{\beta}_{m}-\bm{r}_{m})+\frac{\mu}{2}\|\tilde{\bm{y}}_{m}-\tilde{\bm{X}}_{m}\bm{\beta}_{m}-\bm{r}_{m}\|_{2}^{2}\right.$
		$\displaystyle\left.+\ \bm{e}_{m}^{\top}(\bm{\beta}_{m}-\bm{\beta})+\frac{\mu}{2}\|\bm{\beta}_{m}-\bm{\beta}\|_{2}^{2}\right]+{P}_{\lambda}(|\bm{\beta}|),$		(2.6)

where $\bm{d}_{m}$ and $\bm{e}_{m}$ are dual variables corresponding to the linear constraints. Given $(\bm{r}_{m}^{0},\bm{\beta}_{m}^{0},\bm{d}_{m}^{0},\bm{e}_{m}^{0})$, the iterative scheme of parallel ADMM for problem (2.2) is as follows,

$$\left\{\begin{array}[]{l}\bm{\beta}^{k+1}\ \leftarrow\mathop{\arg\min}\limits_{\bm{\beta}}\left\{L_{\mu}(\bm{\beta},\bm{r}_{m}^{k},\bm{\beta}_{m}^{k},\bm{d}_{m}^{k},\bm{e}_{m}^{k})\right\};\\ \bm{r}_{m}^{k+1}\ \leftarrow\mathop{\arg\min}\limits_{\bm{r}_{m}}\left\{L_{\mu}(\bm{\beta}^{k+1},\bm{r}_{m},\bm{\beta}_{m}^{k},\bm{d}_{m}^{k},\bm{e}_{m}^{k})\right\};\\ \bm{\beta}_{m}^{k+1}\ \leftarrow\mathop{\arg\min}\limits_{\bm{\beta}_{m}}\left\{L_{\mu}(\bm{\beta}^{k+1},\bm{r}^{k+1},\bm{\beta}_{m},\bm{d}_{m}^{k},\bm{e}_{m}^{k})\right\};\\ \bm{d}_{m}^{k+1}\ \leftarrow\bm{d}_{m}^{k}-\mu(-\tilde{\bm{y}}_{m}+\tilde{\bm{X}}_{m}\bm{\beta}_{m}^{k+1}+\bm{r}_{m}^{k+1});\\ \bm{e}_{m}^{k+1}\ \leftarrow\bm{e}_{m}^{k}-\mu(-\bm{\beta}_{m}^{k+1}+\bm{\beta}^{k+1}).\end{array}\right.$$

(2.7)

This distribution enables efficient parallelization of the algorithm. Each local machine independently solves its own subproblem (e.g., $\bm{r}_{m},\bm{\beta}_{m},\bm{u}_{m},\bm{v}_{m}$), while the central machine consolidates the results by updating $\bm{\beta}$ and coordinates updates among the local machines. This framework facilitates parallel processing and effective management of large datasets. In subsection 3.2 of Yu et al. (2017), it is indicated that for convex penalties $P_{\lambda}$, the QPADM converges to the solution of (2.2). Here, we will briefly explain the proof approach. For ease of description, we will only discuss the case where $M=2$; the case where $M>2$ is similar.

When $M=2$, the constraint form of (2.2) is written in matrix form as

$$\bm{A}_{1}\bm{\beta}+\bm{A}_{2}\bm{r}_{1}+\bm{A}_{3}\bm{r}_{2}+\bm{B}_{1}\bm{\beta}_{1}+\bm{B}_{2}\bm{\beta}_{2}=\bm{e},$$

where

$$\bm{A}_{1}=[-\bm{I}_{p},-\bm{I}_{p},\bm{0}^{\top},\bm{0}^{\top}]^{\top},\quad\bm{A}_{2}=[\bm{0}^{\top},\bm{0}^{\top},\bm{I}_{n_{1}},\bm{0}^{\top}]^{\top},\quad\bm{A}_{3}=[\bm{0}^{\top},\bm{0}^{\top},\bm{0}^{\top},\bm{I}_{n_{2}}]^{\top},$$

$$\bm{B}_{1}=[\bm{I}_{p},\bm{0}^{\top},\bm{X}_{1}^{\top},\bm{0}^{\top}]^{\top},\quad\bm{B}_{2}=[\bm{0}^{\top},\bm{I}_{p},\bm{0}^{\top},\bm{X}_{2}^{\top}]^{\top},\quad\bm{e}=[\bm{0}^{\top},\bm{0}^{\top},\bm{y}_{1}^{\top},\bm{y}_{2}^{\top}]^{\top}.$$

Note that $\bm{A}_{1}$, $\bm{A}_{2}$, and $\bm{A}_{3}$ are mutually orthogonal, as are $\bm{B}_{1}$ and $\bm{B}_{2}$. Therefore, $\bm{\beta}$, $\bm{r}_{1}$, and $\bm{r}_{2}$ can be considered as a single variable, and $\bm{\beta}_{1}$ and $\bm{\beta}_{2}$ can also be considered as a single variable. Consequently, the QPADM degrades to a traditional two-block ADMM algorithm in the parallel case. A similar conclusion applies to non-parallel QPADM, provided that $\bm{\beta}$ and $\bm{r}$ are treated as a single block, and $\bm{z}$ is treated as a separate block. The convergence of the two-block traditional ADMM algorithm has been well-studied, thus QPADM is convergent as well. According to He and Yuan (2015), it exhibits linear convergence rate.

However, this does not apply to QPADM-slack in Guan et al. (2020) and Fan et al. (2021) because the introduction of slack variables prevents it from being formulated as a two-block ADMM. We will provide a detailed discussion on this in subsequent sections, which also serves as the motivation for this paper.

Simulation results from Yu et al. (2017) indicated that QPADM may require a large number of iterations to converge for penalized quantile regressions. This poses challenges for their practical application to big data, particularly in distributed environments where communication costs are high. To address this issue, Fan et al. (2021) proposed using two sets of nonnegative slack variables $\bm{\xi}=(\xi_{1},\xi_{2},\ldots,\xi_{n})^{\top}\in\mathbb{R}^{n}$ and $\bm{\eta}=(\eta_{1},\eta_{2},\ldots,\eta_{n})^{\top}\in\mathbb{R}^{n}$ to represent the quantile loss. Thus, (2.2) can be written as

		$\displaystyle\arg\min_{\bm{\beta},\bm{\xi}_{m},\bm{\eta}_{m},\bm{\beta}_{m}}\quad\sum_{m=1}^{M}\left[\tau\bm{1}_{n_{m}}^{\top}\bm{\xi}_{m}+(1-\tau)1_{n_{m}}^{\top}\bm{\eta}_{m}\right]+{P}_{\lambda}(|\bm{\beta}|),$
		$\displaystyle\text{s.t.}\ \bm{\beta}_{m}=\bm{\beta},\ \tilde{\bm{y}}_{m}-\tilde{\bm{X}}_{m}\bm{\beta}_{m}=\bm{\xi}_{m}-\bm{\eta}_{m},\bm{\xi}_{m}\geq 0,\bm{\eta}_{m}\geq 0,$		(3.1)

where $\bm{\xi}=(\bm{\xi}_{1}^{\top},\bm{\xi}_{2}^{\top}.\dots,\bm{\xi}_{M}^{\top})^{\top}$ and $\bm{\eta}=(\bm{\eta}_{1}^{\top},\bm{\eta}_{2}^{\top}.\dots,\bm{\eta}_{M}^{\top})^{\top}$. When $\tau=1$, the above equation is used by Guan et al. (2020) to complete the classification task.

The augmented Lagrangian form of (3.1) is

	$\displaystyle\small L_{\mu}(\bm{\beta},$	$\displaystyle\bm{\xi}_{m},\bm{\eta}_{m},\bm{\beta}_{m},\bm{d}_{m},\bm{e}_{m})=\sum_{m=1}^{M}\left[\tau\bm{1}_{n_{m}}^{\top}\bm{\xi}_{m}+(1-\tau)\bm{1}_{n_{m}}^{\top}\bm{\eta}_{m}+\bm{d}_{m}^{\top}(\bm{\beta}_{m}-\bm{\beta})+\frac{\mu}{2}\|\bm{\beta}_{m}-\bm{\beta}\|_{2}^{2}\right.$
		$\displaystyle\left.+\bm{e}_{m}^{\top}(\tilde{\bm{y}}_{m}-\tilde{\bm{X}}_{m}\bm{\beta}_{m}-\bm{\xi}_{m}+\bm{\eta}_{m})+\frac{\mu}{2}\|\tilde{\bm{y}}_{m}-\tilde{\bm{X}}_{m}\bm{\beta}_{m}-\bm{\xi}_{m}+\bm{\eta}_{m}\|_{2}^{2}\right]+{P}_{\lambda}(|\bm{\beta}|).$		(3.2)

Given $(\bm{\xi}_{m}^{0},\bm{\eta}_{m}^{0},\bm{\beta}_{m}^{0},\bm{d}_{m}^{0},\bm{e}_{m}^{0})$, the iterative scheme of parallel ADMM for problem (3.1) is as follows,

$$\left\{\begin{array}[]{l}\bm{\beta}^{k+1}\ \leftarrow\mathop{\arg\min}\limits_{\bm{\beta}}\left\{L_{\mu}(\bm{\beta},\bm{\xi}_{m}^{k},\bm{\eta}_{m}^{k},\bm{\beta}_{m}^{k},\bm{d}_{m}^{k},\bm{e}_{m}^{k})\right\};\\ \bm{\xi}_{m}^{k+1}\ \leftarrow\mathop{\arg\min}\limits_{\bm{\xi}_{m}\geq\bm{0}}\left\{L_{\mu}(\bm{\beta}^{k+1},\bm{\xi}_{m},\bm{\eta}_{m}^{k},\bm{\beta}_{m}^{k},\bm{d}_{m}^{k},\bm{e}_{m}^{k})\right\};\\ \bm{\eta}_{m}^{k+1}\ \leftarrow\mathop{\arg\min}\limits_{\bm{\eta}_{m}\geq\bm{0}}\left\{L_{\mu}(\bm{\beta}^{k+1},\bm{\xi}_{m}^{k+1},\bm{\eta}_{m},\bm{\beta}_{m}^{k},\bm{d}_{m}^{k},\bm{e}_{m}^{k})\right\};\\ \bm{\beta}_{m}^{k+1}\ \leftarrow\mathop{\arg\min}\limits_{\bm{\beta}_{m}}\left\{L_{\mu}(\bm{\beta}^{k+1},\bm{\xi}_{m}^{k+1},\bm{\eta}_{m}^{k+1},\bm{\beta}_{m},\bm{d}_{m}^{k},\bm{e}_{m}^{k})\right\};\\ \bm{d}_{m}^{k+1}\ \leftarrow\bm{d}_{m}^{k}-\mu(-\bm{\beta}_{m}^{k+1}+\bm{\beta}^{k+1});\\ \bm{e}_{m}^{k+1}\ \leftarrow\bm{e}_{m}^{k}-\mu(-\tilde{\bm{y}}_{m}+\tilde{\bm{X}}_{m}\bm{\beta}_{m}+\bm{\xi}_{m}^{k+1}-\bm{\eta}_{m}^{k+1}).\end{array}\right.$$

(3.3)

Clearly, $\bm{\beta}^{k+1}$ will be updated on the central machine, while $\bm{\xi}_{m}^{k+1},\bm{\eta}_{m}^{k+1},\bm{\beta}_{m}^{k+1},\bm{d}_{m}^{k+1}$ and $\bm{e}_{m}^{k+1}$ will be updated in parallel on the sub machines. In section 3 of Fan et al. (2021), they provided closed-form solutions for the ($\bm{\beta},\bm{\xi}_{m},\bm{\eta}_{m},\bm{\beta}_{m}$)-subproblems, which significantly facilitated the implementation of the parallel algorithm.

The iteration process described above cannot be reduced to a two-block ADMM algorithm; it can only be reduced to a three-block ADMM algorithm. Next, we will briefly explain this point using mathematical expressions. When $M\geq 2$, the constraint form of (3.1) is written in matrix form as

$\displaystyle\bm{A}_{1}\bm{\beta}+\sum_{m=2}^{M+1}\bm{A}_{m}\bm{\xi}_{m-1}+\sum_{m=1}^{M}\bm{B}_{m}\bm{\eta}_{m}+\sum_{m=1}^{M}\bm{C}_{m}\bm{\beta}_{m}=\bm{e},$

(3.4)

where $\bm{A}_{m}$, $\bm{B}_{m}$, and $\bm{C}_{m}$ are block matrices partitioned into $2M$ blocks by rows, and $\bm{e}$ is the corresponding column vector partitioned into $2M$ blocks. For example, for $M=2$, $\bm{A}_{1}=[-\bm{I}_{p},-\bm{I}_{p},\bm{0}^{\top},\bm{0}^{\top}]^{\top}$, $\bm{A}_{2}=[\bm{0}^{\top},\bm{0}^{\top},\bm{I}_{n_{1}},\bm{0}^{\top}]^{\top}$, $\bm{A}_{3}=[\bm{0}^{\top},\bm{0}^{\top},\bm{0}^{\top},\bm{I}_{n_{2}}]^{\top}$, $\bm{B}_{1}=[\bm{0}^{\top},\bm{0}^{\top},-\bm{I}_{n_{1}},\bm{0}^{\top}]^{\top}$, $\bm{B}_{2}=[\bm{0}^{\top},\bm{0}^{\top},\bm{0}^{\top},-\bm{I}_{n_{2}}]^{\top}$, $\bm{C}_{1}=[\bm{I}_{p},\bm{0}^{\top},\bm{X}_{1}^{\top},\bm{0}^{\top}]^{\top}$, $\bm{C}_{2}=[\bm{0}^{\top},\bm{I}_{p},\bm{0}^{\top},\bm{X}_{2}^{\top}]^{\top}$, and $\bm{e}=[\bm{0}^{\top},\bm{0}^{\top},\bm{y}_{1}^{\top},\bm{y}_{2}^{\top}]^{\top}$. The above six matrices cannot be divided into two sets of mutually orthogonal matrices like QPADM. Indeed, these six matrices can be partitioned into three mutually orthogonal groups in the following order: $\bm{A}=[\bm{A}_{1},\bm{A}_{2},\bm{A}_{3}],\bm{B}=[\bm{B}_{1},\bm{B}_{2}]$, and $\bm{C}=[\bm{C}_{1},\bm{C}_{2}]$. The same operation can also be implemented when $M>2$, that is,

$\displaystyle\bm{A}=[\bm{A}_{1},\bm{A}_{2},\dots,\bm{A}_{M+1}],\ \bm{B}=[\bm{B}_{1},\bm{B}_{2},\dots,\bm{B}_{M}],\ \text{and}\ \bm{C}=[\bm{C}_{1},\bm{C}_{2},\dots,\bm{C}_{M}].$

(3.5)

Hence, $\bm{\beta}$ and $\bm{\xi}_{m}$ can be considered as the first variable, $\bm{\eta}_{m}$ as the second variable, and $\bm{\beta}_{m}$ as the third variable. Thus, QPADM-slack is actually a three-block ADMM algorithm in terms of optimization form.

Chen et al. (2016) demonstrated that directly extending the ADMM algorithm for convex optimization with three or more separable blocks may not guarantee convergence, and they provided an example of divergence. As a result, QPADM-slack does not have a theoretical guarantee of convergence even in the convex case.

He et al. (2012) pointed out that although the ADMM algorithm with three blocks cannot be theoretically proven to converge, the algorithm can be corrected for some iterative solutions through simple Gaussian back substitution, thereby making the algorithm convergent. Because the correction matrix is an upper triangular block matrix, it is called Gaussian back substitution. The technique of Gaussian back substitution was extensively applied to various statistical and machine learning domains, including but not limited to the works of Ng et al. (2013), He and Yuan (2013), He et al. (2017), and He and Yuan (2018).

Let

$\displaystyle\bm{a}=(\bm{\beta},\bm{\xi}_{1},\dots,\bm{\xi}_{M}),\ \bm{b}=(\bm{\eta}_{1},\dots,\bm{\eta}_{M})\ \text{and}\ \bm{c}=(\bm{\beta}_{1},\dots,\bm{\beta}_{M}),$

(3.6)

and $\bm{d}=(\bm{d}_{1},\dots,\bm{d}_{M},\bm{e}_{1},\dots,\bm{e}_{M})$. To facilitate the description of the correction process (Gaussian back substitution), we define the sequence of $k+1$ iterations $({\bm{b}}^{k+1},{\bm{c}}^{k+1})$, generated by (3.3), as $(\tilde{\bm{b}}^{k},\tilde{\bm{c}}^{k})$. Then, the Gaussian back substitution of QPADM-slack is defined as

$\displaystyle\begin{bmatrix}\bm{b}^{k+1}\\ \bm{c}^{k+1}\end{bmatrix}=\begin{bmatrix}{\bm{b}}^{k}\\ {\bm{c}}^{k}\end{bmatrix}-\begin{bmatrix}v\bm{I}&-v(\bm{B}^{\top}\bm{B})^{-1}\bm{B}^{\top}\bm{C}\\ \bm{0}&v\bm{I}\end{bmatrix}\begin{bmatrix}{\bm{b}}^{k}-\tilde{\bm{b}}^{k}\\ {\bm{c}}^{k}-\tilde{\bm{c}}^{k}\end{bmatrix},$

(3.7)

where $\nu\in(0,1)$. Note that $\bm{B}^{\top}\bm{B}$ is an identity matrix, and thus

$$(\bm{B}^{\top}\bm{B})^{-1}\bm{B}^{\top}\bm{C}=\bm{B}^{\top}\bm{C}=\begin{bmatrix}-\bm{X}_{1}&0&\cdots&0\\ 0&-\bm{X}_{2}&\cdots&0\\ 0&0&\cdots&0\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\cdots&-\bm{X}_{M}\end{bmatrix},$$

it follows that

$$\begin{bmatrix}\bm{b}^{k+1}\\ \bm{c}^{k+1}\end{bmatrix}=\begin{bmatrix}\bm{\eta}_{m}^{k+1}\\ \bm{\beta}_{m}^{k+1}\end{bmatrix}=\begin{bmatrix}(1-\nu){\bm{\eta}}_{m}^{k}+\nu\tilde{\bm{\eta}}_{m}^{k}-\nu\bm{X}_{m}(\bm{\beta}_{m}^{k}-\tilde{\bm{\beta}}_{m}^{k})\\ (1-\nu){\bm{\beta}}_{m}^{k}+\nu\tilde{\bm{\beta}}_{m}^{k}.\end{bmatrix}.$$

(3.8)

We summarize the process of QPADM-slack with Gaussian back substitution (QPADM-slack(GB)) in Algorithm 1. Its parallel algorithm diagram is shown in Figure 1. Algorithm 1 differs from QPADM-slack only in the correction steps performed in each local machine, specifically steps 3 and 4 in Local machines. This correction process involves only linear operations, with the heaviest computational load being the matrix-vector multiplication in the $\bm{\eta}_{m}^{k+1}$ correction step. Next, we will establish the global convergence of QPADM-slack(GB) and its linear convergence rate. The proof of the following theorem is provided in Appendix A.

In recent years, nonconvex penalized quantile regression has been extensively studied, both theoretically and algorithmically, see Wang et al. (2012), Fan et al. (2014), Gu et al. (2018), Fan et al. (2021), Wang and He (2024). However, there has been little research on quantile loss SVMs with nonconvex regularization terms. In this subsection, we propose the nonconvex quantile loss SVMs, defined as

$$\begin{array}[]{l}\mathop{\arg\min}\limits_{\bm{\beta},{\xi}\geq\bm{0}}\sum\limits_{i=1}^{n}{{{\xi}_{i}}}+P_{\bm{\lambda}}(|\bm{\beta}|)\\ \text{s.t.}\;\;\sum\limits_{j=1}^{p}y_{i}x_{ij}\beta_{j}\geq 1-{1\over\tau}{{\xi}_{i}},\quad i=1,2,\cdots n,\\ \ \ \ \ \ \ \sum\limits_{j=1}^{p}y_{i}x_{ij}\beta_{j}\leq 1+{1\over(1-\tau)}{{\xi}_{i}},\ i=1,2,\cdots n.\end{array}$$

(3.10)

When $\tau=1$, the second inequality is always satisfied, and thus the aforementioned SVM simplifies to the nonconvex hinge SVM addressed in Zhang et al. (2016) and Guan et al. (2018). It is not difficult to see that (3.10) can be included in (1.2). In this paper, we mainly consider two popular nonconvex regularizations, SCAD penalty and MCP penalty.

In high-dimensional penalized linear regression and classification models, convex regularization terms like adaptive LASSO and elastic net ensure global optimality and computational efficiency, while non-convex regularization terms may provide better estimation and prediction performance but pose computational challenges due to the lack of global optimality. Fortunately, for quantile regression models and SVM models with nonconvex penalties, Fan et al. (2014) and Zhang et al. (2016) respectively pointed out that these problems can be uniformly solved by combining the local linear approximation (LLA, Zou and Li (2008)) method with an effective solution for the weighted $\ell_{1}$ penalized quantile regression and SVM.

LLA involves using the first-order Taylor expansion of the nonconvex regularization term to replace the original nonconvex regularization term, that is

$\displaystyle P_{\bm{\lambda}}(|\bm{\beta}|)\approx P_{\bm{\lambda}}(|\bm{\beta}^{l}|)+\nabla P_{\bm{\lambda}}(|\bm{\beta}^{l}|)^{\top}(|\bm{\beta}|-|\bm{\beta}^{l}|),\ \text{for}\ \bm{\beta}\approx\bm{\beta}^{l},$

(3.11)

where $\bm{\beta}^{l}$ is the solution from the last iteration, and $\nabla P_{\bm{\lambda}}(|\bm{\beta}^{l}|)=(\nabla P_{\bm{\lambda}_{1}}(|\bm{\beta}_{1}^{l}|),\nabla P_{\bm{\lambda}_{2}}(|\bm{\beta}_{2}^{l}|),\dots,\nabla P_{\bm{\lambda}_{p}}(|\bm{\beta}_{p}^{l}|))^{\top}$.

$\bullet$ For SCAD, we have

$\displaystyle\nabla P_{\bm{\lambda}_{j}}(|\bm{\beta}_{j}|)=\begin{cases}\bm{\lambda}_{j},&\text{{if }}|\bm{\beta}_{j}|\leq\bm{\lambda}_{j},\\ \frac{a\bm{\lambda}_{j}-|\bm{\beta}_{j}|}{a-1},&\text{{if }}\bm{\lambda}_{j}<|\bm{\beta}_{j}|<a\bm{\lambda}_{j},\\ 0,&\text{{if }}|\bm{\beta}_{j}|\geq a\bm{\lambda}_{j}.\\ \end{cases}$

(3.12)

$\bullet$ For MCP, we have

$\displaystyle\nabla P_{\bm{\lambda}_{j}}(|\bm{\beta}_{j}|)=\begin{cases}\bm{\lambda}_{j}-\frac{|\bm{\beta}_{j}|}{a},&\text{{if }}|\bm{\beta}_{j}|\leq a\bm{\lambda}_{j},\\ 0,&\text{{if }}|\bm{\beta}_{j}|>a\bm{\lambda}_{j}.\\ \end{cases}$

(3.13)

Then, the nonconvex penalized quantile regressions and SVMs can be written as

$\displaystyle\arg\min_{\bm{\beta}}\left\{\sum_{i=1}^{n}\rho_{\tau}\left(\tilde{y_{i}}-\tilde{\bm{x}}_{i}^{\top}\bm{\beta}\right)+P_{\bm{\lambda}}(|\bm{\beta}|)\right\}.$

(3.14)

By substituting equation (3.11) into equation (3.14), we can obtain the following optimized form in a weighted manner,

$\displaystyle\bm{\beta}^{l+1}=\mathop{\arg\min}\limits_{\bm{\beta}}\left\{\sum_{i=1}^{n}\rho_{\tau}\left(\tilde{y_{i}}-\tilde{\bm{x}}_{i}^{\top}\bm{\beta}\right)+\sum_{j=1}^{p}\nabla P_{\bm{\lambda}_{j}}(|\beta_{j}^{l}|)|\beta_{j}|\right\}.$

(3.15)

Note that we only need to make a small change to solve this weighted optimization form using Algorithm 1. This change only requires replacing $\bm{\lambda}$ with $\nabla P_{\bm{\lambda}}(|\bm{\beta}^{l}|)$.

To address nonconvex penalized quantile regressions and SVMs through the LLA algorithm, it is crucial to identify a suitable initial value. Following the guidance of Zhang et al. (2016) and Gu et al. (2018), we can utilize the solution derived from $P_{\bm{\lambda}}(|\bm{\beta}|)=\|\bm{\lambda}\odot\bm{\beta}\|_{1}$ in (3.14) as the initial starting point. Subsequently, the solution to (3.14) is obtained by iteratively solving a series of weighted $\ell_{1}$ penalized quantile regressions and SVMs. The comprehensive iterative steps of this approach are systematically outlined in Algorithm 2.

The preceding discussion reveals that resolving nonconvex penalized quantile regressions and SVMs requires multiple iterations of weighted $\ell_{1}$ penalized models. Theoretically, as shown by Fan et al. (2014) and Zhang et al. (2016), only two to three iterations suffice to find the Oracle solution for (3.14). In implementing Algorithm 2, we utilize the warm-start technique (Friedman et al. (2010)), initializing $\bm{\beta}^{l+1}$ with $\bm{\beta}^{l}$, significantly reducing step 2.2’s iteration count, often leading to convergence within just a few to a dozen steps.

To address this issue, we suggest changing the variable update order of QPADM-slack in (4.1) to

$\displaystyle\cdots\to\bm{\beta}_{m}\to\bm{\xi}_{m}\to\bm{\eta}_{m}\to\bm{\beta}\to\cdots.$

(4.2)

We will explain the reason for adopting this order below. From the discussion in Section 3.2, it is evident that the first block among the three blocks of primal variables to be updated does not require correction. To avoid matrix operations in the correction step, $\bm{\beta_{m}}$ should be placed in the first block because only the constraint matrix corresponding to $\bm{\beta_{m}}$ is not composed of identity and zero matrices. During the Gaussian back substitution process, we need to compute the inverse of $\bm{B}^{\top}\bm{B}$. To keep this inverse as simple as possible, we must place either $\xi_{m}$ or $\eta_{m}$ in the second block. Here, we choose to place $\xi_{m}$ in the second block in sequence. The remaining $\eta_{m}$, $\bm{\beta}$ form the third block.

Thus, with $(\bm{\xi}_{m}^{0},\bm{\eta}_{m}^{0},\bm{\beta}^{0},\bm{d}_{m}^{0},\bm{e}_{m}^{0})$, the iterative scheme of modified QPADM-slack (M-QPADM-slack) for problem (3.1) is as follows,

$$\left\{\begin{aligned} \bm{\beta}_{m}^{k+1}&\leftarrow\mathop{\arg\min}\limits_{\bm{\beta}_{m}}\left\{L_{\mu}(\bm{\beta}_{m},\bm{\xi}_{m}^{k},\bm{\eta}_{m}^{k},\bm{\beta}^{k},\bm{d}_{m}^{k},\bm{e}_{m}^{k})\right\};\\ \tilde{\bm{\xi}}_{m}^{k}&\leftarrow\mathop{\arg\min}\limits_{\bm{\xi}_{m}\geq\bm{0}}\left\{L_{\mu}(\bm{\beta}_{m}^{k+1},\bm{\xi}_{m},\bm{\eta}_{m}^{k},\bm{\beta}^{k},\bm{d}_{m}^{k},\bm{e}_{m}^{k})\right\};\\ \tilde{\bm{\eta}}_{m}^{k}&\leftarrow\mathop{\arg\min}\limits_{\bm{\eta}_{m}\geq\bm{0}}\left\{L_{\mu}(\bm{\beta}_{m}^{k+1},\tilde{\bm{\xi}}_{m}^{k},\bm{\eta}_{m},\bm{\beta}^{k},\bm{d}_{m}^{k},\bm{e}_{m}^{k})\right\};\\ \tilde{\bm{\beta}}^{k}&\leftarrow\mathop{\arg\min}\limits_{\bm{\beta}}\left\{L_{\mu}(\bm{\beta}_{m}^{k+1},\tilde{\bm{\xi}}_{m}^{k},\tilde{\bm{\eta}}_{m}^{k},\bm{\beta},\bm{d}_{m}^{k},\bm{e}_{m}^{k})\right\};\\ \bm{d}_{m}^{k+1}&\leftarrow\bm{d}_{m}^{k}-\mu\left(-\bm{\beta}_{m}^{k+1}+\tilde{\bm{\beta}}^{k}\right);\\ \bm{e}_{m}^{k+1}&\leftarrow\bm{e}_{m}^{k}-\mu\left(-\tilde{\bm{y}}_{m}+\tilde{\bm{X}}_{m}\bm{\beta}_{m}^{k+1}+\tilde{\bm{\xi}}_{m}^{k}-\tilde{\bm{\eta}}_{m}^{k}\right);\end{aligned}\right.$$

(4.3)

where $L_{\mu}$ is defined as in (3.1). $\bm{\beta}_{m}^{k+1}$, $\tilde{\bm{\xi}}_{m}^{k}$, $\tilde{\bm{\eta}}_{m}^{k}$, $\bm{d}_{m}^{k+1}$ and $\bm{e}_{m}^{k+1}$ will be updated in parallel by each sub machine, while $\tilde{\bm{\beta}}^{k}$ will be updated on the central machine. However, this update sequence will cause a new problem, that is, the updates of local sub machines become incoherent, necessitating that the central machine first completes the update of $\tilde{\bm{\beta}}^{k}$ before proceeding with the update of $\bm{d}_{m}^{k+1}$. This will result in an additional round of communication between the central machine and the local machines, thereby reducing the efficiency of the algorithm. This issue will not occur in (3.3) because the updates of variables on its local sub machines are coherent. Fortunately, since the updates of each sub problem in (4.3) do not strictly depend on the previously updated variables, this issue can be resolved by adjusting the update positions of the dual variables $\bm{d}_{m}$ and $\bm{e}_{m}$.

For the convenience of describing the correction steps, we need to define

$\displaystyle\bm{a}=(\bm{\beta}_{1},\dots,\bm{\beta}_{M}),\ \bm{b}=(\bm{\xi}_{1},\dots,\bm{\xi}_{M})\ \text{and}\ \bm{c}=(\bm{\eta}_{1},\dots,\bm{\eta}_{M},\bm{\beta}).$

(4.4)

Correspondingly,

$$\bm{A}=\begin{bmatrix}\bm{I}_{p}&&\\ &\ddots&\\ &&\bm{I}_{p}\\ \bm{X}_{1}&&\\ &\ddots&\\ &&\bm{X}_{M}\end{bmatrix},\bm{B}=\begin{bmatrix}\bm{0}&&\\ &\ddots&\\ &&\bm{0}\\ -\bm{I}_{n_{1}}&&\\ &\ddots&\\ &&-\bm{I}_{n_{M}}\end{bmatrix}\text{and}\ \bm{C}=\begin{bmatrix}&&&-\bm{I}_{p}\\ &&&\vdots\\ &&&-\bm{I}_{p}\\ \bm{I}_{n_{1}}&&&\\ &\ddots&&\\ &&\bm{I}_{n_{M}}&\end{bmatrix}.$$

It is obvious that the matrices $\bm{A}$ and $\bm{B}$ are $2M\times M$ partitioned, while $\bm{C}$ is $2M\times(M+1)$ partitioned. $\bm{A}_{m}$, $\bm{B}_{m}$ and $\bm{C}_{m}$ correspond to the block matrices of the $m$-th columns of $\bm{A}$, $\bm{B}$ and $\bm{C}$, respectively. We also define the sequence of $k+1$ iterations $({\bm{b}}^{k+1},{\bm{c}}^{k+1})$, generated by (4.3), as $(\tilde{\bm{b}}^{k},\tilde{\bm{c}}^{k})$. According to (3.7), it follows from

$$(\bm{B}^{\top}\bm{B})^{-1}\bm{B}^{\top}\bm{C}=\begin{bmatrix}-\bm{I}_{n_{1}}&&&&\bm{0}\\ &-\bm{I}_{n_{2}}&&&\bm{0}\\ &&\ddots&&\vdots\\ &&&-\bm{I}_{n_{M}}&\bm{0}\end{bmatrix}$$

that