Parameter estimation in linear models

Illustration of the linear model in high-dimensions: a data set consists of a response vector ${\displaystyle Y\in \mathbb {R} ^{n}}$ and a design matrix ${\displaystyle X\in \mathbb {R} ^{n\times p}}$ with ${\displaystyle p\gg n}$. Our goal is to estimate the unknown vector ${\displaystyle \beta =(\beta _{1},\dots ,\beta _{p})\in \mathbb {R} ^{p}}$ of regression coefficients where ${\displaystyle \beta }$ is often assumed to be sparse, in the sense that the cardinality of the set ${\displaystyle S:=\{j:\beta _{j}\neq 0\}}$ is small by comparison with ${\displaystyle p}$.

The most basic statistical model for the relationship between a covariate vector ${\displaystyle x\in \mathbb {R} ^{p}}$ and a response variable ${\displaystyle y\in \mathbb {R} }$ is the linear model

${\displaystyle y=x^{\top }\beta +\epsilon ,}$

where ${\displaystyle \beta \in \mathbb {R} ^{p}}$ is an unknown parameter vector, and ${\displaystyle \epsilon }$ is random noise with mean zero and variance ${\displaystyle \sigma ^{2}}$. Given independent responses ${\displaystyle Y_{1},\ldots ,Y_{n}}$, with corresponding covariates ${\displaystyle x_{1},\ldots ,x_{n}}$, from this model, we can form the response vector ${\displaystyle Y=(Y_{1},\ldots ,Y_{n})^{\top }}$, and design matrix ${\displaystyle X=(x_{1},\ldots ,x_{n})^{\top }\in \mathbb {R} ^{n\times p}}$. When ${\displaystyle n\geq p}$ and the design matrix has full column rank (i.e. its columns are linearly independent), the ordinary least squares estimator of ${\displaystyle \beta }$ is

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

When ${\displaystyle \epsilon \sim N(0,\sigma ^{2})}$, it is known that ${\displaystyle {\hat {\beta }}\sim N_{p}{\bigl (}\beta ,\sigma ^{2}(X^{\top }X)^{-1}{\bigr )}}$. Thus, ${\displaystyle {\hat {\beta }}}$ is an unbiased estimator of ${\displaystyle \beta }$, and the Gauss-Markov theorem tells us that it is the Best Linear Unbiased Estimator.

However, overfitting is a concern when ${\displaystyle p}$ is of comparable magnitude to ${\displaystyle n}$: the matrix ${\displaystyle X^{\top }X}$ in the definition of ${\displaystyle {\hat {\beta }}}$ may become ill-conditioned, with a small minimum eigenvalue. In such circumstances ${\displaystyle \mathbb {E} (\|{\hat {\beta }}-\beta \|^{2})=\sigma ^{2}\mathrm {tr} {\bigl (}(X^{\top }X)^{-1}{\bigr )}}$ will be large (since the trace of a matrix is the sum of its eigenvalues). Even worse, when ${\displaystyle p>n}$, the matrix ${\displaystyle X^{\top }X}$ is singular. (See Section 1.2 and Exercise 1.2 in ^[1].)

It is important to note that the deterioration in estimation performance in high dimensions observed in the previous paragraph is not limited to the ordinary least squares estimator. In fact, statistical inference in high dimensions is intrinsically hard, a phenomenon known as the curse of dimensionality, and it can be shown that no estimator can do better in a worst-case sense without additional information (see Example 15.10^[2]). Nevertheless, the situation in high-dimensional statistics may not be hopeless when the data possess some low-dimensional structure. One common assumption for high-dimensional linear regression is that the vector of regression coefficients is sparse, in the sense that most coordinates of ${\displaystyle \beta }$ are zero. Many statistical procedures, including the Lasso, have been proposed to fit high-dimensional linear models under such sparsity assumptions.

Covariance matrix estimation

Another example of a high-dimensional statistical phenomenon can be found in the problem of covariance matrix estimation. Suppose that we observe ${\displaystyle X_{1},\ldots ,X_{n}\in \mathbb {R} ^{p}}$, which are i.i.d. draws from some zero mean distribution with an unknown covariance matrix ${\displaystyle \Sigma \in \mathbb {R} ^{p\times p}}$. A natural unbiased estimator of ${\displaystyle \Sigma }$ is the sample covariance matrix

${\displaystyle {\widehat {\Sigma }}:={\frac {1}{n}}\sum _{i=1}^{n}X_{i}X_{i}^{\top }.}$

In the low-dimensional setting where ${\displaystyle n}$ increases and ${\displaystyle p}$ is held fixed, ${\displaystyle {\widehat {\Sigma }}}$ is a consistent estimator of ${\displaystyle \Sigma }$ in any matrix norm. When ${\displaystyle p}$ grows with ${\displaystyle n}$, on the other hand, this consistency result may fail to hold. As an illustration, suppose that each ${\displaystyle X_{i}\sim N_{p}(0,I)}$ and ${\displaystyle p/n\rightarrow \alpha \in (0,1)}$. If ${\displaystyle {\widehat {\Sigma }}}$ were to consistently estimate ${\displaystyle \Sigma =I}$, then the eigenvalues of ${\displaystyle {\widehat {\Sigma }}}$ should approach one as ${\displaystyle n}$ increases. It turns out that this is not the case in this high-dimensional setting. Indeed, the largest and smallest eigenvalues of ${\displaystyle {\widehat {\Sigma }}}$ concentrate around ${\displaystyle (1+{\sqrt {\alpha }})^{2}}$ and ${\displaystyle (1-{\sqrt {\alpha }})^{2}}$, respectively, according to the limiting distribution derived by Tracy and Widom, and these clearly deviate from the unit eigenvalues of ${\displaystyle \Sigma }$. Further information on the asymptotic behaviour of the eigenvalues of ${\displaystyle {\widehat {\Sigma }}}$ can be obtained from the Marchenko–Pastur law. From a non-asymptotic point of view, the maximum eigenvalue ${\displaystyle \lambda _{\mathrm {max} }({\widehat {\Sigma }})}$ of ${\displaystyle {\widehat {\Sigma }}}$ satisfies

${\displaystyle \mathbb {P} \left(\lambda _{\mathrm {max} }({\widehat {\Sigma }})\geq (1+{\sqrt {p/n}}+\delta )^{2}\right)\leq e^{-n\delta ^{2}/2},}$

for any ${\displaystyle \delta \geq 0}$ and all choices of pairs of ${\displaystyle n,p}$.^[2]

Again, additional low-dimensional structure is needed for successful covariance matrix estimation in high dimensions. Examples of such structures include sparsity, low rankness and bandedness. Similar remarks apply when estimating an inverse covariance matrix (precision matrix).