I found these posts particularly helpful:

https://stats.stackexchange.com/questions/46151/how-to-derive-the-least-square-estimator-for-multiple-linear-regression

https://stats.stackexchange.com/questions/134282/relationship-between-svd-and-pca-how-to-use-svd-to-perform-pca

http://www.math.miami.edu/~armstrong/210sp13/HW7notes.pdf

If $X$ is an $n \times p$ matrix then the matrix $X(X^TX)^{-1}X^T$ defines a *projection* onto the column space of $X$.  Intuitively, you have an overdetermined system of equations, but still want to use it to define a linear map $\mathbb{R}^p \rightarrow \mathbb{R}$ that will map $X$ to something resembling $y$.  So we settle for sending $X$ to the closest thing to $y$ that we can be expressed as a linear combination of your features (the columns of $X$).  

As far as an interpretation of $(X^TX)^{-1}$, I don't have an amazing answer yet.
I know you can think of $(X^TX)$ as basically being the covariance matrix of the dataset.