Pith. sign in

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 →

arxiv 2502.09502 v3 pith:G6GZ45XH submitted 2025-02-13 cs.LG math.OC

classification cs.LGmath.OC MSC 90C2590C2690C27
keywords k-sparseGLMoptimalitycertificationperspectiverelaxationbranchandboundFISTAproximaloperatorpooledadjacentviolatorsalgorithmisotonicregression
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

This paper targets the bottleneck in certifying optimality for $k$-sparse generalized linear models: at every branch-and-bound node one must compute a tight lower bound, and existing conic solvers are either expensive or slow to converge. The authors reformulate the perspective relaxation of the cardinality-constrained problem as an unconstrained convex composite program and prove (Theorem 3.2) that the proximal operator of the nonsmooth component's conjugate can be evaluated exactly by a customized pooled-adjacent-violators algorithm in $O(p\log p)$ time, replacing a second-order cone program. A value-based restart strategy makes the accelerated proximal-gradient method converge linearly in practice, and the per-iteration work reduces to matrix-vector products, so the computation maps onto GPUs. On synthetic and real data the method reports dual-bound speedups of one to two orders of magnitude and certifies optimality for large instances, including a logistic regression problem with about 90,000 features. The certification arguments depend on the safe lower bound of Section 3.4; the weakest-assumption note records a sign issue that appears to affect the logistic-loss formula as printed.

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.

Watch

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

Editorial extensions of the paper, not claims the author makes directly.

  • 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.
Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

