Scalable First-order Method for Certifying Optimal k-Sparse GLMs
Authors: Jiachang Liu
Soroosh Shafiee
Andrea Lodi
Venue: ICML 2025 poster (2025)
Keywords: generalized linear models
sparse learning
proximal method
mixed-integer programming
Primary area: optimization->large_scale_parallel_and_distributed

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.

Scalable First-order Method for Certifying Optimal k-Sparse GLMs
Jiachang Liu' Soroosh Shafiee
Andrea Lodi 2
Abstract
is critical in high-stakes applications. Recently, using the to norm to induce sparsity has gained significant attention. This approach provides distinct advantages over traditional convex relaxation methods, such as replacing &u with 4, particularly in cases involving highly correlated features.
This paper investigates the problem of certifying optimality for sparse generalized linear models (GLMs), where sparsity is enforced through an to cardinality constraint. While branch-and- aanuchhand-buun and-bound (BnB) frameworks can certify optimality by prun- ing nodes using dual bounds, existing methods for computing these bounds either computa- tionally intensive exhibit slow convergence, limiting their scalability to large-scale problems. To address this challenge, we propose first-order proximal gradient algorithm designed to solve the perspective relaxation of the problem within a BnB framework. Specifically, we formulate the relaxed problem composite optimization prob- lem and demonstrate that the proximal operator of the non-smooth component can be computed exactly in log-linear time complexity, climinat- ing the need to solve a computationally expen- sive second-order cone program. Furthermore. we introduce simple restart strategy that en- hances convergence speed while maintaining low per-iteration complexity. Extensive experiments 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.
In this paper, we aim to solve
min
f(XB, + X2 1312
BER
s.t.
|| S M,, 11 S k,
where X R and y E R" denote the matrix of features and the vector of labels, respectively, while the parameter M >> O can be either user-defined based on prior knowledge or estimated from the data (Park & Klabjan, 2020). The GLM loss function, denoted by f : R" x R + R. is assumed to be Lipschitz smooth, the parameter controls the number of nonzero coefficients, and A2 0 is small Tikhonov regularization coefficient to address collinearity. However, problem (1) is NP-hard (Natarajan, 1995). As a result, most existing methods rely on heuristics that deliver high-quality approximations but lack guarantees of optimality. This limitation is particularly problematic in high-stakes applications like healthcare, where ensuring accuracy, reliability, and safety is essential. Therefore, we emphasize the pursuit of certifiably optimal solutions.
N
A naive approach to solve (1) to optimality is to reformulate it as mixed-integer programming (MIP) problem and use commercial MIP solvers. However, these solvers struggle with scalability, especially for large datasets or nonlinear objectives. A key bottleneck is the computation of tight lower bounds at each branch-and-bound (BnB) node, which is essential for effective pruning and solver performance. Existing methods for computing lower bounds typically rely on linear or conic optimization. However, these approaches either generate loose bounds that reduce pruning efficiency or result in high computational costs per iteration. Addition- ally, they are difficult to parallelize, limiting their ability to leverage modern hardware accelerators like GPUs.
1. Introduction
Sparse generalized linear models (GLMs) are essential tools in machine learning (ML), widely applied in fields like healthcare, finance, engineering, and science. These models provide flexible framework for capturing relationships between variables while ensuring interpretability, which
School Operations Research and Information Engineering, Corell University, Ithaca, NY, USA "Jacobs Technion-Cornell Institute, Cornell Tech and Technion-IIT, New York, NY, USA. Correspondence to: Jiachang Liu <jiachang. liu @ comell.edu>> Andrea Lodi
To address these challenges, we propose scalable first- order method for efficiently calculating lower bounds within the BnB framework. We begin with a perspective refor- mulation of (1) and derive its continuous relaxation. The resulting formulation is then expressed as an unconstrained optimization problem, characterized by a convex compos-
Soroosh Shafice cshaticc@comell..cte cal748@comnel..duu
Proceedings of the 42nd International Conference Machine Learning, Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Scalable First-order Method for Certifying Optimal k-Sparse GLMs
ite objective function, which enables the application of the Fast Iterative Marinkage-Threshollee Algorithm (FISTA), well-known first-order method (Beck& & Teboulle, 2009), to compute lower bounds. The successful implementation of FISTA, however, relies on efficient computation of the proximal operator, which requires solving a second order cone program (SOCP) problem. To the best of our knowl- edge, the efficient computation of this proximal operator has not been previously addressed in the literature. Therefore, we propose customized pooled-adjacent-vvolaioo algo- rithm (PAVA) that evaluates the proximal operator exactly with log-linear time complexity. ensuring the scalability of our FISTA approach for large problem instances. A major advantage of our approach is its computational efficiency, in which instead of solving costly linear systems, it only relies on matrix-vector multiplication, which is highly amenable to GPU acceleration. This capability addresses key limita- tion of existing approaches that struggle to parallelize their computations on modern hardware.
1.2. Related Works
MIP for ML. MIP has been successfully applied in med- ical scoring systems (Ustun & Rudin, 2016; 2019; Liu et al., 2022), portfolio optimization (Bienstock. 1996; Wei et al., 2024), nonlinear identification systems (Bertsimas & Gurnee, 2023; Liu et al., 2024), decision trees (Bertsimas & Dunn, 2017; Hu et al., 2019), survival analysis (Zhang et al., 2023; Liu et al., 2025), hierarchical models (Bert- simas & Van Parys, 2020a), regression and classification models (Atamturk & Gómez, 2020; Bertsimas & Van Parys, 2020a; Bertsimas et al., 2020; Bertsimas & Van Parys, 2020b; Hazimeh & Mazumder, 2020; Xie & Deng, 2020; Atamtürk et al., 2021; Dedieu et al., 2021; Hazimeh et al., 2022; Liu et al., 2024; Guyard et al., 2024), graphical mod- els (Manzour et al.. 2021; Kucukyavuz et al., 2023), and outlier detection (Gómez, 2021: Gómez & Neto, 2023). The primary focus of these works is on obtaining high-quality feasible solutions, with only small subset addressing the certification of optimality. Our work aims to contribute to this literature. with a strong focus on enhancing the com- putational scalability of certifying optimality for solving sparse GLM problems.
To accelerate the performance of the FISTA algorithm, we introduce a restart heuristic. This leads to an empirical linear convergence rate, result not previously achieved by other first-order methods for this type of problem. Empirically, our method demonstrates substantial speedups in computing dual bounds- often by 1-2 orders of magnitude- compared to existing techniques. These improvements significantly enhance the overall efficiency of the BnB process, enabling the certification of large-scale instances of (1) that were previously intractable using commercial MIP solvers. All omitted proofs are provided in Appendix A. The experimen- tal setup is detailed in Appendix B. Additional numerical results are reported in Appendix C.
Perspective Formulations. The application of perspec- tive functions to derive convex relaxations for (1) dates back to the seminal work of Ceria & Soares (1999). Per- spective formulations have been developed for separable functions in (Günlük & Linderoth, 2010; Xie & Deng, 2020; Wei et al., 2022; Cacciola et al., 2023; Bacci et al., 2024; Shafice & Kilıng-Karzan, Karzan, 2024) and for rank-one functions in (Frangioni et al., 2020; Wei et al., 2020; 2022; Atamtürk & Gómez, 2023; Han & Gómez, 2024; Shafiee & Kiling- Karzan, 2024) under various conditions. Our work uses perspective formulation of separable functions that appear in (1) as the Tikhonov regularization function. This formu- lation provides tighter relaxations compared to the & bbaadd methods proposed in (Mhenni et al., 2020).
1.1. Contributions
The key contributions of this paper are summarized below.
We propose a FISTA-based first-order method to enhance the scalability of solving (I), with a focus on efficient lower-bound computation within the BnB framework.
Lower Bound Calculation. A key aspect of certifying optimality in MIP problems is the efficient computation of tight lower bounds. Commercial MIP solvers typically iter- atively lincarize the objective function using the celebrated outer approximation method (Kelley, 1960) (via cutting planes) and solve the resulting linear program (Schrijver, 1998; Wolsey, 2020). However, this approach often produce loose lower bounds, especially when high-quality linear cuts are not generated. Alternatively, solvers may use conic convex relaxations and solve them with the interior-point method (IPM) (Dikin, 1967; Renegar, 2001; Nesterov & Ne- mirovskii, 1994). While this approach often yields tighter lower bounds, IPM does not scale well due to its reliance on second- coder information and because- differently from the linear case - effectively warm- carm--tariigg IPMs is not possible.
The proximal operator in the FISTA method is computed using a customized PAVA that leverages hidden mathe- matical structures and enjoys log-linear time complexity, ensuring scalability for large-scale problems.
Besides achieving fast convergence rates (via a restart strategy) and low per-iteration computational complexity, our method can be easily parallelized on GPUs, some- thing not currently achievable by MIP methods.
We validate the practical efficiency of our approach on both synthetic and real- world datasets, demonstrating sub- stantial speedups in computing dual bounds and certifying optimal solutions for large-scale sparse GL.Ms
2

---

Scalable First-order Method for Certifying Optimal k-Sparse GLMs
Recent attempts are based on first-order methods, including subgradient descent (Bertsimas et al., 2020), ADMM (Liu et al., 2024), and coordinate descent (Hazimeh et al., 2022). Our work builds on this, offering faster convergence, low computational complexity, and significant GPU acceleration. We also observe that our proposed FISTA method achieves linear convergence rates empirically, a result not previously achieved by other first-order methods for this problem.
1970, $13), the optimal value of (2) remains unchanged if we replace S with el conv(S), where cl conv(S) denotes the closed convex hull of S. However, the exact descrip- tion of cl conv(S) requires exponentially many (nonlinear) constraints, which leads to the NP-hardness of (1).
We thus explore other options for a convex relaxation of (I). It turns out that tractable convex hull can be obtained if the objective function only includes the Tikhonov regularization term || using the perspective function. The perspective function of the quadratic function h(B) = _ is 1 (3, 2) 8 /2 if > 0, = 0 if 3 = 2 = O, and = oo otherwise. For simplicity, we write 8²/z instead of h (3, 2) even if 2 = 0. The following lemma provides an exact perspective formulation of the convex hull when S does not include fx3 y). This result extends (Günlük & Linderoth, 2010, Lemma 6) by incorporating sparsity constraints, while also extending (Shafiee & Krlmng-Karzan 2024, Theorem 2) to account for too box constraint on B.
GPU Acceleration. Recently, there have been some promising works on using GPUs to accelerate contin- uous optimization problems, including linear program- ming (Applegate et al., 2021: Lu et al., 2023), quadratic programming (Lu & Yang, 2023), and semidefinite program- ming (Han et al., 2024). A natural way to leverage GPUs for discrete problems is to implement greedy heuristics on GPUs (Blanchard & Tanner, 2013). Alternatively, one can use GPU-based LPs within MIP solvers, as demonstrated by De Rosa et al. (2024) for solving clustering problems. However, in (De Rosa et al., 2024), the challenge is to ap- proximate the original objective function with a potentially exponential number of cutting planes. In contrast, we de- velop a customized FISTA method that directly handles the nonlinear objective function, while the computation can be easily parallelized since it only involves matrix-veet cetrii- mul-
Lemma 2.1. The closed convex hull of the set
|1 S M,
# E (0, 13 1 z S k,
(T, 2,)
3, 2,) = 0 Vj E [p].
ΣjEløl $3 ST
tiplication. Other first-order methods, such as ADMM (Liu et al., 2024) and coordinate descent (Hazimeh et al., 2022), are unsuitable for GPUs: ADMM requires solving linear systems, while coordinate descent is inherently sequential.
is given by the set
- Maj S 5 S M2j S v lp).
(T, (T,
= E [0, 1P, 1" k.
Σjeb 03/25 ST
2. Problem Formulation
The convex hull formulation presented in Lemma 2.1 is a second-order conic set. Specifically, the epigraph of the sum of perspective functions in the last line satisfies
In this preliminary section, we introduce some backgrounds on how to obtain lower bound (which will be used for the branch-and-bound process to prune nodes) for the optimal value of Problem (I) by solving an associated convex relax- ation problem. First, note that we can cast problem (1) as
[ S Vj vje (p).
Σ 8312, ST + EE s.t. sElpl
min fr (r, B, 2) E St,
(2)
which is second order cone representable. Motivated by Lemma 2.1, we immediately see that the extended feasi- ble set S defined in (3) admits the following perspective representation
where the extended feasible set is defined as
|| S M, z E 13,, 1T k, 3, 24) 0 Vje (p] f(X3,y) + B221 ST
S=
(T,
- -22 S 33 < <ay E Ip),
S = (T,
2 2 E 13, 1 f(XB,y) + do E 33/2 ST
and [p] {1. .p) stands for the set of all integers up to p E N. Put differently, each binary variable 2j indicates whether a continuous variable 3y is zero not by requiring 33 = 0 when zj = 0 and allowing 33 to take any value when Zj 1. Mcanwhile, the objective function is linearized us- ing the epigraph reformulation technique, which allows us to interpret the optimal value of (2) as the evaluation of the support function of S at (0, 0, 1). By virtue of (Rockafellar,
Plutting in this new perspective representation into Prob- lem (2), we can reformulate (1) as follows
min
f(X + E 63/23
B.zER=
PMUP
s.t.
2E (0,1), 1T -==5k
< k.
(4)
-Mz, S 33 S Maj E (p).
3

---

Scalable First-order Method for Certifying Optimal k-Sparse GLMs
By relaxing the binary variables 2y to the interval [0, I], we obtain the following strong convex relaxation of (4)
3.1. Conjugate function g"
Recall that the conjugate of g is defined as
min 3
f(X,3, y) + A Eeip $3/2,
g'(a) = sup BERP
a T B (33).
P
S.I.
= [0, IJ,, 1" S k;
(5)
- 22 S 83 S Mz, [p)
The following lemma gives a closed-form expression for g¹, where TopSum (-) denotes the sum of the top largest elements, and Has : R + R is the Huber loss function defined
Although this is not the convex hull formulation due to the term f(X3, y). unlike in Lemma P Poom still provides a lower bound for Problem (1).
We can solve (5) using standard conic optimization solvers like Mosek and Gurobi, which rely on IPMs for solving such subproblems in the BuB framework. However, IPMs are computationally expensive and do not scale well for large datasets. Alternatively, first-order conic solvers such as SCS (O' donoghue et al., 2016), based on ADMM, can be used. While these methods are more scalable, they suffer from slow convergence rates and require solving linear sys- tems at each iteration, which can also be computationally intensive for large instances. The main goal of the paper is to introduce efficient and scalable first-order method to address these limitations.
Har(az) = 22 MMoyl 1 if |a > M
if |f S M
(8)
For notational simplicity, we use the shorthand notation HMr (a) to denote H = (Har (a1)......... H (ap)). Lemma 3.1. The conjugate of g is given by g²(a)= = (Ha(()))
(9)
This closed-form expression enables us to compute the prox- imal of g". Note that while the proximal operators of both TopSum, and H functions are known (see, for example, (Beck, 2017, Examples 6.50 & 6,54)), the conjugate func- tion g" is defined as the composition of these two functions. However, there is no general formula to derive the proximal operator of composition of two functions based on the proximal operators of the individual functions. In the next section, see how to bypass this compositional difficulty.
3. Methodology
We begin with reformulating (5) as the following uncon- strained optimization problem
min f333, y) + 2Xx g(B),
(6)
3
where the implicit function g E. + [o0) is defined as
3.2. Proximal operator of g
min
1/2
j( | 83/23
Recall that the proximal operator of g with weight parame- ter p > 0 is defined as
#ER"
((B)
S.L.
= 00, 1l, 1
(7)
M22 S B, S Mzj lp].
proX (1)) argmin
2 llla- |1-11122 + (ax).
(10)
Here, we follow the standard convention that an infeasible minimization problem is assigned value of +oo. Note that g is convex as convexity is preserved under partial mini- mization over a convex set (Rockafellar, 1970, Theorem 5,3). Furthermore, as is assumed to be Lipschitz smooth and g non-smooth, problem (6) is unconstrained optimization problem with a convex composite objective function. As such, it is amenable to be solved using the FISTA algorithm proposed in (Beck& & Teboulle, 2009).
CEER
The evaluation of proxpg" involves a minimization prob- lem that can be reformulated as convex SOCP problem. Generic solvers based IPM and ADMM require solving systems of linear equations. This results in cubic time com- plexity per iteration, making them computationally expen- sive, particularly for large-scale problems. These methods also cannot return exact solutions. The lack of exactness can affect the stability and reliability of the proximal op- erator, which is crucial for the convergence of the FISTA algorithm. Inspired by (Busing, 2022), we present Algo- rithm 1, customized pooled adjacent violators algorithm that provides an exact evaluation of proxpg' in linear time after performing simple 1D sorting step.
In the following, we first analyze the conjugate of g. We then propose an efficient numerical method to compute the proximal operator of g". This, in turn, enables us to compute the proximal operator of g. leading an efficient implementation of the FISTA algorithm. To further enhance the performance of FISTA, we present an efficient approach to solve the minimization problem (7), which guides us in developing an effective restart procedure. Finally, we conclude this section by providing efficient lower bounds for each step of the BnB framework.
Theorem 3.2. For any μ E RP, Algorithm 1 returns the exact evaluation of prox pg (μ) in O(plogp).
The proof relies on several auxiliary lemmas. We start with the following lemma, which uncovers a close connection
4

---

Scalable First-order Method for Certifying Optimal k-Sparse GLMs
Algorithm 1 Customized PAVA to solve prox pg" (#4
PAVA variant designed compute prox pg* exactly. The following lemma shows that the vector generated by Algo- rithm is an exact solution to (II). Intuitively, Algorithm I iteratively merges adjacent blocks until no isotonic con- straint violations remain, at which point the resulting vector is guaranteed to be the optimal solution to (11).
Input: vector μ. scalar multiplier p. and threshold M of the Huber loss function Has()
1: Initialize P E IR" with Pj = p if E (1,2, k} and
Pj = 0 otherwise.
2: PSort il such that |ual 2 |u2| 2 ... 2 Mppl: π is the sorting order.
Lemma 3.4. The vector i in Algorithm I solves (11) ex- actly.
3: μ. = SpecialSort((1)
4: DSTEP 1: Initialize pool of blocks with start and end indices; each block initially has length equal to 1 5: J = =[[1,1], [2, 2).,, [p. pl}
Finally, the merging process in Algorithm I can be executed efficiently. Intuitively, each element of ş is visited at most twice; once during its initial processing and once when it is included in merged block. This ensures that the process achieves linear time complexity.
6: DSTEP 2: Initialize by in each block by ignoring the isotonic constraints
7: for - 1,,,..... .......dd
Lemma 3.5. The merging step (lines 11-14) in Algorithm I can be performed in linear time complexity O(p).
8: D3 = argmin, 2v 11,1) + ()H))())
9: end for
10: >SEEP 3: Whenever there is an isotonic constraint violation between two adjacent blocks, merge the two blocks by setting all values to be the minimizer of the objective function restricted to this merged block; use Algorithm 4 in Appendix A
Armed with these lemmas, one can easily prove Theo- rem 3.2. Details are provided in Appendix A.
3.3. FISTA algorithm with restart
A critical computational step in FISTA is the efficient evalu- ation of the proximal operator of g. By the extended Moreau decomposition theorem (Beck, 2017, Theorem 6.45), for any weight parameter > and any point μ E TRP, the proximal operators of and 9 satisfies
11: while 3aa, a2); (a2 + 1, 1, E I s.t. Dai < Das do
12: J = J |2]] {la2 + 1, a,aa]} U (aa aal}
13:
3aas = argminn Σ [ 10, + py (v)]
==
14: end while
proxp-1g μa - p prox pg* (pn).
(12)
15: >Return i with the inverse sorting order
16: Return sgn(p) @ (D)
Hence, together with Theorem 3.2, we can compute exactly prox--ig using Algorithm 1 with log-linear time complex- ity. This enables an efficient implementation of the FISTA algorithm. We note that when M = oo, the implicit function 9 closely related the k-support norm. For this special case, Argyriou et al. (2012) and McDonald al. (2016) have developed efficient algorithms that compute the proxi- mal operator of g in the primal space with log-linear time complexity. Our work extends these results by enabling the incorporation of big-M constraints through a dual analysis.
between the proximal operator of g" and the generalized isotonic regression problems.
Lemma 3.3. For any μ E RP, we have
prox pg* (p) = sgn(H) 3 v",
where ② denotes the Hadamard (element-wise) product, v² is the unique solution of the following optimization problem
Besides, the vanilla FISTA algorithm is prone to oscilla- tory behavior, which results in a sub-linear convergence rate of (11/T2) after T iterations. In the following, we further accelerate the empirical convergence performance of the FISTA algorithm by incorporating a simple restart strategy based on the function value, originally proposed in (O' d'ooogg & Candes, 2015). In simple terms, the restart strategy operates as follows: if the objective func- tion increases during the iterative process, the momentum coefficient reset to its initial value. The effectiveness of the restart strategy hinges the efficient computation of the loss function. This task essentially reduces to evalu- ating the implicit function defined in (7), which would involve solving a SOCP problem. However, the value of g can be computed efficiently by leveraging the majorization technique (Kim et al., 2022), as shown in Algorithm 2.
min VER" 2 Σsem 1p,1))² + Eçer Ha((s)
(11)
s.t.
13 if if ml Vj E lp),
and I the set of indices of the top largest elements of [nyl.3 E (p).
Problem (11) replaces the TopSun, in (10) from the conju- gate function 9" (as shown in Lemma 3) with linear con- straints. While this may appear computationally complex, it actually converts the problem into an instance of isotonic regression (Best & Chakravarti, 1990). Such problems can be solved exactly in linear time after performing a simple sorting step. The procedure is known as PAVA (Busing, 2022). Specifically, Algorithm implements a customized
5

---

Scalable First-order Method for Certifying Optimal k-Sparse GLMs
Algorithm 2 Algorithm to compute g(33)
tive algorithm, FISTA yields only an approximate solution 3 to (5). Consequently, while we can calculate the objec-
Input: vector 3 E R" from Step 2 Line 8 in Algorithm 3.
tive function for 3 efficiently, this value is not necessarily lower bound of the original problem-on only the optimal value of the relaxed problem (5) serves as a guaranteed lower bound. To get a safe lower bound, we rely on the weak duality theorem, in which for any proper, lower semi- continuous, and convex functions F R + E U {oo} and G: R R U {oo}, we have
1: Initialize: & = 0 E R²
2: Sort B partially such that
13:1 2 |82 2 ... 2 18xl 2 max (188))
k+I.....
3: 8 = Σ3
4: DCompute the majorization vector ub; see Appendix 3.6 for its definition and connection with (3
5: for 1, ... do
inf F(X3)+ + G(B) > BER
sup CER"
- --
6: 8 = =//(k-
1j+1 +1
7: if 2 |3,| then d. = 3; break else 3/ = 13,
2 (-F - () V E R.
8: == |,31|
9: end for
where F and G denote the conjugates of F and G. re- spectively, while the second inequality follows from the definition of the supremum operator. Letting F(X (3 3(X v). G(B) = (x) and Ç V X (),, where & is the output of the FISTA Algorithm 3, we arrive at the safe lower bound
10: return
2 Σー= 133
Theorem 3.6. For any B E RP, Algorithm 2 computes the exact value of g(3), defined (7), in Opp p plog k). Theorem 3.6 guarantees that Algorithm 2 can efficiently compute the value of g(3). which is crucial for our value- based restart strategy to be effective in practice. Empirically, we observe that the function bbsed restart strategy can accelerate FISTA from the sub-linear convergence rate of O(1/T3) to a linear convergence rate in many empirical results. To the best of our knowledge, this is the first linear convergence result of using a first-order method in the MIP context when calculating the lower bounds in the BnB tree. The FISTA algorithm is summarized in Algorithm 3.
MMP -F -
(13)
where the inequality follows from the relaxation bound MIP 2 Paov and the weak duality theorem. For conve- nience, we provide a list of F*).) for different GLM loss functions in D.I. The readers are also referred to D.2, where we derive the safe lower bound for the linear regression problem with pegseeperrpeetive relaxation an example.
4. Experiments
3.4. Safe Lower Bounds for GLMs
We evaluate our proposed methods using both synthetic and real-world datasets to address three key empirical questions:
We conclude this section by commenting how to use Algorithm 3 in the BnB tree to prune nodes. As an itera-
How fast is our customized PAVA algorithm in evaluating prox, compared to existing solvers?
How fast is our proposed FISTA method in calculating the lower bounds compared to existing solvers?
Algorithm 3 Main algorithm to solve problem (5) Input: number of iterations T, coefficient A for the & regularization, and step size T. (Lipschitz-continitty parameter of VF(3))
How fast is our customized BnB algorithm compared to existing solvers?
1: Initialize: g 3 - 0. 0=
We implement our algorithms python. For baselines, we compare with the following -tate-of-the-art commercial and open-source SOCP solvers: Gurobi (Gurobi Optimiza- tion, LLC, 2025), MOSEK (ApS, 2025), SCS (O' Donogguue et al., 2023), and Clarabel (Goulart & Chen, 2024), with the python package cvxpy (Diamond & Boyd, 2016) as the interface to these solvers.
2: Let: p === L/(222), f(B
3: for 1,,,3..... T do
4: >Step 1: momentum acceleration
5:
7'= ='=B +++++ 23' - 31-1)
6:
Step 2: proximal gradient descent; use Algorithm
7:
y²= ==1 7' F(7)(r*)
8:
g++1- - pp-pp-a prox pg- (pr')
4.1. How Fast Can We Evaluate prox,-i,(-)?
9:
Step 3: restart; use Algorithm 2
In this subsection, we demonstrate the computational ef- ficiency of using our PAVA algorithm for evaluating the proximal operators. We conduct the comparisons in two ways - evaluating both a) the proximal operator of the orig- inal function and b) the proximal operator of its conjugate
10:
C++l= = f( yy + (x2g(8*
11:
if [* 2 L then & = 1 else x 0 +
12: end for
13: return +1+1
6

---

Scalable First-order Method for Certifying Optimal k-Sparse GLMs
Proximal Operator of
Proximal Operator of g
Solver Time Elapsed
Solver Time Elapsed
10
10a
10-1
10"1
10-3
10-3
0
20000 40000 60000 80000 100000
0
20000
40000 60000 80000 100000
P
gurobi
mosek
SCS
clarabel
ours
Figure 1. Running time comparison 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.
Linear Loss
Logistic Loss
10a
10a
Solver Time Elapsed
Solver Time Elapsed
102
102
10
10
109
10
10 1
10-1
2500
5000
7500 10000
12500 15600
2500
5000
7500 10000 12500 15600
gurobi
masek
+ scs
clarabel
ours
Figure 2. Running time comparison of solving Problem (5), the perspective relaxation the original MIP problem. We set M 2.0,
A 1.0, and rt-to-p ratio to be Gurobi cannot solve the cardinality constrained logistic regression problem.
Table 1. 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).
to the approximate solutions returned by the baselines. This combination of precision and efficiency constitutes a crit- ical advantage for our first-order optimization framework over generic conic programming solvers, demonstrated in subsequent sections.
P
1k
2k
4k
8k
16k
0.19
0.48
1.54
4.80
19.52
ours CPU
(0.01)
(0.05)
(0.21)
(0.57)
(1.27)
4.2. How Fast Can We Calculate the Lower Bound?
ours GPU
0.29
0.19
0.26
0.59
2.09
We next benchmark the computational speed and scalability of our method against the state-of-the-arr solvers (Gurobi, MOSEK, SCS, and Clarabel) for solving the perspective relaxation of the original MIP problem. Evaluations per- formed both linear and logistic regression tasks. Experi- mental configurations are detailed in Appendix B.2. Addi- tional perturbation studies, such as on the sample-to-feature (n-to-p) ratio, box constraint M, and la regularization co- efficient Ag. are provided Appendix C.I. All solvers are terminated upon achieving an optimality gap tolerance of 10-6 exceeding a runtime limit of 1800 seconds.
(0.04)
(0.02)
(0.02)
(0.08)
(0.11)
0" Detailed experimental configurations, including parame- ter specifications and synthetic data generation process, are provided in Appendix B.I.
The results shown in Figure highlight the superiority of our method. Our algorithm achieves a computational speedup of approximately two orders of magnitude compared to con- ventional SOCP solvers. This performance gain is largely due to our customized PAVA implementation in Algorithm 1. For instance, in high-dimensional settings (p 10"), base- line methods require several seconds to minutes to evaluate the proximal operators, whereas our approach completes the same task in 0.01 seconds. Additionally, our method guaran- tees exact solutions to the optimization problem, in contrast
The results, shown in Figure 2, demonstrates that our method outperforms the fastest conic solver (MOSEK) by over one order of magnitude. For the largest tested instances (n = 16000 and p = 16000), our approach attains the target

---

Scalable First-order Method for Certifying Optimal k-Sparse GLMs
Iteration vs. Optimality Gap
Iteration vs. Optimality Gap
Optima lity Gap
(primal loss optimal loss)/aptimall loss
(optimal loss dual loss)/oplimal loss
(primal loss dual loss//primal loss
Optimality Gap
10
Optimality Gap
10
Optimality Gap
10a
10-3
10
10-3
10-3
100-
10-5
10-3
200 Iteration
300
100
200
400
500
200
300
Iteration
Iteration
PGD
FISTA
Restarted FISTA
Figure 3. Empirical convergence of restarted FISTA (compared with PGD, the proximal gradient method, and FISTA) solving the perspective relaxation in Problem (5) with the logistic loss, n 16000, p = 16000, k = 10, p = 0.5, A2 = 1.0, and M = 2.0.
Table 2. Certifying optimality on large-scale and real-world datasets.
ours
Gurobi
MOSEK
time (8)
opt. gap (%)
time (s)
opt. (%)
time (s)
opt. gup (%)
synthetic 10, 2)
79
0.0
1800
1915
Linear Regression
(n=16k, p=16k, seed=0)
Cancer Drug Response (l: 5, =
0.0
1800
188
0.0
(n=822, p=2300)
Synthetic (k: = 10, 2)
626
0.0
N/A
N/A
2446
Logistic Regression
(n=16k, p=16k, sced=0)
DOROTHEA (A'
10, . 2)
230
0.0
N/A
N/A
1814
0.63
(n=2300, p=89989)
tolerance (10-6) in under 100 seconds across regression and classification datasets, whereas most baselines fail to converge within the 1800-second threshold.
commericial MIP solvers. Integrating our lower-hound com- putation into a minimalist branch-and-bound (BnB) frame- work, we prioritize node pruning via lower bound calcula- tions while intentionally omitting advanced MIP heuristics (e.g., cutting planes, presolve routines) to evaluate the im- pact of our method. Experimental configurations, including dataset descriptions and BnB implementation details, are provided in Appendix B.3. We benchmark our approach against Gurobi and MOSEK, reporting both runtime and fi- nal optimality gaps. Note that the real-world datasets, we have used small k's, which are selected based on 5-fold cross validation (see Appendix B.3). Small A are sufficient for accurate prediction and can help avoid overfitting. They also improve interpretability of the model.
There are two factors driving this speedup. First, our ef- ficient proximal operator evaluation reduces per-iteration complexity. Second, our efficient method to compute g(3) (in Algorithm 2) exactly enables integration of the value- based restart technique within FISTA, significantly improv- ing convergence. Figure illustrate this enhancement: while PGD (the proximal gradient descent method, also known as ISTA) and FISTA exhibit slow convergence rates, FISTA with restarts achieves fast linear convergence on both dual loss and primal-dual gap metrics. To the best of our knowl- edge, this marks the first empirical demonstration of linear convergence for first-order method applied to solving the convex relaxation of this MIP class. Finally, our method permits GPU acceleration because most computation- ally intensive component is matrix-vector multiplications. As shown in Table 1, GPU implementation reduces runtime by an additional order of magnitude on high-dimensional instances.
Results in Table show that our method certifies optimality for of the four tested datasets around 1 minute, the third around 10 minutes, and the fourth around 4 minutes. In con- trast, Gurobi and MOSEK either exceed the time limit (1800 seconds) during the presolve stage or require significantly longer runtimes to achieve zero or small gaps. Crucially, this efficiency stems from our efficient lower-bound compu- tations and dynamic early termination criteria. Specifically, we avoid waiting for full convergence by leveraging two key rules: (1) if the primal loss falls below the incumbent solution's 's loss, we terminate early and proceed to branching; (2) if the dual loss exceeds the incumbent's 's loss, we halt
4.3. How Fast Can We Certify Optimality?
Finally, we demonstrate how our method's ability to com- pute tight lower bounds enables efficient optimality certifica- tion for large-scale datasets, outperforming state-of-the-art
8

---

Scalable First-order Method for Certifying Optimal k-Sparse GLMs
computation and prune the node immediately. This adaptive approach eliminates unnecessary iterations while ensuring we prune the search space effectively.
Argyriou, A., Foygel, R., and Srebro, N. Sparse prediction with the k-support norm. Advances in Neural Information Processing Systems, 25, 2012.
Asuncion, A. and Newman, D. UCI machine learning repos- itory. 2007.
5. Conclusion
introduce a first-order proximal algorithm to solve the perspective relaxation of aurdinality-constrimme GLM prob- lems. By leveraging the problem's unique mathematical structure, we design customized PAVA efficiently eval- uate the proximal operator, ensuring scalability to high- dimensional settings. Further acceleration is achieved through an efficient value-based restart strategy and compat- ibility with GPUs, which collectively enhance convergence rates and computational speed. Extensive empirical results demonstrate that our method outperforms state-of-the-art solvers by 1-2 orders of magnitude, establishing it as practical, pigh-performance component for integration into next-generation MIP solvers.
Atamturk, A. and Gómez, A. Safe screening rules for to- regression from perspective relaxations. In International Conference on Machine Learning, pp. 421-430, 2020. 2020.
Atamtürk, A. and Gómez, A. Supermodularity and valid
inequalities for quadratic optimization with indicators. Mathematical Programming, 201:295-338, 2023.
Atamtürk, A., Gómez, A., and Han, Sparse and smooth signal estimation: Convexification of fo-formulations. Journal of Machine Learning Research, 22:52-1, 2021.
Bacci, T. Frangioni, A., Gentile, C., and Tavlaridis- Gyparakis, K. New mixed-integer nonlinear program- ming formulations for the unit commitment problems with ramping constraints. Operations Research, 72(5): 2153-2167, 2024.
Code Availability
Implementations discussed in this paper are available https:: hittps://githu.. http://gttttt.ccooc aachanggl //tku//kK
Beck, A. First-Order Methods in Optimization. SIAM, 2017.
Acknowledgments
Beck, A. and Teboulle, M. fast iterative shrinkage- thresholding algorithm for linear inverse problems. SIAM Journal Imaging Sciences, 2(1):183-2022 2009.
This work used the Delta system at the National Cen- ter for Supercomputing Applications through allocation CIS250029 from the Advanced Cyberinfrastructure Coordi- nation Ecosystem: Services Support (ACCESS) program, which is supported by National Science Foundation grants #2138259, #2138286, #2138307, #2137603, and #2138296.
Bertsimas, D. and Dunn, J. Optimal classification trees. Machine Learning, 106: 166:1039-1000 2017.
Bertsimas, D. and Gurnee, W. Learning sparse nonlinear dynamics via mixed-integer optimization. Nonlinear Dy- namics, 111(7):6585-660) 2023.
Impact Statement
This paper presents work whose goal advance the field of Machine Learning. There are many potential societal consequences of our work, none which feel must be specifically highlighted here.
Bertsimas, D. and Van Parys, B. Sparse hierarchical re- gression with polynomials. Machine Learning, 109(5): 973-997, 2020a.
Bertsimas, and Van Parys, B. Sparse nigh-dimensional re- gression. The Annals of Statistics, 48(1):300-323, 2020b.
References
Ahuja, R. K. and Orlin, J. B. A fast scaling algorithm for minimizing separable convex functions subject to chain constraints. Operations Research, 4985)7784-789 2001.
Bertsimas, D., Pauphilet, J., and Van Parys, B. Sparse regression. Statistical Science, 35(4):555-578, 2020.
Applegate, D., Díaz, M., Hinder, O., Lu, H., Lubin, M., O'Donoghue, B., and Schudy, W. Practical large-scale linear programming using primal-dual hybrid gradient. In Advances in Neural Information Processing Systems, 2021.
Best, M. J. and Chakravarti, N. Active set algorithms for isotonic regression; unifying framework. Mathematical Programming, 47(1):425-433 1990.
Best, M. J., Chakravarti, N., and Ubhaya, V. A. Minimizing separable convex functions subject to simple chain con- straints. SIAM Journal on Optimization, 10(3):658-66) 672, 2000.
pp.
ApS, M. The MOSEK optimization toolbox for MATLAB manual. Version 11.0.4 2025.

---

Scalable First-order Method for Certifying Optimal k-Sparse GLMs
Bienstock, D. Computational study of family of mixed- integer quadratic programming problems. Mathematical Programming,
Günlük, O. and Linderoth, J. Perspective reformulations of mixed integer nonlinear programs with indicator vari- ables. Mathematical Programming, 124()):183-205 2010.
1996.
Blanchard, J. D. and Tanner, J. GPU accelerated greedy algorithms for compressed sensing. Mathematical Pro- gramming Computation. 5(3):277-34,, 2013.
Gurobi Optimization, LLC. Gurobi Optimizer Reference Manual, 2025.
Guyard, T., Herzet, C., Elvira, C., and Arslan, A.-N. new branch-and and-bo-dd- pruning framework for fn-regularized problems. In International Conference on Machine Learn- ing, 2024.
Boerner, T. J., Deems, S., Furlani, T. R., Knuth, S. L., and
Towns, J. Access: Advancing innovation: NSF's ad- vanced cyberinfrastructure coordination ecosystem: Ser- vices & support. In Practice and Experience in Advanced Research Computing, pp. 173-176, 2023.
Han, Q. Lin, Z., Liu, H., Chen, C., Deng, Q. Ge. D., and Ye, Y. Accelerating low-rank factorization-baser semidefinite programming algorithms on GPU, arXiv:2407.15049, 2024.
Busing, F. M. Monotone regression: A simple and fast (nn) PAVA implementation. Journal of Statistical Software, 102:1-25, 2022.
Han, S. and Gómez, A. Compact extended formulations for low-rank functions with indicator variables. Mathematics of Operations Research (in press), 2024.
Cacciola, M., Frangioni, A., Li, X., and Lodi, A. Deep neural networks pruning via the structured perspective regularization. SIAM Journal on Mathematics of Data Science, 5(4):1051-107) 2023.
Hazimeh, H. and Mazumder, R. Fast best subset selec- tion: Coordinate descent and local combinatorial opti- mization algorithms. Operations Research, 68(5):1517- 1537, 2020.
Ceria, S. and Soares, J. Convex programming for disjunctive convex optimization. Mathematical Programming, 86(3): 595-614, 1999.
Hazimeh, H., Mazumder, R., and Saab, A. Sparse regres- sion at scale: Branch-and-boumd rooted in first-order opti- mization. Mathematical Programming. 196(():3477333) 2022.
De Rosa, A., Khajavirad, A., and Wang, Y. On the power of linear programming for K-means clustering. arXiv:2402.01061 2024.
Hu, X., Rudin, C.,, and Seltzer, M. Optimal sparse decision trees. In Advances in Neural Information Processing Systems, pp. 7267-7275, 2019.
Dedieu, A., Hazimeh, H., and Mazumder, R. Learning sparse classifiers: Continuous and mixed integer opti- mization perspectives. Journal of Machine Learning Research, 22(1):6008-60) 6054, 2021.
Kelley, Jr, J. E. The cutting-plane method for solving con- vex programs. Journal of the society for Industrial and Applied Mathematics, 8(4):703-712, 1960.
Diamond, S. and Boyd, S. CVXPY: A Python-embedded modeling language for convex optimization. Journal of Machine Learning Research, 17(83):1 2016.
Kim, J., Tawarmalani, M., and Richard, J.-P. P. Convexi- fication of cormuutation-nnvarian sets and an application sparse principal component analysis. Mathematics of Operations Research, 47(4):2547-2584, 2022.
Dikin, I. Iterative solution of problems of linear and quadratic programming. Doklady Akademii Nauk, 174 (4):747-748, 1967.
Kucukyavuz, S., Shojaie, A., Manzour, H., Wei, I.., and Wu, H.-H. Consistent second-order conic integer program- ming for learning bayesian networks. Journal of Machine Learning Research, 24(322):1-38, 2023.
Frangioni, A., Gentile, C., and Hungerford, Decom- positions of semidefinite matrices and the perspective reformulation of nonseparable quadratic programs. Math- ematics of Operations Research, 45(1):15-33, 2020.
Lei, T., Bai, J., Brahma, S., Ainslie, J., Lee, K., Zhou, Y., Du, N., Zhao, V., Wu, Y., Li, B.. al. Conditional adapters: Parameter-efficient transfer learning with fast inference. In Advances in Neural Information Processing Systems, pp. 8152 8172, 2023.
Gómez, A. Outlier detection in time series via nixed-integer conic quadratic optimization. SIAM Journal on Optimiza- tion, 3177)1187-11955, 2021.
Gómez, A. and Neto, J. Outlier detection in regression: Conic quadratic formulations. arXin:2307.059755 2023.
Liu, J., Zhong, C., Li, B., Seltzer, M., and Rudin, C. Faster- risk: Fast and accurate interpretable risk scores. In Ad- vances in Neural Information Processing Systems, pp. 17760-17773, 2022.
Goulart, P. J. and Chen, Y. Clarabel: An interior-point solver for conic programs with quadratic objectives, 2024.
10

---

Scalable First-order Method for Certifying Optimal k-Sparse GLMs
Liu, J., Rosen, S., Zhong, C., and Rudin, C. Okridge: Scalable optimal k-sparse ridge regression. In Advances in Neural Information Processing Systems, pp. 41076- 41258, 2024.
Park, Y. W. and Klabjan, D. Subset selection for multiple
linear regression via optimization. Journal of Global Optimization, 77(3):543-574, 2020.
Renegar, J. A Mathematical View of Interior-Poine Methods
Liu, J., Zhang, R., and Rudin, C. Fastsurvival: Hidden com- putational blessings in training cox proportional hazards models. Advances in Neural Information Processing Systems (in press), 2025.
in Convex Optimization. SIAM, 2001.
Rockafellar, R. T. Convex Analysis. Princeton University
Press, 1970.
Schrijver, A. Theory of Linear and Integer Programming. Wiley, 1998.
Liu, Q., Hu, Z., Jiang, R., and Zhou, M. DeepCDR: A hybrid graph convolutional network for predicting cancer drug response. Bioinformatics, 36(Supplement 2):1911- 1918, 2020.
Shafice, S. and Kihnç- Karzan, Constrained optimization of rank-one functions with indicator variables. Mathe- matical Programming, pp. 1-47, 2024.
Lu, H. and Yang, J. A practical and optimal first-order method for large-scale convex quadratic programming. arXiv:2311.07710. 2023.
Ustun, B. and Rudin, C. Supersparse linear integer models for optimized medical scoring systems. Machine Learn- ing, 102:349-391, 2016.
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. mxiv:2331.148322 2023.
Ustun, B. and Rudin, C. Learning optimized risk scores. Journal of Machine Learning Research, 20(150):1-75, 2019.
Manzour, H., Küçükyavuz, S., Wu, H.-H., and Shojaie, A. Integer programming for learning directed acyclic graphs from continuous data. INFORMS Journal on Optimization, 3(1):46-73, 2021.
Wei, L., Gómez, A., and Küçükyavuz, S. On the convexi- fication of constrained quadratic optimization problems with indicator variables. International Conference on Integer Programming and Combinatorial Optimization, pp. 433-447, 2020.
McDonald, A. M., Pontil, M., and Stamos, D. New perspec- tives on k-support and cluster norms. Journal of Machine Learning Research, 17(155):1-38, 2016.
Wei, L., Gómez, A., and Küçükyavuz, S. Ideal formula- tions for constrained convex optimization problems with indicator variables. Mathematical Programming, 192(1):
Mhenni, R. B., Bourguignon, S., Mongeau, M., Ninin, J., and Carfantan, H. Sparse branch and bound for exact optimization of no-morm penalized least squares. In- ternational Conference on Acoustics, Speech and Signal Processing, pp. 5735-5739, 2020.
57-88, 2022.
Wei, I.., Atamtürk, A., Gómez, A., and Küçükyavuz, S. On the convex hull of convex quadratic optimization prob- lems with indicators. Mathematicul Programming, 204: 703-737, 2024.
Natarajan, B. K. Sparse approximate solutions to linear systems. SIAM Journal Computing, 24(2):227-234, 1995.
Wolsey, L. A. Integer Progrumming. Wiley, 2020.
Xie, W. and Deng. X. Scalable algorithms for the sparse ridge regression. SIAM Journal on Optimization, 30(4): 3359-3386, 2020.
Nesterov, Y. and Nemirovskii, A. Interior-Poim Polynomial Algorithms in Convex Programming. SIAM, 1994.
Xic, Y., Dai, H., Chen, M.. Dai, B., Zhao, T. Zha, H.. Wei, W., and Pfister, T. Differentiable top-k with optimal transport. Advances in neural information processing systems, 00:20500220333 2020.
O' O'Donoghue, B., Chu, E.,, Parikh, N., and Boyd, S. SCS: Splitting conic solver, version 3.2.4, 2023.
O' d'donoguue, B. and Candes, E. Adaptive restart for accel- erated gradient schemes. Foundations of Computational Mathematics, 15:715-732, 2015.
Zhang, R., Xin, R., Seltzer, M., and Rudin, C. Optimal sparse regression trees. In AAAI Conference on Artificial Intelligence, pp. 11270- 11279, 2023.
O' d'oongguue, B., Chu, E., Parikh, N., and Boyd, S. Conic optimization via operator splitting and homogeneous self- dual embedding. Journal of Optimization Theory and Applications, 169:1042-1068. 2016.

---

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
(T, 11 - M, E {0,1}, 1 s 8,11 2z) = vj E [p], Σ, Eebll 03 5
is given by the set
{ (T, M2j S B3 s M2, Vj E pp), = E [0, 1P, 1 k, Eje(n B /
|M
Proof. Let T represent the first set mentioned the statement of the lemma. Using the definition of the perspective function and applying the big-M formulation technique, we have
T= (T, M2, s 3, S M2j Vj E lp. (0, 1p, 1 EjEw
- $
As the epigraph of perspective function constitutes cone (Shafiee & Kiling-Karzan, 2024, Lemma 1 & 2), we may write T Proj(r.A..)(T) where
7 T= 7
tj BJ]
(T, z)
= 2 10, 1 k, A,
+ Bjzg E Ky Io]
admits mixed-binary conic representation with
1
0
0
0
O
Ky =
R2
Vj (p]
0
M
-1
M
Here, L, E 3 denotes the rotated second order cone, that is, L4 = {(t, 3, E ER x : st2). Thus, using (Shafiee Kiling-Karzan, 2024, Lemma 4), the set satisfies all the requirements of (Shafiee & Kiling-Karzan, 2024, Theorem 1). and therefore, its continuous relaxation gives the closed convex hull of T. that is,
{ (T)))
t B
cl conv
r, z [0. 1]" 1 A,
+ Byz Bosje,
The prove concludes by applying Fourier-Motzkin elimination method project out the variable t.
A.2. Proof of Lemma 3.1
Lemma 3.1. The conjugate of g is given by
g" (a) = TopSum, (Hv(a)).
Proof. Using the definition of the implicit function g in (7), we have
max a T 1 Σ 83/2z
ga))) = s.t. 3 RD, = [0. 10, 1 < h.
(14)
-Mz, < 3, < M Vje lp]
12

---

Scalable First-order Method for Certifying Optimal k-Sparse GLMs
For any fixed feasible 2, the maximization problem over /3 is simple constrained quadratic problem, that can be solved analytically by the vector 3 whose j'th element is given by / = sgn(ay) min(lay), M)z.. Substituting the optimizer [3', the objective function of the maximization problem in (14) simplifies
1 Σ
(sgm (az) min (layl, 22)
a T
3r
12
Σ sElp]
aj sgn(a;) min([a)), M)),
2
22j
jElp1
Σ
laylmin(lay]: M) 2 (aaaaa.MM) zg
1
jE(p]
= Σ
Har(g)ajs,
jC(M)
where the second equality holds as is binary vector, and the last equality follows from the definition of the Huber loss function:
122²
if | S
M
HM
2 1
if 2
We thus arrive at
g'(a)= max =E[0,1]²
(Eeebl 1T TopSum, (Ha(aa))
This completes the proof.
Note that recent works (Xie al., 2020; Lei al., 2023) have discussed the smooth approximation of the TopSum, (-) However, such smooth approximation is not suitable for this work. The efféctiivnesss of the proximal algorithm relies on the exact evaluation of the proximal operator, which we will talk about next. Replacing TopSum, (-) would lead to solving different problem rather than the proximal operator evaluation. This will not help us to use FISTA to solve the perspective relaxation. As a result, this approach would not guarantee valid lower bounds necessary for optimality certification.
A.3. Proof of Lemma 3.3
Lemma 3.3. For any μ E IRP, we have
prox pg- (10) = sgn(()) @ v², where 2 denotes the Hadamard (element-wii wiset product, v" is the unique solution of the following optimization problem
min MEXP
2 &EjEm (v; - |m51)² E3Eo Hv)
(15)
3 14 2 mel Vj, E lp),
and I is the of indices of the top k largest elements of [ngli (p).
Proof. For simplicity, let ax prox (1t), that is,
a" argmin
~l1aa 11l2 + (a))
(16)
&ER
We first show that sgn(a") sgn(μt) (step 1) and then establish that for every j, [p] with lngl 2 |pr] have lazl 2 loii (step 2). We then conclude the proof using these observations.
property through a proof by contradiction. For the sake of contradiction, suppose that there exists some E [p] such that sgn(a²) o sgn(u1)) Hence, we can construct a new a' by flipping the sign aç, i.e. a; aa and keeping the rest of the elements the same as ax. Now under the assumption that sgn(a;) * sgn(H;), we have |a; - H41 > ||aÿ| - |n,|| = lay - My so the j-th term in the first summation of the objective function will decrease while everything else remains the same. This leads to a smaller objective value for a' than a" which contradicts the optimality of a Thus, the claim follows.
Step 1. We prove the sign-preservin
13

