# COMBSS: Best Subset Selection via Continuous Optimization

Sarat Moka1,∗, Benoit Liquet1, Houying Zhu1, and Samuel Muller1,2

1 School of Mathematical and Physical Sciences, Macquarie University, Sydney,
Australia
2 School of Mathematics and Statistics, University of Sydney, Sydney,
Australia
∗ Corresponding author<EMAIL_ADDRESS>
###### Abstract

We consider the problem of best subset selection in linear regression, where
the goal is to find for every model size $k$, that subset of $k$ features that
best fit the response. This is particularly challenging when the total
available number of features is very large compared to the number of data
samples. We propose COMBSS, a novel continuous optimization based method that
identifies a solution path, a small set of models of varying size, that
consists of candidates for the best subset in linear regression. COMBSS turns
out to be very fast, making subset selection possible when the number of
features is well in excess of thousands. Simulation results are presented to
highlight the performance of COMBSS in comparison to existing popular methods
such as Forward Stepwise, the Lasso and Mixed-Integer Optimization. Because of
the outstanding overall performance, framing the best subset selection
challenge as a continuous optimization problem opens new research directions
for feature extraction for a large variety of regression models.

##  1 Introduction

jrssRecent rapid developments in information technology have enabled the
collection of high-dimensional complex data, including in engineering,
economics, finance, biology, and health sciences [8]. High-dimensional means
that the number of features is large and often far higher than the number of
collected data samples. In many of these applications, it is desirable to find
a small best subset of predictors so that the resulting model has desirable
prediction accuracy [9, 21]. In this article, we recast the challenge of best
subset selection in linear regression as a continuous optimization problem. We
show that this reframing has enormous potential and substantially advances
research into larger dimensional and exhaustive feature selection in
regression, making available technology that can reliably select significant
variables when $p$ is well in excess of 1000s.

Consider the linear regression model
$\mathbf{y}=X\boldsymbol{\beta}+\boldsymbol{\epsilon},$ where
$\mathbf{y}=(y_{1},\dots,y_{n})^{{}^{\mkern-1.5mu\mathsf{T}}}$ is an
$n$-dimensional known response vector, $X$ is a known design matrix of
dimension $n\times p$ with $x_{i,j}$ indicating the $i$-th observation of the
$j$-th explanatory variable,
$\boldsymbol{\beta}=(\beta_{1},\dots,\beta_{p})^{{}^{\mkern-1.5mu\mathsf{T}}}$
is the $p$-dimensional vector of unknown regression coefficients, and
$\boldsymbol{\epsilon}=(\epsilon_{1},\dots,\epsilon_{n})^{{{}^{\mkern-1.5mu\mathsf{T}}}}$
is a vector of unknown errors, unless otherwise specified, assumed to be
independent and identically distributed (iid).

Best subset selection is a classical problem that aims to solve,

$\displaystyle\begin{aligned}
\operatorname*{minimize\,}_{\boldsymbol{\beta}\in\mathbb{R}^{p}}\,\,\frac{1}{n}\|\mathbf{y}-X\boldsymbol{\beta}\|_{2}^{2}\quad\operatorname*{subject\,\,to\,}\|\boldsymbol{\beta}\|_{0}=k,\end{aligned}$
(1)

for a given $k$, where $\|\cdot\|_{2}$ is the usual $\mathcal{L}_{2}$-norm,
and $\|\boldsymbol{\beta}\|_{0}=\sum_{j=1}^{p}I(\beta_{j}\neq 0)$, the number
of non-zero elements in $\boldsymbol{\beta}$. For ease of reading, in this
article, we assume that all columns of $X$ are subject to selection.

In many practical settings, it is assumed that only a small number of
explanatory variables contribute in the data generating model to the response,
say there are at most $q\ll\min(n,p)$ non-zero regression coefficients. Thus,
best subset selection aims to solve (1) for every $k\leq q$ to select indices
of the non-zero regression coefficients while taking all other regression
coefficients to be zero.

When the additional aim is to select the best model among all the best subset
models of size $k=0,\ldots,q$, this is often achieved by using some
performance measure such as criteria that minimize the prediction error or
information criteria. Examples for the latter include the AIC [1] and BIC
[24], we refer to [22] for a review on information criteria in linear
regression.

Exact methods for solving (1) are typically executed by first writing
solutions for low-dimensional problems and then selecting the best subset
solution over these low dimensional problems. To see this, for any binary
vector
$\mathbf{s}=(s_{1},\dots,s_{p})^{{}^{\mkern-1.5mu\mathsf{T}}}\in\\{0,1\\}^{p}$,
let $X_{[\mathbf{s}]}$ be the matrix of size $n\times|\mathbf{s}|$ created by
keeping only columns $j$ of $X$ for which $s_{j}=1$, where $j=1,\ldots,p$.
Then, for every $k\leq q$, in exact best subset selection, the optimal
$\mathbf{s}$ can be found by solving the problem,

$\displaystyle\begin{aligned}
\operatorname*{minimize\,}_{{\mathbf{s}}\in\\{0,1\\}^{p}}\,\,\frac{1}{n}\|\mathbf{y}-{X}_{[\mathbf{s}]}\widehat{\boldsymbol{\beta}}_{[\mathbf{s}]}\|_{2}^{2}\quad\operatorname*{subject\,\,to\,}|{\mathbf{s}}|=k,\end{aligned}$
(2)

where $\widehat{\boldsymbol{\beta}}_{[\mathbf{s}]}$ is the low-dimensional
least squares estimate of elements of $\boldsymbol{\beta}$ with indices
corresponding to non-zero elements of $\mathbf{s}$, i.e.

$\displaystyle\widehat{\boldsymbol{\beta}}_{[\mathbf{s}]}=(X_{[\mathbf{s}]}^{{}^{\mkern-1.5mu\mathsf{T}}}{X}_{[\mathbf{s}]})^{-1}{X}_{[\mathbf{s}]}^{{}^{\mkern-1.5mu\mathsf{T}}}\mathbf{y},$
(3)

whenever the matrix inverse
$({X}_{[\mathbf{s}]}^{{}^{\mkern-1.5mu\mathsf{T}}}{X}_{[\mathbf{s}]})^{-1}$
exists, which is possible when the matrix
${X}_{[\mathbf{s}]}^{{}^{\mkern-1.5mu\mathsf{T}}}{X}_{[\mathbf{s}]}$ is of
full rank $k$. The solution of (2) is identical to the solution of (1),
because $\widehat{\boldsymbol{\beta}}_{[\mathbf{s}]}$ is the least squares
solution when constrained so that $I(\beta_{j}\neq 0)=s_{j}$ for all
$j=1,\dots,p$.

It is well-known that solving the exact optimization problem in (1) is NP-hard
[23] in general. For instance, a popular exact method called leaps and bounds
[10] is practically useful only for values of $p$ smaller than $30$, as for
example is the case when using the R package leaps [19] or mplot [26].

In this paper, we design COMBSS, a novel continuous optimization method
towards best subset selection. To the best of our knowledge, COMBSS is the
first attempt to exploit continuous optimization for finding significant
variables in regression for high-dimensional data.

In a nutshell, COMBSS can be described as follows: Instead of the binary
vector space $\\{0,1\\}^{p}$ as in the exact methods, we consider the hyper-
cube $[0,1]^{p}$ and define for each $\mathbf{t}\in[0,1]^{p}$ a new estimate
$\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}$ so that we have the following
well-defined continuous extension of the exact problem (2):

$\displaystyle\operatorname*{minimize\,}_{{\mathbf{t}}\in[0,1]^{p}}\,\,\frac{1}{n}\|\mathbf{y}-X_{\mathbf{t}}\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}\|_{2}^{2}+\lambda\sum_{j=1}^{p}t_{j},$
(4)

where $X_{\mathbf{t}}$ is obtained from $X$ by multiplying the $j$-th column
of $X$ by $t_{j}$ for all $j=1,\dots,p$, and the tuning parameter $\lambda$
controls the sparsity of the solution obtained, analogous to selecting the
best $k\leq q$ in the exact optimization.

Although, COMBSS aims to find sets of models that are candidates for the best
subset of variables, an important property is that it has no discrete
constraints, unlike the exact optimization problem (2). As a consequence,
COMBSS can take advantage of standard continuous optimization methods, such as
the basic gradient descent method, by starting at an interior point on the
hypercube $[0,1]^{p}$ and iteratively moving towards a corner that minimizes
the objective function. In the implementation, we move the box constrained
problem (4) to an equivalent unconstrained problem so that the gradient
descent method can run without any boundary issues.

The following list provides some key features of COMBSS that make it a novel,
fast, and efficient approach towards finding good and/or best subsets:

  * •

Unlike in the exact method, where $X_{[\mathbf{s}]}$ is constructed by
removing columns from the design matrix $X$ (and hence $X_{[\mathbf{s}]}$ and
$X$ are of different sizes), $X_{\mathbf{t}}$ in COMBSS is constructed by
multiplying the $j$-th column of $X$ by $t_{j}$ for every $j$. Thus, both
$X_{\mathbf{t}}$ and $X$ are of the same size.

(a)

(b)

Figure 1:
$\frac{1}{n}\|\mathbf{y}-{X}_{[\mathbf{s}]}\widehat{\boldsymbol{\beta}}_{[\mathbf{s}]}\|_{2}^{2}$
vs
$\frac{1}{n}\|\mathbf{y}-X_{\mathbf{t}}\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}\|_{2}^{2}$
for the data used later in Example 1. Intuitively, if we roll a ball (i.e.,
running a gradient descent algorithm) on the surface of
$\frac{1}{n}\|\mathbf{y}-X_{\mathbf{t}}\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}\|_{2}^{2}$
shown in Panel (b), it will reach to the minimum point, which is
$(1,1)^{{}^{\mkern-1.5mu\mathsf{T}}}$ for this example. By adding the penalty
$\lambda\sum_{j=1}^{p}t_{j}$ and changing $\lambda$, we can control the path
of the ball so that it reaches a sparse $\mathbf{t}$.

  * •

Our construction of the new coefficient estimator
$\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}$ guarantees that at the corner
points $\mathbf{s}$ of the hyper-cube,
$\|\mathbf{y}-X_{\mathbf{s}}\widetilde{\boldsymbol{\beta}}_{\mathbf{s}}\|_{2}=\|\mathbf{y}-{X}_{[\mathbf{s}]}\widehat{\boldsymbol{\beta}}_{[\mathbf{s}]}\|_{2}$,
where $\widehat{\boldsymbol{\beta}}_{[\mathbf{s}]}$ exists, while making sure
that the new objective function
$\frac{1}{n}\|\mathbf{y}-X_{\mathbf{t}}\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}\|_{2}^{2}$
has smooth transition from one corner to another. For instance, see Figure 1
for an illustration.

  * •

By increasing the value of $\lambda$, we can increase the sparsity of the
solution of the optimization problem (4), because the penalty
$\lambda\sum_{j=1}^{p}t_{j}$ encourages sparsity in $\mathbf{t}$.

  * •

Solving the best subset selection problem using an exhaustive search is
currently infeasible when $p$ is of a moderate size as the number of visited
models increases exponentially. However, COMBSS enables us to visit this huge
space by exploiting the continuity of the objective function. Thus, a gradient
descent algorithm can converge in a few iterations to the solution. While
single runs of COMBSS do not necessarily identify the best subset, in our
empirical research we have found that it finds very competitive solution
paths.

  * •

We show that during the execution of the gradient descent, if some $t_{j}$
converges to a small value close to zero, the algorithm pushes $t_{j}$ to be
zero. We also show that once $t_{j}$ becomes zero, the corresponding $j$-th
column of the design matrix play no role in the gradient descent update. As a
consequence of these two observations, we propose some algorithmic
improvements that provide further substantial speed improvements to COMBSS.

We provide sample code that makes available COMBSS in Python and R at

  * •

https://github.com/saratmoka/COMBSS-Python-VIGNETTE and

  * •

https://github.com/benoit-liquet/COMBSS-R-VIGNETTE, respectively.

The code includes examples where $p$ is as large as of order $10,000$. This
code further allows to replicate our empirical research in this article. In
ongoing technological development, we are currently working towards adding
further functionality to this sample code. We aim to release COMBSS as an R
software package on CRAN and as a Python software package on PyPI as soon as
possible. We will update the above links accordingly when that work is
finalised and will also add references to follow-up papers, where more details
are provided on the computational and other aspects of COMBSS.

Popular Existing Methods: Even though the best subset selection problem is
classical, for many years, it is generally believed that the problem is not
feasible for large values of $p$. Relatively recent work [4] formulates (1) as
a mixed integer optimization (MIO) problem and demonstrates that the best
subset selection problem can be solved for much larger $p$ using modern MIO
solvers such as in the commercial software Gurobi [12], which is not freely
available (except for an initial short period). In particular, the best subset
selection problem is reformulated as the following MIO problem, containing
both continuous constraints and discrete constraints:

$\displaystyle\begin{aligned}
&\operatorname*{minimize\,}_{\boldsymbol{\beta},\,\mathbf{z}\in\mathbb{R}^{p}}\,\,\frac{1}{n}\|\mathbf{y}-X\boldsymbol{\beta}\|_{2}^{2}\\\
&\operatorname*{subject\,\,to\,}-b\,z_{j}\leq\beta_{j}\leq b\,z_{j},\quad
j=1,\dots,p,\\\ &\hskip 54.06006ptz_{j}\in\\{0,1\\},\quad j=1,\dots,p,\\\
&\hskip 54.06006pt\sum_{j=1}^{p}z_{j}=k,\end{aligned}$ (5)

where $b$ is a constant such that if $\boldsymbol{\beta}^{*}$ minimizes
(LABEL:eqn:MIO_general), then $|\beta^{*}_{j}|\leq b$ for all $j=1,\dots,p$.

Despite the fact that the MIO approach directly aims to solve the best subset
selection problem and is faster than the exact method for larger values of
$p$, its implementation via Gurobi is significantly slow from a practical
point of view. As a result, some of the indirect methods, such as Forward
Stepwise selection and the least absolute shrinkage and selection operator
(Lasso) and their variants, are still very common in practice. Forward
Stepwise selection follows a greedy approach which, starting with an empty
model (or intercept only model), iteratively adds the variable that is best
suited [7, 16]. On the other hand, instead of solving the highly non-convex
best subset selection problem, Lasso introduced in [28] solves a convex
relaxation of it by replacing the (discrete) $\mathcal{L}_{0}$ norm
$\|\boldsymbol{\beta}\|_{0}$ in (1) with the $\mathcal{L}_{1}$ norm
$\|\boldsymbol{\beta}\|_{1}$. Due to this clever relaxation, Lasso is a convex
optimization problem. Thus, it is fast (and much faster than the MIO solvers).
However, it is known that Lasso solutions do not normally lead to a best
subset [15, 30] solution.

The rest of the paper is organized as follows: In Section 2, we describe the
mathematical framework required for stating the proposed optimization method
COMBSS. In Section 3, we first establish the continuity of the objective
functions involved in COMBSS, and then we derive expressions for their
gradients. These gradients are later exploited for conducting continuous
optimization. A detailed algorithm for COMBSS is given in Section 4. In
Section 5, we discuss how to choose the tuning parameters that control the
surface shape of the objective functions and the sparsity of the solutions
obtained. Section 6 provides steps for efficient implementation of COMBSS
using some linear algebra techniques. In Section 7 we present a simulation
study to show the good behaviour and the performance of COMBSS compared to
traditional methods in low and high dimensional context. We conclude the paper
in Section 8 with a discussion on possible future improvements for COMBSS.

##  2 Continuous Extension of the Best Subset Selection Problem

In this section, we provide the mathematical setup for COMBSS, which solves a
continuous extension of the exact best subset selection optimization problem
(2).

For $\mathbf{t}=(t_{1},\dots,t_{p})^{{}^{\mkern-1.5mu\mathsf{T}}}\in[0,1]^{p}$
define $T_{\mathbf{t}}=\mathsf{Diag}(\mathbf{t})$, the diagonal matrix with
the diagonal elements being $t_{1},\dots,t_{p}$, and

$X_{\mathbf{t}}=XT_{\mathbf{t}}.$

For a fixed constant $\delta>0$, define

$\displaystyle
L_{\mathbf{t}}=L_{\mathbf{t}}(\delta)=\frac{1}{n}\left[X_{\mathbf{t}}^{{}^{\mkern-1.5mu\mathsf{T}}}X_{\mathbf{t}}+\delta\left(I-T_{\mathbf{t}}^{2}\right)\right]$
(6)

