Why does regularization wreck orthogonality of predictions and residuals in linear regression? Following up on this question...
In ordinary least squares, the predictions and residuals are orthogonal. $$\sum_{i=1}^n\hat{y}_i (y_i - \hat{y}_i) = 0$$
If we estimate the regression coefficients using some other method but the same model, such as using regularization, why, intuitively, should that wreck the orthogonality?
 A: An image might help. In this image, we see a geometric view of the fitting.

*

*Least squares finds a solution in a plane  that has the closest distance to the observation.
(more general a higher dimensional plane for multiple regressors and a curved surface for non-linear regression)
In this case, the vector between observation and solution is perpendicular to the plane (a space spanned be the regressors), and perpendicular to the regressors.


*Regularized regression finds a solution in a restricted set inside the the plane that has the closest distance to the observation.
In this case, the vector between observation and solution is not anymore perpendicular to te plane and not anymore perpendicular to the regressors.
But, there is still some sort of perpendicular relation, namely the vector of the residuals is in some sense perpendicular to the edge of the circle (or whatever other surface that is defined by te regularization)

The model of $\hat{y}$
Our model gives estimates of the observations, $\hat{y}$, the observations as function of parameters $\beta_i$.
$$\hat{y} = f(\beta)$$
In our image this is a linear function with two parameters $\beta_0$ and $\beta_1$
(you can of course generalize this to a large size of coefficients and observations, for simplicity we regard three observations and two coefficients such that we can plot it)
$$\begin{bmatrix} \hat{y}_{1} \\  \hat{y}_{2} \\  \hat{y}_{3} \end{bmatrix} = \beta_0 \begin{bmatrix} 1 \\  1 \\ 1 \end{bmatrix}+
\beta_1 \begin{bmatrix} 0 \\  1 \\ 2 \end{bmatrix}$$
The possible solutions of the model, defined by this linear sum, is represented by the red plane in the image.
Note that this plane in the image relates to the possible solutions of $y_i = \beta_0 + \beta_1 x_i$, when $[x_1,x_2,x_3] = [0,1,2]$. So we plotted the space of all possible $y_i$ (which is a 3D-space, and more generally a n-dimensional space) and the possible solutions that the model allows is a plane inside this space.
Finding the best model with least squares
The model allows any solution in the plane spanned by the model (in the image this is the 2D red plane, in general this can be a higher dimensional plane, and also it does not need to be linear).
The least-squares method will select the 'solution' $\hat{y} = \hat\beta_0 + \hat\beta_1 x_1 $ that has the lowest difference in terms of the squares of the residuals.
In geometric terms, this is equal to finding the point in the plane that has the smalles euclidian distance to the observed value. This smallest difference is achieved when the vector of residuals  is orthogonal to the plane.
Finding the best model with ridge regression (or other regularization)
When we apply a penalty then this is similar to applying some constraint like 'the sum of vectors can not be above some value'. In the image this is represented by the purple drawing.
The solution is still inside the plane, but also inside the circle. Now the estimate solution is still a representing a shortest distance between the space of solutions and the observation. But the optimal solution is not anymore orthogonal projection onto the red plane. It is instead the shortes distance to the purple circles.
A: Think in geometrical terms: the OLS fit is the projection of $Y$ on the space spanned by the columns of $X$, hence the residual vector is orthogonal to that space. If you regularize, performing ridge regression or otherwise, you will in general move your fit away from the projection and destroy orthgonality.
A book which I learned most of this from and is (in my opinion) difficult to surpass is Seber, G.A.F. Linear Regression Analysis, Wiley. I used the ca. 1980 edition, but there is a newer version. Grab a copy if you can.
A: One way to derive the least squares estimate of $\beta$ (the vector of regression coefficients) is that it is the one and only value of $\beta$ (*) that would make the error vector orthogonal to every predictor, and hence orthogonal to linear combinations of the predictors, which is what the predicted values $\hat{y}$ are. See, for example, section 2 of these course notes.
From that perspective, we can see that any estimate of $\beta$ other than the least squares estimate (*) -- for example, any estimate that has been regularized toward 0 to some extent -- will not have this orthogonality property.
(*) Aside from $\hat{\beta} = 0$, which leads to $\hat{y}=0$, which is always orthogonal to any other vector
A: I wrote a comprehensive explanation on this question in my site.
It might be useful for readers.

I'll talk about the ridge regularization here because it can be shown to neatly use the same equations used to derive the OLS solution (see this answer).
The coefficients in ridge regression (with penalty weighting $\lambda$) are simply:
$$\beta = (X^TX+\lambda\mathbb I)^{-1}X^Ty$$
The solution to the OLS can be obtained just as well by setting $\lambda = 0$.
The use of the normal equations to the ridge problem can be recovered from and correspond to an augmentation of $X$.
Concatenating new virtual samples formed by a identity matrix:
$$
\matrix{
X_\text{new}=\left[\matrix{
X_\text{old} \\ \sqrt{\lambda}\mathbb I_{p\times p}
}\right]
\qquad
Y_\text{new}=\left[\matrix{
Y_\text{old} \\ \mathbf 0_{p\times1}
}\right]
}$$
If we do that, it can be quite straightforwardly shown that:
$$\beta = (X_\text{old}^TX_\text{old}+\lambda\mathbb I)^{-1}X^T_\text{old} y_\text{old} = (X_\text{new}^TX_\text{new})^{-1}X_\text{new}^T y_\text{new}$$
Thus, since we are using the normal equations to derive the solution to ridge regression, the property of orthogonal residuals and predictions is kept intact.
But notice that, now, predictions involve these virtual samples.
That's why, when looking only at the real samples, this orthogonality is not guaranteed: you are missing part of the puzzle by not taking into account these "virtual" samples.
