# Fitted value shrinkage

Daeyoung Ham & Adam J. Rothman
School of Statistics, University of Minnesota
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We propose a penalized least-squares method to fit the linear regression model
with fitted values that are invariant to invertible linear transformations of
the design matrix. This invariance is important, for example, when
practitioners have categorical predictors and interactions. Our method has the
same computational cost as ridge-penalized least squares, which lacks this
invariance. We derive the expected squared distance between the vector of
population fitted values and its shrinkage estimator as well as the tuning
parameter value that minimizes this expectation. In addition to using cross
validation, we construct two estimators of this optimal tuning parameter value
and study their asymptotic properties. Our numerical experiments and data
examples show that our method performs similarly to ridge-penalized least-
squares.

Keywords: invariance, penalized least squares, high-dimensional data

##  1 Introduction

We will introduce a new shrinkage strategy for fitting linear regression
models, which assume that the measured response for $n$ subjects is a
realization of the random vector

$Y=X\beta+\varepsilon,$ (1)

where $X\in\mathbb{R}^{n\times p}$ is the nonrandom known design matrix with
ones in its first column and with values of the predictors in its remaining
columns; $\beta\in\mathbb{R}^{p}$ is an unknown vector of regression
coefficients; and $\varepsilon$ has iid entries with mean zero and unknown
variance $\sigma^{2}\in(0,\infty)$.

We will consider fitting (1) in both low and high-dimensional settings, where
the second scenario typically has ${\rm rank}(X)<p$. If ${\rm rank}(X)<p$,
then it is well known that $\beta$ is not identifiable in (1), i.e. there
exists a $\tilde{\beta}\neq\beta$ such that $X\beta=X\tilde{\beta}$.
Similarly, if ${\rm rank}(X)<p$, then there are infinitely many solutions to
the least-squares problem: $\operatorname*{arg\
min}_{b\in\mathbb{R}^{p}}\|Y-Xb\|^{2}$. Given this issue (which is unavoidable
in high dimensions), our inferential target is $X\beta$, which is the expected
value of the response for the $n$ subjects.

To describe least squares estimators whether ${\rm rank}(X)<p$ or ${\rm
rank}(X)=p$, we will use the reduced singular value decomposition of $X$. Let
$q={\rm rank}(X)$. Then $X=UDV^{\prime}$, where $U\in\mathbb{R}^{n\times q}$
with $U^{\prime}U=I_{q}$; $V\in\mathbb{R}^{p\times q}$ with
$V^{\prime}V=I_{q}$; and $D\in\mathbb{R}^{q\times q}$ is diagonal with
positive diagonal entries. The Moore–Penrose generalized inverse of $X$ is
$X^{-}=VD^{-1}U^{\prime}$ and a least-squares estimator of $\beta$ is
$\hat{\beta}=X^{-}Y$. The vector of fitted values is
$X\hat{\beta}=XX^{-}Y=P_{X}Y$, where $P_{X}=XX^{-}=UU^{\prime}$. If ${\rm
rank}(X)=p$, then $P_{X}=X(X^{\prime}X)^{-1}X^{\prime}$.

A nice property of this least-squares method is that its fitted values are
invariant to invertible linear transformations of the design matrix. Suppose
that we replace $X$ by $X_{\bullet}=XT$, where $T\in\mathbb{R}^{p\times p}$ is
invertible. Then $X=X_{\bullet}T^{-1}$. So (1) is

$Y=X\beta+\varepsilon=X_{\bullet}T^{-1}\beta+\varepsilon=X_{\bullet}\beta_{\bullet}+\varepsilon,$

where $\beta_{\bullet}=T^{-1}\beta$. We estimate
$X\beta=X_{\bullet}\beta_{\bullet}$ with $P_{X}Y=P_{X_{\bullet}}Y$, so the
fitted values did not change by changing $X$ to $X_{{\bullet}}$.

Fitting (1) by penalized least-squares has been studied by many scholars.
Well-studied penalties include the ridge penalty (Hoerl and Kennard, 1970),
the bridge/lasso penalty (Frank and Friedman, 1993; Tibshirani, 1996), the
adaptive lasso penalty (Zou, 2006), the SCAD penalty (Fan and Li, 2001), and
the MCP penalty (Zhang, 2010). Unfortunately, these methods’ fitted values are
not invariant to invertible linear transformations of $X$. This is
particularly problematic when categorical variables (with three or more
categories) and their interactions are encoded in $X$ because a change in the
coding can change the fit. The Group Lasso (Yuan and Lin, 2006) and its
variants with different standardization (Choi et al., 2012; Simon and
Tibshirani, 2012) could be used when there are categorical predictors, but
these methods still lack invariance and our simulation and data examples
illustrate their instability. This lack of invariance is also present in
principal components regression (Hotelling, 1933), and partial least squares
(Wold, 1966).

Invariance to invertible linear transformations is important in any
application that includes categorical predictors and their interactions: the
model fit should not change by changing the coding scheme used to represent
these variables in the design matrix. We illustrate that a change in the
reference level coding of categorical predictors (see Section 6.2) and that a
rotational transformation of the design matrix that removes zeros from the
regression coefficient vector (see Section 6.3 and 6.5) leads the regression
shrinkage methods that lack this invariance to perform worse. In contrast, our
method performed the same before and after these transformations.

##  2 A new shrinkage method for linear regression with invariance

###  2.1 Method description

To preserve the invariance to invertible linear transformations of the design
matrix discussed in the previous section, we will use penalties that can be
expressed as a function of the $n$-dimensional vector $Xb$, where $b$ is the
optimization variable corresponding to $\beta$. To construct our estimator, we
start with the following penalized least squares optimization:

$\displaystyle\operatorname*{arg\
min}_{b\in\mathbb{R}^{p}}\left\\{\|Y-Xb\|^{2}+\lambda\|Xb-\bar{Y}1_{n}\|^{2}\right\\},$
(2)

where $\bar{Y}=1_{n}^{\prime}Y/n$;
$1_{n}^{\prime}=(1,\dots,1)\in\mathbb{R}^{n}$; and $\lambda\in[0,\infty)$ is a
tuning parameter. As $\lambda$ increases, the fitted values are shrunk towards
the intercept-only model’s fitted values $\bar{Y}{1_{n}}$. Let
$\gamma=1/(1+\lambda)$. We can express the optimization in (2) as

$\hat{\beta}^{(\gamma)}\in\operatorname*{arg\
min}_{{b}\in\mathbb{R}^{p}}\left\\{\gamma\|{Y}-Xb\|^{2}+(1-\gamma)\|Xb-\bar{Y}{1_{n}}\|^{2}\right\\},$
(3)

where $\gamma\in[0,1]$. The $\in$ is used because there are infinitely many
global minimizers for the optimization in (3) when ${\rm rank}(X)<p$. To
preserve invariance and obtain the minimum 2-norm solution, we define our
estimator as

$\hat{\beta}^{(\gamma)}={X}^{-}\left\\{\gamma{Y}+(1-\gamma)\bar{Y}{1_{n}}\right\\},$

which is a global minimizer of (3) that uses the Moore-Penrose generalized
inverse of $X$. Let
${P_{1}}={1_{n}}({1_{n}}^{\prime}{1_{n}})^{-1}{1_{n}}^{\prime}$. The estimator
of ${X}{\beta}$ is

$X\hat{\beta}^{(\gamma)}=\gamma P_{X}{Y}+(1-\gamma){P_{1}}{Y},$ (4)

which is simply a convex combination of the least-squares fitted values
${P_{X}}{Y}$ and the intercept-only model’s fitted values
${P_{1}}{Y}=\bar{Y}{1_{n}}$. Since $P_{X}=P_{X_{\bullet}}$, where
$X_{\bullet}=XT$ with $T\in\mathbb{R}^{p\times p}$ invertible,
$X\hat{\beta}^{(\gamma)}$ is invariant to invertible linear transformations of
$X$.

Due to its computational simplicity, $\hat{\beta}^{(\gamma)}$ is a natural
competitor to ridge-penalized least squares, which lacks this invariance
property. Both methods can be computed efficiently when $p$ is much larger
than $n$ by using the reduced singular value decomposition of $X$ (Hastie and
Tibshirani, 2004). Specifically, they both cost $O(n\ {\rm rank}^{2}(X))$
floating-point operations. If $\gamma=0$ for our method and
$\lambda\rightarrow\infty$ for ridge-penalized least squares (without
intercept penalization), then both procedures fit the intercept-only model.

We will derive an optimal value of $\gamma$ that minimizes
$\mathbb{E}\|X\hat{\beta}^{(\gamma)}-X\beta\|^{2}$ and propose two estimators
of it: one for low dimensions and one for high dimensions. We also explore
using cross validation to select $\gamma$ when the response and predictor
measurement pairs are drawn from a joint distribution. Conveniently, our
results generalize to shrinkage towards a submodel’s fitted values
$P_{X_{0}}Y$, where $X_{0}$ is a matrix with a proper subset of the columns of
$X$.

###  2.2 Related work

Copas (1997) proposed to predict a future value of the response for the $i$th
subject with a convex combination of its fitted value (from ordinary least
squares) and $\bar{Y}$. Although our methods are related, Copas (1997) used a
future-response-value prediction paradigm and did not establish a theoretical
analysis of his approach.

Azriel and Schwartzman (2020) study the estimation of $a^{\prime}\beta$ when
$p>n$ and $\beta$ is not sparse. They establish a condition for which
$a^{\prime}\beta$ is identifiable and show that using least-squares (with a
pseudoinverse) has an optimal property when the errors are Gaussian. We prove
that our shrinkage method outperforms least squares in same-$X$ prediction and
we illustrate it in our numerical examples. Zhao et al. (2023) also study the
estimation of $a^{\prime}\beta$, but they use procedures that are not
invariant to invertible linear transformations of $X$.

##  3 Theoretical properties of the method

Given the lack of identifiability of $\beta$ in high dimensions, we
investigate the estimation of the $n$-dimensional vector $X\beta$ with
$X\hat{\beta}^{(\gamma)}$. This is an example of same-X prediction (Rosset and
Tibshirani, 2020). It is related to predicting near $X$ when $p>n$ (Cook,
Forzani, and Rothman, 2013, Proposition 3.4). In contrast to a random-design
analysis, our fixed-design analysis allows a practitioner to control a subset
of the columns of $X$, which would be essential for designed experiments.

Suppose that the linear regression model specified in (1) is true (this model
did not specify an error distribution, just that they are iid mean 0 and
variance $\sigma^{2}\in(0,\infty)$). We define $\mu=X\beta$ and we assume that
${\rm rank}(X)\geq 2$ throughout the paper, so $X$ cannot correspond to an
intercept-only model. Then we have the following result:

######  Proposition 1.

For all $(n,p)\in\\{1,2,\ldots\\}\times\\{1,2,\ldots\\}$,

$\mathbb{E}\|X\hat{\beta}^{(\gamma)}-X{\beta}\|^{2}=\sigma^{2}(\gamma^{2}r+1-\gamma^{2})+(1-\gamma)^{2}\|\mu-
P_{1}\mu\|^{2}.$

The proof of Proposition 1 is in Section .1 of the Supplementary material.
When $\gamma=1$, which is least squares,
$\mathbb{E}\|X\hat{\beta}^{(1)}-X\beta\|^{2}=\sigma^{2}{\rm rank}(X)$.

The right side of the equality in Proposition 1 is minimized when
$\gamma=\gamma_{\rm opt}$, where

$\displaystyle\gamma_{\rm opt}=\frac{\|\mu-P_{1}\mu\|^{2}}{\sigma^{2}({\rm
rank}(X)-1)+\|\mu-P_{1}\mu\|^{2}}.$ (5)

So the best our procedure could do is when $\|\mu-P_{1}\mu\|^{2}=0$ (that is,
the intercept-only model is correct), in which case $\gamma_{\rm opt}=0$ and
$\mathbb{E}\|X\hat{\beta}^{(0)}-X\beta\|^{2}=\sigma^{2}$. The expression for
$\gamma_{\rm opt}$ enables us to construct a sample-based one-step estimator
of it, which is more computationally efficiency than cross validation.

##  4 Selection of $\gamma$

###  4.1 Low-dimensional case

Let $\hat{\sigma}^{2}=\|Y-P_{X}Y\|^{2}/(n-{\rm rank}(X))$, which is an
unbiased estimator of $\sigma^{2}$. To construct an estimator of $\gamma_{\rm
opt}$, we use the ratio of $\|P_{X}Y-P_{1}Y\|^{2}-\hat{\sigma}^{2}({\rm
rank}(X)-1)$, which is an unbiased estimator of $\gamma_{\rm opt}$’s
numerator, to $\|P_{X}Y-P_{1}Y\|^{2}$, which is an unbiased estimator of
$\gamma_{\rm opt}$’s denominator. This ratio estimator can be expressed as

$\displaystyle\frac{\|P_{X}Y-P_{1}Y\|^{2}-\hat{\sigma}^{2}({\rm
rank}(X)-1)}{\|P_{X}Y-P_{1}Y\|^{2}}$
$\displaystyle=1-\frac{\hat{\sigma}^{2}({\rm
rank}(X)-1)}{\|P_{X}Y-P_{1}Y\|^{2}}$ $\displaystyle=1-1/F,$

where $F$ is the F statistic that compares the intercept-only model to the
full model:

$\displaystyle F$
$\displaystyle=\frac{(\|Y-P_{1}Y\|^{2}-\|Y-P_{X}Y\|^{2})/({\rm
rank}(X)-1)}{\|Y-P_{X}Y\|^{2}/(n-{\rm rank}(X))}$
$\displaystyle=\frac{\|P_{X}Y-P_{1}Y\|^{2}}{\hat{\sigma}^{2}({\rm
rank}(X)-1)}.$ (6)

Since $F$ could be realized less than one (which corresponds to a fail-to-
reject the intercept-only model situation), we define our estimator of
$\gamma_{\rm opt}$ to be

$\displaystyle\hat{\gamma}=(1-1/F)\cdot 1(F>1)$ (7)

If the regression errors in (1) are Normal, $n>{\rm rank}(X)$, and ${\rm
rank}(X)>1$, then $F$ has a non-central F-distribution with degrees of freedom
parameters ${\rm rank}(X)-1$ and $n-{\rm rank}(X)$; and noncentrality
parameter $\|\mu-P_{1}\mu\|^{2}/\sigma^{2}$. Larger realizations of $F$
correspond to worse intercept-only model fits compared to the full model,
which makes $\hat{\gamma}$ closer to 1.

We also explore two additional estimators of $\gamma_{\rm opt}$:

$\displaystyle\hat{\gamma}_{90}=(1-1/F)\cdot 1(F\geq f_{0.9}),$ (8)
$\displaystyle\hat{\gamma}_{95}=(1-1/F)\cdot 1(F\geq f_{0.95}),$ (9)

where $f_{90}$ and $f_{95}$ are the $0.9$ and $0.95$ quantiles of the central
F-distribution with degrees of freedom ${\rm rank}(X)-1$ and $n-{\rm
rank}(X)$. These estimators may perform better when $\gamma_{\rm opt}$ is near
zero because they have a greater probability of estimating $\gamma_{\rm opt}$
as zero than $\hat{\gamma}$ has.

Interestingly, Copas (1997) proposed to predict a future response value for
the $i$th subject with $(1-\rho)\hat{\beta}^{\prime}x_{i}+\rho\bar{Y}$, where
$\rho\in[0,1]$ is estimated and $\hat{\beta}$ is the ordinary least-squares
estimator. They derived $1/F$ as an estimator of $\rho$ from the normal
equations for the regression of $Y_{{\rm new},i}$ on
$\hat{\beta}^{\prime}x_{i}$, $(i=1,\ldots,n)$, where $Y_{{\rm new},i}$ is an
independent copy of $Y_{i}$. They also discussed using truncation to ensure
their estimator of $\rho$ is in $[0,1]$.

###  4.2 Consistency and the convergence rate of $\hat{\gamma}$

We analyze the asymptotic performance of $\hat{\gamma}$ when the data are
generated from (1) and $n$ and $p$ grow together. Define $r={\rm rank}(X)$ and
$\delta^{2}=\|\mu-P_{1}\mu\|^{2}$. The optimal tuning parameter value is a
function of $r$ and $\delta^{2}$, so its value in the limit will depend on
these sequences.

######  Proposition 2.

Assume that the data-generating model in (1) is correct, that the errors have
a finite fourth moment, and that $r\geq 2$. If $p/n\rightarrow\tau\in[0,1)$
and either $r\rightarrow\infty$ or $\delta^{2}\rightarrow\infty$, then
$\hat{\gamma}-\gamma_{\rm opt}\rightarrow_{P}0$ as $n\rightarrow\infty$.

The proof of Proposition 2 is in Section .2 of the Supplementary material. We
see that consistency is possible whether the design matrix rank $r$ grows. If
$r$ is bounded, then consistency requires $\delta^{2}=\|\mu-
P_{1}\mu\|^{2}\rightarrow\infty$, which is reasonable even when the intercept-
only model is a good approximation because $n$ is growing. One can also show
consistency of $\hat{\gamma}_{90}$ and $\hat{\gamma}_{95}$ with
$\delta^{2}=o(r)$ added to the assumptions for Proposition 2.

Next, we establish a bound on the rate of convergence of $\hat{\gamma}$ with
further assumptions on the design matrix $X$ and the error $\varepsilon$.

######  Proposition 3.

Suppose that the assumptions of Proposition 2 hold, that the errors in (1) are
Gaussian, and that $r\geq 6$ is nondecreasing as $n\rightarrow\infty$. Then

$\hat{\gamma}-\gamma_{\rm opt}=\begin{cases}O_{P}(r^{-1/2})&\text{if
$\delta^{2}=O(r)$}\\\
O_{P}((\delta^{2}/r)^{-3/4})+O_{P}(n^{-1/2}(\delta^{2}/r)^{-1/2})&\text{if
$r\rightarrow\infty$ and $r=o(\delta^{2})$}\\\
O_{P}((\delta^{2})^{-3/4})+O_{P}(n^{-1/2}(\delta^{2})^{-1/2})&\text{if
$r=O(1)$.}\\\ \end{cases}$ (10)

The proof of Proposition 3 is in Section .3 of the Supplementary material.
From the definition of $\gamma_{\rm opt}$, we know that $\gamma_{\rm
opt}\rightarrow 0$ when $\delta^{2}=o(r)$, in which case
$\hat{\gamma}-\gamma_{\rm opt}=O_{P}(n^{-1/2})$ provided that $r\asymp n$. On
the other hand, $\gamma_{\rm opt}\rightarrow 1$ when $r=o(\delta^{2})$. For
example, $\hat{\gamma}-\gamma_{\rm opt}=O_{P}(n^{-1/2})$ provided that
$\delta^{2}\asymp n^{2/3}$ and $r$ is bounded. When $\delta^{2}=o(r)$,
$\hat{\gamma}_{90},\hat{\gamma}_{95}$ and $\hat{\gamma}$ all have the same
convergence rate bound.

###  4.3 Tuning parameter selection in high dimensions

Estimating the unknown parameters in $\gamma_{\rm opt}$ is challenging when
$p>n$ and ${\rm rank}(X)=n$. For example, it is impossible to estimate the
regression’s error variance $\sigma^{2}$ without assuming something extra
about $\mu=X\beta$. This is because the data-generating model in (1) reduces
to

$Y=\mu+\varepsilon,$

where $\mu$ has $n$ unknown free parameters and $\varepsilon$ has iid entries
with mean zero and variance $\sigma^{2}$. So we have a sample size of 1 to
estimate each $\mu_{i}$, which is not enough if we also want to estimate
$\sigma^{2}$. An anonymous referee mentioned that if there are replicated rows
in $X$ (which implies ${\rm rank}(X)<n$), then one could estimate $\sigma^{2}$
using the measured responses corresponding to these the repeated rows.

We explore using cross-validation to choose a value of $\gamma$ that minimizes
the total validation squared error in our numerical experiments. This cross-
validation procedure implicitly assumes that the response and predictor
measurement pairs for each subject are drawn from a joint distribution. As an
alternative, we derive a high-dimensional estimator of $\gamma_{\rm opt}$ that
estimates $\sigma^{2}$ with an assumption about $\mu$. The following
paragraphs introduce this estimator, which is not invariant to invertible
linear transformations of $X$.

Recall that $\gamma_{\rm opt}={\delta^{2}}/({\sigma^{2}({\rm
rank}(X)-1)+\delta^{2}})$, where $\delta^{2}=\|\mu-P_{1}\mu\|^{2}$. Since
$P_{X}$ is an identity operator when ${\rm rank}(X)=n$,

$\displaystyle\mathbb{E}\|Y-P_{1}Y\|^{2}=\mathbb{E}\|P_{X}Y-P_{1}Y\|^{2}=\sigma^{2}({\rm
rank}(X)-1)+\delta^{2}.$

So given an estimator $\check{\sigma}^{2}$ of $\sigma^{2}$, we study the
following plug-in estimator of $\gamma_{\rm opt}$:

$\displaystyle\bar{\gamma}(\check{\sigma}^{2})={\rm
max}\Bigg{(}0,\frac{\|Y-P_{1}Y\|^{2}-\check{\sigma}^{2}({\rm
rank}(X)-1)}{\|Y-P_{1}Y\|^{2}}\Bigg{)},$ (11)

where the truncation at $0$ is necessary to ensure that
$\bar{\gamma}\in[0,1]$. We continue by describing the estimator of
$\sigma^{2}$ that we will use in (11).

