Understanding multiple linear regression residuals I'm a bit confused whether residual vectors in OLS are orthogonal to every vector in X.
The problem I have is: Which of the following are true statements for our multiple linear regression with datapoints as rows of design matrix X, label vector Y, and using the model Ŷ = X̂.
I can understand that the residual vector e = Y - Ŷ is orthogonal to Y. And to me, it would make sense for the residual vector e = Y - Ŷ to also be orthogonal to some cases of vector in X. But what about every?
Any help would be appreciated.
I reference this diagram:

 A: Its been a long time since I've had to do some linear algebra, so forgive me if I've forgotten some of it.  Let's begin by assuming that $X$ is full rank and $\operatorname{rnk}(X)=p<n$.
Minimizing the loss means
$$ \nabla R = 0 \>.$$
Note that $\nabla R$ in terms of $X$ and $y$ is
$$ \nabla R = -2X^T (y- X\theta) $$
Note that $e = y - X\theta$.  This means that $y-X\theta$ is in the null space of $X^T$, so the residual is orthogonal to the columns of $X$.  Those columns span some subspace of $\mathbb{R}^n$ (the blue plane in your image).  So since the residual is orthogonal to the columns, and the columns are linearly independent, it is orthogonal to all linear combinations of those columns. In particular, it is orthogonal to the prediction (the closest vector in the plane to $y$).
A: Under OLS the residual is orthogonal to every column in the design matrix
To understand this issue, it is worth understanding the concept of the column space of the design matrix $\mathbf{x}$.  This is the space spanned by linear combinations of the column vectors in the design matrix, so it represents every possible vector you can form via the linear transformation $\mathbf{x} \boldsymbol{\theta}$.  In your diagram, the column space is the blue region.
When you fit the regression model using OLS estimation, this method is equivalent to forming the predicted response as the projection of $\mathbf{y}$ onto the column space of $\mathbf{x}$.  In fact, the hat matrix $\mathbf{h} = \mathbf{x} (\mathbf{x}^\text{T} \mathbf{x})^{-1} \mathbf{x}^\text{T}$ is the projection matrix that maps the actual response vector onto the predicted response vector.  This orthogonal projection is shown in the diagram in your question.  Since the fitted model is formed by an orthogonal projection, the residual vector is orthogonal to the column space of $\mathbf{x}$ which means it is orthogonal to every column of $\mathbf{x}$ (and therefore also orthonogal to the predicted response vector $\hat{\mathbf{y}}$).
A: Assuming $X$ is full rank, you can see it in terms of the hat matrix $H = X(X^\top X)^{-1}X^\top$. This matrix is an orthogonal projection matrix onto the column space of $X$.
The estimate of $y$ is $\hat{y}=X\hat{\theta}=Hy$. The residual is $e = y - \hat{y} = (I-H)y$. The matrix $(I-H)$ is an orthogonal projection matrix onto the orthogonal complement of the column space of $X$, so the residual vector $e$ is orthogonal to every vector in the column space of $X$.