and we suppress $\delta$ for ease of reading. Intuitively, $L_{\mathbf{t}}$
can be seen as a ‘convex combination’ of the matrices
$(X^{{}^{\mkern-1.5mu\mathsf{T}}}X)/n$ and $\delta I/n$, because
$X_{\mathbf{t}}^{{}^{\mkern-1.5mu\mathsf{T}}}X_{\mathbf{t}}=T_{\mathbf{t}}{X^{{}^{\mkern-1.5mu\mathsf{T}}}X}T_{\mathbf{t}}$
and thus

$L_{\mathbf{t}}=T_{\mathbf{t}}\left(\frac{X^{{}^{\mkern-1.5mu\mathsf{T}}}X}{n}\right)T_{\mathbf{t}}+\left(I-T_{\mathbf{t}}\right)\,\left(\frac{\delta}{n}I\right)\,\left(I-T_{\mathbf{t}}\right).$

Now let

$\displaystyle\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}$
$\displaystyle:=L_{\mathbf{t}}^{+}\left(\frac{X_{\mathbf{t}}^{{}^{\mkern-1.5mu\mathsf{T}}}\mathbf{y}}{n}\right),\quad\mathbf{t}\in[0,1]^{p},$
(7)

where $L^{+}_{\mathbf{t}}$ is the pseudo-inverse (or, the Moore-Penrose
inverse) of $L_{\mathbf{t}}$. We need $L_{\mathbf{t}}^{+}$ in (7) so that
$\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}$ is defined for all
$\mathbf{t}\in[0,1]^{p}$. However, from the way we conduct optimization, we
need to compute $\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}$ only for
$\mathbf{t}\in[0,1)^{p}$. We later show in Theorem 1 that for all
$\mathbf{t}\in[0,1)^{p}$, $L_{\mathbf{t}}$ is invertible and thus in the
implementation of COMBSS, $\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}$ always
takes the form
$\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}=L_{\mathbf{t}}^{-1}\left(X_{\mathbf{t}}^{{}^{\mkern-1.5mu\mathsf{T}}}\mathbf{y}/n\right)$,
eliminating the need to compute any computationally expensive pseudo-inverse.

With the support of these observations, the following is an immediate well
defined generalization of best subset problem (1):

$\displaystyle\begin{aligned}
\operatorname*{minimize\,}_{\mathbf{t}\in[0,1]^{p}}\,\,\frac{1}{n}\|\mathbf{y}-X_{\mathbf{t}}\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}\|_{2}^{2}\quad\operatorname*{subject\,\,to\,}\sum_{j=1}^{p}t_{j}\leq
q.\end{aligned}$ (8)

The optimization problem (8) is a constrained problem. Instead of solving
this, by defining a Lagrangian function

$\displaystyle f_{\lambda}(\mathbf{t})$
$\displaystyle=\frac{1}{n}\|\mathbf{y}-X_{\mathbf{t}}\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}\|_{2}^{2}+\lambda\sum_{j=1}^{p}t_{j},$
(9)

for a tunable parameter $\lambda>0$, we aim to solve instead

$\displaystyle\operatorname*{minimize\,}_{\mathbf{t}\in[0,1]^{p}}f_{\lambda}(\mathbf{t}).$
(10)

We can reformulate the box constrained problem (10) into an equivalent
unconstrained problem by considering the mapping
$\mathbf{t}=\mathbf{t}(\mathbf{w})$ given by

$\displaystyle t_{j}(w_{j})=1-\exp(-w_{j}^{2}),\quad j=1,\dots,p,$ (11)

and then rewriting (10) as,

$\displaystyle\operatorname*{minimize\,}_{\mathbf{w}\in\mathbb{R}^{p}}f_{\lambda}\left(\mathbf{t}(\mathbf{w})\right).$
(12)

The unconstrained optimization problem (12) is equivalent to the box
constrained problem (10), because for any $u,v\in\mathbb{R}$,

$1-\exp(-u^{2})<1-\exp(-v^{2})\quad\text{if and only if}\quad u^{2}<v^{2}.$

In summary, by defining,

$\displaystyle
g_{\lambda}(\mathbf{w})=f_{\lambda}\left(\mathbf{t}(\mathbf{w})\right),$ (13)

we solve,

$\displaystyle\operatorname*{minimize\,}_{\mathbf{w}\in\mathbb{R}^{p}}g_{\lambda}(\mathbf{w}),$
(14)

via some continuous optimization method, such as a gradient descent method.

We recall that our ultimate goal is to find a solution of the best subset
selection problem (2). Suppose that $\mathbf{w}^{*}$ denotes an approximate
solution of (14) obtained by a gradient descent method. Then, the
corresponding solution $\mathbf{t}^{*}$ of (10) is obtained through the
mapping (11). Further using this $\mathbf{t}^{*}$, we select an appropriate
binary vector as a solution to (2); see Section 4 for complete details.

The following example shows the surface plots of the objective functions
$f_{\lambda}(\mathbf{t})$ and $g_{\lambda}(\mathbf{w})$ for a simple case.

######  Example 1.

We choose $p=2$ and $n=100$, and generate each row of $X$ from a multivariate
Gaussian with mean $\mu_{x}$ and covariance $\Sigma_{x}$ given by

$\mu_{x}=\begin{bmatrix}1\\\
1\end{bmatrix}\quad\text{and}\quad\Sigma_{x}=\begin{bmatrix}3&1\\\
1&2\end{bmatrix}.$

The noise $\epsilon$ is generated from a standard normal and the true model
has coefficient vector
$\boldsymbol{\beta}=(2,0)^{{}^{\mkern-1.5mu\mathsf{T}}}$. Suppose $k=1$ in the
best subset selection problem (1). Since $p$ is small, we can easily solve the
exact problem (2) to find its optimal solution
$\mathbf{s}^{*}=(1,0)^{{}^{\mkern-1.5mu\mathsf{T}}}$. Refer to Figure 2 to see
the surface plots of $f_{\lambda}(\mathbf{t})$ and $g_{\lambda}(\mathbf{w})$
for $\lambda=0$ and $\lambda=1.5$. Observe that when $\lambda$ is positive,
the surface of $f_{\lambda}(\mathbf{t})$ is exhibiting a downwards slope
towards $\mathbf{s}^{*}$ from all the directions. Similarly, for positive
$\lambda$, the surface of $g_{\lambda}(\mathbf{w})$ is exhibiting a valley
around the line $w_{2}=0$, and this helps the gradient descent algorithm to
traverse towards optimal points $(\pm\infty,0)^{{}^{\mkern-1.5mu\mathsf{T}}}$.

(a) $\lambda=0$

(b) $\lambda=0$

(c) $\lambda=1.5$

(d) $\lambda=1.5$

Figure 2: Surface plots of the objective functions $f_{\lambda}(\mathbf{t})$
and $g_{\lambda}(\mathbf{w})$ for $\lambda=0$ and $\lambda=1.5$. To see the
effect of changing $\lambda$, compare the values of $f_{\lambda}(\mathbf{t})$
in Panel (a) and Panel (c) at $\mathbf{t}=(1,0)^{{}^{\mkern-1.5mu\mathsf{T}}}$
and $\mathbf{t}=(1,1)^{{}^{\mkern-1.5mu\mathsf{T}}}$.

∎

##  3 Continuity and Gradients of the Objective Function

In this section, we first prove that the objective function
$g_{\lambda}(\mathbf{w})$, defined in (13), is continuous on $\mathbb{R}^{p}$
and then we derive its gradients. En-route, we also establish the relationship
between $\widehat{\boldsymbol{\beta}}_{\mathbf{s}}$ and
$\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}$ which are respectively defined
by (3) and (7). This relationship is useful in understanding the relationship
between our method COMBSS and the exact optimization (2). The detailed proofs
of all the results in this section are available in Appendix A.

Our first result, Theorem 1, shows that for all $\mathbf{t}\in[0,1)^{p}$, the
matrix $L_{\mathbf{t}}$, which is defined in (6), is symmetric positive-
definite and hence invertible.

######  Theorem 1.

For any $\mathbf{t}\in[0,1)^{p}$, $L_{\mathbf{t}}$ is symmetric and positive-
definite and

$\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}=L_{\mathbf{t}}^{-1}\left(\frac{X_{\mathbf{t}}^{{}^{\mkern-1.5mu\mathsf{T}}}\mathbf{y}}{n}\right).$

Recall that $\widehat{\boldsymbol{\beta}}_{[\mathbf{s}]}$ exists only on some
corner points $\mathbf{s}$ of the hyper-cube $[0,1]^{p}$, where the matrix
$X_{[\mathbf{s}]}^{{}^{\mkern-1.5mu\mathsf{T}}}X_{[\mathbf{s}]}$ is
invertible. Theorem 2 establishes a relationship between
$\widetilde{\boldsymbol{\beta}}_{\mathbf{s}}$ and
$\widehat{\boldsymbol{\beta}}_{[\mathbf{s}]}$ at these corner points
$\mathbf{s}\in\\{0,1\\}^{p}$.

######  Theorem 2.

For any point $\mathbf{s}\in\\{0,1\\}^{p}$,
$\widehat{\boldsymbol{\beta}}_{[\mathbf{s}]}$ exists if and only if
$L_{\mathbf{s}}$ is invertible and in that case,

$X_{[\mathbf{s}]}\widehat{\boldsymbol{\beta}}_{[\mathbf{s}]}=X_{\mathbf{s}}\widetilde{\boldsymbol{\beta}}_{\mathbf{s}}.$

As an immediate result of Theorem 2, we have

$f_{\lambda}(\mathbf{s})=\frac{1}{n}\|\mathbf{y}-X_{[\mathbf{s}]}\widehat{\boldsymbol{\beta}}_{[\mathbf{s}]}\|_{2}^{2}+\lambda|\mathbf{s}|,$

whenever $\widehat{\boldsymbol{\beta}}_{[\mathbf{s}]}$ exists. That is, the
objective function of the exact optimization problem (2) is identical to the
objective function of its extended optimization problem (8).

Our next result, Theorem 3, shows that $f_{\lambda}(\mathbf{t})$ is a
continuous function on $[0,1]^{p}$.

######  Theorem 3.

The function $f_{\lambda}(\mathbf{t})$ defined in (9) is continuous over
$[0,1]^{p}$ in the sence that for any sequence
$\mathbf{t}^{(1)},\mathbf{t}^{(2)},\dots\in[0,1)^{p}$ converging to
$\mathbf{t}\in[0,1]^{p}$, the limit
$\lim_{l\to\infty}f_{\lambda}(\mathbf{t}^{(l)})$ exists and

$f_{\lambda}(\mathbf{t})=\lim_{l\to\infty}f_{\lambda}(\mathbf{t}^{(l)}).$

The next result, Corollary 1, establishes the continuity of $g_{\lambda}$ on
$\mathbb{R}^{p}$. This is a simple consequence of Theorem 3, because from the
definition (13),
$g_{\lambda}(\mathbf{w})=f_{\lambda}\left(\mathbf{t}(\mathbf{w})\right)$ with
$\mathbf{t}(\mathbf{w})=\mathbf{1}-\exp(-\mathbf{w}\odot\mathbf{w})$, here
(and afterwards) $\mathbf{1}$ denotes a vector of all $1$’s of appropriate
dimension, $\odot$ denotes the element-wise product of two vectors, and the
function $\exp$ is applied element-wise.

######  Corollary 1.

The objective function $g_{\lambda}(\mathbf{w})$ is continuous at every point
$\mathbf{w}\in\mathbb{R}^{p}$.

As mentioned earlier, COMBSS uses a gradient descent method to solve the
problem (14). Towards that we need to obtain the gradients of
$g_{\lambda}(\mathbf{w})$. Theorem 4 provides an expression of the gradient
$\nabla g_{\lambda}(\mathbf{w})$ for every $\mathbf{w}\in\mathbb{R}^{p}$.

######  Theorem 4.

For every $\mathbf{w}\in\mathbb{R}^{p}$, with
$\mathbf{t}=\mathbf{t}(\mathbf{w})$ is given by (11),

$\nabla
f_{\lambda}(\mathbf{t})=\boldsymbol{\zeta}_{\mathbf{t}}+\lambda\mathbf{1},\quad\mathbf{t}\in(0,1)^{p},$

and

$\nabla
g_{\lambda}(\mathbf{w})=\left(\boldsymbol{\zeta}_{\mathbf{t}}+\lambda\mathbf{1}\right)\odot\left(2\mathbf{w}\odot\exp(-\mathbf{w}\odot\mathbf{w})\right),\quad\mathbf{w}\in\mathbb{R}^{p},$

where

$\displaystyle\boldsymbol{\zeta}_{\mathbf{t}}=2\left({\widetilde{\boldsymbol{\beta}}}_{\mathbf{t}}\odot\left({\mathbf{a}}_{\mathbf{t}}-{\mathbf{d}}_{\mathbf{t}}\right)\right)-2\left({\mathbf{b}}_{\mathbf{t}}\odot{\mathbf{c}}_{\mathbf{t}}\right),$
(15)

with

$\displaystyle\mathbf{a}_{\mathbf{t}}$
$\displaystyle=\left(\frac{X^{{}^{\mkern-1.5mu\mathsf{T}}}X}{n}\right)\left({\mathbf{t}}\odot\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}\right)-\left(\frac{X^{{}^{\mkern-1.5mu\mathsf{T}}}\mathbf{y}}{n}\right),$
$\displaystyle\mathbf{b}_{\mathbf{t}}$
$\displaystyle=\mathbf{a}_{\mathbf{t}}-\frac{\delta}{n}\left({\mathbf{t}}\odot\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}\right),$
$\displaystyle\mathbf{c}_{\mathbf{t}}$
$\displaystyle=L^{-1}_{\mathbf{t}}\left({\mathbf{t}}\odot\mathbf{a}_{\mathbf{t}}\right),\quad\text{and}$
$\displaystyle\mathbf{d}_{\mathbf{t}}$
$\displaystyle=\left(\frac{X^{{}^{\mkern-1.5mu\mathsf{T}}}X}{n}-\frac{\delta}{n}I\right)\left({\mathbf{t}}\odot{\mathbf{c}}_{\mathbf{t}}\right).$

##  4 The COMBSS Algorithm

In this section, we state the pseudocode for COMBSS and provide a discussion
on how to execute the algorithm.

Overall, Algorithm 1 consists of two stages. First stage, Step 1, calls a
gradient descent method of choice to obtain a solution $\mathbf{w}^{*}$ of the
unconstrained optimization problem (14). In the second stage, Steps 1 to 1, a
solution $\mathbf{t}^{*}$ of the box constrained optimization problem (10) is
obtained from $\mathbf{w}^{*}$ using the map
$\mathbf{t}(\mathbf{w})=1-\exp(-\mathbf{w}\odot\mathbf{w})$ and then mapping
small (respectively, large) values of $\mathbf{t}^{*}$ to $0$ (respectively,
$1$) (shown in Step 1), a solution of the discrete constrained best subset
selection problem (2) is obtained.

Input: Dataset $(X,\mathbf{y})$

Tuning parameters $\delta$ and $\lambda$

The initial vector $\mathbf{w}^{(0)}$

Threshold value $\tau\in[0,1]$

Termination condition TermCond

1

Output: A binary vector $\mathbf{s}$ as a solution to the optimization problem
(2)

2

3$\mathbf{w}\leftarrow\mathsf{GradientDescent}\left(\mathbf{w}^{(0)},\mathsf{TermCond},\nabla
g_{\lambda}\right)$

4$\mathbf{t}\leftarrow 1-\exp(-\mathbf{w}\odot\mathbf{w})$

5

6for  _$j=1,\dots,p$_ do

7 $s_{j}=I(t_{j}>\tau)$

8

return ${\mathbf{s}=(s_{1},\dots,s_{p})^{{}^{\mkern-1.5mu\mathsf{T}}}}$

Algorithm 1 COMBSS

The subroutine
$\mathsf{GradientDescent\left(\mathbf{w}^{(0)},\mathsf{TermCond},\nabla
g_{\lambda}\right)}$ calls a gradient descent method for minimizing the
objective function $g_{\lambda}(\mathbf{w})$, takes $\mathbf{w}^{(0)}$ as the
initial point for the gradient descent method, uses the gradient function
$\nabla g_{\lambda}$ for updating the vector $\mathbf{w}$ in each iteration of
the gradient descent method, and terminates to return the final value of
$\mathbf{w}$ when the given termination condition TermCond is satisfied. There
are several well-known gradient descent methods in the literature. For the
reader’s convenience, in Appendix B, we briefly review two such methods,
namely, the basic gradient descent method and the adaptive moment estimation
algorithm (simply known as Adam optimizer). The latter method is popular in
deep learning for solving highly non-convex optimization problems.