If one ignores invariance and assumes Gaussian errors, then one could
simultaneously estimate $\beta$ and $\sigma^{2}$ by penalized likelihood with
the same penalty used in (2). However, this joint optimization is not convex.
We avoid this nonconvexity by modifying a reparametrized penalized Gaussian
likelihood optimization problem proposed by Zhu (2020). Let $\eta=\sigma^{-1}$
and $\beta^{*}=\beta\eta$. We estimate these parameters with

$\displaystyle(\hat{\beta}^{*},\hat{\eta})=\operatorname*{arg\
min}_{(\beta^{*},\eta)\in\mathbb{R}^{p}\times(0,\infty)}\left\\{\frac{1}{2n}\|Y\eta-X\beta^{*}\|^{2}-{\rm
log}(\eta)+\alpha\|\beta^{*}_{-1}\|^{2}\right\\},$ (12)

where
$\beta^{*}=(\beta^{*}_{1},\dots,\beta^{*}_{p})=(\beta^{*}_{1},\beta^{*}_{-1})$;
and $\alpha\geq 0$ be the tuning parameter for the Ridge penalty. This choice
of $\alpha$ was motivated by Liu et al. (2020), who verified that ridge
regression (with tuning parameter $\alpha$) can be used to consistently
estimate $\sigma^{2}$ provided that $\alpha\|\beta\|^{2}=o(1)$. However, Liu
et al. (2020) use a different estimator of $\sigma^{2}$ than the transformed
solution to (12). We also examine other choices for $\alpha$ in the
simulations (see section 6.4).

The reparametrized optimization problem in (12) is strongly convex with the
following global minimizer:

$\displaystyle\hat{\eta}$
$\displaystyle=\Big{(}n^{-1}Y^{\prime}(I-X(X^{\prime}X+2n\alpha
M)^{-1}X^{\prime})Y\Big{)}^{-1/2},$ (13) $\displaystyle\hat{\beta}^{*}$
$\displaystyle=\hat{\eta}(X^{\prime}X+2n\alpha M)^{-1}X^{\prime}Y,$

where $M={\rm diag}(0,1,1,\dots,1)\in\mathbb{R}^{p\times p}$. Since
$\eta=\sigma^{-1}$, which is estimated using (13), the corresponding estimator
of $\sigma^{2}$ is

$\displaystyle\check{\sigma}^{2}=n^{-1}Y^{\prime}(I-K)Y,$ (14)

where $K=X(X^{\prime}X+2n\alpha M)^{-1}X^{\prime}$. Using $\check{\sigma}^{2}$
in (11), we propose a high-dimensional estimator of $\gamma_{\rm pot}$ which
is given by

$\displaystyle\bar{\gamma}=\frac{Y^{\prime}(I-P_{1}-\frac{{\rm
rank}(X)-1}{{\rm rank}(X)}(I-K))Y}{Y^{\prime}(I-P_{1})Y},$ (15)

where truncation at $0$ is not required since the quadratic term on the
numerator of (15) is non-negative definite. Liu et al. (2020) proposed a bias
corrected version of $\check{\sigma}^{2}$, since the uncorrected estimator
does not converge to $\sigma^{2}$ under the assumptions for Theorem 1 of Liu
et al. (2020). Their corrected estimator is
$\check{\sigma}^{2}_{c}=C^{-1}\check{\sigma}^{2}$, where $C=1-{\rm tr}(K)/{\rm
rank}(X)$. Using $\check{\sigma}^{2}_{c}$ in (11), we also propose alternative
high-dimensional estimator of $\gamma_{\rm opt}$ by

$\displaystyle\bar{\gamma}_{c}={\rm
max}\Bigg{(}0,\frac{Y^{\prime}(I-P_{1}-\frac{{\rm rank}(X)-1}{{\rm
rank}(X)-{\rm tr}(K)}(I-K))Y}{Y^{\prime}(I-P_{1})Y}\Bigg{)}.$ (16)

First, we have the following consistency result for the corrected estimator
$\bar{\gamma}_{c}$.

######  Proposition 4.

Assume that the data-generating model in (1) is correct, that error
distribution has a finite fourth moment, $\delta^{2}=o(n)$, and
$d_{2}(n\alpha)^{-1}=o_{P}(1)$, where $d_{2}$ is the second-largest eigenvalue
of $XX^{\prime}$. Then $\bar{\gamma}_{c}-\gamma_{\rm opt}\rightarrow_{P}0$ as
$n\rightarrow\infty$.

The proof is in Section .4 of the Supplementary material. We see that
$\bar{\gamma}_{c}$ converges to $\gamma_{\rm opt}$ when $\gamma_{\rm
opt}\rightarrow 0$. The assumption that $d_{2}(n\alpha)^{-1}=o_{P}(1)$ is met
when $\alpha=n^{3/2}\left(2\|Y-P_{1}Y\|^{2}\right)^{-1}$, $\delta^{2}=o(n)$,
and $d_{2}=O(n)$ because

$\displaystyle d_{2}(n\alpha)^{-1}$
$\displaystyle=2n^{-5/2}d_{2}\|Y-P_{1}Y\|^{2}$
$\displaystyle=2n^{-5/2}d_{2}(\delta^{2}+2(\mu-
P_{1}\mu)^{\prime}\varepsilon+\varepsilon^{\prime}(I-P_{1})\varepsilon)\rightarrow_{P}0,$

since $\|\varepsilon\|^{2}=O_{P}(n)$.

We also established the following result for the uncorrected estimator
$\bar{\gamma}$.

######  Proposition 5.

Assume that the data-generating model in (1) is correct, that the error
distribution has a finite fourth moment, and at least one of the followings
holds:

  * •

$\delta^{2}=o(n)$ and $d_{2}(n\alpha)^{-1}=o_{P}(1)$,

  * •

$n=o(\delta^{2})$ and $n\alpha(d_{n})^{-1}=o_{P}(1)$,

where $d_{2}$ (resp. $d_{n}$) is the secondly largest (smallest) eigenvalue of
$XX^{\prime}$. Then $\bar{\gamma}-\gamma_{\rm opt}\rightarrow_{P}0$ as
$n\rightarrow\infty$.

The proof is in Section .5 of the Supplementary material. Using the
uncorrected error variance estimator $\check{\sigma}^{2}$ enables us to prove
consistency of $\bar{\gamma}$ when either $\gamma_{\rm opt}\rightarrow 0$ or
$\gamma_{\rm opt}\rightarrow 1$. When, $\gamma_{\rm opt}\rightarrow 0$, we
have consistency provided that $d_{2}=O(n)$ and
$\alpha=n^{t}\left(2\|Y-P_{1}Y\|^{2}\right)^{-1}$ with any $t\in(1,\infty)$.
We test $t=2,3$ in Section 6.4. On the other hand, when $\gamma_{\rm
opt}\rightarrow 1$, we have consistency provided that $d_{n}\asymp n$ and
$\alpha=n\left(2\|Y-P_{1}Y\|^{2}\right)^{-1}$.

######  Remark 1.

Azriel (2019) proved that “general” consistent estimation of the error
variance $\sigma^{2}$ in high-dimensional linear regression is impossible
without further assumptions in fixed-X settings. Their proof is based on a
Bayes risk approach that exploits a prior that requires $\delta^{2}\sim n$,
which is not what we assume in Proposition 5 ($\delta^{2}=o(n)$ or
$n=o(\delta^{2})$).

##  5 Shrinking toward the fitted values of a submodel

The previous sections developed our new shrinkage method for linear regression
with fitted values that are invariant to invertible linear transformations of
the design matrix $X$. Its shrinkage target was the fitted values of the
intercept-only model $P_{1}Y$. We can generalize this so that the shrinkage
target is $P_{X_{0}}Y$, where $X_{0}$ is the design matrix for a submodel
formed from a proper subset of the columns of $X$, e.g. $X_{0}=1_{n}$ (the
first column of $X$) as it was previously. We consider the following
generalization of (1):

$\displaystyle\operatorname*{arg\
min}_{b\in\mathbb{R}^{p}}\left\\{\gamma\|Y-Xb\|^{2}+(1-\gamma)\|Xb-
P_{X_{0}}Y\|^{2}\right\\},$

where $\gamma\in[0,1]$. This generalization could be useful when $\|\mu-
P_{X_{0}}\mu\|^{2}$ is substantially smaller than $\|\mu-P_{1}\mu\|^{2}$,
where $\mu=X\beta$. For example, a practitioner may know that some predictors
have large effect sizes, so these predictors would be encoded in $X_{0}$.

All of the results obtained for the special case that $X_{0}=1_{n}$ also hold
in this general case with $P_{1}Y$ replaced by $P_{X_{0}}Y$; and with ${\rm
rank}(X)-1$ replaced by ${\rm rank}(X)-{\rm rank}(X_{0})$. The proofs follow
by making these replacements in the proofs from the special case that
$X_{0}=1_{n}$.

We continue by stating these generalized results. The generalized fitted-value
shrinkage method’s fitted values are

$X\tilde{\beta}^{(\gamma)}=\gamma P_{X}Y+(1-\gamma)P_{X_{0}}Y.$

######  Proposition 6.

For all $(n,p)\in\\{1,2,\ldots\\}\times\\{1,2,\ldots\\}$,

$\mathbb{E}\|X\tilde{\beta}^{(\gamma)}-X\beta\|^{2}=\sigma^{2}\left\\{\gamma^{2}{\rm
rank}(X)+(1-\gamma^{2}){\rm rank}(X_{0})\right\\}+(\gamma-1)^{2}\|\mu-
P_{X_{0}}\mu\|^{2}.$ (17)

The right side of (17) is minimized when $\gamma=\tilde{\gamma}_{\rm opt}$,
where

$\displaystyle\tilde{\gamma}_{\rm opt}=\frac{\|\mu-
P_{X_{0}}\mu\|^{2}}{\sigma^{2}({\rm rank}(X)-{\rm rank}(X_{0}))+\|\mu-
P_{X_{0}}\mu\|^{2}}.$ (18)

We can estimate $\tilde{\gamma}_{\rm opt}$ in low dimensions with

$\displaystyle\tilde{\gamma}=(1-1/\tilde{F})\cdot 1(\tilde{F}>1),$ (19)

where $\tilde{F}$ is the F-statistic for comparing the submodel $X_{0}$ to the
full design matrix model $X$:

$\displaystyle\tilde{F}$
$\displaystyle=\frac{(\|Y-P_{X_{0}}Y\|^{2}-\|Y-P_{X}Y\|^{2})/({\rm
rank}(X)-{\rm rank}(X_{0}))}{\|Y-P_{X}Y\|^{2}/(n-p)}$
$\displaystyle=\frac{\|P_{X}Y-P_{X_{0}}Y\|^{2}}{\hat{\sigma}^{2}({\rm
rank}(X)-{\rm rank}(X_{0}))}.$ (20)

Recall that $r={\rm rank}(X)$ and let $r_{0}={\rm rank}(X_{0})$.

######  Proposition 7.

Assume that the data-generating model in (1) is correct, that the errors
follows a distribution with finite fourth moment, and $r-r_{0}\geq 1$. If
$p/n\rightarrow\tau\in[0,1)$ as $n\rightarrow\infty$ and either
$(r-r_{0})\rightarrow\infty$ or $\tilde{\delta}^{2}\rightarrow\infty$, then
$\tilde{\gamma}-\tilde{\gamma}_{\rm opt}\rightarrow_{P}0$.

######  Proposition 8.

Suppose that the assumptions of Proposition 7 hold, that the errors in (1) are
Gaussian, and that $r-r_{0}\geq 5$ is nondecreasing as $n\rightarrow\infty$.
Then,

$\tilde{\gamma}-\tilde{\gamma}_{\rm
opt}=\begin{cases}O_{P}((r-r_{0})^{-1/2})&\text{if
$\tilde{\delta}^{2}=O(r-r_{0})$}\\\
O_{P}((\tilde{\delta}^{2}/(r-r_{0}))^{-3/4})+O_{P}(n^{-1/2}(\tilde{\delta}^{2}/(r-r_{0}))^{-1/2})&\text{if
$(r-r_{0})\rightarrow\infty$ and $r-r_{0}=o(\tilde{\delta}^{2})$}\\\
O_{P}((\tilde{\delta}^{2})^{-3/4})+O_{P}(n^{-1/2}(\tilde{\delta}^{2})^{-1/2})&\text{if
$r-r_{0}=O(1)$.}\\\ \end{cases}$ (21)

##  6 Simulation studies

###  6.1 Low-dimensional experiments

We conducted a lower-dimensional simulation study in which the data were
generated from the linear regression subjects model (1) with $n=300$ and
$p\in\\{75,150,250\\}$. Also, $\epsilon_{1},\ldots,\epsilon_{n}$ are iid
$N(0,1)$. The design matrix $X$ has ones in its first column and independent
draws from $N_{p-1}(0,\Sigma)$ in the remaining entries on each row, where
$\Sigma_{jk}=0.5^{|j-k|}$. We randomly generated the regression coefficient
vector with the following equation:

$\beta={X}^{-}(1_{p}+\tau Z),$

where
$\tau\in\\{0,10^{-4},10^{-2},10^{-1},10^{-0.5},1,10^{0.5},10^{1},10^{1.5},10^{2}\\}$;
and $Z$ is $N_{p}(0,I)$. Then

$\displaystyle\mu-P_{1}\mu=(P_{X}-P_{1})(1_{p}+\tau Z)=\tau(P_{X}-P_{1})Z.$

So $\tau$ controls the size of $\delta^{2}=\|\mu-P_{1}\mu\|^{2}$.

We used 50 independent replications in each setting. In each replication, we
measured the performance of each estimator $\hat{\beta}_{\rm est}$ using the
same-X loss:

$\displaystyle n^{-1}\|X\beta-X\hat{\beta}_{\rm est}\|^{2}.$ (22)

The candidate estimators $\hat{\beta}_{\rm est}$ that we considered were the
following:

  * •

2n-G: $L_{2}$-squared penalty with 10-fold cross validation for $\gamma$ (3).

  * •

2n-Or: $L_{2}$-squared penalty using the oracle $\gamma_{\rm opt}$ in (5).

  * •

2n-Es: $L_{2}$-squared penalty using $\hat{\gamma}$ in (7).

  * •

2n-Es90: $L_{2}$-squared penalty using $\hat{\gamma}_{90}$ in (8).

  * •

2n-Es95: $L_{2}$-squared penalty using $\hat{\gamma}_{95}$ in (9).

  * •

2n-Rep: $L_{2}$-squared penalty using $\check{\sigma}^{2}$ in (14),
$\alpha=n(2\|Y-P_{1}Y\|^{2})^{-1}$ in (12), and the corresponding
$\bar{\gamma}_{\rm low}={\rm max}(0,1-1/F_{\rm rep})$, where $F_{\rm
rep}=\|P_{X}Y-P_{1}Y\|^{2}/({\check{\sigma}^{2}({\rm rank}(X)-1)})$.

  * •

O: Ordinary least square (OLS) estimator given by

$\displaystyle\hat{\beta}^{(1)}={X}^{-}{Y}.$ (23)

  * •

R: Ridge-penalized least squares (Hoerl and Kennard, 1970)

$\displaystyle\hat{\beta}_{\rm Ridge}={\rm
argmin}_{b\in\mathbb{R}^{p}}\|Y-Xb\|^{2}+\lambda\|b_{-1}\|^{2},$ (24)

where $b=(b_{1},b_{-1})$ with $b_{-1}=\mathbb{R}^{p-1}$; 10-fold cross
validation for the selection of $\lambda$.

  * •

L: Lasso-penalized least squares (Tibshirani, 1996)

$\displaystyle\hat{\beta}_{\rm LASSO}={\rm
argmin}_{b\in\mathbb{R}^{p}}\|Y-Xb\|^{2}+\lambda\|b_{-1}\|_{1},$ (25)

where 10-fold cross validation is used for the selection of $\lambda$.

For the methods that require cross validation, $\lambda$ and $\gamma$ were
selected from $\\{10^{-7+0.25j}:j=0,1,\cdots,44\\}$ and
$\\{\frac{k}{99}:k=0,1,\cdots,99\\}$, respectively. To facilitate the fairest
comparison between our invariant methods and the ridge/lasso methods, we used
the following standardization process, which is the default process used by
the R package glmnet: the ridge/lasso shrunken coefficient estimates are
computed using the standardized design matrix $X_{\bullet}$ defined by

$\displaystyle X_{\bullet}$
$\displaystyle=X\begin{bmatrix}1&-\bar{X}_{2}S_{2}^{-1}&-\bar{X}_{3}S_{3}^{-1}&\dots&-\bar{X}_{p-1}S_{p-1}^{-1}&-\bar{X}_{p}S_{p}^{-1}\\\
0&S_{2}^{-1}&0&\dots&0&0\\\ &0&S_{3}^{-1}&&\vdots&\vdots\\\
\vdots&\vdots&0&\ddots&0&\vdots\\\ 0&0&\vdots&&S_{p-1}^{-1}&0\\\
0&0&0&\dots&0&S_{p}^{-1}\\\ \end{bmatrix}$ $\displaystyle=XT,$

where $\bar{X}_{j}=n^{-1}\sum_{i=1}^{n}X_{ij}$ and
$S_{j}^{2}=(n-1)^{-1}\sum_{i=1}^{n}(X_{ij}-\bar{X}_{j})^{2}$ for
$j\in\\{2,\dots,p\\}$. Let $\hat{\beta}_{\bullet}$ be the shrinkage estimator
of the standardized coefficients. Since all the $S_{j}$’s will be positive, we
invert $T$ to estimate the original $\beta$ with
$T^{-1}\hat{\beta}_{\bullet}$. Our proposed fitted-value shrinkage procedures
are invariant to this standardizing transformation of $X$.

We display side-by-side boxplots of the same-X losses from the 50 replications
when $p=75$ in Figure 1. Further numerical summaries of these results are in
Table 4 in Section .2 of the Supplementary material. Without surprise, our
fitted-value shrinkage with oracle tuning 2n-Or performed the best among these
candidates. Our proposed estimator 2n-Es and its two variants 2n-Es95 and
2n-Es90 followed and generally outperformed OLS, Ridge, and Lasso when
$\tau\geq 10^{-0.5}$. Of the fitted-value shrinkage estimators, 2n-Es
outperformed 2n-Rep. On the other hand, the modified thresholds 2n-Es90 and
2n-Es95 performed better than 2n-Es for smaller values of $\tau$ (that
correspond to smaller values of $\delta^{2}$). Also, 2n-G outperformed Ridge,
Lasso, 2n-ES, and 2n-Es95 for smaller values of $\tau$. Ridge and Lasso
slightly outperformed 2n-Es when $\delta^{2}$ was small, but performed worse
when $\delta^{2}$ was larger.

In Figure 8 in Section .1 of the Supplementary material, we graphed the
average same-X loss values over the 50 replications as a function of
$\lambda$:

$\displaystyle f_{\rm FVS}(\lambda)=\|X\beta-X\hat{\beta}_{\lambda}\|^{2}/n,$
(26) $\displaystyle f_{\rm Ridge}(\lambda)=\|X\beta-X\hat{\beta}_{{\rm
Ridge};\lambda}\|^{2}/n,$ (27)

where $\hat{\beta}_{\lambda}$ and $\hat{\beta}_{{\rm Ridge};\lambda}$ are
solutions for (2) and (24); and
$\lambda\in\\{10^{-7+0.01k}:j=0,1,\ldots,1150\\}$. This allows a further
comparison of fitted-value shrinkage (3) to Ridge regression (24). Further
numerical summaries of these results are in Table 10 in Section .2 of the
Supplementary material. The minimum average same-X loss for fitted value
shrinkage (2) was either less than or nearly equal to that of Ridge (24).
However, the range of values of $\lambda$ that corresponded to average losses
near the minima was much narrower for fitted-value shrinkage than it was for
Ridge when medium to large values of $\tau$ were used.

We also display side-by-side boxplots of the observed same-X losses from the
50 replications when $p=150$ in Figure 2 and when $p=250$ in Figure 7 in
Section .1 of the Supplementary material. In addition, graphs of the average
loss values as a function of $\lambda$ over
$\lambda\in\\{10^{-7+0.01k}:j=0,1,\ldots,1100\\}$ when $p=150$ and $p=250$ are
in Figure 9–10 in Section .1 of the Supplementary material. Further numerical
summaries for this simulation with $p\in\\{150,250\\}$ are in Table 5–6, and
Table 11–12 in Section .2 of the Supplementary material. These results with
$p\in\\{150,250\\}$ are similar to results when $p=75$: the oracle method
2n-Or was the best and our proposed estimators 2n-Es90, 2n-Es95, 2n-Es were
the most competitive when $\tau\geq 1$. However, the performance gap between
our procedure with non-oracle tuning 2n-G and our procedure using oracle
tuning 2n-Es when $\tau\geq 1$ has increased (Figure 2(c), 2(d)). We expect
this is related to the narrower valley observed in the graph of the average
same-X loss as a function of the tuning parameter. As it was when $p=75$, 2n-G
was the most competitive for smaller values of $\tau$.

In Tables 7–8 in Section .2 of the Supplementary material, we report 95%
simulation-based confidence intervals for the expected squared same-X loss
difference (fitted value shrinkage minus Ridge). The confidence intervals are
based on 50 replications for each pair of $(\tau,p)$. These differences were
statistically insignificant when $\tau$ was small, but were significant when
$\tau$ was larger (with fitted value shrinkage outperforming Ridge). In Table
9 in Section .2 of the Supplementary material, we also report the average
realization of $|\hat{\gamma}-\gamma_{\rm opt}|^{2}$ based on 100 independent
replications from the same simulation setting. This quantity increased with
$p$, but remained stable.

