REVIEW 1 major objections 6 minor 13 references
Scalable First-order Method for Certifying Optimal k-Sparse GLMs
T0 review · 1 major / 6 minor · reviewed 2026-08-07 · deepseek-v4-flash
Pith's one-line read This paper claims an exact O(p log p) proximal evaluation that makes branch-and-bound optimality certification of k-sparse GLMs orders of magnitude faster, including logistic problems with about 90,000 features.
desk verdict A clean exact PAVA for the box-constrained perspective conjugate plus big speedups, but the logistic safe-lower-bound derivation as printed is vacuous because of a sign error that must be fixed before the certification claims hold. read the letter →
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
What carries the argument
The perspective relaxation is obtained by writing each squared coefficient as a perspective term $\beta_j^2/z_j$ and relaxing the binary indicators to $[0,1]$; its core machinery is the proximal operator of $g^*$, the conjugate of the implicit perspective-sum function $g(\beta)=\min_z \frac12\sum_j \beta_j^2/z_j$ subject to $\sum_j z_j\le k$, $z\in[0,1]^p$, and $-M z_j\le \beta_j\le M z_j$. Lemma 3.1 identifies $g^*(\alpha)=\mathrm{TopSum}_k(H_M(\alpha))$, where $H_M$ is the Huber loss clipped at $M$; Algorithm 1, a customized pooled adjacent violators algorithm, solves $\mathrm{prox}_{\rho g^*}(\mu)$ exactly by sorting entries of $|\mu|$, grouping them into isotonic blocks, and merging blocks when order constraints are violated. This replaces the second-order cone program that a generic proximal evaluation would require, and the same block machinery plus a majorization characterization of the permutation-invariant convex hull computes $g(\beta)$ exactly in $O(p+p\log k)$ for the restart test. The FISTA iteration then needs only matrix-vector products, and Moreau decomposition turns the exact conjugate prox into the primal prox needed by the gradient step.
What would settle it
Evaluate bound (13) on one reported logistic dataset: take $\hat\beta$ from Algorithm 3, set $\zeta=\nabla F(X\hat\beta)$, and compute $F^*(-\zeta)$ with the appendix's logistic formula. If any $\zeta_i/y_i$ lies outside $[0,1]$, the bound is $-\infty$ and the certificate is vacuous; recomputing with $\zeta=-\nabla F(X\hat\beta)$ gives a finite value. The comparison on a single instance distinguishes a typo from a substantive gap.
Extended reading notes
Core claim
The paper's central claim is that the perspective relaxation of the $k$-sparse GLM problem can be solved fast enough to make optimality certification practical. Its key theorem states that for any $\mu\in\mathbb{R}^p$, Algorithm 1 returns the exact value of $\mathrm{prox}_{\rho g^*}(\mu)$ in $O(p\log p)$, where $g$ is the implicit function $g(\beta)=\min_z \frac12\sum_{j=1}^p \beta_j^2/z_j$ over $z\in[0,1]^p$, $\sum_j z_j\le k$, and $-M z_j\le \beta_j\le M z_j$. The proof identifies the conjugate $g^*(\alpha)=\mathrm{TopSum}_k(H_M(\alpha))$ with the $k$ largest entries of a Huber loss vector, converts the proximal problem into an isotonic regression, and solves it exactly with a customized PAVA. Using Moreau decomposition, this yields an exact proximal step for FISTA; an efficient exact routine for $g(\beta)$ enables value-based restarts; and weak duality converts the approximate FISTA output into a safe lower bound for the original mixed-integer program. The paper reports that the resulting method certifies zero-gap optimality on datasets where commercial solvers time out, with empirically linear convergence.
Load-bearing premise
Everything rests on the safe lower bound (13) being finite and valid; as printed, the logistic-loss version of that bound can evaluate to $-\infty$ because the stated choice $\zeta=\nabla F(X\hat\beta)$ falls outside the domain of the appendix's conjugate formula, so a sign correction such as $\zeta=-\nabla F(X\hat\beta)$ is needed for the certification claims to hold.
Editorial extensions
If this is right
- Exact prox evaluation in $O(p\log p)$ removes the repeated SOCP solves that made proximal methods for this relaxation impractical.
- The method produces valid lower bounds for pruning, so a minimal branch-and-bound can certify optimality on instances where commercial solvers exceed time limits.
- With per-iteration cost limited to matrix-vector products, GPU implementations yield additional order-of-magnitude speedups in high dimensions.
- The perspective relaxation is tighter than the $\ell_1$ relaxation while remaining first-order tractable, so certificates do not require expensive conic relaxations.
- Restarted FISTA exhibits empirically linear convergence on the perspective relaxation, a first for first-order lower-bound computation in this mixed-integer class.
Reading between the lines
- Beyond the paper, the exact PAVA solver for $\mathrm{prox}_{\rho g^*}$ could be reused as a module for other composite penalties that mix top-$k$ selection with Huber-type regularization, such as sparse grouped or monotone-structured models.
- Beyond the paper, the value-based restart's empirical linear convergence is likely provable under a quadratic-growth condition on the perspective objective; a convergence-rate theorem would strengthen the practical claim.
- Beyond the paper, the printed sign of $\zeta$ in Section 3.4 appears to make the logistic conjugate $F^*(-\zeta)$ equal $-\infty$; with $\zeta=-\nabla F(X\hat\beta)$, or an equivalent correction, the reported logistic certificates would be validated.
- Beyond the paper, the GPU-friendly formulation suggests that distributed matrix-vector products could push certification to $p\gg n$ regimes, where coordinate-descent and ADMM baselines become sequential bottlenecks.
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. The paper proposes a first-order proximal-gradient method for computing dual lower bounds in branch-and-bound certificates for k-sparse GLMs. The method works with the perspective relaxation, rewrites it as an unconstrained composite minimization problem, and derives a closed-form conjugate of the implicit regularizer g as well as a customized PAVA procedure that evaluates the relevant proximal operator exactly in O(p log p). A restart strategy and an exact O(p + p log k) evaluation of g complete the algorithm, and a safe lower bound for pruning is derived from weak duality. Experiments report large speedups over SOCP solvers and successful optimality certification on synthetic and real datasets, including DOROTHEA with about 90k features.
Significance. The exact proximal-operator result (Theorem 3.2) and the conjugate characterization (Lemma 3.1) are clean, self-contained, and likely to be useful beyond this specific application; the proofs are detailed and the code is publicly available. If the logistic certificate issue is corrected, the paper constitutes a solid algorithmic contribution to practical MIP certification for sparse GLMs. The empirical linear-convergence claim is interesting but is not supported by a theorem, so it should be presented as an observed property rather than a theoretical novelty.
major comments (1)
- [Section 3.4, Eq. (13), Appendix D.1] The safe lower-bound computation is vacuous for the logistic loss exactly as stated. With ζ=∇F(Xβ̂) and F(Xβ)=Σ_i log(1+exp(-y_i (Xβ)_i)), one gets ∇F_i = -y_i / (1+exp(y_i (Xβ̂)_i)), so ζ_i/y_i ∈ (-1,0). The conjugate formula given in D.1 for F^*(-ζ) is Σ_i [(1 - ζ_i/y_i) log(1 - ζ_i/y_i) + (ζ_i/y_i) log(ζ_i/y_i)], which is finite only when ζ_i/y_i ∈ [0,1]; for ζ_i/y_i < 0, the extended-real value is +∞. Consequently, the right-hand side of Eq. (13) is -∞, giving no certificate. The natural fix is to choose ζ = -∇F(Xβ̂), which places ζ_i/y_i ∈ (0,1] and makes the bound finite; this also matches the linear-regression example in D.2, where ζ̂ = -∇F(β̂). Please correct the sign in Section 3.4, state the logistic conjugate coherently, and verify that the BnB implementation in Appendix B.3 used the corrected choice, since the logistic certification results in Table 2 rely on this derivation.
minor comments (6)
- [Appendix A.2] In the proof of Lemma 3.1, the sentence "the second equality holds as z is a binary vector" is incorrect: z is continuous in [0,1] in the definition of g. The algebra actually holds for every z_j ∈ [0,1], so the formula remains valid, but the stated reason should be fixed.
- [Section 3.4] The notation ζ̂ appears in Eq. (13) without being explicitly defined; it presumably denotes ζ=∇F(Xβ̂), but the manuscript should use a single symbol consistently to avoid confusion about which vector is substituted into the weak-duality bound.
- [Section 3.3] The statement that restarted FISTA achieves "the first linear convergence result of using a first-order method in the MIP context" is an empirical observation, not a theorem. Please temper this claim or state the precise conditions under which linear convergence is expected, and cite relevant prior work on restarted FISTA.
- [Table 2, Section 4.2] Table 2 reports MOSEK runtimes of 1915 and 2446 seconds for datasets that Section 4.2 states are run with a 1800-second time limit. Clarify the actual time limit or explain why these entries exceed it.
- [Algorithm 3] The input line for Algorithm 3 refers to the Lipschitz-continuity parameter of ∇F(β) without defining F. For clarity, state that F(β) = f(Xβ,y) and specify how the step size L is estimated or set in the experiments.
- [Section 3] There is a minor typo in the reference to Rockafellar (1970): "Theorem 5,3" should read "Theorem 5.3".
Circularity Check
No circularity: the proximal operator, conjugate, and safe lower bound are derived from definitions and external isotonic-regression results; the self-citations are comparative and not load-bearing.
full rationale
The derivation chain is self-contained rather than circular. Lemma 3.1 computes g* directly from the definition of g in (7) by solving the inner maximization over β for fixed z and then optimizing over z; the result TopSum_k(H_M(α)) is derived, not assumed. Theorem 3.2 rests on Lemmas 3.3–3.5, which reduce prox_{ρg*} to an isotonic regression problem; the PAVA correctness is argued from standard isotonic-regression optimality properties (Best et al., 2000; Ahuja & Orlin, 2001; Busing, 2022), and the O(p log p) bound follows from sorting plus a linear merging argument. The safe lower bound (13) is obtained from weak duality with ζ = ∇F(Xβ̂); this is a textbook dual bound, and no fitted parameter is relabeled as a prediction. Lemma 2.1 cites Shafiee & Kilinc-Karzan (2024), a prior paper sharing an author, but that citation is to a published, parameter-free convex-hull characterization and is used as a standard building block; it does not smuggle in the paper's own claims, and the central contributions (proximal operator, restart, BnB bounds) do not depend on a self-citation chain. The logistic conjugate domain issue in §3.4/D.1 (ζ_i/y_i < 0 making F*(−ζ) infinite) is a correctness or sign concern, not a circularity: it does not make any output equal to its input. No step was found where a definition, fitted value, or self-citation forces the claimed result.
Assumptions & free parameters
assumptions (4)
- domain assumption The GLM loss f is Lipschitz smooth.
- standard math Standard convex analysis: conjugate calculus, Moreau decomposition, weak duality.
- standard math Optimality and single-valuedness properties of PAVA for separable convex isotonic regression.
- domain assumption The majorization-based convex hull representation (19) of the sign- and permutation-invariant set S0.
Cite this review
Pith. "Pith review of Scalable First-order Method for Certifying Optimal k-Sparse GLMs." pith.science (2026). https://pith.science/paper/G6GZ45XH
@misc{pith2026250209502,
author = {Pith},
title = {Pith review of: Scalable First-order Method for Certifying Optimal k-Sparse GLMs},
year = {2026},
howpublished = {\url{https://pith.science/paper/G6GZ45XH}},
note = {Machine review of arXiv:2502.09502}
}
abstract
This paper investigates the problem of certifying optimality for sparse generalized linear models (GLMs), where sparsity is enforced through an $\ell_0$ cardinality constraint. While branch-and-bound (BnB) frameworks can certify optimality by pruning nodes using dual bounds, existing methods for computing these bounds are either computationally intensive or exhibit slow convergence, limiting their scalability to large-scale problems. To address this challenge, we propose a first-order proximal gradient algorithm designed to solve the perspective relaxation of the problem within a BnB framework. Specifically, we formulate the relaxed problem as a composite optimization problem and demonstrate that the proximal operator of the non-smooth component can be computed exactly in log-linear time complexity, eliminating the need to solve a computationally expensive second-order cone program. Furthermore, we introduce a simple restart strategy that enhances convergence speed while maintaining low per-iteration complexity. Extensive experiments on synthetic and real-world datasets show that our approach significantly accelerates dual bound computations and is highly effective in providing optimality certificates for large-scale problems.
Figures
Figures from the paper (10 more)
Reference graph
Works this paper leans on
-
[1]
We then conclude the proof using these observations
and then establish that for every j, l∈[p] with |µj| ≥ |µl|, we have |α⋆ j | ≥ |α⋆ l |(step 2). We then conclude the proof using these observations. ⋄ Step 1.We prove the sign-preserving property through a proof by contradiction. For the sake of contradiction, suppose that there exists some j∈[p] such that sgn(α⋆ j )̸= sgn(µ j). Hence, we can construct a ...
work page 2000
-
[5]
Lu, H. and Yang, J. A practical and optimal first-order method for large-scale convex quadratic programming. arXiv:2311.07710,
-
[7]
Proofs This section contains all omitted proofs in the paper
11 Scalable First-order Method for Certifying Optimal k-Sparse GLMs Appendix Appendix A. Proofs This section contains all omitted proofs in the paper. A.1. Proof of Lemma 2.1 Lemma 2.1.The closed convex hull of the set n (τ,β,z) ∥β∥∞ ≤M,z∈ {0,1} p,1 ⊤z≤k, β j(1−z j) = 0∀j∈[p], P j∈[p] β2 j ≤τ o is given by the set n (τ,β,z) −M zj ≤β j ≤M zj ∀j∈[p],z∈[0,1]...
work page 2024
-
[9]
3:ν prev = ˆν1,j= 2 4:whilej≤ndo 5:b=b+ 1 6:P b =ρ j,S b =|µ j|,N b = 1,ν= ˆν j 7:▷ If the value for the current singleton block is greater that of the previous block (isotonic violation), merge the current block with the previous block 8:ifν > v prev then 9:b=b−1 10:P b =P b +ρ j,S b =S b +|µ j|,N b =N b + 1,ν=prox Pb Nb HM ( Sb Nb ) 11:▷Look forward: ke...
work page 2022
-
[10]
For this experiment, we vary the feature dimension p∈ {1000,2000,4000,8000,16000}
=P and P(yi =−1) = 1−P . For this experiment, we vary the feature dimension p∈ {1000,2000,4000,8000,16000} . We control the sample size by using a parameter called n-to-p ratio, or sample to feature ratio. For the results in the main paper, we set n-to-p ratio to be 1.0, the box constraint M to be 2, the number of nonzero coefficients k (also the cardinal...
work page 2000
-
[11]
The intuition is that such a variable is important and should be branched early in the BnB framework. B.4. Computing Platforms When investigating how much GPU can accelerate our computation, we ran the experiments with both CPU and GPU implementations on the Nvidia A100s. For everything else, we ran the experiments with the CPU implementation on AMD Milan...
work page 2023
-
[13]
⋄Squared Hinge Loss:For binary classification with labelsy i ∈ {−1,+1}, F(Xβ) = X i∈[n] max(0,1−y i(Xβ) i)2 &F ∗(−ζ) = X i∈[n] h(−yiζi), whereh(z) =z+ z2 4 ifz≤0andh(z) =∞ifz >0. 26 Scalable First-order Method for Certifying Optimal k-Sparse GLMs ⋄ Multinomial Logistic Regression:For multiclass classification with K classes with coefficients β∈R p×K, let ...
work page 2024
-
[500]
24 Scalable First-order Method for Certifying Optimal k-Sparse GLMs C.2. Comparison between Perspective Relaxation andℓ 1-Relaxation Besides the perspective relaxation discussed in Section 2, there is another common relaxation for the sparse learning problem, which is theℓ 1-relaxation. Recall that the perspective relaxation problem is formulated as: P ⋆ ...
work page 2020
Show all 13 references
-
[1999]
On the power of linear programming for K-means clustering
De Rosa, A., Khajavirad, A., and Wang, Y . On the power of linear programming for K-means clustering. arXiv:2402.01061,
-
[2020]
Outlier detection in time series via mixed-integer conic quadratic optimization.SIAM Journal on Optimiza- tion, 31(3):1897–1925,
G´omez, A. Outlier detection in time series via mixed-integer conic quadratic optimization.SIAM Journal on Optimiza- tion, 31(3):1897–1925,
1925
-
[2021]
and Neto, J
G´omez, A. and Neto, J. Outlier detection in regression: Conic quadratic formulations.arXiv:2307.05975,
-
[2023]
cuPDLP-C: A strengthened implementation of cuPDLP for linear programming by C language.arXiv:2312.14832,
Lu, H., Yang, J., Hu, H., Huangfu, Q., Liu, J., Liu, T., Ye, Y ., Zhang, C., and Ge, D. cuPDLP-C: A strengthened implementation of cuPDLP for linear programming by C language.arXiv:2312.14832,
-
[2024]
Accelerating low-rank factorization-based semidefinite programming algorithms on GPU.arXiv:2407.15049,
Han, Q., Lin, Z., Liu, H., Chen, C., Deng, Q., Ge, D., and Ye, Y . Accelerating low-rank factorization-based semidefinite programming algorithms on GPU.arXiv:2407.15049,
Reviewed August 7, 2026 · model on record in the stance chip above.
Discussion (0). Continue with ORCID to comment.