Enhanced Cyclic Coordinate Descent Methods for Elastic Net Penalized Linear Models

Reviewer 3

Thank you for your time and detailed comments. We provide responses to weaknesses and questions below.

W1.

We thank the reviewer for suggesting large-scale benchmarks. To ensure a fair comparison, we ran both glmnet and ECCD with identical convergence criteria and path settings:

• Tolerance: tol = 1e-7.
• Lambda path: lambda.min.ratio = if n < p then 0.1 else 1e-4, nlambda = 100; 0.1 used as $p$ is extremely large in this response's big dataset experiment.
• Standardization: standardize = TRUE.
• Logistic Hessian: type.logistic = "newton" (exact Hessian) and type.logistic = "modified-newton" (diagonal approximation) both tested, with intercept = TRUE, as GLMNet suggests the latter can be faster.

Synthetic Datasets Generation Process

The synthetic datasets are generated through the following steps. We begin by sampling a design matrix $X \in \mathbb{R}^{n \times p}$ from a multivariate normal distribution. Next, we impose sparsity on the true coefficient vector $\boldsymbol{\beta}$ with the intercept $\beta_0$, so the dataset fit in with elastic net problems. Then, the linear predictor is computed by adding Gaussian noise $\boldsymbol{\varepsilon}$ drawn from a multivariate normal distribution:

Finally, each $\eta_i$ is transformed using the logistic link function: $p_i = \frac{1}{1 + e^{-\eta_i}},$ and the binary response variables are sampled as $y_i \sim \mathrm{Bernoulli}(p_i)$.

Using these settings, we ran wall-clock timings on an AMD EPYC 7352 CPU across three synthetic regimes where the number of observations is less than, greater than, or equal to the number of features. The results (reported in the main text) demonstrate ECCD's consistent speedups even when GLMNet runs for minutes or longer.

Interpretation:

• Feature-dominated ($n<p$): ECCD cuts runtime by $\sim\\!1.8\text{–}2.2\times$ versus glmnet's Newton solver and show competitive, even superior performance for majority cases, versus glmnet's modified.Newton solver, even when $p$ reaches millions (e.g. $5\,000\times1\,000\,000$).

• Sample-dominated ($n>p$): Speedups range from $\sim \\!1.5\times$ on small-$p$ problems (e.g. $10\,000\times50$) up to $\sim \\!1.8\times$ versus Newton on moderate-$p$ ($10\,000\times5\,000$), and as much as $\sim \\!7\times$ versus modified.Newton at that scale.

• Balanced ($n=p$): ECCD consistently achieves $\sim \\!3\text{–}4\times$ faster fits than glmnet, e.g. cutting the 5 000$\times$5 000 case from 57.2 s down to 13.8 s, and still delivering a 1.5-1.8$\times$ speedup at $70\,000\times 70\,000$.

ECCD scales gracefully across all GLM regimes—feature-dominated, sample-dominated, and balanced—outperforming both glmnet's default cyclic solver and a Modified-Newton variant on very large problems (up to $70{,}000\times70{,}000$ and $1{,}000{,}000\times5{,}000$), shaving off tens to hundreds of seconds on a single thread.

With ECCD, the $1{,}000{,}000\times5{,}000$ problem runtime drops from 2 633.4 s ($\approx$43.9 min) to 1 632.2 s ($\approx$27.2 min), saving 1 001.2 s ($\approx$16.7 min). This dramatic reduction directly addresses real-world GLM runtime concerns. A shared-memory parallel version of ECCD—currently in development—will drive even greater speedups and make even larger datasets tractable.

We will include these results in the revised manuscript to highlight ECCD's impact on large-scale GLM problems.