(a) $\tau=10^{-2}$

(b) $\tau=10^{-1}$

(c) $\tau=1$

(d) $\tau=10^{0.5}$

Figure 1: Boxplots of the observed same-X loss values from the 50 replications
when $n=300$ and $p=75$. We suppress “2n” from the third to the last
estimators to simplify notation.

(a) $\tau=10^{-2}$

(b) $\tau=10^{-1}$

(c) $\tau=1$

(d) $\tau=10^{0.5}$

Figure 2:  Boxplots of the observed same-X loss values from the 50
replications when $n=300$ and $p=150$. See the caption of Figure 1 for more
details.

###  6.2 Impact of invariance I: In the presence of categorical variables with
different reference level coding schemes

We generated $X\in\mathbb{R}^{n\times p}$ to have $25$ numerical predictors
and $3$ categorical predictors, where each categorical predictor had 5 levels,
and the 25th numerical predictor had interactions with the three categorical
predictors. The design matrix $X$ had $p=1+25+3\times 4\times 2=50$ columns
including the intercept. The $n=100$ observations of the $25$ numerical
predictors were independent draws from $N_{25}(0,\Sigma)$, where
$\Sigma_{jk}=0.5^{|j-k|}$. Observations of each categorical predictor were
independently drawn from a 5-category multinomial distribution with equal
category probabilities. As they were before,
$\epsilon_{1},\ldots,\epsilon_{n}$ are iid $N(0,1)$.

We organize the regression coefficient vector as
$\beta=(\beta_{c},\beta_{f})\in\mathbb{R}^{50}$, where $\beta_{c}$ corresponds
to the first $26$ columns of $X$, which have the numerical predictors; and
$\beta_{f}$ corresponds to the remaining 24 columns of $X$, which have the
categorical predictors and their interaction with the 25th numerical
predictor. Let $X_{c}\in\mathbb{R}^{100\times 26}$ be the matrix with the
first 26 columns of $X$. We generated
$\beta_{c}=X_{c}^{-}(1_{26}+\tau_{c}Z_{c})$; where $Z_{c}\sim N_{26}(0,I)$. We
also generated

$\displaystyle\beta_{f}=(\tau_{f}X_{f}^{-}Z_{f})*v,$ $\displaystyle
v=(w^{\prime},w^{\prime})^{\prime}\in\mathbb{R}^{24},$ $\displaystyle
w=(1,2,0,0,0,2,1,0,0,0,2,1)^{\prime}\in\mathbb{R}^{12},$

where $*$ denotes the elementwise product of two vectors of same dimension;
$Z_{f}\sim N_{24}(0,I)$; $X_{f}\sim N_{24}(0,\Sigma)$ with
$\Sigma_{jk}=0.5^{|j-k|}$; $Z_{c}$ is independent of $Z_{f}$; and $X_{c}$ is
independent of $X_{f}$. We consider $(\tau_{c},\tau_{f})$ in

$\displaystyle\\{(10^{-0.5},10^{-0.5}),(1,10^{-0.5}),(10^{-0.5},1),(1,1),(10^{0.5},1),(1,10^{0.5}),(10^{0.5},10^{0.5})\\}.$

The first 12 entries of $\beta_{f}$ correspond to the $3*4$ main effects
columns for the 3 categorical predictors. The remaining 12 entries of
$\beta_{f}$ correspond to interactions between the 3 categorical predictors
and the 25th numerical predictor.

The first level of each categorical predictor was coded as the reference level
in this design matrix $X$. We call this “Coding-1”. To illustrate the effects
of a coding change, we will also use a “Coding-2” design matrix $X_{\bullet}$,
which is an invertible linear transformation of $X$ that uses the second,
third, and fifth levels of the three categorical predictors as the reference
levels for the first, second, and third categorical predictor, respectively.

The competitors were 2n-Es, 2n-Es95, 2n-Or, 2n-G, 2n-Rep, O, R, L, and Group
Lasso and its two variants. The Group lasso competitors were the following:

  * •

GL: Group Lasso estimator proposed by Yuan and Lin (2006).

$\displaystyle\hat{\beta}_{\rm GL}=\operatorname*{arg\
min}_{b\in\mathbb{R}^{p}}\|Y-Xb\|^{2}+\lambda\sum_{i=1}^{L}\sqrt{p_{i}}\|b^{(i)}\|_{2},$

where $i\in\\{1,\ldots,L\\}$ is the index of each group except for one
intercept column; $p_{i}$ is the size of the $i$-th group; and $b^{(i)}$ is
the corresponding subvector.

  * •

MGL: Modified Group Lasso proposed by Choi et al. (2012).

$\displaystyle\hat{\beta}_{\rm MGL}=\operatorname*{arg\
min}_{b\in\mathbb{R}^{p}}\|Y-Xb\|^{2}+\lambda\sum_{i=1}^{L}\|b^{(i)}\|_{2}.$

  * •

SGL: Standardized Group Lasso proposed by Simon and Tibshirani (2012).

$\displaystyle\hat{\beta}_{\rm SGL}=\operatorname*{arg\
min}_{b\in\mathbb{R}^{p}}\|Y-Xb\|^{2}+\lambda\sum_{i=1}^{L}\sqrt{p_{i}}\|X^{(i)}b^{(i)}\|_{2},$

where $X^{(i)}$ is the $i$-th submatrix.

For the estimators that require tuning parameter selection, 10-fold cross
validation was used. For GL, MGL, and SGL, we formed 31 groups: 25 single-
member groups corresponding to the numerical predictors, 3 groups
corresponding to the categorical predictors’ main effects, and 3 groups
corresponding to the interaction between the 25th numerical predictor and the
three categorical predictors.

We display side-by-side boxplots of the same-X loss values from 50
replications in Figure 3. Further numerical summaries of these results are in
Table 13 in Section .3 of the Supplementary material. Except for
$(\tau_{c},\tau_{f})=(10^{0.5},10^{0.5})$ and $(10^{0.5},10^{-0.5})$, 2n-Or
was the best in Coding-2. Generally, 2n-Es was second-best in Coding-2. The
methods R, L, GL, MGL and SGL generally performed well in Coding-1, which is
not surprising considering the sparsity in the coefficient vector; however, in
Coding-2, these methods were generally worse than 2n-Es, which is invariant to
this change in coding. This illustrates that the invariance of our proposed
estimators is advantageous.

(a) $(10^{-0.5},1,1)$

(b) $(10^{-0.5},1,2)$

(c) $(1,1,1)$

(d) $(1,1,2)$

(e) $(10^{0.5},1,1)$

(f) $(10^{0.5},1,2)$

Figure 3: Boxplots of the observed same-X loss values from the 50 replications
in the simulation settings of Section 6.2. We suppress “2n” from the third to
the last estimators for notational simplicity. Each plot is labeled with the
$(\tau_{c},\tau_{f},{\rm Coding~{}number})$ value used.

###  6.3 Impact of invariance II: General full-rank transformation of the
model matrix (Low dimension)

In this example, the design matrix $X\in\mathbb{R}^{n\times p}$ has ones in
its first column and its $n$ row vectors (excluding the first entry) are drawn
independently from $N_{p-1}(0,\Sigma)$ where $\Sigma_{jk}=0.5^{|j-k|}$ and
$(n,p)=(300,150)$. We generated the regression coefficient vector $\beta$ as
$\beta=u*v$, where each element of $u\in\mathbb{R}^{p}$ is an independent draw
from the Uniform distribution on $(2^{-\psi-1},2^{-\psi})$ with
$\psi\in\\{0,1\\}$; and $v=(1,v_{-1})\in\mathbb{R}^{p}$, where each element of
$v_{-1}\in\mathbb{R}^{p-1}$ is an independent draw from ${\rm Ber}(s)$ with
$s\in\\{0.025,0.05,0.1,0.2,0.3\\}$. We generated $Y$ by (1).

We denote the model fitting scheme with this $X$ as “Coding-1”. We also use a
transformed design matrix $X_{\bullet}=XT$, where $T$ is a Gram-Schmidt
orthogonalization of a matrix with all of its entries drawn independently from
$N(0,1)$. We refer the transformed design matrix $X_{\bullet}$ as “Coding-2”.

In Figure 4, we display side-by-side boxplots of the observed same-X losses
from the 50 replications from representative settings when $\psi=0$. Further
numerical summaries of the results are in Table 14 in Section .4 of the
Supplementary material. When $\tau=1$, Lasso was the best in Coding-1 with
Ridge following next. However, in Coding-2, 2n-Or was the best and 2n-Es was
the second-best regardless of $\psi$ and $s$. Ridge and Lasso were
significantly worse than our proposed estimator as well as OLS under Coding 2.
This further illustrates the benefits of invariance.

(a) $(0.1,1)$

(b) $(0.1,2)$

(c) $(0.2,1)$

(d) $(0.2,2)$

(e) $(0.3,1)$

(f) $(0.3,2)$

Figure 4: Boxplots of the observed same-X loss values from the 50 replications
in the simulation settings of Section 6.3 when $\psi=0$. We suppress “2n” from
the third to the last estimators for notational simplicity. Each plot is
labeled with the $(s,{\rm Coding~{}number})$ value used.

###  6.4 High-dimensional experiments

We used the same data generating model as in Section 6.1 except that $n=200$,
$p=300$, and $\sigma\in\\{2,3\\}$. Since 2n-Es is not applicable in high
dimensions, we tested variants of 2n-G that used 5, 10, and $n$-fold cross
validation. However, except for a few cases, the number of folds used for
tuning parameter selection for 2n-G did not have a significant impact on the
performance. So we only present results from 10-fold cross validation with the
same label 2n-G as in the previous sections. In addition, we tried different
$\alpha$ values that control the matrix $K$ in (16) for the following variants
of 2n-Rep:

  * •

2n-Rep1: Same as original 2n-Rep with $\alpha=n(2\|Y-P_{1}Y\|^{2})^{-1}$.

  * •

2n-Rep2: $\bar{\gamma}_{c}$ (16) using $\check{\sigma}^{2}$ with
$\alpha=n^{3/2}(2\|Y-P_{1}Y\|^{2})^{-1}$.

  * •

2n-Rep3: Same as original 2n-Rep with $\alpha=n^{2}(2\|Y-P_{1}Y\|^{2})^{-1}$.

  * •

2n-Rep4: Same as original 2n-Rep with $\alpha=n^{3}(2\|Y-P_{1}Y\|^{2})^{-1}$.

In Figure 5, we display side-by-side boxplots of the observed same-X losses
from the 50 replications when $n=200$ and $p=300$. There are additional
boxplots displayed in Figure 11 in Section .1 of the Supplementary material.
Further numerical summaries of the results are in Table 15–17 in Section .5 of
the Supplementary material. In general, when $\delta^{2}=\|\mu-P_{1}\mu\|^{2}$
was large relative to $n\sigma^{2}$, which corresponds to the situation
$\gamma_{\rm opt}\rightarrow 1$, 2n-Rep1 performed substantially better than
Ridge, Lasso, 2n-G, 2n-Rep2, and 2n-Rep4 (Figure 5(b), 5(d)). The method
2n-Rep1 also outperformed 2n-Rep3, but the same-X prediction difference was
relatively smaller than the others. Furthermore, larger $\delta^{2}$ led to
improved tuning-parameter selection for 2n-Rep1. However, 2n-G, Ridge, and
Lasso performed better when $\delta^{2}$ was small relative to $n\sigma^{2}$.
In this setting, which corresponds to $\gamma_{\rm opt}\rightarrow 0$,
2n-Rep2, 2n-Rep3, and 2n-Rep4 performed as well as or better than Ridge,
Lasso, and 2n-G (see Figure 5(a), 5(c), 5(e)). On the other hand, 2n-Rep1
struggled for this case. The method 2n-Rep3 was relatively more stable than
the other 2n-Rep versions, and it consistently outperformed Ridge and Lasso
regardless of $\delta^{2}$.

In Figure 12 (see Section .1 of the Supplementary material), we display
average loss values over the 50 replications as a function of $\lambda$ over
$\lambda\in\\{10^{-7+0.01k}:j=0,1,\ldots,1100\\}$ when $n=200$ and $p=300$
(see (26) and (27)). Further numerical summaries of these results are in Table
18 in Section .5 of the Supplementary material. These results look similar to
lower-dimensional results displayed in Figure 8–10 (see Section .1 of the
Supplementary material), except the curve valleys are narrower for fitted-
value shrinkage.

(a) $(1,10^{-0.5})$

(b) $(1,10^{0.5})$

(c) $(2,10^{-0.5})$

(d) $(2,10^{0.5})$

(e) $(3,10^{-0.5})$

(f) $(3,10^{0.5})$

Figure 5:  Boxplots of the observed same-X loss values from the 50
replications when $(n,p)=(200,300)$. Each plot is labeled with the
$(\sigma,\tau)$ value used.

###  6.5 Impact of (approximate) invariance III: General full-rank
transformation of the design matrix (High dimension)

We again generated $X\in\mathbb{R}^{n\times p}$ to have ones in its first
column and and its $n$ row vectors (excluding the first entry) were drawn
independently from $N_{p-1}(0,\Sigma)$ where $\Sigma_{jk}=0.5^{|j-k|}$ and
$(n,p)=(200,300)$. We used (1) again for the data generating procedure, where
we consider $\sigma\in\\{1,3\\}$.

We generated $\beta=u*v$, where each element of $u\in\mathbb{R}^{p}$ is an
independent draw from the Uniform distribution on $(1/4,1/2)$ when $\sigma=1$;
and on $(1/6,1/3)$ when $\sigma=3$; $v=(1,v_{-1})\in\mathbb{R}^{p}$, where
each element of $v_{-1}\in\mathbb{R}^{p-1}$ is an independent draw from ${\rm
Ber}(s)$. When $\sigma=1$, we vary $s\in\\{0.025,0.05,0.1,0.2,0.3\\}$, and
consider 2n-Rep1 for the proposed estimator. On the other hand, when
$\sigma=3$, we vary $s\in\\{0.01,0.02,0.03,0.04,0.05\\}$, and consider 2n-Rep3
for the proposed estimator as well as 2n-Rep2, which is known to be efficient
in the same setting that 2n-Rep3 is (see Proposition 4). We dropped 2n-Rep4,
since it showed relatively unstable performance in comparison to 2n-Rep3 in
Section 6.4.

As in Section 6.3, we denote the model fitting scheme with this $X$ as
“Coding-1”. We also use a transformed design matrix $X_{\bullet}=XT$, where
$T$ through a Gram-Schmidt orthogonalization of a matrix with all of its
entries drawn independently from $N(0,1)$. We refer the transformed design
matrix $X_{\bullet}$ as “Coding-2”.

In Figure 6, we display side-by-side boxplots of the observed same-X losses
from the 50 replications in few representative settings. Further numerical
summaries of the simulation results are in Table 19–20 in Section .6 of the
Supplementary material. When $\tau=1$, Lasso was generally the best in
Coding-1 and Ridge was second best. However, in Coding-2, 2n-Or was the best
and 2n-Rep1 was second best, regardless of $\psi$ and $s$. Ridge and Lasso
performed significantly worse than 2n-Rep1 after transformation.

On the other hand, when $\sigma=3$, 2n-Or achieved the best performance for
both Coding 1 and 2 except for $s=0.05$ where Ridge was the best for both
Coding 1 and 2. For $s\in\\{0.01,0.02,0.03,0.04\\}$, 2n-Rep3 was competitive,
so were 2n-G and 2n-Rep2, which is similar to the pattern observed in Section
6.4 when $\delta^{2}$ was low.

Although the 2n-Rep methods lack exact invariance in high dimensions, their
performance difference between Coding 1 and Coding 2 was not significant.

(a) $(1,0.1,1)$

(b) $(1,0.1,2)$

(c) $(1,0.2,1)$

(d) $(1,0.2,2)$

(e) $(3,0.03,1)$

(f) $(3,0.03,2)$

(g) $(3,0.04,1)$

(h) $(3,0.04,2)$

Figure 6: Boxplots of the observed same-X loss values from the 50 replications
in the simulation settings of Section 6.5. We suppress “2n” from the third to
the last estimators for notational simplicity. Each plot is labeled with the
$(\sigma,s,{\rm Coding~{}number})$ value used.

##  7 Data examples

###  7.1 Low dimensional data experiments

We compared our proposed fitted-value shrinkage procedures to competitors on
three data sets. We used the same non-oracle estimators as the previous
section except we excluded 2n-Es90 because it performed similarly to 2n-Es95.
Each data example was analyzed using the following procedure: For 50
independent replications, we randomly selected 70% of the subjects for the
training set and used the remaining subjects as the test set. Tuning parameter
selection was done using the training set and prediction performance was
measured using squared error loss on the test set. The following is the short
description of the three low-dimensional data set we examined.

  1. (FF)

: The Forest Fire (FF) data are from from Cortez and Morais (2007) and are
stored at the UCI Machine learning repository via
https://archive.ics.uci.edu/dataset/162/forest+fires. There are 517
observations corresponding to forest fires in Portugal from 2000 to 2003. The
response is the total burned area (in $ha$) from the fire, which was
transformed with $x\mapsto{\rm ln}(x+1)$, which was suggested by Cortez and
Morais (2007). There were originally 13 attributes. However, in the pre-
processing step, since we are not focusing on spatio-temporal methods, we
excluded time, date and location coordinates. After this processing, the full-
data design matrix had $(n,p)=(517,9)$ with 8 numerical-variable columns and
one intercept column.

  2. (GDP)

: The GDP data (GDP) are from Barro and Lee (1994). These data consist of 161
observations of GDP growth rates for the two periods 1965-1975 and 1975-1985.
The data are also in the R package quantreg (Koenker, 2022). The response is
Annual change per capita GDP. There are 13 numerical predictors, e.g. Initial
per capita GDP, Life expectancy. We also added a quadratic term for the
predictor Black Market Premium. After processing, the full-data design matrix
has $(n,p)=(161,15)$.

  3. (FC)

: The Forecast data set (FC) is from Cho et al. (2020) for the purpose of bias
correction for the Local Data Assimilation and Prediction System (LDAPS),
which is a numerical weather report model used by Korea Administration (KMA),
Seoul, South Korea. It has a public access through
https://archive.ics.uci.edu/dataset/514/bias+correction+of+numerical+prediction+model+temperature+forecast.
The data are regional observations from 2013 to 2017, from which we randomly
selected 500. We used the true maximal temperature of the next day as the
response and removed date, station ID, and true minimal temperature of the
next day. The full-data design matrix had $(n,p)=(500,20)$.

In Table 1, we display mean squared prediction errors averaged over 50
training/test set splits for the three data examples (FF), (GDP), and (FC).
Our fitted-value shrinkage estimators performed similarly to Ridge and Lasso,
which both lack invariance to invertible linear transformations of the design
matrix.

Performance comparison table
---
Data set | 2n-Rep | 2n-G | 2n-Es | 2n-Es95 | OLS | Ridge | LASSO
Forest Fire | 2.0031 | 2.0160 | 2.0017 | 1.9792 | 2.1834 | 2.0189 | 2.0572
(FF) | (0.0301) | (0.0316) | (0.0300) | (0.0293) | (0.0551) | (0.0322) | (0.0411)
GDP growth | 3.106e-04 | 3.125e-04 | 3.139e-04 | 3.139e-04 | 3.231e-04 | 3.112e-04 | 3.149e-04
(GDP) | (7.465e-06) | (7.615e-06) | (7.706e-06) | (7.706e-06) | (8.120e-06) | (8.126e-06) | (8.307e-06)
Forecast | 1.9045 | 1.8997 | 1.9028 | 1.8933 | 2.0489 | 1.8990 | 1.8978
(FC) | (0.0327) | (0.0325) | (0.0327) | (0.0318) | (0.0371) | (0.0321) | (0.0318)

Table 1: The performance comparison table for three data examples in Section
7.1. The values are the mean squared prediction errors averaged over 50
training/test set splits. The numbers in parentheses are normalized sample
standard deviations. The column labels are defined in Section 6.1. Boldface
indicates the best model. Underlined is our main proposed estimator.

###  7.2 Low dimensional data analyses with categorical variables and their
interactions

We analyzed two data sets from existing R packages to illustrate the
performance of our estimators when categorical variables with interactions are
present in the model. The competitors and setup for the data experiments are
nearly identical to the previous Section 6.2, except we added three new
fitted-value shrinkage estimators that shrink toward the submodel without
interactions instead of the intercept-only model. These new submodel shrinkage
methods are labeled 2n-Repsb, 2n-Gsb, 2n-Essb, 2n-Es95sb, and they
respectively correspond to 2n-Rep, 2n-G, 2n-Es, and 2n-Es95.

The following is a description of the data examples:

  1. (Dia-1)

: The Diamonds data set is from Diamonds data frame in the R package Stat2Data
(Cannon et al., 2019), and it was obtained from https://awesomegems.com/. The
are $n=351$ subjects and the response is the price of the diamond (in
dollars). We divided price per 1000 and used it for the response variable. We
further removed total price from the predictors. There are 2 categorical
predictors: color (with levels D to J) and clarity (with levels IF, VVS1,
VVS2, VS1, VS2, SI1, SI2, and SI3). We divided color into 5 levels (D, E, F,
G, and (H,I,J)), and categorized the clarity into 3 levels ((IF, VVS1, VVS2),
(VS1, VS2), (SI1, SI2, SI3)). We used reference-level coding in the design
matrix, where (H,I,J) was the reference level for color; and (SI1, SI2, SI3)
was the reference level for clarity. Interactions between color and depth as
well as clarity and depth were added. The full design matrix has
$(n,p)=(351,15)$ and the submodel with linear terms only has $p=9$.

  2. (Dia-2)

