# Robust Information Criterion for Model Selection in Sparse High-Dimensional
Linear Regression Models

Prakash B. Gohain, Magnus Jansson This research was supported in part by the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme, grant agreement No. 742648. The authors are
with the Division of Information Science and Engineering, KTH Royal Institute
of Technology, Stockholm SE-10044, Sweden e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Model selection in linear regression models is a major challenge when dealing
with high-dimensional data where the number of available measurements (sample
size) is much smaller than the dimension of the parameter space. Traditional
methods for model selection such as Akaike information criterion, Bayesian
information criterion (BIC) and minimum description length are heavily prone
to overfitting in the high-dimensional setting. In this regard, extended BIC
(EBIC), which is an extended version of the original BIC and extended Fisher
information criterion (EFIC), which is a combination of EBIC and Fisher
information criterion, are consistent estimators of the true model as the
number of measurements grows very large. However, EBIC is not consistent in
high signal-to-noise-ratio (SNR) scenarios where the sample size is fixed and
EFIC is not invariant to data scaling resulting in unstable behaviour. In this
paper, we propose a new form of the EBIC criterion called EBIC-Robust, which
is invariant to data scaling and consistent in both large sample size and
high-SNR scenarios. Analytical proofs are presented to guarantee its
consistency. Simulation results indicate that the performance of EBIC-Robust
is quite superior to that of both EBIC and EFIC.

###### Index Terms:

High-dimension, linear regression, data scaling, statistical model selection,
subset selection, sparse estimation, scale-invariant, variable selection.

##  I Introduction

Selecting the true or best set of covariates from a large pool of potential
covariates is a fundamental requirement in many applications of science,
engineering and biology. In this paper, our primary focus is on model
selection in high-dimensional linear regression models associated with the
maximum likelihood (ML) method of parameter estimation where the number of
measurements, $N$, is quite small compared to the model space or parameter
dimension, $p$, i.e., $N<p$. High-dimensional datasets are a common phenomena
in many fields of scientific studies, and as such model selection is a central
element of data analysis and statistical inference [1].

Consider the linear model

$\mathbf{y}=\mathbf{A}\mathbf{x}+\mathbf{e},$ (1)

where $\mathbf{y}\in\mathbb{R}^{N}$ is the measurement vector and
$\mathbf{A}\in\mathbb{R}^{N\times p}$ is the known design matrix. We are
considering a high-dimensional setting, hence $p>N$. Also, $p$ can be linked
to $N$ as $p=N^{d}$, where $d>0$ is a real value.
$\mathbf{e}\in\mathbb{R}^{N}$ is the associated noise vector whose elements
are assumed to be i.i.d. following a Gaussian distribution, i.e.,
$\mathbf{e}\sim\mathcal{N}(\mathbf{0},\sigma^{2}\mathbf{I}_{N})$ where
$\sigma^{2}$ is the unknown true noise power. $\mathbf{x}\in\mathbb{R}^{p}$ is
the unknown parameter vector. Here, $\mathbf{x}$ is assumed to be sparse,
which implies that very few of the elements of $\mathbf{x}$ are non-zero. We
denote $\mathcal{S}$ as the true support of $\mathbf{x}$, i.e.,
$\mathcal{S}=\\{i:x_{i}\neq 0\\}$ having cardinality
$card(\mathcal{S})=k_{0}\ll N$ and $\mathbf{A}_{\mathcal{S}}$ as the set of
columns of $\mathbf{A}$ corresponding to the support $\mathcal{S}$. The goal
of model selection is estimating $\mathcal{S}$ given $\mathbf{y}$ and
$\mathbf{A}$.

A popular approach for model selection is using information theoretic criteria
[2, 3, 4, 5]. A typical information criterion based model selection rule picks
the best model that minimizes some statistical metric as shown below

$\hat{\mathcal{S}}=\underset{\mathcal{I}\in\mathcal{J}}{\arg\min}\\{f(\mathcal{M}_{\mathcal{I}})+\mathcal{P}(\mathcal{I})\\},$
(2)

where $\hat{\mathcal{S}}$ is the model estimate, $\mathcal{J}$ is the set of
candidate models under consideration and $\mathcal{M}_{\mathcal{I}}$ denotes
the model with support $\mathcal{I}$. The statistical metric consists of two
parts: (1) $f(\mathcal{M}_{\mathcal{I}})$ representing the goodness of fit of
model $\mathcal{M}_{\mathcal{I}}$ and (2) $\mathcal{P}(\mathcal{I})$ is the
penalty term that compensates for overparameterization. The literature on
model selection is quite extensive. Some of the popular classical model
selection rules include Akaike information criterion [6], Bayesian information
criterion (BIC)[7], minimum description length (MDL)[8], gMDL[9], nMDL[10],
and penalizing adaptively the likelihood (PAL) [11]. However, these classical
methods in their current form fail to handle the large dimension cases and
tend to overfit the final model [12, 13].

Among the classical methods of model selection, BIC has been quite successful
due to its simplicity and consistent performance in many fields. BIC is
asymptotically consistent in selecting the true model as $N$ grows very large
given that $p$ and the true noise variance $\sigma^{2}$ is fixed. However, its
performance in high-dimensional settings when $p>N$ is not satisfactory and it
has a tendency to select more co-variates than required, thus overfitting the
model [12]. To handle the large-$p$ small-$N$ scenario, the authors in [12]
proposed a novel extension to the original BIC called extended BIC (EBIC),
that takes into account both the number of unknown parameters and the
complexity of the model space. EBIC adds dynamic prior model probabilities to
each of the models under consideration that is inversely proportional to the
model set dimension. This eliminates the earlier assumption of assigning
uniform prior to all models irrespective of their sizes, which goes against
the principle of parsimony. Under a suitable asymptotic identifiability
condition, EBIC is consistent such that it selects the true model as $N$ tends
to infinity [12]. However, the consistent behaviour of EBIC fails when $N$ is
small and fixed and $\sigma^{2}$ tends to zero [13]. This new consistency
requirement was first introduced in [14], where the authors highlighted that
the original BIC is also inconsistent for fixed $N$ and decreasing noise
variance scenarios where $N>p$.

To overcome the drawbacks of EBIC, the authors in [13] proposed a novel
criterion called extended Fisher information criterion (EFIC) that is inspired
by EBIC and the model selection criteria with Fisher information [15]. The
authors analyzed the performance of EFIC in the high-dimensional setting for
two key cases: (1) when $\sigma^{2}$ is fixed and $N$ tends to infinity; (2)
when $N$ is fixed and $\sigma^{2}$ tends to zero. In each case, it was shown
that EFIC selects the true model with a probability approaching one. However,
as indicated in our simulations, EFIC is not invariant to data scaling and it
tends to suffer from overfitting issues (and sometimes underfitting) in
practical sizes of $N$ when the data is scaled. This scaling problem is a
result of the data dependent penalty design that may blow the penalty to
extremely small or large values depending on how the data is scaled.

Apart from the criteria mentioned above, there are other non-information
theoretic methods available for model selection. One such popular method is
cross-validation (CV) [16]. However, the performance of CV is quite poor in
sample scarce scenarios with large parameter dimensions and even though CV is
unbiased, it can have high variance [17]. Recent additions to the list of
model selection methods for high-dimensional data are residual ratio
thresholding (RRT)[18] and multi-beta-test (MBT) [19]. Both are non-
information theoretic methods based on hypothesis testing using a test
statistic. They operate along with a greedy variable selection method such as
orthogonal matching pursuit (OMP) [20] and involve a tuning parameter
$\in[0,1]$, that controls the false-discovery rate. However, there is no
optimal way to set it and as such their behaviour may tend to overfit or
underfit depending on the chosen tuning parameter value. Moreover, in their
current form, they can only be used with algorithms that generate monotonic
sequences of support estimates such as OMP, which restricts their usability.

In this paper, we propose a modified criterion for model selection in high-
dimensional linear regression models called EBIC-Robust or EBIC${}_{\text{R}}$
in short, where the subscript R stands for robust. EBIC${}_{\text{R}}$ is a
scale-invariant and consistent criterion. To guarantee the consistency, we
provide analytical proofs to show that under a suitable asymptotic
identifiability condition, EBIC${}_{\text{R}}$ selects the true model with a
probability approaching one as $N\to\infty$ as well as when $\sigma^{2}\to 0$.
Some preliminary results have been published in [21].

Throughout the paper, boldface letters denote matrices and vectors. The
notation $(\cdot)^{T}$ stands for transpose. $\mathbf{A}_{\mathcal{I}}$
denotes a sub-matrix of the full matrix $\mathbf{A}$ formed using the columns
indexed by the support set $\mathcal{I}$.
$\boldsymbol{\Pi}_{\mathcal{I}}=\mathbf{A}_{\mathcal{I}}(\mathbf{A}_{\mathcal{I}}^{T}\mathbf{A}_{\mathcal{I}})^{-1}\mathbf{A}_{\mathcal{I}}^{T}$
denotes the orthogonal projection matrix on the span of
$\mathbf{A}_{\mathcal{I}}$ and
$\boldsymbol{\Pi}_{\mathcal{I}}^{\perp}=\mathbf{I}_{N}-\boldsymbol{\Pi}_{\mathcal{I}}$
denotes the orthogonal projection matrix on the null space of
$\mathbf{A}_{\mathcal{I}}^{T}$ and $\mathbf{I}_{N}$ is an $N\times N$ identity
matrix. The notation $\big{\lvert}\mathbf{X}\big{\rvert}$ denotes the
determinant of the matrix $\mathbf{X}$ and $\lVert\cdot\rVert_{2}$ denotes the
Euclidean norm. $X\sim\mathcal{N}(0,1)$ denotes a normal distributed random
variable with mean $0$ and variance $1$. $X\sim\chi^{2}_{k}$ is a central chi-
squared distributed random variable with $k$ degrees of freedom,
$X\sim\chi^{2}_{k}(\lambda)$ is a noncentral chi-squared distributed random
variable with $k$ degrees of freedom and non-centrality parameter $\lambda$.

##  II Background

Given the linear model (1), the entire process of model selection or in other
words estimating the true support set $\mathcal{S}$ involves two major steps:
(i) Predictor/subset selection, which includes finding a competent set of
candidate models out of all the ($2^{p}-1$) possible models. In our work, we
consider the set of competing models as the collection of all plausible
combinatorial models up to a maximum cardinality $K$, under the assumption
that $k_{0}\leq K\ll N$; (ii) estimating the true model among the candidate
models using a suitable model selection criterion. For a candidate model with
support $\mathcal{I}$ having cardinality $card(\mathcal{I})=k$, the linear
model in (1) can be reformulated as follows

$\mathcal{H}_{\mathcal{I}}:\mathbf{y}=\mathbf{A}_{\mathcal{I}}\mathbf{x}_{\mathcal{I}}+\mathbf{e}_{\mathcal{I}},$
(3)

where $\mathcal{H}_{\mathcal{I}}$ denotes the hypothesis that the data
$\mathbf{y}$ is truly generated according to (3),
$\mathbf{A}_{\mathcal{I}}\in\mathbb{R}^{N\times k}$ is the sub-design matrix
consisting of columns from the known design matrix $\mathbf{A}$ with support
$\mathcal{I}$, $\mathbf{x}_{\mathcal{I}}\in\mathbb{R}^{k}$ is the
corresponding unknown parameter vector and
$\mathbf{e}_{\mathcal{I}}\in\mathbb{R}^{N}$ is the associated noise vector
following
$\mathbf{e}_{\mathcal{I}}\sim\mathcal{N}(\mathbf{0},\sigma^{2}_{\mathcal{I}}\mathbf{I}_{N})$
where $\sigma^{2}_{\mathcal{I}}$ is the unknown noise variance corresponding
to the hypothesis $\mathcal{H}_{\mathcal{I}}$.

###  II-A Bayesian Framework for Model Selection

To motivate the proposed criterion we start by describing the Bayesian
framework that leads to the maximum a-posteriori (MAP) estimator, which in
turn forms the backbone for deriving BIC and its extended versions, viz.,
EBIC, EFIC, as well as the proposed criterion EBIC${}_{\text{R}}$. Now, for
the considered model in (3), the probability density function (pdf) of the
data vector $\mathbf{y}$ is given as

$p(\mathbf{y}\lvert\boldsymbol{\theta}_{\mathcal{I}},\mathcal{H}_{\mathcal{I}})=\frac{\exp\\{-\lVert\mathbf{y}-\mathbf{A}_{\mathcal{I}}\mathbf{x}_{\mathcal{I}}\rVert^{2}_{2}/2\sigma^{2}_{\mathcal{I}}\\}}{(2\pi\sigma^{2}_{\mathcal{I}})^{N/2}},$
(4)

where
$\boldsymbol{\theta}_{\mathcal{I}}=[\mathbf{x}^{T}_{\mathcal{I}},\sigma^{2}_{\mathcal{I}}]^{T}$
comprises of all the parameters of the model. Under hypothesis
$\mathcal{H}_{\mathcal{I}}$, the maximum likelihood estimates (MLEs) of
$\hat{\boldsymbol{\theta}}_{\mathcal{I}}=[\hat{\mathbf{x}}^{T}_{\mathcal{I}},\hat{\sigma}^{2}_{\mathcal{I}}]^{T}$
are obtained as [22]

$\hat{\mathbf{x}}_{\mathcal{I}}=\left(\mathbf{A}^{T}_{\mathcal{I}}\mathbf{A}_{\mathcal{I}}\right)^{-1}\mathbf{A}^{T}_{\mathcal{I}}\mathbf{y}\quad\&\quad\hat{\sigma}^{2}_{\mathcal{I}}=\frac{\mathbf{y}^{T}\mathbf{\Pi}^{\perp}_{\mathcal{I}}\mathbf{y}}{N}.$
(5)

Let $p(\boldsymbol{\theta}_{\mathcal{I}}|\mathcal{H}_{\mathcal{I}})$ denote
the prior pdf of the parameter vector $\boldsymbol{\theta}_{\mathcal{I}}$
under $\mathcal{H}_{\mathcal{I}}$. Then we have the joint probability

$p(\mathbf{y},\boldsymbol{\theta}_{\mathcal{I}}|\mathcal{H}_{\mathcal{I}})=p(\mathbf{y}|\boldsymbol{\theta}_{\mathcal{I}},\mathcal{H}_{\mathcal{I}})p(\boldsymbol{\theta}_{\mathcal{I}}|\mathcal{H}_{\mathcal{I}})$
(6)

and the marginal distribution of $\mathbf{y}$ is

$p(\mathbf{y}|\mathcal{H}_{\mathcal{I}})=\int
p(\mathbf{y}|\boldsymbol{\theta}_{\mathcal{I}},\mathcal{H}_{\mathcal{I}})p(\boldsymbol{\theta}_{\mathcal{I}}|\mathcal{H}_{\mathcal{I}})d\boldsymbol{\theta}_{\mathcal{I}}.$
(7)

The posterior probability $\Pr(\mathcal{H}_{\mathcal{I}}|\mathbf{y})$ is given
by

$\Pr(\mathcal{H}_{\mathcal{I}}|\mathbf{y})=\frac{p(\mathbf{y}|\mathcal{H}_{\mathcal{I}})\Pr\left(\mathcal{H}_{\mathcal{I}}\right)}{p(\mathbf{y})},$
(8)

where $\Pr(\mathcal{H}_{\mathcal{I}})$ is the prior probability of the model
with support $\mathcal{I}$. The MAP estimator picks the model with the largest
posterior probability $\Pr(\mathcal{H}_{\mathcal{I}}|\mathbf{y})$. However,
note that the $p(\mathbf{y})$ is a normalizing factor and independent of
$\mathcal{I}$. Hence, the MAP estimate of $\mathcal{S}$ is equivalently given
by

$\hat{\mathcal{S}}_{\text{MAP}}=\underset{\mathcal{I}}{\arg\max}\ \Big{\\{}\ln
p(\mathbf{y}|\mathcal{H}_{\mathcal{I}})+\ln\Pr\left(\mathcal{H}_{\mathcal{I}}\right)\Big{\\}}.$
(9)

