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.

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.