: The setting is identical to that of (Dia-1), except we used the category
(VS1, VS2) as the reference level for coding the categorical predictor
clarity.

  3. (NG-1)

: The NaturalGas data is from Baltagi (2002), and is in the R package AER
(Kleiber and Zeileis, 2008). There are 138 observations on 10 variables. We
removed state name and year and added an interaction between state code and
heating degree days. We set the response as consumption divided by 10000. The
reference level for state code, which is the only categorical predictor, was
set to 35 (NY). The full-data design matrix had $(n,p)=(138,17)$ and the
submodel without interactions had $p=12$.

  4. (NG-2)

: This is the same as (NG-1), except the reference level for state code was
set to 5 (CA).

We display mean squared prediction errors averaged over 50 training/test set
splits for these data examples in Table 2. Our proposed estimators performed
similarly or better than Ridge and Lasso. We also notice that changing the way
that categorical predictors were encoded in the design matrix changes the
performance of Ridge and Lasso, which lack invariance. Furthermore, Group
Lasso and its two standardized versions had unstable performance when the base
level coding was changed, which does not happen to our invariant methods.

Performance comparison table
---
Data set | 2n-Rep | 2n-Repsb | 2n-G | 2n-Gsb | 2n-Es | 2n-Essb | 2n-Es95 | 2n-Es95sb | OLS | Ridge | LASSO | GL | MGL | SGL
Diamonds | 1.2538 | 1.2557 | 1.2545 | 1.2501 | 1.2557 | 1.2507 | 1.2557 | 1.2599 | 1.2587 | 1.2668 | 1.2675 | 1.2699 | 1.2693 | 1.2569
(Dia-1) | (0.0303) | (0.0300) | (0.0298) | (0.0293) | (0.0295) | (0.0294) | (0.0295) | (0.0294) | (0.0293) | (0.0297) | (0.0297) | (0.0290) | (0.0290) | (0.0301)
Diamonds | 1.2538 | 1.2557 | 1.2545 | 1.2501 | 1.2557 | 1.2507 | 1.2557 | 1.2599 | 1.2587 | 1.2615 | 1.2624 | 1.2697 | 1.2677 | 1.2592
(Dia-2) | (0.0303) | (0.0300) | (0.0298) | (0.0293) | (0.0295) | (0.0294) | (0.0295) | (0.0294) | (0.0293) | (0.0296) | (0.0299) | (0.0286) | (0.0287) | (0.0296)
Natural gas | 5.2558 | 5.2485 | 5.2899 | 5.2711 | 5.2559 | 5.2523 | 5.2559 | 5.2523 | 5.2545 | 5.2913 | 5.1254 | 5.2972 | 5.3010 | 5.1584
(NG-1) | (0.2285) | (0.2262) | (0.2282) | (0.2217) | (0.2285) | (0.2257) | (0.2285) | (0.2257) | (0.2279) | (0.2195) | (0.2285) | (0.2241) | (0.2242) | (0.2216)
Natural gas | 5.2558 | 5.2486 | 5.2899 | 5.2711 | 5.2559 | 5.2523 | 5.2559 | 5.2523 | 5.2545 | 5.3563 | 5.2889 | 5.2261 | 5.2010 | 6.7549
(NG-2) | (0.2285) | (0.2262) | (0.2282) | (0.2217) | (0.2285) | (0.2257) | (0.2285) | (0.2257) | (0.2279) | (0.2195) | (0.2285) | (0.2171) | (0.2168) | (0.2877)

Table 2: The performance comparison table for two data sets in Section 7.2
with two different coding strategies. The values are mean squared prediction
errors averaged over 50 training/test splits. The numbers in parentheses are
normalized sample standard deviations. The labels are defined in Section 6.1
except that last three are illustrated in Section 7.2. The suffix ”sb” for
each estimator corresponds to shrinking towards the submodel without
interactions instead of shrinking towards the intercept-only model (see
Section 5 for details). Boldface indicates the best model. Underlined are our
main estimator and its submodel shrinkage variant.

###  7.3 High dimensional data experiments

For high-dimensional data examples, we randomly selected subjects from
existing data sets so that there were fewer subjects than predictors. We used
the same splitting and evaluation procedure that we used in Sections 7.1 and
7.2. The competitors are same as those considered in Section 6.4 except that
we excluded 2n-Rep4 which had nearly identical performance to 2n-Rep3 in the
simulation study. The following is a description of the examples:

  1. (mtp)

: The data set mtp comes from Karthikeyan et al. (2005) and is available at
the OpenML repository via
https://www.openml.org/search?type=data&status=active&id=405. There are 4450
subjects with 203 numerical measurements. The response is oz203. We randomly
selected 120 subjects and removed the 23 predictors that had fewer than 30
distinct values, which ensured that there were no constant columns in the
120-row design matrix other than the intercept column. The full-data design
matrix had $(n,p)=(120,180)$.

  2. (topo)

: The topo.2.1 data set is from Feng et al. (2003) and is available through
the OpenML repository via
https://www.openml.org/search?type=data&sort=runs&id=422&status=active. There
are 8885 subjects with 267 numerical measurements. The response is oz267. We
randomly selected 180 subjects and removed the 22 predictors that had fewer
than 30 distinct values. After this, there were 34 constant columns (other
than the intercept) that were also removed. The R code for this processing is
in Section .0.1 of the Supplementary material. The full-data design matrix has
$(n,p)=(180,214)$.

  3. (tecator)

: The tecator data set comes from Thodberg (2015), and is also available from
OpenML repository via
https://www.openml.org/search?type=data&status=active&sort=runs&order=desc&id=505.
We randomly selected 100 subjects and removed 22 principal components. The
response is fat content. The full-data design matrix had $(n,p)=(100,103)$.

In Table 3, we report the mean squared prediction errors averaged over 50
training/test set splits for these data examples. We see that 2n-Rep2, 2n-Rep3
had similar prediction performance compared to 2n-Rep1. In contrast to its
same-X loss performance in simulations, the cross validation version of our
method 2n-G gave reasonable out-of-sample prediction performance. Generally,
2n-Rep1 and 2n-Rep2 performed competitively compared with Ridge and Lasso, and
2n-Rep3 followed the next within 2n-Reps.

Performance comparison table
---
Data set | OLS | 2n-G | Ridge | LASSO | 2n-Rep1 | 2n-Rep2 | 2n-Rep3
mtp | 0.5075 | 0.0266 | 0.0262 | 0.0227 | 0.0257 | 0.0261 | 0.0258
(mtp) | (0.0592) | (0.0010) | (0.0027) | (0.0017) | (0.0008) | (0.0009) | (0.0008)
topo.2.1 | 0.07524 | 0.00085 | 0.00096 | 0.00087 | 0.00084 | 0.00083 | 0.00084
(topo) | (2.72e-02) | (2.75e-05) | (7.71e-05) | (3.27e-05) | (3.08e-05) | (2.73e-05) | (3.08e-05)
tecator | 2.4400 | 2.4592 | 2.7251 | 2.7120 | 2.4266 | 2.4318 | 2.4733
(tecator) | (0.2359) | (0.2348) | (0.1913) | (0.1940) | (0.2202) | (0.2188) | (0.2148)

Table 3: The performance comparison table for three high dimensional real data
sets in Section 7.3. The values are mean squared prediction errors averaged
over 50 replications. The numbers in parentheses are normalized sample
standard deviations. The labels are defined in Section 6.4. Boldface indicates
the best model.

##  8 Future directions

The fitted-value shrinkage idea presented here to fit linear regression models
with invariance can be extended to more complicated settings. For example, one
could fit a logistic regression model by minimizing the negative loglikelihood
plus the penalty $\lambda\|X\beta-\hat{w}1_{n}\|^{2}$, where $\hat{w}$ is the
sample log-odds that the response takes its first category. We are currently
developing this procedure.

A Bayesian formulation of our method may also be interesting. In addition, one
could study methods that combine our proposed invariant shrinkage penalty with
regular shrinkage penalties like the lasso or ridge.

## Acknowledgement

The authors thank the editor, associate editor, and referees for their helpful
comments and suggestions.

Supplementary material for Fitted value shrinkage

## Technical details

###  .1 Proof of Proposition 1

###### Proof.

We start with the following standard decomposition:

$\displaystyle\mathbb{E}\|X\hat{\beta}^{(\gamma)}-X\beta\|^{2}={\rm
tr}\Big{(}{\rm
var}\Big{(}X\hat{\beta}^{(\gamma)}-X\beta\Big{)}\Big{)}+\Big{\|}\mathbb{E}\Big{(}X\hat{\beta}^{(\gamma)}-X\beta\Big{)}\Big{\|}^{2}.$

From (4),

$\displaystyle{\rm tr}\Big{(}{\rm
var}\Big{(}X\hat{\beta}^{(\gamma)}-X\beta\Big{)}\Big{)}$ $\displaystyle={\rm
tr}\Big{(}\sigma^{2}\Big{(}\gamma^{2}P_{X}+(1-\gamma^{2})P_{1}\Big{)}\Big{)},$
$\displaystyle=\sigma^{2}(\gamma^{2}r+1-\gamma^{2})$
$\displaystyle\mathbb{E}\Big{(}X\hat{\beta}^{(\gamma)}-X\beta\Big{)}$
$\displaystyle=(1-\gamma)(\mu-P_{1}\mu),$

where $\mu=X\beta$. Combining above two equations, we conclude that the
statement holds. ∎

###  .2 Proof of Proposition 2

######  Lemma 1.

When a $n$-variate random variable $V=(V_{1},\dots,V_{n})$ has i.i.d. elements
$V_{i}$ following a distribution that has mean $\mu$, variance $\sigma^{2}$,
and has finite $4$-th moment with $\mathbb{E}(V_{i}^{4})\leq M$, for all
$i=1,\dots,n$. Let an $(i,j)$-th element of $A$ be $a_{ij}$, then, for a
symmetric non-negative definite $A\in\mathbb{R}^{n\times n}$,

$\displaystyle{\rm var}(V^{\prime}AV)\leq
2(M-3\sigma^{4})\sum_{i=1}^{n}a_{ii}^{2}+4\sigma^{4}{\rm
tr}(A^{2})+8\sigma^{2}\mu^{\prime}A^{2}\mu$

holds.

###### Proof.

For the simplest case, we first consider when $\mu=0$. Then, by the i.i.d
property of $\\{V_{i}\\}$ and the moment conditions,

$\displaystyle\mathbb{E}\Big{[}(V^{\prime}AV)^{2}\Big{]}$ $\displaystyle\leq
M\sum_{i=1}^{n}a_{ii}^{2}+\sigma^{4}\sum\sum_{1\leq i\neq j\leq
n}a_{ii}a_{jj}+2\sigma^{4}\sum\sum_{1\leq i\neq j\leq n}a_{ij}a_{ji}$
$\displaystyle=M\sum_{i=1}^{n}a_{ii}^{2}+\sigma^{4}\Big{(}({\rm
tr}(A))^{2}-\sum_{i=1}^{n}a_{ii}^{2}\Big{)}+2\sigma^{4}\Big{(}{\rm
tr}(A^{2})-\sum_{i=1}^{n}a_{ii}^{2}\Big{)}$
$\displaystyle=(M-3\sigma^{4})\sum_{i=1}^{n}a_{ii}^{2}+\sigma^{4}\Big{(}({\rm
tr}(A))^{2}+2{\rm tr}(A^{2})\Big{)}$

Together with, $\mathbb{E}(V^{\prime}AV)=\sigma^{2}{\rm tr}(A)$, we know that,

$\displaystyle{\rm
var}\Big{[}(V^{\prime}AV)^{2}\Big{]}\leq(M-3\sigma^{4})\sum_{i=1}^{n}a_{ii}^{2}+2\sigma^{4}{\rm
tr}(A^{2})$

For general $V$ with $\mu\neq 0$, we denote $V=V_{0}+\mu$, then,

$\displaystyle{\rm var}(V^{\prime}AV)$ $\displaystyle={\rm
var}(V_{0}^{\prime}AV_{0}+2\mu^{\prime}AV_{0})$ $\displaystyle\leq 2{\rm
var}(V_{0}^{\prime}AV_{0})+8{\rm var}(\mu^{\prime}AV_{0}).$

This completes the proof. ∎

Now we turn to the proof of Proposition 2.

###### Proof.

Let $\delta^{2}=\|\mu-P_{1}\mu\|^{2}$ and define
$A_{n}=\|P_{X}Y-P_{1}Y\|^{2}/(r-1),B_{n}=\|Y-P_{X}Y\|^{2}/(n-r)$ and their
expected values $a_{n}=\sigma^{2}+{\delta^{2}}/({r-1}),b_{n}=\sigma^{2}$,
respectively. Then we can express $\gamma_{\rm opt}=1-{1}/{f_{n}}$ with its
sample counterpart $1-{1}/{F_{n}}$ where
$F_{n}={A_{n}}/{B_{n}},f_{n}={a_{n}}/{b_{n}}$. At first glance, considering
that

$f_{n}=\frac{a_{n}}{b_{n}}=1+\frac{\delta^{2}}{\sigma^{2}({\rm rank}(X)-1)},$

it is obvious that $f_{n}\geq 1$. By the assumption, we have the error
distribution to be $\mathbb{E}(e_{i}^{4})\leq M$, for $M>0,i=1,\dots,n$.

For the first step, knowing that $\mathbb{E}(B_{n})=b_{n}$,

$\displaystyle\mathbb{E}(B_{n}-b_{n})$ $\displaystyle=0$ $\displaystyle{\rm
var}(B_{n}-b_{n})$ $\displaystyle={\rm var}(B_{n})$
$\displaystyle=\frac{1}{(n-r)^{2}}{\rm var}\|(I-P_{X})Y\|^{2}$
$\displaystyle=\frac{1}{(n-r)^{2}}{\rm
var}\left\\{Y^{\prime}(I-P_{X})Y\right\\}$
$\displaystyle\leq\frac{2(M-\sigma^{4})}{n-r},$

which converges to 0 because $n\rightarrow\infty$ and
$r/n\rightarrow\tau\in[0,1)$. The last inequality is due to Lemma 1, with the
fact that $I-P_{X}$ is an idempotent matrix with rank $n-r$.

So $B_{n}-b_{n}\rightarrow_{q.m}0$, which implies that
$B_{n}/b_{n}\rightarrow_{q.m.}1$ because $b_{n}=\sigma^{2}$ is constant.

Next, knowing that $\mathbb{E}(A_{n})=a_{n}$, similarly,

$\displaystyle\mathbb{E}\Big{(}\frac{A_{n}}{a_{n}}-1\Big{)}$ $\displaystyle=0$
$\displaystyle{\rm var}\Big{(}\frac{A_{n}}{a_{n}}-1\Big{)}$
$\displaystyle=\frac{1}{a_{n}^{2}}{\rm var}(A_{n})$
$\displaystyle=\frac{1}{a_{n}^{2}(r-1)^{2}}{\rm
var}\left\\{\|(P_{X}-P_{1})Y\|^{2}\right\\}$
$\displaystyle=\frac{1}{a_{n}^{2}(r-1)^{2}}{\rm
var}\left\\{Y^{\prime}(P_{X}-P_{1})Y\right\\}$
$\displaystyle\leq\frac{2(M-\sigma^{4})(r-1)+8\delta^{2}\sigma^{2}}{(\sigma^{2}(r-1)+\delta^{2})^{2}}$
(28)
$\displaystyle=\frac{2(M-\sigma^{4})(r-1)+8\delta^{2}\sigma^{2}}{\sigma^{4}(r-1)^{2}+\delta^{4}+2\sigma^{2}\delta^{2}(r-1)},$
(29)

where we used Lemma 1 for (28) with the fact that $P_{X}-P_{1}$ is an
idempotent matrix with rank $r-1$. If $r\rightarrow\infty$ holds, since (29)
is bounded above by $4(r-1)^{-1}+{2(M-\sigma^{4})}({\sigma^{4}(r-1)})^{-1}$,
it converges to $0$.

On the other hand, under the second condition, we have the same convergence
result. This is because we can have an upper bound of (29) as
$({M-\sigma^{4}})({2\sigma^{2}\delta^{2}})^{-1}+{8\sigma^{2}}({\delta^{2}})^{-1}$,
it converges to $0$ with $\delta^{2}\rightarrow\infty$.

Combining two convergence results,
$B_{n}/b_{n}\rightarrow_{q.m}1,A_{n}/a_{n}\rightarrow_{q.m}1$, it implies
$B_{n}/b_{n}\rightarrow_{p}1,A_{n}/a_{n}\rightarrow_{p}1$, and further yields
$F_{n}/f_{n}=\frac{A_{n}/a_{n}}{B_{n}/b_{n}}\rightarrow_{p}1$ by Slutsky’s
theorem.

Additionally, let us denote that $F^{0}_{n}={\rm max}(1,F_{n})$, then we have
$\hat{\gamma}=1-{1}/{F^{0}_{n}}$. Finally, we get

$\displaystyle\mathbb{P}(|\hat{\gamma}-\gamma_{\rm opt}|\geq\epsilon)$
$\displaystyle=\mathbb{P}\Big{(}\Big{|}1-\frac{1}{F^{0}_{n}}-1+\frac{1}{f_{n}}\Big{|}\geq\epsilon\Big{)}$
$\displaystyle=\mathbb{P}\Big{(}\Big{|}\frac{1}{F^{0}_{n}}-\frac{1}{f_{n}}\Big{|}\geq\epsilon\Big{)}$
$\displaystyle=\mathbb{P}\Big{(}\frac{|F^{0}_{n}-f_{n}|}{F^{0}_{n}f_{n}}\geq\epsilon\Big{)}$
$\displaystyle\leq\mathbb{P}\Big{(}\frac{|F_{n}-f_{n}|}{F_{n}f_{n}}\geq\epsilon,F_{n}\geq
1\Big{)}+\mathbb{P}\Big{(}{|1-f_{n}|}\geq\epsilon,F_{n}<1\Big{)}$
$\displaystyle\leq\mathbb{P}\Big{(}\frac{|F_{n}-f_{n}|}{f_{n}}\geq\epsilon,F_{n}\geq
1\Big{)}+\mathbb{P}\Big{(}f_{n}\geq 1+\epsilon,F_{n}<1\Big{)}$
$\displaystyle\leq\mathbb{P}\Big{(}\Big{|}\frac{F_{n}}{f_{n}}-1\Big{|}\geq\epsilon\Big{)}+\mathbb{P}\Big{(}\Big{|}\frac{F_{n}}{f_{n}}-1\Big{|}\geq\frac{\epsilon}{1+\epsilon}\Big{)}\rightarrow
0.$

∎

###  .3 Proof of Proposition 3

######  Lemma 2.

If a random variable $X$ follows a chi-squared distribution with a degree of
freedom $K\geq 5$, and a non-central parameter $\lambda>0$ which is denoted as
$\chi^{2}(K;\lambda)$, the followings holds.

$\displaystyle\mathbb{E}(X^{-1})\leq{\rm
min}\Big{\\{}\frac{1}{K-2},\frac{1}{2\sqrt{\lambda(K-2)}}\Big{\\}}$ (30)
$\displaystyle\mathbb{E}(X^{-2})\leq{\rm
min}\Big{\\{}\frac{1}{2(K-4)},\frac{(K-4+\lambda)/(2\sqrt{\lambda(K-4)})-1}{2\lambda}\Big{\\}}$
(31)

###### Proof.

Xie (1988) verified the following three results for $K\geq 5$.

$\displaystyle\mathbb{E}^{1}_{K}=\int_{0}^{1}s^{K-3}e^{\frac{\lambda(s^{2}-1)}{2}}ds,$
$\displaystyle\mathbb{E}^{1}_{K}=1/\lambda-((K-4)/\lambda)\mathbb{E}^{1}_{K-2},$
$\displaystyle\mathbb{E}^{n}_{K}-\mathbb{E}^{n}_{K+2}=2n\mathbb{E}^{n+1}_{K+2},$

where $\mathbb{E}^{n}_{K}=\mathbb{E}(1/(\chi^{2}(K;\lambda))^{n})$. Hence

$\displaystyle\mathbb{E}(X^{-2})=\mathbb{E}^{2}_{K}$
$\displaystyle=\frac{\mathbb{E}^{1}_{K-2}-\mathbb{E}^{1}_{K}}{2}$
$\displaystyle=\frac{(K-4+\lambda)\mathbb{E}^{1}_{K-2}-1}{2\lambda}$

holds. Furthermore,

$\displaystyle\mathbb{E}(X^{-1})=\mathbb{E}^{1}_{K}$
$\displaystyle=\int_{0}^{1}s^{K-3}e^{\frac{\lambda(s^{2}-1)}{2}}ds$
$\displaystyle\leq\int_{0}^{1}s^{K-3}ds=\frac{1}{K-2}.$

On the other hand, from the same equation,

