GPU-friendly and Linearly Convergent First-order Methods for Certifying Optimal $k$-sparse GLMs

Jiachang Liu; Andrea Lodi; Soroosh Shafiee

math.OCMar 1, 2026

GPU-friendly and Linearly Convergent First-order Methods for Certifying Optimal $k$ -sparse GLMs

Authors: Jiachang Liu, Andrea Lodi, Soroosh Shafiee

Organizations: Center for Data Science for Enterprise and Society, Cornell University, Ithaca, USA · School of Operations Research and Information Engineering, Cornell University, Ithaca, USA · Jacobs Technion-Cornell Institute, Cornell Tech and Technion–IIT, New York, USA

Abstract

We investigate the problem of certifying optimality for sparse generalized linear models (GLMs), where sparsity is enforced through a cardinality constraint. While Branch-and-Bound (BnB) frameworks can certify optimality using perspective relaxations, existing methods for solving these relaxations are computationally intensive, limiting their scalability. To address this challenge, we reformulate the relaxations as composite optimization problems and develop a unified proximal framework that is both linearly convergent and computationally efficient. Under specific geometric regularity conditions, our analysis links primal quadratic growth to dual quadratic decay, yielding error bounds that make the Fenchel duality gap a sharp proxy for progress towards the solution set. This leads to a duality gap-based restart scheme that upgrades a broad class of sublinear proximal methods to provably linearly convergent methods, and applies beyond the sparse GLM setting. For the implicit perspective regularizer, we further derive specialized routines to evaluate the regularizer and its proximal operator exactly in log-linear time, avoiding costly generic conic solvers. The resulting iterations are dominated by matrix--vector multiplications, which enables GPU acceleration. Experiments on synthetic and real-world datasets show orders-of-magnitude faster dual-bound computations and substantially improved BnB scalability on large instances.

Figures & tables

Panel of Figure 1 — Figure 1: Illustration of quadratic growth for the primal objective $\Phi(\bm{\beta})$ and quadratic decay for the dual objective $\Psi(\bm{\zeta})$ .

Algorithm	Rate Function $r(t)$	Complexity with Restart
PGD	$\displaystyle\frac{L}{2t}$	$\displaystyle\mathcal{O}\left(\frac{\kappa L}{\alpha\sigma}\log\left(\frac{1}{\epsilon}\right)\right)$
FISTA	$\displaystyle\frac{2L}{\left(t+1\right)^{2}}$	$\displaystyle\mathcal{O}\left(\sqrt{\frac{\kappa L}{\alpha\sigma}}\log\left(\frac{1}{\epsilon}\right)\right)$
AC-FGM	$\displaystyle\frac{15cL}{2t\left(t+1\right)\rho\left(1-\rho\right)}$	$\displaystyle\mathcal{O}\left(\sqrt{\frac{\kappa L}{\alpha\sigma}}\log\left(\frac{1}{\epsilon}\right)+\log_{c}\left(\frac{L}{L_{0}}\right)\right)$

Table 1: Standard rate functions

r(t)

(when only the smoothness of

F

is assumed) and the corresponding complexity to reach

\epsilon

-accuracy after applying the duality-gap-based restart scheme.

Figure 2: 3D plots of $g_{\mathcal{N}}$ (left) and $g_{\mathcal{N}}^{*}$ (right) at the root node ( $p=2$ , $k=1$ , $M=1$ ). Left: red diamonds mark the boundary $\lvert\beta_{1}\rvert+\lvert\beta_{2}\rvert=M$ of $\dom(g_{\mathcal{N}})$ (top at $g_{\mathcal{N}}(\bm{\beta})=\frac{1}{2}M^{2}$ , bottom projection). Right: green squares mark the transition set $\max(\lvert\alpha_{1}\rvert,\lvert\alpha_{2}\rvert)=M$ (top at $g_{\mathcal{N}}^{*}(\bm{\alpha})=\frac{1}{2}M^{2}$ , bottom projection); $g_{\mathcal{N}}^{*}$ grows quadratically inside and linearly outside.