To compute the MAP estimate, we need to evaluate the integral in (7).
Traditionally, under the assumption that $N$ and/or SNR are large, we can
obtain an approximation of $\ln p(\mathbf{y}|\mathcal{H}_{\mathcal{I}})$ using
a second order Taylor series expansion, which gives (see [23, 24] for details)

$\begin{split}\ln p(\mathbf{y}|\mathcal{H}_{\mathcal{I}})\approx\ln
p(\mathbf{y}|\hat{\boldsymbol{\theta}}_{\mathcal{I}},\mathcal{H}_{\mathcal{I}})+\ln
p(\hat{\boldsymbol{\theta}}_{\mathcal{I}}|\mathcal{H}_{\mathcal{I}})\\\
+\frac{k+1}{2}\ln(2\pi)-\frac{1}{2}\ln\big{|}\hat{\mathbf{F}}_{\mathcal{I}}\big{|},\end{split}$
(10)

where $k={card}(\mathcal{I})$ and $\hat{\mathbf{F}}_{\mathcal{I}}$ is the
sample Fisher information matrix under $\mathcal{H}_{\mathcal{I}}$ given as
[22]

$\hat{\mathbf{F}}_{\mathcal{I}}=-\frac{\partial^{2}\ln
p(\mathbf{y}|{\boldsymbol{\theta}}_{\mathcal{I}},\mathcal{H}_{\mathcal{I}})}{\partial\boldsymbol{\theta}_{\mathcal{I}}\partial\boldsymbol{\theta}_{\mathcal{I}}^{T}}\bigg{\lvert}_{{\boldsymbol{\theta}}_{\mathcal{I}}=\hat{\boldsymbol{\theta}}_{\mathcal{I}}}.$
(11)

Evaluating (11) using (4) and (5) we get [23]

$\hat{\mathbf{F}}_{\mathcal{I}}=\begin{bmatrix}\frac{1}{\hat{\sigma}^{2}_{\mathcal{I}}}\mathbf{A}_{\mathcal{I}}^{T}\mathbf{A}_{\mathcal{I}}&\mathbf{0}\\\
\mathbf{0}&\frac{N}{2\hat{\sigma}^{4}_{\mathcal{I}}}\end{bmatrix}.$ (12)

Now, for the considered linear model we have

$-2\ln
p(\mathbf{y}|\boldsymbol{\hat{\theta}}_{\mathcal{I}},\mathcal{H}_{\mathcal{I}})=N\ln\hat{\sigma}^{2}_{\mathcal{I}}+\text{const}.$
(13)

Therefore, using (13), we can rewrite (10) as

$\begin{split}-2\ln p(\mathbf{y}|\mathcal{H}_{\mathcal{I}})\approx
N\ln\hat{\sigma}^{2}_{\mathcal{I}}+\ln\big{|}\hat{\mathbf{F}}_{\mathcal{I}}\big{|}-2\ln
p(\hat{\boldsymbol{\theta}}_{\mathcal{I}}|\mathcal{H}_{\mathcal{I}})\\\ -k\ln
2\pi+\text{const.}\end{split}$ (14)

Furthermore, it is assumed that the prior term in (10), i.e., $\ln
p(\hat{\boldsymbol{\theta}}_{\mathcal{I}}|\mathcal{H}_{\mathcal{I}})$ is flat
and uninformative, and hence disregarded from the analysis. Thus, dropping the
constants and the terms independent of the model dimension $k$, we can
equivalently reformulate the MAP based model estimate as

$\hat{\mathcal{S}}_{\text{MAP}}=\underset{\mathcal{I}}{\arg\min}\Big{\\{}N\ln\hat{\sigma}^{2}_{\mathcal{I}}+\ln\big{|}\hat{\mathbf{F}}_{\mathcal{I}}\big{|}-k\ln
2\pi-2\ln\Pr\left(\mathcal{H}_{\mathcal{I}}\right)\Big{\\}}.$ (15)

###  II-B BIC

The BIC can be obtained from the MAP estimator in (15). The term $-k\ln 2\pi$
is ignored as it weakly depends on the model dimension $k$ and hence is
typically much smaller than the dominating terms. Moreover, the prior
probability of each candidate model is assumed to be equiprobable. Hence, the
$-2\ln\Pr(\mathcal{H}_{\mathcal{I}})$ term is dropped as well. Now, expanding
the $|\hat{\mathbf{F}}_{\mathcal{I}}|$ term of (15) using (12) we have

$\ln\big{|}\hat{\mathbf{F}}_{\mathcal{I}}\big{|}=\ln(N/2)-(k+2)\ln\hat{\sigma}^{2}_{\mathcal{I}}+\ln\left\lvert\mathbf{A}_{\mathcal{I}}^{T}\mathbf{A}_{\mathcal{I}}\right\rvert.$
(16)

Here, the following property of the design matrix $\mathbf{A}$ is assumed [23,
25]

$\lim_{N\to\infty}\left\\{N^{-1}(\mathbf{A}^{T}_{\mathcal{I}}\mathbf{A}_{\mathcal{I}})\right\\}=\mathbf{M}_{\mathcal{I}}=\mathcal{O}(1),$
(17)

where $\mathbf{M}_{\mathcal{I}}$ is a $k\times k$ positive definite matrix and
bounded as $N\to\infty$. The assumption in (17) is true in many applications
but not all (see [26] for more details). Using (17), it is possible to show
that for large $N$

$\ln\big{|}\mathbf{A}_{\mathcal{I}}^{T}\mathbf{A}_{\mathcal{I}}\big{|}=\ln\bigg{|}N\cdot{N^{-1}(\mathbf{A}_{\mathcal{I}}^{T}\mathbf{A}_{\mathcal{I}})}\bigg{|}=k\ln
N+\mathcal{O}(1).$ (18)

Furthermore, $\hat{\sigma}^{2}_{\mathcal{I}}$ is considered to be of
$\mathcal{O}(1)$ as well since it does not grow with $N$. As such, the
$\mathcal{O}(1)$ term, $(k+2)\ln\hat{\sigma}^{2}_{\mathcal{I}}$ and $\ln(N/2)$
(a constant) are ignored from (16). This leads to the final form of the BIC

$\text{BIC}(\mathcal{I})=N\ln\hat{\sigma}^{2}_{\mathcal{I}}+k\ln N.$ (19)

BIC is consistent when $p$ is fixed and $N\to\infty$. However, it is
inconsistent when $N$ is fixed and $\sigma^{2}\to 0$ [27, 24] as well as when
$p>N$ and $p$ grows exponentially with $N$ [12].

###  II-C EBIC

The authors in [12] proposed an extended version of the BIC, i.e., EBIC, to
mitigate the drawbacks of BIC for large-$p$ small-$N$ scenarios. EBIC can be
derived from the MAP estimator in (15), using the same assumptions as in BIC,
except for the prior probability term $\Pr(\mathcal{H}_{\mathcal{I}})$. In
EBIC, the idea of equiprobable models is discredited and instead a prior
probability is assigned that is inversely proportional to the size of the
model space. Thus, a model with dimension $k$ is assigned prior probability of
$\Pr(\mathcal{H}_{\mathcal{I}})\propto{p\choose k}^{-\gamma}$, where
$0\leq\gamma\leq 1$ is a tuning parameter. Thus, the EBIC is

$\text{EBIC}(\mathcal{I})=N\ln\hat{\sigma}^{2}_{\mathcal{I}}+k\ln
N+2\gamma\ln{p\choose k}.$ (20)

When $\gamma=0$, EBIC boils down to BIC (19). Moreover, unlike BIC, EBIC is
consistent in selecting the true model for $p\gg N$ cases where $p$ grows
exponentially with $N$. However, it has been observed in [13] that EBIC is
inconsistent when $N$ is fixed and $\sigma^{2}\to 0$.

###  II-D EFIC

To circumvent the shortcomings of EBIC in high-SNR cases, the authors in [13]
proposed EFIC. In EFIC, the assumptions imposed on the sample FIM (16) are
removed and the entire structure is included as it is in the criterion except
for the constant term $\ln(N/2)$. Some further simplifications are involved:

$\displaystyle
N\ln\hat{\sigma}^{2}_{\mathcal{I}}=N\ln\big{\lVert}\mathbf{\Pi}^{\perp}_{\mathcal{I}}\mathbf{y}\big{\rVert}^{2}_{2}-N\ln
N$ (21)
$\displaystyle(k+2)\ln\hat{\sigma}^{2}_{\mathcal{I}}=(k+2)\left[\ln\big{\lVert}\mathbf{\Pi}^{\perp}_{\mathcal{I}}\mathbf{y}\big{\rVert}^{2}_{2}-\ln
N\right].$ (22)

The $-N\ln N$ and $-2\ln N$ term of (21) and (22) respectively are independent
of the model dimension $k$ and hence ignored. Similar to EBIC the prior
probability term is assumed to be proportional to the model space, hence
$\Pr(\mathcal{H}_{\mathcal{I}})\propto{p\choose k}^{-c}$, where $c>0$ is a
tuning parameter. Furthermore, under the large-$p$ approximation and since
$k\leq K\ll p$, the $\ln{p\choose k}$ term is approximated as

$\ln{p\choose k}=\sum_{i=0}^{k-1}\ln(p-i)-\ln(k!)\approx k\ln p.$ (23)

Hence, for large-$p$ case, we can set $-2\ln
p(\mathcal{H}_{\mathcal{I}})\approx 2ck\ln p.$ Thus, the EFIC is given as

$\begin{split}\text{EFIC}(\mathcal{I})=N\ln\big{\lVert}\mathbf{\Pi}^{\perp}_{\mathcal{I}}\mathbf{y}\big{\rVert}^{2}_{2}+k\ln
N+\ln\big{|}\mathbf{A}^{T}_{\mathcal{I}}\mathbf{A}_{\mathcal{I}}\big{|}\\\
-(k+2)\ln\big{\lVert}\mathbf{\Pi}^{\perp}_{\mathcal{I}}\mathbf{y}\big{\rVert}^{2}_{2}+2ck\ln
p.\end{split}$ (24)

EFIC is consistent in both large-$N$ and high-SNR scenarios [13]. However,
EFIC suffers from a data scaling problem due to the inclusion of the data
dependent penalty term and as such the performance of EFIC is not invariant to
data scaling. See further in Section III-A.

##  III Proposed criterion: EBIC-Robust (EBIC${}_{\text{R}}$)

In this section, we present the necessary steps for deriving
EBIC${}_{\text{R}}$. EBIC${}_{\text{R}}$ can be seen as a natural extension of
BIC${}_{\text{R}}$ [24] for performing model selection in large-$p$ small-$N$
scenarios. Below, we provide a detailed derivation and establish the
connection to BIC${}_{\text{R}}$. A similar approach as in [23] is considered,
but here we perform normalization of $\hat{\mathbf{F}}_{\mathcal{I}}$ under
both large-$N$ and high-SNR assumption. It is possible to factorize the
$\ln\big{\lvert}\hat{\mathbf{F}}_{\mathcal{I}}\big{\rvert}$ term in (15) in
the following manner

$\displaystyle\ln{\big{|}\hat{\mathbf{F}}_{\mathcal{I}}\big{|}}=$
$\displaystyle\ln\left[\big{|}\mathbf{L}\big{|}\left\lvert\mathbf{L}^{-1/2}\hat{\mathbf{F}}_{\mathcal{I}}\mathbf{L}^{-1/2}\right\rvert\right]$
$\displaystyle=$
$\displaystyle\ln|\mathbf{L}|+\ln\underbrace{\Big{|}\mathbf{L}^{-1/2}\mathbf{\hat{F}}_{\mathcal{I}}\mathbf{L}^{-1/2}\Big{|}}_{\text{T}}.$
(25)

The goal here is to choose a suitable $\mathbf{L}$ matrix that normalizes the
sample FIM $\hat{\mathbf{F}}_{\mathcal{I}}$ such that the T term in (25) is
$\mathcal{O}(1)$, i.e., in this case T should be bounded as $N\to\infty$
and/or $\sigma^{2}\to 0$. To accomplish this objective, we choose the
following $\mathbf{L}^{-1/2}$ matrix

$\mathbf{L}^{-1/2}=\begin{bmatrix}\sqrt{\frac{1}{N}}\sqrt{\frac{\hat{\sigma}_{\mathcal{I}}^{2}}{\hat{\sigma}_{0}^{2}}}\mathbf{I}_{k}&\mathbf{0}\\\
\mathbf{0}&\sqrt{\frac{1}{N}}\frac{\hat{\sigma}_{\mathcal{I}}^{2}}{\hat{\sigma}_{0}^{2}}\end{bmatrix},$
(26)

where $\hat{\sigma}^{2}_{0}=\lVert\mathbf{y}\rVert^{2}_{2}/N$. The factor,
$\hat{\sigma}^{2}_{0}$, is used in $\mathbf{L}^{-1/2}$ in order to neutralize
the data scaling problem and is motivated by the fact that given (17), when
the SNR is a constant, we have

$\mathbb{E}[\hat{\sigma}^{2}_{0}]\to\text{const.}\quad\&\quad\text{Var}[\hat{\sigma}^{2}_{0}]\to
0$ (27)

as $N\to\infty$. Furthermore, from the considered generating model in (1),
when $N$ is fixed, (27) is also satisfied as $\sigma^{2}\to 0$ (see Appendix B
for details on $\hat{\sigma}^{2}_{0}$). Now using (12), (26) and the
assumptions in (17), (27) it is possible to show that

$\Big{|}\mathbf{L}^{-1/2}\mathbf{\hat{F}}_{\mathcal{I}}\mathbf{L}^{-1/2}\Big{|}=\begin{vmatrix}{\frac{1}{\hat{\sigma}_{0}^{2}}}{\frac{\mathbf{A}^{T}_{\mathcal{I}}\mathbf{A}_{\mathcal{I}}}{N}}&\mathbf{0}\\\
\mathbf{0}&\frac{1}{2\hat{\sigma}_{0}^{4}}\end{vmatrix}=\mathcal{O}(1),$ (28)

and therefore may be discarded without much effect on the criterion.
Furthermore, the $\ln\big{\lvert}\mathbf{L}\big{\rvert}$ term can be expanded
as follows

$\displaystyle\ln|\mathbf{L}|$
$\displaystyle=\ln\begin{vmatrix}N\left(\frac{\hat{\sigma}_{0}^{2}}{\hat{\sigma}^{2}_{\mathcal{I}}}\right)\mathbf{I}_{k}&\mathbf{0}\\\
\mathbf{0}&N\left(\frac{\hat{\sigma}^{2}_{0}}{\hat{\sigma}^{2}_{\mathcal{I}}}\right)^{2}\end{vmatrix}$
$\displaystyle=(k+1)\ln
N+(k+2)\ln\left(\frac{\hat{\sigma}_{0}^{2}}{\hat{\sigma}_{\mathcal{I}}^{2}}\right).$
(29)

Therefore, using (28) and (29) we can rewrite (25) as

$\ln{\big{|}\hat{\mathbf{F}}_{\mathcal{I}}\big{|}}=k\ln
N+(k+2)\ln\left(\frac{\hat{\sigma}_{0}^{2}}{\hat{\sigma}_{\mathcal{I}}^{2}}\right)+\mathcal{O}(1)+\ln
N.$ (30)

Next, for the model prior probability term
$-2\ln\Pr(\mathcal{H}_{\mathcal{I}})$ in (15), a similar proposition is taken
as in EBIC such that $\Pr(\mathcal{H}_{\mathcal{I}})\propto{p\choose
k}^{-\zeta}$, where $\zeta\geq 0$ is a tuning parameter. For large-$p$, we
follow a similar approach as in EFIC by employing the following approximation
$\ln{p\choose k}\approx k\ln p$ . This gives

$-2\ln\Pr(\mathcal{H}_{\mathcal{I}})=2\zeta k\ln p+\text{const}.$ (31)