---

Scalable First-order Method for Certifying Optimal k-Sparse GLMs
Step 2. We prove the relative prgsevude-preserviing property through a proof by contradiction. For the sake of construct new a' by swapping oj and aj, i.e., az = of and af = aj. and keeping the rest of the elements the same as a². Under the assumption that |را| 2 |ue| but |a31 |aï], we have |až رy + |aj pi| > laj - My + |az - -ا la - Myl + |a/ pul, so the sum of the j-th and /-th terms in the first summation the objective function will decrease while everything else remains the same. This leads to a smaller objective value for a' than a², which contradicts the optimality of a*. Thus, the claim follows.
contradiction, suppose that there exists some j, [p] such that [ayl 2 |ma| but | < |of]- Then,
Using these two observations, we are ready to prove that = sgn(p) ② v². We first reparametrize the minimization problem (16) by substinuting the decision variable with new variable E Ry satisfying sgg( v. By the sign-preserving property in step 1, it is easy to show the equivalence between the optimization problem in (16) and the following optimization problem
min MER
1 EEEEp - 1ayl)² + TopSum, (Hv(v)).
By the relative magnitude-pessrringg property in step 2, we can further set the equivalence between the minimization problem in (16) and the following optimization problem
min VER" 2 |Ej 1ms122 + Ejes Har(yy),
s.t. 14 2N 1 if M4 2 lmd.
Lastly, the nonnegative constraint on can be removed as the second summation term in the objective function implies that vy 2 Thus, we have shown that any feasible point or in the minimization problem (16) can be reconstructed by any feasible point in the minimization problem in the statement of lemma, while maintaining the same objective value. Hence, we may conclude that a sgn( @ n², required.
A.4. Proof of Lemma 3.4
Lemma 3.4. The vector ò in Algorithm solves (11) exactly.
Proof. The minimization problem (11) is an instance of generalized isotonic regression problem taking the form
min
hy(2))
s.t. 1/1 2 v2
2 VJJ
(17)
where hy(v)= ] (5-w11)²2 (,) -
J and Pj = 0 otherwise, and the set J is the set of indices of top k largest elements of [pyl. as defined in the statement of Lemma 3.3. Thanks to (Best et al., 2000: Ahuja & Orlin, 2001), the
optimizer of (17) satisfies two key properties:
Property 1: Optimal solution for a merged block is single-valued. Suppose we have two adjacent blocks [a1, ag] and Jag 1, ag such that the optimal solution of each block is single-valued, that is, the minimization problems
min Maring s.t.
Ej hz(1z)
min Mayting s.t. Mast1 2 2 Vag
Ej hT(w,)
and
Va, 2 2 Vuz
are solved by = ma and vis+ Lses with 2/, optimal solution for the merged block (a1; 03) is single-valued that is, the minimization problem
vä, and viz+1 =... vž respectively. If v2 S vost1 then the
min Mayis
Ejan, hT(2))
s.t.
Vaj 2 2 Vag
is solved by ti Vases
with V "s"
14

