# Predictive Criteria for Prior Selection Using Shrinkage in Linear Models

Dean Dustin 11footnotemark: 1 Bertrand Clarke22footnotemark: 2 Jennifer
Clarke33footnotemark: 3 Department of Statistics, University of Nebraska-
Lincoln, 340 Hardin Hall North, PO Box 830963, Lincoln, NE, 68583-0963.
<EMAIL_ADDRESS>Same<EMAIL_ADDRESS>Same<EMAIL_ADDRESS>
###### Abstract

Choosing a shrinkage method can be done by selecting a penalty from a list of
pre-specified penalties or by constructing a penalty based on the data. If a
list of penalties for a class of linear models is given, we provide
comparisons based on sample size and number of non-zero parameters under a
predictive stability criterion based on data perturbation. These comparisons
provide recommendations for penalty selection in a variety of settings. If the
preference is to construct a penalty customized for a given problem, then we
propose a technique based on genetic algorithms, again using a predictive
criterion. We find that, in general, a custom penalty never performs worse
than any commonly used penalties but that there are cases the custom penalty
reduces to a recognizable penalty. Since penalty selection is mathematically
equivalent to prior selection, our method also constructs priors.

The techniques and recommendations we offer are intended for finite sample
cases. In this context, we argue that predictive stability under perturbation
is one of the few relevant properties that can be invoked when the true model
is not known. Nevertheless, we study variable inclusion in simulations and, as
part of our shrinkage selection strategy, we include oracle property
considerations. In particular, we see that the oracle property typically holds
for penalties that satisfy basic regularity conditions and therefore is not
restrictive enough to play a direct role in penalty selection. In addition,
our real data example also includes considerations merging from model mis-
specification.

###### keywords:

[class=MSC] [Primary ]62F15 [; secondary ]62C10, 60G25

###### keywords:

prediction , penalized methods , shrinkage , oracle property , prior selection
, genetic algorithm , evolutionary computation

††journal:

, ,

##  1 Shrinkage and Prediction

In the context of linear models, inference problems in which the number of
parameters $p$ is bigger than the sample size $n$, i.e., with $p>n$, are ill-
posed and require some form of regularization to be solved. That is, some
extra information must be added to ensure the existence of a reasonable
solution. The earliest form of this is called Tikhonov regularization and was
used initially for matrix inversion. In Statistics, ridge regression is
probably the first occurrence of Tikhonov regularization, see [1]. By the
early 1990’s, $L^{2}$ regularization was common in a neural networks context,
see [2], not just to ensure that a solution existed but also to reduce
variance. An important step forward was replacing the $L^{2}$ penalty with an
$L^{1}$ penalty, see [3]. This shrinkage method is called the least absolute
shrinkage and selection operator (LASSO). It provides a form of regularization
that does variable selection as well as variance reduction while ensuring
solutions exist. Another insight was the concept of an ‘oracle property’ (OP)
first proved for a penalty called the smoothly clipped absolute deviation
(SCAD) in the context of linear models; see [4]. The OP meant that
asymptotically, as $n\rightarrow\infty$, the parameter estimates from the SCAD
penalty behaved asymptotically as if the correct $p$ were known, i.e., the
parameter estimates either went to zero or were consistent, asymptotically
normal, and efficient according to whether the variables they represented were
not or were in the true model.

Over the last two decades, numerous shrinkage methods have been proposed and
studied as individual methods, representing individual priors and penalties,
for the purposes of parameter inference. Here, we want to choose a penalty or
prior from a class of penalties or priors. One way to do this is to fix a list
of shrinkage methods, find a basis for comparing them, and choose the best.
Alternatively, we can take a ‘build-your-own’ approach and construct a
shrinkage method from part of the data and use it on the rest of the data. In
this way we can choose the optimal shrinage method adaptively. Here, we
provide both a comparison of ‘off the shelf’ methods and a build-your-own
technique (based on genetic algorithms) under a predictive optimality
criterion. The build-your-own shrinkage methods never perform worse that the
off-the-shelf methods, but may return an off-the shelf method as optimal. We
give an example of each case.

More formally, write a linear model (LM) of the form $Y=X\beta+\epsilon$ where
$Y=(Y_{1},\ldots,Y_{n})^{T}$, $X$ is an $n\times p$ matrix of design values,
and $\epsilon=(\epsilon_{1},\ldots,\epsilon_{n})^{T}$ is random error
$\epsilon_{i}\sim N(0,\sigma^{2})$ for some $\sigma>0$ and assume we have a
data set ${\cal{D}}=\\{(y_{i},x_{i})\mid i=1,\ldots,n\\}$. A nonadaptive
shrinkage method gives parameter estimates

$\displaystyle\hat{\beta}=\arg\min_{\lambda,\beta}\sum_{i=1}^{n}L_{1}(y_{i}-x_{i}^{T}\beta)+\lambda\sum_{j=1}^{p}L_{2}(\beta_{j})$
(1.1)

where $L_{1}$ and $L_{2}$ are loss functions, $x_{i}$ is the $i$-th row of
$X$, and $\lambda$ is the decay parameter. The term shrinkage arises from the
fact that as $\lambda\rightarrow\infty$, each $\beta_{j}\rightarrow 0$.
Sometimes the first term on the right hand side of (1.1)is replaced by a log-
likelihood; in this case and if the likelihood is normal, $L_{1}$ corresponds
to squared error. If $L_{2}$ is also squared error we get the optimization
that gives ridge regression.

By contrast an adaptive shrinkage method typically gives estimates of the form

$\displaystyle\hat{\beta}=\arg\min_{\lambda,w^{p},\beta}\sum_{i=1}^{n}L_{1}(y_{i}-x_{i}^{T}\beta)+\lambda\sum_{j=1}^{p}w_{j}L_{2}(\beta_{j})$
(1.2)

where $w^{p}=(w_{1},\ldots,w_{p})^{T}$ and the $w_{j}$’s are weights on the
individual $\beta_{j}$’s. Again, the first term on the right hand side of
(1.2) may be replaced by a log-likelihood. Also, the dependence on $w_{j}$ in
the second term may be more complicated; we have represented the adaptivity of
the constraint to the data as multiplicative in the penalty term for the sake
of convenience. We regard the SCAD and minimax concave (MCV) penalties ([5])
as adaptive because their penalty terms or constraints are data driven even
though they only introduce one extra non-multiplicative parameter (and have
the OP). This is not simply parameter counting. Note that the elastic net
(EN), [6], introduces two parameters and is not adaptive while the adaptive EN
(AEN) is adaptive and has the OP; see [7].

As written, expression (1.2) introduces $p$ new parameters, the $w_{j}$’s.
When $n$ is not large, it is unclear how to estimate the $w_{j}$’s. Two main
techniques have been proposed. One, due to [8] (see also [9]), is to set
$\hat{w}_{j}=1/|\hat{\beta}_{j}|^{\gamma}$ where $\hat{\beta}_{j}$ is any
$\sqrt{n}$-consistent estimator of $\beta_{j}$ (e.g., from SCAD that only adds
one extra parameter or from linear models when $n$ is large enough) and choose
$\gamma$ by a cross-validation criterion. Another method, due to [10], sets
$\hat{w}_{j}=SE(\beta_{j,OLS})^{\gamma}/|\hat{\beta}_{j,OLS}|^{\gamma}$, using
the ordinary least squares (OLS) estimator. This requires $n>p$ unless some
auxiliary technique is used to find the SE. [10] argued that their method
works well in cases where there is high collinearity. In practice, for both
methods, $\gamma=1$ is used to avoid extra computation. Here, we have
exclusively used the Zou method since it does not require $n>p$ and has a nice
interpretation: As $\hat{\beta}_{j}\rightarrow 0$,
$\hat{w}_{j}\rightarrow\infty$ forcing $\beta_{j}=0$ in (1.2).

Formally, the OP is the following. Write $\beta=(\beta_{1},\beta_{2})$ where
$\beta_{2}$ represents the zero components of $\beta$. Under various
conditions, $\sqrt{n}$ consistent local minimizers
$\hat{\beta}=(\hat{\beta}_{1},\hat{\beta}_{2})$ from certain shrinkage
criteria (such as SCAD) satisfy

$\displaystyle\hat{\beta}_{2}\rightarrow
0\quad\hbox{and}\quad\sqrt{n}I(\beta_{1},0)(\hat{\beta}_{1}-\beta_{1})\rightarrow
N(0,1).$

Roughly, shrinkage methods segregate into those that are nonadaptive, i.e.,
introduce exactly one ‘decay’ parameter and usually do not satisfy the OP, and
those that introduce two or more decay parameters and often satisfy the OP. We
will see that, contrary to initial impressions, the oracle property is not
rare.

We note that (1.1) corresponds to a joint density $\rho$ on the data and
$\beta$. Indeed, exponentiating gives that $\hat{\beta}$ corresponds to the
mode of the posterior

$\displaystyle\rho(\beta\mid Y^{n})\propto
e^{-\lambda\sum_{j=1}^{p}L_{2}(\beta_{j})}e^{-\sum_{i=1}^{n}L_{1}(y_{i}-x_{i}^{T}\beta)}$

in which $L_{2}$ defines a prior on $\beta$ with hyperparameter $\lambda$. A
similar manipulation can be applied to (1.2). This means that penalty
selection is mathematically equivalent to prior selection. However, the
equivalence is only mathematical because the class of reasonable penalties is
a proper subset of the class of reasonable priors. In particular, most
reasonable penalties are convex – our method in Sec. 5.2 assumes convecity,
for instance – but priors do not have to be log-convex.

Even though shrinkage methods were originally introduced as a way to solve the
$n<p$ problem, the OP uses $n\rightarrow\infty$. This is partially ameliorated
by results that give analogs of the OP where $p$ increases much faster than
$n$ does; see [11]. However, the cost of amerlioration is often artificial
conditions on the parameter space and/or design matrix. Moreover, many of the
original examples given to verify that shrinkage methods were effective
actually had $n>p$ and took $p=8$, see [4], [8], and [9]. To the best of our
knowledge no systematic comparison of shrinkage methods for different relative
sizes of $n$ and $p$ has been done. This paper seeks to fill this gap and
provide recommendations for when different methods work well.

Our comparison is in terms of predictive stability and accuracy of model
selection. That is, after finding $\hat{\beta}$ for a given shrinkage method
we define the predictor

$\displaystyle\hat{Y}(x)=x^{T}\hat{\beta}$ (1.3)

for $Y(x)$ at some new value $x$. Then we evaluate how well $\hat{Y}$ predicts
when the data are perturbed. We perturb the data using the technique of [12].
The idea is to add $N(0,\tau^{2})$ noise to the $y_{i}$’s in ${\cal{D}}$ and
then use part of the data to form a predictor and the rest of the data to
evaluate the predictor. We do this by generating ‘instability’ curves,
basically $L^{2}$ errors as a function of $\tau$. A good predictor will have
instability curves with low values, that are smooth, and increase slowly with
$\tau$.

We also use more conventional accuracy measures for variable selection. We
argue that pairing the two provides an assessment that captures analogues of
both variance (from the instability curves) and bias (from the accuracy
measures). We regard instability as more important because even if a variable
is incorrectly included its contribution may be small if its coefficient is
near zero and if it is incorrectly excluded the bias should show up in the
instability curve.

A first contribution of this paper is the use of predictive stability to
compare ten shrinkage methods for four values of $n$ when $p=100$. Very
roughly, for $n<<p$ we find that ridge regression (RR) works comparatively
better but that, as sparsity increased, EN was often best. As $n$ increased
but remained less than $p$, EN improved relative to RR and was often the best.
When $n>p$, LM’s were usually best with little or no perturbation of the data
but performance rapidly deteriorated as perturbations or sparsity increased.
Also, a version of SCAD often performed best. For penalty selection in
practice, however, these conclusions must be tempered by the more detailed
analyses given in Sec. 3 where we include considerations from variable
selection and the intuition developed in Sec. 4 to address model mis-
specification

A second contribution of this paper is to observe that the OP is actually
relatively common; see Sec. 2. So, we needn’t limit ourselves to the specific
shrinkage methods that have been studied. When $n>p$, we search a class of
shrinkage methods that have the OP; when $n<p$ we search a class of shrinkage
methods that do not necessarily have the OP; see Subsec. 5.2. In either case,
we generate priors (or penalties) that are optimal for a given data set, again
in a predictive instability sense. In practice, this means we can use a
technique such as genetic algorithms (GAs) to optimize a fitness function
corresponding to good prediction to find optimal priors. We find that GA
optimization for use in prior selection benefits from choosing point estimates
for the $\beta_{j}$’s when $p<n$. These estimates can always be assumed to
exist and our theory from Sec. 2 carries over to this setting. That is, we
‘shrink’ to estimated values in the James-Stein sense to get shrinkage
estimates in the penalized likelihood sense.

Thus, a third contribution of this paper is the use of GAs and part of the
data to form a prior that we then apply predictively. As seen in Sec. 5, using
both GAs and shrinkage (when we can) together gives better predictive
performance. By optimizing over the prior in this way we are effectively
optimizing over the penalty and hence choosing the optimal shrinkage method.
This optimum may or may not coincide with an established shrinkage method but
will still have the OP when $n$ is large enough. Thus, we have chosen our
shrinkage method to be predictively optimal for our data. In fact, we are de
facto using GA’s to obtain an approxiamtion to the posterior based on the
first part of the data for use with the second part of the data. We verify in
examples that this is predictively better than simply using all the data to
make predictions.

The structure of this paper is as follows. Sec. 2 gives general conditions for
the OP to hold for a wide range of adaptive penalties. We present our
comparisons of existing shrinkage methods in Sec. 3 and provide general
recommendations on which to use in various settings. We compare our
recommendations to the results for a benchmark data in Sec. 4. In Sec. 5, we
present our GA optimization verifying that the theory in Sec. 2 holds and the
result of the GA is optimal given the data. We summarize our overall findings
and intuition in Sec. 6.

##  2 Theoretical Results

Our proofs are motivated by the techniques in [4] and [9]. We assume an
adaptive setting and allow different penalty functions on different
parameters. This is different from [11] who did not treat different penalty
functions on different parameters or adaptivity. Their main result concerned
the asymptotic equivalence of penalized methods and permitted $p$ to increase
with $n$. Like other papers that allow $p$ to increase with $n$, some of the
conditions appear artificial. For example, aside from truncations of the
parameter space, one must assume that there is a sequence of explanatory
variables such that if one of the early variables is correct and accidentally
not included it can be reconstructed from later explanatory variables in the
sequence, thereby sacrificing identifiability. Our result, like many others,
assumes either fixed $p$ or $p$ increasing so slowly with $n$ that the
required convergences hold.

###  2.1 General Penalized Log Likelihood

Write the linear model as

$Y_{i}=x_{i}^{\prime}\beta+\epsilon_{i},\hskip 18.06749pti=1,\ldots,n,$ (2.1)

where $x_{i}=(x_{i1},\ldots,x_{ip})^{\prime}$ is the p-dimensional covariate,
$\beta=(\beta_{i1},\ldots,\beta_{ip})^{\prime}$ is the vector of associated
regression coefficients, and $\epsilon_{i}$ are IID random errors with median
0. Assume that $\beta_{j}\neq 0$ for $j\leq p_{0}$ and $\beta_{j}=0$ for
$j>p_{0}$ for some $p_{0}\geq 0$. Here, we regard the $x_{i}^{\prime}s$ as
deterministic. When needed we write $\beta\in\Omega\subset\mathbb{R}^{p}$
where $\Omega$ is open and $\Omega=\bar{\Omega}^{0}$.

Let $(x_{i},Y_{i})$ for $i=1,\ldots,n$ each have density
$\rho(Y_{i}|x_{i},\beta)$ (with respect to a fixed dominating measure) such
that the six regularity conditions stated below are satisfied. Let
$L(\beta|x^{n})$ be the log-likelihood function of the observations
$(x_{1},Y_{1}),\ldots,(x_{n},Y_{n})$ and denote the penalized log-likelihood
objective function as

$Q(\beta)=L(\beta|x^{n})+n\sum^{n}_{i=1}\lambda_{j}f_{j}(\beta_{j}).$

Here we write $x^{n}$ to mean $x_{1},\ldots,x_{n}$ for ease of notation. Let

$\displaystyle a_{n}=\text{max}\\{\lambda_{j}(\mathcal{D}_{n}):1\leq j\leq
p_{0}\\},$ $\displaystyle
b_{n}=\text{min}\\{\lambda_{j}(\mathcal{D}_{n}):p_{0}<j\leq p\\}.$

Below we state our six regularity conditions.

Condition 1

Each Fisher information matrix
$I(\beta|x_{i})=-E\left(\frac{\partial^{2}}{\partial^{2}\beta}ln\;\rho(Y_{i}|x_{i},\beta)\right)$
exists and is positive definite uniformly in $i$, i.e. is bounded above and
below. Also for some $B,b>0$, we have $BI_{p\times p}\geq I(\beta|x_{i})\geq
bI_{p\times p}$.

Condition 2

The log likelihood of $\rho(Y_{i}|x_{i},\beta)$ has a convergent second order
Taylor expansion. That is, for all $j,\ell=1,\ldots,p$, we have that
$\forall\beta\in\Omega$, $\forall x_{i}$ and as $\eta\rightarrow 0$,

$E\left[\sup_{\beta\in\cal{B}(\beta_{0},\eta)}\left|\frac{\partial^{2}}{\partial\beta_{j}\partial\beta_{\ell}}ln\;\rho(Y_{i}|x_{i},\beta)\right|\right]\rightarrow
I(\beta_{0}|x^{n}),$

where $\cal{B}(\beta_{0},\eta)$ is the ball centered at $\beta_{0}$ with
radius $\eta>0$ in Euclidean distance.

Condition 3

For any $\epsilon$, $\exists N$ such that $\forall n\geq N$

$\sup\left|\frac{1}{n}\sum^{n}_{i=1}I(\beta|x_{i})-I(\beta|x^{\infty})\right|<\epsilon$

where $I(\beta|x^{\infty})$ is positive definite.

Condition 4

There exists an increasing sequence of compact sets $C=C_{n}$ in the parameter
space and constants $M=M_{n}\in\mathbb{R}^{+}$ such that for all $n$,
$\sup_{\beta_{j}\in C_{n}}|f_{j}^{\prime}(\beta_{j})|\leq M_{n}$. That is, the
penalty term has a uniformly bounded first derivative in $j$.

Condition 5

For $\epsilon$ in a neighborhood around zero,
$\sup_{j}\sup_{|\epsilon|\rightarrow 0}f^{\prime}_{j}(\epsilon)=0$.
Consequently, if $\beta_{j}=0$, then we have
$f^{\prime}_{j}(0)=\sup_{|\epsilon|\rightarrow 0}f^{\prime}_{j}(\epsilon)=0$
uniformly in $j$.

