Scalable First-order Method for Certifying Optimal k-Sparse GLMs

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.

We thank the reviewer for thoughtful feedbacks.

1. k is not large enough

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

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

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

Linear regression:

Logistic regression:

2. Novelty is unclear: perspective formulation, BnB algorithm, APG algorithm, and PAVA are not new

We clarify the distinction between prior work and our contributions:

a. BnB: Our key innovation lies not in using BnB but in developing an efficient 1st order method for computing the BnB lower bounds.

b. Perspective relaxation: Again, the novelty is not using the perspective relaxation but in solving it.

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

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

Additional contributions include:

a. Lemma 3.1 shows $g^*$ is a simple composition of the Huber loss and the top-$k$ norm.

b. Our Algorithm 2 computes $g(\beta ^{t+1})$ exactly and efficiently, crucial for restart.

c. Our lower bound computation is GPU-friendly

3. Justification for perspective formulation is not clear

Our main goal is to certify optimality, which necessitates using a BnB framework, whose efficiency depends strongly on obtaining tight lower bounds. Our perspective relaxation produces tighter lower bounds than traditional relaxation approaches (such as the $\ell_1$ relaxation in our Q1 to Reviewer enX8).

4. Greedy/heuristic algorithms

We completely agree that heuristics are very effective. The 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). However, this is not the focus of this work and orthogonal to our goal, which is to certify optimality. Our contributions lies in efficiently computing the lower bounds at each node.

5. Meaningful real-world applications for certifying optimality

Please see Section 1.2. Related Works for more applications. We highlight a few below:

a. Sparse identification of dynamical systems (Bertsimas & Gurnee, 2023): The datasets contain highly-correlated features, on which solving to optimality can greatly outperform heuristic methods.

b. Medical scoring systems (Ustun & Rudin, 2019): Certifying optimality ensures that physicians provide the best disease diagonostic system to the patients.

c. Portfolio optimization (Bienstock, 1996): Solving the problem to optimality ensures that companies are picking the optimal financial trading strategy.

We sincerely appreciate your thoughtful questions and engagement with our work. We hope to have answered to your concerns and questions and we would be happy to provide further explanations.

We thank the reviewer for their thoughtful feedback. We address their concerns below:

1. Connection to Argyriou et al. and McDonald et al.

Thank you for these valuable references--we will cite and discuss both papers. We fully agree with your observations. We clarify some key differences:

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:

a. We are the first to study the Fenchel dual of $g$

b. We derive an explicit closed-form expression for $g^*$ (Equation 9)

c. Our analysis captures bounded variable constraints ($M < \infty$)

Difference 2 From a computational perspective, ours differs fundamentally from prior works:

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)

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.

c. Implementation Advantages: The dual perspective provides a more flexible computational framework, enabling efficient handling even for constrained optimization scenarios

2. Does $\tilde{O}(p)$ mean $O(p \log(p))$?

Yes, this is correct. In Algorithm 1, the most expensive step is the sorting step (line 3), which has computational cost $O(p \log(p))$. All other steps have computational costs $O(p)$. Thus, the overall computational cost is $O(p \log(p))$.

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

To properly address this question, let us first provide some context about BnB algorithms. Problem (1) is both nonconvex and NP-hard. Thus, 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. 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.

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

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

We sincerely appreciate your thoughtful questions and engagement with our work. We hope to have answered to your concerns and questions and we would be happy to provide further explanations.

We thank the reviewer for the thoughtful feedback.

We briefly restate our contribution to avoid any misunderstanding. Our paper addresses certifying optimality or quantifying the optimality gap for the $\ell_0$-constrained GLMs, without making any assumptions on the data. Our 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. However, the problem is NP-hard, necessitating a BnB algorithm. This BnB framework consists of several key components, including heuristics, branching, bounding (calculating lower bounds), pre-solving, among many others. The effectiveness of BnB methods strongly depends on the quality and speed of solving its continuous relaxation (lower bound). In 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.

With this explanation and context, we can now answer the questions raised by the reviewer.

1. Confusion about the problem to solve

Our 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. We hope that by adjustments to the paper (and the title, see below) this confusion is removed.

2. Optimality is supported without theoretical justification

It should be now clear that optimality is guaranteed by the use of the BnB algorithm.

3. Title is misleading

Hopefully we have clarified the reviewer's misunderstanding and explained why the title represents our goal and contribution well. However, 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".

4. Some assumptions are not clearly stated or justified

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

5. Cite Blanchard and Tanner

Thank you for this reference. We will cite this work in revision. While both works use GPU acceleration, they address fundamentally different problems. The cited work focuses on accelerating greedy/heuristic methods through GPUs. There is no requirement of proving optimality. In contrast, our work harnesses GPU specifically for computing lower bounds in BnB, which is crucial for pruning nodes and certifying optimality. By 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.

6. $k$ is not large enough

Due to space limit, please see our reply (Q1) to Reviewer X5XP.

7. Poisson regression

We run the setup in Figure 2 (main paper) with Poisson regression for MOSEK, ours, and ours+GPU. The running time (in seconds) are reported below. "***" denotes reaching the time limit (1800s).

The GPU version of our method scales very well on large-scale instances.

8. Cite more papers when $M=\infty$

Due to space limit, please see our reply (Q1) to Reviewer 3khc. In revision, we will cite those two papers and discuss the key differences.

9. More background information on FISTA

Thank you for the suggestion. In revision, we will include a section in the appendix covering the FISTA algorithm and the Nesterov acceleration method.

10. usage of alas

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

We sincerely appreciate your thoughtful questions and engagement with our work. We hope to have answered to your concerns and questions and we would be happy to provide further explanations.

We thank the reviewer for their thoughtful feedback. We address their concerns below:

1. Connection to Mhenni et al. and experimental comparison

Thank 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. Specifically, Lemma 2.1 establishes that the closed convex hull of

$

\{ ( \beta, z, t): ||\beta||_2^2 \leq t, z \in \{0, 1\}^p, 1^T z \leq k, ||\beta_j|| _{\infty}\leq M \}

$

is

$

\{ ( \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] \}.

$

This is the tightest possible relaxation for the $\ell_2$ term.

The perspective relaxation we want to solve within the BnB algorithm is:

$

\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

$

In contrast, the $\ell_1$-based relaxation proposed by Mhenni et al. (Equation 2 in their paper and subsequent discussion) will solve:

$

\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

$

As Lemma 2.1 indicates, our perspective relaxation provides tighter lower bounds, crucial for pruning nodes and certifying optimallity. Below we provide some numerical comparison (same setup as Figure 2 in the paper):

Linear regression (the higher the better):

Logistic regression (the higher the better):

2. Does Algorithm 2 only work when $M=\infty$?

Algorithm 2 is valid for any $M > 0$. In Algorithm 3, we first compute the proximal update (line 8), before evaluating $g(\beta ^{t+1})$ in Algorithm 2. Since $\beta ^{t+1}$ is the output of the proximal operator, it inherently satisfies all constraints in Equation (7). Namely, 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]$. This 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})$. We will clarify and emphasize this point in revision.

3. Comparison with ISTA

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

4. Connection to optimal transport/deep learning literature for soft-top-k operator

Thank you; we will cite these two papers.

However, the smooth approximation of the top-k operator is not suitable here. The 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. 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.

We hope to have answered to your concerns and questions, and we would be happy to provide further explanations.