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

From Sequential Nodes to GPU Batches: Parallel Branch and Bound for Optimal $k$ -Sparse GLMs
May 21, 2026Jiachang Liu, Andrea LodiBranch and Bound GPU Kernel Optimization
∗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
Reinforcement Learning to Accelerate Primal-Dual Hybrid Gradient for Linear Programming
Oct 1, 2026Jinhwan Sul, Alex Oshin, Evangelos A. TheodorouPrimal-Dual Optimization Reinforcement Learning
Georgia Institute of Technology, USA
Generalized Geometry Block Proximal Linearized Method for Multiblock Nonconvex and Nonsmooth Optimization
Sep 30, 2026Weifeng YangNonconvex Optimization Constrained Optimization

May 21, 2026J/K move · Enter open · S save

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

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

Open paper AskPDF arXiv

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

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