$\displaystyle\int_{0}^{1}s^{K-3}e^{\frac{\lambda(s^{2}-1)}{2}}ds$
$\displaystyle\leq
e^{-\lambda/2}\Big{(}\int_{0}^{1}s^{2K-5}ds\Big{)}^{1/2}\Big{(}\int_{0}^{1}se^{\lambda
s^{2}}ds\Big{)}^{1/2}$
$\displaystyle=e^{-\lambda/2}\frac{1}{\sqrt{2K-4}}\frac{\sqrt{e^{\lambda}-1}}{\sqrt{2\lambda}}$
$\displaystyle\leq\frac{1}{2\sqrt{\lambda(K-2)}}.$

This completes the proof. ∎

With Lemma 2, we prove the statement of Proposition 3.

###### Proof.

From the definition of $\hat{\gamma}$ (4), we obtain

$\displaystyle|\hat{\gamma}-\gamma_{\rm opt}|$
$\displaystyle=\Big{|}I(F_{n}\geq 1)(1-1/F_{n})-(1-1/f_{n})\Big{|}$
$\displaystyle\leq\Big{|}\frac{1}{f_{n}}-\frac{1}{F_{n}}\Big{|},$ (32)

where $\gamma_{\rm opt}=1-1/f_{n}$,
$f_{n}=\frac{\sigma^{2}(r-1)+\delta^{2}}{\sigma^{2}(r-1)}=1+\frac{\delta^{2}}{\sigma^{2}(r-1)}\in[1,\infty)$.

Let us follow the same notations in Proposition 2 as
$F_{n}={A_{n}}/{B_{n}},f_{n}={a_{n}}/{b_{n}}$, where
$A_{n}=\|P_{X}Y-P_{1}Y\|^{2}/(r-1),B_{n}=\|Y-P_{X}Y\|^{2}/(n-r),a_{n}=\sigma^{2}+{\delta^{2}}/({r-1}),b_{n}=\sigma^{2}$.
Since we have normal errors, $A_{n},B_{n}$ are independent. In addition,
$(r-1)A_{n}/\sigma^{2}\sim\chi^{2}(r-1;\delta^{2}/\sigma^{2})$,
$(n-r)B_{n}/\sigma^{2}\sim\chi^{2}(n-r)$, and we denote
$A^{\prime}_{n}:=A_{n}/\sigma^{2},B^{\prime}_{n}:=B_{n}/\sigma^{2}$. First, we
consider the first case where $\delta^{2}=O(r)$. Since $F_{n}$ follows non-
central $F$ distribution,

$\displaystyle\mathbb{E}(F_{n})$
$\displaystyle=\frac{(n-r)(r-1+\delta^{2}/\sigma^{2})}{(r-1)(n-r-2)}$
$\displaystyle=\Big{(}\frac{2}{n-r-2}+1\Big{)}f_{n},$

when $n-r>2$. And, with $n-r>4$,

$\displaystyle{\rm
var}(F_{n})=2\frac{(n-r-2)(n-r)^{2}}{(n-r-4)(n-r-2)^{2}}\Big{[}\frac{(r-1+\delta^{2}/\sigma^{2})^{2}}{(n-r-2)(r-1)^{2}}+\frac{r-1+2\delta^{2}/\sigma^{2}}{(r-1)^{2}}\Big{]}.$
(33)

Then,

$\displaystyle\mathbb{E}(r(F_{n}-f_{n})^{2})$
$\displaystyle=2r\frac{(n-r-2)(n-r)^{2}}{(n-r-4)(n-r-2)^{2}}\left[\frac{(r-1+\delta^{2}/\sigma^{2})^{2}}{(n-r-2)(r-1)^{2}}\right.$
$\displaystyle\qquad+\left.\frac{r-1+2\delta^{2}/\sigma^{2}}{(r-1)^{2}}\right]+\frac{4r}{(n-r-2)^{2}}f_{n}^{2},$

which asymptotically bounded by the assumptions. Furthermore, since we know
that for arbitrary small $\epsilon>0$,

$\displaystyle\hat{\gamma}-\gamma_{\rm opt}=-\gamma_{\rm opt}I(F_{n}\leq
1-\epsilon)+I(F_{n}\geq 1-\epsilon)(\hat{\gamma}-\gamma_{\rm opt})$ (34)

holds, and the second term is bounded by (32). Indeed, on the event $F_{n}\geq
1-\epsilon$, by Cauchy-Schwarz inequality and the fact that $f_{n}\rightarrow
f_{\infty}\in[1,\infty)$,

$\displaystyle\Big{(}\mathbb{E}|1/f_{n}-1/F_{n}|\Big{)}^{2}\leq\frac{1}{f_{n}^{2}}\mathbb{E}|F_{n}-f_{n}|^{2}\mathbb{E}\Big{(}\frac{1}{F_{n}^{2}}\Big{)}$
$\displaystyle\leq\frac{1}{(1-\epsilon)^{2}f_{n}^{2}}\mathbb{E}|F_{n}-f_{n}|^{2}\mathbb{E}\Big{(}\frac{1}{F_{n}^{2}}\Big{)}$
$\displaystyle\asymp\mathbb{E}|F_{n}-f_{n}|^{2}.$

On the other hand, we show that the first term in (34) is negligible. Indeed,
$P(F_{n}\leq 1-\epsilon)$ is at its largest when $f_{n}\rightarrow 1$, since
$F_{n}/f_{n}\rightarrow 1$. Hence, it suffices to show that the tail bound of
$P(F_{n}\leq 1-\epsilon)$ in the case of $f_{n}\rightarrow 1$ is negligible.
Indeed,

$\displaystyle P(F_{n}\leq 1-\epsilon)$
$\displaystyle=P\Big{(}\frac{A^{\prime}_{n}}{B^{\prime}_{n}}\leq
1-\epsilon,A^{\prime}_{n}>1-\frac{\epsilon}{2}\Big{)}+P\Big{(}\frac{A^{\prime}_{n}}{B^{\prime}_{n}}\leq
1-\epsilon,A^{\prime}_{n}\leq 1-\frac{\epsilon}{2}\Big{)}$
$\displaystyle=P\Big{(}B^{\prime}_{n}\geq\frac{1-\frac{\epsilon}{2}}{1-\epsilon}\Big{)}+P\Big{(}A^{\prime}_{n}\leq
1-\frac{\epsilon}{2}\Big{)}.$ (35)

Lemma 1 from Laurent and Massart (2000) with $x=nt/10$ yields $P(Z\geq
2nt)\leq{\rm exp}({-nt/10})$ for all $t\geq 1$, where a random variable $Z$
follows $\chi^{2}(n)$. This further implies that the first term in (35) is
controlled by

$\displaystyle{\rm exp}({-(n-r)(1-\frac{\epsilon}{2})/(20(1-\epsilon))}),$

which has exponential decay in $n$. Furthermore, Theorem 7 of Zhang and Zhou
(2020) yields that the second term in (35) is bounded above by

$\displaystyle{\rm
exp}({-\epsilon^{2}(r-1)^{2}/(16(r-1+\delta^{2}/\sigma^{2})})).$

Based on the case that $\delta^{2}=o(r)$, the upper bound has exponential
decay in $r$. Consequently, $\mathbb{E}|1/f_{n}-1/F_{n}|=O(r^{-1/2})$ for
$\delta^{2}=O(r)$.

However, when $r=o(\delta^{2})$, $f_{n}\rightarrow\infty$, hence the same
reasoning cannot be applied. We use Lemma 2 instead. Indeed,

$\displaystyle\mathbb{E}|1/f_{n}-1/F_{n}|$
$\displaystyle=\mathbb{E}\Big{|}\frac{b_{n}A_{n}-B_{n}a_{n}}{A_{n}a_{n}}\Big{|}$
$\displaystyle\leq\mathbb{E}\frac{b_{n}}{A_{n}}\Big{|}\frac{A_{n}}{a_{n}}-1\Big{|}+\mathbb{E}\frac{1}{A_{n}}|B_{n}-b_{n}|$
$\displaystyle\leq\sigma^{2}\Big{[}\mathbb{E}\Big{(}\frac{A_{n}}{a_{n}}-1\Big{)}^{2}\mathbb{E}\Big{(}\frac{1}{A_{n}}\Big{)}^{2}\Big{]}^{1/2}+\Big{[}\mathbb{E}(B_{n}-b_{n})^{2}\Big{]}^{1/2}\mathbb{E}\Big{(}\frac{1}{A_{n}}\Big{)},$

by Cauchy-Schwarz inequality and independence of $A_{n},B_{n}$. Since

$\displaystyle\mathbb{E}\Big{(}\frac{1}{A_{n}}\Big{)}=\frac{(r-1)}{\sigma^{2}}\mathbb{E}(1/\chi^{2}(r-1;\delta^{2}))$
$\displaystyle\leq\frac{(r-1)}{\sigma^{2}}{\rm
min}\Big{\\{}\frac{1}{r-3},\frac{1}{2\sqrt{\delta^{2}(r-3)}}\Big{\\}}$
$\displaystyle:=c_{1}(r,\delta^{2})/\sigma^{2},$

and,

$\displaystyle\mathbb{E}\Big{(}\frac{1}{A_{n}}\Big{)}^{2}$
$\displaystyle=\frac{(r-1)^{2}}{\sigma^{4}}\mathbb{E}(1/\chi^{2}(r-1;\delta^{2}))^{2}$
$\displaystyle\leq\frac{(r-1)^{2}}{\sigma^{4}}{\rm
min}\Big{\\{}\frac{1}{2(r-5)},\frac{(r-5+\delta^{2})/(2\sqrt{\delta^{2}(r-5)})-1}{2\delta^{2}}\Big{\\}}$
$\displaystyle:=c_{2}(r,\delta^{2})/\sigma^{4}$

from Lemma 2, and combining results from the Proposition 2,

$\displaystyle\mathbb{E}|1/f_{n}-1/F_{n}|$
$\displaystyle\leq\frac{\sqrt{4\sigma^{4}(r-1)+8\delta^{2}\sigma^{2}}}{\sigma^{2}(r-1)+\delta^{2}}\sqrt{c_{2}(r,\delta^{2})}+\frac{\sqrt{2}}{\sqrt{n-r}}c_{1}(r,\delta^{2})$
$\displaystyle\leq\frac{2\sqrt{2}\sigma\sqrt{c_{2}(r,\delta^{2})}}{\sqrt{\sigma^{2}(r-1)+\delta^{2}}}+\frac{\sqrt{2}}{\sqrt{n-r}}c_{1}(r,\delta^{2})$
(36)

If $\delta^{2}/r=o(1)$ holds with growing $r$, the first term in (36) is of
order $O_{P}((\delta^{2}/r)^{-3/4})$, and the second term is of
$O_{P}(n^{-1/2}(\delta^{2}/r)^{-1/2})$. Similarly we can have analogous rates
for the case where $r$ is not divergent, and this completes the proof. ∎

###  .4 Proof of Proposition 4

###### Proof.

We denote $C_{n}=\|Y-P_{1}Y\|^{2},c_{n}=\delta^{2}+\sigma^{2}(n-1)$. Referring
to Lemma 1 with an rank $n-1$ idempotent matrix $A=I-P_{1}$, we have

$\displaystyle\frac{{\rm
var}C_{n}}{c_{n}^{2}}\leq\frac{2(M-\sigma^{4})(n-1)+8\sigma^{2}\delta^{2}}{(\delta^{2}+\sigma^{2}(n-1))^{2}}\rightarrow
0.$

Thus, with $\mathbb{E}C_{n}=c_{n}$, we obtain $C_{n}/c_{n}\rightarrow_{P}1$.
Now,

$\displaystyle\frac{Y^{\prime}(I-P_{1}-\frac{n-1}{nC}(I-K))Y}{Y^{\prime}(I-P_{1})Y}-\frac{\delta^{2}}{\delta^{2}+\sigma^{2}(n-1)}$
$\displaystyle=\frac{C_{n}-(n-1)\check{\sigma}^{2}_{c}}{C_{n}}-\frac{c_{n}-\sigma^{2}(n-1)}{c_{n}}$
$\displaystyle=\frac{(n-1)C_{n}\sigma^{2}-(n-1)c_{n}\check{\sigma}^{2}_{c}}{C_{n}c_{n}}$
$\displaystyle=\frac{(n-1)\sigma^{2}(C_{n}-c_{n})-(n-1)c_{n}(\sigma^{2}-\check{\sigma}^{2}_{c})}{C_{n}c_{n}}.$
(37)

Since $C_{n}/c_{n}\rightarrow_{P}1$, ${(n-1)\sigma^{2}}/{c_{n}}\leq 1$, the
first term in the (37) converges to 0 in probability. Next, we prove the
consistency of $\check{\sigma}^{2}_{c}$. Indeed, since $K$ and $P_{1}$ are
simultaneously diagonalizable and $KP_{1}=P_{1}K=P_{1}$,

$\displaystyle\mathbb{E}\frac{1}{nC}Y^{\prime}(I-K)Y$
$\displaystyle=\frac{1}{nC}(\sigma^{2}{\rm tr}(I-K)+(\mu-
P_{1}\mu)^{\prime}(I-K)(\mu-P_{1}\mu))$
$\displaystyle=\sigma^{2}+\frac{1}{nC}(\mu-P_{1}\mu)^{\prime}(I-K)(\mu-
P_{1}\mu).$

The eigenvalues of $I-K$ are $\frac{2n\alpha}{d_{j}+2n\alpha}$, for
$j=2,\dots,n$, and thus, $nC={\rm
tr}(I-K)=\sum_{j=2}^{n}\frac{2n\alpha}{d_{j}+2n\alpha}$. Hence,

$\displaystyle\frac{1}{nC}(\mu-P_{1}\mu)^{\prime}(I-K)(\mu-P_{1}\mu)$
$\displaystyle\leq\frac{\frac{1}{d_{n}+2n\alpha}}{\sum_{j=2}^{n}\frac{1}{d_{j}+2n\alpha}}\delta^{2}$
$\displaystyle\leq\frac{\frac{1}{d_{n}+2n\alpha}}{\frac{n-1}{d_{2}+2n\alpha}}\delta^{2}=\frac{d_{2}+2n\alpha}{d_{n}+2n\alpha}\frac{\delta^{2}}{n-1}\rightarrow
0,$

since the assumptions $d_{2}(n\alpha)^{-1}=o_{P}(1)$, and $\delta^{2}=o(n)$,
forces $\lim_{n\rightarrow\infty}\frac{d_{2}+2n\alpha}{d_{n}+2n\alpha}=1$.
This yields $\mathbb{E}\check{\sigma}^{2}_{c}-\sigma^{2}\rightarrow 0$. In
addition, by Lemma 2, letting $\phi_{j}=\frac{1}{d_{j}+2n\alpha}$,

$\displaystyle{\rm var}\Big{(}\frac{1}{nC}Y^{\prime}(I-K)Y\Big{)}$
$\displaystyle\leq\frac{4\sigma^{4}}{n^{2}C^{2}}{\rm
tr}(I-K)^{2}+\frac{8\sigma^{2}}{n^{2}C^{2}}(\mu-
P_{1}\mu)^{\prime}(I-K)^{2}(\mu-P_{1}\mu)$ $\displaystyle\leq
4\sigma^{4}\frac{\phi_{2}^{2}+\dots+\phi_{n}^{2}}{(\phi_{2}+\dots+\phi_{n})^{2}}+8\sigma^{2}\delta^{2}\frac{\phi_{n}^{2}}{(\phi_{2}+\dots+\phi_{n})^{2}}$
$\displaystyle\leq
4\sigma^{4}\frac{\phi_{n}^{2}}{n\phi_{2}^{2}}+8\sigma^{2}\delta^{2}\frac{\phi_{n}^{2}}{n^{2}\phi_{2}^{2}}\rightarrow
0.$

Consequently, $\check{\sigma}^{2}_{c}-\sigma^{2}\rightarrow_{P}0$, and this
further yields that the second term in (37) converges to 0 in probability.
Finally, due to $\gamma_{\rm opt}\geq 0$,
$\Big{|}\frac{Y^{\prime}(I-P_{1}-\frac{n-1}{nC}(I-K))Y}{Y^{\prime}(I-P_{1})Y}-\gamma_{\rm
opt}\Big{|}\geq|\bar{\gamma}_{c}-\gamma_{\rm opt}|$ holds. Thus, this
concludes that $\bar{\gamma}_{c}-\gamma_{\rm opt}\rightarrow_{P}0$. ∎

###  .5 Proof of Proposition 5

###### Proof.

Similarly to the proof of Proposition 4, it suffices to show

$\displaystyle\frac{(n-1)(\check{\sigma}^{2}-\sigma^{2})}{c_{n}}\rightarrow_{P}0.$

Since, analogous to the proof of Proposition 4, we have

$\displaystyle\mathbb{E}n^{-1}Y^{\prime}(I-K)Y$
$\displaystyle=\frac{1}{n}(\sigma^{2}{\rm tr}(I-K)+(\mu-
P_{1}\mu)^{\prime}(I-K)(\mu-P_{1}\mu))$
$\displaystyle\leq\frac{\sigma^{2}}{n}\sum_{i=2}^{n}\frac{2n\alpha}{d_{i}+2n\alpha}+\frac{\delta^{2}}{n}\frac{2n\alpha}{d_{n}+2n\alpha}.$

This yields,

$\displaystyle\mathbb{E}n^{-1}Y^{\prime}(I-K)Y-\sigma^{2}$
$\displaystyle\leq\frac{\sigma^{2}}{n}+\frac{\sigma^{2}}{n}\sum_{i=2}^{n}\frac{d_{i}}{d_{i}+2n\alpha}+\frac{\delta^{2}}{n}\frac{2n\alpha}{d_{n}+2n\alpha}$
$\displaystyle\leq\frac{\sigma^{2}}{n}+\frac{(n-1)\sigma^{2}}{n}\frac{d_{2}}{d_{2}+2n\alpha}+\frac{\delta^{2}}{n}\frac{2n\alpha}{d_{n}+2n\alpha}.$

Hence, up to a constant,

$\displaystyle\frac{n-1}{c_{n}}\Big{(}\mathbb{E}n^{-1}Y^{\prime}(I-K)Y-\sigma^{2}\Big{)}\lesssim\frac{n-1}{\sigma^{2}(n-1)+\delta^{2}}\Big{(}\sigma^{2}\frac{d_{2}}{d_{2}+2n\alpha}+\frac{\delta^{2}}{n}\frac{2n\alpha}{d_{n}+2n\alpha}\Big{)}.$
(38)

Under the first condition of $\delta^{2}=o(n)$ and
$d_{2}(n\alpha)^{-1}=o_{P}(1)$, right hand side of (38) goes to $0$. On the
other hand, under $n=\delta^{2}$ and $n\alpha(d_{n})^{-1}=o_{P}(1)$, (38) also
converges to $0$ as well. In addition, analogously to Proposition 4, recalling
that $c_{n}=(n-1)\sigma^{2}+\delta^{2}$, up to a constant, we have

$\displaystyle\frac{(n-1)^{2}}{c_{n}^{2}}{\rm
var}\Big{(}\frac{1}{n}Y^{\prime}(I-K)Y\Big{)}$ (39)
$\displaystyle\qquad\leq\frac{(n-1)^{2}}{c_{n}^{2}}\Big{(}\frac{4\sigma^{4}}{n^{2}}{\rm
tr}(I-K)^{2}+\frac{8\sigma^{2}}{n^{2}}(\mu-P_{1}\mu)^{\prime}(I-K)^{2}(\mu-
P_{1}\mu)\Big{)}$
$\displaystyle\qquad\lesssim\frac{4n\sigma^{4}}{c_{n}^{2}}\Big{(}\frac{2n\alpha}{d_{n}+2n\alpha}\Big{)}^{2}+\frac{8\sigma^{2}\delta^{2}}{c_{n}^{2}}\Big{(}\frac{2n\alpha}{d_{n}+2n\alpha}\Big{)}^{2}.$
(40)

One can easily check that (40) converges to $0$ if any of two conditions
holds. Consequently, $\check{\sigma}^{2}_{c}-\sigma^{2}\rightarrow_{P}0$, and
this completes the proof. ∎

## Additional plots and tables for numerical experiments

###  .1 Additional plots for simulations

Additional supplement plots are in Figure 7–12.

(a) $\tau=10^{-2}$

(b) $\tau=10^{-1}$

(c) $\tau=1$

(d) $\tau=10^{0.5}$

Figure 7:  Boxplots of the observed same-X loss values from the 50
replications when $n=300$ and $p=250$. See the caption of Figure 1 for more
details.

(a) $\tau=10^{-0.5}$

(b) $\tau=1$

(c) $\tau=10^{0.5}$, magnified

(d) $\tau=10^{0.5}$

Figure 8: Additional displays to Section 6.1 when $(n,p)=(300,75)$. Graphs of
the average same-X loss values over the 50 replications as a function of
$\lambda$ comparing our proposed method (2) to Ridge (24). The label FVS2 is
fitted value shrinkage with $L_{2}$-squared penalty (2), and the horizontal
lines, min(FVS2) and min(Ridge), are the minimal average same-X loss values
obtained from each estimator within the set of $\lambda$.

(a) $\tau=10^{-0.5}$

(b) $\tau=1$

(c) $\tau=10^{0.5}$, magnified

(d) $\tau=10^{0.5}$

Figure 9: Additional displays to Section 6.1 when $(n,p)=(300,150)$. Graphs of
the average same-X loss values over the 50 replications as a function of
$\lambda$ comparing our proposed method (2) to Ridge (24). The label FVS2 is
fitted value shrinkage with $L_{2}$-squared penalty (2), and the horizontal
lines, min(FVS2) and min(Ridge), are the minimal average same-X loss values
obtained from each estimator within the set of $\lambda$.