Now, substituting (30), (31) in (15) and dropping the $\mathcal{O}(1)$, the
$\ln N$ term (independent of $k$), the constant and the
$p(\hat{\boldsymbol{\theta}}_{\mathcal{I}}|\mathcal{H}_{\mathcal{I}})$ term we
arrive at the EBIC${}_{\text{R}}$:

$\begin{split}\text{EBIC}_{\text{R}}(\mathcal{I})=&N\ln\hat{\sigma}^{2}_{\mathcal{I}}+k\ln\left(\frac{N}{2\pi}\right)\\\
&+(k+2)\ln\left(\frac{\hat{\sigma}_{0}^{2}}{\hat{\sigma}_{\mathcal{I}}^{2}}\right)+2k\zeta\ln
p.\end{split}$ (32)

The true model is estimated as

$\hat{\mathcal{S}}_{\text{EBIC}_{\text{R}}}=\underset{\mathcal{I}\in\mathcal{J}}{\arg\min}\big{\\{}\text{EBIC}_{\text{R}}(\mathcal{I})\big{\\}},$
(33)

where $\mathcal{J}$ denotes the set of candidate models.

It can be observed from (32) that the penalty of EBIC${}_{\text{R}}$ is a
function of the number of measurements $N$, the ratio
$(\hat{\sigma}_{0}^{2}/\hat{\sigma}_{\mathcal{I}}^{2})$ and the parameter
dimension $p$. Notice that the ratio
$(\hat{\sigma}_{0}^{2}/\hat{\sigma}_{\mathcal{I}}^{2})$ is always greater than
$1$ and independent of the scaling of $\mathbf{y}$. Furthermore, when
$\mathcal{S}\not\subset\mathcal{I}$, the ratio
$(\hat{\sigma}_{0}^{2}/\hat{\sigma}_{\mathcal{I}}^{2})\approx\mathcal{O}(1)$
and for $\mathcal{S}\subset\mathcal{I}$ we have
$(\hat{\sigma}_{0}^{2}/\hat{\sigma}_{\mathcal{I}}^{2})\approx\mathcal{O}(\text{SNR}+1)$.
Hence, the behaviour of the penalty can be summarized as follows: (i) For
fixed $p$ and SNR, as $N\to\infty$ the penalty grows as $\mathcal{O}(\ln N)$;
(ii) If $N$ and $p$ are constant, as SNR $\to\infty$, the penalty grows as
$\mathcal{O}\left(\ln(\text{SNR}+1)\right)$ for all
$\mathcal{I}\supset\mathcal{S}$; (iii) when SNR is a constant and given that
$p$ grows with $N$, then as $N\to\infty$ the penalty grows as $\mathcal{O}(\ln
N)+\mathcal{O}(\ln p)$.

###  III-A Scaling Robustness as Compared to EFIC

In this section, we elaborately discuss the data scaling problem. Ideally, any
model selection criterion should be invariant to data scaling, which means
that if $\mathbf{y}$ is scaled by any arbitrary constant $C>0$, the equivalent
penalty for each of the models $\mathcal{I}$ should not change. This property
is necessary because otherwise the behaviour of the model selection criterion
will be unreliable and may suffer from overfitting or underfitting issues when
the data is scaled. As mentioned before, the penalty of EFIC is not invariant
to data scaling. This can be observed from the following analysis. Let
$\Delta={card}(\mathcal{I})-{card}(\mathcal{S})$. Now, consider the difference
assuming $\mathcal{I}\neq\mathcal{S}$

$\displaystyle\text{EFIC}(\mathcal{I})-\text{EFIC}(\mathcal{S})$
$\displaystyle=\
(N-2)\ln\frac{\big{\lVert}\boldsymbol{\Pi}^{\perp}_{\mathcal{I}}\mathbf{y}\big{\rVert}^{2}_{2}}{\big{\lVert}\boldsymbol{\Pi}^{\perp}_{\mathcal{S}}\mathbf{y}\big{\rVert}^{2}_{2}}+\ln\frac{\big{|}\mathbf{A}^{T}_{\mathcal{I}}\mathbf{A}_{\mathcal{I}}\big{|}}{\big{|}\mathbf{A}^{T}_{\mathcal{S}}\mathbf{A}_{\mathcal{S}}\big{|}}-k\ln\big{\lVert}\boldsymbol{\Pi}^{\perp}_{\mathcal{I}}\mathbf{y}\big{\rVert}^{2}_{2}$
$\displaystyle\ \
+k_{0}\ln\big{\lVert}\boldsymbol{\Pi}^{\perp}_{\mathcal{S}}\mathbf{y}\big{\rVert}^{2}_{2}+\Delta\left(\ln
N+2c\ln p\right)=D_{\text{EFIC}}\text{ (say)}.$ (34)

Ideally, for correct model selection, $D_{\text{EFIC}}>0$ for all
$\mathcal{I}\neq\mathcal{S}$. Now, if we scale the data $\mathbf{y}$ by a
constant $C>0$, the data dependent term becomes
$\ln\lVert\boldsymbol{\Pi}^{\perp}_{\mathcal{I}}C\mathbf{y}\rVert^{2}_{2}=\ln
C^{2}+\ln\lVert\boldsymbol{\Pi}^{\perp}_{\mathcal{I}}\mathbf{y}\rVert^{2}_{2}$
and the difference becomes

$\text{EFIC}(\mathcal{I})-\text{EFIC}(\mathcal{S})=D_{\text{EFIC}}-\Delta\ln
C^{2}.$ (35)

It is evident that (34) and (35) are unequal and the difference after scaling
contains an additional term $-\Delta\ln C^{2}$. This implies that scaling the
data changes the EFIC score difference between any arbitrary model
$\mathcal{I}$ and the true model $\mathcal{S}$. Hence, depending on the $C$
value ($C<1$ or $C\geq 1$) and $\Delta>0$ or $\Delta<0$, the difference in
(35) may become negative leading to a false model selection. Thus, EFIC is not
invariant to data scaling. On the contrary, consider the difference for
EBIC${}_{\text{R}}$,

$\displaystyle\text{EBIC}_{\text{R}}(\mathcal{I})-\text{EBIC}_{\text{R}}(\mathcal{S})$
$\displaystyle=(N-2)\ln\left(\frac{\hat{\sigma}^{2}_{\mathcal{I}}}{\hat{\sigma}^{2}_{\mathcal{S}}}\right)-k\ln\hat{\sigma}^{2}_{\mathcal{I}}+k_{0}\ln\hat{\sigma}^{2}_{\mathcal{S}}+\Delta\ln\hat{\sigma}^{2}_{0}$
$\displaystyle\ \ +\Delta\left(\ln(N/2\pi)+2\zeta\ln
p\right)=D_{\text{EBIC}_{\text{R}}}\text{(say)}$ (36)

Now, scaling $\mathbf{y}$ by $C$, scales the noise variance estimates
$\hat{\sigma}^{2}_{\mathcal{I}}$, $\hat{\sigma}^{2}_{\mathcal{S}}$ and
$\hat{\sigma}^{2}_{0}$ by $C^{2}$, however, the difference remains the same,
i.e., $D_{\text{EBIC}_{\text{R}}}$. This is because in this case the
$-\Delta\ln C^{2}$ term is cancelled by $+\Delta\ln C^{2}$ generated by
$\Delta\ln\hat{\sigma}^{2}_{0}$. Hence, EBIC${}_{\text{R}}$ is invariant to
data scaling, which is a desired property of any model selection criterion.

##  IV Consistency of EBIC${}_{\text{R}}$

In this section, we provide the necessary proofs to show that
EBIC${}_{\text{R}}$ is a consistent criterion. Generally speaking, a model
selection criterion with $\hat{\mathcal{S}}$ as its estimate of the true model
$\mathcal{S}$ is consistent if it satisfies the following conditions [13]

$\lim_{\sigma^{2}\to
0}\Pr\\{\hat{\mathcal{S}}=\mathcal{S}\\}=1\quad\&\quad\lim_{N\to\infty}\Pr\\{\hat{\mathcal{S}}=\mathcal{S}\\}=1.$
(37)

Let us define the set of all overfitted models of dimension $k$ as
$\mathcal{I}_{o}^{k}=\left\\{\mathcal{I}:{card}(\mathcal{I})=k,\mathcal{S}\subset\mathcal{I}\right\\}$
and the set of all misfitted models of dimension $k$ as
$\mathcal{I}^{k}_{m}=\left\\{\mathcal{I}:{card}(\mathcal{I})=k,\mathcal{S}\not\subset\mathcal{I}\right\\}$.
Furthermore, let $\mathbb{O}$ denote the set of all $\mathcal{I}^{k}_{o}$ for
$k=k_{0}+1,\ldots,K$, and let $\mathbb{M}$ denote the set of all
$\mathcal{I}^{k}_{m}$ for $k=1,\ldots,K$, i.e.,

$\mathbb{O}=\bigcup_{k=k_{0}+1}^{K}\mathcal{I}^{k}_{o}\quad\text{and}\quad\mathbb{M}=\bigcup_{k=1}^{K}\mathcal{I}^{k}_{m}\
,$ (38)

where $K$ is some upper bound for $k_{0}$ and $k_{0}\leq K\ll N$. In practice,
EBIC${}_{\text{R}}$ picks the true model $\mathcal{S}$, if the following
conditions are satisfied:

$\displaystyle\mathcal{C}_{1}:\text{EBIC}_{\text{R}}(\mathcal{S})<\text{EBIC}_{\text{R}}(\mathcal{I})\quad\forall\
\mathcal{I}\in\mathbb{O}$ (39)
$\displaystyle\mathcal{C}_{2}:\text{EBIC}_{\text{R}}(\mathcal{S})<\text{EBIC}_{\text{R}}(\mathcal{I})\quad\forall\
\mathcal{I}\in\mathbb{M}.$ (40)

###  IV-A Asymptotic Identifiability of the Model

In general, the model is identifiable if no model of comparable size other
than the true submodel can predict the noise free response almost equally well
[12]. In the context of linear regression, this is equivalent to say
$\mathbf{y}=\mathbf{A}_{\mathcal{S}}\mathbf{x}_{\mathcal{S}}\neq\mathbf{A}_{\mathcal{I}}\mathbf{x}_{\mathcal{I}}$
for
$\big{\\{}\mathcal{I}:{card}(\mathcal{I})\leq{card}(\mathcal{S}),\mathcal{I}\neq\mathcal{S}\big{\\}}$.
The identifiability of the true model in the high-dimensional linear
regression setup is uniformly maintained if the minimal eigenvalue of all
restricted sub-matrices,
$\mathbf{A}^{T}_{\mathcal{I}}\mathbf{A}_{\mathcal{I}}$ for
$\\{\mathcal{I}:{card}(\mathcal{I})\leq 2K\\}$, is bounded away from zero
[13]. A sufficient assumption on the design matrix $\mathbf{A}$ to prove the
consistency of EBIC${}_{\text{R}}$ is the sparse Riesz condition [28]:

$\lim\limits_{N\to\infty}\left\\{N^{-1}\left(\mathbf{A}_{\mathcal{I}}^{T}\mathbf{A}_{\mathcal{I}}\right)\right\\}=\mathbf{M}_{\mathcal{I}},\quad\forall\
card(\mathcal{I})\leq 2K,$ (41)

where $\mathbf{M}_{\mathcal{I}}$ denotes a bounded positive definite matrix.

###  IV-B Consistency as $\sigma^{2}\to 0$ or $\text{SNR}\to\infty$ for fixed
$N$

In this subsection, we examine whether EBIC${}_{\text{R}}$ selects the true
model $\mathcal{S}$ as $\sigma^{2}$ goes vanishingly small (or equivalently
SNR$\to\infty$) under the assumption that $N$ is fixed. We formulate this into
a theorem as follows:

###### Theorem 1

Assume that $N$ and $p$ are fixed and the matrix $\mathbf{A}$ satisfies the
condition given by (41). If $K\geq k_{0}$, then
$\Pr\left\\{\text{EBIC}_{\text{R}}(\mathcal{S})<\text{EBIC}_{\text{R}}(\mathcal{I})\right\\}\to
1$ as $\sigma^{2}\to 0$ for all $\mathcal{I}\neq\mathcal{S}$ and
$card(\mathcal{I})=1,\ldots,K$.

Proof. The proof consists of two parts. In part (a) we show that the
probability of overfitting
$(\mathcal{S}\subset\hat{\mathcal{S}}_{\textrm{EBIC}_{\textrm{R}}})$ tends to
$0$ as $\sigma^{2}\to 0$, which in this case is equivalent to showing
$\lim_{\sigma^{2}\to 0}\Pr(\mathcal{C}_{1})=1$, cf. (39). In part (b) we show
that the probability of misfitting
$(\mathcal{S}\not\subset\hat{\mathcal{S}}_{\textrm{EBIC}_{\textrm{R}}})$ also
tends to $0$ as $\sigma^{2}\to 0$, which is equivalent to $\lim_{\sigma^{2}\to
0}\Pr(\mathcal{C}_{2})=1$, cf. (40).

$(a)$ Over-fitting case
$(\mathcal{S}\subset\hat{\mathcal{S}}_{\textrm{EBIC}_{\textrm{R}}})$: Consider
the set of overfitted subsets having cardinality $k$, which we have denoted as
$\mathcal{I}^{k}_{o}$. Let $\mathcal{I}_{j}$ denote the $j$th subset in the
set $\mathcal{I}^{k}_{o}$. The total number of subsets in
$\mathcal{I}^{k}_{o}$ is ${p-k_{0}}\choose{\Delta}$ where $\Delta=k-k_{0}$ .
For any overfitted subset $\mathcal{I}_{j}\in\mathcal{I}^{k}_{o}$, consider
the following inequality

$\text{EBIC}_{\text{R}}(\mathcal{S})<\text{EBIC}_{\text{R}}(\mathcal{I}_{j}),\quad\mathcal{I}_{j}\in\mathcal{I}^{k}_{o},$
(42)

where $j=1,\ldots,{{p-k_{0}}\choose{\Delta}}$. Using the relation $p=N^{d}$
and after some straightforward rearrangement of (42) we get

$\displaystyle(N-k_{0}-2)\ln\hat{\sigma}^{2}_{\mathcal{S}}$
$\displaystyle-(N-k-2)\ln\hat{\sigma}^{2}_{\mathcal{I}_{j}}$
$\displaystyle-\Delta(1+2\zeta d)$ $\displaystyle\ln
N-\Delta\ln\hat{\sigma}^{2}_{0}+\Delta\ln 2\pi<0.$ (43)

Let us define a random variable
$X_{\mathcal{I}_{j}}={\hat{\sigma}^{2}_{\mathcal{I}_{j}}}/{\sigma^{2}}$, then

$N\cdot X_{\mathcal{I}_{j}}\sim\chi^{2}_{N-k},\quad\forall\
\mathcal{I}_{j}\in\mathcal{I}^{k}_{o}.$ (44)

This implies that the variables $X_{\mathcal{I}_{j}}$ are independent of
$\sigma^{2}$. Now, we can express

$(N-k-2)\ln\hat{\sigma}^{2}_{\mathcal{I}_{j}}=\ln
X_{\mathcal{I}_{j}}^{N-k-2}+(N-k-2)\ln\sigma^{2},$ (45)

and similarly by defining
$X_{\mathcal{S}}=\hat{\sigma}^{2}_{\mathcal{S}}/\sigma^{2}$ we get

$(N-k_{0}-2)\ln\hat{\sigma}^{2}_{\mathcal{S}}=\ln
X_{\mathcal{S}}^{N-k_{0}-2}+(N-k_{0}-2)\ln\sigma^{2}.$ (46)

Using (45) and (46) in (43) and after exponentiation we get

