Parameter estimate in linear regression model with singular matrix $X^TX$ 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.

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.