(a) $\tau=10^{-0.5}$

(b) $\tau=1$

(c) $\tau=10^{0.5}$, magnified

(d) $\tau=10^{0.5}$

Figure 10: Additional displays to Section 6.1 when $(n,p)=(300,250)$. Graphs
of the average same-X loss values over the 50 replications as a function of
$\lambda$ comparing our proposed method (2) to Ridge (24). The label FVS2 is
fitted value shrinkage with $L_{2}$-squared penalty (2), and the horizontal
lines, min(FVS2) and min(Ridge), are the minimal average same-X loss values
obtained from each estimator within the set of $\lambda$.

(a) $(1,10^{-1})$

(b) $(1,10)$

(c) $(2,10^{-1})$

(d) $(2,10)$

(e) $(3,10^{-1})$

(f) $(3,10)$

Figure 11: Additional displays to Figure 5. Boxplots of the observed same-X
loss values from the 50 replications when $(n,p)=(200,300)$.

(a) $(1,10^{-0.5})$

(b) $(1,1)$

(c) $(2,1)$

(d) $(2,10^{0.5})$

(e) $(3,1)$

(f) $(3,10^{0.5})$

Figure 12: Additional displays to Section 6.4 when $(n,p)=(200,300)$. Graphs
of the average same-X loss values labeled with $(\sigma,\tau)$ over the 50
replications as a function of $\lambda$ comparing our proposed method (2) to
Ridge (24). The label FVS2 is fitted value shrinkage with $L_{2}$-squared
penalty (2), and the horizontal lines, min(FVS2) and min(Ridge), are the
minimal average same-X loss values obtained from each estimator within the set
of $\lambda$.

###  .2 Additional tables for Section 6.1

Additional supplement tables are in Table 4–12.

Performance comparison table when $(n,p)=(300,75)$
---
$\tau$ | OLS | Lasso | Ridge | 2n-G | 2n-Rep | 2n-Or | 2n-Es | 2n-Es95
0 | 0.25402 | 0.01053 | 0.00807 | 0.00689 | 0.01076 | 0.00410 | 0.00883 | 0.00706
| (0.00619) | (0.00179) | (0.0013) | (0.00098) | (0.00141) | (0.00071) | (0.00133) | (0.00136)
$10^{-4}$ | 0.25624 | 0.01554 | 0.01039 | 0.00666 | 0.01049 | 0.00349 | 0.00807 | 0.00681
| (0.00625) | (0.00334) | (0.00214) | (0.00107) | (0.00153) | (0.00056) | (0.00145) | (0.00149)
$10^{-2}$ | 0.24498 | 0.00511 | 0.00536 | 0.00583 | 0.00836 | 0.00332 | 0.00676 | 0.00557
| (0.00626) | (0.00128) | (0.00100) | (0.00095) | (0.00120) | (0.00063) | (0.00114) | (0.00111)
$10^{-1}$ | 0.25250 | 0.00923 | 0.01054 | 0.00873 | 0.01286 | 0.00657 | 0.01096 | 0.00990
| (0.00628) | (0.00119) | (0.00150) | (0.00099) | (0.00163) | (0.00068) | (0.00153) | (0.00154)
$10^{-0.5}$ | 0.24418 | 0.03563 | 0.03294 | 0.03053 | 0.03180 | 0.02890 | 0.03162 | 0.03281
| (0.00585) | (0.00157) | (0.00108) | (0.00076) | (0.00100) | (0.00072) | (0.00092) | (0.00087)
$1$ | 0.26053 | 0.16830 | 0.14650 | 0.13044 | 0.12520 | 0.12140 | 0.12749 | 0.12843
| (0.00628) | (0.00452) | (0.00363) | (0.00316) | (0.00261) | (0.00242) | (0.00278) | (0.00320)
$10^{0.5}$ | 0.25099 | 0.24172 | 0.24065 | 0.23429 | 0.24686 | 0.22584 | 0.22714 | 0.22714
| (0.00548) | (0.00526) | (0.00505) | (0.00487) | (0.00513) | (0.00474) | (0.00472) | (0.00472)
$10$ | 0.24962 | 0.24876 | 0.24805 | 0.24849 | 0.26145 | 0.24688 | 0.24709 | 0.24709
| (0.00552) | (0.00554) | (0.00550) | (0.0054) | (0.00567) | (0.00537) | (0.00538) | (0.00538)
$10^{1.5}$ | 0.23895 | 0.23917 | 0.23931 | 0.23895 | 0.24147 | 0.23895 | 0.23897 | 0.23897
| (0.00570) | (0.00570) | (0.00568) | (0.00570) | (0.00584) | (0.00573) | (0.00572) | (0.00572)
$10^{2}$ | 0.25366 | 0.25340 | 0.25336 | 0.25366 | 0.25322 | 0.25346 | 0.25346 | 0.25346
| (0.00638) | (0.00638) | (0.00641) | (0.00638) | (0.00642) | (0.00639) | (0.00639) | (0.00639)

Table 4: The performance comparison table for each simulation setting in
Section 6.1 with $p=75$. The values are the mean squared same-X errors
averaged over 50 independent replications. The numbers in parentheses are
normalized sample standard deviations. The column labels are defined in
Section 6.1. Boldface indicates the best model. Underlined are our main
proposed estimator 2n-Es.

Performance comparison table when $(n,p)=(300,150)$
---
$\tau$ | OLS | Lasso | Ridge | 2n-G | 2n-Rep | 2n-Or | 2n-Es | 2n-Es95
0 | 0.49592 | 0.00643 | 0.00602 | 0.00495 | 0.0189 | 0.00410 | 0.00947 | 0.00592
| (0.00870) | (0.00106) | (0.00081) | (0.00072) | (0.00181) | (0.00071) | (0.00201) | (0.00191)
$10^{-4}$ | 0.50425 | 0.01029 | 0.00612 | 0.00468 | 0.01963 | 0.00349 | 0.00875 | 0.00576
| (0.00778) | (0.00215) | (0.00081) | (0.00069) | (0.00187) | (0.00056) | (0.00175) | (0.00169)
$10^{-2}$ | 0.48596 | 0.00990 | 0.00672 | 0.00406 | 0.01634 | 0.00335 | 0.00589 | 0.00335
| (0.00725) | (0.00195) | (0.00104) | (0.00067) | (0.00161) | (0.00063) | (0.00102) | (0.00063)
$10^{-1}$ | 0.49670 | 0.01597 | 0.01338 | 0.01085 | 0.02410 | 0.00920 | 0.01506 | 0.01191
| (0.00717) | (0.00239) | (0.00137) | (0.00098) | (0.00211) | (0.00069) | (0.00219) | (0.00214)
$10^{-0.5}$ | 0.49954 | 0.05884 | 0.05137 | 0.04916 | 0.05695 | 0.04651 | 0.05255 | 0.05373
| (0.00725) | (0.00304) | (0.00114) | (0.00085) | (0.00162) | (0.00083) | (0.00153) | (0.00139)
$1$ | 0.50286 | 0.44147 | 0.36465 | 0.34512 | 0.28957 | 0.28137 | 0.28761 | 0.28761
| (0.00767) | (0.01307) | (0.00617) | (0.00664) | (0.00341) | (0.00342) | (0.00365) | (0.00365)
$10^{0.5}$ | 0.49542 | 0.47971 | 0.47440 | 0.52221 | 0.52846 | 0.44556 | 0.44695 | 0.44695
| (0.00826) | (0.00773) | (0.00779) | (0.00997) | (0.00898) | (0.00764) | (0.00760) | (0.00760)
$10$ | 0.51361 | 0.51451 | 0.51391 | 0.52422 | 0.58301 | 0.51103 | 0.51110 | 0.51110
| (0.00782) | (0.00784) | (0.00785) | (0.00829) | (0.00875) | (0.00777) | (0.00780) | (0.00780)
$10^{1.5}$ | 0.48673 | 0.48670 | 0.48671 | 0.48762 | 0.49737 | 0.48617 | 0.48626 | 0.48626
| (0.00837) | (0.00852) | (0.00852) | (0.00826) | (0.00886) | (0.00841) | (0.00841) | (0.00841)
$10^{2}$ | 0.49052 | 0.49088 | 0.49100 | 0.49052 | 0.49097 | 0.49041 | 0.49040 | 0.49040
| (0.00858) | (0.00860) | (0.00857) | (0.00858) | (0.00856) | (0.00857) | (0.00857) | (0.00857)

Table 5: The performance comparison table for each simulation setting in
Section 6.1 with $p=150$. The values are the mean squared same-X errors
averaged over 50 independent replications. The numbers in parentheses are
normalized sample standard deviations. The column labels are defined in
Section 6.1. Boldface indicates the best model. Underlined are our main
proposed estimator 2n-Es.

Performance comparison table when $(n,p)=(300,250)$
---
$\tau$ | OLS | Lasso | Ridge | 2n-G | 2n-Rep | 2n-Or | 2n-Es | 2n-Es95
0 | 0.82705 | 0.00904 | 0.00897 | 0.00434 | 0.05467 | 0.00410 | 0.02362 | 0.01692
| (0.0113) | (0.00165) | (0.00141) | (0.00072) | (0.00307) | (0.00071) | (0.00589) | (0.00562)
$10^{-4}$ | 0.84005 | 0.00856 | 0.00847 | 0.00376 | 0.05760 | 0.00349 | 0.02110 | 0.01159
| (0.00935) | (0.00142) | (0.00141) | (0.00057) | (0.00244) | (0.00056) | (0.00521) | (0.00502)
$10^{-2}$ | 0.81920 | 0.00775 | 0.00867 | 0.00356 | 0.05246 | 0.00339 | 0.01662 | 0.00707
| (0.01102) | (0.00131) | (0.00194) | (0.00065) | (0.00305) | (0.00063) | (0.00368) | (0.00266)
$10^{-1}$ | 0.82471 | 0.01577 | 0.01542 | 0.01291 | 0.05834 | 0.01280 | 0.02765 | 0.01750
| (0.00893) | (0.00126) | (0.00097) | (0.00068) | (0.00251) | (0.00069) | (0.00381) | (0.00324)
$10^{-0.5}$ | 0.83966 | 0.08383 | 0.07713 | 0.07461 | 0.10639 | 0.06991 | 0.09452 | 0.08907
| (0.01143) | (0.00270) | (0.00169) | (0.00077) | (0.00293) | (0.00082) | (0.00535) | (0.00509)
$1$ | 0.84641 | 0.73637 | 0.71328 | 0.68847 | 0.41691 | 0.40367 | 0.42248 | 0.44553
| (0.01109) | (0.00679) | (0.00271) | (0.00592) | (0.00397) | (0.00393) | (0.00563) | (0.01446)
$10^{0.5}$ | 0.83136 | 6.20342 | 5.83755 | 3.30717 | 1.24991 | 0.75871 | 0.76393 | 0.76393
| (0.00822) | (0.40666) | (0.31791) | (0.09648) | (0.01594) | (0.00788) | (0.00825) | (0.00825)
$10$ | 0.84872 | 0.84674 | 0.84765 | 2.16844 | 1.46047 | 0.83994 | 0.84025 | 0.84025
| (0.01086) | (0.01077) | (0.01085) | (0.07311) | (0.01548) | (0.01067) | (0.01074) | (0.01074)
$10^{1.5}$ | 0.81350 | 0.81372 | 0.81371 | 1.07344 | 1.02267 | 0.81182 | 0.81190 | 0.81190
| (0.00955) | (0.00954) | (0.00954) | (0.03146) | (0.01159) | (0.00950) | (0.00950) | (0.00950)
$10^{2}$ | 0.82234 | 0.82403 | 0.82401 | 0.84015 | 0.84084 | 0.82251 | 0.82254 | 0.82254
| (0.01043) | (0.01057) | (0.01057) | (0.02032) | (0.01077) | (0.01044) | (0.01044) | (0.01044)

Table 6: The performance comparison table for each simulation setting in
Section 6.1 with $p=250$. The values are the mean squared same-X errors
averaged over 50 independent replications. The numbers in parentheses are
normalized sample standard deviations. The column labels are defined in
Section 6.1. Boldface indicates the best model. Underlined are our main
proposed estimator 2n-Es.

Two sample test table
---
$p$ $\backslash$ $\tau$ | 0 | $10^{-4}$ | $10^{-2}$ | $10^{-1}$ | $10^{-0.5}$
$p=75$ | (-0.00097, 0.00249) | (-0.00429, -0.00034) | (-0.00017, 0.00297) | (-0.00113, 0.00196) | (-0.00234, -0.00031)
$p=150$ | (-0.00060, 0.00749) | (-0.00079. 0.00604) | (-0.00262, 0.00096) | (-0.00193, 0.00528) | (-0.00100, 0.00337)
$p=250$ | (0.00312, 0.02617) | (0.00213, 0.02313) | (0.00011, 0.01577) | (0.00458, 0.01987) | (0.00747, 0.02731)

Table 7: The table includes paired t-test 95% confidence intervals of
$($squared error of our estimator$-$squared error of Ridge$)$ from the
simulation setting in Section 6.1 when $n=300$ and $\tau\leq 10^{-0.5}$. The
confidence intervals are constructed upon paired 50 independent replications.
Boldface indicates a case where FVS outperforms Ridge with statistical
significance.

Two sample test table
---
$p$ $\backslash$ $\tau$ | 1 | $10^{0.5}$ | $10$ | $10^{1.5}$ | $10^{2}$
$p=75$ | (-0.02336, -0.01468) | (-0.01650, -0.01050) | (-0.00258, 0.00067) | (-0.00091, 0.00022) | (-0.00047, 0.00068)
$p=150$ | (-0.08603, -0.06804) | (-0.03273, -0.02217) | (-0.00454, -0.00109) | (-0.00143, -0.00054) | (-0.00091, -0.00029)
$p=250$ | (-0.30338, -0.27822) | (-5.71148, -4.43574) | (-0.00975, -0.00506) | (-0.00277, -0.00086) | (-2.982e-03, 3.625e-05)

Table 8: The table includes paired t-test 95% confidence intervals of
$($squared error of our estimator$-$squared error of Ridge$)$ from the
simulation setting in Section 6.1 when $n=300$ and $\tau\geq 1$. The
confidence intervals are constructed upon paired 50 independent replications.
Boldface indicates a case where FVS outperforms Ridge with statistical
significance.

Estimation accuracy table
---
$p$ $\backslash$ $\tau$ | 0 | $10^{-4}$ | $10^{-2}$ | $10^{-1}$ | $10^{-0.5}$ | 1 | $10^{0.5}$ | $10$ | $10^{1.5}$ | $10^{2}$
$p=75$ | 0.012765 | 0.015367 | 0.012857 | 0.015675 | 0.010584 | 0.006501 | 8.4549e-05 | 8.3909e-07 | 6.86e-09 | 6.3e-11
| (0.004007) | (0.003858) | (0.003706) | (0.004259) | (0.001927) | (0.001748) | (1.45013e-05) | (1.85515e-07) | (1.296e-09) | (1.5e-11)
$p=150$ | 0.011479 | 0.011239 | 0.009588 | 0.011539 | 0.012060 | 0.005970 | 0.000219 | 9.80702e-07 | 9.625e-09 | 2.3e-10
| (0.003712) | (0.003248) | (0.002850) | (0.003683) | (0.002737) | (0.001387) | (4.54872e-05) | (1.93014e-07) | (1.969e-09) | (4.6e-11)
$p=250$ | 0.022469 | 0.015587 | 0.020570 | 0.020587 | 0.018257 | 0.015929 | 0.000351 | 4.32618e-06 | 4.5542e-08 | 5.84e-10
| (0.005912) | (0.004096) | (0.005552) | (0.005919) | (0.004500) | (0.003126) | (6.79379e-05) | (1.11565e-06) | (7.827e-09) | (9.9e-11)

Table 9: The values are the mean squared estimation errors
$|\hat{\gamma}-\gamma_{\rm opt}|^{2}$ averaged over 100 independent
replications from the same simulation setting in Section 6.1. The numbers in
parentheses are normalized sample standard deviations. The sample size is
$n=300$, and $p$ varies in $\\{75,150,250\\}$.

Performance comparison table $(n,p)=(300,75)$
---
$p$ $\backslash$ $\tau$ | 0 | $10^{-4}$ | $10^{-2}$ | $10^{-1}$ | $10^{-0.5}$ | 1 | $10^{0.5}$ | $10$ | $10^{1.5}$ | $10^{2}$
FVS2 | 0.00320 | 0.00293 | 0.00332 | 0.00607 | 0.02776 | 0.11462 | 0.22182 | 0.25013 | 0.25000 | 0.25799
| (31622.78) | (31622.78) | (5623.413) | (128.8250) | (8.70964) | (1.09648) | (0.10965) | (0.01148) | (0.00066) | (2.13796e-05)
Ridge | 0.00323 | 0.00297 | 0.00337 | 0.00610 | 0.02945 | 0.12852 | 0.23381 | 0.25079 | 0.25007 | 0.25799
| (31622.78) | (31622.78) | (31622.78) | (31622.78) | (5370.318) | (245.4709) | (12.58925) | (5.49541) | (0.67608) | (0.10233)

Table 10: The values are the minimal average same-X loss values obtained from
each estimator (see (26) and (27)) within the set of
$\lambda\in\\{10^{-7+0.01k}:j=0,1,\ldots,1150\\}$ from the simulation settings
in Section 6.1 when $(n,p)=(300,75)$. The numbers in parentheses are the
values of $\lambda$ that achieve the minimum for each setting and each
estimator. Boldface indicates the best model for each $\tau$.

Performance comparison table $(n,p)=(300,150)$
---
$p$ $\backslash$ $\tau$ | 0 | $10^{-4}$ | $10^{-2}$ | $10^{-1}$ | $10^{-0.5}$ | 1 | $10^{0.5}$ | $10$ | $10^{1.5}$ | $10^{2}$
FVS2 | 0.00305 | 0.00310 | 0.00470 | 0.00918 | 0.04589 | 0.27579 | 0.44741 | 0.50089 | 0.49639 | 0.50113
| (31622.78) | (31622.78) | (31622.78) | (109.648) | (10.23293) | (0.79433) | (0.11482) | (0.00708) | (0.00083) | (6.60693e-05)
Ridge | 0.00313 | 0.00320 | 0.00479 | 0.00922 | 0.04796 | 0.33086 | 0.47470 | 0.50273 | 0.49678 | 0.50103
| (31622.78) | (31622.78) | (31622.78) | (31622.78) | (5754.399) | (138.0384) | (12.30269) | (1.77828) | (0.18197) | (0.43652)

Table 11: The values are the minimal average same-X loss values obtained from
each estimator (see (26) and (27)) within the set of
$\lambda\in\\{10^{-7+0.01k}:j=0,1,\ldots,1150\\}$ from the simulation settings
in Section 6.1 when $(n,p)=(300,150)$. The numbers in parentheses are the
values of $\lambda$ that achieve the minimum for each setting and each
estimator. Boldface indicates the best model for each $\tau$.

Performance comparison table $(n,p)=(300,250)$
---
$p$ $\backslash$ $\tau$ | 0 | $10^{-4}$ | $10^{-2}$ | $10^{-1}$ | $10^{-0.5}$ | 1 | $10^{0.5}$ | $10$ | $10^{1.5}$ | $10^{2}$
FVS2 | 0.02187 | 0.02105 | 0.02609 | 0.01134 | 0.07150 | 0.41713 | 0.90093 | 0.90040 | 0.92051 | 0.83885
| (31622.78) | (31622.78) | (31622.78) | (323.5937) | (11.48154) | (1.07152) | (0.00047) | (1.07152e-05) | (3.38844e-07) | (2.63027e-05)
Ridge | 0.02235 | 0.02123 | 0.02626 | 0.01138 | 0.15322 | 0.62548 | 0.90505 | 0.90242 | 0.84360 | 0.84018
| (10000) | (31622.78) | (31622.78) | (31622.78) | (7762.471) | (245.4709) | (0.16982) | (0.00115) | (1.02329) | (1e-07)

Table 12: The values are the minimal average same-X loss values obtained from
each estimator (see (26) and (27)) within the set of
$\lambda\in\\{10^{-7+0.02k}:j=0,1,\ldots,1150\\}$ from the simulation settings
in Section 6.1 when $(n,p)=(300,250)$. The numbers in parentheses are the
values of $\lambda$ that achieve the minimum for each setting and each
estimator. Boldface indicates the best model for each $\tau$.

###  .3 Additional table for Section 6.2