$\displaystyle\left(\frac{X^{N-k_{0}-2}_{\mathcal{S}}}{X^{N-k-2}_{\mathcal{I}_{j}}}\right)\left(\frac{1}{N}\right)^{\Delta(1+2\zeta
d)}\left(\frac{2\pi}{\hat{\sigma}^{2}_{0}}\right)^{\Delta}<\left(\frac{1}{\sigma^{2}}\right)^{\Delta}.$
(47)

Let $E^{k}_{\mathcal{I}_{j}}$ denote the entire left hand-side and let
$\eta_{k}$ denote the right-hand side of the inequality in (47). Let
$\mathcal{I}^{*}\in\mathcal{I}^{k}_{o}$ denote the subset that produces the
maximum value of $E_{\mathcal{I}_{j}}^{k}$ among all such subsets
$\mathcal{I}_{j}\in\mathcal{I}^{k}_{o}$ . Then, let us denote

$E^{k}_{\mathcal{I}^{*}}=\underset{\mathcal{I}_{j}\in\mathcal{I}^{k}_{o}}{\max}\left\\{E^{k}_{\mathcal{I}_{j}}\right\\},\quad
j=1,2,\ldots,{{p-k_{0}}\choose{\Delta}}.$ (48)

The condition $\mathcal{C}_{1}$ in (39) is satisfied as $\sigma^{2}\to 0$
under the event $E^{k}_{\mathcal{I}^{*}}<\eta_{k},$ for all
$k=k_{0}+1,\ldots,K.$ Now, we can express the probability that
$E^{k}_{\mathcal{I}^{*}}<\eta_{k}$ as follows

$\displaystyle\Pr\left(E^{k}_{\mathcal{I}^{*}}<\eta_{k}\right)$
$\displaystyle=\Pr\left\\{\bigcap_{j=1}^{{p-k_{0}}\choose{\Delta}}\left(E^{k}_{\mathcal{I}_{j}}<\eta_{k}\right)\right\\}$
$\displaystyle=1-\Pr\left\\{\bigcup_{j=1}^{{p-k_{0}}\choose{\Delta}}\left(E^{k}_{\mathcal{I}_{j}}>\eta_{k}\right)\right\\}$
$\displaystyle\geq
1-{{p-k_{0}}\choose{\Delta}}\Pr\left(E^{k}_{\mathcal{I}_{j}}>\eta_{k}\right)$
$\displaystyle\implies\Pr\left(E^{k}_{\mathcal{I}^{*}}>\eta_{k}\right)$
$\displaystyle\leq{{p-k_{0}}\choose{\Delta}}\Pr\left(E^{k}_{\mathcal{I}_{j}}>\eta_{k}\right),$
(49)

where the inequality follows from the union bound. Now consider the following
probability $\Pr\left\\{E^{k}_{\mathcal{I}_{j}}>\eta_{k}\right\\}$ for any
arbitrary subset $\mathcal{I}_{j}\in\mathcal{I}^{k}_{o}$, which can be
expressed as

$\displaystyle\Pr\left\\{\left(\frac{X^{N-k_{0}-2}_{\mathcal{S}}}{X^{N-k-2}_{\mathcal{I}_{j}}}\right)\left(\frac{1}{N}\right)^{\Delta(1+2\zeta
d)}\left(\frac{2\pi}{\hat{\sigma}^{2}_{0}}\right)^{\Delta}>\left(\frac{1}{\sigma^{2}}\right)^{\Delta}\right\\}.$
(50)

Let $W={X^{N-k_{0}-2}_{\mathcal{S}}}/{X^{N-k-2}_{\mathcal{I}_{j}}}$. Notice
that the random variable $W$ is independent of the noise variance $\sigma^{2}$
and since $N$ is fixed $W$ is bounded as $\sigma^{2}\to 0$. Furthermore,
$\lim\limits_{\sigma^{2}\to 0}\hat{\sigma}^{2}_{0}=c$ (see Appendix B) and the
right-hand side of the inequality in (50) grows unbounded as $\sigma^{2}\to
0$. Thus, we have

$\displaystyle\lim_{\sigma^{2}\to
0}\Pr\left\\{E^{k}_{\mathcal{I}_{j}}>\eta_{k}\right\\}=0.$ (51)

Therefore, using (49) and the result in (51), we have

$\displaystyle\lim_{\sigma^{2}\to
0}\Pr\left(E^{k}_{\mathcal{I}^{*}}>\eta_{k}\right)=0,\qquad\forall\
k=k_{0}+1,\ldots,K.$ (52)

Finally, using the union bound, and the result in (52), we get

$\displaystyle\Pr\left\\{\mathcal{C}_{1}\right\\}=$
$\displaystyle\Pr\left\\{\bigcap_{k=k_{0}+1}^{K}E^{k}_{\mathcal{I}^{*}}<\eta_{k}\right\\}$
$\displaystyle\geq$ $\displaystyle
1-\sum_{k=k_{0}+1}^{K}\Pr\left\\{E^{k}_{\mathcal{I}^{*}}>\eta_{k}\right\\}\to
1,$ (53)

as $\sigma^{2}\to 0$.

$(b)$ Misfitting case
$(\mathcal{S}\not\subset\hat{\mathcal{S}}_{\textrm{EBIC}_{\textrm{R}}})$: Let
$\mathcal{I}_{j}$ be any arbitrary $j$th subset belonging to the set of
misfitted subsets of dimension $k$, i.e., $\mathcal{I}^{k}_{m}$. We consider
the following inequality

$\text{EBIC}_{\text{R}}(\mathcal{S})<\text{EBIC}_{\text{R}}(\mathcal{I}_{j}),\quad\mathcal{I}_{j}\in\mathcal{I}^{k}_{m},$
(54)

where $j=1,\ldots,t$. Here, $t$ denotes the total number of subsets in the set
$\mathcal{I}^{k}_{m}$ and $t={{p}\choose{k}}$ if $k<k_{0}$, otherwise
$t={{p}\choose{k}}-{{p-k_{0}}\choose{\Delta}}$ if $k\geq k_{0}$, where
$\Delta=k-k_{0}$. Denoting
$X_{\mathcal{S}}=\hat{\sigma}^{2}_{\mathcal{S}}/\sigma^{2}$, rearranging and
applying exponentiation we can express (54) as

$\displaystyle\left(\frac{X^{N-k_{0}-2}_{\mathcal{S}}}{(\hat{\sigma}^{2}_{\mathcal{I}_{j}})^{N-k-2}}\right)\left(\frac{1}{N}\right)^{\Delta(1+2\zeta
d)}\left(\frac{2\pi}{\hat{\sigma}^{2}_{0}}\right)^{\Delta}<\left(\frac{1}{\sigma^{2}}\right)^{N-k_{0}-2}.$
(55)

Similar to the overfitting case, let $E^{k}_{\mathcal{I}_{j}}$ denote the
entire left-hand side and $\eta$ the right-hand side of (55). Also, let
$E^{k}_{\mathcal{I}^{*}}=\underset{\mathcal{I}_{j}\in\mathcal{I}_{m}^{k}}{\max}\left\\{E^{k}_{\mathcal{I}_{j}}\right\\}$
for $j=1,\ldots,t$, where $\mathcal{I}^{*}$ is the subset that leads to the
maximum value of $E^{k}_{\mathcal{I}_{j}}$ among all such subsets of dimension
$k$. The condition $\mathcal{C}_{2}$ in (40) is satisfied as $\sigma^{2}\to 0$
under the event $E^{k}_{\mathcal{I}^{*}}<\eta,$ for all $k=1,\ldots,K.$ Now,
we can express the probability that $E^{k}_{\mathcal{I}^{*}}<\eta$ as

$\displaystyle\Pr\left(E^{k}_{\mathcal{I}^{*}}<\eta\right)=$
$\displaystyle\Pr\left\\{\bigcap_{j=1}^{t}\left(E^{k}_{\mathcal{I}_{j}}<\eta\right)\right\\}$
$\displaystyle\implies\Pr\left(E^{k}_{\mathcal{I}^{*}}>\eta\right)\leq$
$\displaystyle t\Pr\left(E^{k}_{\mathcal{I}_{j}}>\eta\right),$ (56)

where the inequality follows from the union bound. Now consider the following
probability for any arbitrary subset $\mathcal{I}_{j}\in\mathcal{I}^{k}_{m}$

$\begin{split}\Pr\left(E^{k}_{\mathcal{I}_{j}}>\eta\right)=\Pr\bigg{\\{}\left(\frac{X^{N-k_{0}-2}_{\mathcal{S}}}{(\hat{\sigma}^{2}_{\mathcal{I}_{j}})^{N-k-2}}\right)\left(\frac{1}{N}\right)^{\Delta(1+2\zeta
d)}\\\
\times\left(\frac{2\pi}{\hat{\sigma}^{2}_{0}}\right)^{\Delta}>\left(\frac{1}{\sigma^{2}}\right)^{N-k_{0}-2}\bigg{\\}}.\end{split}$
(57)

Here, $X^{N-k_{0}-2}_{\mathcal{S}}$ is independent of $\sigma^{2}$ and $N$ is
fixed, therefore $X^{N-k_{0}-2}_{\mathcal{S}}$ is bounded as $\sigma^{2}\to
0$. Also
$\hat{\sigma}^{2}_{\mathcal{I}_{j}}\to\lVert\boldsymbol{\Pi}^{\perp}_{\mathcal{I}_{j}}\mathbf{A}_{\mathcal{S}}\mathbf{x}_{\mathcal{S}}\rVert^{2}_{2}/N$
in probability as $\sigma^{2}\to 0$ and since we are in the misfitting
scenario, from Lemma 4 in Appendix D we have
$\lVert\boldsymbol{\Pi}^{\perp}_{\mathcal{I}_{j}}\mathbf{A}_{\mathcal{S}}\mathbf{x}_{\mathcal{S}}\rVert^{2}_{2}/N>0$.
Furthermore, $\lim\limits_{\sigma^{2}\to 0}\hat{\sigma}^{2}_{0}=\text{const.}$
(see Appendix B) and the right-hand side of the inequality in (57) grows
unbounded as $\sigma^{2}\to 0$. Hence,

$\displaystyle\lim_{\sigma^{2}\to
0}\Pr\left\\{E^{k}_{\mathcal{I}_{j}}>\eta\right\\}=0.$ (58)

Using (56) and the result in (58) we get

$\displaystyle\lim_{\sigma^{2}\to
0}\Pr\left\\{E^{k}_{\mathcal{I}^{*}}>\eta\right\\}=0,\quad\forall\
k=1,\ldots,K.$ (59)

Finally, using the union bound and the result in (59), we get

$\Pr\left\\{\mathcal{C}_{2}\right\\}\geq
1-\sum_{k=1}^{K}\Pr\left\\{E^{k}_{\mathcal{I}^{*}}>\eta\right\\}\to
1\quad\text{as}\quad\sigma^{2}\to 0.$ (60)

From (53) and (60) we can conclude that EBIC${}_{\text{R}}$ is consistent as
$\sigma^{2}\to 0$, which proves Theorem 1.

###  IV-C Consistency as $N\to\infty$ when $\sigma^{2}$ is fixed
$(0<\sigma^{2}<\infty)$

In this section, we prove the consistency of EBIC${}_{\text{R}}$ as the sample
size $N\to\infty$ given that $\sigma^{2}$ is fixed and under the setting
$p=N^{d}$ for some $d>0$. This is a common setting in the model selection
literature (see, e.g., [12, 13, 29]). This leads to the following theorem.

###### Theorem 2

Assume that $p=N^{d}$ for some constant $d>0$, the SNR is fixed and the matrix
$\mathbf{A}$ satisfies (41). If $K\geq k_{0}$, then
$\Pr\left\\{\text{EBIC}_{\text{R}}(\mathcal{S})<\text{EBIC}_{\text{R}}(\mathcal{I})\right\\}\to
1$ as $N\to\infty$ for all $\mathcal{I}\neq\mathcal{S}$ and
$card(\mathcal{I})=1,\ldots,K$ under the condition $\zeta>1-1/2d$.

Proof. As in the previous section, we have two parts of the proof. Part $(a)$
is the overfitting case where we show that $\Pr(\mathcal{C}_{1})\to 1$ as
$N\to\infty$ and part $(b)$ is the misfitting case where we show that
$\Pr(\mathcal{C}_{2})\to 1$ as $N\to\infty$.
$(a)$ Overfitting case
$(\mathcal{S}\subset\hat{\mathcal{S}}_{\text{EBIC}_{\text{R}}})$: Let
$\mathcal{I}_{j}\in\mathcal{I}^{k}_{o}$ be any overfitted subset of dimension
$k$. Consider the following inequality

$\text{EBIC}_{\text{R}}(\mathcal{I}_{j})>\text{EBIC}_{\text{R}}(\mathcal{S}),\quad\mathcal{I}_{j}\in\mathcal{I}^{k}_{o}.$
(61)

Denoting $\Delta=k-k_{0}$ and rearranging (61) we get

$\displaystyle{\left(N-k-2\right)}\ln\left(\frac{\hat{\sigma}^{2}_{\mathcal{I}_{j}}}{\hat{\sigma}^{2}_{\mathcal{S}}}\right)+\Delta(1+2\zeta
d)\ln N$
$\displaystyle+\Delta\ln\left(\frac{\hat{\sigma}^{2}_{0}}{\hat{\sigma}^{2}_{\mathcal{S}}}\right)-\Delta\ln
2\pi>0.$ (62)

Let $E^{k}_{\mathcal{I}_{j}}$ denote the entire left side of the inequality
(62) and $\mathcal{I}^{*}$ denote the subset that leads to the minimum value
of $E^{k}_{\mathcal{I}_{j}}$ among all such subsets of dimension $k$. Hence,

$E^{k}_{\mathcal{I}^{*}}=\underset{\mathcal{I}_{j}\in\mathcal{I}^{k}_{o}}{\min}\left\\{E^{k}_{\mathcal{I}_{j}}\right\\},\quad
j=1,2,\ldots,{{p-k_{0}}\choose{\Delta}}.$ (63)

The condition $\mathcal{C}_{1}$ in (39) is satisfied as $N\to\infty$ under the
event $E^{k}_{\mathcal{I}^{*}}>0,$ for all $k=k_{0}+1,\ldots,K$. Expanding the
ratio we have

$\displaystyle\ln\left(\frac{\hat{\sigma}^{2}_{\mathcal{I}_{j}}}{\hat{\sigma}^{2}_{\mathcal{S}}}\right)$
$\displaystyle=\ln\left(\frac{\mathbf{e}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{I}_{j}}\mathbf{e}}{\mathbf{e}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{S}}\mathbf{e}}\right)$
$\displaystyle=\ln\left[\frac{\mathbf{e}^{T}\left(\mathbf{I}-\boldsymbol{\Pi}_{\mathcal{I}_{j}}+\boldsymbol{\Pi}_{\mathcal{S}}-\boldsymbol{\Pi}_{\mathcal{S}}\right)\mathbf{e}}{\mathbf{e}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{S}}\mathbf{e}}\right]$
$\displaystyle=\ln\left(\frac{\mathbf{e}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{S}}\mathbf{e}-\mathbf{e}^{T}\boldsymbol{\Pi}_{{\mathcal{I}_{j}}\setminus\mathcal{S}}\mathbf{e}}{\mathbf{e}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{S}}\mathbf{e}}\right)$
$\displaystyle=\ln\left(1-\frac{\mathbf{e}^{T}\boldsymbol{\Pi}_{{\mathcal{I}_{j}}\setminus\mathcal{S}}}{\mathbf{e}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{S}}\mathbf{e}}\right),$
(64)

where
$\boldsymbol{\Pi}_{{\mathcal{I}_{j}}\setminus\mathcal{S}}=\boldsymbol{\Pi}_{\mathcal{I}_{j}}-\boldsymbol{\Pi}_{\mathcal{S}}$.
Now we can write

