# On Strenghtenings for the Feature Selection Problem

Renaud Chicoisne, Pierre Latouche, Rui Ma

► To cite this version:

Renaud Chicoisne, Pierre Latouche, Rui Ma. On Strenghtenings for the Feature Selection Problem.
2025. hal-04927551

HAL Id: hal-04927551

https://hal.science/hal-04927551v1

Preprint submitted on 3 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# ON STRENGHTENINGS FOR THE FEATURE SELECTION PROBLEM

Renaud CHICOISNE<fnref n="1" /> & Pierre LATOUCHE <fnref n="2" />,<fnref n="4" /> & Rui MA <fnref n="3" />

**Résumé.** Le principe d'optimisation parcimonieuse, visant à résoudre le problème de sélection des caractéristiques, suscite un intérêt croissant dans les domaines des statistiques, de l'apprentissage automatique et de l'apprentissage profond. Dans cet article, nous nous concentrons sur l'optimisation parcimonieuse exacte, en intégrant la régularisation $L_0$ dans notre formulation. Nous introduisons une relaxation renforcée issue d'une convexification partielle dans l'espace original et proposons un algorithme général d'optimisation par plans coupants, adapté aux problèmes d'optimisation quadratique en nombres entiers mixtes. Notre méthode est évaluée sur des données simulées pour lesquelles des résultats sont présentés.

**Mots-clés.** Optimisation Parcimonieuse, Norme $L_0$, Optimisation Quadratique, Programmation en Nombres Entiers, Plans Coupants.

**Abstract.** The sparse optimization problem, which aims at addressing the feature selection problem, has garnered increasing attention in the fields of statistics, machine and deep learning. In this paper, we focus on exact sparse optimization by incorporating $L_0$ regularization into our formulation. We introduce a strong relaxation derived from partial convexifications in the original space and propose a general cutting-plane optimization algorithm tailored for quadratic mixed-integer optimization problems. Our algorithm is evaluated on simulated data, and the experimental results are presented.

**Keywords.** Feature Selection, $L_0$ - norm, Quadratic Optimization, Integer Programming, Cutting-plane Methods.

# 1 Introduction

The mean square error is one of the most widely used error function for measuring the difference between a regressor and a response vector. When facing high-dimensional data, obtaining sparse regressors is crucial as it helps to lower the dimension of the data, and gives steep insights on its structure, making it more interpretable.

<footnote n="1">¹ Université Clermont-Auvergne, CNRS, Mines de Saint-Etienne, Clermont-Auvergne-INP, LIMOS, France renaud.chicoisne@uca.fr</footnote>

<footnote n="2">² Université Clermont-Auvergne, CNRS, laboratoire LMBP, France, pierre.latouche@uca.fr</footnote>

<footnote n="3">³ Université Clermont-Auvergne, CNRS, Mines de Saint-Etienne, Clermont-Auvergne-INP, LIMOS, France rui.ma@doctorant.uca.fr</footnote>

<footnote n="4">⁴ Institut Universitaire de France (IUF)</footnote>
In practice, to get an approximate sparse solution in reasonable time, the most common way to get a sparse regressor is to apply heuristic methods based on convex - and thus tractable - approximations. In the works of Latouche et al. (2016); Bouveyron et al. (2018, 2020), a simple relaxation along with an expectation-maximization algorithm were employed for inference. Besides, Lasso (Tibshirani, 1996) and Elastic net (Zou and Hastie, 2005) regularizations convexly penalize dense potential regressors. However, such methods can lead to biased regression results (Bertsimas et al., 2016). In contrast, the $L_0$ regularization - which directly controls the sparsity - is less used in the statistical community. More formally, given a set of $m$ observations in an $n$-dimensional feature space represented by an $m \times n$ matrix $X$, a response vector $y \in \mathbb{R}^m$ and a parameter $\alpha > 0$, the $L_0$-regularized regression problem can be written as follows:

$$ \min_{\beta} \frac{1}{2} \|X\beta - y\|_2^2 + \alpha \|\beta\|_0, \qquad (1) $$