Performance comparison table
---
$(\tau_{c},\tau_{f})$,Coding | OLS | 2n-Or | 2n-Es | 2n-Es95 | 2n-G | 2n-Rep | Ridge | LASSO | GL | MGL | SGL
($10^{-0.5}$,$10^{-0.5}$) | 0.50186 | 0.11723 | 0.14239 | 0.15445 | 0.12866 | 0.12664 | 0.12925 | 0.15168 | 0.16374 | 0.15945 | 0.15307
Coding-1 | (0.01530) | (0.00261) | (0.00467) | (0.00354) | (0.00309) | (0.00319) | (0.00540) | (0.00645) | (0.00812) | (0.00579) | (0.00547)
($10^{-0.5}$,$10^{-0.5}$) | 0.50186 | 0.11723 | 0.14239 | 0.15445 | 0.12866 | 0.12663 | 0.13380 | 0.15266 | 0.16712 | 0.15874 | 0.15489
Coding-2 | (0.01530) | (0.00261) | (0.00467) | (0.00354) | (0.00309) | (0.00317) | (0.00499) | (0.00708) | (0.00862) | (0.00693) | (0.00658)
(1,$10^{-0.5}$) | 0.50927 | 0.17630 | 0.19666 | 0.22610 | 0.20258 | 0.18252 | 0.22405 | 0.24681 | 0.25424 | 0.24891 | 0.24874
Coding-1 | (0.01303) | (0.00430) | (0.00513) | (0.00656) | (0.00517) | (0.00443) | (0.00589) | (0.00651) | (0.00562) | (0.00774) | (0.00697)
(1,$10^{-0.5}$) | 0.50927 | 0.17630 | 0.19666 | 0.22610 | 0.20258 | 0.18265 | 0.22813 | 0.24960 | 0.25910 | 0.24812 | 0.24918
Coding-2 | (0.01303) | (0.00430) | (0.00513) | (0.00656) | (0.00517) | (0.00443) | (0.00485) | (0.00722) | (0.00702) | (0.00725) | (0.00646)
($10^{-0.5}$,1) | 0.49276 | 0.22100 | 0.24620 | 0.27809 | 0.27201 | 0.23677 | 0.21132 | 0.22176 | 0.34283 | 0.27566 | 0.24439
Coding-1 | (0.01510) | (0.00419) | (0.00735) | (0.01182) | (0.00739) | (0.00554) | (0.00819) | (0.01228) | (0.01101) | (0.01400) | (0.01359)
($10^{-0.5}$,1) | 0.49276 | 0.22100 | 0.24620 | 0.27809 | 0.27201 | 0.23757 | 0.28954 | 0.29632 | 0.38464 | 0.34066 | 0.34019
Coding-2 | (0.01510) | (0.00419) | (0.00735) | (0.01182) | (0.00739) | (0.00555) | (0.00869) | (0.01005) | (0.00724) | (0.01048) | (0.01039)
(1,1) | 0.50877 | 0.38063 | 0.39377 | 0.39377 | 0.49954 | 0.44260 | 0.34082 | 0.37666 | 0.40414 | 0.39063 | 0.38946
Coding-1 | (0.01491) | (0.01185) | (0.01306) | (0.01306) | (0.01955) | (0.01393) | (0.01153) | (0.01203) | (0.01342) | (0.01355) | (0.01303)
(1,1) | 0.50877 | 0.38063 | 0.39377 | 0.39377 | 0.49954 | 0.46147 | 0.46169 | 0.47966 | 0.47640 | 0.44753 | 0.44525
Coding-2 | (0.01491) | (0.01185) | (0.01306) | (0.01306) | (0.01955) | (0.01426) | (0.02263) | (0.01306) | (0.03130) | (0.02385) | (0.01600)
($10^{0.5}$,1) | 0.49117 | 0.43134 | 0.43714 | 0.43714 | 0.54382 | 0.51857 | 0.42296 | 0.42204 | 0.41071 | 0.41779 | 0.41682
Coding-1 | (0.01107) | (0.00918) | (0.00937) | (0.00937) | (0.01375) | (0.01127) | (0.00959) | (0.01019) | (0.01088) | (0.00985) | (0.00939)
($10^{0.5}$,1) | 0.49117 | 0.43134 | 0.43714 | 0.43714 | 0.54382 | 0.52616 | 0.46032 | 0.47328 | 0.47515 | 0.46817 | 0.46450
Coding-2 | (0.01107) | (0.00918) | (0.00937) | (0.00937) | (0.01375) | (0.01135) | (0.01040) | (0.01158) | (0.01182) | (0.01196) | (0.01110)
(1,$10^{0.5}$) | 0.49584 | 0.46276 | 0.46526 | 0.46526 | 0.54069 | 0.62294 | 0.40551 | 0.32561 | 0.40242 | 0.37382 | 0.36664
Coding-1 | (0.01566) | (0.01462) | (0.01486) | (0.01486) | (0.01945) | (0.01823) | (0.01282) | (0.01191) | (0.01405) | (0.01260) | (0.01275)
(1,$10^{0.5}$) | 0.49584 | 0.46276 | 0.46526 | 0.46526 | 0.54069 | 0.65369 | 0.49049 | 0.47572 | 0.47647 | 0.46569 | 0.47035
Coding-2 | (0.01566) | (0.01462) | (0.01486) | (0.01486) | (0.01945) | (0.01840) | (0.01506) | (0.01511) | (0.01465) | (0.01399) | (0.01488)
($10^{0.5}$,$10^{0.5}$) | 0.48446 | 0.47790 | 0.47878 | 0.47878 | 0.51832 | 0.56642 | 0.47364 | 0.44926 | 0.46486 | 0.46730 | 0.47053
Coding-1 | (0.01323) | (0.01278) | (0.01282) | (0.01282) | (0.01632) | (0.01528) | (0.01344) | (0.01239) | (0.01290) | (0.01303) | (0.01308)
($10^{0.5}$,$10^{0.5}$) | 0.48446 | 0.47790 | 0.47878 | 0.47878 | 0.51832 | 0.62073 | 0.47599 | 0.47412 | 0.47585 | 0.47288 | 0.47398
Coding-2 | (0.01323) | (0.01278) | (0.01282) | (0.01282) | (0.01632) | (0.01692) | (0.01341) | (0.01325) | (0.01318) | (0.01302) | (0.01295)

Table 13: The performance comparison table for each simulation setting in
Section 6.2 characterized by a tuple of $(\tau_{c},\tau_{f})$ and a Coding
number. The values are the mean squared same-X errors averaged over 50
independent replications. The numbers in parentheses are normalized sample
standard deviations. The column labels are defined in Section 6.2. Boldface
indicates the best model. Underlined is our main proposed estimator.

###  .4 Additional table for Section 6.3

Performance comparison table
---
$(\psi,s)$,Coding | OLS | 2n-Or | 2n-Es | 2n-Es95 | 2n-G | 2n-Rep | Ridge | LASSO
(0,0.025) | 0.50839 | 0.39358 | 0.39762 | 0.39762 | 0.47646 | 0.40001 | 0.31489 | 0.06404
Coding-1 | (0.00846) | (0.00662) | (0.00649) | (0.00649) | (0.01003) | (0.00664) | (0.00569) | (0.00375)
(0,0.025) | 0.50839 | 0.39358 | 0.39762 | 0.39762 | 0.47646 | 0.42485 | 0.60588 | 0.68829
Coding-2 | (0.00846) | (0.00662) | (0.00649) | (0.00649) | (0.01003) | (0.00736) | (0.01615) | (0.03787)
(0,0.05) | 0.49586 | 0.46544 | 0.46678 | 0.46678 | 0.52128 | 0.47463 | 0.40975 | 0.11051
Coding-1 | (0.00794) | (0.00769) | (0.00778) | (0.00778) | (0.00896) | (0.00809) | (0.00650) | (0.00441)
(0,0.05) | 0.49586 | 0.46544 | 0.46678 | 0.46678 | 0.52128 | 0.47599 | 0.68096 | 0.69950
Coding-2 | (0.00794) | (0.00769) | (0.00778) | (0.00778) | (0.00896) | (0.00812) | (0.03665) | (0.04346)
(0,0.1) | 0.49187 | 0.47031 | 0.47119 | 0.47119 | 0.51185 | 0.48061 | 0.44077 | 0.17147
Coding-1 | (0.00753) | (0.00713) | (0.00713) | (0.00713) | (0.00900) | (0.00733) | (0.00664) | (0.00481)
(0,0.1) | 0.49187 | 0.47031 | 0.47119 | 0.47119 | 0.51185 | 0.48354 | 0.58174 | 0.63748
Coding-2 | (0.00753) | (0.00713) | (0.00713) | (0.00713) | (0.00900) | (0.00740) | (0.04227) | (0.04410)
(0,0.2) | 0.51313 | 0.50418 | 0.50468 | 0.50468 | 0.52139 | 0.50593 | 0.48666 | 0.26380
Coding-1 | (0.00779) | (0.00776) | (0.00778) | (0.00778) | (0.00873) | (0.00791) | (0.00739) | (0.00769)
(0,0.2) | 0.51313 | 0.50418 | 0.50468 | 0.50468 | 0.52139 | 0.50602 | 1.08212 | 1.76470
Coding-2 | (0.00779) | (0.00776) | (0.00778) | (0.00778) | (0.00873) | (0.00791) | (0.01940) | (0.05839)
(0,0.3) | 0.49183 | 0.48747 | 0.48770 | 0.48770 | 0.50323 | 0.49037 | 0.46678 | 0.27960
Coding-1 | (0.00657) | (0.00667) | (0.00667) | (0.00667) | (0.00713) | (0.00681) | (0.00653) | (0.00569)
(0,0.3) | 0.49183 | 0.48747 | 0.48770 | 0.48770 | 0.50323 | 0.49046 | 4.19051 | 3.41197
Coding-2 | (0.00657) | (0.00667) | (0.00667) | (0.00667) | (0.00713) | (0.00681) | (0.05049) | (0.08670)
(1,0.025) | 0.50711 | 0.33304 | 0.33710 | 0.33710 | 0.39516 | 0.33499 | 0.25947 | 0.07922
Coding-1 | (0.00784) | (0.00516) | (0.00522) | (0.00522) | (0.00717) | (0.00519) | (0.00400) | (0.00418)
(1,0.025) | 0.50711 | 0.33304 | 0.33710 | 0.33710 | 0.39516 | 0.33688 | 0.56775 | 0.62388
Coding-2 | (0.00784) | (0.00516) | (0.00522) | (0.00522) | (0.00717) | (0.00524) | (0.01124) | (0.01830)
(1,0.05) | 0.48570 | 0.31444 | 0.31995 | 0.31995 | 0.38558 | 0.31870 | 0.25929 | 0.09830
Coding-1 | (0.00845) | (0.00497) | (0.00527) | (0.00527) | (0.00736) | (0.00498) | (0.00591) | (0.00593)
(1,0.05) | 0.48570 | 0.31444 | 0.31995 | 0.31995 | 0.38558 | 0.32063 | 0.46405 | 0.46705
Coding-2 | (0.00845) | (0.00497) | (0.00527) | (0.00527) | (0.00736) | (0.00505) | (0.00982) | (0.01272)
(1,0.1) | 0.51035 | 0.45242 | 0.45569 | 0.45569 | 0.53442 | 0.45991 | 0.38374 | 0.18685
Coding-1 | (0.00893) | (0.00731) | (0.00729) | (0.00729) | (0.01030) | (0.00723) | (0.00619) | (0.00620)
(1,0.1) | 0.51035 | 0.45242 | 0.45569 | 0.45569 | 0.53442 | 0.46254 | 0.87261 | 1.01879
Coding-2 | (0.00893) | (0.00731) | (0.00729) | (0.00729) | (0.01030) | (0.00727) | (0.02527) | (0.05128)
(1,0.2) | 0.49424 | 0.45459 | 0.45662 | 0.45662 | 0.50326 | 0.45741 | 0.38720 | 0.23477
Coding-1 | (0.00814) | (0.00793) | (0.00791) | (0.00791) | (0.00995) | (0.00816) | (0.00664) | (0.00653)
(1,0.2) | 0.49424 | 0.45459 | 0.45662 | 0.45662 | 0.50326 | 0.45783 | 0.63929 | 0.81334
Coding-2 | (0.00814) | (0.00793) | (0.00791) | (0.00791) | (0.00995) | (0.00819) | (0.01402) | (0.03587)
(1,0.3) | 0.50315 | 0.48846 | 0.48878 | 0.48878 | 0.51896 | 0.49062 | 0.41947 | 0.29917
Coding-1 | (0.00816) | (0.00810) | (0.00815) | (0.00815) | (0.00918) | (0.00826) | (0.00729) | (0.00752)
(1,0.3) | 0.50315 | 0.48846 | 0.48878 | 0.48878 | 0.51896 | 0.49078 | 0.67030 | 0.82268
Coding-2 | (0.00816) | (0.00810) | (0.00815) | (0.00815) | (0.00918) | (0.00827) | (0.01589) | (0.04973)

Table 14: The performance comparison table for each simulation setting in
Section 6.3 characterized by a tuple of $(\psi,s)$ and a Coding number. The
values are the mean squared same-X errors averaged over 50 independent
replications. The numbers in parentheses are normalized sample standard
deviations. The column labels are defined in Section 6.1. Boldface indicates
the best model. Underlined is our main proposed estimator.

###  .5 Additional tables for Section 6.4

Additional supplement tables are in Table 15–18.

Performance comparison table when $(n,p,\sigma)=(200,300,1)$
---
$\tau$ | Ridge | Lasso | 2n-G | 2n-Or | 2n-Rep1 | 2n-Rep2 | 2n-Rep3 | 2n-Rep4
0 | 0.01786 | 0.02242 | 0.01469 | 0.00593 | 0.28577 | 0.00601 | 0.00616 | 0.00595
| (0.00339) | (0.00542) | (0.00300) | (0.00129) | (0.00780) | (0.00128) | (0.00129) | (0.00129)
$10^{-4}$ | 0.01432 | 0.01468 | 0.01450 | 0.00689 | 0.26896 | 0.00696 | 0.00711 | 0.00692
| (0.00266) | (0.00290) | (0.00275) | (0.00165) | (0.00749) | (0.00164) | (0.00165) | (0.00165)
$10^{-2}$ | 0.02166 | 0.02623 | 0.01219 | 0.00513 | 0.27957 | 0.00524 | 0.00536 | 0.00515
| (0.00430) | (0.00763) | (0.00255) | (0.00084) | (0.00757) | (0.00083) | (0.00083) | (0.00084)
$10^{-1}$ | 0.02276 | 0.03052 | 0.01604 | 0.01172 | 0.26729 | 0.01180 | 0.01175 | 0.01174
| (0.00293) | (0.00521) | (0.00136) | (0.00061) | (0.00624) | (0.00060) | (0.00061) | (0.00061)
$10^{-0.5}$ | 0.11405 | 0.12083 | 0.10753 | 0.10100 | 0.31219 | 0.11014 | 0.10811 | 0.11011
| (0.00300) | (0.00350) | (0.00146) | (0.00118) | (0.00684) | (0.00120) | (0.00115) | (0.00116)
$1$ | 0.91282 | 0.96185 | 0.97451 | 0.49907 | 0.53213 | 0.97599 | 0.93679 | 0.97061
| (0.00330) | (0.00612) | (0.00338) | (0.00546) | (0.00679) | (0.00175) | (0.00105) | (0.00060)
$10^{0.5}$ | 8.49899 | 9.75743 | 8.42798 | 0.88561 | 0.89875 | 9.80771 | 8.24792 | 9.73334
| (0.00707) | (0.02054) | (0.10022) | (0.01310) | (0.01264) | (0.01598) | (0.01343) | (0.00107)
$10$ | 73.96427 | 106.5722 | 106.5722 | 0.97842 | 1.01806 | 106.5722 | 42.61824 | 104.4773
| (0.04330) | (0.00082) | (0.00082) | (0.01307) | (0.01357) | (0.00082) | (0.10719) | (0.00446)
$10^{1.5}$ | 361.32690 | 791.61740 | 711.54760 | 1.00270 | 1.00436 | 681.97860 | 46.79847 | 728.59430
| (0.17178) | (0.30844) | (0.79395) | (0.01140) | (0.01138) | (1.35906) | (0.15520) | (0.06574)
$10^{2}$ | 2377.538 | 10718.830 | 10718.830 | 0.99812 | 0.99848 | 10522.250 | 14.32732 | 5754.502
| (0.57873) | (0.00150) | (0.00150) | (0.01184) | (0.01183) | (6.38851) | (0.07117) | (0.92518)

Table 15: The performance comparison table for each simulation setting in
Section 6.4 with $(n,p)=(200,300)$ and $\sigma=1$. The values are the mean
squared same-X errors averaged over 50 independent replications. The numbers
in parentheses are normalized sample standard deviations. The column labels
are defined in Section 6.4. Boldface indicates the best model. Underlined is
the best model within our proposed estimators 2n-Rep1 to 2n-Rep4.

Performance comparison table when $(n,p,\sigma)=(200,300,2)$
---
$\tau$ | Ridge | Lasso | 2n-G | 2n-Or | 2n-Rep1 | 2n-Rep2 | 2n-Rep3 | 2n-Rep4
0 | 0.07935 | 0.05637 | 0.06283 | 0.03277 | 2.34612 | 0.03383 | 0.04006 | 0.03288
| (0.01238) | (0.01358) | (0.01605) | (0.00934) | (0.04869) | (0.00934) | (0.00936) | (0.00934)
$10^{-4}$ | 0.06238 | 0.03690 | 0.04443 | 0.01775 | 2.2782 | 0.01927 | 0.02505 | 0.01786
| (0.01024) | (0.00682) | (0.01024) | (0.00341) | (0.05627) | (0.00347) | (0.00342) | (0.00341)
$10^{-2}$ | 0.07281 | 0.06067 | 0.04072 | 0.02118 | 2.30044 | 0.02258 | 0.02851 | 0.02129
| (0.01027) | (0.01262) | (0.00533) | (0.00358) | (0.04880) | (0.00350) | (0.00351) | (0.00358)
$10^{-1}$ | 0.08896 | 0.08272 | 0.05622 | 0.03965 | 2.21296 | 0.04086 | 0.04557 | 0.03967
| (0.01202) | (0.0142) | (0.00805) | (0.00547) | (0.0513) | (0.00549) | (0.00545) | (0.00547)
$10^{-0.5}$ | 0.14131 | 0.14596 | 0.12302 | 0.11325 | 2.22741 | 0.11519 | 0.11446 | 0.11492
| (0.00513) | (0.01019) | (0.00500) | (0.00316) | (0.04056) | (0.00311) | (0.00323) | (0.00312)
$1$ | 0.86595 | 0.96481 | 0.87980 | 0.76578 | 2.57110 | 0.91357 | 0.85542 | 0.92523
| (0.00976) | (0.02339) | (0.01148) | (0.00663) | (0.04396) | (0.00637) | (0.00545) | (0.00549)
$10^{0.5}$ | 7.65122 | 9.14910 | 9.12699 | 2.73597 | 3.20867 | 8.98522 | 7.32603 | 9.15122
| (0.04529) | (0.07148) | (0.04857) | (0.03793) | (0.04368) | (0.07203) | (0.03633) | (0.00588)
$10$ | 73.88077 | 109.19530 | 89.77449 | 3.88429 | 3.92254 | 88.15042 | 39.3269 | 107.6866
| (0.08510) | (0.29530) | (0.64223) | (0.06352) | (0.06398) | (0.50894) | (0.18335) | (0.00819)
$10^{1.5}$ | 481.4546 | 1002.2620 | 1002.2620 | 3.91530 | 3.91610 | 1002.2620 | 63.5840 | 911.6336
| (0.31973) | (0.00385) | (0.00385) | (0.05169) | (0.05206) | (0.00385) | (0.31341) | (0.16549)
$10^{2}$ | 2359.036 | 10597.59 | 10597.59 | 3.99555 | 3.99651 | 10597.59 | 18.45940 | 5602.045
| (1.05413) | (0.00368) | (0.00368) | (0.06336) | (0.06341) | (0.00368) | (0.16666) | (2.19424)

Table 16: The performance comparison table for each simulation setting in
Section 6.4 with $(n,p)=(200,300)$ and $\sigma=2$. The values are the mean
squared same-X errors averaged over 50 independent replications. The numbers
in parentheses are normalized sample standard deviations. The column labels
are defined in Section 6.4. Boldface indicates the best model. Underlined is
the best model within our proposed estimators 2n-Rep1 to 2n-Rep4.

Performance comparison table when $(n,p,\sigma)=(200,300,3)$
---
$\tau$ | Ridge | Lasso | 2n-G | 2n-Or | 2n-Rep1 | 2n-Rep2 | 2n-Rep3 | 2n-Rep4
0 | 0.22832 | 0.16547 | 0.08650 | 0.05316 | 6.61394 | 0.06214 | 0.11587 | 0.05343
| (0.04175) | (0.04596) | (0.01490) | (0.00988) | (0.10294) | (0.01082) | (0.01077) | (0.00988)
$10^{-4}$ | 0.17094 | 0.13428 | 0.06448 | 0.04548 | 6.63878 | 0.04895 | 0.10504 | 0.04575
| (0.02052) | (0.03632) | (0.00990) | (0.00954) | (0.11971) | (0.00950) | (0.00990) | (0.00954)
$10^{-2}$ | 0.18817 | 0.13054 | 0.09701 | 0.04566 | 6.43301 | 0.05337 | 0.10362 | 0.04592
| (0.01896) | (0.03597) | (0.01584) | (0.00636) | (0.10094) | (0.00645) | (0.00648) | (0.00636)
$10^{-1}$ | 0.16236 | 0.16727 | 0.08413 | 0.05843 | 6.22776 | 0.06196 | 0.10845 | 0.05859
| (0.01630) | (0.03973) | (0.01257) | (0.01052) | (0.12413) | (0.01047) | (0.01117) | (0.01052)
$10^{-0.5}$ | 0.28054 | 0.27263 | 0.16329 | 0.14208 | 6.45716 | 0.14896 | 0.18881 | 0.14232
| (0.02537) | (0.06909) | (0.01635) | (0.01475) | (0.12854) | (0.01469) | (0.01497) | (0.01476)
$1$ | 1.05281 | 1.09944 | 1.01384 | 0.95356 | 6.77452 | 1.01920 | 0.95611 | 1.04263
| (0.02046) | (0.02191) | (0.00997) | (0.00905) | (0.10774) | (0.00877) | (0.00849) | (0.00715)
$10^{0.5}$ | 8.37061 | 9.48056 | 9.67063 | 4.71207 | 7.36306 | 9.67729 | 7.55671 | 9.57846
| (0.02177) | (0.06041) | (0.01659) | (0.05336) | (0.13437) | (0.01424) | (0.03817) | (0.00842)
$10$ | 72.86179 | 105.5150 | 99.00799 | 8.31751 | 8.80645 | 97.58522 | 39.24465 | 104.3057
| (0.12036) | (0.37105) | (0.82643) | (0.12924) | (0.13857) | (0.90360) | (0.25660) | (0.01449)
$10^{1.5}$ | 501.0324 | 995.3871 | 996.3600 | 8.88848 | 8.93730 | 996.3600 | 71.83093 | 913.9395
| (0.40598) | (0.97373) | (0.00678) | (0.13554) | (0.13665) | (0.00678) | (0.42595) | (0.21334)
$10^{2}$ | 2015.788 | 8484.242 | 7050.887 | 9.07261 | 9.07689 | 7228.744 | 21.54672 | 4691.845
| (1.37168) | (1.37785) | (12.25646) | (0.12175) | (0.12132) | (13.23414) | (0.24329) | (2.66192)