Condition 6

There exists uniform Taylor expandability for $f_{j}(\beta_{j})$. That is, for
$h\in\mathbb{R}$,
$f_{j}(\beta_{0}+h)-f_{j}(\beta_{0})=f_{j}^{\prime}(\beta_{0})h+o_{j}(1)$
uniformly in $j$ where $\sup_{j}o_{j}(1)\rightarrow 0$.

Our main result for penalized likelihoods is the following.

######  Theorem 2.1.

(Oracle Property) Assume Conditions 1–6 hold. Suppose
$\sqrt{n}a_{n}\rightarrow 0$ in which $a_{n}$ satisfies
$a_{n}=\frac{1}{h(n)\sqrt{n}}$ where $\frac{1}{h(n)\sqrt{n}}\rightarrow 0$ and
$h(n)\rightarrow\infty$, and suppose $\sqrt{n}b_{n}\rightarrow\infty$ in which
$b_{n}=\frac{g(n)}{\sqrt{n}}$ where $\frac{g(n)}{\sqrt{n}}\rightarrow 0$ and
$g(n)\rightarrow\infty$. Then, the estimator
$\hat{\beta}=(\hat{\beta_{1}}^{\prime},\hat{\beta_{2}}^{\prime})^{\prime}$
satisfies $P(\hat{\beta_{2}}=0)\rightarrow 1$ and
$\sqrt{n}\hat{I}_{1}(\beta_{1}|x^{n})(\hat{\beta_{1}}-\beta_{1})\rightarrow\mathcal{N}(0,I_{1}(\beta_{1}|x^{n}))$,
where $\beta_{1}$ is the component of $\beta$ containing all the nonzero
elements.

For proof see Supplement A, Sec. 7.1.

Therefore, we can construct an oracle procedure by choosing any
$f_{j}(\beta_{j})$ that satisfy Conditions 4 and 5. We can use the same
function $f_{j}(\beta_{j})$ for each $j$ , if we want, or we can have a
different penalty for each parameter. Thus, we are also not restricted, at
least theoretically, in our choice of the shrinkage methods except perhaps by
requiring the OP.

###  2.2 Penalized Empirical Risk

We now extend the OP results for penalized likelihoods to penalized empirical
risks. This is arguably a more general setting because we allow a wider range
for the function of the data as well as an arbitrary penalty function.

Consider the same regression scenario as Section 2.1 and a distance
$d(y_{i}-x_{i}^{\prime}\beta)$. Let
$R(\beta|x^{n})=\frac{1}{n}\sum^{n}_{i=1}d(y_{i}-x_{i}^{\prime}\beta)$ be the
empirical risk of the observations $(x_{1},Y_{1}),\ldots,(x_{n},Y_{n})$ where
$x^{n}$ is shorthand notation for $x_{1},\ldots,x_{n}$ and denote the
penalized empirical risk objective function as

$Q(\beta)=R(\beta|x^{n})+n\sum^{p}_{j=1}\lambda_{j}f_{j}(\beta_{j}).$

We specify some similar but slightly different conditions before getting into
the proofs. These conditions are similar to those in Section 2.1, but they are
different and necessary for the cases we discuss in this section.

Condition 7

The empirical risk $\frac{1}{n}\sum^{n}_{i=1}d(y_{i}-x_{i}^{\prime}\beta)$ has
a convergent second order Taylor expansion. That is, for all
$j,\ell=1,\ldots,p$ we have that $\forall\beta\in\Omega$ and as
$\eta\rightarrow 0$,

$E\left[\sup_{\beta\in
B(\beta_{0},\eta)}\left|\frac{\partial^{2}}{\partial\beta_{j}\partial\beta_{\ell}}\frac{1}{n}\sum^{n}_{i=1}d(y_{i}-x_{i}^{\prime}\beta)\right|\right]\rightarrow
E_{L}\left[R^{\prime\prime}(\beta|x^{\infty})\right]=I^{*}(\beta|x^{\infty})$

where

$R^{\prime\prime}(\beta|x^{n})=\frac{\partial^{2}}{\partial^{2}\beta}R(\beta|x^{n})=\frac{\partial^{2}}{\partial^{2}\beta}\left(\frac{1}{n}\sum^{n}_{i=1}d(y_{i}-x_{i}^{\prime}\beta)\right)$

and for some $B,b>0$, we have $BI_{p\times p}\geq I^{*}(\beta|x^{\infty})\geq
bI_{p\times p}$. In other words, $I^{*}(\beta|x^{\infty})$ is positive
definite and we abbreviated it to $I^{*}(\beta)$.

We see that $I^{*}(\beta)$ is a sort of analog to the Fisher information
matrix in the context of empirical risks. It is not necessarily the Fisher
information matrix; however, if our empirical risk happens to correspond to a
log-likelihood, then it is possible that $I^{*}(\beta)=I(\beta)$.

Condition 8

There exists an increasing sequence of compact sets $C=C_{n}$ in the parameter
space and constants $M=M_{n}\in\mathbb{R}$ such that for all $n$,
$\sup_{\beta_{j}\in C_{n}}|f_{j}^{\prime}(\beta_{j})|\leq M_{n}$. That is, the
penalty term has a uniformly bounded first derivative in $j$.

Condition 9

For $\epsilon$ in a neighborhood around zero,
$\sup_{j}\sup_{|\epsilon|\rightarrow 0}f_{j}^{\prime}(\epsilon)=0$.
Consequently, if $\beta_{j}=0$, then we have
$f_{j}^{\prime}(0)=\sup_{|\epsilon|\rightarrow 0}f_{j}^{\prime}(\epsilon)=0$
uniformly in $j$.

Condition 10

There exists uniform Taylor expandability for $f_{j}(\beta_{j})$. That is, for
$h\in\mathbb{R}$,
$f_{j}(\beta_{0}+h)-f_{j}(\beta_{0})=f_{j}^{\prime}(\beta_{0})h+o_{j}(1)$
uniformly in $j$ where $\sup_{j}o_{j}(1)\rightarrow 0$.

Again the notion of uniform Taylor expandability in $j$ here requires the
error terms $o_{j}(1)$ in Condition 10 go to zero at the same time.

The following lemma identifies the class of distance functions we consider.

######  Lemma 1.

Let $d(u)$ be an even distance function with a unique minimum at 0. If $u$
comes from some distribution with pdf $f_{U}(u)$ that is symmetric about zero
with support $[-a,a]$ for $a\in\mathbb{R}$, then $E_{U}[d^{\prime}(u)]=0$.

###### Proof.

Since $d(u)$ is an even function, we know $d^{\prime}(u)$ is an odd function.
By definition of expectation,

$E_{U}(d^{\prime}(u))=\int^{a}_{-a}d^{\prime}(u)f_{U}(u)du.$

Note that since $f_{U}(u)$ is a symmetric distribution, $f_{U}(u)=f_{U}(-u)$,
so $f_{U}(u)$ is an even function. Let $g(u)=d^{\prime}(u)f_{U}(u)$. Then
$g(u)$ is an odd function, and we have $\int^{a}_{-a}g(u)du=0,$ so

$E_{U}(d^{\prime}(u))=\lim_{a\rightarrow 0}\int^{a}_{-a}g(u)du=0.$

∎

Note that Lemma 1 is true for any even distance function $d(u)$, and this is
useful for us because we focus on the distance function
$d(y_{i}-x_{i}^{\prime}\beta)$ which is even due to the symmetry of
$(y_{i}-x_{i}^{\prime}\beta)$.

Our main result for penalized empirical risks is the following.

######  Theorem 2.2.

(Oracle Property) Assume Conditions 7–10 and Lemma 3 hold. Suppose
$\sqrt{n}a_{n}\rightarrow 0$ where $a_{n}$ satisfies
$a_{n}=\frac{1}{h(n)\sqrt{n}}$ such that $\frac{1}{h(n)\sqrt{n}}\rightarrow 0$
and $h(n)\rightarrow\infty$, and suppose $\sqrt{n}b_{n}\rightarrow\infty$
where $b_{n}$ satisfies $b_{n}=\frac{g(n)}{\sqrt{n}}$ such that
$\frac{g(n)}{\sqrt{n}}\rightarrow 0$ and $g(n)\rightarrow\infty$. Then the
estimator
$\hat{\beta}=(\hat{\beta_{1}}^{\prime},\hat{\beta_{2}}^{\prime})^{\prime}$
satisfies $P(\hat{\beta_{2}}=0)\rightarrow 1$ and
$\sqrt{n}\hat{I^{*}}_{1}(\tilde{\beta}_{10}|x^{n})(\hat{\beta_{1}}-\beta_{10})\rightarrow\mathcal{N}(0,I_{1}(\beta_{10}|x^{n}))$
.

For proof see Supplement A, Sec. 7.2.

Taken together, Theorems 2.1 and 2.2 give us insight into how large is the
class of oracle procedures. Previously it seemed oracle procedures were rare,
isolated choices of priors. Now, even though we have not characterized the
class of all oracle procedures, we can see that the conditions for a procedure
to have the oracle property are quite general, allowing a large range of
likelihoods, distances, and priors.

Since there are obviously many oracle procedures, the key question becomes
which one to choose in a given context. We propose choosing a ’best’ oracle
procedure by optimizing a stability criterion over both a class of penalty
functions and a class of distance functions. Allowing for $f_{j}(\beta_{j})$
to be a different penalty for each parameter, as long as they are uniformly
Taylor expandable and have similar properties, is powerful because we can
choose variable-dependent penalties. This is a more general sense of
adaptivity than each $\beta_{j}$ merely having its own shrinkage parameter
$\lambda_{j}$. We believe a natural choice for $f_{j}(\beta_{j})$ might be
something like $|\beta_{j}|^{1-|corr(y_{i},x_{i})|}$ because it is data
driven. We have not done this here, and we leave it for future work. Here we
focus on comparing existing methods and then using genetic algorithms to find
optimal choices for the $f_{j}(\beta_{j})$’s.

##  3 Computational Comparisons

Every shrinkage method for linear models generates a predictor of the form
(1.3) ,i.e., $\hat{Y}(x_{n+1})=x_{n+1}^{T}\hat{\beta}$ where the estimate
$\hat{\beta}$ of $\beta$ comes from choosing $\lambda$ and the $w_{j}$’s (or
other data-driven parameters). It is well-known that many shrinkage methods
(LASSO, EN, etc) zero-out coefficients $\beta_{j}$ and so do variable
selection as well as estimation. Here, we look only at the instability of
predictive error and the accuracy of variable selection.

Following [12] we add random $N(0,\tau^{2})$ noise to the $y_{i}$’s in
${\cal{D}}$ and denote the partition of the perturbed data by
${\cal{D}}_{n}(\tau)={\cal{D}}_{train}(\tau)\cup{\cal{D}}_{test}(\tau)$. For
any predictor $\hat{Y}$, we define its instability to be

$S(\hat{Y})_{\tau}=\sqrt{\frac{1}{n_{test}}\sum_{i\in{\cal{D}}_{test}(\tau)}(y_{i}-\hat{Y}_{\tau}(x_{i}))^{2}}$

where $\hat{Y}_{\tau}$ means we have formed a predictor using
${\cal{D}}_{train}(\tau)$. In the computations we present in this section we
used $\tau_{k}=.2k$, $k=1,\ldots,8$, to generate instability curves of the
form $(k,S(\hat{Y})_{\tau_{k}})$, and looked for patterns.

Intuitively, perturbing the $Y$’s by adding normal noise should only increase
$S(\hat{Y})_{\tau_{k}}$, i.e., in curves should generally increase with
$\tau$. Of course, we prefer instability curves that are small – less
instability upon perturbation suggests a better predictor. However, if a
instability curve decreases with $\tau$ then perturbation of $Y$ is making the
predictor more stable. We take this to mean the predictor is discredited for
some reason. We suggest this behavior arises when the predictor has omitted or
included terms incorrectly or has poorly chosen coefficients. We prefer
predictors with instability curves that are lower than the instability curves
of competing predictors and smoothly increase slowly with $\tau$.

Our basic computational procedure is as follows. Fix a number $K$ of values of
$\tau$ to form the points on the instability curve and a (large) number $L$
for the number of iterations to be averaged. For each $k=1,\ldots,K$, let
$\ell=1,\ldots,L$. Instability curves for a given predictor can generically be
formed by the following steps.

  1. 1.

For each $\ell$, randomly split ${\cal{D}}_{n}$ in to ${\cal{D}}_{train}$ and
${\cal{D}}_{test}$.

  2. 2.

Perturb the $y$-values in ${\cal{D}}_{train}$ and ${\cal{D}}_{test}$ using
$N(0,\tau^{2}_{k})$ noise. Call the results ${\cal{D}}_{train,\tau_{k}}$ and
${\cal{D}}_{test,\tau_{k}}$, respectively.

  3. 3.

Using ${\cal{D}}_{train,\tau_{k}}$ form the competing predictors denoted
$\hat{Y}_{\tau_{k}}$.

  4. 4.

Using ${\cal{D}}_{test,\tau_{k}}$ obtain $S(\hat{Y}_{\tau_{k}})$ for each
predictor.

  5. 5.

Let $S(\hat{Y}_{\tau_{k},\ell})$ be the $\ell$-th value of
$S(\hat{Y}_{\tau_{k}})$.

  6. 6.

For each predictor and $k$, find the sample mean:

$S_{\tau_{k}}(\hat{Y})=\frac{1}{L}\sum_{\ell=1}^{L}S(\hat{Y}_{\tau_{k},\ell}).$

  7. 7.

Plot $S_{\tau_{k}}(\hat{Y})$ as a function of $\tau_{k}$ for each predictor.

In this section all simulated data come from

$Y=X\beta+\epsilon$

where the rows $X_{i}=x_{i}$ of $X$, for $i=1,\ldots,n$, are $MVN_{p}(0,M)$
for various choices of variance matrix $M$ or are IID $\sim t_{3}$ to see the
effect of heavier tails. The parameter values for $\beta\in\mathbb{R}^{p}$ are
IID $N_{p}(4,1)$. So, $\dim(Y)=n$, $\dim(X)=n\times p$, and
$\dim(\epsilon)=n$. We use two choices for the distribution of $\epsilon_{i}$,
$N(0,1)$ and $t_{3}$, to represent light and heavy tails in the error,
respectively. We are concerned mainly with the case $p<n$, but include cases
$p>n$ for completeness.

We compared predictors from ten different shrinkage methods as well as a full
linear model. Seven of the shrinkage methods had the OP, namely, LAD-LASSO
(uses $L^{1}$, in place of $L^{2}$ in ALASSO, Wang et al. 2007), ALASSO (Zou
2006), SCAD1 ($L^{1}$ with a truncated $L^{1}$ penalty, Fan and Li 2001), AEN
(Zou and Zhang 2009), ASCAD1 (Dustin 2020: $\|Y-X\beta\|_{1}+\sum
p_{\lambda_{j}}(|\beta_{j}|)$) i.e., an adaptive SCAD penalty on each
$\beta_{j}$), SCAD2 (replace $L^{1}$ with $L^{2}$ in SCAD1) ], and MCP (Zhang
2010). Three of them did not have the OP, namely ridge regression (RR), LASSO
and EN.

Finding the instability curves for the ten shrinkage methods required three
different pieces of software. First, for LASSO, RR, EN, ALASSO, AEN, we used
the glmnet package (see [13]) in RStudio Ver. 1.2.5033. Second, for LAD-LASSO,
SCAD1, and ASCAD1 we used the rqPen package in R, see [14]. The specific
function used (cv.rq.pen() ) requires $X$ to be nonsingular so we only
implement these methods (as well as LM) in the $n>p$ case. Third, for MCP and
SCAD2 we used code that implemented a local linear approximation (LLA).
444Development of the LLA code was by L. Xue and reported in [15]. We are
grateful to these authors for giving us their code.

Because we use different packages, we had to split the data differently for
different methods. The LLA code requires an explicit separation of
${\cal{D}}_{train}$ into two sets say
${\cal{D}}_{train}={\cal{D}}_{train,\lambda}\cup{\cal{D}}_{train,\beta}$ so
that a small portion of the data is set aside for estimating $\lambda$.
However, both glmnet and rqPen approximate $\lambda$ internally and we use the
default methods of the procedures, so the splitting of ${\cal{D}}_{train}$ is
done internally to the programs. Below we list the choices we made for the
explicit splitting.

As noted in Sec. 1, we followed [8] for the adaptive methods and chose
$\hat{w}_{j}=1/|\hat{\beta}_{j,OLS}|$ for $p<n$. When we do not have enough
data to implement OLS, i.e. $p>n$, we used $1/|\hat{\beta}_{j,SCAD2}|$. In
fact, the computations can be done, at least in principle, using any
$\sqrt{n}$-consistent estimator of the $\beta_{j}$’s for the $w_{j}$’s. The
reasoning behind this choice comes from the assumptions on $a_{n}$ and $b_{n}$
in Theorems 2.1 and 2.2.

We focus here on a high dimensional setting of $p$ relative to $n$. Thus, we
set $p=100$ and consider three sparsity levels,10%, 50%, and 90%, and four
sample sizes $n=40,75,150,500$. For each $n$ we let $L=100$ for the
instability computations. That is we average over 100 datasets to get a
instability value for each perturbation level. We use a training data set to
form the predictor and a testing data set to evaluate its performance, and for
the various sample sizes we set $({\cal{D}}_{train},{\cal{D}}_{test})$ =
$(28,12)$, $(60,15)$, $(120,30)$, and $(350,150)$, respectively. For the $LLA$
methods we allocate 10% of the ${\cal{D}}_{train}$ data to
${\cal{D}}_{train,\lambda}$ and the remaining $90\%$ to
${\cal{D}}_{train,\beta}$. Our choices respect the fact that for the methods
to be comparable, they must be trained with exactly the same data, and
evaluated on a hold out set that has not been used in the training process.

The next four subsections present our simulation results for the four sample
sizes and three sparsity levels. The final subsection summarizes our
recommendations before we turn to a benchmark data example in Sec. 4.

###  3.1 Sample Size $n=40$

Our first example uses the value $n=40$ that is relatively small compared to
$p=100$. Here, we examine how good the variable selection, estimates, and thus
predictions are for various sparsity levels where sparsity means the
percentage of zero coefficients in the model used to generate the data.

We use two assessments for this. First, we generate instability curves. Then
we also look ‘inside’ the predictor to see which variables were included
correctly or incorrectly. These two assessments are roughly analogous to
variance and bias. However, as noted in Sec. 1,we think of instability curves
as more important because they reflect ‘variability’ and partially reflect
bias: If a predictor includes a variable incorrectly it can still be
downweighted by small parameter estimates and if the predictor fails to
include a correct variable (that is important) the instability curve can
increase to detect it.