Panel of Figure 3 — Figure 3: A vector $\bm{\omega}$ majorizes $\bm{\beta}$ if and only if $\bm{\beta}\in\conv(\mathrm{perm}(\bm{\omega}))$ . Left ( $p=2$ , $k=1$ ): the hull is the segment between $(2,0)$ and $(0,2)$ , containing $\bm{\beta}=(1.4,0.6)$ . Right ( $p=3$ , $k=2$ ): the hull is the triangle with vertices $(1.5,1.5,0)$ , $(1.5,0,1.5)$ , and $(0,1.5,1.5)$ , containing $\bm{\beta}=(0.6,1.2,1.2)$ . Given $\bm{\beta}$ , Algorithm 1 constructs such a $k$ -sparse vector $\bm{\omega}$ .

Figure 4: Running time comparison of evaluating

g_{\mathcal{N}}

and the proximal operators of

g_{\mathcal{N}}

and

g_{\mathcal{N}}^{*}

, where

\mathcal{N}

is the root node. The baselines solve the corresponding SOCPs directly.

Figure 5: Running time comparison of solving Problem ( 3 ), the perspective relaxation at the root node. We set

M=2.0

\lambda_{2}=1.0

, and

n

-to-

p

ratio to be 1. Gurobi cannot solve the relaxation of the cardinality constrained logistic regression problem.

Figure 6: Convergence speed comparison between accelerated methods (with and without restart), on solving the perspective relaxation in Problem 3 with the linear regression loss,

n=16000

p=16000

k=10

\rho=0.5

\lambda_{2}=1.0

, and

M=2.0

Figure 7: Convergence speed comparison between accelerated methods (with and without restart), on solving the perspective relaxation in Problem 3 with the logistic regression loss,

n=16000

p=16000

k=10

\rho=0.5

\lambda_{2}=1.0

, and

M=2.0

$p$	1000	2000	4000	8000	16000
ours CPU	0.070	0.127	0.447	5.239	12.491
ours CPU	(0.001)	(0.002)	(0.075)	(0.632)	(1.914)
ours GPU	0.296	0.282	0.266	0.444	0.858
ours GPU	(0.003)	(0.002)	(0.001)	(0.002)	(0.003)

Table 2: GPU acceleration of our method on the linear regression task. Top and bottom rows correspond to the mean and standard deviation of running times (seconds).

$p$	1000	2000	4000	8000	16000
ours CPU	0.069	0.124	0.629	4.566	12.075
ours CPU	(0.006)	(0.002)	(0.028)	(0.129)	(0.824)
ours GPU	0.275	0.256	0.249	0.434	0.997
ours GPU	(0.013)	(0.015)	(0.013)	(0.017)	(0.029)

Table 3: GPU acceleration of our method on the logistic regression task. Top and bottom rows correspond to the mean and standard deviation of running times (seconds).

$p$	time(s)	gap(%)	# nodes	time(s)	gap(%)	# nodes	time(s)	gap(%)	# nodes
	ours			gurobi			mosek
1000	9 (1)	0.00	105	146	0.00	237	112	0.00	113
2000	6 (1)	0.00	59	562	0.00	100	470	0.00	57
4000	12 (6)	0.00	31	3146	0.00	34	2909	0.00	31
8000	25 (18)	0.00	21	TL	100.00	1	TL	21.00	4
16000	72 (53)	0.00	21	TL	100.00	1	OOM	OOM	OOM

Table 4: BnB results on synthetic (linear) with

k=10

\rho=0.5

\lambda_{2}=1.0

M=2.0

n/p=1.0

. Parentheses show lower-bound time. OOM = out of memory, TL = time

>7200

$p$	time(s)	gap(%)	# nodes	time(s)	gap(%)	# nodes	time(s)	gap(%)	# nodes
	ours			gurobiOA			mosek