$\displaystyle\underset{1\leq j\leq
T}{\min}\Bigg{\\{}(N-k-2)\ln\left(\frac{\hat{\sigma}^{2}_{\mathcal{I}_{j}}}{\hat{\sigma}^{2}_{\mathcal{S}}}\right)\Bigg{\\}}=$
$\displaystyle\ \ (N-k-2)\ln\left[1-\frac{\underset{1\leq j\leq
T}{\max}\left\\{\left(\mathbf{e}^{T}\boldsymbol{\Pi}_{{\mathcal{I}_{j}}\setminus\mathcal{S}}\mathbf{e}\right)/\sigma^{2}\right\\}}{\left(\mathbf{e}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{S}}\mathbf{e}\right)/\sigma^{2}}\right],$
(65)

where $T={p-k_{0}\choose\Delta}$. Now the term,
$(\mathbf{e}^{T}\boldsymbol{\Pi}_{{\mathcal{I}_{j}}\setminus\mathcal{S}}\mathbf{e})/\sigma^{2}\sim\chi^{2}_{\Delta}$
(see Appendix C). Then from Lemma 2 in Appendix D we have the following upper
bound

$\underset{1\leq j\leq
T}{\max}\Big{\\{}(\mathbf{e}^{T}\boldsymbol{\Pi}_{{\mathcal{I}_{j}}\setminus\mathcal{S}}\mathbf{e})/\sigma^{2}\Big{\\}}\leq\Delta+2\sqrt{\Delta\psi\ln
T}+2\psi\ln T,$ (66)

with probability approaching one as $N\to\infty$ if $\psi>1$. Now, for
sufficiently large $p=N^{d}$ we can write $\ln
T=\ln{p-k_{0}\choose\Delta}\approx\Delta d\ln N$. This gives

$\displaystyle\small\underset{1\leq j\leq
T}{\max}\Big{\\{}(\mathbf{e}^{T}\boldsymbol{\Pi}_{{\mathcal{I}_{j}}\setminus\mathcal{S}}\mathbf{e})/\sigma^{2}\Big{\\}}\leq\Delta+2\Delta\sqrt{\psi
d\ln N}+2\psi\Delta d\ln N$ $\displaystyle\qquad\qquad\qquad=2\psi\Delta d\ln
N\left(1+\frac{1}{\sqrt{\psi d\ln N}}+\frac{1}{2\psi d\ln N}\right)$
$\displaystyle\qquad\qquad\qquad\approx 2\psi\Delta d\ln N,$ (67)

as $N$ grows large. Furthermore, the term in the denominator in (65),
$(\mathbf{e}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{S}}\mathbf{e})/\sigma^{2}\sim\chi^{2}_{N-k_{0}}$
and based on the law of large numbers tends to $N-k_{0}\approx N$. Therefore,
using (67) in (65) and $(N-k-2)\approx N$ under the large-$N$ approximation we
get

$\displaystyle\underset{1\leq j\leq
T}{\min}\Bigg{\\{}N\ln\left(\frac{\hat{\sigma}^{2}_{\mathcal{I}_{j}}}{\hat{\sigma}^{2}_{\mathcal{S}}}\right)\Bigg{\\}}$
$\displaystyle\geq N\ln\left(1-\frac{2\Delta\psi d\ln N}{N}\right)$
$\displaystyle\approx-2\Delta\psi d\ln N,$ (68)

where the last approximation follows by linearization of the logarithm for
small $2\Delta\psi d\ln N/N$ value. Thus, we can write

$\displaystyle E^{k}_{\mathcal{I}^{*}}$ $\displaystyle\geq-2\Delta\psi d\ln
N+\Delta(1+2\zeta d)\ln
N+\Delta\ln\left(\frac{\hat{\sigma}^{2}_{0}}{\hat{\sigma}^{2}_{\mathcal{S}}}\right)$
$\displaystyle\hskip 142.26378pt-\Delta\ln 2\pi$
$\displaystyle=\Delta\left(1+2\zeta d-2\psi d\right)\ln
N+\Delta\ln\left(\frac{\hat{\sigma}^{2}_{0}}{\hat{\sigma}^{2}_{\mathcal{S}}}\right)-\Delta\ln
2\pi.$ (69)

Since $\lim_{N\to\infty}\hat{\sigma}^{2}_{0}=\text{const.}>0$ (see Appendix B)
and $\lim_{N\to\infty}\hat{\sigma}^{2}_{\mathcal{S}}=\sigma^{2}$ (see Appendix
C), $E^{k}_{\mathcal{I}^{*}}\to\infty$ as $N\to\infty$ for all
$k=k_{0}+1,\ldots,K$ under the condition $1+2\zeta d-2\psi d>0$ for any
$\psi>1$. Hence, the lower bound on $\zeta$ becomes

$\boxed{\zeta>1-\frac{1}{2d}}.$ (70)

From the above analysis we can say that

$\lim_{N\to\infty}\Pr\left\\{E^{k}_{\mathcal{I}^{*}}<0\right\\}=0,\quad\forall\
k=k_{0}+1,\ldots,K.$ (71)

Finally, using the union bound and the result in (71) we can express the
probability of $\mathcal{C}_{1}$ (39) happening as

$\displaystyle\Pr\left\\{\mathcal{C}_{1}\right\\}=$
$\displaystyle\Pr\left\\{\bigcap_{k=k_{0}+1}^{K}E^{k}_{\mathcal{I}^{*}}>0\right\\}$
$\displaystyle\geq
1-\sum_{k=k_{0}+1}^{K}\Pr\left\\{E^{k}_{\mathcal{I}^{*}}<0\right\\}\to 1$ (72)

as $N\to\infty$.

$(b)$ Misfitting case
$(\mathcal{S}\not\subset\hat{\mathcal{S}}_{\text{EBIC}_{\text{R}}})$: Let
$\mathcal{I}_{j}\in\mathcal{I}^{k}_{m}$ be any misfitted subset of dimension
$k$. Consider the following inequality

$\text{EBIC}_{\text{R}}(\mathcal{I}_{j})>\text{EBIC}_{\text{R}}(\mathcal{S}),\quad\mathcal{I}_{j}\in\mathcal{I}^{k}_{m}.$
(73)

Denoting $\Delta=k-k_{0}$ and rearranging (73) we get

$\begin{split}(N-k-2)\ln\left(\frac{\hat{\sigma}^{2}_{\mathcal{I}_{j}}}{\hat{\sigma}^{2}_{\mathcal{S}}}\right)+(1+2\zeta
d)\Delta\ln N\\\
+\Delta\ln\left(\frac{\hat{\sigma}^{2}_{0}}{\hat{\sigma}^{2}_{\mathcal{S}}}\right)+\Delta\ln\left(\frac{1}{2\pi}\right)>0.\end{split}$
(74)

Let $E^{k}_{\mathcal{I}_{j}}$ denote the entire left hand side of the
inequality in (74) and $\mathcal{I}^{*}$ denote the subset that generates the
minimum value of $E^{k}_{\mathcal{I}_{j}}$ among all such subsets of dimension
$k$. Then we have

$E^{k}_{\mathcal{I}^{*}}=\underset{\mathcal{I}_{j}\in\mathcal{I}^{k}_{m}}{\min}\left\\{E^{k}_{\mathcal{I}_{j}}\right\\},\quad\quad
j=1,2,\ldots,T,$ (75)

where $T={{p}\choose{k}}$ if $k<k_{0}$ otherwise
$T={{p}\choose{k}}-{{p-k_{0}}\choose{\Delta}}$ if $k\geq k_{0}$. The condition
$\mathcal{C}_{2}$ in (40) is satisfied as $N\to\infty$ under the event
$E^{k}_{\mathcal{I}^{*}}>0,$ for all $k=1,\ldots,K.$ Now, let
$\mathbf{u}=\mathbb{E}[\mathbf{y}]=\mathbf{A}_{\mathcal{S}}\mathbf{x}_{\mathcal{S}}$.
Using this, the ratio
$\frac{\hat{\sigma}^{2}_{\mathcal{I}_{j}}}{\hat{\sigma}^{2}_{\mathcal{S}}}$
can be expanded as

$\displaystyle\frac{\hat{\sigma}^{2}_{\mathcal{I}_{j}}}{\hat{\sigma}^{2}_{\mathcal{S}}}$
$\displaystyle=\frac{\mathbf{y}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{I}_{j}}\mathbf{y}}{\mathbf{y}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{S}}\mathbf{y}}=\frac{(\mathbf{u}+\mathbf{e})^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{I}_{j}}(\mathbf{u}+\mathbf{e})}{\mathbf{e}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{S}}\mathbf{e}}$
$\displaystyle=\frac{\mathbf{u}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{I}_{j}}\mathbf{u}+2\sigma\sqrt{\mathbf{u}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{I}_{j}}\mathbf{u}}\cdot
Z_{j}+\mathbf{e}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{I}_{j}}\mathbf{e}}{\mathbf{e}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{S}}\mathbf{e}},$
(76)

where

$Z_{j}=\frac{\mathbf{u}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{I}_{j}}\mathbf{e}}{\sigma\sqrt{\mathbf{u}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{I}_{j}}\mathbf{u}}}\sim\mathcal{N}(0,1).$
(77)

Now

$\displaystyle\underset{1\leq j\leq
T}{\min}\big{\\{}\hat{\sigma}^{2}_{\mathcal{I}_{j}}/\hat{\sigma}^{2}_{\mathcal{S}}\big{\\}}=$
$\displaystyle\small\underset{1\leq j\leq
T}{\min}\bigg{\\{}\mathbf{u}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{I}_{j}}\mathbf{u}+2\sigma\sqrt{\mathbf{u}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{I}_{j}}\mathbf{u}}\cdot
Z_{j}+\mathbf{e}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{I}_{j}}\mathbf{e}\bigg{\\}}\bigg{/}\mathbf{e}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{S}}\mathbf{e}$
$\displaystyle\geq\bigg{[}\underset{1\leq j\leq
T}{\min}\big{\\{}\mathbf{u}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{I}_{j}}\mathbf{u}\big{\\}}+\sigma^{2}\underset{1\leq
j\leq
T}{\min}\big{\\{}\mathbf{e}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{I}_{j}}\mathbf{e}/\sigma^{2}\big{\\}}$
$\displaystyle\quad-2\sigma\sqrt{\underset{1\leq j\leq
T}{\max}\big{\\{}\mathbf{u}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{I}_{j}}\mathbf{u}\big{\\}}}\cdot\underset{1\leq
j\leq
T}{\max}\big{\\{}Z_{j}\big{\\}}\bigg{]}\bigg{/}\mathbf{e}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{S}}\mathbf{e}.$
(78)

In the misfitting scenario we have two cases: (i) $k<k_{0}$ (ii) $k\geq
k_{0}$. We consider case (i) in our further analysis, which also encapsulates
case (ii). For $k<k_{0}$ we have $\ln T=\ln\binom{p}{k}\approx kd\ln N$.
Therefore, using the result in Lemma 2 we have the following lower bound under
large-$N$ approximation

$\displaystyle\underset{1\leq j\leq
T}{\min}\bigg{\\{}{\mathbf{e}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{I}_{j}}\mathbf{e}/\sigma^{2}}\bigg{\\}}=$
$\displaystyle\ {\mathbf{e}^{T}\mathbf{e}/\sigma^{2}-\underset{1\leq j\leq
T}{\max}\bigg{\\{}\mathbf{e}^{T}\boldsymbol{\Pi}_{\mathcal{I}_{j}}\mathbf{e}/\sigma^{2}\bigg{\\}}}$
$\displaystyle\geq\ N-2\psi^{\prime}kd\ln N,$ (79)

where $\psi^{\prime}>1$ and $\mathbf{e}^{T}\mathbf{e}/\sigma^{2}\approx N$ for
large-$N$. Furthermore, from the result in Lemma 3 we have the following upper
bound

$\underset{1\leq j\leq T}{\max}\\{Z_{j}\\}\leq\sqrt{2\psi^{\prime}kd\ln N},$
(80)

where $\psi^{\prime}>1$. Now, let $C_{\text{min}}=\underset{1\leq j\leq
T}{\min}\big{\\{}\mathbf{u}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{I}_{j}}\mathbf{u}\big{\\}}$
and $C_{\text{max}}=\underset{1\leq j\leq
T}{\max}\big{\\{}\mathbf{u}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{I}_{j}}\mathbf{u}\big{\\}}$.
Also as $N\to\infty$ we can approximate $(N-k-2)\approx N$ and
$\mathbf{e}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{S}}\mathbf{e}\approx\sigma^{2}N$.
Using this, and the results in (79) and (80) we get

$\displaystyle\underset{1\leq j\leq
T}{\min}\left\\{N\ln\left(\frac{\hat{\sigma}^{2}_{\mathcal{I}_{j}}}{\hat{\sigma}^{2}_{\mathcal{S}}}\right)\right\\}=N\ln\left[\underset{1\leq
j\leq
T}{\min}\bigg{\\{}\frac{\hat{\sigma}^{2}_{\mathcal{I}_{j}}}{\hat{\sigma}^{2}_{\mathcal{S}}}\bigg{\\}}\right]$
$\displaystyle\qquad\geq
N\ln\bigg{[}\bigg{\\{}C_{\text{min}}-2\sigma\sqrt{C_{\text{max}}}\cdot\sqrt{2\psi^{\prime}kd\ln
N}$ $\displaystyle\hskip 71.13188pt+\sigma^{2}\left(N-2\psi^{\prime}kd\ln
N\right)\bigg{\\}}\bigg{/}\sigma^{2}N\bigg{]}.$ (81)

Now, observe that
$C_{\text{min}}=\mathbf{u}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{I}^{*}}\mathbf{u}=\mathbf{x}_{\mathcal{S}}^{T}\mathbf{A}^{T}_{\mathcal{S}}\boldsymbol{\Pi}^{\perp}_{\mathcal{I}^{*}}\mathbf{A}_{\mathcal{S}}\mathbf{x}_{\mathcal{S}}$.
Since, we are in the misfitting scenario, from Lemma 4, in Appendix D, we can
express $C_{\min}=Nb_{\min}$ where $b_{\min}=\mathcal{O}(1)>0$. Similarly,
$C_{\max}=Nb_{\max}$ where $b_{\max}=\mathcal{O}(1)>0$ and $0<b_{\min}\leq
b_{\max}$. Hence, we can rewrite (81) as

$\displaystyle\underset{1\leq j\leq
T}{\min}\left\\{N\ln\left(\frac{\hat{\sigma}^{2}_{\mathcal{I}_{j}}}{\hat{\sigma}^{2}_{\mathcal{S}}}\right)\right\\}\geq$
$\displaystyle
N\ln\left(1+\frac{b_{\text{min}}}{\sigma^{2}}-\frac{2\sqrt{b_{\text{max}}}}{\sigma}\sqrt{\frac{2\psi^{\prime}kd\ln
N}{N}}-\frac{2\psi^{\prime}kd\ln N}{N}\right)$ $\displaystyle\approx
N\ln\left(1+\frac{b_{\text{min}}}{\sigma^{2}}\right)$ (82)

as $N$ grows large. For $k<k_{0}$, we get $\Delta<0$, therefore, in this case
we have

$\begin{split}E^{k}_{\mathcal{I}^{*}}\geq
N\ln\left(1+\frac{b_{\text{min}}}{\sigma^{2}}\right)-|\Delta|(1+2\zeta d)\ln
N\\\
-|\Delta|\ln\left(\frac{\hat{\sigma}^{2}_{0}}{2\pi\hat{\sigma}^{2}_{\mathcal{S}}}\right)\to\infty\end{split}$
(83)

as $N\to\infty$ for all $k=1,\ldots,K$, since
$N\ln(1+b_{\text{min}}/\sigma^{2})$ is the dominating term as it tends to
infinity much faster than the $\ln N$ term and
$\lim_{N\to\infty}\hat{\sigma}^{2}_{0}=\text{const.}>0$ (see Appendix B) and
$\lim_{N\to\infty}\hat{\sigma}^{2}_{\mathcal{S}}=\sigma^{2}$ (see Appendix C).
From the above analysis we can say that

$\lim_{N\to\infty}\Pr\left\\{E^{k}_{\mathcal{I}^{*}}<0\right\\}=0,\quad\forall\
k=1,\ldots,K.$ (84)