where the (nonconvex) $L_0$ term $\| \cdot \|_0$ “counts” the number of nonzero coefficients of the vector it takes as an argument. It is well-established that solving such a mixed-integer optimization problem (MIP) is NP-hard (Natarajan, 1995). In the optimization community, two main categories of algorithms are used to solve MIPs: exact and heuristic methods. Exact methods leverage the computational power of modern optimization solvers, solving MIPs through branch and bound algorithms. Additionally, when the matrix $Q = X^\top X$ is sparse, exact methods (Bertsimas et al., 2016; Hastie et al., 2015; Wainwright, 2009) tend to perform better. On the other hand, heuristic methods do not guarantee an optimal regressor but are significantly faster than exact methods. In practice, well-implemented packages (Friedman et al., 2023; Hazimeh and Mazumder, 2020) are readily available and can give reasonable regressors $\beta$, having, however, no guarantee of being optimal solutions to (1).

In this work, we focus on the exact solution of (1): Defining $a := -X^\top y$ and $b := \alpha \cdot 1$, it can be rewritten - up to a constant - as the following extended MIP:

$$ \begin{array}{ll} \min_{\beta,z} & \frac{1}{2}\beta^\top Q\beta + a^\top\beta + b^\top z \\ \text{s.t.} & \beta_j(1 - z_j) = 0, \quad z_j \in \{0, 1\}, \quad \forall j \in [n], \end{array} \qquad (2) $$

where $[n] = \{1, \dots, n\}$. The complementarity constraints $\beta_j(1 - z_j) = 0$ can be enforced using big-M constraints, which - unfortunately- are notoriously hard to deal with in MIPs, as they induce serious numerical issues and provide extremely weak relaxation bounds. A classic approach to strengthen big-M formulations consists in generating strong cutting planes and reoptimize until we obtain a satisfying bound. This approach is unfortunately not adequate here, as reoptimizing nonlinear MIPs such as (2) is not efficient in practice. On the other hand, linear MIPs are particularly easy to reoptimize: In this paper, we build up on two recent results of Wei et al. (2024): 1) their reformulation (MILO) of (2) that is even weaker than classical nonlinear big-M formulations but is a linear MIP, and 2) the computationally expensive convexification that led to MILO to generate strong cutting planes.

In Section 2 we generalize their convexification to lower dimensions to obtain a tradeoff between strength and computational efficiency and then use it to generate strong cutting planes for the linear MIP MILO, which can be reoptimized efficiently after being strengthened
by these new cuts. We present in Section 3 a comparison of the relaxation bounds obtained by our approach, the classical nonlinear big-M based MIP and MILO, on medium-sized synthetic instances. If this paper is accepted, extra results on large instances and on real data will be presented. In this work, we focus exclusively on the case where $Q$ is positive definite. To secure this property, we incorporate an $L_2$ regularization term with a coefficient of 1 in the objective function.

## 2 Optimization

The main difficulty of (2) resides in the fact that the following set is nonconvex:

$$ \mathcal{S} := \{(\beta, z, t) \in \mathbb{R}^n \times \{0, 1\}^n \times \mathbb{R} : t \ge \beta^\top Q \beta, \quad \beta_j(1 - z_j) = 0, \forall j \in [n]\}, $$

making (2) equivalent to $\min_{t, \beta, z} \{t/2 + a^\top \beta + b^\top z : (\beta, z, t) \in \mathcal{S}\}$. As optimizing a linear functional over a set is equivalent to do so over its convex hull, any information of said convex hull is of paramount importance: In Wei et al. (2024), they prove that $(\beta, z, t) \in \text{conv} \mathcal{S}$ if and only if there exists a symmetric matrix $W \in \mathbb{S}^n$ such that:

$$ \begin{gathered} \begin{pmatrix} t & \beta^\top \\ \beta & W \end{pmatrix} \in \mathbb{S}_+^{n+1} \\ (z, W) \in \text{conv} \left\{ (z^\mathcal{N}, Q^\mathcal{N}), \mathcal{N} \subset [n] \right\} \end{gathered} \quad (3) $$

where $\mathbb{S}_+^{n+1}$ is the cone of $(n+1) \times (n+1)$ symmetric semidefinite positive matrices, each vector $z^\mathcal{N} \in \{0, 1\}^n$ is the indicator of $\mathcal{N}$ and each $n \times n$ matrix $Q^\mathcal{N}$ is the inverse of the submatrix of $Q$ restricted to $\mathcal{N}$, padded with zeroes everywhere else.

**Partial Convexification and Separation** While the latter result gives a nice mathematical representation of $\text{conv} \mathcal{S}$, it is intractable in practice due to the presence of the convex hull (3). We now propose to partially convexify $\mathcal{S}$ in the following Theorem, for which the proof is left to the full length paper:

**Theorem 1.** Given $\mathcal{J} \subset [n]$ of cardinality $q$, $(\beta, z, t)$ belongs to the convex hull of

$$ \mathcal{S}_{\mathcal{J}} := \{(\beta, z, t) \in \mathbb{R}^n \times \{0, 1\}^n \times \mathbb{R} : t \ge \beta^\top Q \beta, \quad \beta_j(1 - z_j) = 0, \forall j \in \mathcal{J}\} $$

if and only if, there exists $\mu \in \mathbb{R}$ and a matrix $W \in \mathbb{S}^q$ such that

$$ \begin{gathered} t \ge \beta^\top \begin{pmatrix} T^\top S^{-1} T & T^\top \\ T & S \end{pmatrix} \beta + \mu, \quad \begin{pmatrix} \mu & \beta(\mathcal{J})^\top \\ \beta(\mathcal{J}) & W \end{pmatrix} \in \mathbb{S}_+^{q+1} \\ (z(\mathcal{J}), W) \in \text{conv} \left\{ (z^\mathcal{N}, (P - T^\top S^{-1} T)^\mathcal{N}), \mathcal{N} \subset \mathcal{J} \right\} \end{gathered} \quad (4) $$

where $z(\mathcal{J})$ and $\beta(\mathcal{J})$ are respectively the vectors $z$ and $\beta$ restricted to their components in $\mathcal{J}$ and $Q$ is block-decomposed as follows, where wlog the first rows and columns correspond to the indices in $\mathcal{J}$:

$$ Q = \begin{pmatrix} P & T^\top \\ T & S \end{pmatrix} $$
Theorem 1 allows us to partially strengthen any formulation for our main problem (2). For reasonable sizes of $|\mathcal{J}| = q$, it is tractable to partially convexify $\mathcal{S}$, as (4) can be completely defined by enumerating all $2^q$ extreme points, as opposed to the $2^n$ in (3). By using duality results, we now present a way to separate an incumbent from conv $\mathcal{S}_J$:

**Corollary 1.** $(\beta, z, t) \in \text{conv} \mathcal{S}_J$ if and only if

$$ 
\begin{array}{l}
t \ge \beta^\top \begin{pmatrix} T^\top S^{-1}T & T^\top \\ T & S \end{pmatrix} \beta + \max_{\pi, \lambda, \lambda_0} \pi^\top \beta(\mathcal{J}) + \lambda^\top z(\mathcal{J}) + \lambda_0 \\
\text{s.t.} \quad \frac{1}{4} \pi^\top (P - T^\top S^{-1}T)^\mathcal{N} \pi + \lambda_0 + \sum_{j \in \mathcal{N}} \lambda_j \le 0, \quad \forall \mathcal{N} \subset \mathcal{J}
\end{array}
$$

Given a (fractional) incumbent $(\bar{\beta}, \bar{z}, \bar{t})$, we can generate a cut separating it from conv $\mathcal{S}_J$ by solving the (moderately large) maximization problem above, which can be done in polynomial time via interior point methods (Potra and Wright, 2000). Letting $(\pi^*, \lambda^*, \lambda_0^*)$ be its optimal solution, the following cut is valid for (2):

$$ t \ge \beta^\top \begin{pmatrix} T^\top S^{-1}T & T^\top \\ T & S \end{pmatrix} \beta + \beta(\mathcal{J})^\top \pi^* + z(\mathcal{J})^\top \lambda^* + \lambda_0^* $$

In practice, the selection of $\mathcal{J}$ is crucial. We discuss this in more detail in Section 3. As mentioned earlier, while these cuts are valid, adding them iteratively and reoptimize many times cannot be done efficiently in this *nonlinear* setting. We now show that we can translate these partial enumeration cuts for the *linear* MILO formulation from Wei et al. (2024).

**Linear Reformulation and Equivalent Cuts** Wei et al. (2024) recast the strong, but exponential-size representation (3) into the following weak, but polynomial-size linear MIP that is also devoid of variables $\beta$:

$$ 
\begin{array}{ll}
\min_{z,W,\eta} & \eta \\
\text{s.t.} & \eta \ge -\frac{1}{2} a^\top W a + b^\top z \\
& \displaystyle \sum_{k=1}^{n} Q_{ik} W_{ki} = z_i, \qquad \forall i \in [n] \quad (5) \\
(\text{MILO}) & -M(1-z_i) \le \sum_{k=1}^{n} Q_{ik} W_{kj} \le M(1-z_i), \qquad \forall i \ne j \quad (6) \\
& -\lambda_{\max}(Q^{-1})z_i \le W_{ij} \le \lambda_{\max}(Q^{-1})z_i, \qquad \forall (i,j) \in [n]^2 \quad (7) \\
& -\lambda_{\max}(Q^{-1})z_j \le W_{ij} \le \lambda_{\max}(Q^{-1})z_j, \qquad \forall (i,j) \in [n]^2 \quad (8) \\
& z \in \{0,1\}^n, \quad W \in \mathbb{S}^n
\end{array}
$$

Where $M := \lambda_{\max}(Q^{-1}) \cdot \max_{i \in [n]} ||Q_{[n],i}||_2$. Importantly, it is shown that the constraints (5)-(9) precisely define the exponentially defined set (3), and that we can retrieve the optimal regressor with $\beta^* = -W^*a$ from an optimal solution $(z^*, W^*, \eta^*)$ of MILO. We now state our main result, which is the direct translation of the *nonlinear* separation problem from Corollary 1, now in terms of the MILO formulation:
**Corollary 2.** For any set $\mathcal{J} \subset [n]$, the following inequality is valid for MILO:

$$ \eta \geq -\frac{1}{2}\bar{a}^{\top}S^{-1}\bar{a} + \max_{\lambda, \lambda_0} \lambda^{\top}z(\mathcal{J}) + \lambda_0 \\ \text{s.t. } \lambda_0 + \sum_{j \in \mathcal{J}, z_j^{\mathcal{N}}=1} \lambda_j + \frac{1}{2}(\hat{a} - T^{\top}S^{-1}\bar{a})^{\top}\hat{W}^{\mathcal{N}}(\hat{a} - T^{\top}S^{-1}\bar{a}) \leq 0, \forall \mathcal{N} \subset \mathcal{J}, $$

where $\hat{a} := a(\mathcal{J})$, $\bar{a} := a(\mathcal{J} \setminus [n])$, $\hat{W} := P - T^{\top}S^{-1}T$ and $P, T, S$ the blocks of $Q$ wrt $\mathcal{J}$.

Notice that the maximization problem above is linear and generates the following cuts:

$$ \eta \geq -\frac{1}{2}\bar{a}^{\top}S^{-1}\bar{a} + z(\mathcal{J})^{\top}\lambda^{*} + \lambda_0^{*}. $$

# 3 Experiments on simulated data

We now compare the relaxation bounds obtained with the classic nonlinear big-M formulation and our MILO-based, cut-enhanced models: i.e. we compare the strength of the relaxations obtained by each model when the integrality of the $z \in \{0, 1\}^n$ is relaxed to $z \in [0, 1]^n$. We solve our models using CPLEX 12.10 (Cplex, 2009) with default settings on a desktop computer with Intel Core i5-10505 and 16 GB of RAM.

**Simulated Data** We simulated datasets $(X, y)$ with two generators $\mathbb{G}_1$ and $\mathbb{G}_2$. For each experiment size and generator, we sample 10 datasets with identical settings and record the expected value $(\mathbb{E})$ and standard deviation $(\sigma)$ of the improvement achieved (in %) by our methods against the classical nonlinear big-M model, as well as the running time. All the data sets are generated by sampling $\beta \sim \mathcal{N}(0_n, I_n/\tau_1)$ with a noise $\varepsilon \sim \mathcal{N}(0_m, I_m/\tau_2)$ and by computing $y = X\beta + \varepsilon$. The matrix $X$ itself is generated from a Gaussian distribution with zero mean and a covariance matrix $R$. For $\mathbb{G}_1$, $R$ is set to $R = \text{diag}(R^1, R^2)$, a 2-blocks diagonal matrix, where $R^l$ is such that $R_{ii}^l = 1$ and $R_{ij}^l = \rho$ for $i, j \in [n/2]$ and $i \neq j$. For $\mathbb{G}_2$, a matrix $\tilde{R}$ is first built with each element sampled randomly and independently. Second, each element is set to 10 with probability 1/2 and to 0.01 otherwise. $R$ is finally built by symmetrizing $\tilde{R}$ and by setting each diagonal element to $R_{ii} = \sum_{j=1}^n R_{ij} + 1$. For both $\mathbb{G}_1$ and $\mathbb{G}_2$, we set $\tau_1 = \tau_2 = 1$ and $\rho = 0.75$.

