ST 503 – Fundamentals of Linear Models and Regression

Linear algebra is the study of vector spaces, linear transformations, matrices, and inner product spaces.

Definition. A vector space $V$ over a field $F$ is a set of elements $v \in V$, with the operations of addition and scalar multiplication, that satisfy the following axioms.

$\forall\, x, y \in V,\quad x + y = y + x$
$\forall\, x, y, z \in V,\quad (x + y) + z = x + (y + z)$
There exists an element “$0$” $\in V$ s.t. $x + 0 = x$, $\forall\, x \in V$
$\forall\, x \in V\; \exists\, y \in V$ s.t. $x + y = 0$
There exists an element “$1$” $\in F$ s.t. $1 \cdot x = x$, $\forall\, x \in V$
$\forall\, a, b \in F,\; \forall\, x \in V,\quad (ab)x = a(bx)$
$\forall\, a \in F,\; \forall\, x, y \in V,\quad a(x + y) = ax + ay$
$\forall\, a, b \in F,\; \forall\, x \in V,\quad (a + b)x = ax + bx$ $\square$

In statistics contexts, typically $V = \mathbb{R}^n$ and $F = \mathbb{R}$.

The notion of a vector space is truly foundational for data science. A vector space gives us an entire framework and intuition for thinking about data objects (even as data are stored as objects in a computer).

Definition. A collection of vectors $x_1, \ldots, x_p$ are said to be linearly dependent if and only if there exist scalars $c_1, \ldots, c_p$, not all zero, such that

\sum_{i=1}^{p} c_i x_i = 0.

A set of vectors that are *not* linearly dependent are said to be linearly independent. $\square$

Often, to show that a set of vectors are linearly independent, a good strategy is to demonstrate that $c_1 x_1 + \cdots + c_p x_p = 0$ implies that $c_1 = \cdots = c_p = 0$.

Note that if $x_1, \ldots, x_p \in \mathbb{R}^n$, then there can be at most $\min\{n, p\}$ linearly independent of the $p$ vectors.

Definition. A basis, $B$, for a vector space, $V$, is a linearly independent subset of $V$ with $\text{span}(B) = V.$ $\square$

Example. For $V = \mathbb{R}^n$, the set $\{e_1, \ldots, e_n\}$, where $e_i \in \mathbb{R}^n$ has $i$th component with value 1 and 0 for all other components, is a basis. This is commonly called the standard basis, and $e_i$ a standard basis vector. $\square$

Theorem. Let $V$ be a vector space and $B := \{u_1, \ldots, u_n\} \subseteq V$. Then $B$ is a basis for $V$ if and only if every $v \in V$ can be expressed uniquely as

v = a_1 u_1 + \cdots + a_n u_n

for some scalars $a_1, \ldots, a_n$. $\square$

Proof. Left as an exercise. $\square$

Theorem. If a vector space $V$ is generated by some finite set $S$, then some subset of $S$ is a basis for $V$. $\square$

Proof. Left as an exercise. $\square$

Definition. Let $V$ and $W$ be vector spaces over field $\mathbb{R}$. We call a function $T : V \to W$ a linear transformation from $V$ to $W$ if $\forall\, x, y \in V$ and $\forall\, c \in \mathbb{R}$,

$T(x + y) = T(x) + T(y)$
$T(cx) = cT(x).$ $\square$

Example. Let $A \in \mathbb{R}^{p \times q}$ and define $T(x) := Ax$ $\forall\, x \in \mathbb{R}^q$. Then $T : \mathbb{R}^q \to \mathbb{R}^p$ is a linear transformation since

$\forall\, x, y \in \mathbb{R}^q,\quad T(x + y) = A(x + y) = Ax + Ay = T(x) + T(y)$, and
$\forall\, x \in \mathbb{R}^q,\; \forall\, c \in \mathbb{R},\quad T(cx) = A \cdot (cx) = cAx = cT(x).$ $\square$

Example. Let $V = C(\mathbb{R})$, the vector space of continuous real-valued functions defined on $\mathbb{R}$. Then $\forall\, a, b \in \mathbb{R}$ with $a < b$, the transformation $T : V \to \mathbb{R}$ defined by,

T(f) := \int_a^b f(x)\, dx

is linear. $\square$

Definition. Let $V$ and $W$ be vector spaces, and $T : V \to W$ a linear transformation. The set,

\text{null}(T) := \{x \in V : T(x) = 0\},

is called the null space of $T$. The range of $T$ is defined as,

\text{range}(T) := \{y \in W : T(x) = y \text{ for some } x \in V\}.

Further, $\text{rank}(T) := \dim(\text{range}(T))$. $\square$

Theorem. (dimension theorem) Let $V$ and $W$ be vector spaces, and $T : V \to W$ be linear. If $V$ is finite-dimensional, then

\dim(\text{null}(T)) + \text{rank}(T) = \dim(V).

$\square$

Proof. See any linear algebra textbook. $\square$

Recall that if $T(x) := Ax$, then $\text{range}(T)$ is called the column space of $A$, denoted $\text{col}(A)$. Similarly for $\text{range}(T')$, called the row space of $A$, and denoted by $\text{row}(A)$.

Definition. Let $A \in \mathbb{R}^{p \times q}$ and $B \in \mathbb{R}^{q \times m}$. Then the matrix product, denoted by $AB$, has components defined as

(AB)_{ij} := \sum_{k=1}^{q} A_{ik} B_{kj}

for $i \in \{1, \ldots, p\}$ and $j \in \{1, \ldots, m\}$. $\square$

Further notions for matrices.

The transpose of $A$, denoted $A'$, satisfies the property that $(A')_{ij} = A_{ji}$ for all indices $i, j$.
$A \in \mathbb{R}^{p \times q}$ is called square if $p = q$.
A square matrix is called symmetric if $A' = A$.
$A \in \mathbb{R}^{p \times p}$ is said to be invertible if $\exists\, B \in \mathbb{R}^{p \times p}$ s.t. $AB = I_p = BA$. In that case, $A^{-1} := B$.
$A \in \mathbb{R}^{p \times p}$, the trace of $A$ is defined as $\displaystyle\text{tr}(A) := \sum_{i=1}^{p} A_{ii}$.

Selected properties of matrices.

$(AB)' = B'A'$
$A \in \mathbb{R}^{p \times p}$ is invertible if and only if $\text{rank}(A) = p$ (equivalently $\det(A) \neq 0$).
If both $A, B \in \mathbb{R}^{p \times p}$ are invertible, then $(AB)^{-1} = B^{-1}A^{-1}$.
For $A$ and $B$ of appropriate dimensions, $\text{tr}(AB) = \text{tr}(BA)$.