Fig. 1 shows that there is overall more instability with heavy tails, and,
independent of heaviness of the tails, as sparsity increases the methods
become more stable. We observe two clusters of methods. Namely, the non-
adaptive methods in the lower cluster (L,EN,RR) and the adaptive methods in
the upper cluster(AL,AEN,SCAD2, MCP). This is consistent with the fact that we
have far fewer observations than predictors, so estimating the extra
hyperparameters accurately is difficult. Overall, this figure suggests that RR
and EN perform best for sparsity levels 0.1 and 0.5 while EN is best for
sparsity 0.9, although RR is not discredited, at least for light tails.

This is seen for both light and heavy tails. A possible explanation for the
good performance of RR is that there is not enough information in the small
sample relative to the number of predictors to set any coefficients to zero.
Thus RR performs well because it does not exclude any variables. EN, a mixture
of LASSO and RR, also performs well because it may set fewer coefficients to
zero than other non-adaptive methods.

It may seem odd to say RR is not discredited in the upper right panel of Fig.
1. However, all the other curves decrease markedly as the perturbation
increases suggesting they have chosen poor models. Also, RR is smallest at
perturbation zero.

The curves increase after the perturbation is large enough that it outweighs
the poor selection of the models. Also, it is seen that there are some jumps
in the curves. These are usually at the beginning when the perturbation
changes from zero to a positive number. We suggest these jumps indicate a lack
of stability under perturbation meaning that the shrinkage method they
represent is unstable in terms of choosing a good predictor.

---
Figure 1: $n=40$ From left to right the sparsity increases from 10%, to 50%,
to 90%. From top to bottom the heaviness of the tails of $X$ and $\epsilon$
increases from normal to $t_{3}$. The distribution of the 100 IID outcomes for
the parameter $\beta$ is $N(4,1)$. In this simple case, we used the identity
matrix in the normal to generate the matrix $X$.

Next, we look at the predictors themselves, on average, to help us interpret
the instability curves. Table 1 compares each of the methods for both light
(Li) and heavy (H) tails based on the total percentage of coefficients the
methods sets to 0 on average (Tot; ideally equal to the sparsity level), what
percentage of those set to 0 were truly 0 (Tr; ideally 100%) and what
percentage of those set to 0 were truly non-zero (Fa; ideally zero). Note that
we are focused on the sparsity in this setting, so we frame these tables in
terms of the percentage of variables excluded by a method. Instead, if we wish
to describe the variables included, we simply subtract each entry in the table
from 1.

Sparsity | .1 | .5 | .9 |
---|---|---|---|---
|  | Tot | Tr | Fa | Tot | Tr | Fa | Tot | Tr | Fa
Li | L | .95 | .95 | .95 | .92 | .93 | .90 | .95 | .97 | .76
Li | AL | .88 | .85 | .88 | .89 | .90 | .87 | .91 | .95 | .56
Li | EN | .87 | .86 | .87 | .82 | .84 | .80 | .91 | .94 | .56
Li | AEN | .88 | .84 | .88 | .88 | .90 | .80 | .90 | .94 | .55
Li | SCAD2 | .82 | .78 | .82 | .82 | .84 | .80 | .87 | .91 | .49
Li | MCP | .83 | .78 | .83 | .84 | .85 | .82 | .88 | .92 | .50
H | L | .97 | .98 | .97 | .93 | .95 | .90 | .96 | .97 | .88
H | AL | .89 | .92 | .89 | .90 | .93 | .88 | .91 | .99 | .69
H | EN | .93 | .96 | .93 | .88 | .91 | .85 | .88 | .93 | .72
H | AEN | .87 | .90 | .88 | .90 | .93 | .87 | .90 | 1 | .66
H | SCAD2 | .83 | .85 | .82 | .84 | .86 | .81 | .87 | .86 | .59
H | MCP | .84 | .86 | .84 | .84 | .88 | .82 | .88 | .90 | .63
Table 1:  Variable selection performance for $n=40$

From Table 1, it is seen that none of the methods performs well. The values in
the Tot columns for low and medium sparsity are much too high which leads to
high values in the Tr columns, trivially. The high values in the Fa columns
indicate that all methods are typically excluding far too many terms. There
are some very mild exceptions – SCAD2 or MCP – that perform noticeably better
than other methods but even so their performance is poor. For high sparsity
the methods do exclude fewer correct variables than for low sparsity. Table 1
agrees with Fig. 1 in that variable selection is generally worse for heavier
tails. Curiously, the better performing methods in terms of variable selection
are different from the better performing methods in terms of instability
curves. We attribute the difference to how well coefficients are estimated
rather than to how well variables are selected. Indeed, Table 1 shows that all
methods usually selected far fewer variables than necessary, e.g., for L with
.1 sparsity, five variables were included on average but 90 should have been.

We also considered two cases where the $x_{j}$’s have a nontrivial dependence
structure. Specifically, we set $X\sim MVN_{100}(0,M)$ where $M$ is
tridiagonal or Toeplitz, both with light tailed errors. We do not show the
plots for these cases, howevver we give our recommendations in Subsec. 3.5.

As in the independence case, the methods performed generally better as
sparsity increased. For figures analogous to those in Fig. 1 but for
tridiagonal variance matrices, we found that for low and medium sparsity, EN
and RR were typically best and close to each other. For high sparsity we got
the same results as in the upper right panel of Fig. 1. When the variance
matrix of $X$ was Toeplitz, we found that EN is performing best for all low
and medium sparsity levels. For high sparsity SCAD2 and MCP performed nearly
the same as EN. Interestingly, SCAD2 and MCP were amongst the worst performing
methods in the independence case but with the Toeplitz structure they were
only slightly worse than EN. We suggest that the degree of dependence in the
tridiagonal case is not much more than in the independence case. We also
suggest that when there is high enough dependence e.g., the Toeplitz case,
penalties that have some data dependence (but not too much) such as EN, SCAD2,
and MCP are best able to generate predictions.

In terms of variable selection as indicated in Table 1, we found that the
results for the tridiagonal and Toeplitz were similar to the independence case
across all methods and sparsity levels except for EN in the Toeplitz case
where EN performed noticeably better than the other methods across all
sparsity levels. In addition, L does well in the high sparsity case. We
suggest that being able to choose what the penalty looks like as in EN gives
some advantage over the other methods. We return to this point in Sec. 5.

###  3.2 Sample Size $n=75$

Next we consider a second example where $p>n$. That is, we still have fewer
observations than explanatory variables, but $n$ is much closer to $p$ than in
Subsec. 3.1. Examining the instability plots in Fig. 2, we see that as in Fig.
1, the methods become more stable as sparsity increases and less stable as the
tails become heavier.

Similar to the $n=40$ case, we have the same clusters of methods, namely EN,
RR, and L as the lower (better) cluster with AEN, AL, MCP and SCAD2 as the
worse (upper) cluster. EN outperforms all other methods with low and medium
sparsity regardless of the heaviness of the tails. For high sparsity all
methods are roughly comparable except for RR which performs poorly. That is
the $n=40$ and $n=75$ cases are qualitatively similar apart from the light
tail high sparsity case with $n=40$ where model instability seems to dominate.
Specifically, for $n=40$ EN was the ‘best’ although there was good reason to
prefer RR, but for $n=75$ all methods but RR are essentially indistinguishable
and RR performs poorly.

It is seen that the error decreases as both the sparsity and $n$ increase.
Accordingly, as both increase, the predictive error may converge to $\sigma$,
the SD of $\epsilon$.

---
Figure 2: $n=75$: From left to right the sparsity increases from .1, to .5, to .9. From top to bottom the heaviness of the tails of $X$ and $\epsilon$ increases from normal to $t_{3}$. The distribution of the 100 IID outcomes for the parameter $\beta$ is $N(4,1)$. Again, we used the identity matrix in the normal distribution for $X$. Sparsity | .1 | .5 | .9 |
---|---|---|---|---
|  | Tot | Tr | Fa | Tot | Tr | Fa | Tot | Tr | Fa
Li | L | .75 | .89 | .73 | .83 | .92 | .75 | .76 | .84 | 0
Li | AL | .77 | .88 | .76 | .79 | .88 | .71 | .90 | 1 | 0
Li | EN | .60 | .78 | .58 | .70 | .81 | .58 | .76 | .85 | 0
Li | AEN | .76 | .87 | .74 | .79 | .88 | .70 | .90 | 1 | 0
Li | SCAD2 | .67 | .78 | .66 | .71 | .81 | .61 | .90 | 1 | 0
Li | MCP | .69 | .80 | .68 | .72 | .82 | .63 | .89 | .99 | 0
H | L | .90 | .92 | .90 | .85 | .91 | .90 | .74 | .83 | 0
H | AL | .83 | .99 | .82 | .83 | .90 | .76 | .90 | 1 | 0
H | EN | .85 | .81 | .85 | .77 | .84 | .71 | .74 | .82 | 0
H | AEN | .82 | 1 | .81 | .82 | .89 | .74 | .90 | 1 | 0
H | SCAD2 | .72 | .79 | .71 | .74 | .82 | .66 | .89 | .99 | .01
H | MCP | .74 | .80 | .73 | .75 | .82 | .68 | .89 | .98 | .01
Table 2:  Variable selection performance for $n=75$

Table 2 shows that, as with Table 1, no method performs variable selection
well for low to medium sparsity. For high sparsity, all methods were very much
improved. However, L and EN were the worst – the only methods not having the
OP. The other methods, AL, AEN, SCAD2, and MCP, have the OP and perform
roughly equally well. As a generality, the methods performed better for light
tailed than for heavy tailed distributions. As before, the better performing
methods in the instability curves (EN, RR, L) tended to perform worse in terms
of variable selection and we attribute this to better parameter estimation in
the instability curves since the methods generally included more variables
than necessary. That is, a variable may be included with an estimated
coefficient near zero.

When we included dependence via tridiagonal matrices, we found results similar
to the independent case. Namely, the instability curves show that EN performed
best at all sparsity levels, roughly tying with most other methods for high
sparsity. The key difference was that RR tended to perform poorly overall.
When we switched to Toeplitz matrices, across all sparsity levels SCAD2 was
best with MCP being a close second and virtually tying with SCAD2 with high
sparsity. EN, which was overall best for $n=40$, was noticeably worse than
both SCAD2 and MCP across all sparsity levels. We suggest that since we have
more data, the optimality of the methods that have optimality properties
(other than the OP) is effective.

In terms of variable selection, the tridiagonal case was similar to the
independence case but slightly better. This was surprising and difficult to
interpret. In the Toeplitz case, EN, SCAD2 and MCP are noticeably better than
the other methods for low and medium sparsity. For high sparsity, L also
performs well. This is the same as the Toeplitz case with $n=40$. Thus,
overall, the results for dependence cases with $n=75$ are very close to the
corresponding results for $n=40$.

###  3.3 Sample Size $n=150$

This is our first case where $n>p$, making it qualitatively different from the
earlier two subsections. Here we implement LM’s as well as shrinkage
techniques. Note that the “penalty” associated with LM’s is identically zero
and corresponds to a uniform prior.

Fig. 3 parallels Figs. 2 and 1, but introduces four extra methods, LM’s,
SCAD1, ASCAD1, and L-L. For low sparsity, LM’s start out giving the best
results in predictive instability. However, the instability curve for LM’s
rises faster than for its competitors and arguably SCAD1 is equivalent. For
medium sparsity, LM’s is roughly in the middle and SCAD1 is the best. For high
sparsity, SCAD1 ties with the other top methods. LM’s perform worse as
sparsity increases. As in earlier cases, all methods that have nontrivial
penalties perform better for light tails than heavy tails and as sparsity
increases. Also, as the sparsity level increases, the variability amongst the
methods decreases. Note that RR tends to perform poorly since there is enough
data that the sparsity has an effect.

A key difference between the $n=150$ case and the $n=75,40$ cases is that the
adaptive methods are generally performing better than the nonadaptive methods.
For instance, AEN and AL are performing better for small perturbations than EN
or L, respectively. When the perturbations are too high, it makes sense that
the non-adaptive version of a penalty will perform better that their adaptive
versions because they are less affected by the noise; they use fewer
estimators. One can argue that the perturbation level at which the curves for
non-adaptive penalties and their adaptive versions cross represents the
largest reasonable perturbation that should be considered for that penalty.
Moreover, the OP is not a determining factor for performance: Some methods
with the OP perform well and some do not. Some methods that do not have the OP
perform better than other methods that do.

---
Figure 3: $n=150$: From left to right the sparsity increases from .1, to .5, to .9. From top to bottom the heaviness of the tails of $X$ and $\epsilon$ increases from normal to $t_{3}$. The distribution of the 100 IID outcomes for the parameter $\beta$ is $N(4,1)$. We used the identity matrix in the normal distribution for $X$. Sparsity | .1 | .5 | .9 |
---|---|---|---|---
|  | Tot | Tr | Fa | Tot | Tr | Fa | Tot | Tr | Fa
Li | L | .04 | .38 | 0 | .30 | .60 | 0 | .81 | .90 | 0
Li | AL | .11 | 1 | .01 | .50 | 1 | 0 | .89 | .99 | 0
Li | EN | .03 | .30 | 0 | .29 | .59 | 0 | .81 | .90 | 0
Li | AEN | .11 | 1 | .01 | .50 | 1 | 0 | .88 | .98 | 0
Li | ASCAD1 | .11 | .82 | .03 | .48 | .95 | 0 | .89 | .99 | 0
Li | L-L | .01 | 0 | 0 | .04 | .07 | 0 | .83 | .92 | 0
Li | SCAD1 | .03 | .30 | 0 | .42 | .85 | 0 | .89 | .99 | 0
Li | SCAD2 | .12 | .74 | .07 | .51 | 1 | .02 | .90 | 1 | 0
Li | MCP | .11 | .67 | .05 | .50 | .99 | .01 | .89 | .99 | 0
H | L | .04 | .33 | 0 | .24 | .49 | 0 | .83 | .92 | 0
H | AL | .10 | 1 | 0 | .50 | .99 | 0 | .90 | 1 | 0
H | EN | .03 | .28 | .43 | .24 | .49 | 0 | .82 | .92 | 0
H | AEN | .10 | 1 | .50 | .50 | .99 | 0 | .90 | 1 | 0
H | ASCAD1 | .13 | .87 | .05 | .49 | .97 | .01 | .90 | 1 | .01
H | L-L | .01 | .11 | 0 | .10 | .19 | .01 | .87 | .96 | .01
H | SCAD1 | .07 | .72 | 0 | .48 | .95 | 0 | .90 | 1 | 0
H | SCAD2 | .35 | .72 | .30 | .49 | 98 | .01 | .90 | 1 | 0
H | MCP | .33 | .70 | .29 | .48 | .96 | .01 | .90 | 1 | 0
Table 3:  Variable selection performance for $n=150$

For light tails and low sparsity, Table 3 shows that AL and AEN are the best
methods in terms of variable selection. For medium sparsity they remain the
best but MCP and SCAD2 are almost as good. For high sparsity, AL, AEN, ASCAD1,
SCAD1, SCAD2, and MCP are essentially the same. Note that all methods have
zero in the Fa column meaning they never exclude variables that are important
For heavy tails, the results for low sparsity are the same but for medium
sparsity more methods are nearly as good as the AL and AEN. For high sparsity,
the only methods that are lagging in performance are L, EN, and L-L; the
others are essentially equivalent.

Compared to Table 2 we see that all methods improved in variable selection,
which is not surprising, but that the adaptive methods improved more. This is
true for the instability curves as well. LM’s and RR are not included in Table
3 because they don’t do variable selection.

Comparing the conclusions from Fig 3 and Table 3, we see that SCAD1 performs
best predictively although not necessarily in variable selection. In Table 3,
SCAD1 never excludes a variable that should be included. Thus, its lesser
performance in terms of variable selection may be ameliorated by its parameter
estimation.

Again, we considered two dependence cases with light tails, the tridiagonal
and the Toeplitz. For the tridiagonal case, the instability curves and the
variable selection table are qualitatively the same as for the independence
case. For the Toeplitz, LM’s, SCAD1, SCAD2, and MCP are the only methods that
perform well in terms of instability curves. In the corresponding table,
overall SCAD2 and MCP did best, but sometimes another method does best (but
then does poorly in terms of instability). Overall, variable selection in this
case is worse than in the independence case. This is virtually the opposite of
Toeplitz in the $n<p$ case where variable selection was improved. This case is
more in line with intuition because intuition corresponds to $n>p$.

###  3.4 Sample Size $n=500$

For completeness we also consider the case $n=500$ to identify the limiting
behavior of the methods. Fig. 4 and Table 4 have the same general properties
as the earlier Figures and Tables, namely, as sparsity increases instability
decreases. The methods improve as sparsity increases although the improvement
is not as dramatic as in the smaller sample cases. In the heavy tailed cases,
there is more variability.

In fact, many of the methods at this point are indistinguishable via Fig. 4 or
Table 4. So for descriptive purposes it is easier to identify the methods that
perform poorly rather than the ones that perform well. From the top row in
Fig. 4, the only poor methods are RR, ASCAD1,L-L for low and medium sparsity.
For high sparsity RR is clearly worst. From the bottom row in Fig. 4, RR is
clearly the worst in all cases. For low sparsity LM and SCAD1 (and nearly
SCAD2) are best. For medium sparsity only L-L (and RR) performs poorly. For
high sparsity the worst performers are RR and LM’s. Note that L and EN still
perform well even though they don’t have the OP. However, EN is a
generalization of L and L has some consistency properties, see [16], so the
good performance of L and EN is not surprising.

---
Figure 4: $n=500$: From left to right the sparsity increases from .1, to .5,
to .9. From top to bottom the heaviness of the tails of $X$ and $\epsilon$
increases from normal to $t_{3}$. The distribution of the 100 IID outcomes for
the parameter $\beta$ is $N(4,1)$. We used the identity matrix in the normal
distribution for $X$.

Table 4 shows that at $n=500$, most methods are performing variable selection
quite well. In fact, almost all the methods used that have the OP are nearly
perfect, on average, in performing variable selection. The two exceptions are
LL and ASCAD1 which tend to include too many variables. That said, these two
methods do not often exclude variables that are important, which allows them
to perform comparably well predictively. Even the methods that look worse in
terms of variable selection (L, EN) predict well because they retain all the
important variables. The only methods this table suggest could be ruled out
are in the low or medium sparsity case: L, EN, and L-L. have low entries in
the Tot and Tr columns. (They also have near zero entries in the Fa column.)
However, our comparison is predictive not based on variable selection and
these methods perform generally well in a predictive stability sense. It is
likely that the coefficients of the incorrectly included variables are quite
small.

Note that in this case, the errors at zero of the light tailed methods
approach the theoretical lower limit of the predictive error ($\sigma=1$ ).
For the heavier tailed methods, the smallest errors are below two. Since the
error terms here are $t_{3}$ with $\sigma=\sqrt{3}\approx 1.73$, again the
smallest errors are approaching their theoretical lower limit. (This is seen
in earlier cases but usually only for the highest sparsity.)