Finally, using the union bound and the result in (84) we can express the
probability of $\mathcal{C}_{2}$ (40) happening as

$\displaystyle\Pr\left\\{\mathcal{C}_{2}\right\\}$
$\displaystyle=\Pr\left\\{\bigcap_{k=1}^{K}E^{k}_{\mathcal{I}^{*}}>0\right\\}$
$\displaystyle\geq
1-\sum_{k=1}^{K}\Pr\left\\{E^{k}_{\mathcal{I}^{*}}<0\right\\}\to
1\quad\text{as}\quad N\to\infty.$ (85)

From (72) and (85) we can conclude that EBIC${}_{\text{R}}$ is consistent as
$N\to\infty$, which proves Theorem 2.

###  IV-D Discussion on the Hyperparameter $\zeta$

If $\zeta$ is too large, it will lead to underfitting issues. This is evident
from (83) where a large value of $\zeta$ may force the overall sum to become
negative especially for smaller $N$ values. On the contrary, if $\zeta$ is too
small, it will lead to overfitting issues as the penalty may not be
sufficiently large to compensate for the overparameterization due to large
parameter space.

##  V Predictor selection algorithms

In the high-dimensional scenario, when $p$ is large, it is infeasible to
perform model selection in the conventional manner. For a design matrix with
parameter dimension $p$, the number of possible candidate models is $2^{p}-1$.
Hence, the candidate model space grows exponentially with $p$ and we cannot
afford to calculate model score for all possible models. Therefore, to perform
model selection, we combine a model selection criterion with a predictor
selection (support recovery) algorithm such as OMP or LASSO (least absolute
shrinkage and selection operator) [30]. The goal of predictor selection is to
pick a subset of important predictors from the entire set of $p$ predictors.
In this context, the most important predictors refer to the positions of the
nonzero elements of the input signal $\mathbf{x}$. Thus, predictor selection
reduces the cardinality of the candidate model space to some upper bound $K$
such that $k_{0}\leq K\ll N$ under the assumption of a sparse parameter
vector. This enables us to apply the model selection criterion on the smaller
set of candidate models to pick the best model. The OMP algorithm is shown in
Algorithm 1. To perform model selection, we combine OMP with
EBIC${}_{\text{R}}$ as shown in Algorithm 2.

Algorithm 1 OMP with $K$ iterations

Inputs: Design matrix $\mathbf{A}$, measurement vector $\mathbf{y}$.

Initialization: $\lVert\mathbf{a}_{j}\rVert_{2}=1\ \forall j$,
$\mathbf{r}^{0}=\mathbf{y}$, $\mathcal{S}^{0}_{\text{OMP}}=\emptyset$

for $i=1\text{\ to\ }K$ do

Find next column index:
$d^{i}=\underset{j}{\arg\max}\big{|}\mathbf{a}_{j}^{T}\mathbf{r}^{i-1}\big{|}$

Add current index:
$\mathcal{S}^{i}_{\textrm{OMP}}=\mathcal{S}^{i-1}_{\textrm{OMP}}\cup\\{d^{i}\\}$

Update residual:
$\mathbf{r}^{i}=\left(\mathbf{I}_{N}-\mathbf{\Pi}_{\mathcal{S}^{i}_{\textrm{OMP}}}\right)\mathbf{y}$

end for

Output: OMP generated index sequence $\mathcal{S}^{K}_{\text{OMP}}$

Algorithm 2 Model selection combining EBIC${}_{\text{R}}$ with OMP

Run OMP for $K$ iterations to obtain $\mathcal{S}^{K}_{\text{OMP}}$

for $k=1\text{\ to\ }K$ do

$\mathcal{I}=\mathcal{S}^{k}_{\text{OMP}}$

Compute EBIC${}_{\text{R}}(\mathcal{I})$

end for

Estimated true support:
$\hat{\mathcal{S}}_{\text{EBIC}_{\text{R}}}=\underset{\mathcal{I}}{\arg\min}\\{$EBIC${}_{\text{R}}(\mathcal{I})\\}$

LASSO is a shrinkage method for variable selection/estimation in linear
regression models developed by Tibshirani [30]. Given the linear model in (1),
the LASSO solution for $\mathbf{x}$ for a particular choice of the
regularization parameter $\lambda\geq 0$ is obtained as

$\hat{\mathbf{x}}_{\text{lasso}}(\lambda)=\underset{\mathbf{x}\in\mathbb{R}^{p}}{\min}\left\\{\frac{1}{2N}\lVert\mathbf{y}-\mathbf{A}\mathbf{x}\rVert^{2}_{2}+\lambda\lVert\mathbf{x}\rVert_{1}\right\\},$
(86)

where $\lVert\cdot\rVert_{1}$ denotes the $l_{1}$ norm. The parameter
$\lambda$ determines the level of sparsity. When $\lambda\to\infty$ the
objective function in (86) attains the minimum with
$\hat{\mathbf{x}}_{\text{lasso}}(\lambda)$ being a zero vector. As we
gradually lower the $\lambda$ value, the number of non-zero components in
$\hat{\mathbf{x}}_{\text{lasso}}(\lambda)$ starts increasing. Model selection
combining LASSO and EBIC${}_{\text{R}}$ can be performed as shown in Algorithm
3. Gradually decrease $\lambda$ from a high value so that the number of non-
zero components in $\hat{\mathbf{x}}_{\text{lasso}}(\lambda)$ gradually
increases. Therefore, for each decreasing unique value of $\lambda$ say
$\lambda_{i}$, we acquire a different solution
$\hat{\mathbf{x}}_{\text{lasso}}(\lambda_{i})$, with increasing support and
thus obtaining a sequence of candidate models with maximum cardinality $K$.
The value of EBIC${}_{\text{R}}$ is computed for each of the candidate models
and the model corresponding to the smallest EBIC${}_{\text{R}}$ score is
selected as the final model. A most useful method for solving LASSO in our
context is the (modified) least angle regression (LARS) algorithm [31], since
it also provides the required sequence of regularization parameters for which
the support changes.

Algorithm 3 Model selection combining EBIC${}_{\text{R}}$ with LASSO

Compute LASSO estimates
$\\{\hat{\mathbf{x}}_{\text{lasso}}(\lambda_{1}),\ldots,\hat{\mathbf{x}}_{\text{lasso}}(\lambda_{K_{\text{max}}})\\}$
where
$card(supp\left(\hat{\mathbf{x}}_{\text{lasso}}(\lambda_{K_{\text{max}}})\right))=K$

for $i=1\text{\ to \ }K_{\max}$ do

$\mathcal{I}=supp\left(\hat{\mathbf{x}}_{\text{lasso}}(\lambda_{i})\right)$

Compute EBIC${}_{\text{R}}(\mathcal{I})$

end for

Estimated true support:
$\hat{S}_{\text{EBIC}_{\text{R}}}=\underset{\mathcal{I}}{\arg\min}\
\\{$EBIC${}_{\text{R}}(\mathcal{I})\\}$

##  VI Simulation Results

In this section, we provide numerical simulation results to illustrate the
empirical performance of EBIC${}_{\text{R}}$. The performance of
EBIC${}_{\text{R}}$ is compared with the ‘oracle’, EBIC, EFIC and MBT.
However, the performance comparison with the RRT [18] method is dropped since
it behaves quite similar to MBT (see [19] for details). The ‘oracle’ criterion
assumes a priori knowledge of the true cardinality $k_{0}$. Thus, the model
selection performance of the ‘oracle’ provides the upper bound on the model
selection performance that can be achieved using a particular predictor
selection algorithm and for a given set of data settings. Additionally, we
also provide simulation results to highlight the drawbacks of classical
methods for model selection in high-dimensional linear regression models with
a sparse parameter vector.

###  VI-A General Simulation Setup

In the simulations, we consider the model $\mathbf{y=Ax+e}$, where the design
matrix $\mathbf{A}\in\mathbb{R}^{N\times p}$ is generated with independent
entries following normal distribution $\mathcal{N}(0,1)$. Since $\mathbf{x}$
is assumed to be sparse, we choose $k_{0}=5$. Furthermore, without loss of
generality, we assume that the true support is $\mathcal{S}=[1,2,3,4,5]$,
therefore, $\mathbf{x}_{\mathcal{S}}=[x_{1},x_{2},x_{3},x_{4},x_{5}]^{T}$ and
$\mathbf{A}_{\mathcal{S}}=[\mathbf{a}_{1},\mathbf{a}_{2},\mathbf{a}_{3},\mathbf{a}_{4},\mathbf{a}_{5}]$.
This implies that the elements of $\mathbf{x}$ follows $x_{k}\neq 0$ for
$k=1,\ldots,k_{0}$ and $x_{k}=0$ for $k>k_{0}$. The SNR in dB is SNR (dB) =
$10\log_{10}(\sigma_{s}^{2}/\sigma^{2})$, where $\sigma_{s}^{2}$ and
$\sigma^{2}$ denote signal and true noise power, respectively. The signal
power is computed as
$\sigma_{s}^{2}=||\mathbf{A}_{\mathcal{S}}\mathbf{x}_{\mathcal{S}}||^{2}_{2}/N$.
Based on $\sigma_{s}^{2}$ and the chosen SNR (dB), the noise power is set as
$\sigma^{2}=\sigma_{s}^{2}/10^{\textrm{SNR (dB)}/10}$. Using this
$\sigma^{2}$, the noise vector $\mathbf{e}$ is generated following
$\mathcal{N}(\mathbf{0},\sigma^{2}\mathbf{I}_{N})$. The probability of correct
model selection (PCMS) is estimated over $1000$ Monte Carlo trials. To
maintain randomness in the data, a new design matrix $\mathbf{A}$ is generated
at each Monte Carlo trial. OMP is used for predictor selection for its
simplicity and wider range of applicability.

###  VI-B Tuning Parameter Selection

Figure 1: PCMS vs $N$ with $\mathbf{x}_{\mathcal{S}}=[1,1,1,1,1]$, SNR = 5 dB,
$p=N^{d}$ and $d=1.1$.

An important step in model selection is the choice of the tuning parameter. As
mentioned earlier, too small or large values of the tuning parameter can cause
severe performance degradation in certain scenarios. Fig. 1 shows a
performance comparison of EBIC${}_{\text{R}}$ for four different values of
$\zeta$ (0.4, 0.6, 1, and 2). Here, we set $p=N^{d}$ where $d=1.1$. Hence,
from Theorem 2 we require $\zeta>1-1/2d=0.55$ to achieve consistency. From the
figure we see that for $\zeta=0.4$, the performance of EBIC${}_{\text{R}}$
degrades after a certain point with increasing $N$, which justifies the
theory. For all other $\zeta>0.55$, the performances improve with increasing
$N$. For $\zeta=0.6$, which is very close to the lower bound, the convergence
to correct selection probability one is slow and will require a very large
sample size. For, $\zeta=2$, the performance suffers (due to underfitting) in
the low $N$ regime, but do achieve perfect selection as $N$ increases. In this
case, $\zeta=1$ provides a much better overall performance for a broader range
of $N$. A similar trend as in EBIC${}_{\text{R}}$ is observed even in EBIC and
EFIC for different choices of $\gamma$ and $c$. Hence, to maintain fairness,
the following tuning parameter settings are considered for further analysis:
$\zeta=1$ (EBIC${}_{\text{R}}$), $c=1$ (EFIC) and $\gamma=1$ (EBIC). For MBT
[19], $\lim\limits_{N\to\infty}\text{PCMS}\to 1$ as $\beta\to 1$. Hence, we
choose $\beta=0.999$.

###  VI-C Model Selection with Classical Methods in High-Dimensional Setting

Figure 2: The PCMS versus $N$ for SNR = 30 dB with
$\mathbf{x}_{\mathcal{S}}=[5,4,3,2,1]$ and $p=N^{d}$ where $d=1.1$.

In this section, we present simulations results for model selection using
classical methods in high-dimensional linear regression models and compare
their performances with EBIC${}_{\text{R}}$. The purpose of these results is
to highlight the limitations of the classical methods in dealing with
large-$p$ small-$N$ scenarios. The classical methods used here are BIC [7],
$\widetilde{\text{BIC}}_{N,\text{SNR}}$[23], BIC${}_{\text{R}}$[24], gMDL [9],
and PAL [11].

In the simulation, we consider the true parameter vector to be
$\mathbf{x}_{\mathcal{S}}=[5,4,3,2,1]^{T}$. Fig. 2 presents the plot for PCMS
versus $N$ for SNR = 30 dB with $p=N^{d}$ where $d=1.1$. The figure shows that
EBIC${}_{\text{R}}$ ($\zeta=1$) clearly surpasses the classical methods with
huge differences in performance. In general, when $p$ is fixed and
$N\to\infty$, the classical methods are consistent [24]. However, when $p$ is
varying and grows exponentially with $N$, the consistency attribute does not
hold any longer, hence, we see the decreasing performance trend in Fig. 2.

Figure 3: The PCMS versus SNR (dB) for $N=100$, $p=500$ and
$\mathbf{x}_{\mathcal{S}}=[5,4,3,2,1]$.

(a) $\mathbf{x}_{\mathcal{S}}=[0.05,0.04,0.03,0.02,0.01]^{T}$

(b) $\mathbf{x}_{\mathcal{S}}=[50,40,30,20,10]^{T}$

Figure 4: The PCMS versus SNR (dB) for $N=55$ and $p=1000$.

Fig. 3 illustrates the PCMS versus SNR in dB for fixed $N=100$ and $p=500$.
This gives $d=\log(p)/\log(N)\approx 1.35$, hence, $\zeta>1-1/2d\approx 0.63$.
The first major observation from the figure is that EBIC${}_{\text{R}}$
($\zeta=1$) clearly outperforms all the classical methods by a huge margin.
Secondly, for the considered setting, the performances of BIC${}_{\text{R}}$
and gMDL are quite similar followed by
$\widetilde{\text{BIC}}_{N,\text{SNR}}$. The criteria BIC${}_{\text{R}}$, gMDL
and $\widetilde{\text{BIC}}_{N,\text{SNR}}$ do achieve convergence to
detection probability one but at the expense of very high values of SNR. The
performances of PAL and BIC are extremely poor in this case, even in the high-
SNR regions.

###  VI-D Model Selection with the Latest Methods in High-Dimensional Setting

In Section VI-C, we highlighted the drawbacks of classical methods in model
selection under the high-dimensional setting. We observed that the performance
of the classical methods collapses when $p$ grows exponentially with $N$ and
the consistency property breaks down. In this section, we present simulation
results for model selection comparing EBIC${}_{\text{R}}$ to the existing
state-of-the-art methods, designed to deal with the large-$p$ small-$N$
scenarios.

####  VI-D1 Model Selection versus SNR

To highlight the scale-invariant and consistent behaviour of
EBIC${}_{\text{R}}$, we consider two scenarios. In the first scenario, we
assume the true parameter vector to be
$\mathbf{x}_{\mathcal{S}}=[0.05,0.04,0.03,0.02,0.01]^{T}$ and in the second
scenario, we assume $\mathbf{x}_{\mathcal{S}}=[50,40,30,20,10]^{T}$. Note that
in the simulations we compute the noise variance $\sigma^{2}$ based on the
chosen SNR level and the current signal power value
$\sigma^{2}_{s}=\lVert\mathbf{A}_{\mathcal{S}}\mathbf{x}_{\mathcal{S}}\rVert^{2}_{2}\big{/}N$.
To simulate the probability of correct model selection versus SNR in a high-
dimensional setting we fixed $N=55$ and $p=1000$. This gives
$d=\log(p)/\log(N)\approx 1.724$, hence, $\zeta>1-1/2d\approx 0.71$.

Figure 5: The PCMS versus $N$ for SNR = $6$ dB, $p=1000$ and
$\mathbf{x}_{\mathcal{S}}=[50,40,30,20,10]^{T}$.