There are few things to keep in mind while selecting the initial point
$\mathbf{w}^{(0)}$ for the gradient descent. Observe from Theorem 4 that, for
any $j=1,\dots,p$,

  * •

$t_{j}(w_{j})=0$ if and only if $w_{j}=0$, and

  * •

${\frac{\partial g_{\lambda}(\mathbf{w})}{\partial w_{j}}=0}$ if $w_{j}=0$.

Thus, if we start a gradient descent algorithm with $w^{(0)}_{j}=0$ for some
$j$, both $w_{j}$ and $t_{j}$ can continue to take $0$ forever. As a result,
we might not learn the optimal value for $w_{j}$ (or, equivalently for
$t_{j}$). Thus, it is important to select $w^{(0)}_{j}$ away from $0$ for all
$j=1,\dots,p$. In our simulation, we select $\mathbf{w}_{0}$ such that
$\mathbf{t}(\mathbf{w}_{0})=(1/2,\dots,1/2)^{{}^{\mkern-1.5mu\mathsf{T}}}$,
the mid-point on the hyper cube $[0,1]^{p}$. This is intuitively a good choice
as it is at equal distance from all the corner points, of which one is the
(unknown) target solution of the best subset selection problem. However,
future work on how to better choose the initial vector $\mathbf{w}_{0}$ may
further improve the performance of COMBSS. In this paper, we do not further
investigate this aspect.

The termination condition, TermCond, for the gradient descent method is a user
choice depending on the computational budget. One trivial option is to
terminate the algorithm after reaching a predefined maximum number of
iterations or maximum allowed running time. Another option is to terminate
when the change in $\mathbf{w}$ (or, equivalently, in
$\mathbf{t}(\mathbf{w})$) is significantly small over a set of consecutive
iterations. In our simulations, we take the latter approach by terminating
when the max-norm

$\|\mathbf{t}(\mathbf{w}^{(l)})-\mathbf{t}(\mathbf{w}^{(l-1)})\|_{\infty}=\max_{j=1,\dots,p}|t_{j}(w^{(l)}_{j})-t_{j}(w^{(l-1)}_{j})|$

is smaller than a predefined error over a fixed number of consecutive
iterations, where $\mathbf{w}^{(l)}$ denotes the updated value of $\mathbf{w}$
in the $l$-th iteration of the gradient descent method. This particular
TermCond makes the algorithm robust by guaranteeing that further running of
the algorithm possibly can bring only an insignificant change in every
$t_{j}$.

Now, consider the last argument, $\nabla g_{\lambda}$, in the gradient descent
method. From Theorem 4, observe that computing the gradient $\nabla
g_{\lambda}(\mathbf{w})$ involves finding the values of the expression of the
form $L_{\mathbf{t}}^{-1}\mathbf{u}$ twice, first for computing
$\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}$ (using (7)) and then for
computing the vector $\mathbf{c}_{\mathbf{t}}$ (defined in Theorem 4). Since
$L_{\mathbf{t}}$ is a ${p\times p}$ matrix, computing the matrix inversion
$L_{\mathbf{t}}^{-1}$ can be problematic particularly in high-dimensional
cases where $p$ can be very large [11]. Since $L_{\mathbf{t}}^{-1}\mathbf{u}$
is the solution of the linear equation $L_{\mathbf{t}}\mathbf{z}=\mathbf{u}$,
in Section 6.1, we propose a conjugate gradient based approach for solving
this linear equation efficiently without the need to compute
$L_{\mathbf{t}}^{-1}$.

Moreover, again from Theorem 4, notice that the gradient $\nabla
g_{\lambda}(\mathbf{w})$ depends on both the tuning parameters $\delta$ and
$\lambda$. Particularly, $\delta$ is required for computing $L_{\mathbf{t}}$
and $\lambda$ is used in the penalty term $\lambda\|\mathbf{t}\|$ of the
objective function. In Section 5 we provide more detail on the roles of these
two parameters and instructions on how to choose them.

(a) Basic Gradient Descent

(b) Adam Optimizer

Figure 3: Convergence of $\mathbf{t}$ for a dataset during the execution of
the basic gradient descent (Panel (a)) and the Adam optimizer (Panel (b)). The
dataset is generated using Beta-type 1 with $k_{0}=5$, $n=100,p=1000$,
$\rho=0.8$, and $\mathsf{SNR}=5$; see Section 7 for more details on data
generation and selection of optimization parameters. For both the simulations,
$\lambda=0.1$ and $\delta=n$; see Section 5 for more discussion on how to
choose $\lambda$ and $\delta$. We terminated the gradient descent algorithms
when
$\|\mathbf{t}(\mathbf{w}^{(l)})-\mathbf{t}(\mathbf{w}^{(l-1)})\|_{\infty}\leq
0.001$ consecutively over $10$ iterations. In both the figures, solid lines
correspond to $\beta_{j}=0$ and remaining $5$ curves (with line style
$-\cdot-$) correspond to $\beta_{j}\neq 0$. Figure 4: Histogram of the $t_{j}$
values when COMBSS is executed on $500$ independent datasets using the Adam
optimizer. The parameters for data generation and simulation are the same as
given in Figure 3. In the above histogram, each bin is of size $0.01$ and the
height of the $i$-th bar is the fraction of datasets for which the final
output of the Adam optimizer has at least one $t_{j}$ converged to a point in
the $i$-th bin. Observe that since $p=1000$, except for a tiny fraction of the
datasets, each of the $1000$ $t_{j}$’s converges to a value close to either
$0$ or $1$. This indicates that the output of COMBSS is not sensitive to the
value of $\tau$ if it is chosen away from the terminal values $0$ and $1$.

Figure 3 illustrates the convergence of $t_{j}$ (within the gradient descent
algorithm) using an example dataset for both the basic gradient descent method
and the Adam optimizer. One useful observation from Figure 3 is that the Adam
optimizer takes a smaller number of iterations for convergence compared to the
basic gradient descent method. This is expected because the Adam optimizer
uses momentum updates which are not present in the basic gradient descent. In
our simulations in Section 7, we use only the Adam optimizer.

Finally, the threshold parameter $\tau\in(0,1)$ at Step 1 of Algorithm 1 is
useful for mapping all the insignificantly smaller $t_{j}$’s to $0$ and all
other $t_{j}$’s to $1$. Due to the tolerance allowed by the termination
condition TermCond, some $w_{j}$’s in the output of the gradient descent
algorithm can be almost zero but not exactly zero, even though they are meant
to converge to zero values. As a result, the corresponding $t_{j}$’s also take
values close to zero but not exactly zero because of the mapping from
$\mathbf{w}$ to $\mathbf{t}$ (see Step 1 on the algorithm). Therefore, $\tau$
helps in mapping these insignificantly small values to zeros. Figure 4
illustrates that the final output of the algorithm is not sensitive to $\tau$
as long as it is chosen slightly away from $0$. For all our simulations in
Section 7, we fix $\tau=0.5$.

##  5 Roles of Tuning Parameters $\delta$ and $\lambda$

In this section, we provide insights on how the tuning parameters $\delta$ and
$\lambda$ influence the objective function $f_{\lambda}(\mathbf{t})$ (or,
equivalently $g_{\lambda}(\mathbf{w})$) and hence the convergence of the
algorithm.

###  5.1 Controlling the Shape of $f_{\lambda}(\mathbf{t})$ through $\delta$

The normalized cost
$\|\mathbf{y}-X_{\mathbf{t}}\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}\|_{2}^{2}/n$
provides an estimator of the error variance. For any fixed $\mathbf{t}$, we
expect this variance (and hence $f_{\lambda}(\mathbf{t})$) to be almost the
same for all relatively large values of $n$, particularly, in situations where
the errors $\epsilon_{i}$’s are iid. Clearly, this is the case at all the
corner points $\mathbf{s}\in\\{0,1\\}^{p}$ where $L_{\mathbf{s}}$ is
invertible, because at these corner points, from Theorem 2,
$X_{\mathbf{s}}\widetilde{\boldsymbol{\beta}}_{\mathbf{s}}=X_{[\mathbf{s}]}\widehat{\boldsymbol{\beta}}_{[\mathbf{s}]}$,
which is independent of $\delta$. We would like to have a similar behavior at
all the interior points $\mathbf{t}\in(0,1)^{p}$ as well, so that for each
$\mathbf{t}$, the function $f_{\lambda}(\mathbf{t})$ is roughly the same for
all large values of $n$. Such consistent behavior is helpful in guaranteeing
that the convergence paths of the gradient descent method are approximately
the same for all large values of $n$.

(a) $n=100$ and $\delta=100$

(b) $n=100$ and $\delta=10000$

(c) $n=10000$ and $\delta=100$

(d) $n=10000$ and $\delta=10000$

Figure 5: Illustration of effects of $\delta$ on the objective function
$f_{\lambda}(\mathbf{t})$ with $\lambda=1.5$. A dataset consists of $10000$
samples generated from the illustrative model specified in Example 1. For
Figures 5(a) and 5(b), we used $100$ samples from this dataset.

First reconsider the model given in Example 1. Figure 5 shows surface plots of
$f_{\lambda}(\mathbf{t})$ for different values of $n$ and $\delta$. Here,
there are two plots, i.e. Figure 5(a) and Figure 5(d), that correspond to
$\delta=n$. In both these plots, we observe that the shape of the surface of
$f_{\lambda}(\mathbf{t})$ over $[0,1]^{p}$ is very similar.

To make this more explicit, we now show that the function
$f_{\lambda}(\mathbf{t})$, at any $\mathbf{t}$, takes almost the same value
for all large $n$ if we keep $\delta=c\,n$, for a fixed constant $c>0$, under
the assumption that the data samples are iid (this assumption simplifies the
following discussion; however, the conclusion holds more generally).

Observe that

$\frac{1}{n}\|\mathbf{y}-X_{\mathbf{t}}\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}\|_{2}^{2}=\frac{\mathbf{y}^{{}^{\mkern-1.5mu\mathsf{T}}}\mathbf{y}}{n}-2\,\boldsymbol{\gamma}_{\mathbf{t}}^{{}^{\mkern-1.5mu\mathsf{T}}}\left(\frac{X^{{}^{\mkern-1.5mu\mathsf{T}}}\mathbf{y}}{n}\right)+\boldsymbol{\gamma}_{\mathbf{t}}^{{}^{\mkern-1.5mu\mathsf{T}}}\left(\frac{X^{{}^{\mkern-1.5mu\mathsf{T}}}X}{n}\right)\boldsymbol{\gamma}_{\mathbf{t}},$

where

$\boldsymbol{\gamma}_{\mathbf{t}}=T_{\mathbf{t}}L_{\mathbf{t}}^{-1}T_{\mathbf{t}}\left(\frac{X^{{}^{\mkern-1.5mu\mathsf{T}}}\mathbf{y}}{n}\right).$

Since, under the iid assumption,
$\frac{\mathbf{y}^{{}^{\mkern-1.5mu\mathsf{T}}}\mathbf{y}}{n}$,
$\frac{X^{{}^{\mkern-1.5mu\mathsf{T}}}\mathbf{y}}{n}$, and
$\frac{X^{{}^{\mkern-1.5mu\mathsf{T}}}X}{n}$ converge element-wise as $n$
increases and $T_{\mathbf{t}}$ is independent of $n$, we would like to choose
$\delta$ such that $L_{\mathbf{t}}^{-1}$ also converges as $n$ increases. To
achieve this, first observe that

$L_{\mathbf{t}}=T_{\mathbf{t}}\left(\frac{X^{{}^{\mkern-1.5mu\mathsf{T}}}X}{n}\right)T_{\mathbf{t}}+\frac{\delta}{n}\left(I-T_{\mathbf{t}}^{2}\right).$

It is then evident that the choice $\delta=c\,n$ for a fixed constant ${c}$,
independent of $n$, makes $L_{\mathbf{t}}$ converging as $n$ increases, and
specifically, the choice $c=1$ justifies the behavior observed in Figure 5.

###  5.2 Controlling the Sparsity through $\lambda$

Intuitively, it is clear that larger the value of $\lambda$ the sparser the
solution offered by COMBSS, when all other parameters are fixed. We strengthen
this understanding, first visually and then also mathematically.

Visually, from Figure 2, we can understand that by selecting a suitable
positive $\lambda$ we can make COMBSS converge towards a sparse solution
$s^{*}$, which is in fact the optimal solution. Figure 6 shows that, for the
illustrative model in Example 1, if we take the value of $\lambda$ to be very
high, we can push the algorithm to converge to a much more sparse solution, in
this case, that point being $(0,0)^{{}^{\mkern-1.5mu\mathsf{T}}}$.

Figure 6: Illustration of effects of large $\lambda$ on the curvature of
$f_{\lambda}(\mathbf{t})$ and $g_{\lambda}(\mathbf{w})$. Here, the data is
generated from the model considered in Example 1 and $\lambda=100$.

Mathematically, from Theorem 4,

$\nabla
f_{\lambda}(\mathbf{t})=\boldsymbol{\zeta}_{\mathbf{t}}+\lambda\mathbf{1},\quad\mathbf{t}\in(0,1)^{p},$

and

$\nabla
g_{\lambda}(\mathbf{w})=\left(\boldsymbol{\zeta}_{\mathbf{t}}+\lambda\mathbf{1}\right)\odot\left(2\mathbf{w}\odot\exp(-\mathbf{w}\odot\mathbf{w})\right),\,\mathbf{w}\in\mathbb{R}^{p},$

where $\boldsymbol{\zeta}_{\mathbf{t}}$, given by (15), is independent of
$\lambda$. We note the following property of
$\boldsymbol{\zeta}_{\mathbf{t}}$.

######  Proposition 1.

For any $j=1,\dots,p$, if we fix $t_{i}$ for $i\neq j$,

$\lim_{t_{j}\downarrow 0}\boldsymbol{\zeta}_{\mathbf{t}}(j)=0.$

A proof of Proposition 1 is given in Appendix A. This result implies that for
any $j=1,\dots,p$,

$\lim_{t_{j}\downarrow 0}\tfrac{\partial f_{\lambda}(\mathbf{t})}{\partial
t_{j}}=\lambda,$

where $\lim_{t\downarrow 0}$ denotes the existence of the limit for any
sequence of $t$ that is converging to $0$ from right. Since
$\boldsymbol{\zeta}_{\mathbf{t}}$ is independent of $\lambda$, the above limit
implies that there is a window $(0,a_{j})$ such that the slope
$\tfrac{\partial f_{\lambda}(\mathbf{t})}{\partial t_{j}}>0$ for
$t_{j}\in(0,a_{j})$ and also the window size increases (i.e., $a_{j}$
increases) as $\lambda$ increases. As a result, for the function
$g_{\lambda}(\mathbf{w})$, there exists a constant $a^{\prime}_{j}$ such that

$\tfrac{\partial g_{\lambda}(\mathbf{w})}{\partial
w_{j}}\begin{cases}<0,\quad\text{for}\,\,a_{j}^{\prime}<t_{j}<0\\\
>0,\quad\text{for}\,\,0<t_{j}<a_{j}^{\prime}.\end{cases}$

In other words, for positive $\lambda$, there is a ‘valley’ on the surface of
$g_{\lambda}(\mathbf{w})$ along the line $w_{j}=0$ and the valley becomes
deeper and wider as $\lambda$ increases. In summary, the larger the values of
$\lambda$ the more $w_{j}$’s (or, equivalently $t_{j}$) are pushed towards $0$
by the algorithm and then a sparse model is selected (i.e, small number $k$ of
variables chosen). At the extreme a value
$\lambda_{\max}=\left\|y\right\|^{2}_{2}/n$ will force all the $t_{j}$ to be
equal to $0$ and the null model will be selected.

In practice, this tuning parameter $\lambda$ is generally calibrated by
defining a grid of $\lambda$ values and then exploits $K$-fold cross-
validation procedures (e.g., [25]) or model selection criterion (e.g., [1, 24,
5]).

##  6 Efficient Implementation of COMBSS

In this section, we focus on efficient implementation of COMBSS using the
conjugate gradient method, the Woodbury matrix identity, and the Banachiewicz
Inversion Formula.

