• Native close form solution: just $(X'X)^{-1}(X'Y)$. We can always solve that inverse matrix. It works fine for small dataset but to inverse a large matrix might be very computations expensive.
• QR decomposition: this is the default method in R when you run lm(). In short, QR decomposition is to decompose the X matrix to XQR where Q is an orthogonal matrix and R is an upper triangular matrix. Therefore, $X'X \beta = X'Y$ and then $R'Q'QR \beta = R'Q'Y$ then $R'R \beta = R'Q'Y$then $'R \beta = Q'Y$. Because R is a upper triangular matrix then we can get beta directly after computing Q'Y. (More details)