1 major / 6 minor

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)
  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)
  1. [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.
  2. [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.
  3. [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.
  4. [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.
  5. [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.
  6. [Section 3] There is a minor typo in the reference to Rockafellar (1970): "Theorem 5,3" should read "Theorem 5.3".

Circularity Check

0 steps flagged · score 0.0 of 10

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 0 free parameters · 4 assumptions · 0 invented entities

The derivation relies on standard convex analysis and isotonic regression theory, plus the stated smoothness assumption on the loss. The convex hull representation (19) is taken from prior work. No new entities or ad hoc assumptions are introduced.

assumptions (4)
  • domain assumption The GLM loss f is Lipschitz smooth.
    Stated in Section 1; required for FISTA's convergence and step-size selection in Algorithm 3.
  • standard math Standard convex analysis: conjugate calculus, Moreau decomposition, weak duality.
    Used in Section 3 and Appendix D; from Rockafellar (1970) and Beck (2017).
  • standard math Optimality and single-valuedness properties of PAVA for separable convex isotonic regression.
    Underpins Lemmas 3.4 and 3.5; from Best et al. (2000) and Ahuja & Orlin (2001).
  • domain assumption The majorization-based convex hull representation (19) of the sign- and permutation-invariant set S0.
    Invoked in the proof of Theorem 3.6; cited to Kim et al. (2022, Theorem 4).

how reviews work

0 comments
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 reproduced from arXiv: 2502.09502 by the authors.

Figure 1
Figure 1. Running time comparison of evaluating the proximal operators, for both g (left) and g ∗ (right). The baselines evaluate the proximal operators by directly solving the corresponding second-order conic problems (SOCP), respectively. 2500 5000 7500 10000 12500 15000 P 10 1 10 0 10 1 10 2 10 3 Solver Time Elapsed Linear Loss 2500 5000 7500 10000 12500 15000 P 10 1 10 0 10 1 10 2 10 3 Solver Time Elapsed Logistic Loss gu… view at source ↗
Figure 2
Figure 2. Running time comparison of solving Problem (5), the perspective relaxation of the original MIP problem. We set M = 2.0, λ2 = 1.0, and n-to-p ratio to be 1. Gurobi cannot solve the cardinality constrained logistic regression problem [PITH_FULL_IMAGE:figures/full_fig_p007_2.png] view at source ↗
Figure 3
Figure 3. Empirical convergence rate of our restarted FISTA (compared with PGD, the proximal gradient method, and FISTA) on solving the perspective relaxation in Problem (5) with the logistic loss, n = 16000, p = 16000, k = 10, ρ = 0.5, λ2 = 1.0, and M = 2.0 [PITH_FULL_IMAGE:figures/full_fig_p008_3.png] view at source ↗
Figures from the paper (10 more)
Figure 4
Figure 4. Figure 4: Solve the perspective relaxation in Problem (5). We set M = 1.2, λ2 = 1.0, n/p = 1, and k = 10. 2500 5000 7500 10000 12500 15000 P 10 1 10 0 10 1 10 2 10 3 Solver Time Elapsed Linear Loss 2500 5000 7500 10000 12500 15000 P 10 1 10 0 10 1 10 2 10 3 Solver Time Elapsed L…
Figure 5
Figure 5. Figure 5: Solve the perspective relaxation in Problem (5). We set M = 1.5, λ2 = 1.0, n/p = 1, and k = 10. 20 [PITH_FULL_IMAGE:figures/full_fig_p020_5.png]
Figure 6
Figure 6. Figure 6: Solve the perspective relaxation in Problem (5). We set M = 3.0, λ2 = 1.0, n/p = 1, and k = 10. 2500 5000 7500 10000 12500 15000 P 10 1 10 0 10 1 10 2 10 3 Solver Time Elapsed Linear Loss 2500 5000 7500 10000 12500 15000 P 10 0 10 1 10 2 10 3 Solver Time Elapsed Logist…
Figure 7
Figure 7. Figure 7: Solve the perspective relaxation in Problem (5). We set M = 5.0, λ2 = 1.0, n/p = 1, and k = 10. 2500 5000 7500 10000 12500 15000 P 10 0 10 1 10 2 10 3 Solver Time Elapsed Linear Loss 2500 5000 7500 10000 12500 15000 P 10 1 10 0 10 1 10 2 10 3 Solver Time Elapsed Logist…
Figure 8
Figure 8. Figure 8: Solve the perspective relaxation in Problem (5). We set M = 10.0, λ2 = 1.0, n/p = 1, and k = 10. 21 [PITH_FULL_IMAGE:figures/full_fig_p021_8.png]
Figure 9
Figure 9. Figure 9: Solve the perspective relaxation in Problem (5). We set M = 2.0, λ2 = 0.1, n/p = 1, and k = 10. 2500 5000 7500 10000 12500 15000 P 10 1 10 0 10 1 10 2 10 3 Solver Time Elapsed Linear Loss 2500 5000 7500 10000 12500 15000 P 10 1 10 0 10 1 10 2 10 3 Solver Time Elapsed L…
Figure 10
Figure 10. Figure 10: Solve the perspective relaxation in Problem (5). We set M = 2.0, λ2 = 10.0, n/p = 1, and k = 10. 22 [PITH_FULL_IMAGE:figures/full_fig_p022_10.png]
Figure 11
Figure 11. Figure 11: Solve the perspective relaxation in Problem (5). We set M = 2.0, λ2 = 1.0, n/p = 10.0, and k = 10. 2500 5000 7500 10000 12500 15000 P 10 1 10 0 10 1 10 2 10 3 Solver Time Elapsed Linear Loss 2500 5000 7500 10000 12500 15000 P 10 1 10 0 10 1 10 2 10 3 Solver Time Elaps…
Figure 12
Figure 12. Figure 12: Solve the perspective relaxation in Problem (5). We set M = 2.0, λ2 = 1.0, n/p = 0.1, and k = 10. 23 [PITH_FULL_IMAGE:figures/full_fig_p023_12.png]
Figure 13
Figure 13. Figure 13: Solve the perspective relaxation in Problem (5). We set M = 2.0, λ2 = 1.0, n/p = 0.1, and k = 500. 24 [PITH_FULL_IMAGE:figures/full_fig_p024_13.png]

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

13 extracted references · 9 canonical work pages

  1. [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 ...

  2. [5]

    and Yang, J

    Lu, H. and Yang, J. A practical and optimal first-order method for large-scale convex quadratic programming. arXiv:2311.07710,

  3. [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]...

  4. [9]

    Proof of Theorem 3.6 Theorem 3.6.For anyβ∈R p, Algorithm 2 computes the exact value ofg(β), defined in(7), inO(p+plogk)

    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...

  5. [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...

  6. [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...

  7. [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 ...

  8. [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 ⋆ ...

Show all 13 references
  1. [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,

  2. [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,

  3. [2021]

    and Neto, J

    G´omez, A. and Neto, J. Outlier detection in regression: Conic quadratic formulations.arXiv:2307.05975,

  4. [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,

  5. [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,

Pith tools

Reviewed August 7, 2026 · model on record in the stance chip above.