Almost all the computational cost (or, the running time) of COMBSS is dictated
by the gradient descent method used at Step 1 of Algorithm 1, which in-turn
depends on how fast the gradient $\nabla g_{\lambda}(\mathbf{w})$ is computed
in each iteration of the gradient descent algorithm. For instance, in the
basic gradient descent algorithm, the point $\mathbf{w}^{(l)}$ at the $l$-th
iteration is computed by

$\mathbf{w}^{(l)}=\mathbf{w}^{(l-1)}-\alpha\nabla
g_{\lambda}(\mathbf{w}^{(l-1)}),$

using the point $\mathbf{w}^{(l-1)}$ from the previous iteration; see Appendix
B for more details on gradient descent methods.

First, in Section 6.1, we show how to estimate the gradients $\nabla
g_{\lambda}(\mathbf{w})$ using the conjugate gradient, an efficient method for
solving linear equations with positive-definite matrices. Then, in Section
6.2, we specifically focus on the high-dimensional case (i.e., $n<p$) and
propose an approach to reduce the dimension of the linear equation problem
solved by the conjugate gradient. In Section 6.3, we propose another dimension
reduction approach that ignores the columns of the design matrix $X$
correspond to $w_{j}=0$ for computing the gradients. Finally, in Section 6.4,
we put all these ideas together to make COMBSS fast.

###  6.1 Estimation of Gradients using Conjugate Gradient

Recall the expression of $L_{\mathbf{t}}$ from (6):

$\displaystyle L_{\mathbf{t}}$
$\displaystyle=\frac{1}{n}\left[X^{{{}^{\mkern-1.5mu\mathsf{T}}}}_{\mathbf{t}}X_{\mathbf{t}}+\delta\left(I-T_{\mathbf{t}}^{2}\right)\right].$

We have noticed earlier from Theorem 4 that the exact computation of $\nabla
g_{\lambda}(\mathbf{w})$ requires computing $L_{\mathbf{t}}^{-1}$, which can
be expensive. However, fortunately, all we really need for computing $\nabla
g_{\lambda}(\mathbf{w})$ is evaluating matrix-vector products of the form
$L_{\mathbf{t}}^{-1}\mathbf{u}$, which is the unique solution of the linear
equation $L_{\mathbf{t}}\mathbf{z}=\mathbf{u}$ for all
$\mathbf{t}\in[0,1)^{p}$ (refer to Theorem 1 for invertibility of
$L_{\mathbf{t}}$ for all $\mathbf{t}\in[0,1)^{p}$). Furthermore, the
performance of the gradient descent method does not change significantly as
long as a good approximation of the gradient $\nabla g_{\lambda}(\mathbf{w})$
is provided. As a result, it is optimal to seek efficient techniques for
solving $L_{\mathbf{t}}\mathbf{z}=\mathbf{u}$.

Solving linear equations efficiently is one of the important and well studied
problems in the field of linear algebra, and there are several elegant
approaches in the literature to address this problem. The Krylov subspace
methods, a family of iterative methods for solving linear equations, are most
popular in the literature for solving $A\mathbf{z}=\mathbf{u}$ efficiently
just by using the matrix-vector product $A\mathbf{u}$ (or,
$A^{{}^{\mkern-1.5mu\mathsf{T}}}\mathbf{u}$), see for instance, [11]. Among
these methods, the conjugate gradient method is well-known to work efficiently
for the problems where $A$ is symmetric positive-definite and can achieve a
close approximation to $A^{-1}\mathbf{u}$ within $O(p^{2})$ time. Since,
$L_{\mathbf{t}}$ is a symmetric positive-definite matrix, we can simply use
any standard software package for the conjugate gradient method in the
implementation of COMBSS.

###  6.2 Low-dimension vs High-dimension

As noticed in the previous subsection, the running time of the conjugate
gradient method for solving $A\mathbf{z}=\mathbf{u}$ depends on the dimension
of $A$. In particular, since $L_{\mathbf{t}}$ is of dimension $p\times p$, the
conjugate gradient method can return an approximation of
$L_{\mathbf{t}}^{-1}\mathbf{u}$ within $O(p^{2})$ time by fixing the maximum
number of iterations taken by the conjugate gradient method. This is true for
both low-dimensional models (where $p<n$) and high-dimensional models (where
$n<p$).

Below, we specifically focus on high-dimensional models and transform the
problem of solving the $p$-dimensional linear equation
$L_{\mathbf{t}}\mathbf{z}=\mathbf{u}$ to the problem of solving an
$n$-dimensional linear equation problem. This approach is based on a well-
known result in linear algebra called the Woodbury matrix identity. Since we
are calling the gradient descent method for solving a $n$-dimensional problem
(instead of $p$-dimensional), we can achieve a much lower overall
computational complexity for the high-dimensional models.

The following result is a consequence of the Woodbury matrix identity; refer
to Appendix A for a proof.

######  Theorem 5.

For $\mathbf{t}\in[0,1)^{p}$, let $S_{\mathbf{t}}$ be a $p$-dimensional
diagonal matrix with the $j$-th diagonal element being
$\nicefrac{{n}}{{\delta(1-t_{j}^{2})}}$ and

$\widetilde{L}_{\mathbf{t}}:=I+\frac{1}{n}X_{\mathbf{t}}S_{\mathbf{t}}X_{\mathbf{t}}^{{}^{\mkern-1.5mu\mathsf{T}}}.$

Then,

$\displaystyle L_{\mathbf{t}}^{-1}\mathbf{u}$
$\displaystyle=\left(S_{\mathbf{t}}\mathbf{u}\right)-\frac{1}{n}S_{\mathbf{t}}X_{\mathbf{t}}^{{{}^{\mkern-1.5mu\mathsf{T}}}}\widetilde{L}_{\mathbf{t}}^{-1}\left(X_{\mathbf{t}}S_{\mathbf{t}}\mathbf{u}\right).$

The above expression suggests that instead of solving the $p$-dimensional
problem $L_{\mathbf{t}}\mathbf{z}=\mathbf{u}$ directly, we can first solve the
$n$-dimensional problem
$\widetilde{L}_{\mathbf{t}}\mathbf{z}=\left(X_{\mathbf{t}}S_{\mathbf{t}}\mathbf{u}\right)$
and substitute the result in the above expression to get the value of
$L_{\mathbf{t}}^{-1}\mathbf{u}$. As a consequence, for high-dimensional
models, we can achieve smaller running times by this indirect approach for
computing the gradients in the execution of gradient descent algorithm.

###  6.3 A Dimension Reduction Approach

During the execution of the gradient descent algorithm, Step 1 of Algorithm 1,
some of the $w_{j}$’s (and hence $t_{j}$’s) can reach zero. Particularly, for
the basic gradient descent method, once $w_{j}$ reaches zero it continues to
be zero until the algorithm terminates. This is because, as mentioned earlier,
the update of $\mathbf{w}$ in the $l$-th iteration of the basic gradient
descent depends only on the gradient $g_{\lambda}(\mathbf{w}^{(l)})$, whose
$j$-th element

$\displaystyle\frac{\partial g_{\lambda}(\mathbf{w}^{(l)})}{\partial
w_{j}}=0\quad\text{if}\,\,w^{(l)}_{j}=0.$ (16)

Based on the fact that (16) holds, we ask: Can we just focus only on the terms
$\frac{\partial g_{\lambda}(\mathbf{w})}{\partial w_{j}}$ associated with
$w_{j}\neq 0$ in order to reduce the computation cost of computing the
gradient $\nabla g_{\lambda}(\mathbf{w})$? The answer is yes, and below we
provide an approach for achieving this.

To simplify the notation, let $\mathscr{P}=\\{1,\dots,p\\}$ and for any
$\mathbf{t}=(t_{1},\dots,t_{p})^{{}^{\mkern-1.5mu\mathsf{T}}}\in[0,1)^{p}$,
let $\mathscr{Z}_{\mathbf{t}}$ be the set of indices of the zero elements of
$\mathbf{t}$, that is,

$\displaystyle\mathscr{Z}_{\mathbf{t}}=\\{j:t_{j}=0,j=\mathscr{P}\\}.$ (17)

For a vector $\mathbf{u}\in\mathbb{R}^{p}$, we write
$(\mathbf{u})_{\scriptscriptstyle{+}}$ (respectively, $(\mathbf{u})_{0}$) to
denote the vector of dimension ${p-|\mathscr{Z}_{\mathbf{t}}|}$ (respectively,
$|\mathscr{Z}_{\mathbf{t}}|$) constructed from $\mathbf{u}$ by removing all
its elements with the indices in $\mathscr{Z}_{\mathbf{t}}$ (respectively, in
$\mathscr{P}\setminus\mathscr{Z}_{\mathbf{t}}$). Similarly, for a matrix $A$
of dimension $p\times p$, we write $(A)_{\scriptscriptstyle{+}}$
(respectively, $(A)_{0}$) to denote the new matrix constructed from $A$ by
removing its rows and columns with the indices in $\mathscr{Z}_{\mathbf{t}}$
(respectively, in $\mathscr{P}\setminus\mathscr{Z}_{\mathbf{t}}$). Then we
have the following result.

######  Theorem 6.

Suppose $\mathbf{t}\in[0,1)^{p}$. Then,

$\left(L_{\mathbf{t}}\right)_{\scriptscriptstyle{+}}=\frac{1}{n}\left[\left(T_{\mathbf{t}}\right)_{\scriptscriptstyle{+}}\left(X^{{}^{\mkern-1.5mu\mathsf{T}}}X\right)_{\scriptscriptstyle{+}}\left(T_{\mathbf{t}}\right)_{\scriptscriptstyle{+}}+\delta\left(I-\left(T_{\mathbf{t}}\right)_{\scriptscriptstyle{+}}\right)\right].$

Furthermore, we have

  * (i)

$\displaystyle\left(L_{\mathbf{t}}^{-1}\right)_{0}=\frac{n}{\delta}I\quad\text{and}\quad\left(L_{\mathbf{t}}^{-1}\right)_{\scriptscriptstyle{+}}=\left(\left(L_{\mathbf{t}}\right)_{\scriptscriptstyle{+}}\right)^{-1}$;

  * (ii)

$\displaystyle\left(\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}\right)_{0}=\mathbf{0}\quad\text{and}\quad\left(\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}\right)_{\scriptscriptstyle{+}}=\left(\left(L_{\mathbf{t}}\right)_{\scriptscriptstyle{+}}\right)^{-1}\left(({\mathbf{t}})_{\scriptscriptstyle{+}}\odot\left(\frac{X^{{}^{\mkern-1.5mu\mathsf{T}}}\mathbf{y}}{n}\right)_{\scriptscriptstyle{+}}\right)$;

  * (iii)

$\displaystyle\left(\mathbf{c}_{\mathbf{t}}\right)_{0}=\mathbf{0}\quad\text{and}\quad\left(\mathbf{c}_{\mathbf{t}}\right)_{\scriptscriptstyle{+}}=\left(\left(L_{\mathbf{t}}\right)_{\scriptscriptstyle{+}}\right)^{-1}\left(({\mathbf{t}})_{\scriptscriptstyle{+}}\odot\left(\mathbf{a}_{\mathbf{t}}\right)_{\scriptscriptstyle{+}}\right)$.

A proof of Theorem 6 is given in Appendix A. Theorem 6 (i) is based on the
Banachiewicz Inversion Formula (Lemma 1 in Appendix A). It shows that for
every $j\in\mathscr{Z}_{\mathbf{t}}$, all the off-diagonal elements of the
$j$-th row and the $j$-th column of $L^{-1}_{\mathbf{t}}$ are zero while its
$j$-th diagonal element is $n/\delta$, and all other elements of
$L_{\mathbf{t}}^{-1}$ (which constitute the sub-matrix
$\left(L_{\mathbf{t}}^{-1}\right)_{\scriptscriptstyle{+}}$) depend only on
$\left(L_{\mathbf{t}}\right)_{\scriptscriptstyle{+}}$, which can be computed
using only the columns of the design matrix $X$ with indices in
$\mathscr{P}\setminus\mathscr{Z}_{\mathbf{t}}$. As a consequence, (ii) and
(iii) imply that computing $\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}$ and
$\mathbf{c}_{\mathbf{t}}$ is equal to solving
$p_{\scriptscriptstyle{+}}$-dimensional linear equations of the form
$\left(L_{\mathbf{t}}^{-1}\right)_{\scriptscriptstyle{+}}\mathbf{z}=\mathbf{v}$,
where

$p_{\scriptscriptstyle{+}}=p-|\mathscr{Z}_{\mathbf{t}}|.$

Since $p_{\scriptscriptstyle{+}}\leq p$, solving such a
$p_{\scriptscriptstyle{+}}$-dimensional linear equation using the conjugate
gradient can be faster than solving the original $p$-dimensional linear
equation of the form $L_{\mathbf{t}}\mathbf{z}=\mathbf{u}$.

In summary, for a vector $\mathbf{t}\in[0,1)^{p}$ with some elements being
$0$, the values of $f_{\lambda}(\mathbf{t})$ and $\nabla
f_{\lambda}(\mathbf{t})$ do not depend on the columns $j$ of $X$ where
$t_{j}=0$. Therefore, we can reduce the computational complexity by removing
the all columns $j$ of the design matrix $X$ where $t_{j}=0$.

###  6.4 Making COMBSS Fast

In Section 6.3, we noted that when some elements $t_{j}$ of $\mathbf{t}$ are
zero, it is faster to compute the objective functions
($f_{\lambda}(\mathbf{t})$ and $g_{\lambda}(\mathbf{t})$) and their gradients
($\nabla f_{\lambda}(\mathbf{t})$ and $\nabla g_{\lambda}(\mathbf{t})$) by
ignoring the columns $j$ of the design matrix $X$. This is because the
effective dimension of the design matrix is $n\times
p_{\scriptscriptstyle{+}}$.

In Section 5.2, using Proposition 1, we have observed that for any $\lambda>0$
there is a ’valley’ on the surface of $g_{\lambda}(\mathbf{w})$ along
$w_{j}=0$ for all $j=1,\dots,p$. As a consequence, we have also observed that
for any $j$, when $w_{j}$ (or, equivalently, $t_{j}$) is sufficiently small
during the execution of the gradient descent method, it will eventually become
$0$, as it can not escape the valley. For instance, consider the example
illustrated in Figure 7. We observe that many of the $t_{j}$’s reached values
close to zero within the first $60$ iterations, and by the $80$-th iteration,
except one, the values of all $t_{j}$’s correspond to zero $\beta_{j}$’s in
the true model are zero or insignificantly small. This sparsity in
$\mathbf{t}$ is much more pronounced, when moving to $100$ iterations and then
to the final iteration.

(a) Initial state:
$\mathbf{t}=(1/2,1/2,\dots,1/2)^{{}^{\mkern-1.5mu\mathsf{T}}}$.

(b) After the $40$-th iteration.

(c) After the $60$-th iteration.

(d) After the $80$-th iteration.

(e) After $100$ iterations.

(f) After the $119$-th (final) iteration.

Figure 7: Illustration of convergence of $\mathbf{t}$ during the execution of
the basic gradient descent for an example dataset generated using Beta-type 1
with $k_{0}=5$, $p=20$, $n=100$, $\rho=0.5$, and $\mathsf{SNR}=3$; see Section
7 for more details on data generation. Here, we have taken $\lambda=0.2$ and
$\delta=n$. The set of indices of the (thick) blue bars correspond to the non-
zero elements in the true coefficient vector $\boldsymbol{\beta}$, and (thin)
red bars correspond to the zero elements in $\boldsymbol{\beta}$. The height
of the $j$-th bar denotes the value of $t_{j}$.

As a consequence, in the implementation of COMBSS, to reduce the computational
cost of estimating the gradients during the execution of the gradient descent
method, it is wise to map $w_{j}$ (and $t_{j}$) to $0$ when $w_{j}$ is almost
zero. We incorporate this truncation idea as follows. We first fix a small
constant $\eta$, say $\eta=0.001$. As we run the gradient descent algorithm,
when $t_{j}$ becomes smaller than $\eta$ for some $j\in\mathscr{P}$, we take
$t_{j}$ and $w_{j}$ to be zero and we stop updating them any further (that is,
$t_{j}$ and $w_{j}$ will continue to be zero until the gradient descent
algorithm terminates). In each iteration of the gradient descent algorithm,
the design matrix is updated by removing all the columns corresponding to zero
$t_{j}$’s. If the algorithm starts at $\mathbf{w}$ with all non-zero elements,
the effective dimension $p_{\scriptscriptstyle{+}}$, which denotes the number
of columns in the updated design matrix, monotonically decreases starting from
$p$. In an iteration, if $p_{\scriptscriptstyle{+}}>n$, we can use Theorem 5
to reduce the complexity of computing the gradients. On the other hand, when
$p_{\scriptscriptstyle{+}}$ falls below $n$, we directly use conjugate
gradient for computing the gradients without invoking Theorem 5.