Sparsity | .1 | .5 | .9 |
---|---|---|---|---
|  | Tot | Tr | Fa | Tot | Tr | Fa | Tot | Tr | Fa
Li | L | .04 | .41 | 0 | .34 | .68 | 0 | .83 | .93 | 0
Li | AL | .10 | 1 | 0 | .50 | 1 | 0 | .90 | 1 | 0
Li | EN | .04 | .36 | 0 | .34 | .67 | 0 | .83 | .92 | 0
Li | AEN | .10 | 1 | 0 | .50 | 1 | 0 | .90 | 1 | 0
Li | ASCAD1 | .09 | .94 | .05 | .50 | 1 | 0 | .91 | 1 | .06
Li | L-L | .02 | .15 | .1 | .36 | .72 | .01 | .90 | 1 | .06
Li | SCAD1 | .09 | .93 | 0 | .48 | .97 | 0 | .90 | 1 | 0
Li | SCAD2 | .10 | 1 | 0 | .50 | 1 | 0 | .90 | 1 | 0
Li | MCP | .10 | 1 | 0 | .50 | 1 | 0 | .90 | 1 | 0
H | L | .06 | .60 | 0 | .28 | .56 | 0 | .86 | .95 | 0
H | AL | .10 | 1 | 0 | .50 | 1 | 0 | .90 | 1 | 0
H | EN | .06 | .56 | 0 | .27 | .55 | 0 | .86 | .95 | 0
H | AEN | .10 | 1 | 0 | .50 | 1 | 0 | .90 | 1 | 0
H | ASCAD1 | .09 | .96 | 0 | .49 | .97 | 0 | .90 | 1 | 0
H | L-L | .01 | .11 | 0 | .06 | .11 | 0 | .86 | .96 | 0
H | SCAD1 | .09 | .91 | 0 | .49 | .99 | 0 | .90 | 1 | 0
H | SCAD2 | .10 | 1 | 0 | .50 | 1 | 0 | .90 | 1 | 0
H | MCP | .10 | 1 | 0 | .50 | 1 | 0 | .90 | 1 | 0
Table 4:  Variable selection performance for $n=500$

###  3.5 Summary Tables and Recommendations

To specify our recommendations we assume the following, elaborated from Subec.
3.1. First, when $p>n$, the instability graphs are more important that the
tables because the variable selection is so poor for all methods that we must
rely solely on choosing the method that predicts the best (or least bad).
Thus, for $n=40$ we just consider the figures in choosing the best methods.
For $n=75$, unless one specifies a relative weighting among Tot, Tr, and Fa,
there is usually no clearly preferred method. In these cases, we have used the
instability curves to select from amongst the numerically best shrinkage
methods.

Second, for larger sample sizes, and low to medium sparsity cases, we use both
the instability curves and variable selection tables to make recommendations.
That is, we choose the methods that are both selecting variables most
appropriately and have the lowest instability generally. Again, there are
cases where we have to weight Tot, Tr, and Fa. We continue to believe that Fa
is relatively more important than Tot and Tr and that Tot is relatively more
important than Tr.

Two limitations of our recommendations are that in practice i) the level of
sparsity and ii) the variable selection tables cannot be known. Accordingly,
it is only the instability curve that is available. Thus, our recommendations
are based primarily on the instability curves but taking into account the
variable selection tables. Moreover, there are other limitations that a
potential user should take into account e.g., our restriction to linear
models.

We begin our recommendations for small sample sizes with Table 5 (based on our
computations from Subsec. 3.1). Recall that for $n=40$ no method worked
particularly well. However, the overall best performing of these poor methods
was RR.

Sparsity | .1 | .5 | .9
---|---|---|---
Ind, Li | RR/EN | RR/EN | RR/EN
Ind, H | RR/EN | RR/EN | RR/EN
Tri | RR/EN | RR/EN | RR/EN
Toe | RR | RR | RR/SCAD2/MCP
Table 5:  Best performing shrinkage methods for $n=40$.

Recommendations for the next sample size, $n=75$, are given in Table 6 (based
on the computations in Subsec. 3.2). EN was typically the preferred method for
low to medium sparsity. For high sparsity, all methods besides RR performed
well. Interestingly, we found SCAD2 or MCP performed best when there is strong
dependence as in the Toeplitz structure.

Sparsity | .1 | .5 | .9
---|---|---|---
Ind, Li | EN | EN | Not RR
Ind, H | EN | EN | Not: RR
Tri | EN | EN | Not: RR
Toe | SCAD2/MCP | SCAD2/MCP | SCAD2/MCP
Table 6:  Best performing shrinkage methods for $n=75$.

Table 7 summarizes our recommendations for $n=150$ (based on the computations
in Subsec. 3.3). For low sparsity, LM’s were generally good. For medium
sparsity, the conclusions did not follow a strong pattern. Several methods
performed roughly equally well and tridiagonal was, as in the earlier cases,
similar to the independence cases, except here for medium sparsity. With
stronger dependence, the biggest change from the other cases was that
SCAD1/SCAD2/MCP performed best like the medium sparsity case for tridiagonal.
We also see that generally, as sparsity increases, LM’s do relatively worse.

Sparsity | .1 | .5 | .9
---|---|---|---
Ind, Li | LM | Not: RR/LL/ASCAD1 | Not: RR/LM
Ind, H | LM | SCAD1/AL/AEN | Not : RR/LM/LL/ASCAD1
Tri | LM/SCAD1 | SCAD1/SCAD2/ MCP | Not RR/LM
Toe | LM | LM/SCAD1 | SCAD1/ SCAD2/ MCP
Table 7: Best performing shrinkage methods for $n=150$.

Table 8 summarizes our recommendations for $n=500$ (based on the computations
in Subsec. 3.4) where it is safe to assume the asymptotic properties have
kicked in. Thus, we see several methods performing well regardless of whether
they have the OP. Indeed, it is very hard to chose a single best method. There
are, however, several methods that essentially never do well, namely, RR, LL,
and ASCAD1. Further LM also does relatively poorly as sparsity increases. As
such, we recommend not using these methods. Finally, we see from the two
dependence cases that dependence has a very small effect at most when the
sample size is this large. Indeed, it is the low sparsity with heavy tails
that stands out somewhat from the rest.

Sparsity | .1 | .5 | .9
---|---|---|---
Ind, Li | Not: RR/LL/ASCAD1 | Not: RR/LL/ASCAD1 | Not: RR/LM
Ind, H | LM/SCAD1/SCAD2/MCP | Not: RR/LL/ASCAD1 | Not: RR/LM
Tri | Not: RR/LL/ASCAD1 | Not: RR/LL/ASCAD1 | Not: RR/LM
Toe | Not: RR/LL/ASCAD1 | Not: RR/LL/ASCAD1 | Not: RR/LM
Table 8: Best performing of shrinkage methods for $n=500$.

We conclude with some general observations. In the $n>p$ case, we observe LM’s
often perform best except for high sparsity. Fortunately, at this sparsity
level all shrinkage methods except for RR are essentially equivalent. LM’s are
sometimes are not very good in the heavy tailed case but, again, this mainly
occurs for high sparsity as shrinkage methods set zero coefficients to zero
faster than LM does.

We observe that the higher the sparsity, the better the methods exclude
variables that are not relevant and include only the relevant variables. This
is generally true regardless of the sample size. That is, we are observing a
sort of consistency under increasing sparsity that occurs even for sample
sizes that are not asymptotic.

That said, as sample size increases, methods with the oracle property do
indeed emerge as best, if not uniquely so, regardless of the heaviness of the
tails – as long regularity conditions e.g., conditions 2, 3 and 4 in Subsec.
Condition 2 on moments, are satisfied.

##  4 Corroboration on Real Data

As a test of our recommendations in Subsec 3.5, we used the same shrinkage
methods on the data set Superconductivity presented in [17]. This data set has
81 explanatory variables of a physical or chemical nature to explain a
response $Y$ representing temperature measurements (in degrees K) for when a
compound begins to exhibit superconductivity. Initial data analyses suggested
the data were sparse, but it was unclear how sparse. [17] suggests a sparsity
level of about 90% and our techniques here confirm this in the sense that we
find, if a linear model is fit, around 90% of the coefficients will be zero
and this will be nearly best possible from a predictive standpoint. In
addition, [17] implicitly used light tails in the error term, $\epsilon$, and
did not comment on the distribution of the explanatory variables apart from
effectively taking them as independent and not requiring any special treatment
to account for spread. Accordingly, we treated these as coming from a light
tailed distribution.

Furthermore, [17] identified 20 variables of potential importance. Of those
20, we identify only seven of them being important because the variable
importance factors decreased suddenly at the eighth most important variable.
This gives 7/81 $<$ 10%, confirming this case corresponds to the high sparsity
setting. Histograms of the residuals from the full LM suggest this falls into
the light tail case as well. Thus, we compare our computed results in this
section to the recommendations for the light tailed high sparsity cases
treated in Subsec. 3.5.

In fact, the full Superconductivity data set had $n=21263$, so [17] was able
to use a standard (unpenalized) LM as a ‘benchmark model’ and then improved on
it by developing an XGBoosting model – a boosted, penalized tree model in
which the penalty was carefully constructed to be appropriate for trees.

Here, as is common in pratice, especially where a more justifiable methodology
is infeasible, we have used LM’s for their interpretability, Also, when
$n<<p$, XGBoosting often does not perform well. So, it may sometimes be
reasonable to use shrinkage techniques in mis-specified model situations with
small sample sizes.

Since Superconductivity is so much larger than the data sets used in our
simulations, we drew 40, 75, 150, and 500 data points at random so comparisons
with our recommendations here would be fair. We note that many data sets are
much smaller than Superconductivity so our example here is intended to be
suggestive for them too.

We repeated the analyses presented in Sec. 3 for the independent cases with
light tails but replaced the simulated data with the randomly chosen subsets
of Superconductivity. We were able to generate instability curves but not the
variable selection accuracy tables because the true model is unknown. The
instability curves for Superconductivity are given in Fig. 5.

|
---|---
|
Figure 5: Instability curves for the Superconductivity data for $n=40,75$
(top) and $n=150,500$ (bottom)

For each sample size, we compare the best methods from Fig. 5 to the
corresponding recommendations in Subsec. 3.5. For $n=40$, the upper left panel
in Fig. 5 shows that RR and EN are the best shrinkage methods. This is the
same as recommended in Table 5 for sparsity .9 and light, independent tails.
For $n=75$, the upper right panel in Fig. 5 shows EN is the best, followed by
RR and L which are noticeably worse. Table 6 shows only that RR should not be
used with light, independent tails. So, again we see agreement even if the
recommendation was not specific.

By contrast, for $n=150$, the lower left panel in Fig. 5 shows that EN is
best, closely followed by L and RR. Table 7 indicates that RR and LM’s are to
be avoided (for light independent tails). So,the good performance of RR
disagrees our recommendations. Finally, for $n=500$, the lower right panel in
Fig. 5 shows that five methods form a cluster of the best of the 11 methods.
The cluster of top methods is ALASSO, L, AEN, and RR. The recommendation from
Table 8 is not to use RR or LM’s. Again, we have a disgreement on the use of
RR.

We explain these findings by model mis-specification. First, the true model is
almost certainly not a LM. Indeed, [17] ends up proposing a model based on
trees. The agreement between our recommendations and the data analysis for
small values of $n$ probably means that the sample size is too small to detect
the difference between the true model and a LM. However, when the sample size
increases, the model mis-specification matters. RR normally performs well for
non-sparse cases but here is performing well when the true model is sparse.

We conjecture this occurs because a non-sparse linear model may provide a
better approximation to a sparse non-linear model than a sparse linear model
does. The analogy is to imagine representing a single true tree model with a
single linear model. The linear model would have to have many terms to
approximate a tree even one with relatively few nodes. That is, a large enough
LM might provide a good approximation.

As a final point for this section, if we were to redo our simulations using a
different model class, i.e., not LM’s, we could end up with different
recommendations and if the model class contained the true model for the
Superconductivity data we would expect our recommendations to match the data
analyes.

##  5 Optimizing Over the Shrinkage Method

Although previous sections used existing well-studied shrinkage methods, the
results of Sec. 2 show that there are infinitely many other penalties that
could be used to get shrinkage methods with the OP. Recalling that penalties
are special cases of priors, it is clear that the recommendations for choice
of shrinkage method given in Sec. 3 are limited. Here we propose that, rather
than choosing a shrinkage method from a list of options, one should find a
prior by optimizing a predictive optimality criterion using an adaptive search
technique such as a genetic algorithm (GA) or, more exactly, find a best
posterior based on a portion of the data that can be used as a prior for
remaining and forthcoming data.

We begin by extending the results of Sec. 2 to allow for data-dependent shifts
in parameter locations. The main benefit of shifting the location of the
penalty is that it reduces prior-data conflict. That way, when we find an
optimal penalty in the next subsection, it will correspond to putting priors
on the $\beta_{j}$’s that have more of their mass close to the true values of
the parameters, thereby improving inference. If the location shift is not used
in the penalty, our method below still can be used but is not as effective.
especially in small samples. It just makes sense that if the true value of a
$\beta_{j}$ is not zero, then we should not use a prior centered at zero. Then
we present our GA methodology and verify that our methodology seems to achieve
optimality in simulations. This shows that a GA approach is a viable
alternative to pre-selecting a shrinkage method.

###  5.1 Extending the Theory of Sec. 2

Let $\hat{\beta}^{*}$ be a $\sqrt{n}-$consistent estimator of $\beta$. To take
advantage of the fact that shrinkage methods can set $\hat{\beta}_{j}$’s to
zero, it is natural to choose $\hat{\beta}^{*}$ to be from a specific
shrinkage method such as SCAD2 that only requires the estimation of one extra
parameter. Adaptive methods such as ALASSO, AEN, etc., are also viable. The
idea is to use the $\hat{\beta}^{*}_{j}$’s in $\hat{\beta}^{*}$ to adjust the
location of the penalty function. This overuse of the data is standard in
shrinkage methods (in the James-Stein sense) where it is used to improve the
risk performance of decisions.

We state our extensions to Theorems 2.1 and 2.2 as corollaries since we assume
the same hypotheses. For penalized log likelihoods with location shifted
penalties we have the following sufficient conditions for the OP to hold.

######  Corollary 5.1.

Redefine the objective function in Subsec. 2.1 to be

$Q(\beta)=L(\beta|x^{n})+n\sum^{p}_{j=1}\lambda_{j}f_{j}(\beta_{j}-\hat{\beta}^{*}_{j}).$

Then, under the same conditions as in Theorem 2.1, the estimator
$\hat{\beta}=(\hat{\beta_{1}}^{\prime},\hat{\beta_{2}}^{\prime})^{\prime}$
that minimizes $Q(\beta)$ has the OP, i.e.,

$P(\hat{\beta_{2}}=0)\rightarrow 1\quad{\rm
and}\quad\sqrt{n}\hat{I}_{1}(\beta_{1}|x^{n})(\hat{\beta_{1}}-\beta_{1})\rightarrow\mathcal{N}(0,I_{1}(\beta_{1}|x^{n})).$

###### Proof.

The proof of Cor. 5.1 follows directly from the proof of Theorem 2.1. Indeed,
if a true $\beta_{j}=0$, then for large $n$ we have
$P(\hat{\beta}^{*}_{j}=0)\rightarrow 1$. So, in the proof of Lemma 2,the
inequalities at the end are asymptotically unchanged. Then, since
$\hat{\beta}^{*}_{j}$ only appears in the penalty, not the likelihood, and in
the proof of Theorem 2.1 the penalty only has to be controlled in the last
term of (7.1), to get the asymptotic normality it is enough for
$\sqrt{n}\lambda_{j}\rightarrow 0$, as guaranteed by the hypotheses. ∎

For penalized empirical risks with location shifted penalties we have the
following analog to Theorem 2.2.

######  Corollary 5.2.

Redefine the objective function in Subsec. 2.2 to be

$Q(\beta)=R(\beta|x^{n})+n\sum^{p}_{j=1}\lambda_{j}f_{j}(\beta_{j}-\hat{\beta}^{*}_{j}).$

Then, under the same conditions as Theorem 2.2, the estimator
$\hat{\beta}=(\hat{\beta_{1}}^{\prime},\hat{\beta_{2}}^{\prime})^{\prime}$
that minimizes $Q(\beta)$ has the OP, i.e.,

$P(\hat{\beta_{2}}=0)\rightarrow 1\quad{\rm
and}\quad\sqrt{n}\hat{I^{*}}_{1}(\beta_{1}|x^{n})(\hat{\beta_{1}}-\beta_{1})\rightarrow\mathcal{N}(0,I^{*}_{1}(\beta_{1}|x^{n})).$

Cor. 5.2 follows from Theorem 2.2 the same as Cor. 5.1 follows from Theorem
2.1.

The methods motivated by these corollaries continue to allow shrinkage via the
$w_{j}$’s as well as ‘James-Stein’ type shrinkage. So, we are introducing
another $p$ hyperparameters. For this reason we only recommend this ‘double
shrinkage’ approach when $n$ is not too much smaller than $p$ and preferably
$n>p$. For these cases, we continue to set
$\hat{w}_{j}=1/|\hat{\beta}_{j,OLS}|$. Taken together, double shrinkage lets
us set coefficients to zero from the OP on the $\hat{\beta}^{*}_{j}$’s and
from the OP on the $\beta_{j}$’s. Moreover, the weights $\hat{w}_{j}$ also
function to give tighter intervals around $\beta_{j}$’s that have smaller
$|\hat{\beta}_{j,OLS}|$’s. Thus, overall, we do tend to get penalties/priors
that are centered around zero when they should be and not centered around zero
when they shouldn’t be.

###  5.2 Using GAs to Find a Shrinkage Method

Our goal is to find the penalty that leads to the best predictor. To this end,
recall a GA is a computational algorithm that tries to mimic evolution to
optimize a fitness function by using analogs of mutation, crossover, and
selection. Here, we embed a GA in a two stage optimization to find an optimal
penalty/prior. We start with a ‘population’ of penalty functions and then use
gradient descent on the $Q$ from each of them to find a corresponding set of
$\hat{\beta}$’s. We evaluate a fitness function for each $\hat{\beta}$ and
select the vectors of $\beta$’s from the current population that correspond to
the top 20% (say) of fitness values as our ‘elite’ set. Now we use mutation
and crossover on the remaining penalty functions (i.e., non-elite members of
the current population). The resulting penalty functions along with the
‘elite’ set then form the next population. Iterating this process gives a
series of populations of penalty functions with a non-decreasing ‘best’
fitness value and the $\hat{\beta}$ from the ‘best’ penalty optimizes the
fitness function.