W2. We agree that selecting $s$ requires tuning based on system-specific (i.e. hardware parameters) and data-specific (i.e. condition number, sparsity, dimensions) properties. In the paper, we will add a clear disclaimer that practical tuning is necessary. Theorem 3 suggests an asymptotically optimal $s$, however, it is possible to utilize a more detailed performance model which incorporate hardware/data properties (e.g. memory bandwidth/latency, cache size, processor speed, sparsity, dimensions) and obtain a better estimate for $s$. We will incorporate such a performance model in our planned future work on shared-memory parallelization of ECCD.

W3. We acknowledge that the baselines use different convergence criteria as stated by the reviewer. We chose to align our criteria with glmnet as minimizing the objective is our primary goal. For other baselines, we follow their default settings. While the stopping rules differ, they are closely related; under a second-order Taylor approximation, objective decrease is roughly proportional to the squared coefficient change. Thus, the comparison remains fair and meaningful. Our main goal with this paper is to illustrate that a numerically-stable block coordinate update is feasible for GLM problems through algorithmic redesign. We empirically show performance improvements over our target baseline, but also other packages (in these discussions and in the main paper).

W4. We recognize that public code release is essential for reproducibility and community adoption. We will:

• Open-source ECCD's implementation (R and C++), along with all scripts needed to reproduce the experiments of the main paper, appendices, and reviewer discussions.
• Hosting: The code repository will be hosted on GitHub with detailed installation and benchmarking instructions.

This will ensure that users can easily deploy ECCD and verify our results.

Thank you for your time and detailed comments. We provide responses to weaknesses and questions below. We provide benchmarks on representative datasets (real and synthetic), where the data generation process for synthetic data is described in the response to Reviewer NJYD.

We also appreciate the reviewer for pointing out the references, we will add the below citations to the paper.

References

• [1] J. Qian et al., skglm: Scikit-Learn compatible sparse-GLM solver, 2023.
• [2] Y. Wang et al., abess: Best-subset Selection Library, 2023.
• [3] Breheny & Huang., Group descent algorithms for nonconvex penalized linear and logistic regression models with grouped predictors, 2015.

Weaknesses

W1. We thank the reviewer for pointing out skglm and abess. Below, we explain why direct comparison is limited with empirical experiment results.

1. skglm

   • Implements pathwise coordinate-descent for linear regression in Python.
   • Does not support pathwise GLMs (logistic, Poisson, etc.), so it is not directly comparable to ECCD in those settings.
   • We benchmarked ECCD against the single-point logistic regression algorithm from skglm, evaluating both methods along the penalty path derived from ECCD on the Duke dataset, with default error tolerance matched to that used by ECCD. When $s=4$, ECCD ran in 0.013 s versus skglm's 1.375 s.
   • More results are presented in Table 1.

2. abess

   • The best-subset selection problem is characterized by an explicit $\ell_0$ constraint, which directly limits the number of nonzero coefficients. fundamentally different from the $\ell_1$ or $\ell_2$ regularized problems (Lasso, Elastic Net) we study.
   • abess's claimed advantage arises in high-correlation regimes. We tested under $\rho=0.9$ for both small $(n=20,p=2000)$ and moderate $(n=100,p=1000)$ settings and found abess is still slower than ECCD.
   • Under the exact setting, we found that Abess perform poorly on relatively large dataset (i.e. n = 1000)
   • On larger-scale datasets (e.g., $n=100000$), Abess's performance degrades substantially, yielding slower runtimes than ECCD under identical conditions.

Table 1: ECCD vs. skglm (time in seconds)

Table 2: ECCD vs. glmnet vs. abess (time in seconds)

NOTE: we attempt to run n = 100, p = 100000 with rho = 0.1, 0.9 on cluster, but the data generation process took over 12 hours.

W2. Thank you for pointing this out. Equation (14) in Theorem 1 is indeed a simplified bound derived from a Taylor expansion. It is primarily intended to offer theoretical intuition rather than a tight bound in all scenario. Nevertheless, our algorithm is designed to ensure that the relative approximation error remains small in practice. Specifically, we adopt a convergence criterion based on the relative change in the objective function, which is identical to the stopping rule used in glmnet. Regardless of the block size ($s$), the algorithm proceeds until the objective value achieves the specified threshold. This design guarantees that the final objective value remains consistently low for both small and large $s$, which accounts for the stable performance across block sizes observed in our experiments.