Fig. 4 shows the empirical PCMS versus SNR (dB). Fig. 4(a) and Fig. 4(b)
correspond to $\mathbf{x}_{\mathcal{S}}=[0.05,0.04,0.03,0.02,0.01]$ and
$\mathbf{x}_{\mathcal{S}}=[50,40,30,20,10]$, respectively. Both the figures
depict fixed $N$ increasing SNR scenario. Comparing the figures, the first
clear observation is that unlike the other criteria, the behaviour of EFIC is
not identical for the two different $\mathbf{x}_{\mathcal{S}}$ given that the
other parameters viz, $N$, $p$ and $k_{0}$ are constant and the performance is
evaluated for the same SNR range. This illustrates the scaling problem present
in EFIC that leads to either high underfitting or overfitting issues. This
behavior or EFIC can be explained as follows. The data dependent penalty term
(DDPT) of EFIC is
$\text{DDPT}=-(k+2)\ln\lVert\mathbf{\Pi}^{\perp}_{\mathcal{I}}\mathbf{y}\rVert^{2}_{2}$,
whose overall value depends on the value
$\lVert\mathbf{\Pi}^{\perp}_{\mathcal{I}}\mathbf{y}\rVert^{2}_{2}$, which in
turn is influenced by the signal and noise powers $\sigma^{2}_{s}$ and
$\sigma^{2}$, respectively. If
$\lVert\mathbf{\Pi}^{\perp}_{\mathcal{I}}\mathbf{y}\rVert_{2}^{2}\ll 1$, then
$\text{DDPT}\gg 0$, which may blow the overall penalty to a large value
leading to underfitting issues. This is most likely the case when
$\mathbf{x}_{\mathcal{S}}=[0.05,0.04,0.03,0.02,0.01]^{T}$ (Fig. 4(a)). On the
contrary if
$\lVert\mathbf{\Pi}^{\perp}_{\mathcal{I}}\mathbf{y}\rVert^{2}_{2}\gg 1$, then
$\text{DDPT}\ll 0$, thus lowering the overall penalty leading to overfitting
issues (when $\mathbf{x}_{\mathcal{S}}=[50,40,30,20,10]^{T}$, Fig. 4(b)). The
second major observation is that EBIC is inconsistent when SNR is high but $N$
is small and fixed. This behaviour of EBIC is already reported in [13]. In
general, EFIC, MBT (for $\beta\to 1$) and EBIC${}_{\text{R}}$ are consistent
for increasing SNR scenarios given that $N$ is fixed, but while
EBIC${}_{\text{R}}$ and MBT are invariant to data scaling EFIC is not.

Figure 6: The PCMS versus $N$ for SNR = $25$ dB, $p=N^{d}$, $d=1.3$ and
$\mathbf{x}_{\mathcal{S}}=[50,40,30,20,10]^{T}$.

####  VI-D2 Model Selection versus $N$

Fig. 5 illustrates the empirical PCMS versus $N$ for SNR = 6 dB, $p=1000$ and
$\mathbf{x}_{\mathcal{S}}=[50,40,30,20,10]^{T}$. It depicts a low-SNR
increasing $N$ scenario. It is clearly seen that compared to the other
criteria, EFIC suffers from the scaling issue and requires a large sample size
to achieve detection probability one. Among all the criteria, the performance
of EBIC and EBIC${}_{\text{R}}$ are closest to the oracle. Furthermore,
observe that the performance of EBIC${}_{\text{R}}$ and EBIC are more or less
alike for the current setting. This is primarily because the SNR is low (6dB)
hence the $(k+2)\ln(\hat{\sigma}^{2}_{0}/\hat{\sigma}^{2}_{\mathcal{I}})$ term
of EBIC${}_{\text{R}}$ behaves very close to a $\mathcal{O}(1)$ quantity for
$k\geq k_{0}$. Thus, for low SNR scenarios, the penalties of EBIC and
EBIC${}_{\text{R}}$ are similar and as such the behaviour of these two
criteria overlaps in this case. However, note that this is not true in the
high-SNR cases, which will be evident from the discussion following Fig. 6.

The plots shown in Fig. 4 and Fig. 5 represent fixed-$N$ increasing-SNR and
low-SNR increasing-$N$ scenarios, respectively. In Fig. 6, we present a high-
SNR increasing-$N$ case. Here, we consider a varying parameter space such that
$p=N^{d}$ where $d=1.3$. It is clearly observed that for high-SNR scenarios,
EBIC${}_{\text{R}}$ and MBT provide much faster convergence to oracle
behaviour as compared to EBIC that requires higher sample size to achieve
detection probability one. Furthermore, we also notice that EFIC suffers from
a higher false selection error and performs worse than EBIC in a certain
region of the sample size. This clearly shows the effects of scaling in the
behaviour of EFIC.

##  VII Conclusion

In this paper, we provided a new criterion, which is an extension of
BIC${}_{\text{R}}$, to handle model selection in sparse high-dimensional
linear regression models employing sparse methods for predictor selection. The
extended version is named as EBIC${}_{\text{R}}$, where the subscript ‘R’
stands for robust and it is a scale-invariant and consistent model selection
criterion. Additionally, we analytically examined the behaviour of
EBIC${}_{\text{R}}$ as $\sigma^{2}\to 0$ and as $N\to\infty$. In both cases,
it is shown that the probability of detecting the true model approaches one.
The paper further highlighted the data scaling issue present in EFIC, which is
a consistent criterion for both large sample size and high-SNR scenarios.
Extensive simulation results show that the performance of EBIC${}_{\text{R}}$
is either similar or superior to that of EBIC, EFIC and MBT.

## Appendix A

###### Lemma 1

Let $\mathbf{y}$ be a $N\times 1$ dimensional vector following
$\mathbf{y}\sim\mathcal{N}(\boldsymbol{\mu},\sigma^{2}\mathbf{I}_{N})$ and
$\mathbf{\Pi}$ be a $N\times N$ symmetric, idempotent matrix with
rank$(\mathbf{\Pi})=r$. Then the ratio
$\mathbf{y}^{T}\mathbf{\Pi}\mathbf{y}/\sigma^{2}$ has a non-central chi-square
distribution $\chi^{2}_{r}(\lambda)$ with $r$ degrees of freedom and non-
centrality parameter
$\lambda=\boldsymbol{\mu}^{T}\mathbf{\Pi}\boldsymbol{\mu}/\sigma^{2}$ (see,
e.g., Chapter 5 of [32]).

##  Appendix B statistical analysis of the factor $\hat{\sigma}^{2}_{0}$

From the generating model (1), the true data vector follows
$\mathbf{y}\sim\mathcal{N}\left(\mathbf{A}_{\mathcal{S}}\mathbf{x}_{\mathcal{S}},\sigma^{2}\mathbf{I}_{N}\right)$.
Consider the factor $\hat{\sigma}^{2}_{0}$, which is defined as

$\hat{\sigma}^{2}_{0}=\frac{\lVert\mathbf{y}\rVert^{2}_{2}}{N}=\left(\frac{\sigma^{2}}{N}\right)\frac{\mathbf{y}^{T}\mathbf{I}_{N}\mathbf{y}}{\sigma^{2}}.$
(87)

From Lemma 1 in Appendix A we have

$\frac{\mathbf{y}^{T}\mathbf{I}_{N}\mathbf{y}}{\sigma^{2}}\sim\chi^{2}_{N}(\lambda)\
\text{where}\
\lambda=\frac{\lVert\mathbf{A}_{\mathcal{S}}\mathbf{x}_{\mathcal{S}}\rVert^{2}_{2}}{\sigma^{2}}.$
(88)

This implies that
$\left(\frac{N}{\sigma^{2}}\right)\hat{\sigma}^{2}_{0}\sim\chi^{2}_{N}(\lambda).$
Therefore, the mean and variance of $\hat{\sigma}^{2}_{0}$ are:

$\begin{split}&\mathbb{E}[\hat{\sigma}^{2}_{0}]=\frac{\sigma^{2}}{N}(N+\lambda)=\sigma^{2}+\frac{\lVert\mathbf{A}_{\mathcal{S}}\mathbf{x}_{\mathcal{S}}\rVert^{2}_{2}}{N}\\\
&\text{Var}[\hat{\sigma}^{2}_{0}]=2\frac{\sigma^{4}}{N^{2}}(N+2\lambda)=2\frac{\sigma^{4}}{N}+4\frac{\sigma^{2}}{N^{2}}\lVert\mathbf{A}_{\mathcal{S}}\mathbf{x}_{\mathcal{S}}\rVert^{2}_{2}.\end{split}$
(89)

Hence, for a fixed $N$,

$\displaystyle\lim_{\sigma^{2}\to
0}\mathbb{E}[\hat{\sigma}^{2}_{0}]=\frac{\lVert\mathbf{A}_{\mathcal{S}}\mathbf{x}_{\mathcal{S}}\rVert^{2}_{2}}{N}\quad\&\quad\lim_{\sigma^{2}\to
0}\text{Var}[\hat{\sigma}^{2}_{0}]=0.$ (90)

Further, when SNR or $\sigma^{2}$ is fixed, using the assumption
$\lim_{N\to\infty}\left\\{\frac{\mathbf{A}_{\mathcal{S}}^{T}\mathbf{A}_{\mathcal{S}}}{N}\right\\}=\mathbf{M}_{\mathcal{S}}$
we get

$\displaystyle\lim_{N\to\infty}\mathbb{E}[\hat{\sigma}^{2}_{0}]=\sigma^{2}+\mathbf{x}_{\mathcal{S}}^{T}\mathbf{M}_{\mathcal{S}}\mathbf{x}_{\mathcal{S}}\quad\&\quad\lim_{N\to\infty}\text{Var}[\hat{\sigma}^{2}_{0}]=0,$
(91)

where $\mathbf{M}_{\mathcal{S}}$ is a bounded positive definite matrix and as
such
$\mathbf{x}_{\mathcal{S}}^{T}\mathbf{M}_{\mathcal{S}}\mathbf{x}_{\mathcal{S}}=\mathcal{O}(1)$
as $N$ grows large.

##  Appendix C statistical analysis of $\hat{\sigma}^{2}_{\mathcal{I}}$ when
$\mathcal{S}\subseteq\mathcal{I}$

The noise variance estimate under hypothesis $\mathcal{H}_{\mathcal{I}}$ can
be rewritten as

$\displaystyle\hat{\sigma}^{2}_{\mathcal{I}}=\left(\frac{\sigma^{2}}{N}\right)\frac{\mathbf{y}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{I}}\mathbf{y}}{\sigma^{2}}.$
(92)

The true model $\mathbf{u}=\mathbf{A}_{\mathcal{S}}\mathbf{x}_{\mathcal{S}}$
lies in a linear subspace spanned by the columns of
$\mathbf{A}_{\mathcal{S}}$. Consequently, for
$\mathcal{I}\supseteq\mathcal{S}$ we have
$\boldsymbol{\Pi}^{\perp}_{\mathcal{I}}\mathbf{u}=\mathbf{0}$. This implies
that
$\mathbf{y}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{I}}\mathbf{y}=\mathbf{e}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{I}}\mathbf{e}$.
Thus we have,

$\displaystyle\frac{\mathbf{y}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{I}}\mathbf{y}}{\sigma^{2}}=\frac{\mathbf{e}^{T}\boldsymbol{\Pi}^{\perp}_{\mathcal{I}}\mathbf{e}}{\sigma^{2}}\sim\chi^{2}_{N-k}\
\text{(Using Lemma \ref{lemma-1})},$ (93)

where $k=card(\mathcal{I})\geq k_{0}$. This implies that
$\left(\frac{N}{\sigma^{2}}\right)\hat{\sigma}^{2}_{\mathcal{I}}\sim\chi^{2}_{N-k}.$
Therefore, the mean and variance of $\hat{\sigma}^{2}_{\mathcal{I}}$ for
$\mathcal{I}\supseteq\mathcal{S}$ are:

$\mathbb{E}[\hat{\sigma}^{2}_{\mathcal{I}}]=\frac{\sigma^{2}}{N}(N-k)\quad\&\quad\text{Var}[\hat{\sigma}^{2}_{\mathcal{I}}]=2\frac{\sigma^{4}}{N^{2}}(N-k).$
(94)

Hence, when $\sigma^{2}$ is a constant,

$\lim_{N\to\infty}\mathbb{E}[\hat{\sigma}^{2}_{\mathcal{I}}]=\sigma^{2}\quad\&\quad\lim_{N\to\infty}\text{Var}[\hat{\sigma}^{2}_{\mathcal{I}}]=0.$
(95)

## Appendix D

###### Lemma 2

Let $Z_{\text{max}}=\underset{i}{\max}\big{\\{}Z_{i}\big{\\}}_{i=1}^{m}$ where
$Z_{1},Z_{2},\ldots,Z_{m}$ is a sequence of identically distributed random
variables (not necessarily independent) having a Chi-square distribution with
$k$ degrees of freedom where $k<m$. Then $Z_{\text{max}}\leq k+2\sqrt{k\psi\ln
m}+2\psi\ln m$ for some constant $\psi>1$ with probability approaching one as
$m\to\infty$.

Proof: From the union bound we have

$\Pr\left(Z_{\text{max}}\leq\eta\right)\geq 1-m\Pr\left(Z_{i}\geq\eta\right).$
(96)

Since $Z_{i}\sim\chi^{2}_{k}$, then from the Chi-square tail bound (Lemma 1 of
[33]) we have the following result

$\Pr\left(Z_{i}\geq k+2\sqrt{kt}+2t\right)\leq e^{-t}.$ (97)

Setting $t=\psi\ln m$ in (97) where $\psi>1$ we get

$\Pr\left(Z_{i}\geq k+2\sqrt{k\psi\ln m}+2\psi\ln m\right)\leq e^{-\psi\ln
m}=m^{-\psi}.$ (98)

Using (98) in (96) we get

$\Pr\left(Z_{\text{max}}\leq k+2\sqrt{k\psi\ln m}+2\psi\ln m\right)\geq
1-\frac{1}{m^{\psi-1}}.$ (99)

Therefore, $Z_{\text{max}}\leq k+2\sqrt{k\psi\ln m}+2\psi\ln m$ with
probability approaching one as $m\to\infty$ if $\psi>1$.

###### Lemma 3

Let $X_{\text{max}}=\underset{i}{\max}\big{\\{}X_{i}\big{\\}}_{i=1}^{m}$ where
$X_{1},X_{2},\ldots,X_{m}$ is a sequence of identically distributed random
variables (not necessarily independent) having a Gaussian distribution with
zero mean and variance one. Then $X_{\text{max}}\leq\sqrt{2\ln m}$ with
probability approaching one as $m\to\infty$.

Proof: From the union bound we have

$\Pr\left(X_{\text{max}}\leq\eta\right)\geq 1-m\Pr\left(X_{i}\geq\eta\right).$
(100)

Since $X_{i}\sim\mathcal{N}(0,1)$, from the Gaussian tail bound we have

$\Pr\left(X_{i}\geq\eta\right)\leq\frac{1}{\eta}\frac{e^{-\eta^{2}/2}}{\sqrt{2\pi}},$
(101)

for all $\eta>0$. Setting $\eta=\sqrt{2\ln m}$ in (101) we get

$\Pr\left(X_{i}\geq\sqrt{2\ln m}\right)\leq\frac{m^{-1}}{2\sqrt{\pi\ln m}}.$
(102)

Using (102) in (100) we get

$\Pr\left(X_{\text{max}}\leq\sqrt{2\ln m}\right)\geq 1-\frac{1}{2\sqrt{\pi\ln
m}}.$ (103)

Therefore, $X_{\text{max}}\leq\sqrt{2\ln m}$ with probability approaching one
as $m\to\infty$.

###### Lemma 4

For any arbitrary support $\mathcal{I}\in\mathcal{I}^{k}_{m}\in\mathbb{M}$,
under the asymptotic identifiability condition in (41) the following
inequality holds