---

Scalable First-order Method for Certifying Optimal k-Sparse GLMs
No isotonic constraint violation between single-valued blocks implies the solution is optimal. Suppose that we have s blocks [2); (d2 (. (as 1,0+11] (with d1 1 and 0,+1 p) such that the optimal solution for each block is single-valued, that is, Mat+1 optimal solution to (17).
Property 2:
Va Vai- for all ( E [s]. Then, if Poz 2 Pog+1 2 ... .a. then i is the
Using these two properties, it is now easy to see why Algorithm I returns the optimal solution. We start by constructing blocks which have length 1. The initial value restricted to each block is optimal. Then, we iteratively merge adjacent blocks and update the values of 23 's whenever there is a violation of the isotonic constraint. By the first property, the optimal solution for the merged block is single-valued. Therefore, we can compute the optimal solution for the merged block by solving univariate optimization problem. We keep merging blocks until there is no isotonic constraint violation. When this happens, by construction, the solution for each block is single-valued and optimal. By the second property, the final vector the optimal solution to (17), required.
A.5. Proof of Lemma 3.5
Lemma 3.5. The merging step (lines 11-14) in Algorithm can be performed in linear time complexity Opp).
Proof. A detailed implementation of line 11-14 (Step 3) of the PAVA Algorithm that achieves a linear time complexity is presented in Algorithm 4. In the following, we first show that Algorithm 4 accomplishes the objective in lines 1 11-1 of Algorithm 1. We then establish that Algorithm 4 runs in linear time complexity.
To prove the first claim, we show that the parameters Po, So, and Vs amount to
Pb = EjeB(s) Pj. S. = EjeBon 1usb v = prox E Hy (1u5))
for each block index b, where B(b) denoting the set of indices in the b'th block. easy to verify that Algorithm recursively computes P, and Sy. Thus, will focus v/- Note that the computation of the proximal operator in v/ is reduced to solving univariate optimization problem for each b and satisfies
Σ jEB(b)
(210- - 11,1)² + P3 HM (u)
V/ = argmin VER
argmin
Σ
0² ا|uاl + p3 Har(v)
jEB(b)
argmin (5/2
So
|,,1
Po No
Har(v))
argiin
1 2 v
(1-
Sb Nb
2
+
Po No
Har(e)) prox Ps No
Sb No
No
Har
Thus, Algorithm merges two adjacent blocks if the isotonic violation persists, and the output of the proximal operator is the minimizer of the univariate function in the merged block. This is exactly the same as the objective in lines 11-14 of Algorithm 1. Hence, the first claim follows.
To show that the algorithm runs in linear time, notice that in the while loop j p in Algorithm 4, the variable j is incremented by 1 in each iteration, and the loop terminates when j = p. Although there are two while loops inside the main while loop, the total number of iterations in the two inner while loops is at most p. This is because we start with p blocks, and each iteration of the inner while loops either merges two blocks forward or merges two blocks backward. The total number of merging operations is at most - 1. Thus, the total number of iterations in the while loop S is at most p. Lastly, since we can evaluate the proximal operator of the Huber loss function, Hy, in constant time complexity, the total time complexity of Algorithm 4 Opp).
A.6. Proof of Theorem 3.2
Theorem 3.2. For any μ E IR,, Algorithm I returns the exact evaluation of prox (H) in O(p log
Proof. By Lemmas 3.3 and 3.4, the output of Algorithm computes prox exactly. The log linear time complexity statement also holds thanks to Lemma 3.5 and the initial sorting step.
15

---

Scalable First-order Method for Certifying Optimal k-Sparse GLMs
Algorithm 4 Up and Down Block Algorithm for Merging in PAVA
Input: vector E IR,, nonnegative weights E 184 (Pli:k = P, P>Pk+11p 1:p = 0), vector D (0) prox Ma (ngy)) (, integer k E (first k elements subject to Huber penalty), and threshold M > 0 for the Huber loss function.
1: DInitialization for the first block
2: Initialize = 1, P1 = Pr. S = = MI
3: Pproy = ] = 2
4: while j 5 n do
5: b b-b+1
6:
P. Pj. 11)
7:
>ff 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 U/ree then
9:
b b-b-1
P. Po+ S. Sb + 11yl. N₂ + v prox Po HM
Sb))
11:
>Look forward: keep merging the current block with the next block if the isotonic violation persists while j < n and v S bj do
j=j+1 Po = P. Pj. S. = S 1myl No No +1, prox Po end while
14:
Had Ny.
No
16:
>Look backward: keep merging the current block with the previous block if the isotonic violation persists while b > 1 and Pb-1 < v do
17:
18:
b = 1-1 P. = Po Po-1, So = So + St41 No = No + No41. == prox Has
19:
SL
No
20:
end while
21:
end if
22:
DSave the current block's ss value and the index of the last element in the block
us = U, Tb j
DStart fresh the next element
25: Upree = v. =
26: end while
27: DModdffy the output vector to have the same new value for all elements in each block
28: for 1r l, b do
29:
P|n=itim]] 1t
30: end for
return
A.7. Proof of Theorem 3.6
Theorem 3.6. For any B IR,, Algorithm 2 computes the exact value of 9(3), defined in (7), in O(p plog k).
Proof. We first show that the algorithm correctly computes the value of g(3) and then analyze its computational complexity. Define the mixed-binary set
So = (t. Σ}E(n 83 St, || S M, || S k
Using the perspective and big-M MM reformulation techniques, the set So admits the equivalent representation
So { (1, 3z E (0, 1)2 s.t. EjE 83/2, S1, 1T 1T - Mz, S B3 S M2, Vj E [m]
Following the proof of Lemma 1. one can show that the closed convex hull of So is given by
conv(S.) (t. 3z [0, 11" s.t. 1 1 83/2, St,
M2, s B, s M2j Vj lm]
16

---

Scalable First-order Method for Certifying Optimal k-Sparse GLMs
Therefore, the implicit function y can be written the evaluation of the support function of el conv(()) (1,0), that is,
( min(t (t, cl conv(So)).
(18)
Notice that the set So is sign- and permutation cermutation-ivaariants Hence, by (Kim al., 2022, Theorem 4), its closed convex hull admits the following (different) lifted represenation
ز 03 1, ||] 3m Φ, 0 のか S ... S Φ1 S M, @ Qk-2 _... On
cl conv(So)
(t, 2 E IRP s.t.
(19)
where the absolute value operator |-| applied to vector in an element-wise fashion, and the constraint |33| 3m &b denotes that Φ majorizes [3], that is.
|
m -m Φ
Eyem1oy1 Eyem o 1I E (p 1] and Σje(p o;= -
Using this alternative convex hull description of So (19) and the implicit formulation (18), we may conclude that
. =min &cR" {2 E
03 |||| 3 Φ, s de S ... S oo S M, Φk+1 = 2x+2 :.. - =
(20)
In the following we show that Algorithm 2 can efficiently solve the minimization problem in (20). At the first iteration j = of the algorithm, we have
kΦp 2 EjE(K 0j -Ej 0j jj Φ1 2 Σ 181
At the same time, we also need to satisfy [8i] S Φ1 from the first majorization constraint. We now discuss two cases
Case 1: If 8 1831 2 1821. in order to solve the minimization problem in (20), we set 01 k 11
| 131 = 1 E Mz, S M. This leads to Notice that Φ1 S M is automatically satisfied because di E
2 = O& = 1 Σ jE
To see this, for the sake of contradition, assume that Ej E [2,...,k) such that
2 < 1 EEbn Since 6; S Φ1 = - |3, we have ΣjEN) O5 < >]EEN| 2/2 E = E)E() 18yl. which contradicts the majorization constraint.
Case 2: If
$ Eje(n 1881 < [811. we can set é1 = [311. Notice that di I M is automatically satisfied because |a S M21 M. Then are left with k 1 coefficients to and can follow the same argument as we did for 1 with slight difference that the majorization constraints are changed to
E (2, pp. and
We repeat this process until we set all & coefficients $1,..., Øb. as implemented by Algorithm 2. The output of the algorithm coincides with the optimal value of the minimization problem in (20). Hence, the first claim follows.
As for the complexity claim, it is easy to see that Algorithm only requires partial sorting step on Line 2, which has a complexity of O(p log k). The summation step on Line has a complexity of O(p). The for-loop step on Line 4-8 has complexity of O(k), does the final summation step Line 9, Therefore, the overall computational complexity of Algorithm 2 is O(p + p log k). This concludes the proof.
17

---

Scalable First-order Method for Certifying Optimal k-Sparse GLMs
Appendix B. Experimental Setup Details
B.1. Setup for Evaluating Proximal Operators
The synthetic data generation process is as follows. We sample the input vector 7 E R' from the standard multivariate Gaussian distribution, Y ~ (O, Ip), where Ip denotes the identity matrix with dimension p. We vary the dimension .......))) _ 10a and set the cardinality & be 10, the box constraint M to be 1.0, and the weight parameter p to be 1.0. We report the running time for evaluating these proximal operators. To obtain the mean and standard deviation of the running time, we repeat each setting 5 times, each with different random seed.
P E (29; 12021 21 ....,
...,,
B.2. Setup for Solving the Perspective Relaxation
We generate our synthetic datasets in the following procedure. First, we sample each feature vector 2c; RP from Gaussian distribution, 2; ~ (0, E). where the covariance matrix has entries En o The variable E (0, 1) controls the features correlation: if we increase o, feature columns in the design matrix X become more correlated. Throughout the experimental section, we set a = 0.5. Next, we create the sparse coefficient vector 3 with l equally spaced nonzero entries, where 33 = 1 if j mod (p/k) = and 3j otherwise. After these two steps, we build the prediction vector y. If our loss function is squared error loss (regression task), we set y, = x 3 + Ei where is Gaussian random noise with G ~ NO, || [xx=_|| ). and SNR stands for the signal-to-noise ratio. In all our experiments, we choose SNR 5. SNR
If our loss function is logistic loss (classification task), we set y/i Bern(a: 3* + ei), where Bern(P) is Bernoulli random variable with P(y
1) = P and Ply, = 11) 1 P. For this experiment, we vary the feature dimension p E (1000, 2000, 4000, 8000, 160003. 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 cardinality constraint) to be 10, and f2 regularization coefficient A2 to be 1.0. Again, we report and compare the running times, with means and standard deviations calculated based on 5 repeated simulations with different random seeds.
B.3. Setup for Certifying Optimality
Datasets and Preprocessing We run on both synthetic and real-world datasets. For the synthetic datasets, we run on the largest synthetic instances (n = 16000 and p = 16000). For the real-world datasets, we use the dataset cancer drug response (Liu et al., 2020) for linear regression and DOROTHEA (Asuncion & Newman, 2007) for logistic regression.
The cancer drug response dataset has 822 samples and orginally has 34674 features. However, many feature only has single value, so we prune all these features, which result in 2200 features. The DOROTHEA dataset originally has 1150 samples and 100000 features. However, the dataset is highly imbalanced (only 112 samples for the minority class). To alleviate this data imbalance problem and train a meaningful classifier, we apply the resampling preprocessing technique on the DOROTHEA dataset. We sample (with replacement) 1150 samples from each class. Thus, in total, we have 2300 samples. After pruning redundant features, we have 89989 features.
For both the cancer drug response and DOROTHEA dataset, we center each feature to have mean and norm equal to 1.
Choice of Hyperparameters For the cardinality constraint k we set ls 10 for both synthetic datasets. This corresponds to the true number of nonzero coefficients. For the cancer drug response dataset, we set h = 5. For DOROTHEA, we set k = 10. The k values for both real-world datasets are selected based on doing 5 fold cross validation with a heuristic sparse learning algorithm first. We apply the heuristics to get a path of solutions with different number of nonzero coefficients. The values are selected on this path of solutions right before the loss objective just starts to plateau or increase on the test set.
For the Eg regularization coefficient, we set A2 1. For the box constraint, we set M for the synthetic datasets and DOROTHEA. The infinity norm of the final optimal solution less than this value. For the cancer drug response dataset, we set M- 5, which is also bigger than the infinity norm of the final optimal solution.
Branch and Bound We first provide some background on the branch and bound (BnB) framework and then elaborate on the details of our BnB implementation. Problem 1 is both nonconvex and NP-hard. Thus, any method capable of certifying optimality, including heuristics, branching, bounding (calculating lower bounds), presolving, among many others. Our paper focuses on rapidly solving relaxation problems at both the root and node levels. This, in turn, provides fast lower bound certificates for the BnB algorithm.
18

---

Scalable First-order Method for Certifying Optimal k-Sparse GLMs
For our method, we write customized BnB framework. We Algorithm solve the relaxation at each node and Equation (13) calculate the safe lower bound prune the search space. To find feasible solutions, approach called beamsearch (Liu et al., 2022) from the existing literature. For branching, we branch on the feature based on the best feasible solution found by the beamsearch algorithm at each node. For the nonzero coefficients of this solution, we branch on the variable which would lead to the largest loss increase if the coefficient to The intuition that such a variable important and should be branched early the BnB framework.
effective
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 with CPU speed 2.45 Ghz and 8 cores. We conducted the experiments the the Delta system the National Center for Supercomputing Applications from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program (Boerner al., 2023).
19

---

Scalable First-order Method for Certifying Optimal k-Sparse GLMs
Appendix C. Additional Numerical Results
C.1. Perturbation Study regarding Solving the Perspective Relaxation
C.I.I. PERTURBATION STUDY ON M VALUES
Linear Loss
Logistic Loss
10a
10a
Solver Time Elapsed
Solver Time Elapsed
10a
10
10a
10
100
10
10-1
10-1
2500
5000
7500 10000
12500 15000
2500
5000
7500 10000 12500 15000
+ gurobi
mosek
clarabel
ours
Figure 4. Solve the perspective relaxation Problem (5). We
1.2, A2 1.0, n/p = 1, and
Linear Loss
Logistic Loss
103
10a
Solver Time Elapsed
Solver Time Elapsed
1
10
101
10
10a
10
10
10-1
2500
5000
7500
10000
12500 15000
2500
5000 7500 10000 12500 15000
P
gurobi
mosek
+ scs
+
clarabel
ours
Figure 5. Solve the perspective relaxation in Problem (5). We set M 155, A2 1.0. n/p = 1, and &: 10.
20

---

Scalable First-order Method for Certifying Optimal k-Sparse GLMs
Linear Loss
Logistic Loss
103
10a
Solver Time Elapsed
Solver Time Elapsed
10a
10
10
10
10
10
10
10
2500
5000
7500
10000 12500 15000
2500
5000
7500 10000 12500 15000
P
gurabi +
masek
+ scs
+
clarabel
ours
Figure 6. Solve the perspective relaxation in Problem (5). We set M = 3.0, A2 = 1.0. n/p = 1, and l: = 10.
Linear Loss
Logistic Loss
103
Solver Time Elapsed
Solver Time Elapsed
103
10
102
10
101
10
100
10
2500
5000
7500
10000
12500 15000
2500
5000
7500 10000 12500 15000
gurobi
mosek
+ clarabel
Figure 7. Solve the perspective relaxation in Problem (5). We M 5.0, 1.0, n/p 1, and & 10.
Linear Loss
Logistic Loss
1
10a
Solver Time Elapsed
Solver Time Elapsed
102
102
10a
10a
10a
10a
10-1
2500
5000
7500
10000
12500 15000
2500
5000
7500 10000 12500 15000
P
P
gurobi
mosek
scs
clarabel
ours
Figure 8. Solve the perspective relaxation in Problem (5). We set M = 10.0. A2 = 1.0, n/p 1, and &: 10.
21

---

Scalable First-order Method for Certifying Optimal k-Sparse GLMs
C.1.2. PERTURBATION STUDY ON AS VALUES
Linear Loss
Logistic Loss
103
10a
Solver Time Elapsed
Solver Time Elapsed
102
10
10a
10a
10
100
10-1
10-1
2500
5000
7500
10000
12500 15000
2500
5000
7500 10000 12500 15000
P
gurobi
mosek
scs
+
clarabel
ours
Figure 9. Solve the perspective relaxation in Problem (5). We set M 2.0, A2 = 0.1, n/p- 1, and &: 10.
Linear Loss
Logistic Loss
10
10a
Solver Time Elapsed
Solver Time Elapsed
102
10
10
10
10
10
10-1
10
2500
5000
7500
10000
12500 15000
2500
5000 7500 10000 12500 15000
P
gurobi
masek
SCs
+
clarabel
ours
Figure 10. Solve the perspective relaxation in Problem (5). We set M = 2.0, Ag = 10.0, n/p 1, and k: = 10.
22

---

Scalable First-order Method for Certifying Optimal k-Sparse GLMs
C.1.3. PERTURBATION STUDY n-TO-p RATIOS
Linear Loss
Logistic Loss
10a
Solver Time Elapsed
Solver Time Elapsed
10a
10a
10a
10a
10a
100
100
2500
5000
7500
10000
12500
15000
2500
5000
7500
10000 12500 15000
P
P
gurobi
mosek
clarabel
aurs
Figure 11. Solve the perspective relaxation in Problem (5). We set M = 2.0, A2 = 1.0, n/p = 10.0, and A: = 10.
Linear Loss
Logistic Loss
103
103
Solver Time Elapsed
Solver Time Elapsed
102
10
10
10a
10
10
10-1
2500
5000
7500
10000
12500 15000
2500
5000
7500 10000 12500 15000
P
+ gurobi
mosek
+ SCs
+
clarabel
ours
Figure 12. Solve the perspective relaxation in Problem (5). We set M 2.0, A2 1.0, n/p 0.1, and k 10.
23

---

Scalable First-order Method for Certifying Optimal k-Sparse GLMs
C.1.4. PERTURBATION STUDY LARGE k VALUE
Linear Loss
Logistic Loss
104
Solver Time Elapsed
Solver Time Elapsed
10a
10a
10
10a
10a
10a
100
100
2500
5000
7500
10000
12500
15600
2500
5000
7500 10000 12500 15000
gurobi
mosek
clarabel
ours
Figure 13. Solve the perspective relaxation in Problem (5). We set M 2.0, A2 1.0, n/p 0.1, and k 500.
24

---

Scalable First-order Method for Certifying Optimal k-Sparse GLMs
C.2. Comparison between Perspective Relaxation and E1-Relaxation
Besides the perspective relaxation discussed in Section 2, there is another common relaxation for the sparse learning problem. which is the & -erelaxation.
Recall that the perspective relaxation problem is formulated as:
min B
f(XB, y) + A2 jelp 83/2g
P+ conv
s.t.
2 2 [0, 1],, 1 S k,
(21)
- -Mzj S B3 S M Vj E [p].
In contrast, the fı-relaxation problem formulated as (following Mhenni et al. (2020, Equation 2)):
min BR=ERP
f(X3,y) y) + 2211
Poon
(22)
|3|| S kM, || M.
Since Lemma 2.1 establishes that the perspective relaxation is equivalent the convex hull of the Eg-regularization term subject to the cardinality and box constraints, this is the tightest convex relaxation for the term. Thus, the perspective relaxation is a better choice than the é é--relxxation We verify the superiority of the perspective relaxation over the é--relaxation by numerically comparing the two relaxation methods in terms of the objective values.
Table Linear Regression Results: objective values of the perspective relaxation and & relaxxttio (the higher the better).
method
p=1000
p=2000
p=4000
p=8000
p=16000
11
1088.78
2418.70
5518.70
11877.03
25710.79
perspective
1119.53
2449.34
5549.40
11907.39
25741.67
Table 4. Logistic Regression Results: objective values of the perspective relaxation and fi-relaxation (the higher the better).
method
p=1000
p=2000
p=4000
p=8000
p=16000
11
202.12
485.82
1000.84
2216.15 15
4667.79
perspective
234.34
518.99
1032.95
2248.64
4699.78
Lastly, we want to mention that there exist theoretically tigher relaxations in the literature (particularly rank-one SOCP and SDP formulations), but these formulations face fundamental computational limitations. They require interior-point methods, which cannot be warn-started effectively and also scale poorly with problem size. Our perspective relaxation provides the optimal practical balance: delivers sufficiently tight bounds while remaining computationally efficient through first-order methods.
25

---

Scalable First-order Method for Certifying Optimal k-Sparse GLMs
Appendix D. Additional Discussions
We first provide common calculus rules for conjugate functions, whose proof can be found in standard optimization textbooks such (Beck, 2017).
Separable Sum Rule: Let f(x) EjEppi (1zz), where fj : R IR is convex for all j E [p]. Then, the conjugate of is given by (m) = EjEbl (1z))
Scalar Multiplication Rule: Let g RP - IR be convex and a > 0 be scalar. Then, the conjugate of f(æ)) = ag(at) is given by f f'())= = oag" (u/a).
Addition to Affine Function Rule: Let g IRP + R be convex and a, E IR" two vectors. Then, the conjugate of ((x) g(x) a x + is given by f (14) = 9" (n a)
Composition with Invertible Linear Mapping Rule: Let g. R2 R. be convex and R2xx be an invertible matrix.
Then, the convex conjugate of (x) ((Ax:) is given by (1) = 9" (A- μ.).
Infimal Convolution Rule: Let g, h IP - IR be convex. Then, the convex conjugate of f(æ) = infoy g(y) + h(x y) is given by (1) 9" (m) + h (14).
These rules are useful for discussions in D.1 and D.2.
D.I. Convex Conjugate for GLM Loss Functions
The convex conjugates of some of GLM loss functions are summarized bellow.
Linear Regression:
F(XØ) X xx--11 y11 & F*(-() = - y²C.
Logistic Regression:
1 1- S
log (1-5)
Si
Si
73 yi
F(XB)= =
Σ
log(1 + (X & (+ = Σ
+
log
iE(m]
Vi
y/
Poisson Regression:
F(X3) Σ (exp(X33); (x3)), & F-(-C)= Σ h1- + vo),
iE(m)
iE(n)
where h(2) log(z) - if > 0 and h(z) = Of == O.
Gamma Regression:
F(X0)= = Σ
(y: cxp(-(X33)) + (x3))) & F*(-)) = Σ
yihl
Vi
rE(n]
iE|m]
where h()) = 2 <lg(22) - 2 if > 0 and M(2) = 0f O.
Squared Hinge Loss: For binary classification with labels y E (-1,+1}
F(XB)= = Σ
max(O. 1 & F-(-)) -E M M--9(S)) (E(M)
iE(m
where h(2) = + =/= if = S and h(z) = 00 if > 0.
26

---

Scalable First-order Method for Certifying Optimal k-Sparse GLMs
Multinomial Logistic Regression: For multiclass classification with K classes with coefficients 3 R let Dit be binary indicator such that Yık if the i-th sample belongs to class k and Yik 0 otherwise.
F(X Σ
log Σ exp(((X():))
Uek (x
& (-())
Σ
h(vi Si),
vE(M)
sE(m
where h(z) - 2k log(2)) i
and 1" = 1, and h(z) = oo otherwise.
D.2. Safe Lower Bound
The linear regression problem with eigen peeeppeespeettve relaxation is formulated
Ps cig-conv
min 3 Qeigg 2y x, 2B)
BERU
where Qeig
X Dmin(( (, Adi X2 + min (XT X), and (-) denotes minimum eigenvalue of the input matrix. Using the standard version of weak duality theorem, have
MIP&
P Pagg-conv 2 -
where F(B) 3 3e G(B) 2y X3 the following closed form expressions
T3 2geg (B), and Ç VF((3) -2Qeig.. The conjugate functions admit
( Qeiy (-Qcig++x XT )
F-(- = 1 - 3Qa & G²(()= = 2 9"
Nete
where use (-)t denote the gpeudo-inverss of a matrix. We may conclude that
PMIP 3Qag 2deig 9" ( [-Qigg (-Qci0(
Aeig
The above lower bound can be viewed as generalization of the safe lower bound formula from (Liu et al., 2024, Theorem 3.1). Specifically, as M approaches oo, the above lower bound matches the lower bound in in (Liu et al., 2024, Theorem 3.1). Furthermore, Our proof uses a simple weak duality argument and concise, contrast to the lengthy two-page algebraic proof of (Liu al., 2024, Theorem 3.1).
27


Reviews:
[{"review_id": "FCTUT1t5T7", "reviewer": "Reviewer_enX8", "summary": "This paper considers the problem of sparse generalized linear model (GLM) optimization. Specifically, the goal is to fit a generalized linear model under the constraint that at most $k$ features are used. Typical approaches for such problem include LASSO, greedy/local search-based techniques, as well as branch-and-bound approaches.\n\nThe authors focus on developing a fast algorithm for solving the convex relaxation that appears in each step of the branch-and-bound process. While the usual relaxation employed (e.g. in LASSO) replaces the $\\ell_0$ sparsity constraint by an $\\ell_1$ constraint, this relaxation is not necessarily tight. Instead, the authors employ convex duality to study the convex hull of the $\\ell_0$-based discrete optimization set, and end up with a tighter but non-closed form constraint. This constraint essentially regularizes the top (largest-magnitude) features separately from the bottom ones, using an $\\ell_2$ norm for the former but an $\\ell_1$ norm for the latter. The authors show how to evaluate as well as optimize over this regularizer efficiently, leading to an overall efficient algorithm for solving the relaxation. Experiments show that this specialized algorithm outperforms off-the-shelf second-order cone programming solvers.\n\nOverall, I believe this is a good paper that provides a thorough and valid analysis of the sparse GLM problem. I believe it could be used to further move the state of the art by making branch-and-bound-based sparsity solvers more efficient. This would be significant for quality-sensitive applications.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their thoughtful feedback. We address their concerns below:\n\n1. **Connection to Mhenni et al. and experimental comparison**\n\nThank you; we will cite this paper. However, we note that our perspective relaxation is tighter than the $\\ell_1$ relaxation from by Mhenni et al.\nSpecifically, Lemma 2.1 establishes that the closed convex hull of\n\n\\begin{align}\n\\\\{ ( \\beta,  z, t): ||\\beta||_2^2 \\leq t, z \\in \\\\{0, 1\\\\}^p, 1^T z \\leq k,  ||\\beta_j||  _{\\infty}\\leq M  \\\\}\n\\end{align}\nis\n\n\\begin{align}\n\\\\{ ( \\beta, z, t) : \\sum _{j=1}^p \\beta_j^2/z_j \\leq t, z \\in [0,1]^p, 1^\\top z \\leq k, -M z_j \\leq \\beta_j \\leq M z_j ~~ \\forall j \\in [p] \\\\}.\n\\end{align}\nThis is the tightest possible relaxation for the $\\ell_2$ term.\n\nThe perspective relaxation we want to solve within the BnB algorithm is:\n\\begin{align}\n\\min_{\\beta, z} || y - X \\beta || _2^2 + \\lambda_2 \\sum _{j=1}^p \\beta_j^2/z_j \\quad \\text{s.t.} \\quad z \\in [0, 1]^p, 1^T z \\leq k, -M z_j \\leq \\beta_j \\leq M z_j\n\\end{align}\n\nIn contrast, the $\\ell_1$-based relaxation proposed by Mhenni et al. (Equation 2 in their paper and subsequent discussion) will solve:\n\\begin{align}\n\\min_{\\beta, z} || y - X \\beta ||_2^2 + \\lambda_2 ||\\beta||_2^2 \\quad \\text{s.t.} \\quad ||\\beta||_1 \\leq k M, ||\\beta|| _{\\infty} \\leq M\n\\end{align}\n\nAs Lemma 2.1 indicates, our perspective relaxation provides tighter lower bounds, crucial for pruning nodes and certifying optimallity.\nBelow we provide some numerical comparison (same setup as Figure 2 in the paper):\n\nLinear regression (the higher the better):\n| method      | p=1000  | p=2000  | p=4000  | p=8000   | p=16000  |\n| :-----:     | :-----: | :-----: | :-----: | :-----:  | :-----:  |\n| l1          | 1088.78 | 2418.70 | 5518.70 | 11877.03 | 25710.79 |\n| perspective | 1119.53 | 2449.34 | 5549.40 | 11907.39 | 25741.67 |\n\nLogistic regression (the higher the better):\n| method      |  p=1000  | p=2000  | p=4000  | p=8000  | p=16000 |\n| :-----:     |  :-----: | :-----: | :-----: | :-----: | :-----: |\n| l1          |  202.12  | 485.82  | 1000.84 | 2216.15 | 4667.79 |\n| perspective |  234.34  | 518.99  | 1032.95 | 2248.64 | 4699.78 |\n    \n2. **Does Algorithm 2 only work when $M=\\infty$?**\n\nAlgorithm 2 is valid for any $M > 0$.\nIn Algorithm 3, we first compute the proximal update (line 8), before evaluating $g(\\beta ^{t+1})$ in Algorithm 2.\nSince $\\beta ^{t+1}$ is the output of the proximal operator, it inherently satisfies all constraints in Equation (7).\nNamely, there exists some $z \\in [0, 1]^p$ such that $1^T z \\leq k$ and $-M z_j \\leq \\beta_j \\leq M z_j \\; \\forall j \\in [p]$.\nThis property is critical to the proof of Theorem 3.6 (Line 896-897 and Line 902-903), where we show that Algorithm 2 correctly computes $g(\\beta ^{t+1})$.\nWe will clarify and emphasize this point in revision.\n\n3. **Comparison with ISTA**\n\nIf by ISTA the reviewer refers to the non-accelerated version of FISTA, the comparison is already included in Figure 3, as PGD (Proximal Gradient Descent) and ISTA are essentially the same algorithm. We used the more general term PGD since ISTA typically implies an $\\ell_1$ regularizer. We would be happy to rename legend to ''PGD/ISTA'' in revision. Please let us know if we have misinterpreted your meaning regarding ISTA.\n\n4. **Connection to optimal transport/deep learning literature for soft-top-k operator**\n\nThank you; we will cite these two papers.\n\nHowever, the smooth approximation of the top-k operator is not suitable here.\nThe effectiveness of the proximal algorithm relies on the exact evaluation of the proximal operator. Replacing the top-k would lead to solving a different problem rather than the proximal operator evaluation.\nThis will not help us to use FISTA to solve the perspective relaxation. As a result, this approach would not guarantee valid lower bounds necessary for optimality certification.\n\n\nWe hope to have answered to your concerns and questions, and we would be happy to provide further explanations.", "claims_and_evidence": "I found all the claims presented to be valid.", "methods_and_evaluation_criteria": "The method is evaluated on synthetic tasks, using both CPU and GPU, as well as two real linear/logistic regression datasets, which are highly over-parametrized. The experiments are performed multiple times to reduce variance. The results show a significant benefit of using the proposed method (the version with ranodm restarts) compared to off-the-shelf optimizers like gurobi.", "theoretical_claims": "All the theoretical claims look correct to me. I checked the proof of Theorem 3.6 and all the claims in the main body.", "experimental_designs_or_analyses": "N/A", "relation_to_broader_scientific_literature": "The idea of solving sparse optimization with branch-and-bound has been studied before, e.g. \"Sparse Branch and Bound for Exact Optimization of L0-Norm Penalized Least Squares\" by Mhenni et al. However that work uses an $\\ell_1$ relaxation. I am not aware of any work using a tighter relaxation.", "essential_references_not_discussed": "N/A", "other_strengths_and_weaknesses": "The paper is well-written.\n\nThe box constraint parameter $M$ is introduced but then ignored in Algorithm 2 (I guess it assumes $M=\\infty$). It might also be implicitly assumed in other algorithms / theorems. It would be good for the authors to fix this.", "other_comments_or_suggestions": "- Add some comparison with ISTA.\n- An optional suggestion is to take a look at papers in the optimal transport / deep learning literature for similar approaches regarding a soft-top-k operator. https://arxiv.org/pdf/2002.06504 https://arxiv.org/pdf/2304.04947 At a high level they have the same goal of inducing tighter relaxations to the $\\ell_0$ operator.", "questions_for_authors": "- How do the results compare to ISTA (i.e. no acceleration)?\n- Is there any experimental evidence that the relaxation used here is superior to the $\\ell_1$-based relaxation? It would be great to see at least one experiment or some pointers to other works that do that.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "wluTiDCzsS", "reviewer": "Reviewer_3khC", "summary": "The paper proposes a new algorithm for solving lower bounds in the BnB framework. It solves a composite problem using an efficient restarted FISTA algorithm, for which an efficient way to exactly compute the function value and the proximal operator is given. As such, the paper allows to achieve impressive efficiency gains compared to existing methods.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their thoughtful feedback. We address their concerns below:\n\n1. **Connection to Argyriou et al. and McDonald et al.**\n\nThank you for these valuable references--we will cite and discuss both papers.\nWe fully agree with your observations. \nWe clarify some key differences:\n\n**Difference 1** Our function $g(\\beta)$ generalizes the k-support norm. Specifically, when  $M=\\infty$, our formulation reduces exactly to the standard k-support norm. However, our work makes several novel contributions beyond prior works:\n\n   a. We are the first to study the Fenchel dual of $g$\n\n   b. We derive an explicit closed-form expression for $g^*$ (Equation 9)\n\n   c. Our analysis captures bounded variable constraints ($ M < \\infty $) \n\n**Difference 2** From a computational perspective, ours differs fundamentally from prior works:\n\n   a. Algorithmic Framework: Both cited papers compute the proximal operator directly in the primal space (their Algorithm 1's), while our paper employs a dual approach by computing the proximal operator of the dual function $g^*$ and applying the Moreau decomposition theorem (Equation 12)\n\n   b. Constraint Handling: To the best of our knowledge, the primal approach struggles to incorporate box constraints effectively, whereas our dual framework naturally accommodates box constraints while maintaining computational efficiency. \n\n   c. Implementation Advantages: The dual perspective provides a more flexible computational framework, enabling efficient handling even for constrained optimization scenarios\n\n2. **Does $\\tilde{O}(p)$ mean $O(p \\log(p))$?**\n\nYes, this is correct.\nIn Algorithm 1, the most expensive step is the sorting step (line 3), which has computational cost $O(p \\log(p))$.\nAll other steps have computational costs $O(p)$.\nThus, the overall computational cost is $O(p \\log(p))$.\n\n3. **Is our method (BnB with lower bounds computed by our Algorithm 3) better than all other methods in solving the overall Problem 1? Under what conditions?**\n\nTo properly address this question, let us first provide some context about BnB algorithms. Problem (1) is both nonconvex and NP-hard.\nThus, any method capable of certifying optimality, including ours, must ultimately rely on BnB. This BnB framework consists of several key components, including heuristics, branching, bounding (calculating lower bounds), presolving, among many others. \nOur paper focuses on rapidly solving relaxation problems at both the root and node levels. This, in turn, provides fast lower bound certificates for the BnB algorithm.\n\nWith this background information, we can now more directly address your 1st question. The answer is yes in the sense that it is currently the method that achieves the best results and scalability thanks to efficient lower bound computation. The key advantage lies in our formulation's compatibility with first-order optimization methods, which enables warm-starting and maintains computational tractability even for large-scale problems.\nWhile there exist theoretically tighter relaxations in the literature (particularly rank-one SOCP and SDP formulations), these formulations face fundamental computational limitations. They require interior-point methods, which cannot be warm-started effectively and also scale poorly with problem size. Our perspective relaxation provides the optimal practical balance: it delivers sufficiently tight bounds while remaining computationally efficient through first-order methods.\n\nFor your 2nd question, our method performs the best when applied to **large-scale** problems, on which existing lower bound computation methods scale poorly. In such cases, our approach provides substantial improvements in both speed and scalability. For small size problems, we observe that computing lower bounds is generally computationally cheap. In these scenarios, both our branch-and-bound implementation and commercial mixed-integer programming solvers can typically certify optimality within seconds. While these smaller cases serve as useful validation of our method's correctness, they are not the primary focus of our work, as they do not present the same computational challenges that motivate our technical contributions.\n\nWe sincerely appreciate your thoughtful questions and engagement with our work.\nWe hope to have answered to your concerns and questions and we would be happy to provide further explanations.", "claims_and_evidence": "The authors prove their main claims or give appropriate references when necessary (Section 2 and 3), and provide extensive experiments to show the advantage of their method (Section 4)", "methods_and_evaluation_criteria": "The setting of the paper makes a lot of sense: it builds on existing methods (BnB methods using lower bounds based on relaxations, to solve sparse optimization problems), and replaces the part of the lower bound computation by an efficient FISTA algorithm using a novel proximal operator.", "theoretical_claims": "I did not check the correctness of the proofs in details, but I am familiar with similar literature on the computation of proximal operators for sparse penalties, and the results make sense to me at first sight.", "experimental_designs_or_analyses": "The experiments are synthetic experiments and their setting makes sense to me.", "relation_to_broader_scientific_literature": "This paper give a more efficient way to solve sparse optimization problems, which is very useful in many areas of machine learning such as medical scoring systems or portfolio optimization, as described in the paper's introduction.", "essential_references_not_discussed": "This is just a suggestion, but I think it might be interesting to cite the papers [1] or [2] below related to the k-support norm. Indeed, the function $g^{\\*}$ in the paper is very similar to the (half-squared) top k-norm (the top-k norm is the norm of the top-k elements in absolute value), which is the dual function of the (half-squared) k-support norm (since the top-k norm is the dual norm of the k-support norm). More precisely, for a given $\\boldsymbol{\\alpha}$, if all for all $i$, $|\\alpha_i| \\leq M$, then $g^*(\\boldsymbol{\\alpha})$ is exactly the half-squared top-k norm.\nNow, the proximal operator of the k-support norm is known [1] and can be efficiently computed from [2]: therefore I wonder if a similar technique could be used to give an alternative computation for the proximal operator of $g^*$ (related to the one of $g$). Indeed, Algorithm 1 seems to involve different techniques, with the merging of intervals. In any case, I notice that in [2], the complexity of the proximal operator is in $O(p log(p))$, which is the same as the proximal operator in this reviewed paper (assuming $\\tilde{O}(p)$ in the paper means $O(p \\log(p))$, not $O(p   \\text{poly}(\\log(p))$ for some higher order polynomial ? Maybe the authors can confirm). Therefore, I think this reviewed paper gives a very efficient formulation, and therefore this remark is more just for sake of curiosity rather than a suggestion for improving the algorithm.\n\n[1] Sparse Prediction with the k-Support Norm, Argyriou et al.\n\n[2] New Perspectives on k-Support and Cluster Norms, McDonald et al.", "other_strengths_and_weaknesses": "I think this paper uses ingenious techniques to make the computation of lower bounds in BnB much more effective. The paper uses innovative techniques for the computation of the proximal operator, and I think the restart technique is also welcomed, to improve the efficiency of FISTA: that last technique actually makes a big difference as the authors show in their experiments.\nThe overall complexity of the algorithm is shown to be much better than state of the art solvers like Gurobi, which is very encouraging.", "other_comments_or_suggestions": "See questions below.", "questions_for_authors": "I first have just a small question regarding the comparison with the k-support norm, although the authors can omit it as it is just for sake of curiosity rather than a suggestion of improvement.\n\nAnd I would also have another second question regarding the performance of the full algorithm plugged in the BnB framework, to solve (1). Indeed, I think it would be interesting to compare the complexity or running time of solving the full optimization problem (1) with the method in the paper (BnB with lower bounds based on the improved FISTA) vs. with other methods. Can we already say that the method in the paper would also be more efficient than other methods for solving (1) ? More precisely: this paper provides a more efficient first order method to compute lower bounds for BnB. In the introduction, it is said such first order methods are usually slow, but give tighter bounds than methods based e.g. on conic relaxations. Now, since the method in the paper is faster, can we say that overall (i.e. for solving the full (1)), it is a better methods than all other alternatives ? Under which conditions is it ? (Note: I am less familiar with the literature on BnB so feel free to react if my question needs to be reformulated/updated)", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "0Lqtbw47Vb", "reviewer": "Reviewer_UQGX", "summary": "This paper studies sparse generalized linear models with cardinality constraints. Existing branch-and-bound methods are not computationally efficient due to expensive or slow dual bound computations. To overcome this, the authors propose a first-order proximal gradient algorithm to solve the perspective relaxation of the problem. The method eliminates the need for costly second-order cone programming and incorporates a restart strategy to accelerate convergence. Experiments demonstrate improvements in computation.", "questions": "", "limitations": "", "rating": 2, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the thoughtful feedback.\n\nWe briefly restate our contribution to avoid any misunderstanding.\nOur paper addresses certifying optimality or quantifying the optimality gap for the $\\ell_0$-constrained GLMs, \\textbf{without making any assumptions on the data}.\nOur motivation stems from the need for optimal solutions in high-stakes applications, where solving the $\\ell_0$ problem ensures superior quality, interpretability, and trustworthiness compared to heuristics.\nHowever, the problem is NP-hard, necessitating a BnB algorithm.\nThis BnB framework consists of several key components, including heuristics, branching, bounding (calculating lower bounds), pre-solving, among many others.\nThe effectiveness of BnB methods strongly depends on the quality and speed of solving its continuous relaxation (lower bound).\nIn this work, we have proposed a first-order method to solve it efficiently. By doing so, we have significantly expanded the size/scale of the datasets on which the BnB approach can be successfully applied.\n\nWith this explanation and context, we can now answer the questions raised by the reviewer.\n\n\n1. **Confusion about the problem to solve**\n\nOur contribution is tailored at efficiently computing the lower bound for the $\\ell_0$-constrained GLMs as the main workhorse for solving the full problem by BnB.\nWe hope that by adjustments to the paper (and the title, see below) this confusion is removed. \n\n2. **Optimality is supported without theoretical justification**\n\nIt should be now clear that optimality is guaranteed by the use of the BnB algorithm.\n\n3. **Title is misleading**\n\nHopefully we have clarified the reviewer's misunderstanding and explained why the title represents our goal and contribution well.\nHowever, we are open to discuss an alternative title along the line of ``Certifying Optimality of k-Sparse GLMs by Scalable First-order Lower Bound Computation\".\n\n4. **Some assumptions are not clearly stated or justified**\n\nWe do not impose any assumption on our datasets. However, if the reviewer thinks there are specific parts of paper that require more explicit discussion, we are open to follow them up to improve the clarity of the paper.\n\n5. **Cite Blanchard and Tanner**\n\nThank you for this reference. We will cite this work in revision. While both works use GPU acceleration, they address fundamentally different problems.\nThe cited work focuses on accelerating greedy/heuristic methods through GPUs.\nThere is no requirement of proving optimality.\nIn contrast, our work harnesses GPU specifically for computing lower bounds in BnB, which is crucial for pruning nodes and certifying optimality.\nBy relying solely on matrix-vector multiplications rather than solving linear systems, we achieve the first truly GPU-efficient implementation for computing the BnB's lower bounds.\n\n\n6. **$k$ is not large enough**\n\nDue to space limit, please see our reply (Q1) to Reviewer X5XP.\n\n7. **Poisson regression**\n\nWe run the setup in Figure 2 (main paper) with Poisson regression for MOSEK, ours, and ours+GPU.\nThe running time (in seconds) are reported below. \"***\" denotes reaching the time limit (1800s).\n\n| method   |  p=1000   |  p=2000   |  p=4000   |  p=8000   |  p=16000   |\n| :-----:  |  :-----:  |  :-----:  |  :-----:  |  :-----:  |  :-----:   |\n| MOSEK |  2.537    |  15.505    |  136.790   |  919.412   |   ***  |\n| ours CPU  |  2.505   |  12.821    |  95.583   |  572.876  |  ***   |\n| ours GPU |  2.637 | 6.142 |  5.933   |  17.338   |  25.954   |\n\nThe GPU version of our method scales very well on large-scale instances.\n\n8. **Cite more papers when $M=\\infty$**\n\nDue to space limit, please see our reply (Q1) to Reviewer 3khc. In revision, we will cite those two papers and discuss the key differences.\n\n9. **More background information on FISTA**\n\nThank you for the suggestion. In revision, we will include a section in the appendix covering the FISTA algorithm and the Nesterov acceleration method.\n\n10. **usage of  alas**\n\nThank you very much. We will change ''alas'' to make clear that we meant the sentence in adversative form, i.e., we will use ''however'' in the revision.\n\nWe sincerely appreciate your thoughtful questions and engagement with our work. \nWe hope to have answered to your concerns and questions and we would be happy to provide further explanations.", "claims_and_evidence": "The title of the paper is somewhat misleading, as the proposed method primarily relies on a convex relaxation rather than directly addressing the original sparse ridge regression problem (1). Additionally, the claims about the optimality gap are supported only by numerical experiments, without theoretical justification.", "methods_and_evaluation_criteria": "- The theoretical result in the paper only focus on the properties in computation while lacks of results in optimization. For example, what is the quality of solution? \n- What is the empirical performance the method when applying the method to Poisson regression?", "theoretical_claims": "I have reviewed Appendix A and examined the proof of Lemma 3.1 in detail. While the proof appears to be correct, some assumptions are not clearly stated or justified in the paper.", "experimental_designs_or_analyses": "The experimental design does not consider scenarios where $k$ is large (e.g., $p/2$), especially when $p$ is also large.", "relation_to_broader_scientific_literature": "The key contributions of this paper align with the broader literature on optimization with sparsity constraints. Thus, it may be helpful for high-dimensional data analysis.", "essential_references_not_discussed": "Using GPUs to speed up sparse learning is not new. There is existing literature on this topic (see, e.g., [1]). Therefore, I believe some statements in this paper to be overly exaggerated.\n\n- [1] Blanchard, Jeffrey D., and Jared Tanner. \"GPU accelerated greedy algorithms for compressed sensing.\" Mathematical Programming Computation 5.3 (2013): 267-304.", "other_strengths_and_weaknesses": "- Strengths\nThe empirical results show that the proposal is promising for reducing computational time.\n\n- Weaknesses\n  - I believe the special case of Lemma 3.1, where $M=\\infty$, likely appears in existing literature, but this is not discussed in the paper.\n  - The paper is not well-structured:\n    - It is confusing about the problem the authors aim to solve.\n    - Some information about FISTA should be presented before the Methodology section.\n    - The assumptions in the paper are not clearly stated.", "other_comments_or_suggestions": "Typos:\n\n- Table 1: p should be presented in a math form. \n- \"Alas\" --> \"Also\"", "questions_for_authors": "Please refer to **Claims and Evidence**, **Methods and Evaluation Criteria**, and **Weaknesses**.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "tttfr3KuCn", "reviewer": "Reviewer_X5XP", "summary": "This paper explores the use of branch-and-bound (BnB) frameworks to solve sparsity-constrained optimization problems.\n\nThe authors derive a lower bound for the optimization problem using a perspective relaxation formulation.\n\nTo efficiently solve the resulting perspective relaxation, the authors employ a first-order proximal gradient algorithm.\n\nExtensive experiments on both synthetic and real-world datasets demonstrate that the proposed approach accelerates dual bound computations and improves efficiency.", "questions": "", "limitations": "", "rating": 2, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for thoughtful feedbacks. \n\n1. **k is not large enough**\n\na. *Many real-world datasets are naturally sparse*. Small $k$'s are sufficient for accurate prediction and can help avoid overfitting, especially on the validation set. Small $k$ also improves interpretability. This is true for the two real-world datasets used in the paper. We can follow the experimental setup section (Line 973-977) to select the correct $k$ via cross-validation. The best $k$ for Cancer Drug Response is $5$ (already used), and the best $k$ for DOROTHEA is $26$ (in the submission, we used $15$). Below is the BnB runtime comparison when $k=26$ for DOROTHEA. Our method is still significantly faster than MOSEK.\n\n|    | MOSEK  |  ours  |\n|:---:|:---:|:---:|\n| running time (s) | 661.33 | 228.34 |\n\nb. *Some data science practitioners require $k$ to be small*. This is the case for the sparse identification of dynamical systems  (Bertsimas & Gurnee, 2023). Physicists often want the differential equation to be succinct in order to match some physical intuition or prior knowledge.\n\n\nc. Lastly, *our lower bound calculation is still fast when $k$ is large*. We re-ran Fig. 2 (main paper) with $k=500$ (runtime in seconds):. We only compare with MOSEK as other three baselines are already struggling when $p=4000$.\n\nLinear regression:\n| method  |  p=1000   |  p=2000   |  p=4000   |  p=8000   |  p=16000   |\n| :-----: |  :-----:  |  :-----:  |  :-----:  |  :-----:  |  :-----:   |\n| MOSEK   |  0.970    |  3.097    |  13.565   |  115.554  |  573.247   |\n| ours    |  0.327    |  0.876    |  10.623   |  44.638   |  108.475   |\n\nLogistic regression:\n\n| method  |   p=1000   |  p=2000   |  p=4000   |  p=8000   |  p=16000   |\n| :-----: |   :-----:  |  :-----:  |  :-----:  |  :-----:  |   :-----:  |\n| MOSEK   |   1.794    |  9.266    |  49.801   |  284.638  |  1969.305  |\n| ours    |   0.601    |  5.993    |  42.802   |  177.481  |   640.345  |\n\n2. **Novelty is unclear: perspective formulation, BnB algorithm, APG algorithm, and PAVA are not new**\n\nWe clarify the distinction between prior work and our contributions:\n\na. *BnB*: Our key innovation lies not in using BnB but in developing an efficient 1st order method for computing the BnB lower bounds.\n\nb. *Perspective relaxation*: Again, the novelty is not using the perspective relaxation but in solving it.\n\nc. *APG*: Although the APG framework is standard, its direct application to our problem is impossible without our method to evaluate the proximal operator efficiently and exactly.\n\nd. *PAVA*: Novelty is not using PAVA but linking $\\text{prox}_{\\rho^{-1}g}(\\cdot)$ with PAVA via generalized isotonic regression reformulation and then invoking the Moreau decomposition Theorem.\n\nAdditional contributions include:\n\na. Lemma 3.1 shows $g^*$ is a simple composition of the Huber loss and the top-$k$ norm.\n\nb. Our Algorithm 2 computes $g(\\beta ^{t+1})$ exactly and efficiently, crucial for restart.\n\nc. Our lower bound computation is GPU-friendly\n    \n3. **Justification for perspective formulation is not clear**\n\nOur main goal is to certify optimality, which necessitates using a BnB framework, whose efficiency depends strongly on obtaining tight lower bounds.\nOur perspective relaxation produces tighter lower bounds than traditional relaxation approaches (such as the $\\ell_1$ relaxation in our Q1 to Reviewer enX8).\n \n4. **Greedy/heuristic algorithms**\n\nWe completely agree that heuristics are very effective.\nThe heuristic method used in our work is already a variant of OMP (Line 983-989; also see our source code) and can sometimes find the optimal solution even within a few seconds (the whole BnB algorithm can still take over hundreds of seconds though).\nHowever, this is not the focus of this work and orthogonal to our goal, which is to certify optimality.\nOur contributions lies in efficiently computing the lower bounds at each node.\n\n5. **Meaningful real-world applications for certifying optimality**\n\nPlease see Section 1.2. Related Works for more applications.\nWe highlight a few below:\n\na. *Sparse identification of dynamical systems (Bertsimas \\& Gurnee, 2023)*: The datasets contain highly-correlated features, on which solving to optimality can greatly outperform heuristic methods.\n\nb. *Medical scoring systems (Ustun \\& Rudin, 2019)*: Certifying optimality ensures that physicians provide the best disease diagonostic system to the patients.\n\nc. *Portfolio optimization (Bienstock, 1996)*: Solving the problem to optimality ensures that companies are picking the optimal financial trading strategy.\n\nWe sincerely appreciate your thoughtful questions and engagement with our work.\nWe hope to have answered to your concerns and questions and we would be happy to provide further explanations.", "claims_and_evidence": "yes", "methods_and_evaluation_criteria": "The parameter $k$, set to $k \\in $\\{5,10,15\\}$, is too small.", "theoretical_claims": "yes", "experimental_designs_or_analyses": "yes", "relation_to_broader_scientific_literature": "Most of the results are empirical, with a lack of theoretical analysis.", "essential_references_not_discussed": "yes", "other_strengths_and_weaknesses": "### **Strengths:**  \n- **S1.** The proposed method significantly accelerates computations compared to previous BnB algorithms.  \n- **S2.** The authors derive the proximal operator of the function $g(\\beta)$ in problem (7). The authors demonstrate that the customized PAVA algorithm in Algorithm 1 has an $\\tilde{O}(p)$ time complexity.  These results are reasonable to me.\n\n\n### **Weaknesses:**  \n- **W1.** The main novelty of the paper is unclear. The perspective formulation is not new, the BnB framework and the APG algorithm are well-established, and the Pool Adjacent Violators Algorithm (PAVA) has already been used for solving isotonic regression problems and SLOPE models.  \n\n- **W2.** The justification for using the perspective formulation to estimate the lower bound is not well explained.  \n\n- **W3.** The Customized PAVA algorithm in Algorithm 1 appears to be an incremental extension of the standard PAVA algorithm. Its design is primarily based on leveraging the property that the optimal solution of problem (11) maintains the same order as the input data \n$\\mu$. While the algorithm introduces some new ideas, its overall contribution seems limited.\n\n- **W4.** The parameter $k$ used in the experiments is too small, set to $k \\in \\{5,10,15\\}$. With such small values, standard local algorithms such as OMP, iterative hard thresholding, BCD methods would likely perform well. This raises concerns about whether the proposed method offers a meaningful improvement in real-world applications.", "other_comments_or_suggestions": "No", "questions_for_authors": "- **Q1.** How does the scalability of the proposed method change when $k = 0.1p$?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "2lThYJ6v6I", "reviewer": "Reviewer_xg9D", "summary": "In this work, the authors proposed a novel FISTA-based algorithm for computing a lower bound in the Branch-and-Bound method. The new algorithm utilizes several customized components and outperforms universal optimization solvers on both artificial and practical datasets. Please see the following sections for my detailed comments.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for carefully reading our paper! This is indeed a typo. The implementation is correct in our submitted source code. We will fix this typo in the revision.", "claims_and_evidence": "Due to the time limit, I did not check the correctness of the theory, except those briefly mentioned in the main paper. The theoretical claims seem correct by checking the main paper.", "methods_and_evaluation_criteria": "The methods and evaluation criteria make sense to me.", "theoretical_claims": "Due to the time limit, I did not check the correctness of the theory. The theoretical claims and the proofs in the main paper seem correct.", "experimental_designs_or_analyses": "The experiment settings and results are sound and reasonable to me.", "relation_to_broader_scientific_literature": "This paper is related to the topic of ICML conference and should be interesting to audiences from machine learning and optimization fields.", "essential_references_not_discussed": "N/A.", "other_strengths_and_weaknesses": "Please see my comments in other sections.", "other_comments_or_suggestions": "- Algorithm 3, Line 5: should the step size be $\\phi / (\\phi + 3)$?", "questions_for_authors": "Please see my comments in other sections.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: The paper presents a first-order proximal gradient algorithm for optimizing lower bounds within a BnB framework for solving the $k$-sparse GLM problem. The paper received five detailed reviews. Three reviews are in favor of accepting this work, while the other two lean towards rejection. While the opinions expressed by the reviewers are somewhat mixed, the majority opinion is that the proposed FISTA-type lower bound optimizer in BnB is innovative, and it works to considerably reduce the overall computation cost of non-convex sparse GLM.  The reviewers also raised some concerns regarding the inadequate references and missing theoretical/empirical comparisons with some related works, which have been thoroughly addressed by authors in their responses. Based on the reviews and author responses, it was assessed that this submission makes a sufficiently novel and solid contribution to the addressed problem, and thus may be accepted. The authors are encouraged to address the issues raised in the reviews as much as possible in the revised paper.
Decision: Accept (poster)