1000	140 (7)	0.00	1121	TL	54.17	597933	386	0.00	1083
2000	110 (13)	0.00	595	TL	66.26	582581	1188	0.00	583
4000	126 (60)	0.00	219	TL	54.04	255320	2917	0.00	235
8000	180 (106)	0.00	131	TL	46.95	81405	OOM	OOM	OOM
16000	281 (187)	0.00	67	TL	35.04	21479	OOM	OOM	OOM

Table 5: BnB results on synthetic (logistic) with

k=10

\rho=0.5

\lambda_{2}=1.0

M=2.0

n/p=1.0

. Parentheses show lower-bound time. OOM = out of memory, TL = time

>7200

$k$	time(s)	gap(%)	# nodes	time(s)	gap(%)	# nodes	time(s)	gap(%)	# nodes
	ours			gurobi			mosek
6	345 (259)	0.00	1359	TL	100.00	1	TL	0.33	437
7	536 (395)	0.00	2125	TL	100.00	1	TL	0.43	497
8	1206 (908)	0.00	4393	TL	100.00	1	TL	0.36	393
9	3060 (2281)	0.00	10085	TL	100.00	1	TL	0.49	303
10	TL (5213)	0.01	18775	TL	100.00	1	TL	0.52	369

Table 6: BnB results on Santander (linear),

n=4459

p=4735

\lambda_{2}=1.0

M=10

. Parentheses show lower-bound time. OOM = out of memory, TL = time

>7200

$k$	time(s)	gap(%)	# nodes	time(s)	gap(%)	# nodes	time(s)	gap(%)	# nodes
	ours			gurobiOA			mosek
5	19 (15)	0.00	1	683	0.00	938	1075	0.00	0
10	66 (48)	0.00	23	2304	0.00	3118	OOM	OOM	OOM
15	116 (76)	0.00	33	3005	0.00	3209	OOM	OOM	OOM
20	223 (142)	0.00	61	TL	0.07	3680	OOM	OOM	OOM
25	629 (487)	0.00	167	TL	0.12	3635	OOM	OOM	OOM

Table 7: BnB results on Dorothea (logistic),

n=2300

p=89989

\lambda_{2}=1.0

M=10

. Parentheses show lower-bound time. OOM = out of memory, TL = time

>7200

$p$	time(s)	gap(%)	# nodes	time(s)	gap(%)	# nodes
	ours CPU			ours GPU
1000	8 (1)	0.00	105	10 (2)	0.00	105
2000	6 (1)	0.00	59	6 (1)	0.00	59
4000	4 (1)	0.00	31	4 (1)	0.00	31
8000	12 (7)	0.00	21	6 (1)	0.00	21
16000	31 (20)	0.00	21	15 (2)	0.00	21

Table 8: BnB CPU vs GPU results for linear regression with

k=10

\rho=0.5

\lambda_{2}=1.0

M=2.0

n/p=1.0

. Parentheses show lower-bound time.

$p$	time(s)	gap(%)	# nodes	time(s)	gap(%)	# nodes
	ours CPU			ours GPU
1000	132 (6)	0.00	1121	141 (20)	0.00	1123
2000	94 (7)	0.00	595	99 (11)	0.00	595
4000	64 (14)	0.00	219	55 (5)	0.00	219
8000	111 (45)	0.00	131	68 (5)	0.00	131
16000	160 (72)	0.00	67	88 (8)	0.00	67

Table 9: BnB CPU vs GPU results for logistic regression with

k=10

\rho=0.5

\lambda_{2}=1.0

M=2.0

n/p=1.0

. Parentheses show lower-bound time.

$k$	time(s)	gap(%)	# nodes	time(s)	gap(%)	# nodes
	ours CPU			ours GPU
6	169 (107)	0.00	1357	83 (23)	0.00	1363
7	267 (166)	0.00	2149	145 (46)	0.00	2143
8	600 (368)	0.00	4385	337 (111)	0.00	4387
9	1655 (1058)	0.00	10087	838 (267)	0.00	10121
10	5808 (3536)	0.00	34397	3135 (909)	0.00	34379