W3. We thank the reviewer for this observation and will modify our claims in the abstract and paper body to reflect the more realistic speedups of $1 - 3\times$.

W4.1. To investigate, we ran BigLasso on one node of the Cray EX Cluster used in the main paper. We experimented with using 2, 4, 8, 16, 32, 64, and 128 threads, and compared against our single-thread ECCD. Across a range of datasets and both full and short $\lambda$-paths, ECCD's active-set screening and path-wise updates yielded consistent performance advantage even versus BigLasso's multi-threaded solver. We report representative results below.

Table 3: Runtime on glmnet's default path ($\lambda_{\rm max} \to \lambda_{\rm max} \times 10^{-4}$ or $10^{-2}$, 100 points)

Table 4: Runtime on BigLasso's default path ($\lambda_{\rm max} \to \lambda_{\rm max} \times 10^{-1}$, 100 points)

Key observations:

1. Sequential ECCD outperforms BigLasso even when it is parallelized using 2–64 threads.
2. glmnet's default path incurs more work, yet ECCD's optimized screening keeps it faster.
3. BigLasso's default path still favor ECCD for moderate block sizes ($s\le16$).

This confirms that ECCD's pathwise, active-set strategy provides a fundamental performance advantage. We plan to address parallel algorithms design for ECCD in shared- and distributed-memory environments in future work. Once parallelized, we believe that ECCD will continue to outperform BigLASSO.

W4.2. We include additional experiments when $n$ and $p$ are large (e.g. $10^4$) below, comparing three solvers (glmnet, glmnet with Modified Newton, and ECCD) on total runtime:

As $n=p$ increases from 5 000 to 70 000:

• glmnet grows by roughly $30.79\times$ (15.0 s $\to$ 461.6 s),
• glmnet w. Modified Newton grows by roughly $6.73\times$ (57.2 s $\to$ 384.8 s), and
• ECCD grows by roughly $18.7\times$ (13.8 s $\to$ 258.3 s).

Moreover, we observed the relative speedup:

• At $n=p=40,000$, ECCD (69.4 s) is $\sim1.87\times$ faster than glmnet (129.9 s) and outperform modified case (195.2 s).
• Even at the largest size, $n=p=70\,000$, ECCD (258.3 s) maintains a $\sim1.5\times$ speedup over glmnet (384.8 s) and a $\sim1.8\times$ advantage over Modified Newton (461.6 s).

These results show that ECCD's active-set screening and block-coordinate updates scale much more efficiently in large, dimension-balanced problems.

All scaling experiments reported here were conducted on an AMD EPYC 7352 CPU using a single thread, as noted in our response to Reviewer NJYD's W1.

Table 5: Synthetic $n = p$ Benchmarks on Cluster

W5.

ECCD's approach: ECCD does not rely on a diagonal approximation of the Hessian. As shown in Equation (11) of the main paper, we incorporate all entries of the true Hessian to form a second‑order approximation of the loss. This full‑matrix update more faithfully captures curvature, yielding more accurate coordinate steps and leads to faster convergence overall.

glmnet's options: The glmnet package does offer a diagonal Hessian variant (e.g. via type.logistic = "modified-newton" for logistic regression)—which is claimed in the documentation to be a more efficient method. However, we evaluated this variant on large-scale datasets, and our experiments (see Table 5 above) show that when both $n$ and $p$ are large (e.g. $10^4$), this diagonal variant suffers noticeable degradation in both runtime and solution quality. By contrast, ECCD's full Hessian updates remain robust and efficient.

W6. Thanks for point this out, we have corrected this.

Questions