The following list summarizes the key steps we take to make COMBSS fast:

  1. 1.

We use the well-known Adam optimizer to minimize the objective function
$g_{\lambda}$ for any given $\lambda$; refer to Figure 3 to see that the Adam
optimizer converges faster than the basic gradient descent.

  2. 2.

In every iteration of the gradient descent, we truncate all the $t_{j}$ (and
hence the corresponding $w_{j}$) to zero if they are smaller than the pre-
selected truncation parameter $\eta$. That is, for any $j\in\mathscr{P}$, we
assign $0$ to both $t_{j}$ and $w_{j}$ if $t_{j}<\eta$.

  3. 3.

In any iteration of the gradient descent, we invoke the Woodbury matrix
identity via Theorem 5 in computing $L_{\mathbf{t}}^{-1}\mathbf{u}$ if and
only if $p_{+}>n$.

  4. 4.

Whenever we need to compute an expression of the form $A^{-1}\mathbf{u}$, we
use the conjugate gradient method to solve the linear equation
$A\mathbf{z}=\mathbf{u}$ to obtain an approximation of $A^{-1}\mathbf{u}$.

Figure 8 illustrates that COMBSS is substantially much faster when the above
mentioned improvement ideas are incorporated in its implementation.

(a) $p=1000$

(b) $p=5000$

Figure 8: Running time of execution of COMBSS (until convergence with
respected to max-norm) for $\lambda=0.1$ for dataset generated using Beta-type
1 with $k_{0}=5$, $n=100$, $p=1000$, $\rho=0.8$, and $\mathsf{SNR}=5$ (same as
for Figures 3 and 4). The boxplot are based on 300 replications. Here we
compare running times for COMBSS using only conjugate gradient (CG), conjugate
gradient with Woodbury matrix identity (CG + WB), and conjugate gradient with
both Woodbury matrix identity and truncation improvement (CG +WB + Trunc). For
the truncation, $\eta=0.001$.

##  7 Simulation Experiments

In this section, through empirical research, we assess the efficiency and
correctness of our proposed method, COMBSS, and compare it with three main
existing methods, Forward Selection (FS), MIO, and Lasso, both in low-
dimensional and high-dimensional settings. Even though our method was
motivated by the challenge of features extraction for high-dimensional data
sets, it is important and instructive to investigate the behaviour of our
approach in a low dimensional setting as well.

###  7.1 Simulation Setting

The data is generated from the linear model:

$\mathbf{y}=X\boldsymbol{\beta}+\boldsymbol{\epsilon},\quad\text{ where
}\epsilon\sim N_{n}\left(0,\sigma^{2}\mathbb{I}_{n}\right),$ (18)

where each row of the predictor matrix $X$ is generated from a multivariate
normal distribution with zero mean and covariance matrix $\Sigma$ with
diagonal elements $\Sigma_{j,j}=1$ and off-diagonal elements
$\Sigma_{i,j}=\rho^{|i-j|}$, $i\neq j$, for some correlation parameter
$\rho\in(-1,1)$. In order to investigate a challenging situation, we use
$\rho=0.8$ to mimic strong correlation between predictors. For each
simulation, we fix the signal-to-noise ratio (SNR) and compute the variance
$\sigma^{2}$ of the noise $\boldsymbol{\epsilon}$ using

$\sigma^{2}=\frac{\boldsymbol{\beta}^{{}^{\mkern-1.5mu\mathsf{T}}}\Sigma\boldsymbol{\beta}}{\mathsf{SNR}}.$

We consider the following two simulation settings (also used in [4]):

  * •

Beta-type 1: $k_{0}=10$ components of $\boldsymbol{\beta}$ are equal to $1$ at
equally spaced indices between $1$ and $p$, and all other components of
$\boldsymbol{\beta}$ are $0$.

  * •

Beta-type 2: The first $k_{0}=10$ components of $\boldsymbol{\beta}$ are equal
to $1$ and all other components of $\boldsymbol{\beta}$ are equal to $0$.

Beta-type 1 will provide low correlation between active predictors while Beta-
type 2 assumes strong correlation between the active predictors. For both
these types, we investigate the performance of our method in low- and high-
dimensional settings. For the low-dimensional setting, we take $n=100$ and
$p=20$ for $\mathsf{SNR}\in\\{0.5,1,2,\dots,8\\}$, while for the high-
dimensional setting, $n=100$ and $p=1000$ for
$\mathsf{SNR}\in\\{2,3,\dots,8\\}$.

The comparison between the methods are given in terms of prediction error and
variable selection accuracy performances. The prediction error (PE)
performance of a method is defined as

$\text{
PE}=\left\|X\widehat{\boldsymbol{\beta}}-X\boldsymbol{\beta}\right\|_{2}^{2}/\left\|X\boldsymbol{\beta}\right\|_{2}^{2},$

where $\widehat{\boldsymbol{\beta}}$ is the estimated coefficient obtained by
the method. The variable selection accuracy performances used are sensibility
(true positive rate), specificity (true negative rate), accuracy, $F_{1}$
score, and the Mathew’s correlation coefficient (MCC) [6].

###  7.2 Model Tunning

In the low-dimensional setting, FS and MIO were tuned over $k=0,\ldots,20$.
For the high dimensional setting, FS was tuned over $k=0,\ldots,50$, while MIO
was tuned over $k=0,\ldots,30$ only due to time computational constraints. We
re-observed that MIO based on the Gurobi optimizer is quite time consuming for
high dimensional cases (see [14]). In this simulation we ran MIO through the R
package bestsubest offered in [13] and we fixed the default budget to 1 minute
40 seconds (per problem per $k$).

In low and high dimensional settings, the Lasso was tuned for 50 values of
$\lambda$ ranging from
$\lambda_{\max}=\left\|X^{T}\mathbf{y}\right\|_{\infty}$ to a small fraction
of $\lambda_{\max}$ on a log scale, as per the default in bestsubest package.
In both the low- and high-dimensional settings, COMBSS was also tuned using 50
values of $\lambda$ ranging from
$\lambda_{\max}=\left\|\mathbf{y}\right\|^{2}_{2}/n$ to a small fraction of
$\lambda_{\max}$. To be precise the 50 values of $\lambda$ are defined by
$\lambda_{l}=\lambda_{\max}\times 0.8^{(l-1)}$ with $l=1,\dots,50$. Further,
we have set $\eta=0.001$, $\tau=0.5$, and $\delta=n$ (involved in
$L_{\mathbf{t}}$). For all the simulations, we use Adam optimizer with
$\alpha=0.1$, $\xi_{1}=0.9$, and $\xi_{2}=0.999$ (see Appendix B.2).

For each subset size $k$ for FS and MIO and each tuning parameter $\lambda$
for Lasso and COMBSS, we evaluate the mean square error (MSE) using a
validation set of $5000$ samples from the true model defined in (18). Then,
for all the four methods, the best model is the one with the lowest MSE.

###  7.3 Simulation Results

Figure 9 presents the results in the low dimensional setting for both Beta-
type 1 and Beta-type 2. The panels in this figure display the average of MCC,
accuracy and prediction error over 500 replications for the four procedures,
and vertical bars denote one standard error. The other metrics (F1,
sensitivity and specificity) are presented in Figure 11 in Appendix C.

Figure 9: Performance results in terms of MCC, accuracy and prediction error,
for the low-dimensional case, where $n=100$, $p=20$, $\rho=0.8$, and true
active predictor is $10$ with Beta-type 1 (top plots) and Beta-type 2 (bottom
plots).

Overall, COMBSS outperforms FS and MIO procedures in terms of MCC, accuracy,
and prediction error. It also outperforms the Lasso in terms of MCC and
accuracy even for Beta-type 2 (high correlation of the true active signal) and
for small SNR, both the approaches give similar results. We can see that FS
and MIO have similar pattern for Beta-type 1 while FS seems to be preferable
for Beta-type 2, where there the true set of variables are highly correlated.
Note that the Lasso presents lower prediction error in general as it tends to
provide dense model compared to other methods. As a result, the Lasso suffers
with low specificity which is clearly exhibited in Figure 11 in Appendix C.

Figure 10: Performance results in terms of MCC and prediction error, for the
high-dimensional cases: $n=100$, $p=1000$, $\rho=0.8$, and true active
predictor is $10$ with Beta-type 1 (top plots) and Beta-type 2 (bottom plots).

Figure 10 presents the results in the high-dimensional setting for both Beta-
type 1 and Beta-type 2. The panels in this figure display average MCC and
prediction error over $10$ replications for all the four methods, and vertical
bars denote one standard error. The other metrics (F1, sensitivity, and
specificity) are presented in Figure 12 in Appendix C. We do not present
accuracy for the high-dimensional setting, because even a procedure which
always selects the null model will get an accuracy of $990/1000=0.99$). We
restrict this simulation to $10$ replications due to the time constraint of
running MIO in this setting. Even if the other approaches have been ran other
$500$ replications easily, we present here only the results on the same $10$
repetitions.

For the Beta-type 1 scenario, the three best subset selection methods (FS, MIO
and COMBS) give similar results and outperform the Lasso in terms of MCC. In
this situation, the Lasso obtains lower prediction error for lower SNR but
still suffers from selecting dense models and thus gets small specificity (see
12 in Appendix C). Overall, COMBSS outperforms the other procedures in the
Beta-type 2 scenario.

##  8 Discussion

In this paper, we have introduced COMBSS, a novel continuous optimization
method towards best subset selection in linear regression. We developed this
method with a particular focus on high-dimensional settings where the number
of features $p$ can be very large. Key goal of COMBSS is to extend the highly
complex discrete constrained best subset selection problem to an unconstrained
continuous optimization problem. In particular, COMBSS involves extending the
objective function of the best subset selection, which is defined at the
corners of the hypercube $[0,1]^{p}$, to a differentiable function defined on
the whole hypercube. For this extended function, starting from an interior
point, a gradient descent method is executed to find a corner of the hypercube
where the objective function is minimum.

The major technological innovations were made in the coding of COMBSS and we
have made Python and R codes available at https://github.com/saratmoka/COMBSS-
Python-VIGNETTE and https://github.com/benoit-liquet/COMBSS-R-VIGNETTE,
respectively. In this article, we have specifically highlighted some of the
steps that make COMBSS computationally fast, including using conjugate
gradient and other linear algebra techniques. In our implementation of COMBSS,
the primary operations involved are matrix-vector product, vector-vector
element-wise product, and scalar-vector product (noting that the primary
operation for executing a conjugate gradient method for solving a linear
equation $A\mathbf{z}=\mathbf{u}$ is the matrix-vector product $A\mathbf{u}$).
All these operations are known to execute faster on GPU based computers using
parallel programming. Such parallel programming based implementation of COMBSS
could substantially increase its speed.

In the paper, we have conducted several simulation experiments in both low-
dimensional and high-dimensional setups to illustrate the good performance of
COMBSS in comparison to the existing popular methods: Forward Stepwise, Lasso,
and Mixed Integer Optimization. These promising empirical results highlight
the potential of COMBSS for feature extractions. We leave extensive
performance comparisons with Lasso variants such as relaxed Lasso [20] and
best subset selection based on coordinate descent and local combinatorial
optimization [15] to future work.

We highlight that the primary focus of this paper was to introduce the novel
COMBSS approach for feature selection and illustrate its performance. In
ongoing work, we are incorporating additional tricks and functionality to make
COMBSS even faster and more user-friendly. We expect, similarly to the
significant body of work that focuses on the Lasso and on MIO, respectively,
that there are many avenues that can be explored and investigated for building
on the presented COMBSS framework. Particularly, to tackle best subset
selection when problems are ultra high dimensional.

The following is a partial list of considerations and future research, that
goes beyond the scope of this article.

  * •

Despite using the same grid of $\lambda$ values for COMBSS in all our
simulations, it is competitive with all the existing methods considered in
this paper. Future work on more efficient selection of the grid may further
enhance the performance of COMBSS. As with all feature selection methods where
the tuning of parameters is part of the challenge, future research is needed
to appropriately investigate to what extent insights into the tuning of the
Lasso and other regularisation methods translates to the tuning of COMBSS.

  * •

Although, in this paper, COMBSS achieves sparsity via $\sum_{j=1}^{p}t_{j}$,
which turns out to be the $\mathcal{L}_{1}$-norm of $\mathbf{t}$ because
$\|\mathbf{t}\|_{1}=\sum_{j=1}^{p}t_{j}$ for $\mathbf{t}\in[0,1]^{p}$.
However, regularising with the $\mathcal{L}_{1}$-norm is not crucial and in
fact sparsity in the solution of COMBSS can be achieved using other norms; for
example we could replace $\sum_{j=1}^{p}t_{j}$ with
$\|\mathbf{t}\|_{r}^{r}=\sum_{j=1}^{p}t_{j}^{r}$ for $r\geq 1$. Future work
could focus on the selection of an optimal metric for the penalty term in
COMBSS.

  * •

In this paper, we explore first-order continuous optimization methods for
executing COMBSS. The learning rate parameter $\alpha$ of these methods can
influence their performance. Since, it is possible to compute the second-order
derivatives of the objective functions $f_{\lambda}$ and $g_{\lambda}$, future
work can improve the performance of COMBSS by exploiting second-order methods
to avoid the calibration of the learning rate.

  * •

Recall that the purpose of the mapping $t_{j}(w_{j})=1-\exp(-w_{j}^{2})$ is to
convert the boxed constrained optimization problem in (10) to the
unconstrained optimization problem (14). We note that any strictly increasing
function on $\mathbb{R}_{+}$ that starts with $0$ at $0$ and goes to $1$ as
$w_{j}\to\infty$ will also convert the boxed constrained optimization problem
to the unconstrained optimization problem. For instance, we could consider
$t_{j}(w_{j})=\tanh(w_{j}^{2})$. Finding an optimal such function could be a
part of future work.

In conclusion, to the best of our knowledge, COMBSS is the first method that
exploits unconstrained continuous optimization to find sets of models that are
candidates for the best subset of variables in linear regression. As
highlighted in this discussion and throughout the article, COMBSS potentially
opens a new avenue of research in high-dimensional feature selection, going
well beyond the linear regression model.

## Acknowledgement

Muller was supported by the Australian Research Council (DP210100521).

##  Appendix A Proofs

###### Proof of Theorem 1.

Since both $X_{\mathbf{t}}^{{}^{\mkern-1.5mu\mathsf{T}}}X_{\mathbf{t}}$ and
$T_{\mathbf{t}}$ are symmetric, the symmetry of $L_{\mathbf{t}}$ is obvious.
We now show that $L_{\mathbf{t}}$ is positive-definite for
$\mathbf{t}\in[0,1)^{p}$ by establishing

$\displaystyle\mathbf{u}^{{}^{\mkern-1.5mu\mathsf{T}}}L_{\mathbf{t}}\mathbf{u}>0,\quad\text{for
all}\,\,\mathbf{u}\in\mathbb{R}^{p}\setminus\\{\boldsymbol{0}\\}.$ (19)

Trivially, $X_{\mathbf{t}}^{{}^{\mkern-1.5mu\mathsf{T}}}X_{\mathbf{t}}$ is a
positive semi-definite matrix, because

$\mathbf{u}^{{}^{\mkern-1.5mu\mathsf{T}}}X_{\mathbf{t}}^{{}^{\mkern-1.5mu\mathsf{T}}}X_{\mathbf{t}}\mathbf{u}=\|X_{\mathbf{t}}\mathbf{u}\|_{2}^{2}\geq
0.$

In addition, for all $\mathbf{t}\in[0,1)^{p}$,
$\delta\left(I-T_{\mathbf{t}}^{2}\right)$ is a positive-definite matrix due to
the fact that $\delta>0$ and

$\displaystyle\mathbf{u}^{{}^{\mkern-1.5mu\mathsf{T}}}\left(I-T_{\mathbf{t}}^{2}\right)\mathbf{u}$
$\displaystyle=\|\mathbf{u}\|_{2}^{2}-\|T_{\mathbf{t}}\mathbf{u}\|_{2}^{2}$
$\displaystyle=\sum_{j=1}^{p}u_{j+1}^{2}(1-t_{j}^{2}),$ (20)