Definition. Let $A \in \mathbb{R}^{p \times p}$. A nonzero vector $v \in \mathbb{R}^p$ is called an eigenvector of $A$ if $Av = \lambda v$ for some scalar $\lambda$, called an eigenvalue. $\square$

Theorem. A scalar $\lambda$ is an eigenvalue of a matrix $A \in \mathbb{R}^{p \times p}$ if and only if $\det(A - \lambda I_p) = 0$. The polynomial $f(t) := \det(A - tI_p)$ is called the characteristic polynomial of $A$. $\square$

Proof. See any linear algebra textbook. $\square$

For example, if

A = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix},

then $\forall\, v \in \mathbb{R}^3$, $Av = 1 \cdot v$. How about

A = \begin{pmatrix} 1 & 2 & 3 \\ 0 & 4 & 5 \\ 0 & 0 & 6 \end{pmatrix}?

f(t) = \det(A - tI_3) \] \[ = \det \begin{pmatrix} 1-t & 2 & 3 \\ 0 & 4-t & 5 \\ 0 & 0 & 6-t \end{pmatrix} \] \[ = (1-t)\det\begin{pmatrix} 4-t & 5 \\ 0 & 6-t \end{pmatrix} - 0 + 0 \] \[ = (1-t)\left[(4-t)(6-t) - 0\right] \] \[ = (1-t)(4-t)(6-t)

which is a cubic polynomial with roots $t = 1$, $t = 4$, $t = 6$.

How to find the eigenvectors of

A = \begin{pmatrix} 1 & 2 & 3 \\ 0 & 4 & 5 \\ 0 & 0 & 6 \end{pmatrix},

and verify that $1, 4, 6$ are the eigenvalues of $A$?

Observe that the eigenvectors must satisfy the system of linear equations

Av = \lambda v \quad \text{or} \quad (A - \lambda I)v = 0

\begin{aligned} (1 - \lambda)v_1 + 2v_2 + 3v_3 &= 0 \\ (4 - \lambda)v_2 + 5v_3 &= 0 \\ (6 - \lambda)v_3 &= 0 \end{aligned}

For $\lambda = 6$,

\left. \begin{aligned} -5v_1 + 2v_2 + 3v_3 &= 0 \\ -2v_2 + 5v_3 &= 0 \\ 0 &= 0 \end{aligned} \right\} \text{WLOG take } v_3 = 1