Q1. Yes, ECCD can be naturally extended to other penalties that are compatible with coordinate descent based optimization. And prior work by Breheny and Huang (2015) [3] has shown that both SCAD and MCP can be efficiently optimized using coordinate descent frameworks.

While these penalties primarily address nonconvex regularization, our method is designed to handle the nonlinearity in the loss function. ECCD focuses on accelerating coordinate descent by accurately approximating the nonlinear loss via second-order updates, and is agnostic to the specific form of the penalty.

Therefore, as long as the penalty admits coordinate-wise or block-wise updates, our ECCD framework can be applied. In our work, we adopt the $\ell_1$ penalty for simplicity and clarity of exposition, but the method naturally generalizes to other penalties.

Thank you for your time and detailed comments. We provide responses to weaknesses and questions below. We provide benchmarks on representative datasets (real and synthetic), where the data generation process for synthetic data is described in the response to Reviewer NJYD.

Weakness

W1. Theorem 3 (Appendix A.5) provides space complexity analysis which includes low order terms, including the additional memory required by ECCD beyond the storage cost of the design matrix. Specifically, ECCD allocates $\mathcal{O}(ns)$ memory for the working buffer of the active block $X_s$, and $\mathcal{O}(s^2)$ for the corresponding Gram matrix. These bounds can be specialized to glmnet by setting $s = 1$. The additional memory overhead of ECCD remains negligible in high-dimensional settings, with relative cost $\frac{ns}{np} \to 0$ as $p \gg s$. This is also empirically supported, as discussed in our response to the second question.

Importantly, the ideal value of $s$ is a small constant that is approximately independent of $n$ and $p$, and reflects the per-coordinate cost of evaluating the link function, as characterized in Theorem 2. Thus, assuming $p \gg s$ is both theoretically justified and reasonable in practice. If desirable, we can incorporate Theorem 3 from Appendix A.5 into the main text.

W2. To address this, we have rerun our benchmarks and now include peak RSS (resident set size) measurements for ECCD as the block size $s$ varies. These experiments serve to empirically validate our theoretical space-complexity analysis in Theorem 3 (Appendix A.5).

In our analysis, the extra storage required by ECCD consists of $s$ vectors of length $n$ (the working block $X_s$) plus the $s\times s$ Gram matrix, $X_s^\intercal X_s$. In extreme case where $s \approx p$, the Gram matrix storage cost ($O(s^2)$ term) could dominate; however, since we incorporate the active-set method similar to glmnet, this cost is smaller than storing $X$ and the active block $X_s$, itself. Our approach for $s$ is to bound it to $\sqrt{|\mathcal{A}|}$, where $\mathcal{A}$ is the active set in the current iteration. This ensures that further increases in $s$ do not increase memory usage. The memory overhead of ECCD remains negligible compared to the cost of storing the original data matrix.

For lambda path, by reusing the previous lambda solution as a warm start along the path, the RSS memory grows only until $s$ reaches square root of the largest active set in the lambda path.

Table 1 shows Memory vs. $s$ experiments for a synthetic $n=p=5000$ matrix. We instrumented our C++ code to obtain memory stats via getrusage(). The baseline RSS is $754.7$ MiB for $s = 1$ (i.e. equivalent to glmnet); At the largest value of $s=4096$, peak RSS increases by $41.2$ MiB when compared to $s = 1$. These results demonstrate that ECCD's memory footprint increases slowly with the primary cost being the original design matrix storage, even for large block sizes. If $s$ is allowed to exceed the active-set upper bound, the $O(s^2)$ term can dominate and increase memory usage substantially. We show experiments without active set to illustrate the memory trade off.

Table 1: Memory vs. block size ($n=p=5000$, with active-set constraint on $s$, baseline RSS = 754.7 MiB)

Table 1 summarize the case where the block is capped at the squared root of active set size $|\mathcal{A}|$, so memory usage rises only modestly with $s$ and plateaus once $s$ exceeds the $\sqrt{|\mathcal{A}|}$ upperbound when $n$ is moderately significant.

