Parameter estimate in linear regression model with singular matrix $X^TX$

Question

I started to study linear regression this autumn and got stuck at the estimation of regression parameters (in matrix form).

In parameter estimation of the general linear model, the OLS method is used to estimate the parameter $\beta$ and result in least square estimator $\hat{\beta}$ in normal equation $$X^\top X\hat{\beta} = X^\top Y$$

I tried to figure out how to find the unique solution of $\hat{\beta}$ when $X^\top X$ has full rank and is non-singular. However, as I am a bit rusty on linear algebra, I am really confused about the methods to acquire multiple solutions of $\hat{\beta}$ for singular $X^\top X$. Therefore I started to try some simple examples.

For instance, I have a linear regression model likes this(for the $i_{th}$ response): $$y_i = \beta_0 + \beta_1\psi(x_i) +e_i, i = 1, ..., n,$$ where $\psi(x_i)$ is the function of $x_i$ and the $e_i$ are independent and have zero mean.

I am now considering the situation that $x_i$ satisfies the case that $\psi(x_i) = 0 \ \forall \ i$. $$ X^\top X = \begin{pmatrix} n&0 \\ 0&0 \end{pmatrix} $$

The matrix above tends to be singular and I can not get a unique solution. And I am trying to find the possible solutions of $\hat{\beta}$.

I am now consulting methods like Gaussian elimination and Cholesky decomposition. As I am rusty on linear algebra, I still could not get a clue. Could anyone give me some inspiration for the simple case above? Thank you.

Answer 1

The question asks

The matrix above tends to be singular and I can not get a unique solution. And I am trying to find the possible solutions of $\hat{\beta}$.

The solution set to the specific problem in OP's post is to choose $\beta_0= \bar{Y}$ and $\beta_1$ can be any real number. To see why, expand the $X^\top X \beta$ and $X^\top Y$ products. $$\begin{align} X^\top X \beta &= \begin{bmatrix} n&0 \\ 0&0 \end{bmatrix}\beta = X^\top Y\\ &=\begin{bmatrix}n\beta_0 + 0\beta_1 \\ 0\beta_0 + 0\beta_1\end{bmatrix} = X^\top Y\\ &=\begin{bmatrix}n\beta_0 \\ 0\end{bmatrix} = \begin{bmatrix} 1^\top Y \\ 0^\top Y \end{bmatrix} \end{align}$$ You'll see you have two sums that do not depend on each other: we have $0 \beta_1 = 0$ and $n \beta_0 = \sum_{i=1}^n Y_i$, from which the conclusion is obvious.

A ridge or lasso regularization method can provide a unique solution given a choice of $\lambda > 0$. These methods provide solutions to a modified version of the original problem, but do not provide the set of all possible solutions to OP's stated problem.

Likewise, a pseudo-inverse could be computed using SVD. However, this will give the solution $\beta^*$, which is the solution with the smallest Euclidean norm, from among all $\beta$ that can solve the under-determined system of equations exactly. But the least-norm solution $\beta^*$ is not the only $\beta$ that can solve the equation exactly.

You mention Gaussian elimination. In this example, there's essentially no more work to be done: you've found the first column has 1 pivot, and the second column has 0 pivots. In a more complex example of a rank-deficient regression, you might have to undertake more manipulations.

Answer 2

In the classical linear regression, we estimate the regression coefficients $\beta$ by Ordinary Least Squares (OLS) estimator

$$ \underset{\beta \in \mathbb{R}^p}{\operatorname{arg\,min}} \sum_{i=1}^{n} \left(Y_{i} - \sum_{j=1}^{p}X_{ij}\beta_{j}\right)^{2}, \quad(*) $$

which has a closed-form expression

$$ \hat\beta = \left(X^\top X\right)^{-1} X^\top Y. $$ Here $X$ is the design matrix and $Y$ is the vector of response values for the $n$ individuals.

As you correctly point out, OLS is computable only if $X^\top X$ is invertible. If that's not the case, a huge number of approaches are available to deal with this issue.

One way to go would be to find a generalized inverse for $X^\top X$, through the singular value decomposition (see Michael M's comment). However, as highlighted by Sycorax's answer, this will not necessarily find all the possible solutions.

Another approach could be to use a penalized estimation problem, although your solution will depend on a tuning parameter.

In particular, the idea is to replace $(*)$ by the following penalized estimation problem

$$ \underset{\beta \in \mathbb{R}^p}{\operatorname{argmin}} \sum_{i=1}^{n} \left(Y_{i} - \sum_{j=1}^{p}X_{ij}\beta_{j}\right)^{2} \color{}{+ \lambda \sum_{j=1}^{p}\phi(\beta_{j})}. $$

The function $\phi(\beta_j)$, which makes the minimisation problem possible, is the penalty term. There are many possibilities of $\phi$ proposed in the literature, but the two most popular are, perhaps, the lasso penalty $$\phi(\beta_j) = |\beta_j|$$

and the ridge penalty

$$ \phi(\beta_j) = \beta_j^2. $$

$\lambda$ is called regularisation parameter and determines the amount of penalty; it is typically unknown and has to be estimated from data; see e.g. here.

You can learn deeper about these topics also in this website, there are a lot of interesting answers; see for instance here and here among others.