Your question is a bit unclear to me. However, as far I can understand...

> if the generalized linear model is
> $y(x_1,x_2)=b_0+b_1x_1+b_2x_2+b_3x_1x_2+b_4x^2_1+b_5x^2_2$, and there
> are, say, four data points, then

...then your model matrix $X$ is a $4\times 6$ matrix:
$$\begin{bmatrix}
1 & x_{11} & x_{21} & x_{11}x_{21} & x_{11}^2 & x_{11}^2 \\
1 & x_{12} & x_{22} & x_{12}x_{22} & x_{12}^2 & x_{12}^2 \\
1 & x_{13} & x_{23} & x_{13}x_{23} & x_{13}^2 & x_{13}^2 \\
1 & x_{14} & x_{24} & x_{14}x_{24} & x_{14}^2 & x_{14}^2
\end{bmatrix}$$

If $\text{rank}(X)=4$, then $X^TX$ is a $6\times 6$ singular matrix, while $XX^T$ is a $4\times 4$ non singular matrix.

You should need $(X^TX)^{-1}$ to estimate $\beta$, but since $X^TX$ is singular you have to use a _right inverse_, i.e. a matrix $X_R$ such that $XX_R=I$:
$$ X\beta=y,\qquad X\beta=XX_Ry,\qquad \beta=X_Ry$$
A right inverse that always guarantees a solution is: $X_R=X^T(XX^T)^{-1}$.
See [Wikipedia](https://en.wikipedia.org/wiki/Generalized_inverse) and Cherkassky & Mulier, _Learning From Data: Concepts, Theory, and Methods_, John Wiley & Sons, 2007, [Appendix B](https://onlinelibrary.wiley.com/doi/pdf/10.1002/9780470140529.app2).