$\left\lVert\mathbf{\Pi}^{\perp}_{\mathcal{I}}\mathbf{A}_{\mathcal{S}}\mathbf{x}_{\mathcal{S}}\right\rVert^{2}_{2}>0.$

Proof: Let $\mathcal{S}^{\prime}=\\{\mathcal{S}\setminus\mathcal{I}\\}$. The
true support $\mathcal{S}$ can be split into two disjoint subsets as
$\mathcal{S}=\\{\mathcal{S}\cap\mathcal{I}\\}\cup\\{\mathcal{S}\setminus\mathcal{I}\\}$.
Since
$\text{span}(\mathbf{A}_{\mathcal{S}\cap\mathcal{I}})\subset\text{span}(\mathbf{A}_{\mathcal{I}})$
we have

$\displaystyle\lVert\mathbf{\Pi}^{\perp}_{\mathcal{I}}\mathbf{A}_{\mathcal{S}}\mathbf{x}_{\mathcal{S}}\rVert^{2}_{2}=$
$\displaystyle\lVert\mathbf{\Pi}^{\perp}_{\mathcal{I}}\mathbf{A}_{\mathcal{S}^{\prime}}\mathbf{x}_{\mathcal{S}^{\prime}}\rVert^{2}_{2}$
$\displaystyle=$ $\displaystyle
N\mathbf{x}_{\mathcal{S}^{\prime}}^{T}\left(N^{-1}\mathbf{A}_{\mathcal{S}^{\prime}}^{T}\mathbf{\Pi}^{\perp}_{\mathcal{I}}\mathbf{A}_{\mathcal{S}^{\prime}}\right)\mathbf{x}_{\mathcal{S}^{\prime}}.$

Now, consider the matrix
$\mathbf{M}=\begin{bmatrix}\mathbf{A}_{\mathcal{S}^{\prime}}&\mathbf{A}_{\mathcal{I}}\end{bmatrix}$
where $card(\mathcal{S}^{\prime})\leq K$ and $card(\mathcal{I})\leq K$, such
that $card(\mathcal{S}^{\prime}\cup\mathcal{I})\leq 2K$. Under the assumption
(41)

$\displaystyle
N^{-1}\mathbf{M}^{T}\mathbf{M}=N^{-1}\begin{bmatrix}\mathbf{A}^{T}_{\mathcal{S}^{\prime}}\mathbf{A}_{\mathcal{S}^{\prime}}&\mathbf{A}^{T}_{\mathcal{S}^{\prime}}\mathbf{A}_{\mathcal{I}}\\\
\mathbf{A}^{T}_{\mathcal{I}}\mathbf{A}_{\mathcal{S}^{\prime}}&\mathbf{A}^{T}_{\mathcal{I}}\mathbf{A}_{\mathcal{I}}\end{bmatrix}$
(104)

is a bounded positive definite matrix. Then the Schur complement of the block
matrix $\mathbf{A}_{\mathcal{I}}^{T}\mathbf{A}_{\mathcal{I}}$ is

$\displaystyle
N^{-1}\left[\mathbf{A}^{T}_{\mathcal{S}^{\prime}}\mathbf{A}_{\mathcal{S}^{\prime}}-\mathbf{A}^{T}_{\mathcal{S}^{\prime}}\mathbf{A}_{\mathcal{I}}(\mathbf{A}^{T}_{\mathcal{I}}\mathbf{A}_{\mathcal{I}})^{-1}\mathbf{A}^{T}_{\mathcal{I}}\mathbf{A}_{\mathcal{S}^{\prime}}\right]$
$\displaystyle=$ $\displaystyle
N^{-1}\mathbf{A}^{T}_{\mathcal{S}^{\prime}}\mathbf{\Pi}^{\perp}_{\mathcal{I}}\mathbf{A}_{\mathcal{S}^{\prime}}$

is also positive definite and bounded as $N\to\infty$. Let
$\widetilde{\mathbf{M}}=N^{-1}\mathbf{A}_{\mathcal{S}^{\prime}}^{T}\mathbf{\Pi}^{\perp}_{\mathcal{I}}\mathbf{A}_{\mathcal{S}^{\prime}}$,
then,
$\mathbf{x}_{\mathcal{S}^{\prime}}^{T}\widetilde{\mathbf{M}}\mathbf{x}_{\mathcal{S}^{\prime}}=b\text{
(say) }=\mathcal{O}(1)>0$. Hence,
$\lVert\mathbf{\Pi}^{\perp}_{\mathcal{I}}\mathbf{A}_{\mathcal{S}}\mathbf{x}_{\mathcal{S}}\rVert^{2}_{2}=Nb>0$
for all $\mathcal{I}\in\mathcal{I}^{k}_{m}\in\mathbb{M}$.

## References

  * [1] J. Ding, V. Tarokh, and Y. Yang, “Model selection techniques: An overview,” _IEEE Signal Processing Magazine_ , vol. 35, no. 6, pp. 16–34, 2018.
  * [2] P. Stoica and Y. Selen, “Model-order selection: a review of information criterion rules,” _IEEE Signal Processing Magazine_ , vol. 21, no. 4, pp. 36–47, 2004.
  * [3] C. Rao, Y. Wu, S. Konishi, and R. Mukerjee, “On model selection,” _Lecture Notes-Monograph Series_ , pp. 1–64, 2001.
  * [4] D. Anderson and K. Burnham, “Model selection and multi-model inference,” _Second. NY: Springer-Verlag_ , vol. 63, p. 10, 2004.
  * [5] A. Chakrabarti and J. K. Ghosh, “AIC, BIC and recent advances in model selection,” _Philosophy of statistics_ , pp. 583–605, 2011.
  * [6] H. Akaike, “A new look at the statistical model identification,” _IEEE transactions on automatic control_ , vol. 19, no. 6, pp. 716–723, 1974.
  * [7] G. Schwarz _et al._ , “Estimating the dimension of a model,” _Annals of statistics_ , vol. 6, no. 2, pp. 461–464, 1978.
  * [8] J. Rissanen, “Modeling by shortest data description,” _Automatica_ , vol. 14, no. 5, pp. 465–471, 1978.
  * [9] M. H. Hansen and B. Yu, “Model selection and the principle of minimum description length,” _Journal of the American Statistical Association_ , vol. 96, no. 454, pp. 746–774, 2001.
  * [10] J. Rissanen, “MDL denoising,” _IEEE Transactions on Information Theory_ , vol. 46, no. 7, pp. 2537–2543, 2000.
  * [11] P. Stoica and P. Babu, “Model order estimation via penalizing adaptively the likelihood (PAL),” _Signal Processing_ , vol. 93, no. 11, pp. 2865–2871, 2013.
  * [12] J. Chen and Z. Chen, “Extended Bayesian information criteria for model selection with large model spaces,” _Biometrika_ , vol. 95, no. 3, pp. 759–771, 2008.
  * [13] A. Owrang and M. Jansson, “A model selection criterion for high-dimensional linear regression,” _IEEE Transactions on Signal Processing_ , vol. 66, no. 13, pp. 3436–3446, 2018.
  * [14] S. Kay, “Exponentially embedded families-new approaches to model order estimation,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 41, no. 1, pp. 333–345, 2005.
  * [15] H. Bozdogan, “Model selection and Akaike’s information criterion (AIC): The general theory and its analytical extensions,” _Psychometrika_ , vol. 52, no. 3, pp. 345–370, 1987.
  * [16] R. R. Picard and R. D. Cook, “Cross-validation of regression models,” _Journal of the American Statistical Association_ , vol. 79, no. 387, pp. 575–583, 1984.
  * [17] L. de Torrenté and T. Hastie, “Does cross-validation work when $p\gg n$?” 2012\.
  * [18] S. Kallummil and S. Kalyani, “Signal and noise statistics oblivious orthogonal matching pursuit,” in _International Conference on Machine Learning_. PMLR, 2018, pp. 2429–2438.
  * [19] P. B. Gohain and M. Jansson, “Relative cost based model selection for sparse high-dimensional linear regression models,” in _ICASSP IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2020, pp. 5515–5519.
  * [20] T. T. Cai and L. Wang, “Orthogonal matching pursuit for sparse signal recovery with noise,” _IEEE Transactions on Information theory_ , vol. 57, no. 7, pp. 4680–4688, 2011.
  * [21] P. B. Gohain and M. Jansson, “New improved criterion for model selection in sparse high-dimensional linear regression models,” in _ICASSP IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2022, pp. 5692–5696.
  * [22] S. M. Kay, _Fundamentals of statistical signal processing: estimation theory_. Prentice Hall PTR, 1993.
  * [23] P. Stoica and P. Babu, “On the proper forms of BIC for model order selection,” _IEEE Transactions on Signal Processing_ , vol. 60, no. 9, pp. 4956–4961, 2012.
  * [24] P. B. Gohain and M. Jansson, “Scale-invariant and consistent Bayesian information criterion for order selection in linear regression models,” _Signal Processing_ , p. 108499, 2022.
  * [25] D. F. Schmidt and E. Makalic, “The consistency of MDL for linear regression models with increasing signal-to-noise ratio,” _IEEE transactions on signal processing_ , vol. 60, no. 3, pp. 1508–1510, 2011.
  * [26] P. M. Djuric, “Asymptotic MAP criteria for model selection,” _IEEE Transactions on Signal Processing_ , vol. 46, no. 10, pp. 2726–2735, 1998.
  * [27] Q. Ding and S. Kay, “Inconsistency of the MDL: On the performance of model order selection criteria with increasing signal-to-noise ratio,” _IEEE Transactions on Signal Processing_ , vol. 59, no. 5, pp. 1959–1969, 2011.
  * [28] C.-H. Zhang and J. Huang, “The sparsity and bias of the lasso selection in high-dimensional linear regression,” _The Annals of Statistics_ , vol. 36, no. 4, pp. 1567–1594, 2008.
  * [29] N. Meinshausen and P. Bühlmann, “High-dimensional graphs and variable selection with the lasso,” _The annals of statistics_ , vol. 34, no. 3, pp. 1436–1462, 2006.
  * [30] R. Tibshirani, “Regression shrinkage and selection via the lasso,” _Journal of the Royal Statistical Society: Series B (Methodological)_ , vol. 58, no. 1, pp. 267–288, 1996.
  * [31] B. Efron, T. Hastie, I. Johnstone, and R. Tibshirani, “Least angle regression,” _The Annals of statistics_ , vol. 32, no. 2, pp. 407–499, 2004\.
  * [32] A. M. Mathai and S. B. Provost, _Quadratic forms in random variables: theory and applications_. Dekker, 1992\.
  * [33] B. Laurent and P. Massart, “Adaptive estimation of a quadratic functional by model selection,” _Annals of Statistics_ , pp. 1302–1338, 2000.


# Araguaia Medical Vision Lab at ISIC 2017 Skin Lesion Classification
Challenge

Rafael Teixeira Sousa  Universidade Federal
de Mato Grosso
Email<EMAIL_ADDRESS>Larissa Vasconcellos de Moraes  Universidade
Federal
de Mato Grosso

###### Abstract

This paper describes the participation of Araguaia Medical Vision Lab at the
International Skin Imaging Collaboration 2017 Skin Lesion Challenge. We
describe the use of deep convolutional neural networks in attempt to classify
images of Melanoma and Seborrheic Keratosis lesions. With use of finetuned
GoogleNet and AlexNet we attained results of $0.950$ and $0.846$ AUC on
Seborrheic Keratosis and Melanoma respectively.

###### Index Terms:

IEEEtran, journal, LaTeX, paper, template.

##  I Motivation and Overview

This paper is a small overview of Araguaia Medical Vision Lab (AMVL) at the
International Skin Imaging Collaboration (ISIC) 2017 challenge, more
specifically the skin lesion classification task. Our main objective is to
perform an automatic classification of skin lesions on two main tasks, the
Melanoma and Seborrheic Keratosis recognition, using the image dataset
available by ISIC, which was already diagnosed by specialists and used as
ground truth. The algorithm proposed a combination of deep convolutional
neural networks (CNN), GoogleNet[1] and AlexNet[2], fine-tuned [3] with
augmented skin lesion images. The next sessions will describe how was the
training and evaluation process.

##  II Image Pre-processing

The original dataset is composed by 2000 images, with 374 samples of Melanoma
and 254 samples of Seborrheic Keratosis, the other 1372 are defined as Nevus.
The images have different sizes from 1022 x 767 to 6748 x 4499.

The first step was split a validation set with around 20% of images from each
class to evaluate the neural network performance during the training stage.

All train dataset pass through a pre-process filter which applied random
shear, zoom, and vertical and horizontal shift and flip. This step was
necessary to increase the dataset size (around 5 times), make it less
unbalanced and improve the neural network accuracy.

##  III Network Training and Evaluation

To perform all the training and evaluation we used Nvidia Digits interface
running with Caffe [4].

###  III-A Seborrheic Keratosis Task

On Seborrheic classification task we got the best result resizing all training
images to 350 x 350 without losing proportions and training an AlexNet pre-
trained on ImageNet classification task dataset by the Berkeley Vision and
Learning Center (ref). The network was trained over 30 epochs with a
Stochastic Gradient Descent (SGD) using three stages, first 10 epochs with
learning rate of $0.001$, than 10 epochs with $0.0001$, and 10 more with
$0.00001$. The top accuracy on validation set was $98.2\%$, with $97.73\%$ on
Seborrheic class and $98.28\%$ on non-Seborrheic. On validation dataset from
ISIC the network we got $89.3\%$ accuracy and $0.950$ of Area under Roc Curve
(AUC), with sensibility of $0.786$ and specificity of $0.935$.

###  III-B Melanoma Task

On Melanoma task we got our best result doing an average between three
networks, GoogleNet 256, GoogleNet 224 and AlexNet 224, where all was pre-
trained on ImageNet classification dataset.

GoogleNet 256 was retrained over 72 epochs with 256 x 256 resized images with
random crop to 224 x 224, which is the network standard size, the learning was
with SGD and three learning rate steps $0.001$, $0.0001$, $0.00001$, the best
result was got with the last epoch.

GoogleNet 224 was retrained over 50 epochs with images resized to 224 x 224
without the random crop with SGD and learning steps of $0.005$, $0.0025$,
$0.00125$, $0.000625$ and $0.0003125$. Best result was at epoch 15

AlexNet 224 was retrained over 30 epochs with 224 x 224 images, SGD and three
learning rate steps $0.001$, $0.0001$, $0.00001$ over 30 epochs, best result
was at epoch 15.

The result was obtained by the arithmetic mean between the three networks
softmax output, making it more balanced between classes. On validation dataset
from ISIC the network we got $84.7\%$ accuracy and $0.846$ of AUC, with
sensibility of $0.633$ and specificity of $0.900$.

##  IV Conclusion

The CNNs got great results on both tasks, but seem to lack on sensibility,
probably because the unbalanced dataset with few positive images, making it
harder to generalize the visual features of lesions.

## References

  * [1] C. Szegedy, W. Liu, Y. Jia, P. Sermanet, S. E. Reed, D. Anguelov, D. Erhan, V. Vanhoucke, and A. Rabinovich, “Going deeper with convolutions,” _CoRR_ , vol. abs/1409.4842, 2014. [Online]. Available: http://arxiv.org/abs/1409.4842
  * [2] A. Krizhevsky, I. Sutskever, and G. E. Hinton, “Imagenet classification with deep convolutional neural networks,” in _Advances in neural information processing systems_ , 2012, pp. 1097–1105.
  * [3] C.-K. Shie, C.-H. Chuang, C.-N. Chou, M.-H. Wu, and E. Y. Chang, “Transfer representation learning for medical image analysis,” in _Engineering in Medicine and Biology Society (EMBC), 2015 37th Annual International Conference of the IEEE_. IEEE, 2015, pp. 711–714.
  * [4] Y. Jia, E. Shelhamer, J. Donahue, S. Karayev, J. Long, R. Girshick, S. Guadarrama, and T. Darrell, “Caffe: Convolutional architecture for fast feature embedding,” _arXiv preprint arXiv:1408.5093_ , 2014.