Table 17: The performance comparison table for each simulation setting in
Section 6.4 with $(n,p)=(200,300)$ and $\sigma=3$. The values are the mean
squared same-X errors averaged over 50 independent replications. The numbers
in parentheses are normalized sample standard deviations. The column labels
are defined in Section 6.4. Boldface indicates the best model. Underlined is
the best model within our proposed estimators 2n-Rep1 to 2n-Rep4.

Performance comparison table $(n,p)=(200,300)$
---
$p$ $\backslash$ $\tau$ | 0 | $10^{-4}$ | $10^{-2}$ | $10^{-1}$ | $10^{-0.5}$ | 1 | $10^{0.5}$ | $10$ | $10^{1.5}$ | $10^{2}$
$\sigma=1$ | 0.00620 | 0.00237 | 0.00494 | 0.01363 | 0.09847 | 0.48945 | 0.91546 | 1.00967 | 1.15949 | 0.98135
FVS2 | (31622.78) | (31622.78) | (5623.413) | (123.0269) | (9.33254) | (1.02329) | (0.10471) | (0.00977) | (0.00115) | (0.00019)
$\sigma=1$ | 0.00638 | 0.00255 | 0.00513 | 0.01373 | 0.10376 | 0.66386 | 0.96319 | 1.01786 | 1.14341 | 0.98173
Ridge | (31622.78) | (31622.78) | (31622.78) | (31622.78) | (5623.413) | (239.883) | (13.18257) | (1.148154) | (2.13796) | (0.89125)
$\sigma=2$ | 0.01499 | 0.01633 | 0.01702 | 0.03516 | 0.11523 | 0.76342 | 2.79710 | 3.78894 | 4.01031 | 3.95205
FVS2 | (31622.78) | (31622.78) | (27542.29) | (331.1311) | (37.15352) | (4.36516) | (0.42658) | (0.03715) | (0.00355) | (0.00028)
$\sigma=2$ | 0.01783 | 0.01892 | 0.01984 | 0.03759 | 0.11682 | 0.84901 | 3.43498 | 3.86065 | 4.02001 | 3.95274
Ridge | (31622.78) | (31622.78) | (31622.78) | (31622.738) | (31622.78) | (5128.614) | (79.43282) | (14.12538) | (2.13796) | (0.02692)
$\sigma=3$ | 0.04586 | 0.04670 | 0.03703 | 0.06753 | 0.12454 | 1.11509 | 4.63284 | 8.29820 | 8.93890 | 8.80025
FVS2 | (31622.78) | (31622.78) | (8709.636) | (851.1380) | (85.11380) | (9.12011) | (0.93325) | (0.08511) | (0.00977) | (0.00087)
$\sigma=3$ | 0.05893 | 0.05933 | 0.05001 | 0.08110 | 0.13283 | 1.15278 | 6.49951 | 8.72412 | 9.01420 | 8.80203
Ridge | (31622.78) | (31622.78) | (31622.78) | (31622.78) | (31622.78) | (14125.38) | (602.5596) | (27.54229) | (5.01187) | (3.38844)

Table 18: The values are the minimal average same-X loss values obtained from
each estimator (see (26) and (27)) within the set of
$\lambda\in\\{10^{-7+0.01k}:j=0,1,\ldots,1150\\}$ from the simulation settings
in Section 6.4 when $(n,p)=(200,300)$. The numbers in parentheses are the
values of $\lambda$ that achieve the minimum for each setting and each
estimator. Boldface indicates the best model for each $\tau$ and $\sigma$.

###  .6 Additional tables for Section 6.5

Performance comparison table when $\sigma=1$
---
$(s,{\rm Coding~{}\\#})$ | 2n-Or | Ridge | LASSO | 2n-G | 2n-Rep1
$s=0.025$ | 0.41999 | 0.34635 | 0.13778 | 0.49072 | 0.55543
Coding-1 | (0.00502) | (0.00589) | (0.00727) | (0.00962) | (0.00936)
$s=0.025$ | 0.41999 | 1.00610 | 1.14077 | 0.49072 | 0.52888
Coding-2 | (0.00502) | (0.01917) | (0.03211) | (0.00962) | (0.00898)
$s=0.05$ | 0.55923 | 0.46571 | 0.16670 | 0.71573 | 0.65167
Coding-1 | (0.00665) | (0.00816) | (0.00730) | (0.01901) | (0.00759)
$s=0.05$ | 0.55923 | 0.64650 | 0.80326 | 0.71573 | 0.64631
Coding-2 | (0.00665) | (0.01711) | (0.03086) | (0.01901) | (0.00757)
$s=0.1$ | 0.88849 | 0.73842 | 0.50612 | 1.02179 | 0.92288
Coding-1 | (0.01241) | (0.01289) | (0.01148) | (0.03451) | (0.01333)
$s=0.1$ | 0.88849 | 1.57656 | 1.62821 | 1.02179 | 0.92044
Coding-2 | (0.01241) | (0.03464) | (0.04985) | (0.03451) | (0.01331)
$s=0.2$ | 0.94443 | 0.86450 | 0.68938 | 1.00405 | 0.95970
Coding-1 | (0.01276) | (0.01517) | (0.01294) | (0.01984) | (0.01253)
$s=0.2$ | 0.94443 | 1.18902 | 1.31990 | 1.00405 | 0.95948
Coding-2 | (0.01276) | (0.02706) | (0.05416) | (0.01984) | (0.01253)
$s=0.3$ | 0.95895 | 0.87565 | 0.87934 | 0.98095 | 0.96919
Coding-1 | (0.01285) | (0.01319) | (0.01666) | (0.01405) | (0.01297)
$s=0.3$ | 0.95895 | 1.24615 | 1.11530 | 0.98095 | 0.96907
Coding-2 | (0.01285) | (0.02910) | (0.02987) | (0.01405) | (0.01297)

Table 19: The performance comparison table for each simulation setting in
Section 6.5 when $\sigma=1$. The values are the mean squared same-X errors
averaged over 50 independent replications. The numbers in parentheses are
normalized sample standard deviations. The column labels are defined in
Section 6.4. Boldface indicates the best model. Underlined is our main
proposed estimator: 2n-Rep1.

Performance comparison table when $\sigma=3$
---
$(s,{\rm Coding~{}\\#})$ | 2n-Or | Ridge | LASSO | 2n-G | 2n-Rep2 | 2n-Rep3
$s=0.01$ | 0.11933 | 0.24273 | 0.16385 | 0.14141 | 0.12694 | 0.17240
Coding-1 | (0.01243) | (0.02976) | (0.02155) | (0.01606) | (0.01302) | (0.01246)
$s=0.01$ | 0.11933 | 0.39756 | 0.31803 | 0.14141 | 0.12621 | 0.16576
Coding-2 | (0.01243) | (0.03321) | (0.03025) | (0.01606) | (0.01293) | (0.01199)
$s=0.02$ | 0.35239 | 0.39643 | 0.42977 | 0.37503 | 0.35854 | 0.37686
Coding-1 | (0.01035) | (0.01197) | (0.02694) | (0.01212) | (0.01053) | (0.01008)
$s=0.02$ | 0.35239 | 0.50265 | 0.51940 | 0.37503 | 0.35827 | 0.37556
Coding-2 | (0.01035) | (0.02094) | (0.02899) | (0.01212) | (0.01052) | (0.0098)
$s=0.03$ | 0.70534 | 0.71936 | 0.88481 | 0.74575 | 0.72611 | 0.71032
Coding-1 | (0.01334) | (0.02562) | (0.05025) | (0.01639) | (0.01363) | (0.01358)
$s=0.03$ | 0.70534 | 1.05775 | 1.14859 | 0.74575 | 0.72954 | 0.70970
Coding-2 | (0.01334) | (0.05330) | (0.05018) | (0.01639) | (0.01383) | (0.01400)
$s=0.04$ | 0.56631 | 0.59650 | 0.68032 | 0.59733 | 0.57926 | 0.57477
Coding-1 | (0.01107) | (0.02028) | (0.03621) | (0.01139) | (0.01099) | (0.01162)
$s=0.04$ | 0.56631 | 0.74298 | 0.88815 | 0.59733 | 0.58074 | 0.57244
Coding-2 | (0.01107) | (0.02746) | (0.06009) | (0.01139) | (0.01093) | (0.01130)
$s=0.05$ | 1.78976 | 1.60166 | 1.83978 | 1.94710 | 1.91890 | 1.85974
Coding-1 | (0.01647) | (0.05570) | (0.04532) | (0.02877) | (0.02369) | (0.01577)
$s=0.05$ | 1.78976 | 1.72413 | 2.22285 | 1.94710 | 1.92666 | 1.86978
Coding-2 | (0.01647) | (0.04972) | (0.05914) | (0.02877) | (0.02534) | (0.01907)

Table 20: The performance comparison table for each simulation setting in
Section 6.5 when $\sigma=3$. The values are the mean squared same-X errors
averaged over 50 independent replications. The numbers in parentheses are
normalized sample standard deviations. The column labels are defined in
Section 6.4. Boldface indicates the best model. Underlined is our main
proposed estimator: 2n-Rep3.

###  .7 Additional simulations: submodel shrinkage

We conduct additional simulations which is again the linear regression
subjects model (1) with $n=100,\,p=75$. The design matrix $X$ has ones in its
first column and independent draws from $N_{p-1}(0,\Sigma)$ in the remaining
entries on each row, where $\Sigma_{jk}=0.5^{|j-k|}$. However, in this section
we explore, possible scenarios where the shrinkage towards known submodel (or
submatrix) $X_{0}$ has an advantage over the shrinkage towards the intercept
only model (see Section 5 for details). Let $X_{0}$ denote the sub-matrix of
$X$ which consists of the first $p_{0}$ columns of $X$. We denote
$X=(X_{0},X_{1})$. We consider $p_{0}\in\\{5,10,15,20,25\\}$. And, we update
$X_{1}$ by $X_{1}^{*}=Q_{X_{0}}X_{1}=(I-P_{X_{0}})X_{1}$, and further replace
$(X_{0},X_{1})$ with $(X_{0},X_{1}^{*})$ for $X$. We conduct this
orthogonalization procedure to decompose terms easier (see (41) below). Then
we create

$\displaystyle\beta_{0}=X_{0}^{-}(1_{p_{0}}+\tau_{0}Z),$
$\displaystyle\beta_{1}=\tau_{1}(X_{1}^{*})^{-}Z^{\prime},$

where $Z,\,Z^{\prime}\sim N_{n}(0,I)$ and $Z,\,Z^{\prime}$ are independent of
each other. In this setting, we have

$\displaystyle\|\mu-P_{X_{0}}\mu\|^{2}$
$\displaystyle=\|X_{1}^{*}\beta_{1}\|^{2}=\tau_{1}^{2}\|P_{X_{1}^{*}}Z^{\prime}\|^{2},$
$\displaystyle\|\mu-P_{1}\mu\|^{2}$
$\displaystyle=\|(P_{X_{0}}-P_{1})(X_{0}\beta_{0})\|^{2}+\|X_{1}^{*}\beta_{1}\|^{2}$
$\displaystyle=\tau_{0}^{2}\|P_{X_{0}^{*}}Z\|^{2}+\tau_{1}^{2}\|P_{X_{1}^{*}}Z^{\prime}\|^{2},$
(41)

where $X_{0}^{*}=(P_{X_{0}}-P_{1})X_{0}$. We recall that the oracle tuning
parameters (5), (18) are given by

$\displaystyle\gamma_{\rm opt}=\frac{\|\mu-P_{1}\mu\|^{2}}{\sigma^{2}({\rm
rank}(X)-1)+\|\mu-P_{1}\mu\|^{2}},$ $\displaystyle\tilde{\gamma}_{\rm
opt}=\frac{\|\mu-P_{X_{0}}\mu\|^{2}}{\sigma^{2}({\rm rank}(X)-{\rm
rank}(X_{0}))+\|\mu-P_{X_{0}}\mu\|^{2}}.$

Thus, since $\sigma^{2}=1$, we have

$\displaystyle{\mathbb{E}}\left\|X\hat{\beta}^{(\gamma_{\rm
opt})}-X\beta\right\|^{2}$ $\displaystyle=\sigma^{2}+\frac{\sigma^{2}({\rm
rank}(X)-1)\|\mu-P_{1}\mu\|^{2}}{\sigma^{2}({\rm rank}(X)-1)+\|\mu-
P_{1}\mu\|^{2}}$ $\displaystyle=1+\frac{({\rm rank}(X)-1)\|\mu-
P_{1}\mu\|^{2}}{({\rm rank}(X)-1)+\|\mu-P_{1}\mu\|^{2}},$
$\displaystyle{\mathbb{E}}\left\|X\hat{\beta}^{(\tilde{\gamma}_{\rm
opt})}-X\beta\right\|^{2}$ $\displaystyle=\sigma^{2}{\rm
rank}(X_{0})+\frac{\sigma^{2}({\rm rank}(X)-{\rm rank}(X_{0}))\|\mu-
P_{X_{0}}\mu\|^{2}}{\sigma^{2}({\rm rank}(X)-{\rm rank}(X_{0}))+\|\mu-
P_{X_{0}}\mu\|^{2}}$ $\displaystyle={\rm rank}(X_{0})+\frac{({\rm
rank}(X)-{\rm rank}(X_{0}))\|\mu-P_{X_{0}}\mu\|^{2}}{({\rm rank}(X)-{\rm
rank}(X_{0}))+\|\mu-P_{X_{0}}\mu\|^{2}}.$

We control $\tau_{1}$ to satisfy the following equation:

$\displaystyle\|\mu-P_{X_{0}}\mu\|^{2}$
$\displaystyle=\tau_{1}^{2}\|P_{X_{1}^{*}}Z^{\prime}\|^{2}=R_{1}({\rm
rank}(X)-{\rm rank}(X_{0})),$ $\displaystyle\|\mu-P_{1}\mu\|^{2}$
$\displaystyle=\tau_{0}^{2}\|P_{X_{0}^{*}}Z\|^{2}+\tau_{1}^{2}\|P_{X_{1}^{*}}Z^{\prime}\|^{2}$
$\displaystyle=\tau_{0}^{2}\|P_{X_{0}^{*}}Z\|^{2}+R_{1}({\rm rank}(X)-{\rm
rank}(X_{0}))$ $\displaystyle=R_{2}({\rm rank}(X)-1).$

Thus, ${\mathbb{E}}\left\|X\hat{\beta}^{(\gamma_{\rm
opt})}-X\beta\right\|^{2}>{\mathbb{E}}\left\|X\hat{\beta}^{(\tilde{\gamma}_{\rm
opt})}-X\beta\right\|^{2}$ holds when

$\displaystyle\frac{{\rm rank}(X)-1}{R_{2}+1}>{\rm rank}(X_{0})+\frac{{\rm
rank}(X)-{\rm rank}(X_{0})}{R_{1}+1}.$

We vary $R_{1}\in\\{10^{-1},10^{-1/2},1,10^{1/2},10\\}$ and
$R_{2}\in\\{10^{-1/2},1,10^{1/2},10,10^{3/2}\\}$ while keeping $R_{2}>R_{1}$.
This implies that there are 15 pairs of $(R_{1},R_{2})$. Smaller $R_{1}$
implies smaller (expected) realization of $\|\mu-P_{X_{0}}\mu\|^{2}$. Smaller
$R_{2}$ indicates smaller (expected) realization of $\|\mu-
P_{1}\mu\|^{2}-\|\mu-P_{X_{0}}\mu\|^{2}$. We include 2n-Or, 2n-Orsb, 2n-Es,
2n-Es95, 2n-Essb, 2n-Essb95, R, L, R-sb, L-sb as the competitors. In the set
of competitors, 2n-Orsb exploits the oracle tuning parameter for the submodel
shrinkage towards $X_{0}$ (see (18) for details). 2n-Essb95 is based on the
estimator for $\tilde{\gamma}_{\rm opt}$ through

$\displaystyle\hat{\tilde{\gamma}}_{95}=(1-1/\tilde{F})\cdot
1(\tilde{F}\geq\tilde{f}_{0.95}),$

where $\tilde{f}_{0.95}$ is the 0.95 quantile of the central F-distribution
with degrees of freedom ${\rm rank}(X)-{\rm rank}(X_{0})$ and $n-{\rm
rank}(X)$ (see (19) for the definition of $\tilde{F}$). And, R-sb and L-sb
does not shrink the coefficients according to $X_{0}$, which are fitted by
glmnet as for R and L.

The results are shown in Table 21–25. In general, 2n-Orsb was the best among
the candidates except for few cases where 2n-Or was the best. When
$R_{2}=10^{1/2}R_{1}$ for each $R_{1}$ and $p$, R-sb followed the next. And,
2n-Essb95 was much competitive than 2n-Essb. However, except for these cases,
2n-Essb was generally the best estimator among the candidates set without the
oracle submodel shrinkage 2n-Orsb.

Performance comparison table ($p_{0}=5$)
---
$(R_{1},R_{2})$ | OLS | 2n-Or | 2n-Orsb | 2n-Es | 2n-Essb | 2n-Es95 | 2n-Essb95 | R | R-sb | L | L-sb
$(10^{-1},10^{-1/2})$ | 0.75315 | 0.18980 | 0.11566 | 0.22162 | 0.13632 | 0.25021 | 0.13246 | 0.23288 | 0.13411 | 0.20668 | 0.14579
| (0.01806) | (0.00322) | (0.00450) | (0.00632) | (0.00749) | (0.00399) | (0.00708) | (0.00801) | (0.00820) | (0.01131) | (0.01161)
$(10^{-1},1)$ | 0.75391 | 0.37009 | 0.11780 | 0.40247 | 0.14407 | 0.50593 | 0.13321 | 0.49804 | 0.14304 | 0.28169 | 0.14074
| (0.01409) | (0.00671) | (0.00400) | (0.00987) | (0.00858) | (0.02553) | (0.00837) | (0.01675) | (0.01016) | (0.01163) | (0.00779)
$(10^{-1},10^{1/2})$ | 0.74555 | 0.55969 | 0.11745 | 0.58596 | 0.13973 | 0.58596 | 0.12298 | 0.46691 | 0.13354 | 0.23834 | 0.13071
| (0.01941) | (0.01321) | (0.00465) | (0.01470) | (0.00663) | (0.01470) | (0.00450) | (0.01297) | (0.00650) | (0.01203) | (0.00653)
$(10^{-1},10)$ | 0.74243 | 0.66124 | 0.11471 | 0.66965 | 0.14174 | 0.66965 | 0.13262 | 0.47052 | 0.12693 | 0.19610 | 0.13456
| (0.01805) | (0.01599) | (0.00478) | (0.01557) | (0.00858) | (0.01557) | (0.00823) | (0.01317) | (0.00550) | (0.00965) | (0.00643)
$(10^{-1},10^{3/2})$ | 0.75591 | 0.73121 | 0.11318 | 0.73492 | 0.15318 | 0.73492 | 0.13518 | 0.62439 | 0.12744 | 0.23151 | 0.13621
| (0.01977) | (0.01879) | (0.00419) | (0.01881) | (0.00952) | (0.01881) | (0.00868) | (0.01493) | (0.00585) | (0.00887) | (0.00706)
$(10^{-1/2},1)$ | 0.74920 | 0.38666 | 0.22035 | 0.42746 | 0.25591 | 0.53655 | 0.28396 | 0.39515 | 0.26311 | 0.31882 | 0.28338
| (0.01878) | (0.00856) | (0.00644) | (0.01359) | (0.01025) | (0.02627) | (0.00859) | (0.01157) | (0.00701) | (0.01086) | (0.00915)
$(10^{-1/2},10^{1/2})$ | 0.73399 | 0.55839 | 0.20824 | 0.57860 | 0.23963 | 0.57860 | 0.26789 | 0.52743 | 0.26231 | 0.36967 | 0.27402
| (0.01736) | (0.01335) | (0.00431) | (0.01347) | (0.00601) | (0.01347) | (0.00423) | (0.01421) | (0.00519) | (0.01252) | (0.00577)
$(10^{-1/2},10)$ | 0.77899 | 0.70325 | 0.22648 | 0.71316 | 0.26098 | 0.71316 | 0.27987 | 0.66012 | 0.26736 | 0.41049 | 0.27625