Table 10: BnB CPU vs GPU results on Santander (linear),

n=4459

p=4735

\lambda_{2}=1.0

M=10

. Parentheses show lower-bound time.

$k$	time(s)	gap(%)	# nodes	time(s)	gap(%)	# nodes
	ours CPU			ours GPU
5	12 (11)	0.00	1	3 (1)	0.00	1
10	52 (46)	0.00	23	9 (4)	0.00	23
15	82 (69)	0.00	33	21 (6)	0.00	33
20	167 (135)	0.00	61	37 (11)	0.00	63
25	466 (403)	0.00	167	84 (30)	0.00	165

Table 11: BnB CPU vs GPU results on Dorothea (logistic),

n=2300

p=89989

\lambda_{2}=1.0

M=10

. Parentheses show lower-bound time.

Appendix figures & tables20 assets

Supplementary material from the paper’s appendix.

Appendix

Figure D.1: Running time comparison of solving Problem 3 when the box constraint is tightened to

M=1.2

Figure D.2: Running time comparison of solving Problem 3 when the box constraint is relaxed to

M=5.0

Figure D.3: Running time comparison of solving Problem 3 with the

\ell_{2}

regularization reduced to

\lambda_{2}=0.1

Figure D.4: Running time comparison of solving Problem 3 with the

\ell_{2}

regularization increased to

\lambda_{2}=10

Figure D.5: Running time comparison of solving Problem 3 when the

n

-to-

p

ratio is increased to

10.0

Figure D.6: Running time comparison of solving Problem 3 when the

n

-to-

p

ratio is decreased to

0.1

Figure D.7: Convergence speed comparison of FISTA methods (with and without line search) and ACFGM (with duality-gap-based restarted ) for

\eta=e

\eta=e^{2}

, and

\eta=e^{3}

on solving the perspective relaxation in Problem 3 with the linear regression loss,

n=16000

p=16000

k=10

\rho=0.5

\lambda_{2}=1.0

, and

M=2.0

Figure D.8: Convergence speed comparison of FISTA methods (with and without line search) and ACFGM (with duality-gap-based restarted ) for

\eta=e

\eta=e^{2}

, and

\eta=e^{3}

on solving the perspective relaxation in Problem 3 with the logistic regression loss,

n=16000

p=16000

k=10

\rho=0.5

\lambda_{2}=1.0

, and

M=2.0

Figure D.9: Convergence speed comparison between accelerated methods (with and without restart), on solving the

\ell_{1}

-ball constrained problem with the least squares loss.

Figure D.10: Convergence speed comparison between accelerated methods (with and without restart), on solving the

\ell_{1}

-ball constrained problem with the logistic regression loss.

Figure D.11: Convergence speed comparison between accelerated methods (with and without restart), on solving the $\ell_{1}$ -regularized problem with the least squares loss.

Figure D.12: Convergence speed comparison between accelerated methods (with and without restart), on solving the $\ell_{1}$ -regularized problem with the logistic regression loss.

$p$	time(s)	gap(%)	# nodes	time(s)	gap(%)	# nodes	time(s)	gap(%)	# nodes
	ours			gurobi			mosek
1000	31 (3)	0.00	355	227	0.00	653	208	0.00	361
2000	16 (2)	0.00	161	2091	0.00	1544	758	0.00	171
4000	27 (15)	0.00	77	4574	0.00	97	3612	0.00	73
8000	63 (45)	0.00	49	TL	100.00	1	TL	33.57	4
16000	151 (115)	0.00	33	TL	100.00	1	OOM	OOM	OOM

Appendix

Table D.1: BnB results on synthetic (linear) with

k=10

\rho=0.5

\lambda_{2}=1.0

M=3.0

n/p=1.0

. Parentheses show lower-bound time. OOM = out of memory, TL = time

>7200

$p$	time(s)	gap(%)	# nodes	time(s)	gap(%)	# nodes	time(s)	gap(%)	# nodes
	ours			gurobiOA			mosek