which is strictly positive if $\mathbf{t}\in[0,1)^{p}$ and
$\mathbf{u}\in\mathbb{R}^{p}\setminus\\{\boldsymbol{0}\\}$.

Since positive-definite matrices are invertible, we have
$L_{\mathbf{t}}^{+}=L_{\mathbf{t}}^{-1}$, and thus,
$\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}=L_{\mathbf{t}}^{-1}\left(\frac{X_{\mathbf{t}}^{{}^{\mkern-1.5mu\mathsf{T}}}\mathbf{y}}{n}\right)$.
∎

We use Lemma 1 for proving Theorem 2. For a proof of this lemma, we refer to
[27].

######  Lemma 1 (Banachiewicz Inversion Formula).

For a block matrix

$B=\begin{pmatrix}A_{11}&A_{12}\\\ A_{21}&A_{22}\end{pmatrix},$

we have:

  * (i)

If $A_{11}$ is nonsingular, then $B$ is nonsingular if and only if the Schur
compliment ${S=A_{22}-A_{21}A_{11}^{-1}A_{12}}$ is nonsingular.

  * (ii)

Furthermore, if both $A_{11}$ and $A_{22}$ are nonsingular, then

$\displaystyle B^{-1}$
$\displaystyle=\begin{pmatrix}I&-A^{-1}_{11}A_{12}S^{-1}\\\
\mathbf{0}&S^{-1}\end{pmatrix}\begin{pmatrix}A_{11}^{-1}&\mathbf{0}\\\
-A_{21}A_{11}^{-1}&I\end{pmatrix}.$

###### Proof of Theorem 2.

It is well-known that the inverse of a matrix after a permutation of rows
(respectively, columns) is identical to the matrix obtained by applying the
same permutation on columns (respectively, rows) on the inverse of the matrix.
Therefore without loss of generality we assume that all the zero-elements of
$\mathbf{s}\in\\{0,1\\}^{p}$ appear at the end, in the form:

$\mathbf{s}=(s_{1},\dots,s_{m},0,\dots,0),$

where $m$ indicates the number of non-zeros in $\mathbf{s}$. Recall that
${X}_{[\mathbf{s}]}$ is the matrix of size $n\times|\mathbf{s}|$ created by
keeping only columns $j$ of $X$ for which $s_{j}=1$. Then, with the notation
$(A)_{\scriptscriptstyle{+}}$, defined in Section 6.3, $L_{\mathbf{s}}$ is
given by,

$\displaystyle L_{\mathbf{s}}$
$\displaystyle=\frac{1}{n}\left[\begin{pmatrix}\left(X_{\mathbf{s}}^{{}^{\mkern-1.5mu\mathsf{T}}}X_{\mathbf{s}}\right)_{\scriptscriptstyle{+}}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}\end{pmatrix}+\delta\begin{pmatrix}I-\left(T^{2}_{\mathbf{s}}\right)_{\scriptscriptstyle{+}}&\mathbf{0}\\\
\mathbf{0}&I\end{pmatrix}\right]$
$\displaystyle=\frac{1}{n}\begin{pmatrix}X_{[\mathbf{s}]}^{{}^{\mkern-1.5mu\mathsf{T}}}X_{[\mathbf{s}]}&\mathbf{0}\\\
\mathbf{0}&\delta I\end{pmatrix},$ (21)

where we used the fact that for binary vectors $\mathbf{s}$,

$\left(X_{\mathbf{s}}^{{}^{\mkern-1.5mu\mathsf{T}}}X_{\mathbf{s}}\right)_{\scriptscriptstyle{+}}=X_{[\mathbf{s}]}^{{}^{\mkern-1.5mu\mathsf{T}}}X_{[\mathbf{s}]}\quad\text{and}\quad\left(T^{2}_{\mathbf{s}}\right)_{\scriptscriptstyle{+}}=I.$

From Lemma 1, it is evident that $L_{\mathbf{t}}$ is invertible if and only if
$X_{[\mathbf{s}]}^{{}^{\mkern-1.5mu\mathsf{T}}}X_{[\mathbf{s}]}$ is
invertible. Since $\widehat{\boldsymbol{\beta}}_{[\mathbf{s}]}$ exists if and
only if $X_{[\mathbf{s}]}^{{}^{\mkern-1.5mu\mathsf{T}}}X_{[\mathbf{s}]}$ is
invertible, $L_{\mathbf{s}}$ is invertible if and only if
$\widehat{\boldsymbol{\beta}}_{[\mathbf{s}]}$ exists. Thus, using Lemma 1, we
have

$\displaystyle L_{\mathbf{s}}^{-1}$
$\displaystyle=n\begin{pmatrix}I&\mathbf{0}\\\
\mathbf{0}&\frac{1}{\delta}I\end{pmatrix}\begin{pmatrix}\left({X}_{[\mathbf{s}]}^{{}^{\mkern-1.5mu\mathsf{T}}}{X}_{[\mathbf{s}]}\right)^{-1}&\mathbf{0}\\\
0&I\end{pmatrix}$
$\displaystyle=n\begin{pmatrix}\left({X}_{[\mathbf{s}]}^{{}^{\mkern-1.5mu\mathsf{T}}}{X}_{[\mathbf{s}]}\right)^{-1}&\mathbf{0}\\\
0&\frac{1}{\delta}I\end{pmatrix}$ (22)

To show that
$X_{\mathbf{s}}\widetilde{\boldsymbol{\beta}}_{\mathbf{s}}={X}_{[\mathbf{s}]}\widehat{\boldsymbol{\beta}}_{[\mathbf{s}]}$
for $\mathbf{s}\in\\{0,1\\}^{p}$ where
${X}_{[\mathbf{s}]}^{{}^{\mkern-1.5mu\mathsf{T}}}{X}_{[\mathbf{s}]}$ is
invertible, using (22), we obtain

$L_{\mathbf{s}}^{-1}T_{\mathbf{s}}X^{{}^{\mkern-1.5mu\mathsf{T}}}\mathbf{y}=\left({X}_{[\mathbf{s}]}^{{}^{\mkern-1.5mu\mathsf{T}}}{X}_{[\mathbf{s}]}\right)^{-1}{X}_{[\mathbf{s}]}^{{}^{\mkern-1.5mu\mathsf{T}}}\mathbf{y},$

which guarantees
${X}_{[\mathbf{s}]}\widehat{\boldsymbol{\beta}}_{[\mathbf{s}]}=X_{\mathbf{s}}\widetilde{\boldsymbol{\beta}}_{\mathbf{s}}$.
∎

We use Lemma 2 for proving Theorem 3. Here, recall that for any vector
$\mathbf{u}\in\mathbb{R}^{p}$, the max-norm of $\mathbf{u}$ is given by

$\|\boldsymbol{u}\|_{\infty}=\max_{j=1,\dots,p}|u_{j}|.$

######  Lemma 2.

([3, Theorem 16]) Suppose $A\in\mathbb{R}^{p\times p}$ is a non-trivial real-
valued square matrix with eigenvalues
$\boldsymbol{\eta}=(\eta_{1},\dots,\eta_{p})^{{}^{\mkern-1.5mu\mathsf{T}}}$
and $r$ is a real number such that

$0<r<\frac{2}{\displaystyle\|\boldsymbol{\eta}\|_{\infty}^{2}}.$

Then,

$A^{+}=r\sum_{k=0}^{\infty}\left(I-rA^{{}^{\mkern-1.5mu\mathsf{T}}}A\right)^{k}A^{{}^{\mkern-1.5mu\mathsf{T}}}.$

###### Proof of Theorem 3.

Consider a point $\mathbf{t}\in[0,1]^{p}$ and a sequence
$\mathbf{t}^{(1)},\mathbf{t}^{(2)},\dots\in[0,1)^{p}$ converging to
$\mathbf{t}$. Then $L_{\mathbf{t}^{(l)}}$ converges to $L_{\mathbf{t}}$
element-wise as $i\to\infty$. To see this, observe from the definitions of
$L_{\mathbf{t}^{(l)}}$ and $T_{\mathbf{t}^{(l)}}$ that, with
$Z=\frac{X^{{}^{\mkern-1.5mu\mathsf{T}}}X}{n}-\frac{\delta}{n}I$,

$\displaystyle L_{\mathbf{t}^{(l)}}$
$\displaystyle=T_{\mathbf{t}^{(l)}}\,Z\,T_{\mathbf{t}^{(l)}}+\frac{\delta}{n}I=\left(T_{\mathbf{t}}-T_{\Delta\mathbf{t}^{(l)}}\right)\,Z\,\left(T_{\mathbf{t}}-T_{\Delta\mathbf{t}^{(l)}}\right)+\frac{\delta}{n}I,$

where $\Delta\mathbf{t}^{(l)}={\mathbf{t}}-\mathbf{t}^{(l)}$. Expanding the
above expression we obtain

$\displaystyle L_{\mathbf{t}^{(l)}}=L_{\mathbf{t}}-D_{i},$

where the symmetric matrix

$D_{i}=T_{\mathbf{t}}\,Z\,T_{\Delta\mathbf{t}^{(l)}}+T_{\Delta\mathbf{t}^{(l)}}\,Z\,T_{\mathbf{t}}-T_{\Delta\mathbf{t}^{(l)}}\,Z\,T_{\Delta\mathbf{t}^{(l)}}.$

Clearly, element-wise, $D_{i}$ is going to $\mathbf{0}$ as $i\to\infty$. In
other words, $L_{\mathbf{t}^{(l)}}$ converges to $L_{\mathbf{t}}$ element-wise
as $i\to\infty$.

Now suppose that $\boldsymbol{\eta}_{l}$ and $\boldsymbol{\eta}$ are the
vector of the eigenvalues of $L_{\mathbf{t}^{(l)}}$ and $L_{\mathbf{t}}$,
respectively. Then Proposition 5.2.2 (c) of [2] guarantees the convergence of
the eigenvalues of $L_{\mathbf{t}^{(l)}}$ to that $L_{\mathbf{t}}$. Thus,

$\lim_{l\to\infty}\|\boldsymbol{\eta}_{l}\|_{\infty}=\|\boldsymbol{\eta}\|_{\infty}.$

This implies that there exists a large number $N$ and a constant $r$ such that

$0<r<\frac{2}{\displaystyle\|\boldsymbol{\eta}\|_{\infty}^{2}}$

and

$0<r<\frac{2}{\displaystyle\|\boldsymbol{\eta}_{l}\|_{\infty}^{2}},\quad\text{for
all}\,\,l\geq N.$

For such $r$ and $N$, from Lemma 2,

$L_{\mathbf{t}}^{+}=r\sum_{k=0}^{\infty}\left(I-rL_{\mathbf{t}}^{2}\right)^{k}L_{\mathbf{t}},$

and

$L_{\mathbf{t}^{(l)}}^{+}=L_{\mathbf{t}^{(l)}}^{-1}=r\sum_{k=0}^{\infty}\left(I-rL_{\mathbf{t}^{(l)}}^{2}\right)^{k}L_{\mathbf{t}^{(l)}},\quad\text{for
all}\,\,l\geq N,$

where the equality $L_{\mathbf{t}^{(l)}}^{+}=L_{\mathbf{t}^{(l)}}^{-1}$
follows from Theorem 1 since each $\mathbf{t}^{(l)}\in[0,1)^{p}$. Now observe
that

$\displaystyle
L_{\mathbf{t}^{(l)}}^{2}=\left(L_{\mathbf{t}}-D_{l}\right)^{2}=L_{\mathbf{t}}^{2}+G_{l},$

where $G_{l}=D_{l}^{2}-D_{l}L_{\mathbf{t}}-L_{\mathbf{t}}D_{l}$, which
converges to $\mathbf{0}$ element-wise as $l\to\infty$. Since

$L_{\mathbf{t}^{(l)}}^{-1}=r\sum_{k=0}^{\infty}\left(I-rL_{\mathbf{t}}^{2}-rG_{l}\right)^{k}\left(L_{\mathbf{t}}-D_{l}\right),\quad\text{for
all}\,\,l\geq N,$

$L_{\mathbf{t}^{(l)}}^{-1}$ converges to $L_{\mathbf{t}}^{+}$ element-wise.

As a consequence, since $X_{\mathbf{t}^{(l)}}$ converges element-wise to
$X_{\mathbf{t}}$, we can guarantee that

$\displaystyle\lim_{l\to\infty}\widetilde{\boldsymbol{\beta}}_{\mathbf{t}^{(l)}}$
$\displaystyle=\lim_{l\to\infty}L_{\mathbf{t}^{(l)}}^{-1}\left(\frac{X_{\mathbf{t}^{(l)}\mathbf{y}}}{n}\right)$
$\displaystyle=\widetilde{\boldsymbol{\beta}}_{\mathbf{t}},$

and thus we trivially have from the definition of $f_{\lambda}$ that
$\lim_{l\to\infty}f_{\lambda}(\mathbf{t}^{(l)})=f_{\lambda}(\mathbf{t}).$ ∎

We use Lemma 3, which obtains the partial derivatives of
$\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}$ with respect to the elements of
$\mathbf{t}$, for proving Theorem 4.

######  Lemma 3.

For any $\mathbf{t}\in(0,1)^{p}$, the derivatives of
$\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}$ are given by

$\frac{\partial\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}}{\partial
t_{j}}=L_{\mathbf{t}}^{-1}\left[E_{j}-E_{j}ZT_{\mathbf{t}}L^{-1}_{\mathbf{t}}T_{\mathbf{t}}-T_{\mathbf{t}}ZE_{j+1}L^{-1}_{\mathbf{t}}T_{\mathbf{t}}\right]\left(\frac{X^{{}^{\mkern-1.5mu\mathsf{T}}}\mathbf{y}}{n}\right),\quad
j=1,\dots,p,$

where

$Z=\left(\frac{X^{{}^{\mkern-1.5mu\mathsf{T}}}X}{n}-\frac{\delta}{n}I\right),$

and $E_{j}$ is a square matrix of dimension $p\times p$ with $1$ at the
$(j,j)$-th position and $0$ everywhere else.

###### Proof.

Existence of $\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}$ for every
$\mathbf{t}\in(0,1)^{p}$ and $\delta>0$ follows from Theorem 1 which states
that $L_{\mathbf{t}}$ is positive-definite which guarantees the invertibility
of $L_{\mathbf{t}}$.

Since
$\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}=L_{\mathbf{t}}^{-1}T_{\mathbf{t}}\left(\frac{X^{{}^{\mkern-1.5mu\mathsf{T}}}\mathbf{y}}{n}\right)$,
using matrix calculus, for any $j=1,\dots,p$,

$\displaystyle\frac{\partial\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}}{\partial
t_{j}}$
$\displaystyle=\frac{\partial\left(L^{-1}_{\mathbf{t}}T_{\mathbf{t}}\right)}{\partial
t_{j}}\left(\frac{X^{{}^{\mkern-1.5mu\mathsf{T}}}\mathbf{y}}{n}\right)$
$\displaystyle=\left[\frac{\partial L^{-1}_{\mathbf{t}}}{\partial
t_{j}}T_{\mathbf{t}}+L^{-1}_{\mathbf{t}}\frac{\partial
T_{\mathbf{t}}}{\partial
t_{j}}\right]\left(\frac{X^{{}^{\mkern-1.5mu\mathsf{T}}}\mathbf{y}}{n}\right)$
$\displaystyle=\left[L^{-1}_{\mathbf{t}}\frac{\partial
T_{\mathbf{t}}}{\partial t_{j}}-L^{-1}_{\mathbf{t}}\frac{\partial
L_{\mathbf{t}}}{\partial
t_{j}}L^{-1}_{\mathbf{t}}T_{\mathbf{t}}\right]\left(\frac{X^{{}^{\mkern-1.5mu\mathsf{T}}}\mathbf{y}}{n}\right),$

where we used differentiation of an invertible matrix which implies

$\frac{\partial L^{-1}_{\mathbf{t}}}{\partial
t_{j}}=-L^{-1}_{\mathbf{t}}\frac{\partial L_{\mathbf{t}}}{\partial
t_{j}}L^{-1}_{\mathbf{t}}.$

Since

$\frac{\partial T_{\mathbf{t}}}{\partial t_{j}}=E_{j},$

and the fact that
$L_{\mathbf{t}}=T_{\mathbf{t}}ZT_{\mathbf{t}}+\frac{\delta}{n}I$, we get

