What estimation methods other than ordinary least squares guarantee the orthogonality of predictions and residuals? From my question here, it is evident that estimation approaches to linear regression other than ordinary least squares can result in the predictions and residuals lacking orthogonality, despite the model being linear.
What, if any, approaches to estimating the $\beta$ of $y=X\beta+\epsilon$ are there that are not equivalent to ordinary least squares (yield a different answer than $\hat\beta=(X^TX)^{-1}X^Ty$) yet still yield this orthogonality?
Let’s rule out $\hat\beta=\vec 0$. If it happens that $X\hat\beta=X\hat\beta_{ols}$, so be it, but that is not a requirement.

For the benefit of completeness, I'll note that you have ruled out the case where $\hat{\boldsymbol{\beta}}=\mathbf{0}$, which gives the zero vector for the predicted response.  The zero vector is (trivially) orthogonal to any other vector, so this estimator also gives a predicted response vector that is (trivially) orthogonal to the residual vector.
Setting aside this special case, or mixtures with this case, the hat matrix from OLS is the unique projection matrix that gives an orthogonal projection onto the column space of the design matrix.  This projection uniquely defines OLS estimation (equivalently MLE under the Gaussian linear model) and so there are no other estimation methods that use this orthogonal projection that are not equivalent to OLS.  This means that the only way you can get a guarantee of (non-trivial) orthogonality is to use OLS estimation.