1000	771 (37)	0.00	6049	TL	67.29	417432	1871	0.00	5847
2000	630 (81)	0.00	3451	TL	87.29	443958	7031	0.00	3427
4000	384 (168)	0.00	855	TL	82.12	269137	TL	13.97	569
8000	666 (404)	0.00	459	TL	69.86	93806	OOM	OOM	OOM
16000	1031 (677)	0.00	207	TL	56.55	35414	OOM	OOM	OOM

Appendix

Table D.2: BnB results on synthetic (logistic) with

k=10

\rho=0.5

\lambda_{2}=1.0

M=3.0

n/p=1.0

. Parentheses show lower-bound time. OOM = out of memory, TL = time

>7200

$p$	time(s)	gap(%)	# nodes	time(s)	gap(%)	# nodes
	ours CPU			ours GPU
1000	29 (2)	0.00	355	36 (9)	0.00	355
2000	15 (2)	0.00	161	18 (4)	0.00	161
4000	12 (4)	0.00	77	9 (2)	0.00	77
8000	44 (30)	0.00	49	15 (2)	0.00	49
16000	78 (53)	0.00	33	26 (4)	0.00	33

Appendix

Table D.3: BnB CPU vs GPU results for linear regression with

k=10

\rho=0.5

\lambda_{2}=1.0

M=3.0

n/p=1.0

. Parentheses show lower-bound time.

$p$	time(s)	gap(%)	# nodes	time(s)	gap(%)	# nodes
	ours CPU			ours GPU
1000	689 (38)	0.00	6049	758 (108)	0.00	6049
2000	536 (39)	0.00	3451	556 (65)	0.00	3451
4000	251 (62)	0.00	855	203 (21)	0.00	855
8000	412 (176)	0.00	459	246 (18)	0.00	459
16000	533 (276)	0.00	207	282 (23)	0.00	207

Appendix

Table D.4: BnB CPU vs GPU results for logistic regression with

k=10

\rho=0.5

\lambda_{2}=1.0

M=3.0

n/p=1.0

. Parentheses show lower-bound time.

Figure D.13: Beamsearch cross-validation results on the SantanderBank dataset with

\lambda_{2}\in\{0.01,0.1,1.0,2.0,3.0,10.0\}

Figure D.14: Beamsearch cross-validation results on the DOROTHEA dataset with

\lambda_{2}\in\{0.01,0.1\}

Figure D.15: Beamsearch cross-validation results on the DOROTHEA dataset with

\lambda_{2}\in\{1.0,2.0\}

Figure D.16: Beamsearch cross-validation results on the DOROTHEA dataset with

\lambda_{2}\in\{3.0,10.0\}

Constrained Optimization

Constrained optimization finds the best decisions while satisfying equality, inequality, resource, safety, or fairness requirements.

See all 335 papers

May 21, 2026

From Sequential Nodes to GPU Batches: Parallel Branch and Bound for Optimal $k$ -Sparse GLMs

GPUs have significantly accelerated first-order methods for large-scale optimization, especially in continuous optimization. However, this success has not transferred cleanly to problems with discrete variables, combinatorial structure, and nonlinear objectives, such as certifying optimal solutions for cardinality-constrained generalized linear models. Major challenges include the sequential processing of heterogeneous nodes in branch and bound (BnB) and frequent data movement between the CPU and GPU. We propose a simple, generic, and modular CPU--GPU framework that processes multiple BnB nodes in batches on GPUs. The framework is built around a small set of GPU-efficient routines and uses padding together with lightweight custom kernels to handle irregular node data structures. Experiments show one to two orders of magnitude speedups and zero optimality gap on challenging instances. The framework can also be extended to collect the entire Rashomon set, enabling downstream statistical analysis such as variable-importance analysis and model selection under secondary user-specific measures (e.g., AUC in classification).

Jiachang Liu, Andrea Lodi