$\displaystyle\frac{\partial L_{\mathbf{t}}}{\partial t_{j}}$
$\displaystyle=\frac{\partial T_{\mathbf{t}}}{\partial
t_{i}}ZT_{\mathbf{t}}+T_{\mathbf{t}}Z\frac{\partial T_{\mathbf{t}}}{\partial
t_{j}}$ $\displaystyle=E_{j}ZT_{\mathbf{t}}+T_{\mathbf{t}}ZE_{j}.$

Therefore,

$\displaystyle L^{-1}_{\mathbf{t}}\frac{\partial T_{\mathbf{t}}}{\partial
t_{j}}-L^{-1}_{\mathbf{t}}\frac{\partial L_{\mathbf{t}}}{\partial
t_{j}}L^{-1}_{\mathbf{t}}T_{\mathbf{t}}$
$\displaystyle=L^{-1}_{\mathbf{t}}E_{j}-L^{-1}_{\mathbf{t}}E_{j}ZT_{\mathbf{t}}L^{-1}_{\mathbf{t}}T_{\mathbf{t}}-L^{-1}_{\mathbf{t}}T_{\mathbf{t}}ZE_{j}L^{-1}_{\mathbf{t}}T_{\mathbf{t}}$
$\displaystyle=L_{\mathbf{t}}^{-1}\left[E_{j}-E_{j}ZT_{\mathbf{t}}L^{-1}_{\mathbf{t}}T_{\mathbf{t}}-T_{\mathbf{t}}ZE_{j}L^{-1}_{\mathbf{t}}T_{\mathbf{t}}\right].$

This completes the proof Lemma 3. ∎

###### Proof of Theorem 4.

To obtain the gradient $\nabla f_{\lambda}(\mathbf{t})$ for
$\mathbf{t}\in(0,1)^{p}$, let
$\boldsymbol{\gamma}_{\mathbf{t}}=T_{\mathbf{t}}\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}=\mathbf{t}\odot\widetilde{\boldsymbol{\beta}}$.
Then,

$\displaystyle\|\mathbf{y}-X_{\mathbf{t}}\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}\|_{2}^{2}$
$\displaystyle=\|\mathbf{y}-X\boldsymbol{\gamma}_{\mathbf{t}}\|_{2}^{2}$
$\displaystyle=\left(\mathbf{y}-X\boldsymbol{\gamma}_{\mathbf{t}}\right)^{{}^{\mkern-1.5mu\mathsf{T}}}\left(\mathbf{y}-X\boldsymbol{\gamma}_{\mathbf{t}}\right)$
$\displaystyle=\mathbf{y}^{{}^{\mkern-1.5mu\mathsf{T}}}\mathbf{y}-2\,\boldsymbol{\gamma}_{\mathbf{t}}^{{}^{\mkern-1.5mu\mathsf{T}}}\left(X^{{}^{\mkern-1.5mu\mathsf{T}}}\mathbf{y}\right)+\boldsymbol{\gamma}_{\mathbf{t}}^{{}^{\mkern-1.5mu\mathsf{T}}}\left(X^{{}^{\mkern-1.5mu\mathsf{T}}}X\right)\boldsymbol{\gamma}_{\mathbf{t}}.$
(23)

Consequently,

$\displaystyle\frac{\partial f_{\lambda}(\mathbf{t})}{\partial t_{j}}$
$\displaystyle=\frac{1}{n}\frac{\partial}{\partial
t_{j}}\left[\|\mathbf{y}-X_{\mathbf{t}}\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}\|_{2}^{2}\right]+\lambda$
$\displaystyle=-\frac{2}{n}\left(\frac{\partial\boldsymbol{\gamma}_{\mathbf{t}}}{\partial
t_{j}}\right)^{{}^{\mkern-1.5mu\mathsf{T}}}\left(X^{{}^{\mkern-1.5mu\mathsf{T}}}\mathbf{y}\right)+\frac{2}{n}\left(\frac{\partial\boldsymbol{\gamma}_{\mathbf{t}}}{\partial
t_{j}}\right)^{{}^{\mkern-1.5mu\mathsf{T}}}\left(X^{{}^{\mkern-1.5mu\mathsf{T}}}X\right)\boldsymbol{\gamma}_{\mathbf{t}}+\lambda$
$\displaystyle=\frac{2}{n}\left(\frac{\partial\boldsymbol{\gamma}_{\mathbf{t}}}{\partial
t_{j}}\right)^{{}^{\mkern-1.5mu\mathsf{T}}}\left[\left(X^{{}^{\mkern-1.5mu\mathsf{T}}}X\right)\boldsymbol{\gamma}_{\mathbf{t}}-\left(X^{{}^{\mkern-1.5mu\mathsf{T}}}\mathbf{y}\right)\right]+\lambda$
$\displaystyle=2\left(\frac{\partial\boldsymbol{\gamma}_{\mathbf{t}}}{\partial
t_{j}}\right)^{{}^{\mkern-1.5mu\mathsf{T}}}\mathbf{a}_{\mathbf{t}}+\lambda,$
(24)

where

$\mathbf{a}_{\mathbf{t}}=\left(\frac{X^{{}^{\mkern-1.5mu\mathsf{T}}}X}{n}\right)\boldsymbol{\gamma}_{\mathbf{t}}-\left(\frac{X^{{}^{\mkern-1.5mu\mathsf{T}}}\mathbf{y}}{n}\right).$

From the definition of $\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}$ and
$\boldsymbol{\gamma}_{\mathbf{t}}$,

$\displaystyle\frac{\partial\boldsymbol{\gamma}_{\mathbf{t}}}{\partial t_{j}}$
$\displaystyle=\frac{\partial
T_{\mathbf{t}}\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}}{\partial t_{j}}$
$\displaystyle=\frac{\partial T_{\mathbf{t}}}{\partial
t_{j}}\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}+T_{\mathbf{t}}\frac{\partial\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}}{\partial
t_{j}}$
$\displaystyle=E_{j}\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}+T_{\mathbf{t}}L_{\mathbf{t}}^{-1}\left[E_{j}-E_{j}ZT_{\mathbf{t}}L^{-1}_{\mathbf{t}}T_{\mathbf{t}}-T_{\mathbf{t}}ZE_{j}L^{-1}_{\mathbf{t}}T_{\mathbf{t}}\right]\left(\frac{X^{{}^{\mkern-1.5mu\mathsf{T}}}\mathbf{y}}{n}\right),$

where the last equality is obtained using Lemma 3 and the fact that
$\frac{\partial T_{\mathbf{t}}}{\partial t_{j}}=E_{j}$ and
$Z=\left(\frac{X^{{}^{\mkern-1.5mu\mathsf{T}}}X}{n}-\frac{\delta}{n}I\right)$.

Further simplification yields

$\displaystyle\frac{\partial\boldsymbol{\gamma}_{\mathbf{t}}}{\partial t_{j}}$
$\displaystyle=E_{j}\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}+T_{\mathbf{t}}L_{\mathbf{t}}^{-1}\left[E_{j}\left(\frac{X^{{}^{\mkern-1.5mu\mathsf{T}}}\mathbf{y}}{n}\right)-E_{j}Z\boldsymbol{\gamma}_{\mathbf{t}}-T_{\mathbf{t}}ZE_{j}\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}\right]$
$\displaystyle=E_{j}\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}-T_{\mathbf{t}}L_{\mathbf{t}}^{-1}E_{j}\mathbf{b}_{\mathbf{t}}-T_{\mathbf{t}}L_{\mathbf{t}}^{-1}T_{\mathbf{t}}ZE_{j}\widetilde{\boldsymbol{\beta}}_{\mathbf{t}},$
(25)

and we recall that

$\mathbf{b}_{\mathbf{t}}=Z\boldsymbol{\gamma}_{\mathbf{t}}-\left(\frac{X^{{}^{\mkern-1.5mu\mathsf{T}}}\mathbf{y}}{n}\right)=\mathbf{a}_{\mathbf{t}}-\frac{\delta}{n}\boldsymbol{\gamma}_{\mathbf{t}}.$

For a further simplification, recall that

$\displaystyle\mathbf{c}_{\mathbf{t}}=L^{-1}_{\mathbf{t}}\left(\mathbf{t}\odot\mathbf{a}_{\mathbf{t}}\right),$

and

$\displaystyle\mathbf{d}_{\mathbf{t}}=Z\left(\mathbf{t}\odot\mathbf{c}_{\mathbf{t}}\right).$

Then, from (25), the matrix
$\frac{\partial\boldsymbol{\gamma}_{\mathbf{t}}}{\partial{\mathbf{t}}}$ of
dimension $p\times p$ with $j$-th column being
$\frac{\partial\boldsymbol{\gamma}_{\mathbf{t}}}{\partial t_{j}}$ can be
expressed as

$\displaystyle\frac{\partial\boldsymbol{\gamma}_{\mathbf{t}}}{\partial{\mathbf{t}}}=\mathsf{Diag}\left(\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}\right)-T_{\mathbf{t}}L^{-1}_{\mathbf{t}}\mathsf{Diag}\left(\mathbf{b}_{\mathbf{t}}\right)-T_{\mathbf{t}}L_{\mathbf{t}}^{-1}T_{\mathbf{t}}Z\mathsf{Diag}\left(\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}\right).$
(26)

From (24), with $\mathbf{1}$ representing a vector of all ones, $\nabla
f_{\lambda}(\mathbf{t})$ can be expressed as

$\displaystyle\nabla f_{\lambda}(\mathbf{t})$
$\displaystyle=2\mathsf{Diag}\left(\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}\right)\mathbf{a}_{\mathbf{t}}-2\mathsf{Diag}\left({\mathbf{b}}_{\mathbf{t}}\right)L^{-1}_{\mathbf{t}}T_{\mathbf{t}}\mathbf{a}_{\mathbf{t}}-2\mathsf{Diag}\left(\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}\right)ZT_{\mathbf{t}}L_{\mathbf{t}}^{-1}T_{\mathbf{t}}\mathbf{a}_{\mathbf{t}}+\lambda\mathbf{1}$
$\displaystyle=2\left(\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}\odot{\mathbf{a}}_{\mathbf{t}}\right)-2\mathsf{Diag}\left({\mathbf{b}}_{\mathbf{t}}\right){\mathbf{c}}_{\mathbf{t}}-2\mathsf{Diag}\left(\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}\right)ZT_{\mathbf{t}}{\mathbf{c}}_{\mathbf{t}}+\lambda\mathbf{1}$
$\displaystyle=2\left(\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}\odot{\mathbf{a}}_{\mathbf{t}}\right)-2\left({\mathbf{b}}_{\mathbf{t}}\odot{\mathbf{c}}_{\mathbf{t}}\right)-2\left({\widetilde{\boldsymbol{\beta}}}_{\mathbf{t}}\odot{\mathbf{d}}_{\mathbf{t}}\right)+\lambda\mathbf{1}$
$\displaystyle=2\left(\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}\odot\left({\mathbf{a}}_{\mathbf{t}}-{\mathbf{d}}_{\mathbf{t}}\right)\right)-2\left({\mathbf{b}}_{\mathbf{t}}\odot{\mathbf{c}}_{\mathbf{t}}\right)+\lambda\mathbf{1}$
$\displaystyle={\boldsymbol{\zeta}}_{\mathbf{t}}+\lambda\mathbf{1},$

where

$\displaystyle{\boldsymbol{\zeta}}_{\mathbf{t}}=2\left(\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}\odot\left({\mathbf{a}}_{\mathbf{t}}-{\mathbf{d}}_{\mathbf{t}}\right)\right)-2\left({\mathbf{b}}_{\mathbf{t}}\odot{\mathbf{c}}_{\mathbf{t}}\right).$

To solve the optimization problem (14) using any gradient based method, we
need to compute the gradients of $g_{\lambda}(\mathbf{w})$:

$\nabla g_{\lambda}(\mathbf{w})=\left(\frac{\partial
g_{\lambda}(\mathbf{w})}{\partial w_{1}},\dots,\frac{\partial
g_{\lambda}(\mathbf{w})}{\partial
w_{p}}\right)^{{}^{\mkern-1.5mu\mathsf{T}}}.$

Recall that we define
$g_{\lambda}(\mathbf{w})=f_{\lambda}\left(\mathbf{t}(\mathbf{w})\right)$ in
(13) and

$f_{\lambda}(\mathbf{t})=\frac{1}{n}\|\mathbf{y}-X_{\mathbf{t}}\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}\|_{2}^{2}+\lambda\sum_{j=1}^{p}t_{j}.$

Then, from the chain rule of differentiation, for each $j=1,\dots,p$,

$\frac{\partial g_{\lambda}(\mathbf{w})}{\partial w_{j}}=\frac{\partial
f_{\lambda}(\mathbf{t})}{\partial t_{j}}\left(2w_{j}\exp(-w_{j}^{2})\right).$

Alternatively, in short,

$\displaystyle\nabla g_{\lambda}(\mathbf{w})=\nabla
f_{\lambda}(\mathbf{t})\odot\left(2\mathbf{w}\odot\exp(-\mathbf{w}\odot\mathbf{w})\right).$
(27)

∎

###### Proof of Proposition 1.

From Theorem 4, we obtain $\nabla f_{\lambda}(\mathbf{t})$ as follows,

$\nabla f_{\lambda}(\mathbf{t})=\zeta_{\mathbf{t}}+\lambda\mathbf{1},$

where $\zeta_{\mathbf{t}}\in\mathbb{R}^{p}$ (recalling from Theorem 4) is
given by

$\displaystyle\boldsymbol{\zeta}_{\mathbf{t}}=2\left({\widetilde{\boldsymbol{\beta}}}_{\mathbf{t}}\odot\left({\mathbf{a}}_{\mathbf{t}}-{\mathbf{d}}_{\mathbf{t}}\right)\right)-2\left({\mathbf{b}}_{\mathbf{t}}\odot{\mathbf{c}}_{\mathbf{t}}\right),$
(28)

with

$\displaystyle\mathbf{a}_{\mathbf{t}}$
$\displaystyle=\left(\frac{X^{{}^{\mkern-1.5mu\mathsf{T}}}X}{n}\right)\left({\mathbf{t}}\odot\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}\right)-\left(\frac{X^{{}^{\mkern-1.5mu\mathsf{T}}}\mathbf{y}}{n}\right),$
$\displaystyle\mathbf{b}_{\mathbf{t}}$
$\displaystyle=\mathbf{a}_{\mathbf{t}}-\frac{\delta}{n}\left({\mathbf{t}}\odot\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}\right),$
$\displaystyle\mathbf{c}_{\mathbf{t}}$
$\displaystyle=L^{-1}_{\mathbf{t}}\left({\mathbf{t}}\odot\mathbf{a}_{\mathbf{t}}\right),\quad\text{and}$
$\displaystyle\mathbf{d}_{\mathbf{t}}$
$\displaystyle=\left(\frac{X^{{}^{\mkern-1.5mu\mathsf{T}}}X}{n}-\frac{\delta}{n}I\right)\left({\mathbf{t}}\odot{\mathbf{c}}_{\mathbf{t}}\right).$

From the definition of $\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}$, it is
easy to show that for any $j$, if we fix $t_{i}$ for all $i\neq j$, then the
$j$-th component of $\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}$ goes to zero
as $t_{j}\downarrow 0$, that is, $\lim_{t_{j}\downarrow
0}\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}(j)=0.$ Similarly,
$\lim_{t_{j}\downarrow 0}\mathbf{c}_{\mathbf{t}}(j)=0.$ Therefore, from the
expression of $\boldsymbol{\zeta}_{\mathbf{t}}$ in (28), we have
$\lim_{t_{j}\downarrow 0}\boldsymbol{\zeta}_{\mathbf{t}}(j)=0.$ ∎

We need Lemma 4 for proving Theorem 5; we refer to [29] for a proof of Lemma
4.

######  Lemma 4 (Woodbury matrix identity or Duncan Inversion Formula).

For any conformable matrices $A,B_{1}$, and $B_{2},$ and $C$,

$\left(A+B_{1}CB_{2}\right)^{-1}=A^{-1}-A^{-1}B_{1}\left(C^{-1}+B_{2}A^{-1}B_{1}\right)^{-1}B_{2}A^{-1}.$

###### Proof of Theorem 5.

Now recall the expression of $L_{\mathbf{t}}$ from (6):

$\displaystyle L_{\mathbf{t}}$
$\displaystyle=\frac{1}{n}\left[X^{{{}^{\mkern-1.5mu\mathsf{T}}}}_{\mathbf{t}}X_{\mathbf{t}}+\delta\left(I-T_{\mathbf{t}}^{2}\right)\right].$

