I found these posts particularly helpful:
How to derive the least square estimator for multiple linear regression?
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.