**Set $\mathcal{J}$ Selection** Given a fractional incumbent $\bar{z}$, we jointly used three heuristic choices for selecting the set $\mathcal{J}$ on which we will partially enumerate: 1) wrt the correlation of $Q$ in (2), 2) wrt the “fractionality” of each $\bar{z}_j$, or 3) a random selection of indices. We find out that in practice, our cuts are not always improving the bound. Let $T$ be the maximum number of non-contributing cuts that we apply before updating $\mathcal{J}$. We consider two macro-strategies for selecting $\mathcal{J}$: The first approach uses a fixed value of $q$, small enough so that $2^q$ is enumerable, for which we take $q \in \{2, ..., 6\}$ and the selection process for $\mathcal{J}$ terminates when $T$ reaches 20. The second approach allows for a dynamic adjustment of $q$. In this
case, $q$ is incremented by 1, until it reaches $q = 6$. At that point the cutting-plane algorithm terminates. We tested this second strategy with $T \in \{1, 10, 20\}$.

<table><caption>Table 1: Relaxation bound improvement (%) from adding cuts.</caption><thead><tr><th></th><th>(n, m)</th><th></th><th>q = 2</th><th>q = 3</th><th>q = 4</th><th>q = 5</th><th>q = 6</th><th>T = 5</th><th>T = 10</th><th>T = 20</th></tr></thead><tbody><tr><td rowspan="3">G1</td><td>(10,50)</td><td>E</td><td>0.03</td><td>2.46</td><td>4.27</td><td>6.84</td><td>8.21</td><td>7.43</td><td>7.85</td><td>8.05</td></tr><tr><td></td><td>σ</td><td>4.04</td><td>5.32</td><td>6.88</td><td>7.70</td><td>8.49</td><td>8.95</td><td>8.69</td><td>8.32</td></tr><tr><td>(20,100)</td><td>E</td><td>-2.51</td><td>-1.47</td><td>-0.6</td><td>0.6</td><td>1.2</td><td>0.73</td><td>0.79</td><td>0.83</td></tr><tr><td></td><td></td><td>σ</td><td>1.29</td><td>1.41</td><td>1.47</td><td>1.71</td><td>2.23</td><td>2.49</td><td>2.41</td><td>2.54</td></tr><tr><td rowspan="3">G2</td><td>(10,50)</td><td>E</td><td>10.01</td><td>15.23</td><td>15.36</td><td>17.58</td><td>22.57</td><td>14.55</td><td>22.24</td><td>24.03</td></tr><tr><td></td><td>σ</td><td>27.59</td><td>28.32</td><td>28.08</td><td>28.72</td><td>30.44</td><td>27.28</td><td>31.48</td><td>32.22</td></tr><tr><td>(20,100)</td><td>E</td><td>139.34</td><td>140.36</td><td>144.68</td><td>148.93</td><td>153.29</td><td>144.07</td><td>147.14</td><td>155.36</td></tr><tr><td></td><td></td><td>σ</td><td>233.81</td><td>236.42</td><td>237.41</td><td>247.67</td><td>253.58</td><td>234.47</td><td>243.30</td><td>255.04</td></tr></tbody></table>

The values (%) presented in Table 1 are calculated as $100 \times \frac{(\text{MILO\&CUTS bound})-(\text{BigM bound})}{(\text{BigM bound})}$. Unsurprisingly, we observe a general improvement of the relaxation bound as $q$ increases.

<table><caption>Table 2: Running times (in seconds) for each data setting.</caption><thead><tr><th></th><th>(n, m)</th><th></th><th>q = 2</th><th>q = 3</th><th>q = 4</th><th>q = 5</th><th>q = 6</th><th>T = 5</th><th>T = 10</th><th>T = 20</th></tr></thead><tbody><tr><td rowspan="3">G1</td><td>(10,50)</td><td>E</td><td>3.97</td><td>4.43</td><td>4.17</td><td>4.92</td><td>4.76</td><td>5.54</td><td>9.36</td><td>18.09</td></tr><tr><td></td><td>σ</td><td>0.95</td><td>1.12</td><td>0.91</td><td>1.85</td><td>1.68</td><td>0.31</td><td>0.83</td><td>3.65</td></tr><tr><td>(20,100)</td><td>E</td><td>54.30</td><td>53.06</td><td>51.49</td><td>49.21</td><td>51.38</td><td>59.77</td><td>94.1</td><td>157.19</td></tr><tr><td></td><td></td><td>σ</td><td>21.52</td><td>15.97</td><td>15.62</td><td>13.94</td><td>18.39</td><td>15.09</td><td>22.44</td><td>46.92</td></tr><tr><td rowspan="3">G2</td><td>(10,50)</td><td>E</td><td>4.54</td><td>4.54</td><td>4.73</td><td>4.20</td><td>3.02</td><td>6.48</td><td>14.49</td><td>11.72</td></tr><tr><td></td><td>σ</td><td>2.26</td><td>2.44</td><td>2.92</td><td>2.17</td><td>1.93</td><td>4.06</td><td>7.07</td><td>8.90</td></tr><tr><td>(20,100)</td><td>E</td><td>39.8</td><td>37.58</td><td>38.64</td><td>41.94</td><td>43.06</td><td>45.35</td><td>63.6</td><td>96.65</td></tr><tr><td></td><td></td><td>σ</td><td>12.62</td><td>16.78</td><td>14.66</td><td>12.6</td><td>21.34</td><td>17.31</td><td>22.12</td><td>30.81</td></tr></tbody></table>