Since the gradient descent algorithm in COMBSS computes
$g_{\lambda}(\mathbf{w})$ for $\mathbf{w}\in\mathbb{R}^{p}$, from the map
$\mathbf{w}\mapsto\mathbf{t}(\mathbf{w})$, we care about solving
$L_{\mathbf{t}}\mathbf{z}=\mathbf{u}$ only for $\mathbf{t}\in[0,1)^{p}$. For
such a $\mathbf{t}$,

$S_{\mathbf{t}}:=\frac{n}{\delta}(I-T^{2}_{\mathbf{t}})^{-1}$

exists. In particular, $S_{\mathbf{t}}$ is a $p$-dimensional diagonal matrix
with the diagonal elements given by
$(\nicefrac{{n}}{{\delta(1-t_{1}^{2})}},\dots,\nicefrac{{n}}{{\delta(1-t_{p}^{2})}})$,
which is easy to compute. Further, if we take

$A=\frac{\delta}{n}\left(I-T_{\mathbf{t}}^{2}\right),\quad
B_{1}=B_{2}^{{}^{\mkern-1.5mu\mathsf{T}}}=\frac{1}{\sqrt{n}}X^{{{}^{\mkern-1.5mu\mathsf{T}}}}_{\mathbf{t}},\quad\text{and}\quad
C=I$

in Lemma 4, then,

$\displaystyle L_{\mathbf{t}}^{-1}$
$\displaystyle=S_{\mathbf{t}}-\frac{1}{n}S_{\mathbf{t}}X_{\mathbf{t}}^{{{}^{\mkern-1.5mu\mathsf{T}}}}\left(I+\frac{1}{n}X_{\mathbf{t}}S_{\mathbf{t}}X_{\mathbf{t}}^{{}^{\mkern-1.5mu\mathsf{T}}}\right)^{-1}X_{\mathbf{t}}S_{\mathbf{t}}.$

Since $S_{\mathbf{t}}$ is diagonal matrix, it is easy to see that

$\widetilde{L}_{\mathbf{t}}:=I+\frac{1}{n}X_{\mathbf{t}}S_{\mathbf{t}}X_{\mathbf{t}}^{{}^{\mkern-1.5mu\mathsf{T}}}$

is a symmetric positive-definite matrix of dimension $n\times n$ (in
comparison, $L_{\mathbf{t}}$ is of dimension $p\times p$). Thus,
$\widetilde{L}_{\mathbf{t}}^{-1}$ exists and

$\displaystyle L_{\mathbf{t}}^{-1}\mathbf{u}$
$\displaystyle=\left(S_{\mathbf{t}}\mathbf{u}\right)-\frac{1}{n}S_{\mathbf{t}}X_{\mathbf{t}}^{{{}^{\mkern-1.5mu\mathsf{T}}}}\widetilde{L}_{\mathbf{t}}^{-1}\left(X_{\mathbf{t}}S_{\mathbf{t}}\mathbf{u}\right).$

∎

###### Proof of Theorem 6.

For the same reasons mentioned in the proof of Theorem 2, without loss of
generality we can assume that all the zero elements of $\mathbf{t}$ appear at
the end of $\mathbf{t}$; that is, $\mathbf{t}$ is of the form:

$\mathbf{t}=(t_{1},\dots,t_{m},0,\dots,0),$

where $m$ indicates the number of non-zero elements in $\mathbf{t}$. Then,
$L_{\mathbf{t}}$ is given by

$\displaystyle L_{\mathbf{t}}$
$\displaystyle=\frac{1}{n}\left[\begin{pmatrix}\left(X_{\mathbf{t}}^{{}^{\mkern-1.5mu\mathsf{T}}}X_{\mathbf{t}}\right)_{\scriptscriptstyle{+}}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}\end{pmatrix}+\delta\begin{pmatrix}I-\left(T^{2}_{\mathbf{t}}\right)_{\scriptscriptstyle{+}}&\mathbf{0}\\\
\mathbf{0}&I\end{pmatrix}\right]$ (29)
$\displaystyle=\begin{pmatrix}\left(L_{\mathbf{t}}\right)_{\scriptscriptstyle{+}}&\mathbf{0}\\\
\mathbf{0}&\frac{\delta}{n}I\end{pmatrix}.$

Since $\mathbf{t}\in[0,1)^{p}$, from Theorem 1, $L_{\mathbf{t}}$ is a
positive-definite matrix. It is a well-known fact that every principle
submatrix of a positive-definite matrix is also positive-definite. Thus,
$\left(L_{\mathbf{t}}\right)_{\scriptscriptstyle{+}}$ is positive-definite and
hence invertible. Using Lemma 1,

$\displaystyle L^{-1}_{\mathbf{t}}$
$\displaystyle=\begin{pmatrix}\left(\left(L_{\mathbf{t}}\right)_{\scriptscriptstyle{+}}\right)^{-1}&\mathbf{0}\\\
\mathbf{0}&I\end{pmatrix}\begin{pmatrix}I&\mathbf{0}\\\
\mathbf{0}&\frac{n}{\delta}I\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\left(\left(L_{\mathbf{t}}\right)_{\scriptscriptstyle{+}}\right)^{-1}&\mathbf{0}\\\
\mathbf{0}&\frac{n}{\delta}I\end{pmatrix}.$ (30)

This completes the proof of (i).

Now observe that

$\left(\frac{X_{\mathbf{t}}^{{}^{\mkern-1.5mu\mathsf{T}}}\mathbf{y}}{n}\right)=\begin{pmatrix}\left(\frac{X_{\mathbf{t}}^{{}^{\mkern-1.5mu\mathsf{T}}}\mathbf{y}}{n}\right)_{\scriptscriptstyle{+}}\\\
\mathbf{0}\end{pmatrix}=\begin{pmatrix}(\mathbf{t})_{\scriptscriptstyle{+}}\odot\left(\frac{X^{{}^{\mkern-1.5mu\mathsf{T}}}\mathbf{y}}{n}\right)_{\scriptscriptstyle{+}}\\\
\mathbf{0}\end{pmatrix}.$

Since
$\widetilde{\boldsymbol{\beta}}_{\mathbf{t}}=L_{\mathbf{t}}^{-1}\left(\frac{X_{\mathbf{t}}^{{}^{\mkern-1.5mu\mathsf{T}}}\mathbf{y}}{n}\right)$,
using (30), we complete the proof of (ii). Similar arguments prove (iii);
hence for conciseness that proof is omitted here. ∎

##  Appendix B Gradient Descent Methods

In this section, we review two popular gradient descent methods, namely, the
basic gradient descent method and the Adam optimizer. In each iteration, the
basic gradient descent method uses only the current gradient for finding the
descent direction $\Delta\mathbf{w}$. On the other hand, the Adam optimizer
uses the current as well as all the past gradients for finding
$\Delta\mathbf{w}$, in such a way that the algorithm moves with momentum that
can provide acceleration and can help in bypassing saddle points and local
minima of the objective function. See, for example, [18], for more details on
gradient descent methods.

###  B.1 Basic Gradient Descent Method

For each iteration $l=1,2,\dots,$ let $\mathbf{w}^{(l-1)}$ be the value of
$\mathbf{w}$ obtained from the previous iteration; here, the first iteration
starts with the initial vector $\mathbf{w}^{(0)}$ given to the algorithm.
Then, the next point is obtained from the update:

$\mathbf{w}^{(l)}=\mathbf{w}^{(l-1)}-\alpha\,\nabla
g_{\lambda}(\mathbf{w}^{(l-1)}),$

for a positive constant $\alpha$ known as learning rate. That means, the
descent step $\Delta\mathbf{w}$ in the $l$-th iteration is $-\alpha\,\nabla
g_{\lambda}(\mathbf{w}^{(l-1)})$, which is in the negative gradient direction,
i.e., the most rapid decreasing direction. The learning parameter $\alpha$
controls the convergence rate: smaller $\alpha$ would slow convergence and
larger $\alpha$ would lead to divergence. To overcome the difficulty of
selecting the learning parameter and increase the convergence speed, there are
several methods proposed in the existing literature; among them, most popular
one is the Adam optimizer, which is presented below.

###  B.2 The Adam Optimizer

The Adam (short for Adaptive Moment Estimation) optimizer [17] uses the
following exponential smoothing of the gradients and the squared gradients of
the objective function $g_{\lambda}(\mathbf{w})$. For each iteration
$l=1,2,\dots,$ with $\mathbf{w}^{(l-1)}$ being the value of $\mathbf{w}$
obtained from the previous iteration, let

$\displaystyle\mathbf{u}^{(l)}$
$\displaystyle=-(1-\xi_{1})\sum_{k=0}^{l-1}\xi_{1}^{l-k-1}\nabla
g_{\lambda}(\mathbf{w}^{(k)}),\quad(\text{momentum update})$
$\displaystyle\mathbf{v}^{(l)}$
$\displaystyle=(1-\xi_{2})\sum_{k=0}^{l-1}\xi_{2}^{l-k-1}\left(\nabla
g_{\lambda}(\mathbf{w}^{(k)})\odot\nabla
g_{\lambda}(\mathbf{w}^{(k)})\right),$

with $\mathbf{u}^{(0)}=\mathbf{v}^{(0)}=\mathbf{0}$ and
$\xi_{1},\xi_{2}\in[0,1)$. The default choices are $\xi_{1}=0.9$ and
$\xi_{2}=0.999$. Further take,

$\displaystyle\widehat{\mathbf{u}}^{(l)}=\frac{\mathbf{u}^{(l)}}{(1-\xi_{1}^{l})},$

and

$\displaystyle\widehat{\mathbf{v}}^{(l)}=\frac{\mathbf{v}^{(l)}}{(1-\xi_{2}^{l})}.$

Finally, the $\mathbf{w}$ vector is updated using

$\mathbf{w}^{(l)}=\mathbf{w}^{(l-1)}+\alpha\frac{\widehat{\mathbf{u}}^{(l)}}{\sqrt{\widehat{\mathbf{v}}^{(l)}+c}},\quad
l=1,2,\dots,$

where the division, multiplication and addition are element-wise, and $c$ is a
constant of order $1\times 10^{-8}$ to avoid division by zero and $\alpha$ is
a positive constant. This means, the descent step $\Delta\mathbf{w}$ in the
$l$-th iteration is $\alpha\mathbf{u}^{(l)}/\sqrt{\mathbf{v}^{(l)}+c}$, which
depends on all the gradients computed up to the $l$-th iteration.

Observe that the updates $\mathbf{u}$ and $\mathbf{v}$ can be rewritten as the
following easy to implement recursive expressions:

$\displaystyle\begin{aligned}
\mathbf{u}^{(l)}&=\xi_{1}\mathbf{u}^{(l-1)}-(1-\xi_{1})\nabla
g_{\lambda}(\mathbf{w}^{(l-1)}),\\\
\mathbf{v}^{(l)}&=\xi_{2}\mathbf{v}^{(l-1)}+(1-\xi_{2})\left(\nabla
g_{\lambda}(\mathbf{w}^{(l-1)})\odot\nabla
g_{\lambda}(\mathbf{w}^{(l-1)})\right).\end{aligned}$ (31)

##  Appendix C Additional Simulation Results

Figure 11: Performance results in terms of F1-score, Sensibility, and
Specificity, for the low-dimensional cases: $n=100$, $p=20$, $\rho=0.8$, true
active predictor is 10 with Beta-type 1 (top plots) and Beta-type 2 (bottom
plots).

Figure 12: Performance results in terms of F1-score, Sensibility, and
Specificity, for the high-dimensional cases: $n=100$, $p=1000$, $\rho=0.8$,
true active predictor is 10 with Beta-type 1 (top plots) and Beta-type 2
(bottom plots).

## References

  * [1] H. Akaike.  Information theory and an extension of the maximum likelihood principle.  In Second International Symposium on Information Theory (Tsahkadsor, 1971), pages 267–281. 1973.
  * [2] M. Artin.  Algebra.  Pearson Prentice Hall, 2011.
  * [3] A. Ben-Israel and A. Charnes.  Contributions to the theory of generalized inverses.  J. Soc. Indust. Appl. Math., 11:667–699, 1963.
  * [4] D. Bertsimas, A. King, and R. Mazumder.  Best subset selection via a modern optimization lens.  The Annals of Statistics, 44(2):813 – 852, 2016.
  * [5] J. Chen and Z. Chen.  Extended Bayesian information criteria for model selection with large model spaces.  Biometrika, 95(3):759–771, 2008.
  * [6] D. Chicco and G. Jurman.  The advantages of the matthews correlation coefficient (mcc) over f1 score and accuracy in binary classification evaluation.  BMC genomics, 21(1):1–13, 2020.
  * [7] M. A. Efroymson.  Stepwise regression—a backward and forward look.  Presented at the Eastern Regional Meetings of of the Institute of Mathematical Statistics, Florham Park, New Jersey, 1966.
  * [8] J. Fan and R. Li.  Statistical challenges with high dimensionality: feature selection in knowledge discovery.  In International Congress of Mathematicians. Vol. III, pages 595–622. Eur. Math. Soc., Zürich, 2006.
  * [9] J. Fan and J. Lv.  A selective overview of variable selection in high dimensional feature space.  Statist. Sinica, 20(1):101–148, 2010.
  * [10] G. M. Furnival and R. W. Wilson.  Regressions by leaps and bounds.  Technometrics, 42:69–79, 2000.
  * [11] G. H. Golub and C. F. Van Loan.  Matrix computations.  Johns Hopkins Studies in the Mathematical Sciences. Johns Hopkins University Press, Baltimore, MD, third edition, 1996.
  * [12] Gurobi Optimization, LLC.  Gurobi Optimizer Reference Manual, 2022.
  * [13] T. Hastie, R. Tibshirani, and R. Tibshirani.  bestsubset: Tools for best subset selection in regression, 2018\.  R package version 1.0.10.
  * [14] T. Hastie, R. Tibshirani, and R. Tibshirani.  Best subset, forward stepwise or lasso? analysis and recommendations based on extensive comparisons.  Statistical Science, 35(4):579–592, 2020.
  * [15] H. Hazimeh and R. Mazumder.  Fast best subset selection: coordinate descent and local combinatorial optimization algorithms.  Oper. Res., 68(5):1517–1537, 2020.
  * [16] R. R. Hocking and R. N. Leslie.  Selection of the best subset in regression analysis.  Technometrics, 9:531–540, 1967.
  * [17] D. P. Kingma and J. Ba.  Adam: A method for stochastic optimization.  In Y. Bengio and Y. LeCun, editors, 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings, 2015.
  * [18] M. J. Kochenderfer and T. A. Wheeler.  Algorithms for optimization.  MIT Press, Cambridge, MA, 2019.
  * [19] T. Lumley.  leaps: Regression Subset Selection, 2020.  R package version 3.1, based on Fortran code by Alan Miller.
  * [20] N. Meinshausen.  Relaxed Lasso.  Comput. Statist. Data Anal., 52(1):374–393, 2007.
  * [21] A. Miller.  Subset selection in regression, volume 95 of Monographs on Statistics and Applied Probability.  Chapman & Hall/CRC, Boca Raton, FL, 2019.
  * [22] S. Müller and A. H. Welsh.  On model selection curves.  Intnl. Statist. Reviews, 78(2):240–256, 2010.
  * [23] B. K. Natarajan.  Sparse approximate solutions to linear systems.  SIAM J. Comput., 24(2):227–234, 1995.
  * [24] G. Schwarz.  Estimating the dimension of a model.  Ann. Statist., 6(2):461–464, 1978.
  * [25] Y. Takano and R. Miyashiro.  Best subset selection via cross-validation criterion.  TOP, 28:475–488, 2020.
  * [26] G. Tarr, S. Muller, and A. H. Welsh.  mplot: An r package for graphical model stability and variable selection procedures.  J. Statist. Software, 83(9):1–28, 2018.
  * [27] Y. Tian and Y. Takane.  Schur complements and Banachiewicz-Schur forms.  Electron. J. Linear Algebra, 13:405–418, 2005.
  * [28] R. Tibshirani.  Regression shrinkage and selection via the lasso.  J. Roy. Statist. Soc. Ser. B, 58(1):267–288, 1996.
  * [29] M. A. Woodbury.  Inverting modified matrices.  Princeton University, Princeton, N. J., 1950.  Statistical Research Group, Memo. Rep. no. 42,.
  * [30] J. Zhu, C. Wen, J. Zhu, H. Zhang, and X. Wang.  A polynomial algorithm for best-subset selection problem.  Proc. Natl. Acad. Sci. USA, 117(52):33117–33123, 2020.