First, we define our initial class of penalty functions. Cors. 5.1 and 5.2
imply we can use any penalty function within a very large class, as long as
the regularity conditions are met. Since our goal is to find the optimal
penalty and the optimal penalty may be very complicated, we satisfy ourselves
with merely approximating it. Here we represent $f_{j}(\beta_{j})$ using
finitely many polynomials. That is, with mild abuse of notation, we set

$\displaystyle
f_{j}(\beta_{j})=\sum^{p}_{j=1}\sum^{6}_{k=1}\alpha_{k}|\beta_{j}-\hat{\beta}^{*}_{j}|^{k}.$
(5.1)

Obviously, we would get a better approximation to an optimal penalty if we
used more terms but for present purposes sixth order polynomials turned out to
be sufficient. Our initial population of penalty functions is generated from
(5.1) by selecting $M$ values of $\alpha=(\alpha_{1},\ldots,\alpha_{6})$ IID
from a Unif[0,20], say $\alpha_{m}=(\alpha_{1,m},\ldots,\alpha_{6,m})$ for
$m=1,\ldots,M$. The GA will update this initial population denoted
$A^{0}=\\{\alpha_{1}^{0},\ldots,\alpha_{M}^{0}\\}$ of size $M$ over $F$
iterations to a final population
$A^{F}=\\{\alpha_{1}^{F},\ldots,\alpha_{M}^{F}\\}$ also of size $M$ in which
we expect essentially all members to be the same. ([18] p. 75 states that the
algorithm often stops when there is little diversity in the population, as we
detected.)

We start by showing how the typical iteration from $A^{0}$ to $A^{1}$
proceeds. Assume we have data
${\cal{D}}={\cal{D}}_{n}=\\{(y_{i},x_{i})|i=1,\ldots,n\\}$ and $\dim(x_{i})=p$
and the empirical risk

$R(\beta|{\cal{D}}_{n})=\frac{1}{n}\sum^{n}_{i=1}(y_{i}-x^{\prime}_{i}\beta)^{2}.$

In view of Cor. 5.2 we seek

$\hat{\beta}_{\alpha^{0}_{m}}=\arg\inf_{\beta}\left(\frac{1}{n}\sum^{n}_{i=1}(y_{i}-x^{\prime}_{i}\beta)^{2}+\lambda\sum^{p}_{j=1}w_{j}\sum^{6}_{k=1}\alpha_{m,k}^{0}|\beta_{j}-\hat{\beta}^{*}_{j}|^{k}\right)$
(5.2)

for each $\alpha_{m}^{0}\in A^{0}$. We find suitable values of the decay
parameter $\lambda\in\mathbb{R}^{+}$ and the $\hat{\beta}_{j}^{*}$’s based on
the data as described shortly. We will use two versions of (5.2) depending on
the relative sizes of $n$ and $p$. Specifically, if $p\geq n$ or not too much
smaller than $n$, we set all $w_{j}=1$ and all $\hat{\beta}^{*}_{j}=0$. If
$p<n$, we set $w_{j}=1/|\hat{\beta}_{j,OLS}|$ as noted in Subsec. 5.1. We make
this choice because when $n<p$ typically our asymptotic results do not apply.
In the case that $p\geq n$ (5.2) reduces to

$\hat{\beta}_{\alpha^{0}_{m}}=\arg\inf_{\beta}\left(\frac{1}{n}\sum^{n}_{i=1}(y_{i}-x^{\prime}_{i}\beta)^{2}+\lambda\sum^{p}_{j=1}\sum^{6}_{k=1}\alpha_{m,k}^{0}|\beta_{j}|^{k}\right).$
(5.3)

To solve (5.2) or (5.3), we randomly split the data to estimate the various
parameters. We begin by writing
${\cal{D}}={\cal{D}}_{train}\cup{\cal{D}}_{test}$. We reserve
${\cal{D}}_{test}$ for comapring predictors after the entire GA process is
completed. Next we split the training data

${\cal{D}}_{train}={\cal{D}}_{train,\lambda}\cup{\cal{D}}_{train,\beta}\cup{\cal{D}}_{train,\alpha}.$

We use ${\cal{D}}_{train,\lambda}$ to find $\lambda$ and
${\cal{D}}_{train,\beta}$ to find the $\hat{\beta}_{j}^{*}$’s and the
$\hat{w}_{j}$’s ($n>p$). Since $\hat{\beta}_{\alpha_{m}^{0}}$ depends on
$\lambda$ we begin by searching over a list of values $\Lambda$ equally spaced
from
$\lambda_{max}=(1/n_{train,\beta})\max|Y_{train,\beta}^{T}X_{train,\beta}|$ to
$\lambda_{min}=\gamma\lambda_{max}$ for some $0<\gamma<1$. (Here,
$n_{train,\beta}=\\#{\cal{D}}_{train,\beta}$ with corresponding data indicated
by $Y_{train,\beta}$ and $X_{train,\beta}$.) For each fixed $\alpha^{0}_{m}$
and each choice of $\lambda\in\Lambda$, we find
$\hat{\beta}_{\alpha_{m}^{0},\lambda}$ from ${\cal{D}}_{train,\beta}$ and
choose the $\hat{\lambda}^{0}_{m}$ that minimizes
$R(\hat{\beta}_{\alpha_{m}^{0},\lambda}|{\cal{D}}_{train,\lambda})$.

We find $\hat{\beta}_{\alpha^{0}_{m}}$ for each $m$ in (5.2) by sub-gradient
descent since $\alpha^{0}_{m}$, $\lambda=\hat{\lambda}^{0}_{m}$,
$w_{j}=\hat{w}_{j}$ and $\hat{\beta}^{*}_{j}$ can be taken as given. (The
$\hat{w}_{j}$ and $\hat{\beta}^{*}_{j}$ should also have sub- and super-
scripts $m$ and 0; we omit these for convenience.) Recall, the sub-gradient
descent algorithm allows for us to have points of non-differentiability in the
penalty (e.g., a corner as in L or SCAD2), and in cases where the penalty is
differentiable, the sub-gradient is uniquely defined by the gradient. Note
that the objective function is constructed to be convex, so we do indeed have
a minimum. We initialize the gradient descent algorithm at the LASSO solution
for $n>p$ and at the RR solution for $n<p$.

Now define the fitness function for the GA to be

$\displaystyle
f=\sum_{i\in{\cal{D}}_{train,\alpha}}(y_{i}-x^{\prime}_{i}\hat{\beta}_{\alpha^{0}_{m},\hat{\lambda}^{0}_{m}})^{2}.$
(5.4)

We evaluate the fitness for each $\alpha_{m}^{0}$ in $A^{0}$. Note that for
each $\alpha_{m}^{0}$ for $m=1,\ldots,M$ we get a single best choice for
$\hat{\lambda}_{m}^{0}$ and
$\hat{\beta}_{\alpha_{m}^{0},\hat{\lambda}_{m}^{0}}$ and hence a single
fitness value. However, it is possible for different $\alpha_{m}^{0}$’s to
give exactly the same $f$-value because its possible
$\hat{\beta}_{\alpha_{i}^{0},\hat{\lambda}_{m}^{0}}=\hat{\beta}_{\alpha_{j}^{0},\hat{\lambda}_{m}^{0}}$
for some $i\neq j$. Although this would appear to happen with probability
zero, it is observed on a regular basis. This arises because different but
similar penalties may lead to the same solution and because computing only has
limited precision.

Next, by elitism we select off the top 20% of members of $A^{0}$. We fill in
the ‘missing’ 80% by applying crossover and mutation to the bottom 80% of
fitness values to obtain a new generation of size $M$ from the algorithm to go
into the second iteration. Crossing means switching some entries of a genome
$\alpha^{\prime}_{m}$ with entries from another $\alpha^{\dagger}_{m}$ to
generate a ‘new genome’. This is done at random keeping only the ‘child’ until
the population size $M$ is achieved. Mutation means adding a perturbation to
all members of $\alpha$ (here a random number between the user specified
maximum and minimum values for each component in $\alpha$). Mutation does not
change the size of the population, only the specific genomes already in it. In
this way we get a new population $A^{1}$ to which we can apply the same
procedure. Then, we can iterate to get $A^{2}$, $A^{3}$ and so on until
$A^{F}$ contains little diversity on the population.

To see that this is the typical behavior of this sort of GA, we use the
framework of [19]. First, it is easy to see that as we have set it up here,
the GA is a Markov process. That is, the probabilistic behavior in moving from
time $t$ to time $t+1$ depends only on the state at time $t$. Moreover, this
Markov process is homogeneous in the sense that the transition from time step
to time step is the same for any two adjacent time steps. Note that the Markov
process is ‘discrete time’ and has a discrete population (leading to distinct
crosses) but the mutation is continuous because of the uniform distribution.
Thus, there is no transition matrix. Instead, there is a transition kernel,
$K(x,S)$, where $x$ is a population member at time $t$ and $S$ is a set of
possible states to which $x$ may be transformed and $K$ is independent of $t$.
In fact, $K(\cdot,\cdot)$ can be partitioned into a $K_{m}$ and $K_{c}$, a
mutation and crossover kernel. The crossover kernel is a transition matrix
since crossover is discrete. The mutation kernel includes the continuous
mutation phase based on the uniform distribution. So, let $x$ be any state at
time $t$ and suppose an optimum $f^{*}$ exists and the Markov process has
state space $E$. Then, there will be elements of $E$ arbitrarily close to
$f^{*}$. Let $b(x_{t})$ be the best fitness value within the $t$-th population
and let $d(x)=b(x)-f^{*}$. As long as the population is large enough,
$B_{\epsilon}=\\{d(x)<\epsilon\\}$ will have nonzero probability for
$\epsilon>0$ and hence $K_{m}(x,B_{\epsilon})$ will be bounded away from zero.
Now, given that we have used elitism, Theorem 2 in [19] applies to give
convergence of the GA to the global minimum of $f$ within the class of priors
that have the OP as in Cors. 5.1 and 5.2.

The behavior of the GA depends on $M$, the elements of $A^{0}$, the size of
$F$, the choice of $f$, the data, etc. Indeed, a pragmatic check on the
behavior of a GA would be to run it with different initial populations to see
if the GA outputs approximately the same minimum. To ensure convergence of the
GA one should set a large population as well as a large number of generations.
We comment that in the sub-gradient descent phase of our procedure, we have
limited ourselves to convex objective functions. For more general results we
would have to ensure convergence of the gradient based optimization to ensure
convergence of the GA-based optimization. We implemented our GA computations
using genalg, see [20].

Our intuition tells us that this method will be beneficial in low to medium
sparsity cases, as well as heavy dependence or non-asymptotic cases. This is
due to the fact that asymptotically the OP methods are equivalent, and thus a
GA can do no better. Further, in the smaller sample cases with high sparsity
we observe most methods performing roughly the same. In low to medium sparsity
cases, there is more variability between the methods and thus, we should be
able to optimize to find a penalty that in fact does perform better than the
common shrinkage methods. We acknowledge that only using standard basis
expansions may not allow us to approximate some penalties well.

###  5.3 Simulations

Here we present two simulations, one for $p>n$ and one for $n>p$ to show how
implementing the GA performs relative to other shrinkage methods in a
predictive setting. We simulate IID observations from

$Y=X\beta+\epsilon$

where $X\sim MVN_{p}(0,I_{100})$,
$\epsilon=(\epsilon_{1},\ldots,\epsilon_{n})^{T}$ with $\epsilon_{i}\sim
N(0,1)$, and $\beta=(\beta_{1},\beta_{2})^{T}$ and we set $p=100$, as before.
We assume 50% sparsity, so the dimension of both $\beta_{1}$ and $\beta_{2}$
are 50. We take $\beta_{1}\sim MVN_{50}(4,I_{50})$ and set $\beta_{2}=0$. We
consider $n=40$ and $n=150$ and we split the data as described in Sec. 5.2.

The GA will find an optimal penalty as defined by an optimal vector
$\alpha_{opt}=(\alpha_{1,opt},\ldots,\alpha_{6,opt})^{T}$. The entries
$\alpha_{j,opt}=\alpha_{j,opt}({\cal{D}}_{train,\alpha})$ so we are treating
the penalty as a hyperparameter in the prior that would be mathematically
equivalent to it. The difference from actually estimating a hyperparameter
comes from the fact we are only using
${\cal{D}}_{train,\alpha}\subsetneqq{\cal{D}}_{train}$. Given the penalty, we
have a potentially new shrinkage method, dependent on a proper subset of
${\cal{D}}_{train}$. So, we can form a posterior using the prior determined
from the penalty given the rest of the data. This posterior can be used to
generate predictions for ${\cal{D}}_{test}$ that can be compared with the
predictions from the other shrinkage methods used in Sec. 3.

####  5.3.1 GA example $n=40$

Here we split the data so that
$\\#({\cal{D}}_{train})=\\#({\cal{D}}_{train,\lambda}\cup{\cal{D}}_{train,\beta}\cup{\cal{D}}_{train,\alpha})=36$
with corresponding sample sizes $(2,30,4)$ and $\\#({\cal{D}}_{test})=4$. All
other methods we compare use all of the training data to find estimates of
$\beta$ and $\lambda$. Note for comparisons with other methods, those that use
glmnet use all 36 observations in the training data to form the predictor and
for the methods that we implemented with LLA we combined
${\cal{D}}_{train,\lambda}$ and ${\cal{D}}_{train,\alpha}$ to estimate
$\lambda$. Thus, we ensured that each method used all the training data,
providing a fair comparison.

Interestingly, but perhaps not surprisingly, we find
$\hat{\alpha}=(0,1,0,0,0,0)$ which corresponds exactly to RR. This is
consistent with the methods that performed the best in the analogous cases in
Subsec. 3.1. This suggests that when we have few data points relative to
explanatory variables, we do not have enough information to obtain an
informative prior (in terms of its location and variance) so we default to the
prior that makes us retain all the explanatory variables.

The predictive errors for $\\#({\cal{D}}_{test})$ are given in Table 9. Note
that GA, RR,EN, and AEN are the same here, because we set $w_{j}=1,\forall j$.
We comment that because GA’s require a lot of computing time, we have not
averaged over many data sets to get the prediction errors reported in this
table. However, we believe we have used a large enough population and large
enough number of generations that our results are accurate.

GA | L | RR | EN | AEN | AL | SCAD2 | MCP
---|---|---|---|---|---|---|---
22.75 | 45.37 | 22.75 | 22.75 | 22.75 | 46.62 | 36.18 | 36.18

Table 9: MSPE for our new GA method and seven other methods.

This example illustrates that by optimizing over the choice of penalties, we
are not guaranteed to find a penalty that is different from an established
method (although we argue this is the typical case). The guarantee is only
that we will find an optimal penalty for prediction and it is no surprise if
there are settings where a well known technique is optimal. The novelty in our
GA approach is that it can be used in any linear regression problem and, if
properly implemented, will always give the best predictions.

####  5.3.2 GA example $n=150$

When splitting the data in this situation, we must keep more than 100
observations in ${\cal{D}}_{train,\beta}$ to ensure $n>p$. Accordingly, we set

$\\#({\cal{D}}_{train})=\\#({\cal{D}}_{train,\lambda}\cup{\cal{D}}_{train,\beta}\cup{\cal{D}}_{train,\alpha})=135$

with corresponding sample sizes of $(9,113,13)$ respectively, and
$\\#({\cal{D}}_{test})=15$. Parallel to our methodology in Subsec. 5.3.1, the
methods implemented using glmnet and rqPen used all 135 observations in the
training data to form the predictor. Also, as before, for the methods that we
implemented with LLA we combined ${\cal{D}}_{train,\lambda}$ and
${\cal{D}}_{train,\alpha}$ to estimate $\lambda$. Again, this ensured all
methods were being treated fairly.

After running the GA we found $\hat{\alpha}=(2,0,0,20,0,0)$. The associated
prediction error on ${\cal{D}}_{test}$ for each method is given in Table 10.
We observe the penalty selected through GA achieves the best predictive error
among all methods considered. As in Subsec. 5.3.1, we comment that because
GA’s require a lot of computing time, we have not averaged over many data sets
to get the prediction errors reported in this table. However, we believe we
have used a large enough population and large enough number of generations
that our results are accurate.

Since we found a new (and better) penalty, we have graphed it in Fig. 6.
Recall that half the parameter values are zero, half are non-zero, and the
penalty term has $f_{j}$’s, i.e., the terms depend on the index of the
parameter. This means that we allow different penalties on different
parameters. The left hand panel shows a plot of the optimal penalty for one of
the $\beta_{j}$’s that is known to be zero. It is compared with the common
penalties RR, LASSO, and SCAD. The right hand panel shows a plot of the
optimal penalty for one of the $\beta_{j}$’s that is known to be non-zero,
again compared with RR, LASSO, and SCAD. It is obvious that the GA method
described in Subsec. 5.1 gives two sorts of $f_{j}$’s. The training data
forces the $f_{j}$’s corresponding to $\beta_{j}=0$ to concentrate at zero and
forces the $f_{j}$’s that correspond to nonzero $\beta_{j}$’s to concentrate
away from zero. This explains the improvement in prediction error seen in
Table 10.

GA | LM | L | RR | EN | AEN | AL | ASCAD1 | SCAD1 | SCAD2 | LL
---|---|---|---|---|---|---|---|---|---|---
1.12 | 1.52 | 1.78 | 3.66 | 1.78 | 1.50 | 1.50 | 1.36 | 1.53 | 1.18 | 1.67
Table 10: MSPE for our new GA method and for 10 other methods.

To end this section, note that our simulations only show proof of concept; the
priors we found here may not be genuinely optimal for prediction. That is so
because we have not run the GA for many generations with a large population
size so we cannot assume the GA has converged. In fact, in both cases here
($n=40,150$) we only ran a single generation of the GA and we only used a
population size of $M=150$. However, because of the elitism operation, running
the GA longer can never result in a worse predictor and our results show that
it can be relatively easy to find a penalty that is better for prediction than
established penalties – even if they are not optimal within the class of all
penalties with the OP.

---
Figure 6: GA Optimal Penalty vs Standard Penalties for Zero and Nonzero
$\beta_{j}$’s.

##  6 Discussion

This paper has assumed a predictive stability perspective and within that
context shown several results that may be a bit unexpected. First, the OP is
not rare; it is actually rather common. Its proof requires little more than
what most would regard as regularity conditions. Second, for small $n$ and
large $p$ shrinkage methods did not perform very well even if they have the OP
and the true model has reasonable sparsity. Third, on the other hand, if
optimal or near optimal penalties are used they give shrinkage methods that
work noticeably better than the established ones. Our findings indicate that
methods having the OP do not perform particularly well for $n<p$ and for $n>p$
the OP is no guarantee that they perform better than methods not having the
OP. Fourth, our results also indicate that with increasing sparsity the
performance of shrinkage methods improves. This intuition needs to be
developed further because the obvious limit of perfect sparsity gives the
trivial model.