Table 2: Memory vs. block size ($n=p=5000$, without active-set constraint on $s$, baseline RSS = 809.6 MiB)

In Table 2, we removed the active set cap $|\mathcal{A}|$. Without the active set cap $|\mathcal{A}|$, the $\mathcal{O}(s^2)$ component dominates for large $s$, leading to a substantial rise in peak memory as block size increases. At $s=4096$ peak RSS is 1316.5 MiB, compared to 1000.4 MiB at $s = 1$ (an increase of $316.1$ MiB). Without active set, the storage cost of $X_s$ and the Gram matrix $X_s^\intercal X_s$ are not negligible. However, this is an impractical setting for performance reasons as glmnet incorporates active set to improve performance/efficiency. Similarly, the fastest variant of ECCD is one where active set is included.

Table 3: glmnet peak memory

In Table 3 shows peak RSS memory for glmnet using peakRAM R package. Because glmnet is implemented in optimized C++, we did not directly instrument the C++ code similar to ECCD, as this would require source code changes and building glmnet from the instrumented source code. We provide Table 3 to give a general sense in memory usage of glmnet (without instrumented C++ code). ECCD with $s = 1$ is equivalent to glmnet, so we expect the actual peak RSS memory to be similar to $s = 1$. Overall, these results confirm that ECCD's practical memory overhead is negligible—even for large $s$—with the largest memory cost arising from storing the original data matrix.

We provide additional RSS memory usage below under different settings:

Table 4: Memory vs. block size ($n=100$, $p=100\,000$; with active-set constraint on $s$, baseline RSS = 601.9 MiB)

Table 4 shows RSS memory for $p \gg n$ with active set.

Table 5: Memory vs. block size ($n = 100000$, $p = 100$, with active-set constraint on $s$, baseline RSS = 672.3 MiB)

Table 5 shows RSS memory for $n \gg p$ with active set.

Table 6: Memory vs. block size (Duke dataset, with active-set constraint on $s$, base = 510.5 MiB)

Table 6, shows RSS memory for the Duke dataset with active set. Tables 4 - 6 illustrate that memory usage for large $s$ is not significant when using the active set method.

W3. This error has been fixed.

Questions:

Q1. Please see our responses under Weakness 1 & 2 which includes a discussion of other memory costs (low order terms) from Theorem 3 Appendix A.5 and additional experiments to support our analysis.

Q2. We appreciate the reviewer's suggestion to evaluate larger block sizes. In response, we conducted additional experiments on the Duke dataset with block sizes up to $s = p = 7129$; the results are summarized in Table 8 below.

ECCD remains stable and convergent even when $s = p$, underscoring the robustness of our second‑order Taylor correction and its advantage over naive BCD (see Table 3 in the main paper).

Interpretation of Table 8.

• Column 2 (Active Set): We employ an active‑set strategy where updates are restricted to variables with nonzero coefficients, and use KKT conditions to govern active set expansion. When the block size $s$ exceeds the active-set size $|\mathcal{A}|$, we dynamically set $s \leftarrow \sqrt{|\mathcal{A}|}$ to improve efficiency.
• Column 3 (No Active Set): We turn the active-set feature off so that the block size, $s$, remains fixed at the initial setting. We report accuracy using relative objective difference.

The optimal setting for $s$ is typically a small constant determined by the cost of the link function and hardware features (e.g. DRAM memory bandwidth/latency, cache size, processor speed) and, as illustrated in Table 8, not due to numerical stability or accuracy. Our experiments in the main paper yielded $s \in [2, 32]$ as the most performant, so the additional memory footprint is bounded with performance (and not accuracy) being the upper bound on the best value of $s$. The large‑$s$ results will be included in the camera-ready version for completeness.

Table 8. ECCD Convergence Behavior under Fixed and Adaptive Block Sizes ($s$) on the Duke Dataset