v_2 = \tfrac{5}{2} \] \[ v_1 = \tfrac{2}{5}v_2 + \tfrac{3}{5} = \tfrac{8}{5}

To verify,

A \cdot v = \begin{pmatrix} 1 & 2 & 3 \\ 0 & 4 & 5 \\ 0 & 0 & 6 \end{pmatrix} \begin{pmatrix} 8/5 \\ 5/2 \\ 1 \end{pmatrix} = \begin{pmatrix} 48/5 \\ 15 \\ 6 \end{pmatrix} = 6 \cdot \begin{pmatrix} 8/5 \\ 5/2 \\ 1 \end{pmatrix} = \lambda v

In general, for a system of linear equations in row echelon form,

\begin{pmatrix} M_{11} & M_{12} & \cdots & M_{1p} \\ 0 & M_{22} & \cdots & M_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & M_{pp} \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_p \end{pmatrix} = \begin{pmatrix} b_1 \\ b_2 \\ \vdots \\ b_p \end{pmatrix}

Suppose $M_{ij} \neq 0$ $\forall\, j \geq i$. Then

x_p = \frac{b_p}{M_{pp}} \] \[ \vdots \] \[ x_2 = \left(b_2 - \sum_{j=3}^{p} M_{2j} x_j\right) \frac{1}{M_{22}} \] \[ x_1 = \left(b_1 - \sum_{j=2}^{p} M_{1j} x_j\right) \frac{1}{M_{11}}

In general, this is called back substitution

x_i = \left(b_i - \sum_{j=i+1}^{p} M_{ij} x_j\right) \frac{1}{M_{ii}} \qquad \text{for } i \in \{1, \ldots, p\}

Eigenvalues have an important role in matrix decompositions for certain classes of matrices. In particular, consider the following fundamental result for the class of symmetric matrices

Spectral Theorem. For any symmetric matrix $A \in \mathbb{R}^{p \times p}$ there exists an orthogonal matrix $Q$ (i.e., $Q'Q = QQ' = I_p$) such that,

\[ A = QDQ', \]

where $D$ is a diagonal matrix composed of the eigenvalues of $A$. $\square$

Note also that the columns of $Q$ form an orthonormal basis for $\mathbb{R}^p$ consisting of eigenvectors of $A$.

The basis part means that every vector $v \in \mathbb{R}^p$ can be expressed as

v = \sum_{i=1}^{p} Q_i a_i,

for some $a_1, \ldots, a_p \in \mathbb{R}$, and that $Q_1, \ldots, Q_p$ are linearly independent.

Recall the definition of linear independence:

Vectors $u_1, \ldots, u_n \in \mathbb{R}^p$ are said to be linearly independent if

\sum_{i=1}^{n} u_i a_i = 0 \quad \text{implies that} \quad a_1 = \cdots = a_n = 0

A set of vectors are said to be linearly dependent if they are not linearly independent.

Orthonormal means $\forall\, i, j \in \{1, \ldots, p\}$

Q_i' Q_j = \begin{cases} 1 & \text{if } i = j \\ 0 & \text{if } i \neq j \end{cases}.

Next, verify that the columns $Q_1, \ldots, Q_p$ are eigenvectors of $A$ corresponding to eigenvalues given by the diagonal components of $D$. Take any $i \in \{1, \ldots, p\}$:

AQ_i = QDQ'Q_i = QD \begin{pmatrix} 0 \\ \vdots \\ 1 \\ \vdots \\ 0 \end{pmatrix} = Q \begin{pmatrix} 0 \\ \vdots \\ D_{ii} \\ \vdots \\ 0 \end{pmatrix} = D_{ii} \cdot Q_i

Thus, $Q_i$ is an eigenvector of $A$, by definition, and by the same definition $D_{ii}$ is an eigenvalue.

Further properties: $\displaystyle\text{tr}(A) = \sum_{i=1}^{p} \lambda_i$ and $\displaystyle\det(A) = \prod_{i=1}^{p} \lambda_i$.

\text{tr}(A) = \text{tr}(QDQ') = \text{tr}(Q'QD) = \text{tr}(D) = \sum_{i=1}^{p} D_{ii}

\det(A) = \det(QDQ') = \det(Q) \cdot \det(D) \cdot \det(Q') = \det(QQ') \cdot \det(D) = 1 \cdot \prod_{i=1}^{p} D_{ii}

What if a matrix is not symmetric nor even square? How to decompose?

Let $X \in \mathbb{R}^{n \times p}$. Then $X'X \in \mathbb{R}^{p \times p}$ and $XX' \in \mathbb{R}^{n \times n}$, and both are symmetric:

(X'X)' = X'(X')' = X'X \] \[ (XX')' = (X')'\, X' = XX'

Accordingly, we can apply the Spectral theorem to both $X'X$ and $XX'$ to decompose as

X'X = V\Delta_1 V' \qquad \text{and} \qquad XX' = U\Delta_2 U'

for matrices $V \in \mathbb{R}^{p \times p}$ and $U \in \mathbb{R}^{n \times n}$, both orthogonal matrices, and diagonal matrices $\Delta_1$ and $\Delta_2$.

Moreover, it can be shown that $X = U\Sigma V'$, where $\Sigma$ is a diagonal matrix with nonzero diagonal components equal to the square roots of the nonzero diagonal components of both $\Delta_1$ and $\Delta_2$. The diagonal components of $\Sigma$ are called the singular values of $X$. For symmetric matrices $A \in \mathbb{R}^{p \times p}$ with singular values $\sigma_1, \ldots, \sigma_p$ and eigenvalues $\lambda_1, \ldots, \lambda_p$,

\sigma_i = |\lambda_i| \quad \text{for } i \in \{1, \ldots, p\}

SVD and other decompositions are important for data compression and dimension reduction. Consider this for $X \in \mathbb{R}^{n \times p}$ with $\text{rank}(X) = r$,

\begin{aligned} X &= U\Sigma V' = (U_1, \ldots, U_n) \begin{pmatrix} \sigma_1 & & & \\ & \ddots & & \\ & & \sigma_r & \\ & & & 0 \\ & & & & \ddots \\ & & & & & 0 \end{pmatrix} \begin{pmatrix} V_1' \\ \vdots \\ V_p' \end{pmatrix} \\[1.5em] &= (U_1, \ldots, U_r) \begin{pmatrix} \sigma_1 & & \\ & \ddots & \\ & & \sigma_r \end{pmatrix} \begin{pmatrix} V_1' \\ \vdots \\ V_r' \end{pmatrix} \quad \text{“compact SVD”} \end{aligned}

where $U_1, \ldots, U_n \in \mathbb{R}^n$ and $V_1, \ldots, V_r \in \mathbb{R}^p$. Storing raw $X$ requires $n \cdot p$ pieces of information. Storing the SVD requires

n \cdot r + r + p \cdot r = (n + p + 1) \cdot r

pieces of information. For example, suppose $n = 100$, $p = 30$. Then

(n + p + 1) \cdot r = 131 \cdot r < 3000 = n \cdot p \quad \text{if } r < 23.

Furthermore, if $\sigma_k, \ldots, \sigma_r \approx 0$, for some $k \leq r$, then

X \approx (U_1, \ldots, U_k) \begin{pmatrix} \sigma_1 & & \\ & \ddots & \\ & & \sigma_k \end{pmatrix} \begin{pmatrix} V_1' \\ \vdots \\ V_k' \end{pmatrix}

which allows for even less memory requirements.

Least Squares Problems

Suppose that $X \in \mathbb{R}^{n \times p}$ and $y \in \mathbb{R}^n$. Consider the conditions in which the system

\[ Xb = y \]

has a solution for some vector $b \in \mathbb{R}^p$. As we have seen, to determine whether $y \in \text{col}(X)$,

Gaussian elimination to express system in row echelon form.
Back substitution algorithm

In the special case that $p = n$ and $\text{rank}(X) = n$, $b = X^{-1}y$

How about the system when $p < n$?

$X$ is tall and narrow (i.e., more rows than columns)
Likely, $y \notin \text{col}(X)$
This is the context of least squares problems

If the system $Xb = y$ does not have a solution then perhaps the next best thing would be to find the $b \in \mathbb{R}^p$ which makes $Xb$ as “close” as possible to $y$.

Q(b) := \|y - Xb\|^2.

The $\underset{b}{\text{argmin}}\{Q(b)\}$ is called the least squares solution.

Desirable properties of $Q(b)$:

Convex
Differentiable
$\underset{b}{\text{argmin}}\{Q(b)\} = \{b : \nabla_b Q(b) = 0\}$

Recall that for a function $f : \mathbb{R}^p \to \mathbb{R}$,

\nabla_x f(x) := \begin{pmatrix} \dfrac{\partial f(x)}{\partial x_1} \\[6pt] \vdots \\[2pt] \dfrac{\partial f(x)}{\partial x_p} \end{pmatrix} \in \mathbb{R}^p.

Lemma. For $a, b \in \mathbb{R}^p$ and $A \in \mathbb{R}^{p \times p}$, the following properties hold.

$\nabla_b(a'b) = a$
$\nabla_b(b'Ab) = (A + A')b$ $\square$

Proof. See any multivariate calculus textbook. $\square$

Observing the fact that

Q(\beta) = (y - X\beta)'(y - X\beta) = y'y - y'X\beta - \beta'X'y + \beta'X'X\beta,

it follows by the lemma that

\nabla_\beta Q(\beta) = -2X'y + (X'X + X'X)\beta = -2X'(y - X\beta).

Setting $\nabla_\beta Q(\beta) = 0$ yields the “normal equations”

X'X\beta = X'y.

Example. Recall a simple linear regression,

y_i = \beta_0 + \beta_1 x_i + u_i \qquad \text{for } i \in \{1, \ldots, n\}.

In this case,

y = X\beta + u, \qquad \text{where} \quad X := \begin{pmatrix} 1 & x_1 \\ \vdots & \vdots \\ 1 & x_n \end{pmatrix}.

Then

X'X\beta = \begin{pmatrix} n & \displaystyle\sum_i x_i \\[6pt] \displaystyle\sum_i x_i & \displaystyle\sum_i x_i^2 \end{pmatrix} \beta \qquad \text{and} \qquad X'y = \begin{pmatrix} \displaystyle\sum_i y_i \\[6pt] \displaystyle\sum_i x_i y_i \end{pmatrix}.

Accordingly, the normal equations are

n\beta_0 + n\bar{x}_n \beta_1 = n\bar{y}_n \] \[ n\bar{x}_n \beta_0 + \sum_i x_i^2\, \beta_1 = \sum_i x_i y_i,

and so the least squares solution for $\beta_0$ and $\beta_1$ are

\begin{aligned} \hat{\beta}_0 &= \bar{y}_n - \bar{x}_n \hat{\beta}_1 \\[1em] \hat{\beta}_1 &= \frac{\displaystyle\sum_i y_i(x_i - \bar{x}_n)}{\displaystyle\sum_i (x_i - \bar{x}_n)^2}. \end{aligned}

$\square$

Moreover, it can be mathematically proven that there always exists a solution to the normal equations. In particular, when $\text{rank}(X) = p$, the unique least squares solution has the form

\hat{b} = (X'X)^{-1} X'y.

This solution turns out to be the coefficients corresponding to the orthogonal projection of $y$ onto $\text{col}(X)$:

\hat{y} = X\hat{b} = \underbrace{X(X'X)^{-1}X'}_{=:\, \mathcal{P}_X} y

A linear transformation $\mathcal{P} : \mathbb{R}^n \to \mathbb{R}^n$ is said to be a projection onto some subspace $V \subseteq \mathbb{R}^n$ if

$\mathcal{P}v \in V$ for every $v \in \mathbb{R}^n$
$\mathcal{P}v = v$ for every $v \in V.$

Further, $\mathcal{P}$ is called an orthogonal projection onto $V$ if it also satisfies

$(\mathcal{P}v)'(u - \mathcal{P}u) = (v - \mathcal{P}v)'\mathcal{P}u = 0$ $\forall\, v, u \in \mathbb{R}^n$.

Observe that for a projection matrix $\mathcal{P}$, of any $v \in \mathbb{R}^n$

\mathcal{P}^2 v = \mathcal{P}\mathcal{P}v = \mathcal{P}v \qquad \text{since } \mathcal{P}v \in V;

$(\mathcal{P}v)'(v - \mathcal{P}v) = 0$

$\mathcal{P}$ is “idempotent”. In fact, any idempotent matrix is a projection matrix, and any symmetric idempotent matrix is an orthogonal projection matrix.

Verify that $\mathcal{P}_X$ is an orthogonal projection matrix onto $\text{col}(X)$.

Suppose $v \in \mathbb{R}^n$. Then
$\mathcal{P}_X v = X(X'X)^{-1}X'v \in \text{col}(X)$
Suppose $v \in \text{col}(X)$. Then $v = Xa$ for some $a \in \mathbb{R}^p$, and so
$\mathcal{P}_X v = X(X'X)^{-1}X'v = X(X'X)^{-1}X'Xa = Xa = v$
Let $v, u \in \mathbb{R}^n$. Then
$(\mathcal{P}_X v)'(u - \mathcal{P}_X u) = v'(\mathcal{P}_X u - \mathcal{P}_X \mathcal{P}_X u) = 0 \] \[ (v - \mathcal{P}_X v)'\mathcal{P}_X u = (\mathcal{P}_X' v - \mathcal{P}_X' \mathcal{P}_X v)'\, u = 0$
since $\mathcal{P}_X' = \bigl(X(X'X)^{-1}X'\bigr)' = X(X'X)^{-1}X' = \mathcal{P}_X$

Example. Let $X = \mathbf{1}_n$. Then $\text{col}(X)$ is the subspace of vectors that have all components equal to the same value. That being true, projecting onto $\text{col}(X)$ is a projection of a vector to a single value.

\mathcal{P}_X = X(X'X)^{-1}X' = \tfrac{1}{n}\, XX' = \tfrac{1}{n} \begin{pmatrix} 1 & \cdots & 1 \\ \vdots & \ddots & \vdots \\ 1 & \cdots & 1 \end{pmatrix}

Then for any $y \in \mathbb{R}^n$,

\mathcal{P}_X y = \tfrac{1}{n} \begin{pmatrix} \displaystyle\sum_i y_i \\ \vdots \\ \displaystyle\sum_i y_i \end{pmatrix} = \bar{y}_n \cdot \mathbf{1}_n

$\square$

What happens if $\text{rank}(X) < p$ since in this case $(X'X)^{-1}$ does not exist?

Consider the orthogonal diagonalization $XX' = UDU' \in \mathbb{R}^{n \times n}$ and use the fact that $\text{col}(X) = \text{col}(\widetilde{U})$, where if $k := \text{rank}(X)$ then

\widetilde{U} := (U_1, \ldots, U_k)

and $U_i$ for $i \in \{1, \ldots, k\}$ corresponds to the $i$th nonzero eigenvalue of $XX'$. Then

\mathcal{P}_X = \widetilde{U}(\widetilde{U}'\widetilde{U})^{-1}\widetilde{U}' = \widetilde{U}\widetilde{U}'.

Gram-Schmidt Orthonormalization

Here we study a procedure for how to construct a set of mutually orthogonal vectors from a set of linearly independent vectors. Suppose that $x_1, \ldots, x_p \in \mathbb{R}^n$ are linearly independent. The result of the Gram-Schmidt orthonormalization procedure is a set of orthonormal vectors $u_1, \ldots, u_p \in \mathbb{R}^n$ such that

\text{span}\{u_1, \ldots, u_p\} = \text{span}\{x_1, \ldots, x_p\}.

In matrix form with $U := (u_1, \ldots, u_p)$ and $X := (x_1, \ldots, x_p)$, that is

\text{col}(U) = \text{col}(X).

How to construct $u_1$?

Take $u_1 := x_1$

Next, for $u_2$, modify $x_2$ so that it is orthogonal to $u_1$. Such a vector is in the orthogonal complement of the span of $u_1$,

Take $u_2 := (I_n - \mathcal{P}_{u_1})\, x_2$

Verify that

u_2' u_1 = x_2'(I_n - \mathcal{P}_{u_1})\, u_1 = 0.

Continuing on, the third vector must be orthogonal to both $u_1$ and $u_2$. That is,

Take $u_3 := (I_n - \mathcal{P}_{u_1} - \mathcal{P}_{u_2})\, x_3$

so that

\begin{aligned} u_3' u_2 &= x_3'(I_n - \mathcal{P}_{u_1} - \mathcal{P}_{u_2})\, u_2 \\ &= -x_3' \mathcal{P}_{u_1} u_2 \\ &= -x_3' u_1 (u_1' u_1)^{-1} u_1' u_2 \\ &= 0, \end{aligned}

and

\begin{aligned} u_3' u_1 &= x_3'(I_n - \mathcal{P}_{u_1} - \mathcal{P}_{u_2})\, u_1 \\ &= -x_3' \mathcal{P}_{u_2} u_1 \\ &= -x_3' u_2 (u_2' u_2)^{-1} u_2' u_1 \\ &= 0. \end{aligned}

Accordingly,

\begin{aligned} u_1 &= x_1 \\[.5em] u_2 &= x_2 - \frac{u_1 u_1' x_2}{\|u_1\|^2} \\[.5em] u_3 &= x_3 - \frac{u_1 u_1' x_3}{\|u_1\|^2} - \frac{u_2 u_2' x_3}{\|u_2\|^2} \\ &\vdots \end{aligned}

More concisely, for $j \in \{1, \ldots, n\}$,

u_j = x_j - \left(\sum_{k=1}^{j-1} \frac{u_k u_k'}{\|u_k\|^2}\right) x_j.

So now we have an orthogonal set $\{u_1, \ldots, u_p\}$. How to show that,

\text{span}\{u_1, \ldots, u_p\} = \text{span}\{x_1, \ldots, x_p\}\text{ ?}

Since this is equivalent to $\text{col}(U) = \text{col}(X)$, express

X = (u_1, \ldots, u_p) \cdot \underbrace{ \begin{pmatrix} 1 & \dfrac{u_1' x_2}{\|u_1\|^2} & \dfrac{u_1' x_3}{\|u_1\|^2} & \dfrac{u_1' x_4}{\|u_1\|^2} & \cdots & \dfrac{u_1' x_p}{\|u_1\|^2} \\[10pt] 0 & 1 & \dfrac{u_2' x_3}{\|u_2\|^2} & \dfrac{u_2' x_4}{\|u_2\|^2} & \cdots & \dfrac{u_2' x_p}{\|u_2\|^2} \\[10pt] 0 & 0 & 1 & \dfrac{u_3' x_4}{\|u_3\|^2} & \cdots & \dfrac{u_3' x_p}{\|u_3\|^2} \\[10pt] 0 & 0 & 0 & 1 & \cdots & \dfrac{u_4' x_p}{\|u_4\|^2} \\[6pt] \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\[4pt] 0 & 0 & 0 & 0 & \cdots & 1 \end{pmatrix} }_{=:\, S}

= \underbrace{\left(\frac{u_1}{\|u_1\|}, \ldots, \frac{u_p}{\|u_p\|}\right)}_{=:\, Q} \cdot \underbrace{\begin{pmatrix} \|u_1\| & & \\ & \ddots & \\ & & \|u_p\| \end{pmatrix} \cdot S}_{=:\, R}

Note that $X = QR$ is commonly called the “QR decomposition” for an orthogonal matrix $Q$, and an upper triangular matrix $R$. This expression also yields the Cholesky decomposition of $X'X$,

\[ X'X = R'Q'QR = R'R, \]

where $R$ is upper triangular. Think about the uniqueness of these decompositions.

The General Linear Model

Suppose $x_1, \ldots, x_p \in \mathbb{R}$ are covariates/features with a meaningful interpretation in some setting such that a random variable $Y$ is generated as

Y = \beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p + U

for some coefficients $\beta_0, \beta_1, \ldots, \beta_p \in \mathbb{R}$ and some random variable $U$ with $E(U) = 0$. This is regarded as the general linear model, and repeated sampling from it leads to a sample of size $n$, where

\begin{aligned} Y_1 &= \beta_0 + \beta_1 x_{1,1} + \cdots + \beta_p x_{p,1} + U_1 \\[.3em] &\quad\vdots \\[.3em] Y_n &= \beta_0 + \beta_1 x_{1,n} + \cdots + \beta_p x_{p,n} + U_n \end{aligned}

Equivalently,

Y = X\beta + U

where $Y \in \mathbb{R}^n$, $X \in \mathbb{R}^{n \times p}$, $\beta \in \mathbb{R}^p$, and $U \in \mathbb{R}^n$. Note that only the vector $Y$ is commonly understood as the data. The matrix $X$ is referred to as the design matrix, and $\beta$ is regarded as fixed but unknown.

The statistical inference problem is to estimate the unknown parameters $\beta$. A natural candidate for an estimate of $\beta$ is the more general idea of the least squares solution

\underset{\beta}{\operatorname{argmin}}\{Q(\beta)\} = \{\beta : X'X\beta = X'y\}

If we propose this estimator, then we need to study its sampling distribution. For example, if $U \sim N_n(0, \sigma^2 I_n)$, then

\hat{\beta} \sim N_p(\beta,\, \sigma^2(X'X)^{-1})

if $\text{rank}(X) = p$.

Next,

\hat{y} := X\hat{\beta} \in \text{col}(X).

If we decompose

y = X\hat{\beta} + \underbrace{y - X\hat{\beta}}_{=:\, \hat{e}},

it turns out that $\hat{e} \in \text{null}(X')$ since

X'\hat{e} = X'y - X'X\hat{\beta} = 0.

Moreover, note that this decomposition is unique because $\text{col}(X) \perp \text{null}(X')$, and

\begin{aligned} \|y\|^2 &= (X\hat{\beta} + \hat{e})'(X\hat{\beta} + \hat{e}) \\ &= \|X\hat{\beta}\|^2 + 2\hat{\beta}'X'\hat{e} + \|\hat{e}\|^2 \\ &= \|\hat{y}\|^2 + \|\hat{e}\|^2. \end{aligned}

The uniqueness of this decomposition/Pythagorean theorem implies that for the orthogonal projection matrix $P$ onto $\text{col}(X)$,

y = Py + (I - P)y = X\hat{\beta} + \hat{e}

with $Py = X\hat{\beta}$ and $(I - P)y = \hat{e}$.

Goals:

Determine whether certain functions of parameters are “estimable”.
Construct “unbiased estimators” for the “estimable” functions.

Definition. An estimator $t(y)$ is an unbiased estimator for the scalar $\lambda'\beta$ if and only if $E(t(Y)) = \lambda'\beta$, $\forall\, \beta$. $\square$

Definition. An estimator $t(y)$ is a linear estimator in $y$ if and only if $t(y) = c + a'y$, for some constants $a$ and $c$. $\square$

Definition. A function $\lambda'\beta$ is linearly estimable if and only if there exists a linear unbiased estimator of it. If no such estimator exists, then the function is called nonestimable. $\square$

Theorem. $\lambda'\beta$ is linearly estimable if and only if there exists a vector $a$ such that $E(a'Y) = \lambda'\beta$, $\forall\, \beta$. $\square$

Proof. First suppose that $\exists\, a \in \mathbb{R}^n$ such that $E(a'Y) = \lambda'\beta$, $\forall\, \beta \in \mathbb{R}^p$. Then $t(y) = a'y$ is an unbiased linear estimator of $\lambda'\beta$, $\forall\, \beta$, with $c = 0$. Conversely, if $\lambda'\beta$ is linearly estimable, then there exist constants $a$ and $c$ such that $t(y) = c + a'y$, and $\forall\, \beta$,

\lambda'\beta = E(t(Y)) = c + a'E(Y) = c + a'X\beta.

Choosing $\beta = 0$ demonstrates that $c = 0$, so that $\forall\, \beta \in \mathbb{R}^p$, $t(y) = a'y$ is a linear unbiased estimator of $\lambda'\beta$. $\square$

In words, $\lambda'\beta$ is linearly estimable if and only if it is equal to the expected value of a linear combination of the data. Observe the condition that

E(a'Y) = \lambda'\beta, \quad \forall\, \beta

implies that $\lambda = X'a$. So essentially, linear estimability of a function $\lambda'\beta$ is equivalent to the condition that $\lambda \in \text{col}(X')$.

Example. Suppose $Y = \mathbf{1}_n\beta + U$. Then for $\lambda = 1$, $\lambda'\beta = \beta$. Choosing $a = e_1 \in \mathbb{R}^n$,

E(a'Y) = a'\mathbf{1}_n\beta = \beta,

$\forall\, \beta \in \mathbb{R}$. That being true, $t(y) = y_1$ is a linear unbiased estimator of $\beta$. Probably $y_1$ is not a great estimator of $\beta$ due to high variability, and the discard of the information $y_2, \ldots, y_n$, but the theorem only gives the necessary and sufficient conditions for the existence of a linear estimator. A better choice would be $a = \tfrac{1}{n} \cdot \mathbf{1}_n$. Then $t(y) = a'y = \bar{y}_n$, and

E(\bar{y}_n) = E(a'Y) = \tfrac{1}{n}\mathbf{1}_n'\mathbf{1}_n\beta = \beta.

$\square$

Example. Suppose

Y = \begin{pmatrix} 1 & 1 \\ 1 & 1 \\ 1 & 0 \\ 1 & 0 \end{pmatrix} \underbrace{\begin{pmatrix} \mu \\ \alpha \end{pmatrix}}_{=:\, \beta} + U

To estimate $\alpha$, choose $\lambda = (0, 1)'$ since $\lambda'\beta = \alpha$. For $\alpha$ to be estimable, it must be the case that $\lambda \in \text{col}(X')$, so determine if there exists a solution to

\begin{pmatrix} 1 & 1 & 1 & 1 \\ 1 & 1 & 0 & 0 \end{pmatrix} \begin{pmatrix} a_1 \\ a_2 \\ a_3 \\ a_4 \end{pmatrix} = \begin{pmatrix} 0 \\ 1 \end{pmatrix}.

One solution is $a = (1, 0, -1, 0)'$, and so $t(y) = y_1 - y_3$ is a linear unbiased estimator of $\alpha$.

$\square$

Example. Now let’s consider an example where $\lambda'\beta$ is not linearly estimable.

Y = \begin{pmatrix} 1 & 1 & 0 \\ 1 & 1 & 0 \\ 1 & 0 & 1 \\ 1 & 0 & 1 \end{pmatrix} \underbrace{\begin{pmatrix} \mu \\ \alpha_1 \\ \alpha_2 \end{pmatrix}}_{=:\, \beta} + U

In this case, $\lambda = (0, 1, 0)'$ to estimate $\alpha_1$, but

\begin{pmatrix} 1 & 1 & 1 & 1 \\ 1 & 1 & 0 & 0 \\ 0 & 0 & 1 & 1 \end{pmatrix} \begin{pmatrix} a_1 \\ a_2 \\ a_3 \\ a_4 \end{pmatrix} = \begin{pmatrix} 0 \\ 1 \\ 0 \end{pmatrix}

does not have a solution. Thus, there exists no linear unbiased estimator for $\alpha_1$ because $\lambda \neq X'a$ for any $a \in \mathbb{R}^n$.

$\square$

Notice that the set of vectors $\lambda$ for which $\lambda'\beta$ is linearly estimable actually forms an entire subspace,

\{\lambda \in \mathbb{R}^p : \lambda = X'a \text{ for some } a \in \mathbb{R}^n\} = \text{col}(X').

Moreover, if $\text{rank}(X') = p$, then $\text{col}(X') = \mathbb{R}^p$ and so $\lambda'\beta$ is linearly estimable $\forall\, \lambda \in \mathbb{R}^p$. An important corollary of this discussion is that each component of $\beta$ is linearly estimable (i.e., identifiable) if $X$ has full column rank.

With the notion of linear estimability now established, what are some effective strategies for determining if a function $\lambda'\beta$ is linearly estimable?

Method 1. If $\lambda \in \text{col}(X')$, then $\lambda'\beta$ is estimable. Accordingly, for any basis $\{v_1, \ldots, v_k\}$ for $\text{col}(X')$ determine if the system $v_1 a_1 + \cdots + v_k a_k = \lambda$ has a solution. Equivalently, since $\mathcal{P}_{X'} = X'(XX')^{-1}X$ projections onto $\text{col}(X')$, it suffices to show that either $\mathcal{P}_{X'}\lambda = \lambda$

Method 2. Since $\text{col}(X') \oplus \text{null}(X) = \mathbb{R}^n$, $\lambda \in \text{col}(X')$ if and only if $\lambda \perp \text{null}(X)$. That being so, construct a basis for $\text{null}(X)$, $\{w_1, \ldots, w_{p-k}\}$, and show that

\lambda' w_j = 0 \quad \forall\, j \in \{1, \ldots, p-k\}.

Method 3. If $\lambda'\beta$ can be expressed as a linear combination of $E(Y_1), \ldots, E(Y_n)$, then $\lambda'\beta$ is estimable.

Theorem. If $\lambda'\beta$ is estimable, then the least squares estimator $\lambda'\hat{\beta}$ is the same for all solutions to the normal equations. $\square$

Proof. Let $\hat{\beta}_1, \hat{\beta}_2 \in \{X'X\beta = X'y\}$. Then

\begin{aligned} X'X\hat{\beta}_1 &= X'X\hat{\beta}_2 \\ X'X(\hat{\beta}_1 - \hat{\beta}_2) &= 0, \end{aligned}

and so $\hat{\beta}_1 - \hat{\beta}_2 \in \text{null}(X'X) = \text{null}(X)$. For estimable $\lambda'\beta$, $\lambda \perp \text{null}(X)$ so it follows that $\lambda'(\hat{\beta}_1 - \hat{\beta}_2) = 0$. $\square$

Theorem. The least squares estimator $\lambda'\hat{\beta}$ of an estimable function $\lambda'\beta$ is a linear unbiased estimator of $\lambda'\beta$. $\square$

Proof. First observe that for some $a \in \mathbb{R}^n$ and $\lambda \in \mathbb{R}^p$,

\lambda'\hat{\beta} = a'X(X'X)^{-1}X'y = a'\mathcal{P}_X y.

Hence, $\lambda'\hat{\beta} = v'y$ for some $v \in \mathbb{R}^n$. Next,

E(\lambda'\hat{\beta}) = a'\mathcal{P}_X E(y) = a'\mathcal{P}_X X\beta = a'X\beta = \lambda'\beta.

$\square$

Gauss-Markov Model

The motivation here is to understand additional inferences that ensue from making assumptions on the variance of the errors in the general linear model. We begin with the Gauss-Markov assumptions/model,

Y = X\beta + U,

with $E(U) = 0 \in \mathbb{R}^n$, $\text{Var}(U) = \sigma^2 I_n$, and $\sigma \in \mathbb{R}_+$. In particular, $E(U_i) = 0$ for $i \in \{1, \ldots, n\}$, and the errors are homoskedastic; that is,

\text{Cov}(U_i, U_j) = \begin{cases} \sigma^2 & \text{for } i = j \\ 0 & \text{for } i \neq j \end{cases}.

Example. Suppose that $\lambda'\beta$ is an estimable function. Then for $\hat{\beta} \in \{X'X\beta = X'y\}$,

\begin{aligned} \text{Var}(\lambda'\hat{\beta}) &= \lambda'\,\text{Var}\!\left[(X'X)^{-1}X'Y\right]\lambda \\ &= \lambda'(X'X)^{-1}X'\,\text{Var}(Y)\,X(X'X)^{-1}\lambda \\ &= a'X(X'X)^{-1}X'\,\text{Var}(Y)\,X(X'X)^{-1}X'a, \quad \text{for some } a, \text{ since } \lambda \in \text{col}(X') \\ &= \sigma^2 \cdot a'X(X'X)^{-1}X'X(X'X)^{-1}X'a \\ &= \sigma^2 \cdot a'X(X'X)^{-1}X'a \\ &= \sigma^2\,\lambda'(X'X)^{-1}\lambda. \end{aligned}

$\square$

It turns out that the least squares estimate plays an important role under the Gauss-Markov assumptions, as described by the following result.

Theorem. (Gauss-Markov theorem) Under the Gauss-Markov assumptions, if $\lambda'\beta$ is estimable, then $\lambda'\hat{\beta}$ is the best (minimum variance) linear unbiased estimator (BLUE) of $\lambda'\beta$, where $\hat{\beta} \in \{X'X\beta = X'y\}$. $\square$

Proof. Assume that $c + d'y$ is another unbiased estimator of $\lambda'\beta$. Then, $\forall\, \beta \in \mathbb{R}^p$,

\lambda'\beta = E(c + d'Y) = c + d'X\beta.

Choosing $\beta = 0$ demonstrates that $c = 0$. Accordingly, $\lambda = X'd$, and

\begin{aligned} \text{Var}(c + d'Y) &= \text{Var}(d'Y) \\ &= \text{Var}(\lambda'\hat{\beta} + d'Y - \lambda'\hat{\beta}) \\ &= \text{Var}(\lambda'\hat{\beta}) + 2\,\text{Cov}(\lambda'\hat{\beta},\, d'Y - \lambda'\hat{\beta}) + \text{Var}(d'Y - \lambda'\hat{\beta}). \end{aligned}

Next,

\begin{aligned} \text{Cov}(\lambda'\hat{\beta},\, d'Y - \lambda'\hat{\beta}) &= \text{Cov}\!\left(\lambda'(X'X)^{-1}X'Y,\; d'Y - \lambda'(X'X)^{-1}X'Y\right) \\ &= \lambda'(X'X)^{-1}X'\cdot \text{Var}(Y)\cdot\left(d - X(X'X)^{-1}\lambda\right) \\ &= \sigma^2\,\lambda'(X'X)^{-1}X'\left(d - X(X'X)^{-1}\lambda\right) \\ &= 0. \end{aligned}

Hence,

\text{Var}(c + d'Y) = \text{Var}(\lambda'\hat{\beta}) + \text{Var}(d'Y - \lambda'\hat{\beta}) \geq \text{Var}(\lambda'\hat{\beta})

with equality if and only if

\begin{aligned} 0 &= \text{Var}(d'Y - \lambda'\hat{\beta}) \\ &= \left(d' - \lambda'(X'X)^{-1}X'\right)\text{Var}(Y)\left(d' - \lambda'(X'X)^{-1}X'\right)' \\ &= \sigma^2\left\lVert d - X(X'X)^{-1}\lambda\right\rVert_2^2. \end{aligned}

That is, $\text{Var}(c + d'Y) = \text{Var}(\lambda'\hat{\beta})$ if and only if $d = X(X'X)^{-1}\lambda$, in which case

d'y = \lambda'(X'X)^{-1}X'y = \lambda'\hat{\beta},

so that $\lambda'\hat{\beta}$ is the unique BLUE of $\lambda'\beta$. $\square$

Variance Estimation

Recall the unique decomposition,

\begin{aligned} Y &= \mathcal{P}_X Y + (I - \mathcal{P}_X)Y \\ &= X\hat{\beta} + (I - \mathcal{P}_X)(X\beta + U) \\ &= X\hat{\beta} + (I - \mathcal{P}_X)U \end{aligned}

for any $\hat{\beta} \in \{X'X\beta = X'Y\} = \{X\beta = \mathcal{P}_X Y\}$. Since $\mathcal{P}_X Y$ is related to the estimation of $\beta$, in the decomposition, the intuition is that $(I - \mathcal{P}_X)Y$ is related to the estimation of $\sigma^2$.

Theorem. Under the Gauss-Markov assumptions, $\widehat{\sigma^2} := Y'(I - \mathcal{P}_X)Y / (n - r)$ is an unbiased estimator of $\sigma^2$, where $r := \text{rank}(X) \leq p$. $\square$

Proof.

\begin{aligned} E(\widehat{\sigma^2}) &= \frac{1}{n-r}\,E\!\left[(X\beta + U)'(I - \mathcal{P}_X)(X\beta + U)\right] \\ &= \frac{1}{n-r}\,E\!\left[U'(I - \mathcal{P}_X)U\right] \\ &= \frac{1}{n-r}\,E\!\left[\text{tr}\!\left(U'(I - \mathcal{P}_X)U\right)\right] \\ &= \frac{1}{n-r}\,\text{tr}\!\left[(I - \mathcal{P}_X)E(UU')\right] \\ &= \frac{1}{n-r}\,\text{tr}(I - \mathcal{P}_X)\cdot\sigma^2 \\ &= \sigma^2. \end{aligned}

$\square$

Note that we have not shown that $\sqrt{\widehat{\sigma^2}}$ is an unbiased estimator for $\sigma$.

Implications on Model Selection

Here we make a distinction between the true data-generating model and the model being fit. In the context of a linear model, two discrepancies are:

Underfitting (missing covariates)
Overfitting (redundant covariates)

Misspecification from Underfitting

Consider the true model $Y = X\beta + \eta + U$, against the model used by the practitioner,

Y = X\beta + U.

Using the practitioner’s model, $Y = X\beta + U$, under the Gauss-Markov assumptions results in an omitted-variable bias. Consider the effect on the least squares estimator. For $\lambda \in \text{col}(X')$, where $\lambda = X'a$ for some $a \in \mathbb{R}^n$,

\begin{aligned} E(\lambda'\hat{\beta}) &= \lambda'(X'X)^{-1}X'E(Y) \\ &= \lambda'(X'X)^{-1}X'(X\beta + \eta) \\ &= \lambda'\beta + \lambda'(X'X)^{-1}X'\eta \\ &= \lambda'\beta + \underbrace{a'\mathcal{P}_X\eta}_{\text{misspecification bias}} \end{aligned}

Observe that if $\eta \perp \text{col}(X)$, then $\lambda'\hat{\beta}$ is unbiased. How about estimation of $\sigma^2$?

\begin{aligned} E\!\left(Y'(I - \mathcal{P}_X)Y\right) &= (X\beta+\eta)'(I - \mathcal{P}_X)(X\beta+\eta) + \sigma^2\,\text{tr}(I - \mathcal{P}_X) \\ &= \eta'(I - \mathcal{P}_X)\eta + \sigma^2(n-r) \end{aligned}

That is, if $\eta \in \text{col}(X)$ then $\widehat{\sigma^2}$ is unbiased. In the case that $\eta \in \text{col}(X)$, $\exists\, \delta$ such that

Y = X\beta + \eta + U = X(\beta+\delta) + U.

Overfitting and Multicollinearity

Practitioner’s model: $Y = X\beta + W\gamma + U$, where $\gamma = 0$ in the true data-generating model. Assume that the concatenated matrix $(X, W)$ has full column rank. Then under the Gauss-Markov assumptions,

\begin{pmatrix} \hat{\beta} \\ \hat{\gamma} \end{pmatrix} = \begin{pmatrix} X'X & X'W \\ W'X & W'W \end{pmatrix}^{-1}\begin{pmatrix} X'Y \\ W'Y \end{pmatrix}, \quad \text{with} \quad E\begin{pmatrix} \hat{\beta} \\ \hat{\gamma} \end{pmatrix} = \begin{pmatrix} \beta \\ \gamma \end{pmatrix} = \begin{pmatrix} \beta \\ 0 \end{pmatrix}.

It can be worked out that,

\text{Var}(\hat{\beta}) = \sigma^2(X'X)^{-1} + \sigma^2(X'X)^{-1}X'W\left(W'(I-\mathcal{P}_X)W\right)^{-1}W'X(X'X)^{-1},

where the second term can be understood as the increase in variability of the least squares estimator that results from including redundant covariates. However, estimation of $\sigma^2$ remains unbiased:

E\!\left[Y'(I - \mathcal{P}_{(X,W)})Y\right] = \sigma^2\left(n - \text{rank}(X,W)\right).

Notice that if the columns of $W$ are not only linearly independent but are also orthogonal to the columns of $X$, then the variance penalty is zero. In this case, there is no cost to overfitting. More generally, though, if $W$ is not orthogonal to $X$, then a loss of efficiency from overfitting ensues.

A measure to gauge the influence of multicollinearity relates to the SVD of $(X, W)$.

Definition. The mean squared error (MSE) of an estimator $\hat{\theta}$ for some parameter $\theta$ is

\text{MSE}(\hat{\theta}) := E\left[\lVert \hat{\theta} - \theta \rVert_2^2\right].

$\square$

For some unbiased estimator $\hat{\beta}$ of $\beta$,

\begin{aligned} \text{MSE}(\hat{\beta}) &= E\left[\lVert \hat{\beta} - \beta \rVert_2^2\right] \\ &= \text{tr}\!\left(E\left[(\hat{\beta}-\beta)(\hat{\beta}-\beta)'\right]\right) \\ &= \text{tr}\!\left(\text{Var}(\hat{\beta})\right). \end{aligned}

Further, denoting the SVD of the design matrix as $X = U\Sigma V'$,

\begin{aligned} \text{MSE}(\hat{\beta}) &= \text{tr}\!\left(\text{Var}(\hat{\beta})\right) \\ &= \sigma^2\,\text{tr}\!\left[(X'X)^{-1}\right] \\ &= \sigma^2\,\text{tr}\!\left(V\Sigma^{-2}V'\right) \\ &= \sigma^2\sum_{i=1}^{p}\frac{1}{\lambda_i}, \end{aligned}

where $\lambda_i$ is the $i$th eigenvalue of $X'X$. In the case that multicollinearity is severe, $X'X$ is nearly singular. As such, $\min_i\{\lambda_i\}$ will be close to zero, and so $\text{MSE}(\hat{\beta})$ will be excessively large. Another measure of multicollinearity is the condition number,

\kappa := \frac{\max_i\{\lambda_i\}}{\min_i\{\lambda_i\}}.

Alternatively, there is a residual-sum-of-squares approach to assessing the severity of multicollinearity, called variance inflation factors.