As shown in Table 2, our cutting-plane method generally terminates in a reasonable amount of time. In particular, when compared to the time required to solve the the nonlinear big-M model, which takes an average of around 500 seconds, our method demonstrates a significant reduction in computation time.

# 4 Conclusions

In this paper, we present a partial convexification method for solving quadratic mixed-integer optimization problems, and used it to provide a new optimization scheme of the mean square error function under $L_0$ penalization. We also introduced a cutting-plane algorithm based on this method, which strengthens the linear reformulation MILO proposed by Wei et al. (2024) and yields favorable results compared to the conventional big-M formulation within a reasonable time. The work presented in this paper is part of an ongoing work. In the next phase, we will conduct experiments with larger values of $q$ and real-life data.
References

Dimitris Bertsimas, Angela King, and Rahul Mazumder. Best subset selection via a modern optimization lens. *The Annals of Statistics*, 44(2):813–852, 2016.

Charles Bouveyron, Pierre Latouche, and Pierre-Alexandre Mattei. Bayesian variable selection for globally sparse probabilistic pca. *Electronic Journal of Statistics*, 2018.

Charles Bouveyron, Pierre Latouche, and Pierre-Alexandre Mattei. Exact dimensionality selection for bayesian pca. *Scandinavian Journal of Statistics*, 47(1):196–211, 2020.

IBM ILOG Cplex. V12. 1: User's manual for CPLEX. *International Business Machines Corporation*, 46(53):157, 2009.

Jerome Friedman, Trevor Hastie, Rob Tibshirani, Balasubramanian Narasimhan, Kenneth Tay, Noah Simon, Junyang Qian, and James Yang. glmnet: Lasso and elastic-net regularized generalized linear models. *Astrophysics Source Code Library*, pages ascl-2308, 2023.

Trevor Hastie, Robert Tibshirani, and Martin Wainwright. Statistical learning with sparsity. *Monographs on statistics and applied probability*, 143(143):8, 2015.

Hussein Hazimeh and Rahul Mazumder. Fast best subset selection: Coordinate descent and local combinatorial optimization algorithms. *Operations Research*, 68(5):1517–1537, 2020.

Pierre Latouche, Pierre-Alexandre Mattei, Charles Bouveyron, and Julien Chiquet. Combining a relaxed EM algorithm with occam's razor for bayesian variable selection in high-dimensional regression. *Journal of Multivariate Analysis*, 146:177–190, 2016.

Balas Kausik Natarajan. Sparse approximate solutions to linear systems. *SIAM journal on computing*, 24(2):227–234, 1995.

Florian A Potra and Stephen J Wright. Interior-point methods. *Journal of computational and applied mathematics*, 124(1-2):281–302, 2000.

Robert Tibshirani. Regression shrinkage and selection via the Lasso. *Journal of the Royal Statistical Society Series B: Statistical Methodology*, 58(1):267–288, 1996.

Mark Wainwright. *Photosensitisers in biomedicine*. John Wiley & Sons, 2009.

Linchuan Wei, Alper Atamtürk, Andrés Gómez, and Simge Küçükyavuz. On the convex hull of convex quadratic optimization problems with indicators. *Mathematical Programming*, 204(1):703–737, 2024.

Hui Zou and Trevor Hastie. Regularization and variable selection via the elastic net. *Journal of the Royal Statistical Society Series B: Statistical Methodology*, 67(2):301–320, 2005.