So, even though the OP is important, it is not at all clear how important it
is or when it is important in cases where sample sizes are finite. We still
think it’s better to have the OP than not if only because it gives
consistency, asymptotic normality and efficiency. This is especially the case
with high sparsity and large $n$ relative to $p$, but in these cases other
methods often perform comparably. When $n$ is small compared to $p$, the OP is
not a useful property, and thus the adaptive methods that have the OP do not
perform well. A possible explanation for this is a poor bias variance trade
off when $p>n$. It does not seem to be a good idea to use methods that require
estimating $w_{j}$ for each $\beta_{j}$: For $p>n$ we have not seen any
example where the adaptive penalty gives better results that the nonadaptive
version.

Since the OP really requires $n\rightarrow\infty$ whereas $n$ often must be
taken as truly finite, we introduce the notion of instability of predictions
as a criterion for selecting a penalty or prior. Comparing instability curves
is a finite sample check for good predictive performance. Using this approach
we can easily rule out unstable predictors. For instance, with high sparsity
using a linear model by itself is often unstable. In general, quantifying the
variability of variable selection when $p$ is large is difficult, so defining
instability in terms of the prediction errors seems reasonable.

Our simulation studies show that as a generality, shrinkage methods tend to
perform better in terms of variable selection, and thus prediction, as
sparsity increases as well as when $n$ increases. In fact, our simulations
showed that regardless of $n$, as the sparsity increased, the methods seemed
to perform roughly equally well. For instance, recall the $n=75$ simulations
in Subsec. 3.2 . At $90\%$ sparsity, we observed what appeared to be
asymptotic convergence of the method with the OP. Of course this situation is
not asymptotic as $n<p$, but an increase in sparsity is associated with an
increase in efficiency of the methods.

Since there are infinitely many such choices for penalties that have the OP,
we take the subjectivity out of penalty (or prior) selection by using a GA to
find an optimal penalty/prior for prediction. When $n>p$ we use the GA
approach to find a predictively optimal penalty that has the OP. When $n<<p$,
we do not search for methods with the OP because we do not benefit from the
known asymptotic results. Thus, we search over the class of penalties that are
non-adaptive and do not require estimation of many hyperparameters. In
principle, as long as we let the GA run long enough to converge, this approach
can never do worse that simply choosing a standard shrinkage method. In fact,
we have examples where the GA approach does better than others; when the GA
approach selects the best among standard methods, we can infer the standard
method was the right choice.

Another way to look at the procedure in Sec. 5 is that when we find a
penalty/prior based on the data we are producing an approximation to a
predictively optimal posterior given the training data that can then be used
with the log-likelihood. Thus, the predictive improvement comes from the
efficiency of the way the posterior uses the data with an optimal prior. That
is, we are using the data to form a pseudo-posterior based on a de facto
optimal prior only determined by the (pseudo-)posterior it forms.

We close with another heuristic that seems to be borne out by our results.
Namely, we associate corners and other points of non-differentiability in the
penalty with setting parameter values equal to zero in finite samples. Recall
that minimizing

$\sum^{n}_{i=1}L(y_{i}-x^{\prime}_{i}\beta)+n\lambda\sum^{p}_{j=1}w_{j}f_{j}(\beta_{j})$

is equivalent to minimizing $Q=\sum^{n}_{i=1}L(y_{i}-x^{\prime}_{i}\beta)$
subject to the constraint $\sum^{p}_{j=1}w_{j}f_{j}(\beta_{j})\leq
R\in\mathbb{R}$ where $R$ typically decreases as $\lambda$ increases. Denote
the constraint region by

$D=\left\\{\sum^{p}_{j=1}w_{j}f_{j}(\beta_{j})\leq R\right\\}.$

Since $D$ is closed and compact, the Krein-Milman theorem (see [21]) gives
that $D$ is the closed convex hull of its extreme points, i.e.,
$D=CCH(D_{ext})$. For reasonable choices of $f_{j}$, $D$ is defined by the
intersection of regions of the form

$U_{k}(\beta_{1},\ldots,\beta_{p})\leq\alpha_{k},k=1,\ldots,p,$

where the $U_{k}$’s are defined from the $f_{j}$’s and $w_{j}$’s. When our
goal is optimizing $Q$ over $D$, it is often the case that the optima occur at
extreme points of $D$. When the extreme points of $D$ are on the coordinate
axes we will find at least some of $\beta_{j}$’s are zero. Indeed, if $Q$ is
convex and continuous on an open set containing $D$, then $Q$ often attains
its minimum over $D$ on a ‘face ’of $D$ and the exact point where the optima
occur may lie at the intersections of some or all of the $U_{k}=\alpha_{k}$;
this defines a subset of the extreme points of $D$. This is well-established
for the case of linear optimization with linear constraints. Indeed, if $Q$ is
minimized for at least one extreme point of $D$ that lies on a coordinate axis
then at least some $\beta_{j}$’s will be set to zero. This means that any
locally convex penalty with a ‘corner’ on a coordinate axis will perform
nontrivial variable selection if $R$ is small enough. We conjecture a converse
to this statement will hold, too.

##  7  Supplement : Proofs from Sec. 2

###  7.1 Penalized Log-Likelihood Case, Sec. 2.1

######  Theorem 7.3.

Suppose Conditions 1–6 are satisfied and suppose $\sqrt{n}a_{n}\rightarrow 0$
where $a_{n}$ satisfies $a_{n}=\frac{1}{h(n)\sqrt{n}}$ in which
$\frac{1}{h(n)\sqrt{n}}\rightarrow 0$ and $h(n)\rightarrow\infty$. Then, there
exists a local minimizer $\hat{\beta}$ of $Q(\beta)$ such that
$||\hat{\beta}-\beta_{0}||=O_{p}(n^{-\frac{1}{2}}+a_{n})$.

###### Proof.

Step 1. We want to show for any $\epsilon>0$ there exists a large constant $C$
such that

$P\left\\{\inf_{||u||=C}Q(\beta_{0}+\alpha_{n}u)>Q(\beta_{0})\right\\}\geq
1-\epsilon,$ (7.1)

where $\alpha_{n}=n^{-\frac{1}{2}}+a_{n}$. Denote

$\displaystyle D_{n}(u)$ $\displaystyle=Q(\beta_{0}+\alpha_{n}u)-Q(\beta_{0})$
$\displaystyle=L(\beta_{0}+\alpha_{n}u|x^{n})+\sum_{j=1}^{p}\lambda_{j}f_{j}(\beta_{j0}+\alpha_{n}u_{j})-L(\beta_{0}|x^{n})-n\sum_{j=1}^{p}\lambda_{j}(f_{j}(\beta_{j0}))$
$\displaystyle=L(\beta_{0}+\alpha_{n}u|x^{n})-L(\beta_{0}|x^{n})+n\sum_{j=1}^{p}\lambda_{j}[f_{j}(\beta_{j0}+\alpha_{n}u_{j})-f_{j}(\beta_{j0})]$
$\displaystyle=L(\beta_{0}+\alpha_{n}u|x^{n})-L(\beta_{0}|x^{n})+n\sum_{j=1}^{p}\lambda_{j}f_{j}^{\prime}(\tilde{\beta_{j}})\alpha_{n}u_{j}$
$\displaystyle=L(\beta_{0}+\alpha_{n}u|x^{n})-L(\beta_{0}|x^{n})+n\sum_{j=1}^{p_{0}}\lambda_{j}f_{j}^{\prime}(\tilde{\beta_{j}})\alpha_{n}u_{j}+n\sum_{j=p_{0}+1}^{p}\lambda_{j}f_{j}^{\prime}(\tilde{\beta_{j}})\alpha_{n}u_{j}$

where $f_{j}^{\prime}(\tilde{\beta_{j}})\alpha_{n}u_{j}$ is the Taylor
expansion of the penalty term by Conditions 4 and 6 and where
$\tilde{\beta_{j}}$ is on the line joining $\beta_{j}$ to
$\beta_{j}+\alpha_{n}u$. More formally,
$\tilde{\beta_{j}}\in<\beta_{j},\beta_{j}+\alpha_{n}u>$. Thus we have,

$\displaystyle D_{n}(u)$
$\displaystyle=L(\beta_{0}+\alpha_{n}u|x^{n})-L(\beta_{0}|x^{n})+n\alpha_{n}\left(\sum_{j=1}^{p_{0}}\lambda_{j}f_{j}^{\prime}(\tilde{\beta_{j}})u_{j}+n\sum_{j=p_{0}+1}^{p}\lambda_{j}f_{j}^{\prime}(\tilde{\beta_{j}})u_{j}\right)$
$\displaystyle=L(\beta_{0}+\alpha_{n}u|x^{n})-L(\beta_{0}|x^{n})+n(n^{-\frac{1}{2}}+a_{n})\left(\sum_{j=1}^{p_{0}}\lambda_{j}f_{j}^{\prime}(\tilde{\beta_{j}})u_{j}+n\sum_{j=p_{0}+1}^{p}\lambda_{j}f_{j}^{\prime}(\tilde{\beta_{j}})u_{j}\right)$
$\displaystyle=L(\beta_{0}+\alpha_{n}u|x^{n})-L(\beta_{0}|x^{n})+(\sqrt{n}+na_{n})\left(\sum_{j=1}^{p_{0}}\lambda_{j}f_{j}^{\prime}(\tilde{\beta_{j}})u_{j}+n\sum_{j=p_{0}+1}^{p}\lambda_{j}f_{j}^{\prime}(\tilde{\beta_{j}})u_{j}\right)$
$\displaystyle\geq
L(\beta_{0}+\alpha_{n}u|x^{n})-L(\beta_{0}|x^{n})+(\sqrt{n}+na_{n})a_{n}\left(\sum_{j=1}^{p_{0}}f_{j}^{\prime}(\tilde{\beta_{j}})u_{j}+n\sum_{j=p_{0}+1}^{p}f_{j}^{\prime}(\tilde{\beta_{j}})u_{j}\right)$
$\displaystyle=L(\beta_{0}+\alpha_{n}u|x^{n})-L(\beta_{0}|x^{n})+(\sqrt{n}a_{n}+na_{n}^{2})\left(\sum_{j=1}^{p_{0}}f_{j}^{\prime}(\tilde{\beta_{j}})u_{j}+n\sum_{j=p_{0}+1}^{p}f_{j}^{\prime}(\tilde{\beta_{j}})u_{j}\right)$
$\displaystyle=L(\beta_{0}+\alpha_{n}u|x^{n})-L(\beta_{0}|x^{n})+\sqrt{n}a_{n}(1+\sqrt{n}a_{n})\left(\sum_{j=1}^{p_{0}}f_{j}^{\prime}(\tilde{\beta_{j}})u_{j}+n\sum_{j=p_{0}+1}^{p}f_{j}^{\prime}(\tilde{\beta_{j}})u_{j}\right)$
(7.2)

where $p_{0}$ is the number of components in $\beta_{10}$.

Now by Taylor expansion of the log likelihood at $\beta_{0}$,