∗Center for Data Science for Enterprise and Society, Cornell University, Ithaca, USA · †School of Operations Research and Information Engineering, Cornell University, Ithaca, USA · ‡Jacobs Technion-Cornell Institute, Cornell Tech and Technion–IIT, New York, USA

Branch and Bound GPU Kernel Optimization GPU Acceleration High-Performance Computing Combinatorial Optimization Constrained Optimization+1

Ask PDF

Oct 1, 2026

Reinforcement Learning to Accelerate Primal-Dual Hybrid Gradient for Linear Programming

Primal-dual hybrid gradient (PDHG) methods solve large-scale linear programs (LPs) using GPU-friendly matrix-vector products and projections, but their practical performance depends on coordinating algorithm parameters, acceleration, and restarts. We introduce GALLOP, which uses reinforcement learning to jointly learn continuous algorithm parameters and discrete restart decisions without differentiating through the solver. Its generalized accelerated PDHG update combines separate primal and dual extrapolation, history corrections, and restart anchoring with independently adjustable coefficients. We train a dimension-agnostic feedback policy using a groupwise proximal policy optimization objective that clips likelihood ratios separately for different control groups and excludes inactive acceleration controls on restart transitions. We evaluate GALLOP on six LP families and a public item-placement benchmark. On the main evaluation settings across the six families, GALLOP reduces iteration counts by factors of

1.9

5.6

and achieves up to a

16.0\times

speedup in algorithm wall-clock time over MPAX. With one policy trained per family, the learned policies generalize without retraining to within-family LPs

3\times

400\times

larger than the largest training instances, including Transport LPs with

10.24

million variables.

Jinhwan Sul, Alex Oshin, Evangelos A. Theodorou

Georgia Institute of Technology, USA

Primal-Dual Optimization Reinforcement Learning Linear Programming Accelerated Gradient Methods

Ask PDF

Sep 30, 2026

Generalized Geometry Block Proximal Linearized Method for Multiblock Nonconvex and Nonsmooth Optimization

This paper considers a class of multiblock nonconvex and nonsmooth optimization problems arising in many applications. Existing methods construct proximal linearized operators or their variants within standard Euclidean geometry to solve this class of problems, forcing their block variable updates to rely on the standard inner product and its induced norm. Nevertheless, this construction fails to capture the geometric structure of the target problem, leading to low numerical efficiency. To overcome these drawbacks, we propose a generalized geometry proximal linearized operator for updating block variables, and develop the Generalized Geometry Block Proximal Linearized (GGBPL) method based on this operator. Compared with existing proximal linearized operators, the proposed operator allows the block surrogate functions to be constructed using arbitrary inner products and general admissible metrics, thereby enabling the GGBPL method to adapt its updates to the geometric structure of various problems. We also introduce the inertial version of GGBPL, named the inertial GGBPL (iGGBPL) method. We further establish a new unified convergence framework under this generalized geometry, within which we prove that our methods guarantee convergence of the objective function values, establish global convergence of the generated sequence to a critical point, and derive the convergence rate of our methods. We also establish an

\mathcal{O}(\varepsilon^{-2})

iteration complexity bound for obtaining an

\varepsilon

-stationary point. We apply our methods to two nonconvex and nonsmooth problems: sparse nonnegative matrix factorization with

\ell_0

-constraints and sparse nonnegative CP decomposition with

\ell_0

-constraints. Numerical results demonstrate the superior numerical performance of our proposed methods over several state-of-the-art methods.

Weifeng Yang

Nonconvex Optimization Constrained Optimization

Ask PDF

Abstract

Constrained Optimization

Explore similar work

From Sequential Nodes to GPU Batches: Parallel Branch and Bound for Optimal kkk-Sparse GLMs

Reinforcement Learning to Accelerate Primal-Dual Hybrid Gradient for Linear Programming

Generalized Geometry Block Proximal Linearized Method for Multiblock Nonconvex and Nonsmooth Optimization

From Sequential Nodes to GPU Batches: Parallel Branch and Bound for Optimal $k$ -Sparse GLMs