$\displaystyle L(\beta_{0})+\alpha_{n}u)-L(\beta_{0}|x^{n})$
$\displaystyle=\alpha_{n}uL^{\prime}(\beta_{0}|x^{n})+\frac{n}{2}(\alpha_{n}u)^{\prime}\left(I(\beta_{0}|x^{n})\right)(\alpha_{n}u)$
$\displaystyle+\frac{n}{2}(\alpha_{n}u)^{\prime}\left(\frac{L^{\prime\prime}(\tilde{\beta}|x^{n}}{n}-I(\beta_{0}|x^{n})\pm\frac{1}{n}\sum^{n}_{i=1}I(\beta_{0}|x^{n})\pm
I(\beta_{0})\right)(\alpha_{n}u)$
$\displaystyle=\alpha_{n}L^{\prime}(\beta_{0}|x^{n})u-\frac{n}{2}u^{\prime}I_{n}(\beta_{0})u\alpha^{2}_{n}\\{1+o_{p}(1)\\}$
(7.3)

where $I_{n}(\beta_{0})=\frac{1}{n}\sum^{n}_{i=1}I(\beta_{0}|x^{n})$ by
Condition 4 and we have

$E\left[\sup_{\beta\in
B(\beta_{0},\eta)}\left|\frac{L^{\prime\prime}(\tilde{\beta}|x^{n})}{n}-I(\beta_{0}|x^{n})\right|\right]\rightarrow
0\text{ as }\eta\rightarrow 0$

which implies that $I_{n}(\beta_{0})\rightarrow I(\beta_{0})$ by Condition 3.
Using $(\ref{dnu})$ and $(\ref{lnu})$ ,

$\displaystyle D_{n}(u)$ $\displaystyle\geq$
$\displaystyle\alpha_{n}L^{\prime}(\beta_{0}|x^{n})^{\prime}u+\frac{n}{2}u^{\prime}I(\beta_{0}|x^{n})u\alpha^{2}_{n}\\{1+o_{p}(1)\\}+\sqrt{n}a_{n}(1+\sqrt{n}a_{n})\sum_{j=1}^{p_{0}}f_{j}^{\prime}(\tilde{\beta_{j}})u_{j}$
(7.4)
$\displaystyle+\sqrt{n}a_{n}(1+\sqrt{n}a_{n})\sum_{j=p_{0}+1}^{p}f_{j}^{\prime}(\tilde{\beta_{j}})u_{j}.$

We argue the second term on the RHS of (7.4) is dominant. To see this consider
the first term on the RHS of (7.4) using Conditions 1 and 3 We multiply by
$\frac{\sqrt{n}}{\sqrt{n}}$ and obtain
$\sqrt{n}\alpha_{n}\frac{L^{\prime}(\beta_{0}|x^{n})}{\sqrt{n}}u$. Let
$\bar{L}=\frac{L^{\prime}(\beta_{0}|x^{n})}{\sqrt{n}}$, and observe
$\sqrt{n}\alpha_{n}\rightarrow 1$. So, $C^{2}=||u||^{2}=u^{\prime}u>\bar{L}u$
as long as $0<\frac{1}{B}<\bar{L}<B$ with high probability for $C>B$. Then
with high probability $C^{2}\geq B\underline{1}(u)\geq B\bar{L}$ for $||u||=C$
and $\underline{1}=(1,1,\ldots,1)$. Thus we have that by choosing sufficiently
large $C$, the second term of (7.4) is larger in absolute value, with high
probability than the first term uniformly in $||u||=C$.

Also, by hypothesis $\sqrt{n}a_{n}\rightarrow 0$, so the whole third term in
(7.4) goes to zero and the second term of (7.4) is also larger than the third
term. The fourth term of (7.4) also goes to zero because for $p_{0}\leq j\leq
p$, we know that $\beta_{j}=0$, so $\tilde{\beta_{j}}\rightarrow 0$ as
$n\rightarrow 0$.Therefore
$\sum_{j=p_{0}+1}^{p}f_{j}^{\prime}(\tilde{\beta_{j}})u_{j}\rightarrow 0$ by
Condition 5 and since $\sqrt{n}a_{n}\rightarrow 0$, the fourth term on the RHS
of (7.4) goes to 0. Thus since the third term third and fourths term on the
RHS of (7.4) both go to 0, we have that the second term on the RHS of (7.4) is
larger than all of the other terms.

Therefore, by choosing large enough $C$,

$P\left\\{\inf_{||u||=C}D_{n}(u)\geq\frac{n}{2}u^{\prime}I(\beta_{0}|x^{n})u\alpha^{2}_{n}\\{1+o_{p}(1)\\}\right\\}\geq
1-\epsilon.$

Thus, (7.1) is true for sufficiently large $C$.

Step 2. We now show the conclusion of Step 1 holds for any $C^{*}>C$. Let
$C^{\star}>C$. Then,

$\displaystyle\hat{\beta}$
$\displaystyle\in\\{\beta_{0}+\alpha_{n}u:||u||=C^{*}\\}$
$\displaystyle\iff\hat{\beta}\in\cal{B}\left(\beta_{0},\alpha_{n}C^{*}\right)$
$\displaystyle\iff\hat{\beta}-\beta_{0}\in B\left(0,\alpha_{n}C^{*}\right)$
$\displaystyle\iff\alpha_{n}^{-1}(\hat{\beta}-\beta_{0})\in B(0,C^{*})$
$\displaystyle\iff||\alpha_{n}^{-1}(\hat{\beta}-\beta_{0})||<C^{*}$
$\displaystyle\iff||\hat{\beta}-\beta_{0}||^{2}<\alpha_{n}^{2}{C^{*}}^{2}$

It follows that
$\sum^{p}_{j=1}(\hat{\beta}_{j}-\beta_{0j})^{2}\leq\alpha_{n}^{2}{C^{*}}^{2}$
and $\forall j$ ,
$(\hat{\beta}_{j}-\beta_{0j})^{2}\leq\alpha_{n}^{2}{C^{*}}^{2}$. So,

$-\alpha_{n}C^{*}\leq\hat{\beta}_{j}-\beta_{0j}\leq\alpha_{n}C^{*}$

and we can absorb $C^{*}$ into $\alpha_{n}$ because $C^{*}$ is a constant.
Thus we have $\hat{\beta}-\beta_{0}=O_{p}(\alpha_{n})$. Note that

$\alpha_{n}=\frac{1}{\sqrt{n}}+a_{n}=\frac{1+\sqrt{n}a_{n}}{\sqrt{n}}$

and as $n\rightarrow\infty$ we have $\sqrt{n}a_{n}\rightarrow 0$, so
$\hat{\beta}-\beta_{0}=O_{p}\left(\frac{1}{\sqrt{n}}\right)$. ∎

Remark: The argument for Theorem 1 is true for differentiable $f_{j}(\cdot)$,
but not non-differentiable $f_{j}(\cdot)$. In the case of non-differentiable
$f_{j}(\cdot)$, consider a smooth function that approximates $f_{j}(\cdot)$
well. Say, $f_{j}^{\star}(\cdot)$, which differs from $f_{j}(\cdot)$ by a
margin of $\epsilon$ where $\epsilon_{N}\rightarrow 0$ as
$n\rightarrow\infty$. As long as
$\lim_{\beta\rightarrow\beta^{*-}}f_{j}^{\star}(\cdot)=\lim_{\beta\rightarrow\beta^{*+}}f_{j}^{\star}(\cdot)$,
then the above argument holds.

######  Lemma 2.

Assume conditions 1- 6, and the result from Theorem 7.3 holds. If
$\sqrt{n}b_{n}\rightarrow\infty$ and $b_{n}\geq\frac{g(n)}{\sqrt{n}}$ as
$n,g(n)\rightarrow\infty$, then with probability tending to 1, for any given
$\beta_{1}$ satisfying $||\beta_{1}-\beta_{10}||=O_{p}(n^{-\frac{1}{2}})$ and
any constant $C$,

$Q\left\\{\left(\begin{matrix}\beta_{1}\\\
0\end{matrix}\right)\right\\}=\min_{||\beta_{2}||\leq
Cn^{-\frac{1}{2}}}Q\left\\{\left(\begin{matrix}\beta_{1}\\\
\beta_{2}\end{matrix}\right)\right\\}.$ (7.5)

###### Proof.

Consider the objective function
$Q(\beta)=L(\beta|x^{n})+n\sum^{p}_{j=1}\lambda_{j}f_{j}(\beta_{j})$. Note
that

$\frac{\partial Q(\beta)}{\partial\beta_{j}}=\frac{\partial
L(\beta|x^{n})}{\partial\beta_{j}}+n\lambda_{j}f_{j}^{\prime}(\beta_{j}).$
(7.6)

Then by Taylor expanding at $\beta_{0}=\left(\begin{matrix}\beta_{10}\\\
0\end{matrix}\right)$ we have

$\displaystyle\frac{\partial Q(\beta)}{\partial\beta_{j}}$
$\displaystyle=\frac{\partial
L(\beta_{0}|x^{n})}{\partial\beta_{j}}+\sum^{p}_{\ell=1}\frac{\partial^{2}L(\beta_{0}|x^{n})}{\partial\beta_{j}\partial\beta_{\ell}}(\beta_{\ell}-\beta_{\ell
0})+\sum^{p}_{\ell=1}\sum^{p}_{k=1}\frac{\partial^{3}L(\beta^{*}|x^{n})}{\partial\beta_{j}\partial\beta_{\ell}\partial\beta_{k}}(\beta_{\ell}-\beta_{\ell
0})(\beta_{k}-\beta_{k0})$
$\displaystyle+n\lambda_{j}f_{j}^{\prime}(\beta_{j})$

where $\beta^{*}$ lies between $\hat{\beta}$ and $\beta_{0}$. Note that

$\frac{1}{n}\frac{\partial
L(\beta_{0}|x^{n})}{\partial\beta_{j}}=O_{p}(n^{-\frac{1}{2}})$

and due to [22] and by the Law of Large Numbers for non-identically
distributed random variables we have

$\frac{1}{n}\frac{\partial^{2}L(\beta_{0}|x^{n})}{\partial\beta_{j}\partial\beta_{\ell}}=E\left[\frac{\partial^{2}L(\beta_{0}|x^{n})}{\partial\beta_{j}\partial\beta_{\ell}}\right]+o_{p}(1).$

Note that due to Conditions 1–6, and the result in Theorem 1,
$\hat{\beta}-\beta_{0}=O_{p}\left(\frac{1}{\sqrt{n}}\right)$. So, we have

$\displaystyle\frac{\partial Q(\beta)}{\partial\beta_{j}}$
$\displaystyle=n\frac{1}{n}\frac{\partial
L(\beta_{0}|x^{n})}{\partial\beta_{j}}+n\frac{1}{n}\sum^{p}_{\ell=1}\frac{\partial^{2}L(\beta_{0}|x^{n})}{\partial\beta_{j}\partial\beta_{\ell}}(\hat{\beta_{\ell}}-\beta_{\ell
0})$
$\displaystyle+n\frac{1}{n}\sum^{p}_{\ell=1}\sum^{p}_{k=1}\frac{\partial^{3}L(\beta^{*}|x^{n})}{\partial\beta_{j}\partial\beta_{\ell}\partial\beta_{k}}(\hat{\beta_{\ell}}-\beta_{\ell
0})(\hat{\beta_{k}}-\beta_{k0})+n\lambda_{j}f_{j}^{\prime}(\beta_{j})$
$\displaystyle=nO_{p}\left(\frac{1}{\sqrt{n}}\right)+n\sum^{p}_{\ell=1}\left(E\left[\frac{\partial^{2}L(\beta_{0}|x^{n})}{\partial\beta_{j}\partial\beta_{\ell}}\right]+o_{p}(1)\right)O_{p}\left(\frac{1}{\sqrt{n}}\right)$
$\displaystyle+n\sum^{p}_{\ell=1}\sum^{p}_{k=1}E\left[\sup_{\beta\in
N}\left|\frac{\partial^{3}L(\beta^{*}|x^{n})}{\partial\beta_{j}\partial\beta_{\ell}\partial\beta_{k}}\right|+o_{p}(1)\right]O_{p}\left(\frac{1}{\sqrt{n}}\right)O_{p}\left(\frac{1}{\sqrt{n}}\right)+n\lambda_{j}f_{j}^{\prime}(\beta_{j})$
$\displaystyle=n\left[O_{p}\left(\frac{1}{\sqrt{n}}\right)+O_{p}\left(\frac{1}{\sqrt{n}}\right)+O_{p}\left(\frac{1}{\sqrt{n}}\right)\right]+n\lambda_{j}f_{j}^{\prime}(\beta_{j})$
$\displaystyle=n\left[O_{p}\left(\frac{1}{\sqrt{n}}\right)+\lambda_{j}f_{j}^{\prime}(\beta_{j})\right]$
$\displaystyle=\sqrt{n}\left[O_{p}(1)+\sqrt{n}\lambda_{j}f_{j}^{\prime}(\beta_{j})\right]$
$\displaystyle\geq\sqrt{n}\left[O_{p}(1)+\sqrt{n}b_{n}f_{j}^{\prime}(\beta_{j})\right]$

because $p$ is finite and using Slutsky’s Theorem. Since
$\sqrt{n}b_{n}\rightarrow\infty$, the sign of

$\frac{\partial
Q(\beta)}{\partial\beta_{j}}=\sqrt{n}\left[O_{p}(1)+\sqrt{n}\lambda_{j}f_{j}^{\prime}(\beta_{j})\right]$

is determined by the sign of $f_{j}^{\prime}(\beta_{j})$. Thus, since
$f_{j}^{\prime}(0)=0$, we have

$\begin{cases}f_{j}^{\prime}(\beta_{j})<0,&\beta_{j}<0,\\\
f_{j}^{\prime}(\beta_{j})=0,&\beta_{j}=0\\\
f_{j}^{\prime}(\beta_{j})>0,&\beta_{j}>0\end{cases}$

which implies that we have a minimum at $\beta_{j}=0$ and (7.5) is true. ∎

The requirement that $b_{n}$ not go to zero too fast, if it goes to zero at
all, is consistent with the fact that as the $\lambda_{j}$’s increase, the
$\beta_{j}$’s must decrease to obtain an optimal solution.

The proof of the OP for penalized log-likelihood Theorem 2.1 is as follows.

###### Proof.

First, we observe that under Conditions 1–6 and our assumptions on $b_{n}$,
Lemma 2 holds. So,
$\hat{\beta}=(\hat{\beta_{1}}^{\prime},\hat{\beta_{2}}^{\prime})^{\prime}$
satisfies $P(\hat{\beta_{2}}=0)\rightarrow 1$.

Now to prove the asymptotic normality of $\hat{\beta}_{1}$, consider
$\frac{\partial Q(\beta)}{\partial\beta_{j}}$. It can be shown that there
exists a $\hat{\beta_{1}}$ in Theorem 1 that is a $\sqrt{n}$-consistent
minimizer of $Q\left\\{\left(\begin{matrix}\beta_{1}\\\
0\end{matrix}\right)\right\\}$ which is a function of $\beta_{1}$ and
satisfies the likelihood equation

$\frac{\partial
Q(\beta)}{\partial\beta_{j}}\bigg{|}_{\hat{\beta}=\left(\begin{matrix}\hat{\beta_{1}}\\\
0\end{matrix}\right)}=0,$

for $j=1,\ldots,p_{0}$. We see that Taylor expanding at $\beta_{0}$ gives

$\displaystyle 0=\frac{\partial
Q(\beta)}{\partial\beta_{j}}\bigg{|}_{\hat{\beta}=\left(\begin{matrix}\hat{\beta_{1}}\\\
0\end{matrix}\right)}$ $\displaystyle=\frac{\partial
L(\hat{\beta}|x^{n})}{\partial\beta_{j}}+n\frac{\partial}{\partial\beta_{j}}\sum^{p_{0}}_{j=1}\lambda_{j}f_{j}(\hat{\beta_{j}})$
$\displaystyle=\frac{\partial
L({\beta_{0}}|x^{n})}{\partial\beta_{j}}+\sum^{p_{0}}_{j=1}\frac{\partial^{2}{\L(\tilde{\beta}|x^{n})}}{\partial\beta_{\ell}\partial\beta_{j}}(\hat{\beta_{\ell}}-\beta_{\ell
0})+n\lambda_{j}f_{j}^{\prime}(\hat{\beta_{j}})$

where $\tilde{\beta}\in<\beta_{0},\hat{\beta_{1}}>$. Furthermore,

$\displaystyle-\frac{\partial L(\beta_{0}|x^{n})}{\partial\beta_{j}}$
$\displaystyle=\sum^{p_{0}}_{j=1}\frac{\partial^{2}L(\tilde{\beta}|x^{n})}{\partial\beta_{\ell}\partial\beta_{j}}(\hat{\beta_{\ell}}-\beta_{\ell
0})+n\lambda_{j}f_{j}^{\prime}(\hat{\beta_{j}})$ $\displaystyle-\frac{\partial
L(\beta_{0}|x^{n})}{\partial\beta_{j}}$
$\displaystyle=n\sum^{p_{0}}_{j=1}\frac{1}{n}\frac{\partial^{2}L(\tilde{\beta}|x^{n})}{\partial\beta_{\ell}\partial\beta_{j}}(\hat{\beta_{\ell}}-\beta_{\ell
0})+n\lambda_{j}f_{j}^{\prime}(\hat{\beta_{j}})$ $\displaystyle-\frac{\partial
L(\beta_{0}|x^{n})}{\partial\beta_{j}}$
$\displaystyle=-n\left(\sum^{p_{0}}_{j=1}I_{\ell
j}(\beta_{10}|x^{n})\right)(\hat{\beta_{\ell}}-\beta_{\ell
0})+n\lambda_{j}f_{j}^{\prime}(\hat{\beta_{j}})$ $\displaystyle-\frac{\partial
L(\beta_{0}|x^{n})}{\partial\beta_{j}}$
$\displaystyle=-\sqrt{n}\left[\sqrt{n}I_{1}(\beta_{10}|x^{n})(\hat{\beta_{1}}-\beta_{10})+\sqrt{n}\lambda_{j}f_{j}^{\prime}(\hat{\beta_{j}})\right]$
$\displaystyle\frac{1}{\sqrt{n}}\frac{\partial
L(\beta_{0}|x^{n})}{\partial\beta_{j}}$
$\displaystyle=\sqrt{n}I_{1}(\beta_{10}|x^{n})(\hat{\beta_{1}}-\beta_{10})-\sqrt{n}\lambda_{j}f_{j}^{\prime}(\hat{\beta_{j}})$
(7.7)

where the third to last line is due to Condition 2. Here $\hat{I}_{1}$ is like
any $\hat{I}_{n}$ but restricted the upper block of $\hat{I}$.

To finish the proof, we examine the terms in (7.1). By supposition,
$\sqrt{n}a_{n}\rightarrow 0$ and $f_{j}^{\prime}(\hat{\beta_{j}})$ is
uniformly bounded by Condition 4, so
$\sqrt{n}\lambda_{j}f_{j}^{\prime}(\hat{\beta_{j}})\rightarrow 0$ because
$\lambda_{j}<a_{n}$ for $p_{0}<j\leq p$. Then, by the Central Limit Theorem,
Condition 3 and the fact that $E\left(\frac{\partial
L(\beta_{0}|x^{n})}{\partial\beta_{j}}\right)=0$, we have

$\sqrt{n}\left(\frac{1}{n}\frac{\partial
L(\beta_{0}|x^{n})}{\partial\beta_{j}}-E\left(\frac{\partial
L(\beta_{0}|x^{n})}{\partial\beta_{j}}\right)\right)\overset{D}{\longrightarrow}\mathcal{N}(0,I_{1}(\beta_{10})).$

Thus,

$\frac{1}{\sqrt{n}}\frac{\partial
L(\beta_{0}|x^{n})}{\partial\beta_{j}}\overset{D}{\longrightarrow}\mathcal{N}(0,I_{1}(\beta_{10})).$

By equality in (7.1), it follows that

$\sqrt{n}I_{1}(\beta_{10}|x^{n})(\hat{\beta_{1}}-\beta_{10})\overset{D}{\longrightarrow}\mathcal{N}(0,I_{1}(\beta_{10})).$

∎

###  7.2 Penalized Empirical Risk Case, Sec. 2.2

######  Theorem 7.4.

Suppose Conditions 7-10 and Lemma 1 are satisfied. If
$\sqrt{n}a_{n}\rightarrow 0$ where $a_{n}$ satisfies s
$a_{n}=\frac{1}{h(n)\sqrt{n}}$ where $\frac{1}{h(n)\sqrt{n}}\rightarrow 0$ and
$h(n)\rightarrow\infty$ , then there exists a local minimizer $\hat{\beta}$ of
$Q(\beta)$ such that
$||\hat{\beta}-\beta_{0}||=O_{p}(n^{-\frac{1}{2}}+a_{n})$.

###### Proof.

Step 1. We want to show for any $\epsilon>0$ there exists a large constant $C$
such that

$P\left\\{\inf_{||u||=C}Q(\beta_{0}+\alpha_{n}u)>Q(\beta_{0})\right\\}\geq
1-\epsilon.$ (7.8)

Denote

$\displaystyle D_{n}(u)$ $\displaystyle=Q(\beta_{0}+\alpha_{n}u)-Q(\beta_{0})$
$\displaystyle=R(\beta_{0}+\alpha_{n}u|x^{n})+n\sum_{j=1}^{p}\lambda_{j}f(\beta_{j0}+\alpha_{n}u_{j})-R(\beta_{0}|x^{n})-n\sum_{j=1}^{p}\lambda_{j}(f(\beta_{j0}))$
$\displaystyle=R(\beta_{0}+\alpha_{n}u|x^{n})-R(\beta_{0}|x^{n})+n\sum_{j=1}^{p}\lambda_{j}[f_{j}(\beta_{j0}+\alpha_{n}u_{j})-f_{j}(\beta_{j0})]$
$\displaystyle=R(\beta_{0}+\alpha_{n}u|x^{n})-R(\beta_{0}|x^{n})+n\sum_{j=1}^{p}\lambda_{j}f_{j}^{\prime}(\tilde{\beta_{j}})\alpha_{n}u_{j}$
$\displaystyle=R(\beta_{0}+\alpha_{n}u|x^{n})-R(\beta_{0}|x^{n})+n\sum_{j=1}^{p_{0}}\lambda_{j}f_{j}^{\prime}(\tilde{\beta_{j}})\alpha_{n}u_{j}+n\sum_{j=p_{0}+1}^{p}\lambda_{j}f_{j}^{\prime}(\tilde{\beta_{j}})\alpha_{n}u_{j}$

where $f_{j}^{\prime}(\tilde{\beta_{j}})\alpha_{n}u_{j}$ is the Taylor
expansion of the penalty term by Condition 10 and where $\tilde{\beta_{j}}$ is
on the line joining $\beta_{j}$ to $\beta_{j}+\alpha_{n}u$. Thus we have that
the last two terms on the right are the Taylor expansion of the penalty term
by Condition 6. Thus,

$\displaystyle D_{n}(u)$
$\displaystyle=R(\beta_{0}+\alpha_{n}u|x^{n})-R(\beta_{0}|x^{n})+n\alpha_{n}\left(\sum_{j=1}^{p_{0}}\lambda_{j}f_{j}^{\prime}(\tilde{\beta_{j}})u_{j}+n\sum_{j=p_{0}+1}^{p}\lambda_{j}f_{j}^{\prime}(\tilde{\beta_{j}})u_{j}\right)$
$\displaystyle=R(\beta_{0}+\alpha_{n}u|x^{n})-R(\beta_{0}|x^{n})+n(n^{-\frac{1}{2}}+a_{n})\left(\sum_{j=1}^{p_{0}}\lambda_{j}f_{j}^{\prime}(\tilde{\beta_{j}})u_{j}+n\sum_{j=p_{0}+1}^{p}\lambda_{j}f_{j}^{\prime}(\tilde{\beta_{j}})u_{j}\right)$
$\displaystyle=R(\beta_{0}+\alpha_{n}u|x^{n})-R(\beta_{0}|x^{n})+(\sqrt{n}+na_{n})\left(\sum_{j=1}^{p_{0}}\lambda_{j}f_{j}^{\prime}(\tilde{\beta_{j}})u_{j}+n\sum_{j=p_{0}+1}^{p}\lambda_{j}f_{j}^{\prime}(\tilde{\beta_{j}})u_{j}\right)$
$\displaystyle\geq
R(\beta_{0}+\alpha_{n}u|x^{n})-R(\beta_{0}|x^{n})+(\sqrt{n}+na_{n})a_{n}\left(\sum_{j=1}^{p_{0}}f_{j}^{\prime}(\tilde{\beta_{j}})u_{j}+n\sum_{j=p_{0}+1}^{p}f_{j}^{\prime}(\tilde{\beta_{j}})u_{j}\right)$
$\displaystyle=R(\beta_{0}+\alpha_{n}u|x^{n})-R(\beta_{0}|x^{n})+(\sqrt{n}a_{n}+na_{n}^{2})\left(\sum_{j=1}^{p_{0}}f_{j}^{\prime}(\tilde{\beta_{j}})u_{j}+n\sum_{j=p_{0}+1}^{p}f_{j}^{\prime}(\tilde{\beta_{j}})u_{j}\right)$
$\displaystyle=R(\beta_{0}+\alpha_{n}u|x^{n})-R(\beta_{0}|x^{n})+\sqrt{n}a_{n}(1+\sqrt{n}a_{n})\left(\sum_{j=1}^{p_{0}}f_{j}^{\prime}(\tilde{\beta_{j}})u_{j}+n\sum_{j=p_{0}+1}^{p}f_{j}^{\prime}(\tilde{\beta_{j}})u_{j}\right)$
(7.9)

where $p_{0}$ is the number of components in $\beta_{10}$.

Now by Taylor expansion of the empirical risk at $\beta_{0}$,

$\displaystyle R(\beta_{0}+\alpha_{n}u|x^{n})$
$\displaystyle-R(\beta_{0}|x^{n})=\alpha_{n}uR^{\prime}(\beta_{0}|x^{n})+\frac{n}{2}(\alpha_{n}u)^{\prime}\left(I^{*}(\beta_{0}|x^{n})\right)(\alpha_{n}u)$
$\displaystyle+\frac{n}{2}(\alpha_{n}u)^{\prime}\left(R^{\prime\prime}(\tilde{\beta}|x^{n})\pm\frac{1}{n}\sum^{n}_{i=1}I^{*}(\beta|x^{n})\pm
I^{*}(\beta)\right)(\alpha_{n}u)$
$\displaystyle=\alpha_{n}R^{\prime}(\beta_{0}|x^{n})u+\frac{n}{2}u^{\prime}I^{*}(\beta_{0})u\alpha^{2}_{n}\\{1+o_{p}(1)\\}$
(7.10)

because by Condition 7 we have

$E\left[\sup_{\beta\in
B(\beta_{0},\eta)}\left|R^{\prime\prime}(\tilde{\beta}|x^{n})-I^{*}(\beta_{0}|x^{n})\right|\right]\rightarrow
0$

as $\eta\rightarrow 0$ and where
$I^{*}(\beta_{0}|x^{n})=E_{L}(R^{\prime\prime}_{jk}(\beta_{0}|x^{n}))_{jk}$.

Using $(\ref{dnu2})$ and $(\ref{lnu2})$ ,

$\displaystyle D_{n}(u)$ $\displaystyle\geq$
$\displaystyle\alpha_{n}R^{\prime}(\beta_{0}|x^{n})u+\frac{n}{2}u^{\prime}I^{*}(\beta_{0}|x^{n})u\alpha^{2}_{n}\\{1+o_{p}(1)\\}$
(7.11) $\displaystyle+$
$\displaystyle\sqrt{n}a_{n}(1+\sqrt{n}a_{n})\sum_{j=1}^{p_{0}}f_{j}^{\prime}(\tilde{\beta_{j}})u_{j}$
$\displaystyle+$
$\displaystyle\sqrt{n}a_{n}(1+\sqrt{n}a_{n})\sum_{j=p_{0}+1}^{p}f_{j}^{\prime}(\tilde{\beta_{j}})u_{j}.$

Now we argue that the second term on the RHS of (7.11) dominates the others.
To see this, consider the first term on the RHS of (7.11). We multiply by
$\frac{\sqrt{n}}{\sqrt{n}}$ and we have
$\frac{\sqrt{n}}{\sqrt{n}}\alpha_{n}R^{\prime}(\beta_{0}|x^{n})u$. Let
$\bar{R}=\sqrt{n}\frac{\sum^{n}_{i=1}d^{\prime}(y_{i}-x_{i}^{\prime}\beta)}{n}$.
Thus we rewrite
$\frac{\sqrt{n}}{\sqrt{n}}\alpha_{n}R^{\prime}(\beta_{0}|x^{n})u$ as
$\frac{\alpha_{n}}{\sqrt{n}}\bar{R}u$. Now by the Central Limit Theorem,

$\sqrt{n}\left[\frac{\sum^{n}_{i=1}d^{\prime}(y_{i}-x_{i}^{\prime}\beta)}{n}-E[d^{\prime}(y_{i}-x_{i}^{\prime}\beta)]\right]$

converges in distribution to a normal distribution since
$E[d^{\prime}(y_{i}-x_{i}^{\prime}\beta)]=0$ due to Lemma 1 and because
$\frac{\alpha_{n}}{\sqrt{n}}\rightarrow 1$. Thus,
$C^{2}=||u||^{2}=u^{\prime}u>\bar{R}u$ as long as $0<\frac{1}{B}<\bar{R}<B$
for $C>B$ with high probability. Then with high probability $C^{2}\geq
B\underline{1}(u)\geq B\underline{1}u$ for $||u||=C$. Thus we have that by
choosing sufficiently large $C$, the second term of (7.4) is larger in
absolute value, with high probability than the first term uniformly in
$||u||=C$.

Also, by hypothesis $\sqrt{n}a_{n}\rightarrow 1$, so the whole third term in
(7.11) goes to zero and the second term of (7.11) is also larger than the
third term. The fourth term of (7.11) also goes to zero because for $p_{0}\leq
j\leq p$, we know that $\beta_{j}=0$, so $\tilde{\beta_{j}}\rightarrow 0$ as
$n\rightarrow 0$.Therefore
$\sum_{j=p_{0}+1}^{p}f_{j}^{\prime}(\tilde{\beta_{j}})u_{j}\rightarrow 0$ by
Condition 9. It follows that the third term of (7.11) is larger than the
fourth term, which implies that the second term of (7.11) is larger than all
of the other terms.

Therefore, by choosing large enough $C$,

$P\left\\{\inf_{||u||=C}D_{n}(u)\geq\frac{n}{2}u^{\prime}I^{*}(\beta_{0}|x^{n})u\alpha^{2}_{n}\\{1+o_{p}(1)\\}\right\\}\geq
1-\epsilon.$

Thus, (7.8) is true for sufficiently large $C$.

Step 2. We now show the conclusion of Step 1 holds for any $C^{*}>C$. Let
$C^{\star}>C$. Then,

$\displaystyle\hat{\beta}$
$\displaystyle\in\\{\beta_{0}+\alpha_{n}u:||u||=C^{*}\\}$
$\displaystyle\iff\hat{\beta}\in\cal{B}\left(\beta_{0},\alpha_{n}C^{*}\right)$
$\displaystyle\iff\hat{\beta}-\beta_{0}\in\cal{B}\left(0,\alpha_{n}C^{*}\right)$
$\displaystyle\iff\alpha_{n}^{-1}(\hat{\beta}-\beta_{0})\in\cal{B}(0,C^{*})$
$\displaystyle\iff||\alpha_{n}^{-1}(\hat{\beta}-\beta_{0})||<C^{*}$
$\displaystyle\iff||\hat{\beta}-\beta_{0}||^{2}<\alpha_{n}^{2}{C^{*}}^{2}$

It follows that
$\sum^{p}_{j=1}(\hat{\beta}_{j}-\beta_{0j})^{2}\leq\alpha_{n}^{2}{C^{*}}^{2}$
and $\forall j$ ,
$(\hat{\beta}_{j}-\beta_{0j})^{2}\leq\alpha_{n}^{2}{C^{*}}^{2}$. So,

$-\alpha_{n}C^{*}\leq\hat{\beta}_{j}-\beta_{0j}\leq\alpha_{n}C^{*}$

and we can absorb $C^{*}$ into $\alpha_{n}$ because $C^{*}$ is a constant.
Thus we have $\hat{\beta}-\beta_{0}=O_{p}(\alpha_{n})$. Note that

$\alpha_{n}=\frac{1}{\sqrt{n}}+a_{n}=\frac{1+\sqrt{n}a_{n}}{\sqrt{n}}$

and as $n\rightarrow\infty$ we have $\sqrt{n}a_{n}\rightarrow 0$, so
$\hat{\beta}-\beta_{0}=O_{p}\left(\frac{1}{\sqrt{n}}\right)$. ∎

Again, we discuss the way you can get this result for penalty functions that
have finitely many points that are non-differentiable, particularly at the
origin in the following remark.

######  Remark.

The argument for Theorem 1 is true for differentiable $f_{j}(\cdot)$, but not
non-differentiable $f_{j}(\cdot)$. In the case of non-differentiable
$f_{j}(\cdot)$, consider a smooth function that approximates $f_{j}(\cdot)$
well. Say, $f_{j}^{\star}(\cdot)$, which differs from $f_{j}(\cdot)$ by a
margin of $\epsilon$ where $\epsilon_{N}\rightarrow 0$ as
$n\rightarrow\infty$. As long as
$\lim_{\beta\rightarrow\beta^{*-}}f_{j}^{\star}(\cdot)=\lim_{\beta\rightarrow\beta^{*+}}f_{j}^{\star}(\cdot)$,
then the above argument holds.

######  Lemma 3.

Assume conditions 7–10 hold. If $\sqrt{n}b_{n}\rightarrow\infty$ as
$n\rightarrow\infty$ and and $b_{n}$ satisfies $b_{n}=\frac{g(n)}{\sqrt{n}}$
where $\frac{g(n)}{\sqrt{n}}\rightarrow 0$ and $g(n)\rightarrow\infty$, then
with probability tending to 1, for any given $\beta_{1}$ satisfying
$||\beta_{1}-\beta_{10}||=O_{p}(n^{-\frac{1}{2}})$ and any constant $C$,

$Q\left\\{\left(\begin{matrix}\beta_{1}\\\
0\end{matrix}\right)\right\\}=\max_{||\beta_{2}||\leq
Cn^{-\frac{1}{2}}}Q\left\\{\left(\begin{matrix}\beta_{1}\\\
\beta_{2}\end{matrix}\right)\right\\}.$ (7.12)

###### Proof.

Consider the objective function
$Q(\beta)=R(\beta|x^{n})+n\sum^{p}_{j=1}\lambda_{j}f_{j}(\beta_{j})$. Note
that

$\frac{\partial Q(\beta)}{\partial\beta_{j}}=\frac{\partial
R(\beta|x^{n})}{\partial\beta_{j}}+n\lambda_{j}f_{j}^{\prime}(\beta_{j}).$
(7.13)

Then by Taylor expanding at $\beta_{0}=\left(\begin{matrix}\beta_{10}\\\
0\end{matrix}\right)$ we have

$\displaystyle\frac{\partial Q(\beta)}{\partial\beta_{j}}=\frac{\partial
R(\beta_{0}|x^{n})}{\partial\beta_{j}}$
$\displaystyle+\sum^{p}_{\ell=1}\frac{\partial^{2}R(\beta_{0}|x^{n})}{\partial\beta_{j}\partial\beta_{\ell}}(\beta_{\ell}-\beta_{\ell
0})$
$\displaystyle+\sum^{p}_{\ell=1}\sum^{p}_{k=1}\frac{\partial^{3}R(\beta^{*}|x^{n})}{\partial\beta_{j}\partial\beta_{\ell}\partial\beta_{k}}(\beta_{\ell}-\beta_{\ell
0})(\beta_{k}-\beta_{k0})+n\lambda_{j}f_{j}^{\prime}(\beta_{j})$

where $\beta^{*}$ lies between $\hat{\beta}$ and $\beta_{0}$. Note that

$\frac{\partial
R(\beta_{0}|x^{n})}{\partial\beta_{j}}=O_{p}(n^{-\frac{1}{2}})$

and due to [22] and by the Law of Large Numbers for non-identically
distributed random variables we have

$\frac{\partial^{2}R(\beta_{0}|x^{n})}{\partial\beta_{j}\partial\beta_{\ell}}=E\left[\frac{\partial^{2}R(\beta_{0}|x^{n})}{\partial\beta_{j}\partial\beta_{\ell}}\right]+o_{p}(1).$

Also, due to Theorem 1,
$\hat{\beta}-\beta=O_{p}\left(\frac{1}{\sqrt{n}}\right)$, so we have

$\displaystyle\frac{\partial Q(\beta)}{\partial\beta_{j}}$
$\displaystyle=\frac{\partial
R(\beta_{0}|x^{n})}{\partial\beta_{j}}+\sum^{p}_{\ell=1}\frac{\partial^{2}R(\beta_{0}|x^{n})}{\partial\beta_{j}\partial\beta_{\ell}}(\hat{\beta_{\ell}}-\beta_{\ell
0})$
$\displaystyle+\sum^{p}_{\ell=1}\sum^{p}_{k=1}\frac{\partial^{3}R(\beta^{*}|x^{n})}{\partial\beta_{j}\partial\beta_{\ell}\partial\beta_{k}}(\hat{\beta_{\ell}}-\beta_{\ell
0})(\hat{\beta_{k}}-\beta_{k0})+n\lambda_{j}f_{j}^{\prime}(\beta_{j})$
$\displaystyle=O_{p}\left(\frac{1}{\sqrt{n}}\right)+\sum^{p}_{\ell=1}\left(E\left[\frac{\partial^{2}R(\beta_{0}|x^{n})}{\partial\beta_{j}\partial\beta_{\ell}}\right]+o_{p}(1)\right)O_{p}\left(\frac{1}{\sqrt{n}}\right)$
$\displaystyle+\sum^{p}_{\ell=1}\sum^{p}_{k=1}E\left[\sup_{\beta\in
N}\left|\frac{\partial^{3}R(\beta^{*}|x^{n})}{\partial\beta_{j}\partial\beta_{\ell}\partial\beta_{k}}\right|+o_{p}(1)\right]O_{p}\left(\frac{1}{\sqrt{n}}\right)O_{p}\left(\frac{1}{\sqrt{n}}\right)+n\lambda_{j}f_{j}^{\prime}(\beta_{j})$
$\displaystyle=\left[O_{p}\left(\frac{1}{\sqrt{n}}\right)+O_{p}\left(\frac{1}{\sqrt{n}}\right)+O_{p}\left(\frac{1}{\sqrt{n}}\right)\right]+n\lambda_{j}f_{j}^{\prime}(\beta_{j})$
$\displaystyle=O_{p}\left(\frac{1}{\sqrt{n}}\right)+n\lambda_{j}f_{j}^{\prime}(\beta_{j})$
$\displaystyle=\sqrt{n}\left[O_{p}\left(\frac{1}{n}\right)+\sqrt{n}\lambda_{j}f_{j}^{\prime}(\beta_{j})\right]$
$\displaystyle\geq\sqrt{n}\left[O_{p}\left(\frac{1}{n}\right)+\sqrt{n}b_{n}f_{j}^{\prime}(\beta_{j})\right]$

because $p$ is finite and using Slutsky’s Theorem. Since
$\sqrt{n}b_{n}\rightarrow\infty$, the sign of

$\frac{\partial
Q(\beta)}{\partial\beta_{j}}=\sqrt{n}\left[O_{p}\left(\frac{1}{n}\right)+\sqrt{n}\lambda_{j}f_{j}^{\prime}(\beta_{j})\right]$

is determined by the sign of $f_{j}^{\prime}(\beta_{j})$. That is, since
$f_{j}^{\prime}(0)=0$, we have

$\begin{cases}f_{j}^{\prime}(\beta_{j})<0,&\beta_{j}<0,\\\
f_{j}^{\prime}(\beta_{j})=0,&\beta_{j}=0\\\
f_{j}^{\prime}(\beta_{j})>0,&\beta_{j}>0\end{cases}$

which implies that we have a minimum at $\beta_{j}=0$ and (7.12) is true. ∎

The proof of the OP for penalized log-likelihood Theorem 2.2 is as follows.

###### Proof.

First, we observe that the estimator
$\hat{\beta}=(\hat{\beta_{1}}^{\prime},\hat{\beta_{2}}^{\prime})^{\prime}$
satisfies $P(\hat{\beta_{2}}=0)\rightarrow 1$ due to Lemma 3.

Now to prove the asymptotic normality piece, consider $\frac{\partial
Q(\beta)}{\partial\beta_{j}}$. It can be shown that there exists a
$\hat{\beta_{1}}$ in Theorem 3 that is a $\sqrt{n}$ consistent minimizer of
$Q\left\\{\left(\begin{matrix}\beta_{1}\\\ 0\end{matrix}\right)\right\\}$
which is a function of $\beta_{1}$ and satisfies the likelihood equations

$\frac{\partial
Q(\beta)}{\partial\beta_{j}}\bigg{|}_{\hat{\beta}=\left(\begin{matrix}\hat{\beta_{1}}\\\
0\end{matrix}\right)}=0$

for $j=1,\ldots,p_{0}$.

We see that Taylor expanding at $\beta_{0}$ gives

$\displaystyle 0=\frac{\partial
Q(\beta)}{\partial\beta_{j}}\bigg{|}_{\hat{\beta}=\left(\begin{matrix}\hat{\beta_{1}}\\\
0\end{matrix}\right)}$ $\displaystyle=\frac{\partial
R(\hat{\beta}|x^{n})}{\partial\beta_{j}}+\frac{\partial}{\partial\beta_{j}}n\sum^{p_{0}}_{j=1}\lambda_{j}f_{j}(\hat{\beta_{j}})$
$\displaystyle=\frac{\partial
R(\beta_{0}|x^{n})}{\partial\beta_{j}}+\sum^{p_{0}}_{j=1}\frac{\partial^{2}R(\tilde{\beta}|x^{n})}{\partial\beta_{\ell}\partial\beta_{j}}(\hat{\beta_{\ell}}-\beta_{\ell
0})+n\lambda_{j}f_{j}^{\prime}(\hat{\beta_{j}})$

where $\tilde{\beta}\in<\beta_{0},\hat{\beta_{1}}>$. Furthermore,

$\displaystyle-\frac{\partial R(\beta_{0}|x^{n})}{\partial\beta_{j}}$
$\displaystyle=\sum^{p_{0}}_{j=1}\frac{\partial^{2}R(\tilde{\beta}|x^{n})}{\partial\beta_{\ell}\partial\beta_{j}}(\hat{\beta_{\ell}}-\beta_{\ell
0})+n\lambda_{j}f_{j}^{\prime}(\hat{\beta_{j}})$ $\displaystyle-\frac{\partial
R(\beta_{0}|x^{n})}{\partial\beta_{j}}$
$\displaystyle=n\sum^{p_{0}}_{j=1}\frac{1}{n}\frac{\partial^{2}R(\tilde{\beta}|x^{n})}{\partial\beta_{\ell}\partial\beta_{j}}(\hat{\beta_{\ell}}-\beta_{\ell
0})+n\lambda_{j}f_{j}^{\prime}(\hat{\beta_{j}})$ $\displaystyle-\frac{\partial
R(\beta_{0}|x^{n})}{\partial\beta_{j}}$
$\displaystyle=-n\sum^{p_{0}}_{j=1}I^{*}_{\ell
j}(\beta_{10}|x^{n})(\hat{\beta_{\ell}}-\beta_{\ell
0})+n\lambda_{j}f_{j}^{\prime}(\hat{\beta_{j}})$ $\displaystyle-\frac{\partial
R(\beta_{0}|x^{n})}{\partial\beta_{j}}$
$\displaystyle=-\sqrt{n}\left[\sqrt{n}I^{*}_{1}(\beta_{10}|x^{n})(\hat{\beta_{1}}-\beta_{10})+\sqrt{n}\lambda_{j}f_{j}^{\prime}(\hat{\beta_{j}})\right]$
$\displaystyle\frac{1}{\sqrt{n}}\frac{\partial
R(\beta_{0}|x^{n})}{\partial\beta_{j}}$
$\displaystyle=\sqrt{n}I^{*}_{1}(\beta_{10}|x^{n})(\hat{\beta_{1}}-\beta_{10})-\sqrt{n}\lambda_{j}f_{j}^{\prime}(\hat{\beta_{j}})$
(7.14)

where the third to last line is due Condition 7. Here $I^{*}_{1}$ is like any
$I^{*}_{n}$ but restricted the upper block of $I^{*}$.

To finish the proof, we examine the terms in (7.2). By supposition,
$\sqrt{n}a_{n}\rightarrow 0$ and $f_{j}^{\prime}(\hat{\beta})$ is bounded.
Thus, $\sqrt{n}\lambda_{j}f_{j}^{\prime}(\hat{\beta_{j}})\rightarrow 0$ for